text:
function evaluate_results(incidents::DataFrame,
                          sim_data::DataFrame,
                          ressource_flow::Array{Int64,3},
                          simulation_capacity::Array{Int64,2})
# Evaluate the results of the simulation
    incd = innerjoin(incidents, sim_data, on=:incidentid)
    for i = 1:nrow(incd)
        if ismissing(incd[i,:dispatch_minute_first])
            incd[i,:dispatch_minute_first] = ceil(incd[i,:incident_minute] + drop_incident)
            incd[i,:arrival_minute_first]  = ceil(incd[i,:dispatch_minute_first] + max_drive)
        end
        if ismissing(incd[i,:dispatch_minute_all])
            incd[i,:dispatch_minute_all] = ceil(incd[i,:incident_minute] + drop_incident)
            incd[i,:arrival_minute_all]  = ceil(incd[i,:dispatch_minute_all] + max_drive)
        end
    end
    incd[:,:response_time_first]   = incd[:,:arrival_minute_first] - incd[:,:incident_minute]
    incd[:,:response_time_all]     = incd[:,:arrival_minute_all] -  incd[:,:incident_minute]
    incd[:,:dispatch_time_first]   = incd[:,:dispatch_minute_first] - incd[:,:incident_minute]
    incd[:,:dispatch_time_all]     = incd[:,:dispatch_minute_all] - incd[:,:incident_minute]
    incd[:,:driving_time_first]    = incd[:,:arrival_minute_first] - incd[:,:dispatch_minute_first]

# Group the results for evaluation
# 1. group the results according to the location and weekhour
    weekly_location = groupby(dropmissing(select!(copy(incd), Not([:priority]))), [:location_responsible, :weekhour])
    weekly_location = combine(weekly_location, nrow => :total_incidents_response,
                                            :cars => mean => :requested_cars,
                                            :cars_dispatched => mean => :dispatched_cars,
                                            :dispatch_time_first => mean => :dispatch_time_first,
                                            :dispatch_time_all => mean => :dispatch_time_all,
                                            :response_time_first => mean => :response_time_first,
                                            :response_time_all => mean => :response_time_all,
                                            :driving_time_first => mean => :driving_time_first,
                                            :cars_missing => mean => :cars_missing,
                                            :cars_location_responsible => mean => :cars_location_responsible)
    weekly_location = sort!(weekly_location, [:weekhour, :location_responsible])
    weekly_location[!,:exchange_ratio] = 1 .- weekly_location[!,:cars_location_responsible] ./
                                                weekly_location[!,:dispatched_cars]
    weekly_location.location_responsible = convert.(Int64, weekly_location.location_responsible)

# 2. group the results according to the priority and weekhour
    weekly_priority = groupby(select!(copy(incd), Not([:location_dispatched_first])), [:weekhour, :priority])
    weekly_priority = combine(weekly_priority,
                        nrow => :total_incidents_response,
                        :cars => mean => :requested_cars,
                        :cars_dispatched => mean => :dispatched_cars,
                        :dispatch_time_first => mean => :dispatch_time_first,
                        :dispatch_time_all => mean => :dispatch_time_all,
                        :response_time_first => mean => :response_time_first,
                        :response_time_all => mean => :response_time_all,
                        :driving_time_first => mean => :driving_time_first,
                        :cars_missing => mean => :cars_missing,
                        :cars_location_responsible => mean => :cars_location_responsible)
    weekly_priority = sort!(weekly_priority, [:weekhour, :priority])
    weekly_priority[!,:exchange_ratio] = 1 .- weekly_priority[!,:cars_location_responsible] ./
                                                weekly_priority[!,:dispatched_cars]
    weekly_priority[:,:priority] = Int.(weekly_priority[:,:priority])
    weekly_priority.priority = convert.(Int64, weekly_priority.priority)

# 3. group the results for a single output row
    main_results = combine(select!(copy(incd), Not([:location_dispatched_first])),
                                nrow => :total_incidents,
                                :cars => mean => :requested_cars,
                                :cars_dispatched => mean => :dispatched_cars,
                                :dispatch_time_first => mean => :dispatch_time_first,
                                :response_time_first => mean => :response_time_first,
                                :driving_time_first  => mean => :driving_time_first,
                                :dispatch_time_all   => mean => :dispatch_time_all,
                                :response_time_all   => mean => :response_time_all,
                                :cars_missing        => mean => :cars_missing,
                                :cars_location_responsible => mean => :cars_location_responsible)
    main_results[!,:incidents_unfulfilled]      .= size(filter(:cars_dispatched => isequal(0), incd),1)
    main_results[!,:total_incidents_fulfilled]   = size(incd,1) .- main_results[!,:incidents_unfulfilled]
    main_results[!,:cars_total_dispatched]       = main_results[!,:dispatched_cars] .* main_results[!,:total_incidents]
    main_results[!,:exchange_ratio]              = 1 .- main_results[!,:cars_location_responsible] ./ main_results[!,:dispatched_cars]
    main_results[!,:ratio_cars_undispatched]     = main_results[!,:cars_missing] ./ main_results[!,:dispatched_cars]

# 5. create a DataFrame that holds the ressource flow during the simulation
    capacity_status, mcs = ressource_status(ressource_flow::Array{Int64,3}, incd::DataFrame, simulation_capacity::Array{Int64,2})
    main_results = hcat(main_results,mcs)

    return main_results, weekly_location, weekly_priority, capacity_status
end


text:
Approaching the Unit Circle

    x(t) = cospi(t);  y(t) = sinpi(t)
    (x(t), y(t)) circles as t ↦ 0..2

```julia
using CatmullRom, Plots

fx(t) = cospi(t);
fy(t) = sinpi(t);

xs = [fx(t) for t=range(0.0, stop=2.0, length=17)];
ys = [fy(t) for t=range(0.0, stop=2.0, length=17)];
points = collect(zip(xs, ys));

cxs,cys = catmullrom(points, 16);

plot(cxs,cys, linecolor=:grey, linewidth=4, size=(600,600), legend=nothing, xaxis=nothing, yaxis=nothing)
plot!(xs, ys, linecolor=:blue, linewidth=2, size=(600,600), legend=nothing, xaxis=nothing, yaxis=nothing)
```

ext: jl


text:
const ADSUPPORT = (:ForwardDiff, :Zygote)
const ADAVAILABLE = Dict{Module, Function}()

Hamiltonian(metric::AbstractMetric, ℓπ, m::Module) = ADAVAILABLE[m](metric, ℓπ)

function Hamiltonian(metric::AbstractMetric, ℓπ)
    available = collect(keys(ADAVAILABLE))
    if length(available) == 0
        support_list_str = join(ADSUPPORT, " or ")
        error("MethodError: no method matching Hamiltonian(metric::AbstractMetric, ℓπ) because no backend is loaded. Please load an AD package ($support_list_str) first.")
    elseif length(available) > 1
        available_list_str = join(keys(ADAVAILABLE), " and ")
        constructor_list_str = join(map(m -> "Hamiltonian(metric, ℓπ, $m)", available), "\n  ")
        error("MethodError: Hamiltonian(metric::AbstractMetric, ℓπ) is ambiguous because multiple AD pakcages are available ($available_list_str). Please use AD explictly. Candidates:\n  $constructor_list_str")
    else
        return Hamiltonian(metric, ℓπ, first(available))
    end
end

### ForwardDiff

@require ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" begin

import .ForwardDiff, .ForwardDiff.DiffResults

function ∂ℓπ∂θ_forwarddiff(ℓπ, θ::AbstractVector)
    res = DiffResults.GradientResult(θ)
    ForwardDiff.gradient!(res, ℓπ, θ)
    return DiffResults.value(res), DiffResults.gradient(res)
end

# Implementation 1
function ∂ℓπ∂θ_forwarddiff(ℓπ, θ::AbstractMatrix)
    jacob = similar(θ)
    res = DiffResults.JacobianResult(similar(θ, size(θ, 2)), jacob)
    ForwardDiff.jacobian!(res, ℓπ, θ)
    jacob_full = DiffResults.jacobian(res)

    d, n = size(jacob)
    for i in 1:n
        jacob[:,i] = jacob_full[i,1+(i-1)*d:i*d]
    end
    return DiffResults.value(res), jacob
end

# Implementation 2
# function ∂ℓπ∂θ_forwarddiff(ℓπ, θ::AbstractMatrix)
#     local densities
#     f(x) = (densities = ℓπ(x); sum(densities))
#     res = DiffResults.GradientResult(θ)
#     ForwardDiff.gradient!(res, f, θ)
#     return ForwardDiff.value.(densities), DiffResults.gradient(res)
# end

# Implementation 3
# function ∂ℓπ∂θ_forwarddiff(ℓπ, θ::AbstractMatrix)
#     v = similar(θ, size(θ, 2))
#     g = similar(θ)
#     for i in 1:size(θ, 2)
#         res = GradientResult(θ[:,i])
#         gradient!(res, ℓπ, θ[:,i])
#         v[i] = value(res)
#         g[:,i] = gradient(res)
#     end
#     return v, g
# end

function ForwardDiffHamiltonian(metric::AbstractMetric, ℓπ)
    ∂ℓπ∂θ(θ::AbstractVecOrMat) = ∂ℓπ∂θ_forwarddiff(ℓπ, θ)
    return Hamiltonian(metric, ℓπ, ∂ℓπ∂θ)
end

ADAVAILABLE[ForwardDiff] = ForwardDiffHamiltonian

end # @require

### Zygote

@require Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" begin

import .Zygote

function ∂ℓπ∂θ_zygote(ℓπ, θ::AbstractVector)
    res, back = Zygote.pullback(ℓπ, θ)
    return res, first(back(Zygote.sensitivity(res)))
end

function ∂ℓπ∂θ_zygote(ℓπ, θ::AbstractMatrix)
    res, back = Zygote.pullback(ℓπ, θ)
    return res, first(back(ones(eltype(res), size(res))))
end

function ZygoteADHamiltonian(metric::AbstractMetric, ℓπ)
    ∂ℓπ∂θ(θ::AbstractVecOrMat) = ∂ℓπ∂θ_zygote(ℓπ, θ)
    return Hamiltonian(metric, ℓπ, ∂ℓπ∂θ)
end

ADAVAILABLE[Zygote] = ZygoteADHamiltonian

# Zygote.@adjoint

end # @require


text:
function selfie(x::Integer)
	ds = reverse(digits(x))
	if sum(ds) != length(ds) return false end
	for (i, d) in enumerate(ds)
		if d != sum(ds .== i - 1) return false end
	end
	return true
end

@show selfie(2020)
@show selfie(2021)

selfies(x) = for i in 1:x selfie(i) && println(i) end
@time selfies(4000000)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# define constants
using PyPlot

const g = 10.
const H = 10.
const N = 500
const Nt = 3000
const cfl = 0.9
const rho = 1.
const ν = 0.006
const F0 = 500.
const Lx = 1.
const dx = Lx/N
const dxinv = 1./dx
const cph = sqrt(g*H)
const dt = cfl * dx / cph
const RKα = [1/6.,1/3.,1/3.,1/6.]
const RKβ = [0.5,0.5,1.]

# grid
const x_h = dx/2:dx:Lx
const x_u = dx:dx:Lx        # no u-point at x=0 but at x=L (periodicity)

F(x,t) = -F0*sin.(8π*x/Lx + 200*t).*sin.(π*x/Lx).^2/rho

# OPERATORS
function Gux(res,u)
    res[1] = dxinv*(u[1]-u[end])
    res[2:end] = dxinv*(u[2:end]-u[1:end-1])
    return res
end

function GTx(res,η)
    res[1:end-1] = dxinv*(η[2:end]-η[1:end-1])
    res[end] = dxinv*(η[1]-η[end])
    return res
end

function ITu(res,η)
    res[1:end-1] = .5*(η[1:end-1] + η[2:end])
    res[end] = .5*(η[1]+η[end])
    return res
end

function IuT(res,u)
    res[1] = .5*(u[1] + u[end])
    res[2:end] = .5*(u[1:end-1]+u[2:end])
    return res
end

# initial conditions
η0 = zeros(N)
u0 = zeros(N)

function rhs(du,dη,h_u,u_h,dudx,u,η,t)
    h_u = ITu(h_u,η+H)
    u_h = IuT(u_h,u.^2)
    dudx = Gux(dudx,u)

    du = -GTx(du,.5*u_h + g*η - ν*dudx)
    du += F(x_h,t)./h_u

    dη = -Gux(dη,u.*h_u)

    return du,dη
end

function time_integration(Nt,u,η)

    # pre-allocate memory
    u0,η0 = zeros(u),zeros(η)
    u1,η1 = zeros(u),zeros(η)
    du,dη = zeros(u),zeros(η)
    h_u = zeros(u)
    u_h = zeros(η)
    dudx = zeros(η)
    t = 0.

    # for output
    u_out = zeros(Nt+1,length(u))
    η_out = zeros(Nt+1,length(η))

    # store the initial conditions
    u_out[1,:] = u
    η_out[1,:] = η

    for i = 1:Nt
        u1[:] = u
        η1[:] = η

        for RKi = 1:4
            du,dη = rhs(du,dη,h_u,u_h,dudx,u1,η1,t)

            if RKi < 4 # RHS update for the next RK-step
                u1 = u + RKβ[RKi]*dt*du
                η1 = η + RKβ[RKi]*dt*dη
            end

            # Summing all the RHS on the go
            u0 += RKα[RKi]*dt*du
            η0 += RKα[RKi]*dt*dη

        end

        u[:] = u0
        η[:] = η0
        t += dt

        # # output
        u_out[i+1,:] = u
        η_out[i+1,:] = η

    end
    u_out,η_out
end

time_integration(1,u0,η0)
@time u,η = time_integration(Nt,u0,η0)

## PLOTTING
fig,ax = subplots()

l1, = ax[:plot](x_h,η[1,:])
ax[:set_ylim](-2,2)
ax[:set_xlabel]("x")
ax[:set_ylabel]("η")
tight_layout()

for it = 2:5:Nt+1
    pause(0.00001)
    l1[:set_data](x_h,η[it,:])
end

pause(2)
close(fig)

max_stars_repo_name: milankl/misc

max_issues_repo_name: milankl/misc

max_forks_repo_name: milankl/misc


text:
using LatPhysBase
using Test

@testset "LatPhysBase.jl" begin

    # include subfile for abstract type tests
    include("tests_abstract_types.jl")

    # include subfile for default labels
    include("tests_default_labels.jl")

    # include subfile for concrete type tests
    include("tests_concrete_types.jl")

    # include subfile for custom type printing tests
    include("tests_custom_printing.jl")

end


text:
using Statistics, StatsBase
using PyCall
xr = pyimport("xarray")
using Blosc
using Elefridge

path = "/Users/milan/cams"
allvars = ["no2","go3","so2","aermr04","aermr05","aermr06","ch4","co"]
vari = allvars[3]

filelist = filter(x->endswith(x,vari*".grib"),readdir(path))
X = xr.open_dataarray(joinpath(path,filelist[1]),engine="cfgrib").data

# volume and error norms (columns):
# vol, mean, median, 90% and max
# total 5

# compression methods (rows):
# LinQuant24, LogQuant16, Bitgroom16/18/20/24 + Blosc, RoundNearest16/18/20/24 + Blosc
# total 10

# storing array
EV = Array{Float32,2}(undef,10,5)

## LinQuant24
EV[1,1] = 32/24     # compression factor
Xc = Array{Float32}(LinQuant24Array(X))
de = abs.(log10.(X ./ Xc))
EV[1,2] = mean(de)
EV[1,3] = median(de)
EV[1,4] = percentile(vec(de),90)
EV[1,5] = maximum(de)

# LoqQuant16
EV[2,1] = 2
Xc = Array{Float32}(LogQuant16Array(X))
de = abs.(log10.(X ./ Xc))
EV[2,2] = mean(de)
EV[2,3] = median(de)
EV[2,4] = percentile(vec(de),90)
EV[2,5] = maximum(de)

## Bitgroom16 + Blosc
Xc = groom(X,7)
de = abs.(log10.(X ./ Xc))
EV[3,2] = mean(de)
EV[3,3] = median(de)
EV[3,4] = percentile(vec(de),90)
EV[3,5] = maximum(de)
Xcc = compress(Xc,level=9)
EV[3,1] = sizeof(Xc)/sizeof(Xcc)

# Bitgroom18 + Blosc
Xc = groom(X,9)
de = abs.(log10.(X ./ Xc))
EV[4,2] = mean(de)
EV[4,3] = median(de)
EV[4,4] = percentile(vec(de),90)
EV[4,5] = maximum(de)
Xcc = compress(Xc,level=9)
EV[4,1] = sizeof(Xc)/sizeof(Xcc)

# Bitgroom20 + Blosc
Xc = groom(X,11)
de = abs.(log10.(X ./ Xc))
EV[5,2] = mean(de)
EV[5,3] = median(de)
EV[5,4] = percentile(vec(de),90)
EV[5,5] = maximum(de)
Xcc = compress(Xc,level=9)
EV[5,1] = sizeof(Xc)/sizeof(Xcc)

# Bitgroom24 + Blosc
Xc = groom(X,15)
de = abs.(log10.(X ./ Xc))
EV[6,2] = mean(de)
EV[6,3] = median(de)
EV[6,4] = percentile(vec(de),90)
EV[6,5] = maximum(de)
Xcc = compress(Xc,level=9)
EV[6,1] = sizeof(Xc)/sizeof(Xcc)

# RoundNearest16 + Blosc
Xc = round(X,7)
de = abs.(log10.(X ./ Xc))
EV[7,2] = mean(de)
EV[7,3] = median(de)
EV[7,4] = percentile(vec(de),90)
EV[7,5] = maximum(de)
Xcc = compress(Xc,level=9)
EV[7,1] = sizeof(Xc)/sizeof(Xcc)

# RoundNearest18 + Blosc
Xc = round(X,9)
de = abs.(log10.(X ./ Xc))
EV[8,2] = mean(de)
EV[8,3] = median(de)
EV[8,4] = percentile(vec(de),90)
EV[8,5] = maximum(de)
Xcc = compress(Xc,level=9)
EV[8,1] = sizeof(Xc)/sizeof(Xcc)

# RoundNearest20 + Blosc
Xc = round(X,11)
de = abs.(log10.(X ./ Xc))
EV[9,2] = mean(de)
EV[9,3] = median(de)
EV[9,4] = percentile(vec(de),90)
EV[9,5] = maximum(de)
Xcc = compress(Xc,level=9)
EV[9,1] = sizeof(Xc)/sizeof(Xcc)

# RoundNearest24 + Blosc
Xc = round(X,15)
de = abs.(log10.(X ./ Xc))
EV[10,2] = mean(de)
EV[10,3] = median(de)
EV[10,4] = percentile(vec(de),90)
EV[10,5] = maximum(de)
Xcc = compress(Xc,level=9)
EV[10,1] = sizeof(Xc)/sizeof(Xcc)

# just blosc
Xcc = compress(X,level=9)
bloscc = sizeof(X)/sizeof(Xcc)

## theoretical
theosize = [32/i for i in 32:-1:5]
x0 = 1f0
theoerror = [log10(nextfloat(x0,2^i)) for i in 0:27]

## plotting
pygui(true)
n = size(EV)[1]

cm1 = ColorMap("viridis")
cm2 = ColorMap("inferno")

function RGB(cv::Tuple)
    r = UInt8(round(cv[1]*255))
    g = UInt8(round(cv[2]*255))
    b = UInt8(round(cv[3]*255))
    return "#"*repr(r)[3:4]*repr(g)[3:4]*repr(b)[3:4]
end

colours = ["0.5","0.3",RGB(cm2(0.9)),RGB(cm2(0.8)),RGB(cm2(0.6)),RGB(cm2(0.4)),
                    RGB(cm1(0.8)),RGB(cm1(0.6)),RGB(cm1(0.4)),RGB(cm1(0.2))]

names = ["LinQuant24","LogQuant16","Bitgroom16+Blosc","Bitgroom18+Blosc",
            "Bitgroom20+Blosc","Bitgroom24+Blosc","RoundNearest16+Blosc",
            "RoundNearest18+Blosc","RoundNearest20+Blosc","RoundNearest24+Blosc"]

markers = ["s","s","d","d","d","d","o","o","o","o"]

fig,ax = subplots(1,1,figsize=(8,6))
ax.set_xscale("log")

ax.plot(theoerror,theosize,"k")
ax.plot(eps(1f0)/2e0,1,"w",markersize=15,marker="*",lw=0,markeredgewidth=1.4,
            markeredgecolor="k",alpha=1,label="Float32",zorder=10)
ax.plot(eps(1f0)/2e0,bloscc,"w",markersize=10,marker="v",lw=0,markeredgewidth=1.4,
            markeredgecolor="k",alpha=0.7,label="Float32+Blosc",zorder=10)

for i in [1,2,7,8,9,10,3,4,5,6]
    # ax.plot(EV[i,[3,5]],[EV[i,1],EV[i,1]],"k",lw=5.5,zorder=-1)
    ax.plot(EV[i,3:4],[EV[i,1],EV[i,1]],colours[i],
            marker=markers[i],markevery=2,
            markeredgecolor="k",markersize=10,
            lw=5,label=names[i])
    ax.plot(EV[i,4:5],[EV[i,1],EV[i,1]],colours[i],lw=2)
end

ax.set_xlim(eps(1f0)/3,1)
ax.set_ylim(0.6,3.5)
ax.set_xlabel("decimal error")
ax.set_ylabel("compression factor")
ax.legend(loc=4,ncol=3,fontsize=9)
ax.set_title(L"NO$_2$ compression",loc="left")

tight_layout()


text:
@testset "LLVM IR" begin

@testset "JuliaLang/julia#21121" begin
    function foobar()
        weight_matrix = @cuStaticSharedMem(Float32, (16, 16))
        sync_threads()
        weight_matrix[1, 16] *= 2
        sync_threads()
    end

    ir = sprint(io->CUDA.code_llvm(io, foobar, Tuple{}))
    @test !occursin("inttoptr", ir)
end

@testset "CUDA.jl#553" begin
    function kernel(ptr)
       unsafe_store!(ptr, CUDA.fma(unsafe_load(ptr), unsafe_load(ptr,2), unsafe_load(ptr,3)))
       return
    end

    ir = sprint(io->CUDA.code_llvm(io, kernel, Tuple{Ptr{Float32}}))
    @test !occursin("@__nv_fmaf", ir)
end

@testset "assume" begin
    foo(i) = cld(42, i)
    ir = sprint(io->CUDA.code_llvm(io, foo, Tuple{Int}))
    @test occursin("@gpu_report_exception", ir)

    bar(i) = (CUDA.assume(i > 0); cld(42, i))
    ir = sprint(io->CUDA.code_llvm(io, bar, Tuple{Int}))
    @test !occursin("gpu_report_exception", ir)
end

@testset "stripping invariant.load" begin
    function kernel(ptr, x)
        i = CUDA.threadIdx_x()
        @inbounds ptr[] = x[i]
        return
    end

    arr = CuArray(zeros(Float64))

    @cuda kernel(arr, (1., 2., ))
    @test Array(arr)[] == 1.
end

@testset "stripping const TBAA" begin
    # this one is particularly nasty because it occurs in a nested function

    _a = rand(Int, 2, 1)
    b = ((1,9999),(1,9999))

    out = CuArray(zeros(Int, 2,1))
    a = Tuple(_a)

    function kernel(out, a, b)
        i = threadIdx().x
        blockIdx().x
        @inbounds out[i,1] = a[i] + b[i][1]
        return
    end

    @cuda threads=2 kernel(out, a, b)
    @test Array(out) == (_a .+ 1)
end

@testset "ptxas-compatible control flow" begin
    @noinline function throw_some()
        throw(42)
        return
    end

    @inbounds function kernel(input, output, n)
        i = threadIdx().x

        temp = @cuStaticSharedMem(Int, 1)
        if i == 1
            1 <= n || throw_some()
            temp[1] = input
        end
        sync_threads()

        1 <= n || throw_some()
        unsafe_store!(output, temp[1], i)

        return
    end

    function gpu(input)
        output = CuArray(zeros(eltype(input), 2))
        ptr = pointer(output)
        ptr = reinterpret(Ptr{eltype(input)}, ptr)

        @cuda threads=2 kernel(input, ptr, 99)

        return Array(output)
    end

    function cpu(input)
        output = zeros(eltype(input), 2)

        for j in 1:2
            @inbounds output[j] = input
        end

        return output
    end

    input = rand(1:100)
    @test cpu(input) == gpu(input)
end

end

############################################################################################

@testset "PTX" begin

@testset "local memory stores due to byval" begin
    # JuliaGPU/GPUCompiler.jl#92
    function kernel(y1, y2)
        y = threadIdx().x == 1 ? y1 : y2
        @inbounds y[] = 0
        return
    end

    asm = sprint(io->CUDA.code_ptx(io, kernel, NTuple{2,CuDeviceArray{Float32,1,AS.Global}}))
    @test !occursin(".local", asm)
end

end

############################################################################################

@testset "SASS" begin

@testset "basic reflection" begin
    valid_kernel() = return
    invalid_kernel() = 1

    @not_if_memcheck @test CUDA.code_sass(devnull, valid_kernel, Tuple{}) == nothing
    @not_if_memcheck @test_throws CUDA.KernelError CUDA.code_sass(devnull, invalid_kernel, Tuple{})
end

@testset "function name mangling" begin
    @eval @noinline $(Symbol("dummy_^"))(x) = x

    @eval kernel_341(ptr) = (@inbounds unsafe_store!(ptr, $(Symbol("dummy_^"))(unsafe_load(ptr))); nothing)

    @not_if_memcheck CUDA.code_sass(devnull, kernel_341, Tuple{Ptr{Int}})
end

@testset "device runtime" begin
    kernel() = (CUDA.cudaGetLastError(); return)

    @not_if_memcheck CUDA.code_sass(devnull, kernel, Tuple{})
end

end


text:
# build the mock dynamic library
run(`make`)

# put the mock MultiNest library on DL_LOAD_PATH
push!(DL_LOAD_PATH, dirname(@__FILE__()))

using MultiNest
using Base.Test

# write your own tests here
@test 1 == 1

# done with tests, clean up
run(`make clean`)


text:
using Flux: Dense

in_channel = 3
out_channel = 5
N = 4
adj = [0. 1. 0. 1.;
       1. 0. 1. 0.;
       0. 1. 0. 1.;
       1. 0. 1. 0.]

@testset "layer" begin
    @testset "GCNConv" begin
        X = rand(in_channel, N)
        @testset "layer with graph" begin
            gc = GCNConv(adj, in_channel=>out_channel)
            @test size(gc.weight) == (out_channel, in_channel)
            @test size(gc.bias) == (out_channel,)
            @test graph(gc.fg) === adj

            Y = gc(X)
            @test size(Y) == (out_channel, N)
        end

        @testset "layer without graph" begin
            gc = GCNConv(in_channel=>out_channel, cache=false)
            @test size(gc.weight) == (out_channel, in_channel)
            @test size(gc.bias) == (out_channel,)
            @test !has_graph(gc.fg)

            fg = FeaturedGraph(adj, X)
            fg_ = gc(fg)
            @test size(node_feature(fg_)) == (out_channel, N)
            @test_throws AssertionError gc(X)
        end
    end

    @testset "ChebConv" begin
        k = 6
        X = rand(in_channel, N)
        @testset "layer with graph" begin
            cc = ChebConv(adj, in_channel=>out_channel, k)
            @test size(cc.weight) == (out_channel, in_channel, k)
            @test size(cc.bias) == (out_channel,)
            @test graph(cc.fg) === adj
            @test cc.k == k
            @test cc.in_channel == in_channel
            @test cc.out_channel == out_channel

            Y = cc(X)
            @test size(Y) == (out_channel, N)
        end

        @testset "layer without graph" begin
            cc = ChebConv(in_channel=>out_channel, k, cache=false)
            @test size(cc.weight) == (out_channel, in_channel, k)
            @test size(cc.bias) == (out_channel,)
            @test !has_graph(cc.fg)
            @test cc.k == k
            @test cc.in_channel == in_channel
            @test cc.out_channel == out_channel

            fg = FeaturedGraph(adj, X)
            fg_ = cc(fg)
            @test size(node_feature(fg_)) == (out_channel, N)
            @test_throws AssertionError cc(X)
        end
    end

    @testset "GraphConv" begin
        @testset "layer with graph" begin
            gc = GraphConv(adj, in_channel=>out_channel)
            @test graph(gc.fg) == [[2,4], [1,3], [2,4], [1,3]]
            @test size(gc.weight1) == (out_channel, in_channel)
            @test size(gc.weight2) == (out_channel, in_channel)
            @test size(gc.bias) == (out_channel,)

            X = rand(in_channel, N)
            Y = gc(X)
            @test size(Y) == (out_channel, N)
        end

        @testset "layer without graph" begin
            gc = GraphConv(in_channel=>out_channel)
            @test size(gc.weight1) == (out_channel, in_channel)
            @test size(gc.weight2) == (out_channel, in_channel)
            @test size(gc.bias) == (out_channel,)

            X = rand(in_channel, N)
            fg = FeaturedGraph(adj, X)
            fg_ = gc(fg)
            @test size(node_feature(fg_)) == (out_channel, N)
            @test_throws AssertionError gc(X)
        end
    end

    @testset "GATConv" begin
        X = rand(in_channel, N)

        @testset "layer with graph" begin
            @testset "concat=true" begin
                for heads = [1, 6]
                    gat = GATConv(adj, in_channel=>out_channel, heads=heads, concat=true)
                    @test graph(gat.fg) == [[2,4], [1,3], [2,4], [1,3]]
                    @test size(gat.weight) == (out_channel * heads, in_channel)
                    @test size(gat.bias) == (out_channel * heads,)
                    @test size(gat.a) == (2*out_channel, heads, 1)

                    Y = gat(X)
                    @test size(Y) == (out_channel * heads, N)
                end
            end

            @testset "concat=false" begin
                for heads = [1, 6]
                    gat = GATConv(adj, in_channel=>out_channel, heads=heads, concat=false)
                    @test graph(gat.fg) == [[2,4], [1,3], [2,4], [1,3]]
                    @test size(gat.weight) == (out_channel * heads, in_channel)
                    @test size(gat.bias) == (out_channel * heads,)
                    @test size(gat.a) == (2*out_channel, heads, 1)

                    Y = gat(X)
                    @test size(Y) == (out_channel * heads, 1)
                end
            end
        end

        @testset "layer without graph" begin
            fg = FeaturedGraph(adj, X)

            @testset "concat=true" begin
                for heads = [1, 6]
                    gat = GATConv(in_channel=>out_channel, heads=heads, concat=true)
                    @test size(gat.weight) == (out_channel * heads, in_channel)
                    @test size(gat.bias) == (out_channel * heads,)
                    @test size(gat.a) == (2*out_channel, heads, 1)

                    fg_ = gat(fg)
                    Y = node_feature(fg_)
                    @test size(Y) == (out_channel * heads, N)
                    @test_throws AssertionError gat(X)
                end
            end

            @testset "concat=false" begin
                for heads = [1, 6]
                    gat = GATConv(in_channel=>out_channel, heads=heads, concat=false)
                    @test size(gat.weight) == (out_channel * heads, in_channel)
                    @test size(gat.bias) == (out_channel * heads,)
                    @test size(gat.a) == (2*out_channel, heads, 1)

                    fg_ = gat(fg)
                    Y = node_feature(fg_)
                    @test size(Y) == (out_channel * heads, 1)
                    @test_throws AssertionError gat(X)
                end
            end
        end
    end

    @testset "GatedGraphConv" begin
        num_layers = 3
        @testset "layer with graph" begin
            ggc = GatedGraphConv(adj, out_channel, num_layers)
            @test graph(ggc.fg) == [[2,4], [1,3], [2,4], [1,3]]
            @test size(ggc.weight) == (out_channel, out_channel, num_layers)

            X = rand(in_channel, N)
            Y = ggc(X)
            @test size(Y) == (out_channel, N)
        end

        @testset "layer without graph" begin
            ggc = GatedGraphConv(out_channel, num_layers)
            @test size(ggc.weight) == (out_channel, out_channel, num_layers)

            X = rand(in_channel, N)
            fg = FeaturedGraph(adj, X)
            fg_ = ggc(fg)
            @test size(node_feature(fg_)) == (out_channel, N)
            @test_throws AssertionError ggc(X)
        end
    end

    @testset "EdgeConv" begin
        @testset "layer with graph" begin
            ec = EdgeConv(adj, Dense(2*in_channel, out_channel))
            @test graph(ec.fg) == [[2,4], [1,3], [2,4], [1,3]]

            X = rand(in_channel, N)
            Y = ec(X)
            @test size(Y) == (out_channel, N)
        end

        @testset "layer without graph" begin
            ec = EdgeConv(Dense(2*in_channel, out_channel))

            X = rand(in_channel, N)
            fg = FeaturedGraph(adj, X)
            fg_ = ec(fg)
            @test size(node_feature(fg_)) == (out_channel, N)
            @test_throws AssertionError ec(X)
        end
    end
end


text:
##using ITensors

## ######################################################################
## # XXX: Is this needed anywhere?
## # ITensors
## # PR: commontags with no inputs
## #
##
## import ITensors: commontags
##
## commontags() = TagSet()

## ######################################################################
## # XXX: Is this needed anywhere?
## # ITensors
## # PR: improve `rrule` for `setinds`.
## #
##
## using ChainRulesCore
## using ITensors.NDTensors
##
## using ITensors.NDTensors: setinds
##
## import ChainRulesCore: rrule
##
## function rrule(::typeof(setinds), x::Tensor, a...)
##   y = setinds(x, a...)
##   function setinds_pullback(ȳ)
##     x̄ = setinds(ȳ, inds(x))
##     return (NoTangent(), x̄, map(Returns(NoTangent()), a))
##   end
##   return y, setinds_pullback
## end

## ######################################################################
## # XXX: Is this needed anywhere?
## # ITensors
## # PR: define `ITensor(::ITensor)` and related functions.
## #
##
## using ChainRulesCore
## using ITensors.NDTensors
##
## import ITensors: itensor, ITensor, tensor, Tensor
##
## itensor(T::ITensor) = T
## ITensor(T::ITensor) = copy(T)
## tensor(T::Tensor) = T
## Tensor(T::Tensor) = copy(T)
##
## function rrule(::typeof(itensor), x::ITensor)
##   y = x
##   function itensor_pullback(ȳ)
##     x̄ = ȳ
##     return (NoTangent(), x̄)
##   end
##   return y, itensor_pullback
## end
##
## function rrule(::typeof(itensor), x::Tensor)
##   y = itensor(x)
##   function itensor_pullback(ȳ)
##     x̄ = tensor(ȳ)
##     return (NoTangent(), x̄)
##   end
##   return y, itensor_pullback
## end
##
## function rrule(::Type{ITensor}, x::Tensor)
##   y = ITensor(x)
##   function ITensor_pullback(ȳ)
##     x̄ = Tensor(ȳ)
##     return (NoTangent(), x̄)
##   end
##   return y, ITensor_pullback
## end
##
## function rrule(::Type{ITensor}, x::ITensor)
##   y = copy(x)
##   function ITensor_pullback(ȳ)
##     x̄ = copy(ȳ)
##     return (NoTangent(), x̄)
##   end
##   return y, ITensor_pullback
## end
##
## function rrule(::typeof(tensor), x::ITensor)
##   y = tensor(x)
##   function tensor_pullback(ȳ)
##     x̄ = itensor(ȳ)
##     return (NoTangent(), x̄)
##   end
##   return y, tensor_pullback
## end
##
## function rrule(::typeof(tensor), x::Tensor)
##   y = x
##   function tensor_pullback(ȳ)
##     x̄ = ȳ
##     return (NoTangent(), x̄)
##   end
##   return y, tensor_pullback
## end
##
## function rrule(::Type{Tensor}, x::ITensor)
##   y = Tensor(x)
##   function Tensor_pullback(ȳ)
##     x̄ = ITensor(ȳ)
##     return (NoTangent(), x̄)
##   end
##   return y, Tensor_pullback
## end
##
## function rrule(::Type{Tensor}, x::Tensor)
##   y = copy(x)
##   function Tensor_pullback(ȳ)
##     x̄ = copy(ȳ)
##     return (NoTangent(), x̄)
##   end
##   return y, Tensor_pullback
## end

## ######################################################################
## # XXX: This is WIP to improve differentiability of priming/tagging MPS/MPO.
## # ITensors
## # PR: make `map` of `AbstractMPS` differentiable
## #
##
## using ITensors: AbstractMPS, data
##
## import ITensors: MPS, MPO
## import Base: map
## import ChainRulesCore: rrule
##
## function map(f::Function, M::AbstractMPS; set_limits::Bool=true)
##   new_ortho_lims = set_limits ? eachindex(M) : ortho_lims(M)
##   return typeof(M)(map(f, data(M)); ortho_lims=new_ortho_lims)
## end
##
## function rrule(::Type{MPS}, x::Vector{ITensor}, a1::Int, a2::Int)
##   y = MPS(x, a1, a2)
##   function MPS_rrule(ȳ)
##     x̄ = data(ȳ)
##     ā1 = NoTangent()
##     ā2 = NoTangent()
##     return (NoTangent(), x̄, ā1, ā2)
##   end
##   return y, MPS_rrule
## end
##
## function rrule(::Type{MPS}, x::Vector{ITensor}; kwargs...)
##   y = MPS(x; kwargs...)
##   function MPS_rrule(ȳ)
##     x̄ = ȳ.data
##     return (NoTangent(), x̄)
##   end
##   return y, MPS_rrule
## end
##
## function rrule(::Type{MPO}, x::Vector{ITensor}, a1::Int, a2::Int)
##   y = MPO(x, a1, a2)
##   function MPO_rrule(ȳ)
##     x̄ = data(ȳ)
##     ā1 = NoTangent()
##     ā2 = NoTangent()
##     return (NoTangent(), x̄, ā1, ā2)
##   end
##   return y, MPO_rrule
## end
##
## function rrule(::Type{MPO}, x::Vector{ITensor}; kwargs...)
##   y = MPO(x; kwargs...)
##   function MPO_rrule(ȳ)
##     x̄ = data(ȳ)
##     return (NoTangent(), x̄)
##   end
##   return y, MPO_rrule
## end
##
## function rrule(::typeof(getindex), x::MPS, I::Int)
##   @show I
##   y = getindex(x, I)
##   function getindex_pullback(ȳ)
##     x̄ = MPS([ITensor(inds(x[i])) for i in eachindex(x)])
##     x̄[I] = ȳ
##     return (NoTangent(), x̄, NoTangent())
##   end
##   return y, getindex_pullback
## end
##
## # XXX: This may be a bad definition! Maybe should be `x1 + x1`?
## Zygote.accum(x1::MPS, x2::MPS) = x1 .+ x2


text:
# # Minimizing the sum of the k-largest λ
#
# We show how to find the sum of absolute value of the `k` largest eigenvalues of a symmetric matrix $A \in \mathbb{S}^n$. This problem can be solved as a semidefinite program. The primal and dual forms are stated in *Alizadeh* \[1\]:
# $$
# \begin{array}{llll} \text{maximize} &   \text{Tr}(AY) - \text{Tr}(AW) &    \text{minimize} &  kz + Tr(U) + Tr(V)            \\
# \text{subject to} &  \text{Tr}(Y + W) = k                             &   \text{subject to} &   zI + V - A \succeq 0            \\
#                   &  0 \preceq  Y \preceq I                           &  &                    zI + U + A \succeq 0             \\
#                   &  0 \preceq  W \preceq I                           & &                     U, V \succeq 0,
# \end{array}
# $$
# where $Y, W$ are the variables of the primal and $U, V$ are the variables of the dual problem.
#-
using LinearAlgebra, JuMP, COSMO, Random
rng = Random.MersenneTwister(212)

n = 10
A = 5 .* randn(rng, 10, 10)
A = Symmetric(A, :U)

# We are interested in minimizing the sum of absolute values of the `k=3` largest eigenvalues. Let's formulate the problem in `JuMP` with `COSMO` as the backend solver:

k = 3
model = JuMP.Model(optimizer_with_attributes(COSMO.Optimizer, "verbose" => true));
@variable(model, Y[1:n, 1:n], PSD);
@variable(model, W[1:n, 1:n], PSD);

@objective(model, Max, tr(A * Y) - tr(A * W));
@constraint(model, tr(Y + W) == k);
@constraint(model, Symmetric(I - Y) in PSDCone());
@constraint(model, Symmetric(I - W) in PSDCone());
status = JuMP.optimize!(model)

#-
opt_objective = JuMP.objective_value(model)

# Now, we can check the solution by computing the sum of the absolute value of the 3-largest eigenvalues:
k_λ_abs = sum(sort(abs.(eigen(A).values), rev = true)[1:k])

# ### Solve the dual
#
# Alternatively, we can solve the dual problem:
#-
model = JuMP.Model(with_optimizer(COSMO.Optimizer, verbose=true));
@variable(model, V[1:n, 1:n], PSD);
@variable(model, U[1:n, 1:n], PSD);
@variable(model, z);

@objective(model, Min, k * z + tr(V) + tr(U));
@constraint(model, Symmetric(z .* diagm(0 => ones(n)) + V - A) in PSDCone());
@constraint(model, Symmetric(z .* diagm(0 => ones(n)) + U + A) in PSDCone());
status = JuMP.optimize!(model)

#-
opt_objective = JuMP.objective_value(model)

# This gives the same result.
# ## Problem with A as variable
#
# Above problems are mostly helpful for illustrative purpose. It is obviously easier to find the sum of the k-largest eigenvalues by simply computing the eigenvalues of $A$. However, above results become useful if finding $A$ itself is part of the problem. For example, assume we want to find a valid matrix $A$ under the constraints: $C\, \text{vec}(A) = b$ with the minimum sum of absolute values of the k-largest eigenvalues. We can then solve the equivalent problem:
# $$
# \begin{array}{ll} \text{minimize} &  kz + Tr(U) + Tr(V)     \\
#  \text{subject to} &   C \text{vec}(A) = b \\
#                    & zI + V - A \succeq 0            \\
#                    &    zI + U + A \succeq 0             \\
#                    &   U, V \succeq 0.
# \end{array}
# $$
#
# ## References
# [1] Alizadeh - Interior point methods in semidefinite programming with applications to combinatorial optimization (1995)


text:
abstract type GraphvizPoperties end

mutable struct Property{T}
    key::String
    value::T
end

const Properties = Vector{Property}

# return val of attribute:
function val(attributes::Properties, attribute::String)
    if !isempty(attributes)
        for a in attributes
            if a.key == attribute
                return a.value
            end
        end
    end
    return []
end

# return Tuple (bool, pos). bool=true if key exist in Attributes
import Base.haskey
function haskey(attributes::Properties, key::String)
    if !isempty(attributes)
        for i = 1:length(attributes)
            if attributes[i].key == key
                return true, i
            end
        end
    end
    return false, 0
end

# set val to attributeDict
set!(attributes::Properties, prop::Property; override=true) = set!(attributes, prop.key, prop.value; override)

function set!(attributes::Properties, key::String, value; override=true)
    key_exist, idx = haskey(attributes, key)
    if key_exist & (override == true)
        attributes[idx].value = check_value(value)
    else
        push!(attributes, Property(key, check_value(value)))
    end
end

# remove attribute of attributes
rm!(attributes::Properties, prop::Property) = rm!(attributes, prop.key)

function rm!(attributes::Properties, key::String)
    key_exist, idx = haskey(attributes, key)
    (key_exist) ? deleteat!(attributes, idx) : nothing
end

max_stars_repo_name: tragisch/PlotGraphviz

max_issues_repo_name: tragisch/PlotGraphviz


text:
# https://adventofcode.com/2021/day/13
using AdventOfCode

input = readlines("2021/data/day_13.txt")
dots=parse.(input[begin:begin+findfirst(==(""),input)-2])
folds=replace.(input[begin+findfirst(==(""),input):end],"fold along "=>":")
folds=replace.(folds,"="=>"=>").|>parse
axes=Dict(:x=>1,:y=>2)
input=(dots=dots,folds=folds,axes=axes)

function makepaper(input)
    papersize = ((maximum.(input.dots)|>maximum) +1)*2
    paper=zeros(Bool,papersize,papersize)
    paper[CartesianIndex.(input.dots).+Ref(CartesianIndex(1,1))].=1
    paper
end

function fold!(paper, foldselem)
    axis,at = foldselem
    if axis==:x
        paper[1:at,:] .= paper[1:at,:].| paper[1+at.+(at:-1:1),:]
        paper[at+1:end,:].=0
    end
    if axis==:y
        @views paper[:,1:at] .= paper[:,1:at].| paper[:,1+at.+(at:-1:1)]
        paper[:,at+1:end].=0
    end
end
function part_1(input)
    paper=makepaper(input)
    fold!(paper,input.folds[1])
    count(paper)
end

function part_2(input)
    paper=makepaper(input)
    fold!.(Ref(paper),input.folds)
    println.(eachrow((Int.(paper)')[1:6,1:40]).|>join)
    nothing
end


text:
const TIMER_OUTPUTS = true
const being_timed = Ref{Bool}(false)

if TIMER_OUTPUTS
    using TimerOutputs

    macro timer(name, expr)
        :(if being_timed[]
              @timeit $(esc(name)) $(esc(expr))
          else
              $(esc(expr))
          end)
    end

    macro iftimer(expr)
        esc(expr)
    end

else
    macro timer(name, expr)
        esc(expr)
    end

    macro iftimer(expr)
    end
end

using Base: ImmutableDict

# Linked List interface
@inline assoc(d::ImmutableDict, k, v) = ImmutableDict(d, k=>v)

struct LL{V}
    v::V
    i::Int
end

islist(::Any) = false
islist(l::Union{LL, Term, AbstractArray, Tuple}) = true

Base.empty(l::LL) = empty(l.v)
Base.isempty(l::LL) = l.i > length(l.v)

Base.length(l::LL) = length(l.v)-l.i+1
@inline car(l::LL) = l.v[l.i]
@inline cdr(l::LL) = isempty(l) ? empty(l) : LL(l.v, l.i+1)

Base.length(t::Term) = length(arguments(t)) + 1 # PIRACY
Base.isempty(t::Term) = false
@inline car(t::Term) = operation(t)
@inline cdr(t::Term) = arguments(t)

@inline car(v) = first(v)
@inline cdr(v) = isempty(v) ? empty(l) : LL(v, 2)

@inline take_n(ll::LL, n) = isempty(ll) || n == 0 ? empty(ll) : @views ll.v[ll.i:n+ll.i-1] # @views handles Tuple
@inline take_n(ll, n) = @views ll[1:n]

drop_n(ll, n) = n === 0 ? ll : drop_n(cdr(ll), n-1)

@inline drop_n(ll::Term, n) = drop_n(arguments(ll), n-1)
@inline drop_n(ll::AbstractArray, n) = drop_n(LL(ll, 1), n)
@inline drop_n(ll::LL, n) = LL(ll.v, ll.i+n)

pow(x,y) = y==0 ? 1 : y<0 ? inv(x)^(-y) : x^y
pow(x::Symbolic,y) = y==0 ? 1 : Base.:^(x,y)
pow(x, y::Symbolic) = Base.:^(x,y)
pow(x::Symbolic,y::Symbolic) = Base.:^(x,y)

# Simplification utilities
has_trig(x) = false
function has_trig(term::Term)
    fns = (sin, cos, tan, cot, sec, csc)
    op = operation(term)

    if Base.@nany 6 i->fns[i] === op
        return true
    else
        return any(has_trig, arguments(term))
    end
end

fold(x) = x
function fold(t::Term)
    tt = map(fold, arguments(t))
    if !any(x->x isa Symbolic, tt)
        # evaluate it
        return operation(t)(tt...)
    else
        return Term{symtype(t)}(operation(t), tt)
    end
end

### Predicates

sym_isa(::Type{T}) where {T} = @nospecialize(x) -> x isa T || symtype(x) <: T
is_operation(f) = @nospecialize(x) -> (x isa Term) && (operation(x) == f)

isliteral(::Type{T}) where {T} = x -> x isa T
isnumber(x) = isliteral(Number)(x)

_iszero(t) = false
_iszero(x::Number) = iszero(x)
_isone(t) = false
_isone(x::Number) = isone(x)

issortedₑ(args) = issorted(args, lt=<ₑ)

# are there nested ⋆ terms?
function isnotflat(⋆)
    function (args)
        for t in args
            if t isa Term && operation(t) === (⋆)
                return true
            end
        end
        return false
    end
end

function hasrepeats(x)
    length(x) <= 1 && return false
    for i=1:length(x)-1
        if isequal(x[i], x[i+1])
            return true
        end
    end
    return false
end

function merge_repeats(merge, xs)
    length(xs) <= 1 && return false
    merged = Any[]
    i=1

    while i<=length(xs)
        l = 1
        for j=i+1:length(xs)
            if isequal(xs[i], xs[j])
                l += 1
            else
                break
            end
        end
        if l > 1
            push!(merged, merge(xs[i], l))
        else
            push!(merged, xs[i])
        end
        i+=l
    end
    return merged
end

# Numbers to the back
function flatten_term(⋆, args)
    # flatten nested ⋆
    flattened_args = []
    for t in args
        if t isa Term && operation(t) === (⋆)
            append!(flattened_args, arguments(t))
        else
            push!(flattened_args, t)
        end
    end
    Term(⋆, flattened_args)
end

function sort_args(f, args)
    if length(args) < 2
        return Term(f, args)
    elseif length(args) == 2
        x, y = args
        return Term(f, x <ₑ y ? [x,y] : [y,x])
    end
    args = args isa Tuple ? [args...] : args
    Term(f, sort(args, lt=<ₑ))
end


text:
@testset "test Kaczmarz update" begin
    for T in [Float32,Float64,ComplexF32,ComplexF64]
        # set up
	M = 127
	N = 16

	A = rand(T,M,N)
	Aᵀ = transpose(A)
	b = zeros(T,M)
	β = rand(T)
	k = rand(1:N)
	# end set up

	RegularizedLeastSquares.kaczmarz_update!(Aᵀ,b,k,β)
	@test b ≈ β*conj(A[:,k])
    end

    for T in [Float32,Float64,ComplexF32,ComplexF64]
        # set up
        M = 127
        N = 16

        A = rand(T,N,M)
        b = zeros(T,M)
        β = rand(T)
        k = rand(1:N)
        # end set up

	RegularizedLeastSquares.kaczmarz_update!(A,b,k,β)
        @test b ≈ β*conj(A[k,:])
    end
end


text:
using POMCPOW
using Test

using POMDPs
using POMDPModels
using ParticleFilters
using POMDPTesting
using D3Trees
using BeliefUpdaters
using POMDPModelTools

solver = POMCPOWSolver()

pomdp = BabyPOMDP()

test_solver(solver, pomdp, updater=DiscreteUpdater(pomdp))
test_solver(solver, pomdp)

# make sure internal function is type stable
planner = solve(solver, pomdp)
b = initialstate_distribution(pomdp)
B = POMCPOW.belief_type(POMCPOW.POWNodeFilter, typeof(pomdp))
tree = POMCPOWTree{B,Bool,Bool,typeof(b)}(b, 2*planner.solver.tree_queries)
@inferred POMCPOW.simulate(planner, POMCPOW.POWTreeObsNode(tree, 1), true, 10)
# @code_warntype POMCPOW.simulate(planner, POMCPOW.POWTreeObsNode(tree, 1), true, 10)

pomdp = LightDark1D()
solver = POMCPOWSolver(default_action=485)
planner = solve(solver, pomdp)

b = ParticleCollection([LightDark1DState(-1, 0)])
@test @inferred(action(planner, b)) == 485

b = initialstate_distribution(pomdp)
@inferred action(planner, b)

a, info = action_info(planner, b)
# d3t = D3Tree(planner)
@test_throws KeyError d3t = D3Tree(info[:tree])

a, info = action_info(planner, b, tree_in_info=true)
# d3t = D3Tree(planner)
d3t = D3Tree(info[:tree])
# inchrome(d3t)

include("init_node_sr_belief_error.jl")


text:
@testset "structures.jl" begin
    include(joinpath("structures", "cluster_test.jl"))
    include(joinpath("structures", "dataitem_test.jl"))
    include(joinpath("structures", "eva_test.jl"))
    include(joinpath("structures", "fittedeva_test.jl"))
    include(joinpath("structures", "returnlevel_test.jl"))

end


text:
using MinFinVol
using Plots

# defining a mesh
Lx=1.0
Ly=1.0

nx=20
ny=20

xs=LinRange(0.0, Lx, nx)
ys=LinRange(0.0, Ly, ny)

# mesh generation
msh=FinVolMesh{2}()

ptinds=zeros(Int64, nx, ny)
for (i, x) in enumerate(xs)
    for (j, y) in enumerate(ys)
        isedge=(i==1 || i==nx || j==1 || j==ny)
        w=((isedge & i==1) ? 1.1 : 1.0)

        ptinds[i, j]=add_point!(msh, [x, y]; fixed=isedge, weight=w)
    end
end

left_face=Vector{Int64}()
lower_face=Vector{Int64}()
upper_face=Vector{Int64}()
right_face=Vector{Int64}()

volume_inds=zeros(Int64, nx-1, ny-1)
for i=1:(nx-1)
    for j=1:(ny-1)
        newfcs=[
            add_face!(msh, [ptinds[i, j], ptinds[i, j+1]]; univocal=true),
            add_face!(msh, [ptinds[i, j+1], ptinds[i+1, j+1]]; univocal=true),
            add_face!(msh, [ptinds[i+1, j+1], ptinds[i+1, j]]; univocal=true),
            add_face!(msh, [ptinds[i+1, j], ptinds[i, j]]; univocal=true)
        ]

        if i==1
            push!(left_face, newfcs[1])
        elseif i==(nx-1)
            push!(right_face, newfcs[3])
        end

        if j==1
            push!(lower_face, newfcs[4])
        elseif j==(ny-1)
            push!(upper_face, newfcs[2])
        end

        volume_inds[i, j]=add_volume!(msh, newfcs)
    end
end

process_geometry!(msh)
# mesh_smoothing!(msh)

Dirichlet_left=1.0
Dirichlet_right=0.0
Dirichlet_lower=0.0
# Dirichlet_upper=0.0
x_Neumann_upper=zeros(Float64, length(upper_face))
y_Neumann_upper=ones(Float64, length(upper_face))
Neumann_upper=[collect(t) for t in zip(x_Neumann_upper, y_Neumann_upper)]

# pre-allocations
_at_faces=zeros(Float64, nx*(ny-1)+ny*(nx-1))
# _at_volumes=zeros(Float64, (nx-1)*(ny-1))

_qx=zeros(Float64, size(_at_faces)...)
_qy=zeros(Float64, size(_at_faces)...)

function getres!(r, u)
    # variable values at faces
    face_average!(msh, u, _at_faces)

    # applying Dirichlet conditions
    _at_faces[left_face].=Dirichlet_left
    _at_faces[right_face].=Dirichlet_right
    _at_faces[lower_face].=Dirichlet_lower
    # _at_faces[upper_face].=Dirichlet_upper
    Neumann_conditions_impose!(msh, u, _at_faces, (upper_face, Neumann_upper))

    # obtain derivatives
    _qx.=_at_faces
    face_gradient!(msh, u, _qx, 1)

    _qy.=_at_faces
    face_gradient!(msh, u, _qy, 2)

    flux_integral!(msh, r, _qx, _qy)
end

u=zeros(Float64, (nx-1)*(ny-1))
r=zeros(Float64, size(u)...)

dt=1e-4

for i=1:2000
    getres!(r, u)

    u.-=(dt.*r)
end

plotlyjs()

x_plot=[msh.volume_centers[(nx-1)*(i-1)+j][1] for i=1:(nx-1), j=1:(ny-1)]
y_plot=[msh.volume_centers[(nx-1)*(i-1)+j][2] for i=1:(nx-1), j=1:(ny-1)]

p=surface(x_plot, y_plot, u[volume_inds])
display(p)


text:
using Base.Cartesian

function cudims(n::Integer)
  threads = min(n, 256)
  ceil(Int, n / threads), threads
end

cudims(a::AbstractArray) = cudims(length(a))

@inline ind2sub_(a::AbstractArray{T,0}, i) where T = ()
@inline ind2sub_(a, i) = Tuple(CartesianIndices(a)[i])

macro cuindex(A)
  quote
    A = $(esc(A))
    i = (blockIdx().x-1) * blockDim().x + threadIdx().x
    i > length(A) && return
    ind2sub_(A, i)
  end
end

@generated function nindex(i::T, ls::NTuple{N,T}) where {N,T}
  na = one(i)
  quote
    Base.@_inline_meta
    $(foldr((n, els) -> :(i ≤ ls[$n] ? ($n, i) : (i -= ls[$n]; $els)), :($na, $na), one(i):i(N)))
  end
end

@inline function catindex(dim, I::NTuple{N}, shapes) where N
  @inbounds x, i = nindex(I[dim], getindex.(shapes, dim))
  x, ntuple(n -> n == dim ? i : I[n], Val{N})
end

function growdims(dim, x)
  if ndims(x) >= dim
    x
  else
    reshape(x, size.((x,), 1:dim)...)
  end
end

function _cat(dim, dest, xs...)
  function kernel(dim, dest, xs)
    I = @cuindex dest
    @inbounds n, I′ = catindex(dim, Int.(I), size.(xs))
    @inbounds dest[I...] = xs[n][I′...]
    return
  end
  xs = growdims.(dim, xs)
  blk, thr = cudims(dest)
  @cuda blocks=blk threads=thr kernel(dim, dest, xs)
  return dest
end

function Base.cat_t(dims::Integer, T::Type, x::CuArray, xs::CuArray...)
  catdims = Base.dims2cat(dims)
  shape = Base.cat_shape(catdims, (), size.((x, xs...))...)
  dest = Base.cat_similar(x, T, shape)
  _cat(dims, dest, x, xs...)
end

Base.vcat(xs::CuArray...) = cat(xs..., dims=1)
Base.hcat(xs::CuArray...) = cat(xs..., dims=2)

"""
    @sync ex

Run expression `ex` and synchronize the GPU afterwards. This is a CPU-friendly
synchronization, i.e. it performs a blocking synchronization without increasing CPU load. As
such, this operation is preferred over implicit synchronization (e.g. when performing a
memory copy) for high-performance applications.

It is also useful for timing code that executes asynchronously.
"""
macro sync(ex)
    quote
        local e = CuEvent(CUDAdrv.EVENT_BLOCKING_SYNC | CUDAdrv.EVENT_DISABLE_TIMING)
        local ret = $(esc(ex))
        CUDAdrv.synchronize(e)
        ret
    end
end


text:
# Wrapper of Matrix that is callable as `A(v)`
struct MatrixMap{T}
  data::Matrix{T}
end
(A::MatrixMap)(v) = A.data * v
Base.:+(A::MatrixMap, B::MatrixMap) = MatrixMap(A.data + B.data)

max_stars_repo_name: mtfishman/ITensorParallel


text:
"""
doublet_matrix(panels_1, panels_2)

Create the matrix of doublet potential influence coefficients between pairs of panels_1 and panels_2.
"""
doublet_matrix(panels_1, panels_2) = [ doublet_influence(panel_j, panel_i) for panel_i in panels_1, panel_j in panels_2 ]

"""
    kutta_condition(panels)

Create the vector describing Morino's Kutta condition given Panel2Ds.
"""
kutta_condition(panels) = [ 1 zeros(length(panels) - 2)' -1 ]

"""
    wake_vector(wake_panel, panels)

Create the vector of doublet potential influence coefficients from the wake on the panels given the wake panel and the array of Panel2Ds.
"""
wake_vector(woke_panel :: AbstractPanel2D, panels) = doublet_influence.(Ref(woke_panel), panels)

"""
    influence_matrix(panels, wake_panel)

Assemble the Aerodynamic Influence Coefficient matrix consisting of the doublet matrix, wake vector, Kutta condition given Panel2Ds and the wake panel.
"""
influence_matrix(panels, woke_panel :: AbstractPanel2D) =
    [ doublet_matrix(panels, panels)  wake_vector(woke_panel, panels) ;
        kutta_condition(panels)                      1.               ]

"""
    source_matrix(panels_1, panels_2)

Create the matrix of source potential influence coefficients between pairs of `panels_1` and `panels_2`.
"""
source_matrix(panels_1, panels_2) = [ source_influence(panel_j, panel_i) for (panel_i, panel_j) in product(panels_1, panels_2) ]

"""
    source_strengths(panels, freestream)

Create the vector of source strengths for the Dirichlet boundary condition ``\\sigma = \\vec U_{\\infty} \\cdot \\hat{n}`` given Panel2Ds and a Uniform2D.
"""
source_strengths(panels, u) = dot.(Ref(u), panel_normal.(panels))

"""
    boundary_vector(panels, u)

Create the vector for the boundary condition of the problem given an array of Panel2Ds and velocity ``u``.
"""
boundary_vector(panels, u) = [ - source_matrix(panels, panels) * source_strengths(panels, u); 0 ]

boundary_vector(colpoints, u, r_te) = [ dot.(colpoints, Ref(u)); dot(u, r_te) ]

# boundary_vector(panels, u, r_te) = [ dot.(collocation_point.(panels), Ref(u)); dot(u, r_te) ]

function boundary_vector(panels :: Vector{<: AbstractPanel2D}, wakes :: Vector{<: AbstractPanel2D}, u)
    source_panels = [ panels; wakes ]
    [ - source_matrix(panels, source_panels) * source_strengths(source_panels, u); 0 ]
end

"""
    solve_strengths(panels, u, sources, bound)

Solve the system of equations ``[AIC][\\phi] = [\\vec{U} \\cdot \\hat{n}] - B[\\sigma]`` condition given the array of Panel2Ds, a velocity ``\\vec U``, a condition whether to disable source terms (``\\sigma = 0``), and an optional named bound for the length of the wake.
"""
function solve_strengths(panels, u, α, r_te, sources :: Bool; bound = 1e2)
    # Wake
    woke_panel  = wake_panel(panels, bound, α)
    woke_vector = wake_vector(woke_panel, panels)
    woke_matrix = [ -woke_vector zeros(length(panels), length(panels) -2) woke_vector ]

    # AIC
    AIC 	= doublet_matrix(panels, panels) + woke_matrix
    boco 	= dot.(collocation_point.(panels), Ref(u)) - woke_vector * dot(u, r_te)

    # AIC
    # AIC 	= influence_matrix(panels, woke_panel)
    # boco 	= boundary_vector(ifelse(sources, panels, collocation_point.(panels)), u, r_te) - [ woke_vector; 0 ] .* dot(u, r_te)

    AIC \ boco
end

"""
    tangential_velocities(panels, φs, u, sources :: Bool)

Compute the tangential velocities and panel distances given the array of `Panel2D`s, their associated doublet strengths ``\\phi``s, the velocity ``u``, and a condition whether to disable source terms (``\\sigma = 0``).
"""
function tangential_velocities(panels, φs, u, sources :: Bool)
    # Δrs   = midpair_map(panel_dist, panels)
    # Δφs   = -midpair_map(-, φs[1:end-1])

    Δrs   	 = @. panel_dist(panels[2:end], panels[1:end-1])
    Δφs   	 = @. φs[1:end-1] - φs[2:end]
    tangents = @. panel_tangent(panels[2:end])

    vels  = ifelse(sources, panel_velocity.(Δφs, Δrs, Ref(u), tangents), Δφs ./ Δrs)

    vels, Δrs
end

## WAKE VERSIONS
#==========================#

"""
    solve_strengths(panels, u, wakes, bound)

Solve the system of equations ``[AIC][\\phi] = [\\vec{U} \\cdot \\hat{n}] - B[\\sigma]`` condition given the array of Panel2Ds, a velocity ``\\vec U``, a vector of wake `Panel2D`s, and an optional named bound for the length of the wake.
"""
function solve_strengths(panels, u, α, wakes; bound = 1e2)
    AIC  = influence_matrix(panels, wake_panel(panels, bound, α))
    boco = boundary_vector(panels, wakes, u)

    AIC \ boco
end

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
import Glob

function proc_subset(
    lats,
    lons,
    out_dir;
    ref_station = nothing,
    ref_row = nothing,
    ref_col = nothing,
)
    intf = Glob.glob("*.int")[1]  # Dummy file
    i1 = MapImages.MapImage(intf)

    # demrsc = Sario.load("dem.rsc")
    cropped = i1[lats, lons]
    Sario.save(joinpath(out_dir, "dem.rsc"), cropped.demrsc)

    # Now get the range for cropping all the other files
    rowrange, colrange = MapImages.latlon_to_rowcol(i1.demrsc, lats, lons)

    @show rowrange, colrange, lats, lons
    println("New dem:")
    @show cropped.demrsc
    println("Ints:")
    @time InsarTimeseries.loop_over_files(
        x -> x[rowrange, colrange],
        ".",
        ".int",
        ".int",
        out_dir = out_dir,
        overwrite = false,
        do_permute = true,
    )
    println("ccs:")
    @time InsarTimeseries.loop_over_files(
        x -> x[rowrange, colrange],
        ".",
        ".cc",
        ".cc",
        out_dir = out_dir,
        overwrite = false,
        do_permute = true,
    )
    println("geo masks:")
    @time InsarTimeseries.loop_over_files(
        x -> x[rowrange, colrange],
        ".",
        ".geo.mask",
        ".geo.mask",
        out_dir = out_dir,
        overwrite = false,
        do_permute = true,
    )

    cd(out_dir)
    run(`insar process --step 8`)
    @time InsarTimeseries.prepare_stacks(
        ".",
        ref_station = ref_station,
        ref_row = ref_row,
        ref_col = ref_col,
        overwrite = false,
        zero_masked = false,
    )
end


text:
using CSV, DataFrames, Dates

df1 = CSV.read("/mnt/3fbf694d-d8e0-46c0-903d-69d994241206/mintsData/raw/001e06305a6c/2019/11/09/MINTS_001e06305a6c_AS7262_2019_11_09.csv", DataFrame)         # Reading ith CSV data file
df2 = CSV.read("/mnt/3fbf694d-d8e0-46c0-903d-69d994241206/mintsData/raw/001e06305a6c/2019/11/08/MINTS_001e06305a6c_AS7262_2019_11_08.csv", DataFrame)
#dfn = df[1:3,:]
#df.dateTime =  SubString.(string.(df.dateTime), 1, 19)
#df.dateTime = DateTime.(df.dateTime,"yyyy-mm-dd HH:MM:SS")
#DateTime.(df.dateTime[2],"yyyy-mm-dd HH:MM:SS")
#skipmissing(df)
#df1 = df[1, :]
#df2 = df[2, :]

#println(dfn)
println(names(df1))
println(names(df2))

  #=                                          # Removing milisecond part of dateTime
for i in 1:length(df.dateTime)
    # println("========")
    # println(i)
    # print(df[i,:])
    # println("---------")
    # println(df.dateTime[i])
    try
        DateTime.(df.dateTime[i],"yyyy-mm-dd HH:MM:SS")
    catch e
        println("========")
        println(i)
        println(e)
    end

end

    #println(df)
    #println(typeof(df.dateTime[7914]))

    =#


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct LCLFilter <: Filter
        lf::Float64
        rf::Float64
        cf::Float64
        lg::Float64
        rg::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int64
    end

Parameters of a LCL filter outside the converter, the states are in the grid's reference frame

# Arguments
- `lf::Float64`: Series inductance in p.u. of converter filter, validation range: (0, nothing)
- `rf::Float64`: Series resistance in p.u. of converter filter, validation range: (0, nothing)
- `cf::Float64`: Shunt capacitance in p.u. of converter filter, validation range: (0, nothing)
- `lg::Float64`: Series inductance in p.u. of converter filter to the grid, validation range: (0, nothing)
- `rg::Float64`: Series resistance in p.u. of converter filter to the grid, validation range: (0, nothing)
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`
- `n_states::Int64`
"""
mutable struct LCLFilter <: Filter
    "Series inductance in p.u. of converter filter"
    lf::Float64
    "Series resistance in p.u. of converter filter"
    rf::Float64
    "Shunt capacitance in p.u. of converter filter"
    cf::Float64
    "Series inductance in p.u. of converter filter to the grid"
    lg::Float64
    "Series resistance in p.u. of converter filter to the grid"
    rg::Float64
    ext::Dict{String, Any}
    states::Vector{Symbol}
    n_states::Int64
end

function LCLFilter(lf, rf, cf, lg, rg, ext=Dict{String, Any}(), )
    LCLFilter(lf, rf, cf, lg, rg, ext, [:ir_cnv, :ii_cnv, :vr_filter, :vi_filter, :ir_filter, :ii_filter], 6, )
end

function LCLFilter(; lf, rf, cf, lg, rg, ext=Dict{String, Any}(), )
    LCLFilter(lf, rf, cf, lg, rg, ext, )
end

# Constructor for demo purposes; non-functional.
function LCLFilter(::Nothing)
    LCLFilter(;
        lf=0,
        rf=0,
        cf=0,
        lg=0,
        rg=0,
        ext=Dict{String, Any}(),
    )
end

"""Get LCLFilter lf."""
get_lf(value::LCLFilter) = value.lf
"""Get LCLFilter rf."""
get_rf(value::LCLFilter) = value.rf
"""Get LCLFilter cf."""
get_cf(value::LCLFilter) = value.cf
"""Get LCLFilter lg."""
get_lg(value::LCLFilter) = value.lg
"""Get LCLFilter rg."""
get_rg(value::LCLFilter) = value.rg
"""Get LCLFilter ext."""
get_ext(value::LCLFilter) = value.ext
"""Get LCLFilter states."""
get_states(value::LCLFilter) = value.states
"""Get LCLFilter n_states."""
get_n_states(value::LCLFilter) = value.n_states

"""Set LCLFilter lf."""
set_lf!(value::LCLFilter, val::Float64) = value.lf = val
"""Set LCLFilter rf."""
set_rf!(value::LCLFilter, val::Float64) = value.rf = val
"""Set LCLFilter cf."""
set_cf!(value::LCLFilter, val::Float64) = value.cf = val
"""Set LCLFilter lg."""
set_lg!(value::LCLFilter, val::Float64) = value.lg = val
"""Set LCLFilter rg."""
set_rg!(value::LCLFilter, val::Float64) = value.rg = val
"""Set LCLFilter ext."""
set_ext!(value::LCLFilter, val::Dict{String, Any}) = value.ext = val
"""Set LCLFilter states."""
set_states!(value::LCLFilter, val::Vector{Symbol}) = value.states = val
"""Set LCLFilter n_states."""
set_n_states!(value::LCLFilter, val::Int64) = value.n_states = val


text:
#using FixedSizeArrays
CD = CollisionDetection

@test CD.fitsinbox([1,1,1],1.2,[0,0,0],1,2.21) == true
@test CD.fitsinbox([1,1,1],1.2,[0,0,0],1,2.19) == false
@test CD.fitsinbox([-1,-1],1.2,[0,0],1,2.21) == true

d, n = 2, 9
data = SVector{d,Float64}[2*rand(SVector{d,Float64}) for i in 1:n]  # create the list of points for the tree test
push!(data,SVector(2.0,2.0))        # The is the the top right corner
push!(data,SVector(0.0,0.0))        # This is the lower left corner
push!(data, SVector(0.90,0.5))      # This point is at the edgre of sector(0) (-,-)
radii = abs.(zeros(n+2))           # We set raduis for all points to zero to test if they fill in sectors as well
push!(radii,0.2)                  # Now only the test point has a raduis
tree = CD.Octree(data, radii)     # Create an Octree with normal (ratio=1) so that means the point is unmovable to child sector
@test tree.rootbox.data[1]==n+3   # Now we test if the test point is the only unmovable (it should be)
tree = CD.Octree(data, radii,1.2) # We create a new tree but with wider boxes to include points at sectors edge
@test in(n+3,tree.rootbox.children[1].data)==true # now the point should be located at sector zero (+1) of the direct child
push!(data, SVector(0.90,0.5))      # now we add another point but with slightly fatter raduis than the ratio 1.2
push!(radii,0.26)
tree = CD.Octree(data, radii,1.2)
@test in(n+4,tree.rootbox.children[1].data)==false # Now it should go to the child sector
@test tree.rootbox.data[1]==n+4                    # but will be definatly the only one in top level,

max_stars_repo_name: shakisparki/CollisionDetection

max_issues_repo_name: shakisparki/CollisionDetection

max_forks_repo_name: shakisparki/CollisionDetection


text:
# ----------
# Time Integrator Function for MHDFlows
# ----------

function TimeIntegrator!(prob,t₀ :: Number,N₀ :: Int;
                                       usr_dt = 0.0,
                                        diags = [],
                                  loop_number = 100,
							     	     save = false,
							         save_loc = "",
							         filename = "",
							          dump_dt = 0)

  # Check if save function related parameter
  if (save)

  	if length(save_loc) == 0 || length(filename) == 0 || dump_dt == 0
  		error("Save Function Turned ON but save_loc/filename/dump_dt is not declared!\n");
  	end
  end

  # Declare the vars update function and CFL time calclator
  updatevars! = ifelse(prob.flag.b, MHDSolver.MHDupdatevars!,
                                    HDSolver.HDupdatevars!);
  updateCFL!  = ifelse(prob.flag.b, CFL_MHD!, CFL_HD!);

  # Declare the iterator paramters
  Nᵢ = 0;
  t_next_save = prob.clock.t + dump_dt;

  # Check if user is declared a looping dt
  usr_declared_dt = usr_dt != 0.0 ? true : false
  if (usr_declared_dt)
    prob.clock.dt = usr_dt;
  end

  while (N₀ >= Nᵢ) && (t₀ >= prob.clock.t)
    if (!usr_declared_dt)
        #update the CFL condition;
      	updateCFL!(prob)
    end

  	#update the system;
  	stepforward!(prob.sol, prob.clock, prob.timestepper, prob.eqn,
                 prob.vars, prob.params, prob.grid);

    #update the diags
    increment!(diags)

    #Corret v and b if VP method is turned on
    if (prob.flag.vp == true)
      MHDSolver_VP.DivVCorrection!(prob);
      prob.flag.b == true ? MHDSolver_VP.DivBCorrection!(prob) : nothing;
    end

    #Dye Update
    prob.dye.dyeflag == true ? prob.dye.stepforward!(prob) : nothing;

  	#User defined function
  	for foo! ∈ prob.usr_func
  		foo!(prob);
  	end

  	#update the vars
	updatevars!(prob);

	Nᵢ += 1;

	#Save Section
	if (save) && prob.clock.t >= t_next_save;
      KE_ = ProbDiagnostic(prob, Nᵢ; print_ = false);
      isnan(KE_) ? error("detected NaN! Quite the simulation right now.") : nothing;
      savefile(prob, file_number; file_path_and_name = file_path_and_name)
      t_next_save += dump_dt;
      file_number +=1;
    end

    if Nᵢ % loop_number == 0
    	KE_ = ProbDiagnostic(prob, Nᵢ; print_ = true);
        isnan(KE_) ? error("detected NaN! Quite the simulation right now.") : nothing;
    end

  end

end

function CFL_MHD!(prob;Coef = 0.25)
    #Solving the dt of CFL condition using dt = Coef*dx/v

	#Maxmium velocity
    @. prob.vars.nonlin1 *=0;
    v2 = prob.vars.nonlin1;
    v2 += prob.vars.ux.^2 + prob.vars.uy.^2 + prob.vars.uz.^2;
    v2_max = maximum(v2);

	#Maxmium Alfvenic velocity
    @. prob.vars.nonlin1 *=0;
    v2a = prob.vars.nonlin1;
    v2a += prob.vars.bx.^2 + prob.vars.by.^2 + prob.vars.bz.^2;
    v2a_max = maximum(v2a);

    vmax = sqrt(maximum([v2_max,v2a_max]));
    dx = prob.grid.Lx/prob.grid.nx;
    dy = prob.grid.Ly/prob.grid.ny;
    dz = prob.grid.Lz/prob.grid.nz;
    dl = minimum([dx,dy,dz]);
    dt =  Coef*dl/vmax;
    prob.clock.dt = dt;

end

function CFL_HD!(prob;Coef = 0.25)
    #Solving the dt of CFL condition using dt = Coef*dx/v

	#Maxmium velocity
    @. prob.vars.nonlin1 *=0;
    v2 = prob.vars.nonlin1;
    v2 += prob.vars.ux.^2 + prob.vars.uy.^2 + prob.vars.uz.^2;
    vmax = sqrt(maximum(v2));

    dx = prob.grid.Lx/prob.grid.nx;
    dy = prob.grid.Ly/prob.grid.ny;
    dz = prob.grid.Lz/prob.grid.nz;
    dl = minimum([dx,dy,dz]);
    dt =  Coef*dl/vmax;
    prob.clock.dt = dt;

end

function ProbDiagnostic(prob,N; print_=false)
    dV = (2π/prob.grid.nx)^3;
    vx,vy,vz = prob.vars.ux,prob.vars.uy,prob.vars.uz;
    KE =  string(round(sum(vx.^2+vy.^2 + vz.^2)*dV,sigdigits=3));
    tt =  string(round(prob.clock.t,sigdigits=3));
    nn = string(N);
    for i = 1:8-length(string(tt));tt= " "*tt;end
    for i = 1:8-length(string(KE));KE= " "*KE;end
    for i = 1:8-length(string(nn));nn= " "*nn;end

    if (prob.flag.b == true)
        bx,by,bz = prob.vars.bx,prob.vars.by,prob.vars.bz;
        ME =  string(round(sum(bx.^2+by.^2 + bz.^2)*dV,sigdigits=3));
        for i = 1:8-length(string(ME));ME= " "*ME;end
        print_ == true ? println("n = $nn, t = $tt, KE = $KE, ME= $ME") : nothing;
    else
        print_ == true ? println("n = $nn, t = $tt, KE = $KE") : nothing;
    end

    return parse(Float32,KE)
end

function savefile(prob,ρₖ,file_number;file_path_and_name="")
    space_0 = ""
    for i = 1:4-length(string(file_number));space_0*="0";end
    fw = h5open(file_path_and_name*"_t_"*space_0*string(file_number)*".h5","w")
    write(fw, "i_velocity",  Array(prob.vars.ux));
    write(fw, "j_velocity",  Array(prob.vars.uy));
    write(fw, "k_velocity",  Array(prob.vars.uz));
    if (prob.dye.dyeflag == true)
        write(fw, "dye_density",  Array(prob.dye.ρ));
    end
    if (prob.flag.b == true)
        write(fw, "i_mag_field", Array(prob.vars.bx));
        write(fw, "j_mag_field", Array(prob.vars.by));
        write(fw, "k_mag_field", Array(prob.vars.bz));
    end
    write(fw, "time", prob.clock.t);
    close(fw)
end

function MHDintegrator!(prob,t)
	dt = prob.clock.dt;
	NStep = Int(round(t/dt));
	for i =1:NStep
	    stepforward!(prob);
	    MHDupdatevars!(prob);
	end
end

function HDintegrator!(prob,t)
	dt = prob.clock.dt;
	NStep = Int(round(t/dt));
	for i =1:NStep
	    stepforward!(prob);
	    HDupdatevars!(prob);
	end
end

max_stars_repo_name: MHDFlows/MHDFlows

max_issues_repo_name: MHDFlows/MHDFlows

max_forks_repo_name: MHDFlows/MHDFlows


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Frames (Modia3D/Frames/_module.jl)
#

"""
    const Modia3D.Quaternion = SVector{4,Float64}

Describes the rotation from a frame 1 into a frame 2 with a quaternion vector.
If `e` is the (normalized) axis of rotation to rotate frame 1 into frame 2
(either resolved in frame 1 or frame 2) and `angle`
is the rotation angle for this rotation then the quaternion vector
`q::Modia3D.Quaternions` is defined as:

```julia
q = [e*sin(angle/2),
       cos(angle/2]
```
"""
const Quaternion = SVector{4,Float64}

"""
    const Modia3D.NullQuaternion = Quaternion(0,0,0,1)

Constant Quaternion vector of a null rotation (= no rotation from frame 1 to frame 2)
"""
NullQuaternion = Quaternion(0.0, 0.0, 0.0, 1.0)

"""
    Modia3D.assertQuaternion(q::AbstractVector)

Assert that vector `q` has the properties of a `Quaternion` vector
(has 4 elements, `norm(q) = 1`)
"""
function assertQuaternion(q::AbstractVector)
    @assert(length(q) == 4)
    @assert(abs(norm(q) - 1.0) <= 1e-10)
end

"""
    R = Modia3D.from_q(q::Modia3D.Quaternion)

Return RotationMatrix `R` from Quaternion `q`.
"""
from_q(q::Quaternion) = RotationMatrix(2.0 * (q[1] * q[1] + q[4] * q[4]) - 1.0,
                                       2.0 * (q[2] * q[1] - q[3] * q[4]),
                                       2.0 * (q[3] * q[1] + q[2] * q[4]),
                                       2.0 * (q[1] * q[2] + q[3] * q[4]),
                                       2.0 * (q[2] * q[2] + q[4] * q[4]) - 1.0,
                                       2.0 * (q[3] * q[2] - q[1] * q[4]),
                                       2.0 * (q[1] * q[3] - q[2] * q[4]),
                                       2.0 * (q[2] * q[3] + q[1] * q[4]),
                                       2.0 * (q[3] * q[3] + q[4] * q[4]) - 1.0)

const p4limit = 0.1
const c4limit = 4.0 * p4limit * p4limit

"""
    q = Modia3D.from_R(R::Modia3D.RotationMatrix;
                       q_guess = NullQuaternion)

Return `Quaternion q` from `RotationMatrix R`.

From the two possible solutions `q` the one is returned that is closer
to `q_guess` (note, `q` and `-q` define the same rotation).
"""
function from_R(R::RotationMatrix;
                q_guess::Quaternion=NullQuaternion)::Quaternion
    # Based on https://d3cw3dd2w32x2b.cloudfront.net/wp-content/uploads/2015/01/matrix-to-quat.pdf
    tt::Float64 = 0.0
    if R[3,3] < 0
        if R[1,1] > R[2,2]
            tt = 1 + R[1,1] - R[2,2] - R[3,3]
            q = Quaternion(tt, R[1,2] + R[2,1], R[3,1] + R[1,3], R[2,3] - R[3,2])
        else
            tt = 1 - R[1,1] + R[2,2] - R[3,3]
            q = Quaternion(R[1,2] + R[2,1], tt, R[2,3] + R[3,2], R[3,1] - R[1,3])
        end
    else
        if R[1,1] < -R[2,2]
            tt = 1 - R[1,1] - R[2,2] + R[3,3]
            q = Quaternion(R[3,1] + R[1,3], R[2,3] + R[3,2], tt, R[1,2] - R[2,1])
        else
            tt = 1 + R[1,1] + R[2,2] + R[3,3]
            q = Quaternion(R[2,3] - R[3,2], R[3,1] - R[1,3], R[1,2] - R[2,1], tt)
        end
    end
    q = q * 0.5/sqrt(tt)

    return dot(q, q_guess) >= 0 ? q : -q
end

from_R(R::AbstractMatrix; q_guess::AbstractVector=NullQuaternion)::Quaternion =
    from_R(RotationMatrix(R); q_guess=Quaternion(q_guess))

"""
    q = Modia3D.qrot1(angle; q_guess = NullQuaternion)

Return Quaternion `q` that rotates with angle `angle` along the x-axis of frame 1.

From the two possible solutions `q` the one is returned that is closer
to `q_guess` (note, `q` and `-q` define the same rotation).
"""
@inline function qrot1(angle::Number; q_guess::Quaternion=NullQuaternion)::Quaternion
    q = Quaternion(sin(angle / 2), 0.0, 0.0, cos(angle / 2))

    return dot(q, q_guess) >= 0 ? q : -q
end

"""
    q = Modia3D.qrot2(angle; q_guess = NullQuaternion)

Return Quaternion `q` that rotates with angle `angle` along the y-axis of frame 1.

From the two possible solutions `q` the one is returned that is closer
to `q_guess` (note, `q` and `-q` define the same rotation).
"""
@inline function qrot2(angle::Number; q_guess::Quaternion=NullQuaternion)::Quaternion
    q = Quaternion(0.0, sin(angle / 2), 0.0, cos(angle / 2))

    return dot(q, q_guess) >= 0 ? q : -q
end

"""
    q = Modia3D.qrot3(angle; q_guess = NullQuaternion)

Return Quaternion `q` that rotates with angle `angle` along the z-axis of frame 1.

From the two possible solutions `q` the one is returned that is closer
to `q_guess` (note, `q` and `-q` define the same rotation).
"""
@inline function qrot3(angle::Number; q_guess::Quaternion=NullQuaternion)::Quaternion
    q = Quaternion(0.0, 0.0, sin(angle / 2), cos(angle / 2))

    return dot(q, q_guess) >= 0 ? q : -q
end

absoluteRotation(q1::Quaternion, q_rel::Quaternion)::Quaternion =
     (@SMatrix [ q_rel[4]  q_rel[3] -q_rel[2] q_rel[1];
                -q_rel[3]  q_rel[4]  q_rel[1] q_rel[2];
                 q_rel[2] -q_rel[1]  q_rel[4] q_rel[3];
                -q_rel[1] -q_rel[2] -q_rel[3] q_rel[4]]) * q1

"""
    q = Modia3D.qrot123(angle1, angle2, angle3)

Return Quaternion `q` by rotating with angle1 along the x-axis of frame 1,
then with angle2 along the y-axis of this frame and then with angle3 along
the z-axis of this frame.

From the two possible solutions `q` the one is returned that is closer
to `q_guess` (note, `q` and `-q` define the same rotation).
"""
qrot123(angle1::Number, angle2::Number, angle3::Number)::Quaternion = absoluteRotation(absoluteRotation(qrot1(angle1), qrot2(angle2)), qrot3(angle3))

"""
    q = Modia3D.qrot_e(e, angle; q_guess = NullQuaternion)

Return Quaternion `q` that rotates with angle `angle` along unit axis `e`.
This function assumes that `norm(e) == 1`.

From the two possible solutions `q` the one is returned that is closer
to `q_guess` (note, `q` and `-q` define the same rotation).
"""
@inline function qrot_e(e::Vector3D, angle::Number; q_guess::Quaternion=NullQuaternion)::Quaternion
    sa = sin(angle / 2)
    q = Quaternion(e[1] * sa, e[2] * sa, e[3] * sa, cos(angle / 2))

    return dot(q, q_guess) >= 0 ? q : -q
end
qrot_e(e::AbstractVector, angle::Number)::Quaternion = qrot_e(Vector3D(e), convert(Float64, angle))

"""
    q = Modia3D.qrot_nxy(nx, ny)

It is assumed that the two input vectors `nx` and `ny` are resolved in frame 1 and
are directed along the x and y axis of frame 2.
The function returns the Quaternion `q` to rotate from frame 1 to frame 2.

The function is robust in the sense that it returns always a Quaternion `q`,
even if `ny` is not orthogonal to `nx` or if one or both vectors have zero length.
This is performed in the following way:
If `nx` and `ny` are not orthogonal to each other, first a unit vector `ey` is
determined that is orthogonal to `nx` and is lying in the plane spanned by
`nx` and `ny`. If `nx` and `ny` are parallel or nearly parallel to each other
or `ny` is a vector with zero or nearly zero length, a vector `ey` is selected
arbitrarily such that `ex` and `ey` are orthogonal to each other.
If both `nx` and `ny` are vectors with zero or nearly zero length, an
arbitrary Quaternion `q` is returned.

# Example

```julia
using Unitful
import Modia3D

q1 = Modia3D.qrot1(90u"°")
q2 = Modia3D.qrot_nxy([1  , 0, 0], [0  , 0, 1  ])
q3 = Modia3D.qrot_nxy([0.9, 0, 0], [1.1, 0, 1.1])
isapprox(q1,q2)   # returns true
isapprox(q1,q3)   # returns true
```
"""
qrot_nxy(nx, ny)::Quaternion = from_R(rot_nxy(nx, ny))

resolve1(q::Quaternion, v2::Vector3D)::Vector3D =
    2 * ((q[4] * q[4] - 0.5) * v2 + dot(q[1:3], v2) * q[1:3] + q[4] * cross(q[1:3], v2))

resolve2(q::Quaternion, v1::Vector3D)::Vector3D =
    2 * ((q[4] * q[4] - 0.5) * v1 + dot(q[1:3], v1) * q[1:3] - q[4] * cross(q[1:3], v1))

relativeRotation(q1::Quaternion, q2::Quaternion)::Quaternion =
     (@SMatrix [ q1[4]  q1[3] -q1[2] -q1[1];
                -q1[3]  q1[4]  q1[1] -q1[2];
                 q1[2] -q1[1]  q1[4] -q1[3];
                 q1[1]  q1[2]  q1[3]  q1[4]]) * q2

inverseRotation(q::Quaternion)::Quaternion = Quaternion(-q[1], -q[2], -q[3], q[4])


text:
using churrosJulia
using Test

include("helloWorld.jl") # contains a long Hellow World program

fibDirty = "{o}=} {o}=} push 2 1's
{o}==========} {o}==} {=={o} push 10
{==={*} peek stack, exit loop if it is 0
{o}=} push 1
{====={o} put 10 in memory 1
{={*} peek top 2 spots and push their sum
{o}=} push 1
{======{o} pop 1 get memory location 1
{*}=} push -1
{={o} pop 10 -1 push 10+-1
{===={*} exit loop if 10-1!=0
{{o} {======={o} pop and print 2nd"

fibClean = "{o}=} {o}=} {o}==========} {o}==} {=={o} {==={*} {o}=} {====={o}
{={*} {o}=} {======{o} {*}=} {={o} {===={*} {{o} {======={o}"

breakTest = "{o}=} {======={o} {=========={o} {o}=} {======={o}"

skipLoops1 = "{o}} {==={*} {o}=} {===={*} {======={o}"
skipLoops2 = "{o}=} {==={*} {o}} {===={*} {======={o} {======={o}"
skipLoops3 = "{o}} {==={*} {o}=} {===={*} {o}} {======={o} {======={o}"

loopParse1 = "{==={o} {===={o} {===={o} {==={o} {==={o}"
loopParse2 = "{==={o} {===={o} {==={o} {==={o} {===={o}"

invalidMem = "{o}==} {======{o} {======={o}"

@testset "churrosJulia.jl" begin
    @test churro("")==""
    @test churro("{o}====} {======={o} {{==o}")=="4"
    @test churro("{o}====} {o}====} {={o} {======={o}")=="8"
    @test churro(HelloWorld)=="Hello, World!"
    @test churro(fibDirty)=="55"
    @test churro(fibClean)=="55"
    @test churro(breakTest)=="1"
    @test churro(skipLoops1)=="0"
    @test churro(skipLoops2)=="01"
    @test churro(skipLoops3)=="00"
    @test churro(fibClean,diagnostics = true)=="55"
    @test try
        churro(loopParse1)
        "no error"
    catch e
        e.msg
    end == "Too many loops were closed"
    @test try
        churro(loopParse2)
        "no error"
    catch e
        e.msg
    end == "Not enough loops were closed"
    @test churro(invalidMem) == "0"
end


text:
const tbbox = (min_x = 440720.0, min_y = 3.74532e6, max_x = 446720.0, max_y = 3.75132e6)

@testset "GeoUtils" begin
    @testset "bbox" begin
        ga = GeoArrays.read("data/utmsmall.tif")
        @test GeoArrays.bbox(ga) == tbbox
    end

    @testset "bbox!" begin
        ga = GeoArrays.read("data/utmsmall.tif")
        @test GeoArrays.bbox(ga) == tbbox
        GeoArrays.bbox!(ga, GeoArrays.bbox(ga))
        @test GeoArrays.bbox(ga) == tbbox
    end

    @testset "flip upside down" begin
        ga = GeoArrays.read("data/utmsmall.tif")
        ga_identical = GeoArrays.read("data/utmsmall.tif")
        @test ga == ga_identical
        GeoArrays.flipud!(ga)
        @test ga != ga_identical
        @test GeoArrays.bbox(ga) == tbbox
        GeoArrays.flipud!(ga)
        @test ga == ga_identical
    end

    @testset "geotransform should remain the same" begin
        gt = [1000.,1.,0.2,2000.,0.,1.]
        f = GeoArrays.geotransform_to_affine(gt)
        gtnew = GeoArrays.affine_to_geotransform(f)
        @test gt == gtnew
    end
end


text:
module SphericalHarmonics

include("ylm.jl")

include("clebschgordan.jl")

end


text:
using Test
using Algorithm101
const ALGO = Algorithm101

@testset "Test QuickSort" begin

    # Test: 5 elements
    data = [5,2,4,3,1]
    answer = sort(data)
    @test ALGO.quick_sort!(data, 1, length(data)) == answer

    # Test: 100 elements
    data = rand(100)
    answer = sort(data)
    @test ALGO.quick_sort!(data, 1, length(data)) == answer

end


text:
""
function constraint_current_balance(pm::AbstractPowerModel, i::Int; nw::Int=nw_id_default)
    bus_arcs      = _PMs.ref(pm, nw, :bus_arcs, i)
    bus_arcs_dc   = _PMs.ref(pm, nw, :bus_arcs_dc, i)
    bus_arcs_xfmr = _PMs.ref(pm, nw, :bus_arcs_xfmr, i)
    bus_gens      = _PMs.ref(pm, nw, :bus_gens, i)
    bus_loads     = _PMs.ref(pm, nw, :bus_loads, i)
    bus_shunts    = _PMs.ref(pm, nw, :bus_shunts, i)

    bus_gs = Dict(k => _PMs.ref(pm, nw, :shunt, k, "gs") for k in bus_shunts)
    bus_bs = Dict(k => _PMs.ref(pm, nw, :shunt, k, "bs") for k in bus_shunts)

    constraint_current_balance(pm, nw, i, bus_arcs, bus_arcs_dc, bus_arcs_xfmr, bus_gens, bus_loads, bus_gs, bus_bs)
end

""
function constraint_transformer_core_excitation(pm::AbstractPowerModel, t::Int; nw::Int=nw_id_default)
    current_harmonics_ntws = ref(pm, nw, :xfmr, t, "current_harmonics_ntws")

    if nw in current_harmonics_ntws
        int_a  = ref(pm, nw, :xfmr, t, "INT_A")
        int_b  = ref(pm, nw, :xfmr, t, "INT_B")

        constraint_transformer_core_excitation(pm, nw, t, int_a, int_b)
    else
        constraint_transformer_core_excitation(pm, nw, t)
    end
end

""
function constraint_transformer_core_voltage_drop(pm::AbstractPowerModel, t::Int; nw::Int=nw_id_default)
    f_bus = ref(pm, nw, :xfmr, t, "f_bus")
    t_bus = ref(pm, nw, :xfmr, t, "t_bus")
    f_idx = (t,f_bus,t_bus)

    xsc = ref(pm, nw, :xfmr, t, "xsc")

    constraint_transformer_core_voltage_drop(pm, nw, t, f_idx, xsc)
end

""
function constraint_transformer_core_voltage_balance(pm::AbstractPowerModel, t::Int; nw::Int=nw_id_default)
    f_bus = ref(pm, nw, :xfmr, t, "f_bus")
    t_bus = ref(pm, nw, :xfmr, t, "t_bus")
    t_idx = (t,t_bus,f_bus)

    tr = ref(pm, nw, :xfmr, t, "tr")
    ti = ref(pm, nw, :xfmr, t, "ti")

    constraint_transformer_core_voltage_balance(pm, nw, t, t_idx, tr, ti)
end

""
function constraint_transformer_core_current_balance(pm::AbstractPowerModel, t::Int; nw::Int=nw_id_default)
    f_bus = ref(pm, nw, :xfmr, t, "f_bus")
    t_bus = ref(pm, nw, :xfmr, t, "t_bus")
    f_idx = (t,f_bus,t_bus)
    t_idx = (t,t_bus,f_bus)

    tr = ref(pm, nw, :xfmr, t, "tr")
    ti = ref(pm, nw, :xfmr, t, "ti")

    constraint_transformer_core_current_balance(pm, nw, t, f_idx, t_idx, tr, ti)
end

""
function constraint_transformer_winding_config(pm::AbstractPowerModel, t::Int; nw::Int=nw_id_default)
    nh = pm.data["harmonics"]["$nw"]

    f_bus = ref(pm, nw, :xfmr, t, "f_bus")
    t_bus = ref(pm, nw, :xfmr, t, "t_bus")
    bus = [f_bus, t_bus]
    idx = [(t,f_bus,t_bus), (t,t_bus,f_bus)]

    r  = [ref(pm, nw, :xfmr, t, nk) for nk in ["r1","r2"]]
    re = [ref(pm, nw, :xfmr, t, nk) for nk in ["re1","re2"]]
    xe = [ref(pm, nw, :xfmr, t, nk) for nk in ["xe1","xe2"]]

    gnd = [ref(pm, nw, :xfmr, t, nk) for nk in ["gnd1","gnd2"]]

    for w in 1:2
        constraint_transformer_winding_config(pm, nw, nh, bus[w], idx[w], r[w], re[w], xe[w], gnd[w])
    end #                                               nh, i,      idx,    r,  re,     xe,     gnd)
end

""
function constraint_transformer_winding_current_balance(pm::AbstractPowerModel, t::Int; nw::Int=nw_id_default)
    nh = pm.data["harmonics"]["$nw"]

    f_bus = ref(pm, nw, :xfmr, t, "f_bus")
    t_bus = ref(pm, nw, :xfmr, t, "t_bus")
    idx = [(t,f_bus,t_bus), (t,t_bus,f_bus)]

    r = [ref(pm, nw, :xfmr, t, nk) for nk in ["r1","r2"]]
    # b_sh = [ref(pm, nw, :xfmr, t, nk) for nk in ["b_sh1","b_sh2"]] ## TODO
    # g_sh = [ref(pm, nw, :xfmr, t, nk) for nk in ["g_sh1","g_sh2"]]
    b_sh, g_sh = [0.0,0.0], [0.0,0.0]

    cnf = [ref(pm, nw, :xfmr, t, nk) for nk in ["cnf1","cnf2"]]

    for w in 1:2
        constraint_transformer_winding_current_balance(pm, nw, nh, idx[w], r[w], b_sh[w], g_sh[w], cnf[w])
    end
end

""
function constraint_voltage_magnitude_rms(pm::AbstractPowerModel, i::Int; fundamental::Int=1)
    vminrms = ref(pm, fundamental, :bus, i, "vminrms")
    vmaxrms = ref(pm, fundamental, :bus, i, "vmaxrms")

    constraint_voltage_magnitude_rms(pm, i, vminrms, vmaxrms)
end

""
function constraint_voltage_thd(pm::AbstractPowerModel, i::Int; fundamental::Int=1)
    thdmax = ref(pm, fundamental, :bus, i, "thdmax")
    constraint_voltage_thd(pm, i, fundamental, thdmax)
end

""
function constraint_voltage_harmonics_relative_magnitude(pm::AbstractPowerModel, i::Int; nw::Int=nw_id_default, fundamental::Int=1)
    rm = ref(pm, fundamental, :bus, i, "rm")
    hmap = pm.ref[:it][:pm][:harmonics]
    hnumber = hmap["$nw"]
    if haskey(rm, hnumber)
        constraint_voltage_harmonics_relative_magnitude(pm, nw, i, rm[hnumber], fundamental)
    end
end

function constraint_load_power(pm::AbstractPowerModel, i::Int; nw::Int=nw_id_default)
    load = ref(pm, nw, :load, i)
    busi = load["load_bus"]

    if nw == 1
        constraint_load_constant_power(pm, nw, i, busi, load["pd"], load["qd"])
    else
        constraint_load_constant_current(pm, nw, i, busi, load["multiplier"])
    end
end

""
function constraint_vm_auxiliary_variable(pm::AbstractPowerModel, i::Int; nw::Int=nw_id_default)
    constraint_vm_auxiliary_variable(pm, nw, i)
end

function constraint_ref_bus(pm::AbstractPowerModel, i::Int; nw::Int=nw_id_default)
    constraint_ref_bus(pm, nw, i)
end

function constraint_current_limit_rms(pm::AbstractPowerModel, i::Int; fundamental::Int=1)
    branch = ref(pm, fundamental, :branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)

    if haskey(branch, "c_rating_a")
        constraint_current_limit_rms(pm, f_idx, branch["c_rating_a"])
    end
end

function constraint_active_filter(pm::AbstractPowerModel, i::Int; fundamental::Int=1)
    gen = ref(pm, fundamental, :gen, i)

    if haskey(gen, "isfilter") && gen["isfilter"] == 1
        constraint_active_filter(pm, i, fundamental)
    end
end


text:
"""
Determine road maxspeeds given osm way tags dictionary.
"""
function maxspeed(tags::AbstractDict)::DEFAULT_OSM_MAXSPEED_TYPE
    maxspeed = get(tags, "maxspeed", "default")
    U = DEFAULT_OSM_MAXSPEED_TYPE

    if maxspeed != "default"
        if maxspeed isa Integer
            return maxspeed
        elseif maxspeed isa AbstractFloat
            return U(round(maxspeed))
        elseif maxspeed isa String
            if occursin("conditional", maxspeed)
                maxspeed = remove_sub_string_after(maxspeed, "conditional")
            end

            maxspeed = split(maxspeed, COMMON_OSM_STRING_DELIMITERS)

            cleaned_maxspeeds = []
            for speed in maxspeed
                speed = occursin("mph", speed) ? remove_non_numeric(speed) * KMH_PER_MPH : remove_non_numeric(speed)
                push!(cleaned_maxspeeds, speed)
            end

            return U(round(mean(cleaned_maxspeeds)))
        else
            throw(ErrorException("Maxspeed is neither a string nor number, check data quality: $maxspeed"))
        end
    else
        highway_type = get(tags, "highway", "other")
        key = getkey(DEFAULT_MAXSPEEDS, highway_type, "other")
        return U(DEFAULT_MAXSPEEDS[key])
    end
end

"""
Determine number of lanes given osm way tags dictionary.
"""
function lanes(tags::AbstractDict)::DEFAULT_OSM_LANES_TYPE
    lanes = get(tags, "lanes", "default")
    U = DEFAULT_OSM_LANES_TYPE

    if lanes != "default"
        if lanes isa Integer
            return lanes
        elseif lanes isa AbstractFloat
            return U(round(lanes))
        elseif lanes isa String
            lanes = split(lanes, COMMON_OSM_STRING_DELIMITERS)
            lanes = [remove_non_numeric(l) for l in lanes]
            return U(round(mean(lanes)))
        else
            throw(ErrorException("Lanes is neither a string nor number, check data quality: $lanes"))
        end
    else
        highway_type = get(tags, "highway", "other")
        key = getkey(DEFAULT_LANES, highway_type, "other")
        return U(DEFAULT_LANES[key])
    end
end

"""
Determind if way is a roundabout given osm way tags dictionary.
"""
is_roundabout(tags::AbstractDict)::Bool = get(tags, "junction", "") == "roundabout" ? true : false

"""
Determine oneway road attribute given osm way tags dictionary.
"""
function is_oneway(tags::AbstractDict)::Bool
    oneway = get(tags, "oneway", "")

    if oneway in ONEWAY_FALSE
        return false
    elseif oneway in ONEWAY_TRUE
        return true
    elseif is_roundabout(tags)
        return true
    else
        highway_type = get(tags, "highway", "other")
        key = getkey(DEFAULT_ONEWAY, highway_type, "other")
        return DEFAULT_ONEWAY[key]
    end
end

"""
Determine reverseway road attribute given osm way tags dictionary.
"""
is_reverseway(tags::AbstractDict)::Bool = get(tags, "oneway", "") in Set(["-1", -1]) ? true : false

"""
Determine if way is of highway type given osm way tags dictionary.
"""
is_highway(tags::AbstractDict)::Bool = haskey(tags, "highway") ? true : false

"""
Determine if way is of railway type given osm way tags dictionary.
"""
is_railway(tags::AbstractDict)::Bool = haskey(tags, "railway") ? true : false

"""
Determine if way matches the specified network type.
"""
function matches_network_type(tags::AbstractDict, network_type::Symbol)::Bool
    for (k, v) in WAY_EXCLUSION_FILTERS[network_type]
        if haskey(tags, k)
            if tags[k] in Set(v)
                return false
            end
        end
    end
    return true
end

"""
Determine if relation is a restriction.
"""
is_restriction(tags::AbstractDict)::Bool = get(tags, "type", "") == "restriction" && haskey(tags, "restriction") ? true : false

"""
Determine if a restriction is valid and has usable data.
"""
function is_valid_restriction(members::AbstractArray, ways::AbstractDict{T,Way{T}})::Bool where T <: Integer
    role_counts = DefaultDict(0)
    role_type_counts = DefaultDict(0)
    ways_set = Set{Int}()
    ways_mapping = DefaultDict(Vector)
    via_node = nothing

    for member in members
        id = member["ref"]
        type = member["type"]
        role = member["role"]

        if type == "way"
            if !haskey(ways, id) || id in ways_set
                # Cannot process missing and duplicate from/via/to ways
                return false
            else
                push!(ways_set, id)
                push!(ways_mapping[role], id)
            end

        elseif type == "node"
            via_node = id
        end

        role_counts[role] += 1
        role_type_counts["$(role)_$(type)"] += 1
    end

    if !(role_counts["from"] == 1) ||
       !(role_counts["to"] == 1) ||
       !(
           (role_type_counts["via_node"] == 1 && role_type_counts["via_way"] < 1) ||
           (role_type_counts["via_node"] < 1 && role_type_counts["via_way"] >= 1)
        )
        # Restrictions with multiple "from" and "to" members cannot be processed
        # Restrictions with multiple "via" "node" members cannot be processed
        # Restrictions with combination of "via" "node" and "via" "way" members cannot be processed
        return false
    end

    to_way = ways_mapping["to"][1]
    trailing_to_way_nodes = trailing_elements(ways[to_way].nodes)
    from_way = ways_mapping["from"][1]
    trailing_from_way_nodes = trailing_elements(ways[from_way].nodes)

    if via_node isa Integer && (!(via_node in trailing_to_way_nodes) || !(via_node in trailing_from_way_nodes))
        # Via node must be a trailing node on to_way and from_way
        return false
    end

    if haskey(ways_mapping, "via")
        try
            # Via way trailing nodes must intersect with all to_ways and from_ways
            via_ways = ways_mapping["via"]
            via_way_nodes_list = [ways[w].nodes for w in via_ways]
            via_way_nodes = join_arrays_on_common_trailing_elements(via_way_nodes_list...)
            trailing_via_way_nodes = trailing_elements(via_way_nodes)

            if isempty(intersect(trailing_via_way_nodes, trailing_to_way_nodes)) ||
            isempty(intersect(trailing_via_way_nodes, trailing_from_way_nodes))
                return false
            end
        catch
            # Restriction cannot be process - via ways cannot be joined on common trailing nodes
            return false
        end
    end

    return true
end

"""
Parse OpenStreetMap data into `Node`, `Way` and `Restriction` objects.
"""
function parse_osm_network_dict(osm_network_dict::AbstractDict, network_type::Symbol=:drive)::OSMGraph
    U = DEFAULT_OSM_INDEX_TYPE
    T = DEFAULT_OSM_ID_TYPE
    W = DEFAULT_OSM_EDGE_WEIGHT_TYPE

    ways = Dict{T,Way{T}}()
    highway_nodes = Set{Int}([])
    for way in osm_network_dict["way"]
        if haskey(way, "tags") && haskey(way, "nodes")
            tags = way["tags"]
            if is_highway(tags) && matches_network_type(tags, network_type)
                tags["maxspeed"] = maxspeed(tags)
                tags["lanes"] = lanes(tags)
                tags["oneway"] = is_oneway(tags)
                tags["reverseway"] = is_reverseway(tags)
                nds = way["nodes"]
                union!(highway_nodes, nds)
                id = way["id"]
                ways[id] = Way(id, nds, tags)
            elseif is_railway(tags) && matches_network_type(tags, network_type)
                tags["rail_type"] = haskey(tags,"railway") ? tags["railway"] : "unknown"
                tags["electrified"] = haskey(tags,"electrified") ? tags["electrified"] : "unknown"
                tags["gauge"] = haskey(tags,"gauge") ? tags["gauge"] : nothing
                tags["usage"] = haskey(tags,"usage") ? tags["usage"] :  "unknown"
                tags["name"] = haskey(tags,"name") ? tags["name"] : "unknown"
                tags["lanes"] = haskey(tags,"tracks") ? tags["tracks"] : 1
                tags["maxspeed"] = maxspeed(tags)
                tags["oneway"] = is_oneway(tags)
                tags["reverseway"] = is_reverseway(tags)
                nds = way["nodes"]
                union!(highway_nodes, nds)
                id = way["id"]
                ways[id] = Way(id, nds, tags)
            end
        end
    end

    nodes = Dict{T,Node{T}}()
    for node in osm_network_dict["node"]
        id = node["id"]
        if id in highway_nodes
            nodes[id] = Node{T}(
                id,
                GeoLocation(node["lat"], node["lon"]),
                haskey(node, "tags") ? node["tags"] : Dict{String,Any}()
            )
        end
    end

    restrictions = Dict{T,Restriction{T}}()
    if haskey(osm_network_dict, "relation")
        for relation in osm_network_dict["relation"]
            if haskey(relation, "tags") && haskey(relation, "members")
                tags = relation["tags"]
                members = relation["members"]

                if is_restriction(tags) && is_valid_restriction(members, ways)
                    restriction_kwargs = DefaultDict(Vector)
                    for member in members
                        key = "$(member["role"])_$(member["type"])"
                        if key == "via_way"
                            push!(restriction_kwargs[Symbol(key)], member["ref"])
                        else
                            restriction_kwargs[Symbol(key)] = member["ref"]
                        end
                    end

                    id = relation["id"]
                    restrictions[id] = Restriction{T}(
                        id=id,
                        tags=tags,
                        type=haskey(restriction_kwargs, :via_way) ? "via_way" : "via_node",
                        is_exclusion=occursin("no", tags["restriction"]) ? true : false,
                        is_exclusive=occursin("only", tags["restriction"]) ? true : false,
                        ;restriction_kwargs...
                    )
                end
            end
        end
    end
    return OSMGraph{U,T,W}(nodes=nodes, ways=ways, restrictions=restrictions)
end

"""
Parse OpenStreetMap data downloaded in `:xml` or `:osm` format into a dictionary consistent with data downloaded in `:json` format.
"""
function parse_xml_dict_to_json_dict(dict::AbstractDict)::AbstractDict
    # This function is needed so dict parsed from xml is consistent with dict parsed from json
    for (type, elements) in dict
        for (i, el) in enumerate(elements)
            !(el isa AbstractDict) && continue
            sub_dict::Dict{String, Any} = el
            if haskey(sub_dict, "tag")
                dict[type][i]["tags"] = Dict{String,Any}(tag["k"] => tag["v"] for tag in sub_dict["tag"] if tag["k"] isa String)
                delete!(dict[type][i], "tag")
            end

            if haskey(sub_dict, "member")
                dict[type][i]["members"] = pop!(dict[type][i], "member") # rename
            end

            if haskey(sub_dict, "nd")
                dict[type][i]["nodes"] = [nd["ref"] for nd in sub_dict["nd"]]
                delete!(dict[type][i], "nd")
            end
        end
    end

    return dict
end

"""
Parse OpenStreetMap data downloaded in `:xml` or `:osm` format into a dictionary.
"""
function osm_dict_from_xml(osm_xml_object::XMLDocument)::AbstractDict
    root_node = root(osm_xml_object)
    dict = xml_to_dict(root_node, OSM_METADATA)
    return parse_xml_dict_to_json_dict(dict)
end

"""
Reorder OpenStreetMap data downloaded in `:json` format so items are groups by type.
"""
function osm_dict_from_json(osm_json_object::AbstractDict)::AbstractDict
    dict = DefaultDict(Vector)

    for el in osm_json_object["elements"]
        push!(dict[el["type"]], el)
    end

    return dict
end

"""
Initialises the OSMGraph object from OpenStreetMap data downloaded in `:xml` or `:osm` format.
"""
function init_graph_from_object(osm_xml_object::XMLDocument, network_type::Symbol=:drive)::OSMGraph
    dict_to_parse = osm_dict_from_xml(osm_xml_object)
    return parse_osm_network_dict(dict_to_parse, network_type)
end

"""
Initialises the OSMGraph object from OpenStreetMap data downloaded in `:json` format.
"""
function init_graph_from_object(osm_json_object::AbstractDict, network_type::Symbol=:drive)::OSMGraph
    dict_to_parse = osm_dict_from_json(osm_json_object)
    return parse_osm_network_dict(dict_to_parse, network_type)
end