text:
using TensorDSL

V4 = VectorSpace(4)
V67 = VectorSpace(67)
V3 = VectorSpace(3)
V1000 = VectorSpace(1000)

A = VariableTensor("A", V4, V4, V3, V67)
B = VariableTensor("B", V4, V4, V3, V67)
C = VariableTensor("C", V67', V3')
D = VariableTensor("D", V1000)

a = FreeIndex(V4, "a")
b = FreeIndex(V4, "b")
x = FreeIndex(V67, "x")
y = FreeIndex(V3, "y")
z = FreeIndex(V67, "z")

expr = componenttensor((A + B)[a, b, y, x] ⊗ C[x', y'], a, b) + componenttensor(B[a, b, 3, 1], a, b)
expr = ((expr - expr)⊗A + (3*expr)⊗B) ⊗ D
expr = +([expr[1,2,3,4,1,1,i] for i in 1:50]...)
expr = +([(D⊗expr)[i] for i in 1:1000]...)
@time assign(expr, B => Tensor(fill(2.4, (4,4,3,67)), V4, V4, V3, V67))
println(expr.shape)
# 28.266366 seconds (11.97 M allocations: 598.007 MiB, 1.44% gc time) for Robin dict
# 27.349857 seconds (8.26 M allocations: 417.624 MiB, 1.05% gc time) for dict
# 84.748977 seconds (80.31 M allocations: 3.421 GiB, 13.89% gc time) for nothing

max_stars_repo_name: Sohva/oolong

max_issues_repo_name: Sohva/oolong

max_forks_repo_name: oolongtensor/oolong


text:
module ForwardKinematics

export Point, forward_kinematics

struct Point
  x::Float64
  y::Float64
end

# angles is an array of angles, in order from base to end, in radians.
# lengths is an array of segment lengths, in order from base to end, in any (consistent) unit.
# returns Point(x, y) of final segment.
function forward_kinematics(angles, lengths)
  if (size(angles, 1) != size(lengths, 1))
    throw(error("parameters angles and lengths must be the same size"))
  end
  x = 0.0
  y = 0.0
  angle_sum = 0.0
  for i = 1:size(angles, 1)
    angle_sum += angles[i]
    x += lengths[i] * cos(angle_sum)
    y += lengths[i] * sin(angle_sum)
  end
  Point(x, y)
end

end # module

ext: jl


text:
using Mimi
using MimiFUND
using Plots
using Statistics
using Distributions

# get the model
m = MimiFUND.get_model()

update_param!(m, :climatedynamics, :climatesensitivity, 4.5)

discounting = [i for i in 0:0.001:0.03]

results_sc = []

for i in eachindex(discounting)
    push!(results_sc, MimiFUND.compute_scc(m, year = 2020, eta = 1., prtp = discounting[i]) * 1.68)
end

results_sc = filter(x -> x >0, results_sc)

## Company Statement
Base.@kwdef mutable struct CompanyResults
    EBITDA::Float64 #Ebitda
    ghg_emissions::Float64 # GHG emissions in milions tonnes of CO2eq
    climate_costs::Union{Nothing,Vector{Float64},Float64} = nothing # Climate-related damages to society from the company, in Billion USD
end

company_a = CompanyResults(EBITDA = 25.2, ghg_emissions =  374.)

## GHG emissions valuation
using Plots

function climate_valuation!(data::CompanyResults, scc::Vector{Float64})::CompanyResults
    data.climate_costs = data.ghg_emissions .* scc * 10^(-3) # to express in billion USD
    return data
end

climate_valuation!(company_a, convert.(Float64,results_sc))

climate_adjusted_revenue = company_a.EBITDA .- company_a.climate_costs

test = fit(Normal, climate_adjusted_revenue)

cdf.(test, climate_adjusted_revenue)

cum_cdf = cdf.(test, [i for i in minimum(climate_adjusted_revenue)*0.5:1:maximum(climate_adjusted_revenue)*1.5])

p = plot([i for i in minimum(climate_adjusted_revenue)*0.5:1:maximum(climate_adjusted_revenue)*1.5], cum_cdf .* 100, label = nothing, ylabel = "Cumulative distribution function (in %)", xlabel = "Climate-adjusted EBITDA, in Billion USD")

plot!(p, [0], [100], seriestype = :bar, label = nothing)

## dynamic

growth_rate_revenue = rand(truncated(Normal(0, 0.05),0,0.05), 21)
growth_rate_emissions = rand(truncated(Normal(0, 0.05),0,0.05), 21)

historical_company_a = []

historical_proba_unsustainability = []
time_x = [i for i in 2000:2020]

for i in eachindex(growth_rate_revenue)
    if i == 1
        push!(historical_company_a, CompanyResults(EBITDA = 25.2, ghg_emissions =  374))
    else
        push!(historical_company_a, CompanyResults(EBITDA = historical_company_a[i - 1].EBITDA * (1 + growth_rate_revenue[i]),
                                                    ghg_emissions = historical_company_a[i - 1].ghg_emissions * (1 + growth_rate_emissions[i])))
    end

    discounting = [i for i in 0:0.001:0.03]

    results_sc = []

    for j in eachindex(discounting)
        push!(results_sc, MimiFUND.compute_scc(m, year = time_x[i], eta = 1., prtp = discounting[j]) * 1.68)
    end

    climate_valuation!(historical_company_a[i], convert.(Float64,results_sc))
    climate_adjusted_revenue = historical_company_a[i].EBITDA .- historical_company_a[i].climate_costs
    proba = fit(Normal, climate_adjusted_revenue)
    push!(historical_proba_unsustainability, cdf.(proba, 0.)) # proba that climate adjusted revenue be negative
    println(cdf.(proba, 0.))
end

plot(time_x, historical_proba_unsustainability .* 100, label = nothing, xlabel = "Year", ylabel = "Probability of Climate Unsustainability (in %)")

max_stars_repo_name: tlorans/SustainableFinance

max_issues_repo_name: tlorans/SustainableFinance

max_forks_repo_name: tlorans/SustainableFinance


text:
import Downloads
Downloads.download("https://www.baidu.com/","tmp/1.html")
# Downloads.download("https://api.ixiaowai.cn/api/api.php","tmp/1.jpg")

r=Downloads.request("https://www.baidu.com/";method="GET")
dump(r)

max_stars_repo_name: Rratic/notes

max_issues_repo_name: Rratic/notes

max_forks_repo_name: Rratic/notes


text:
# Create example data
x = LinRange(-2, 2, 40)
y = x.^3 .+ x.^2 .+ x
# Create a plot with light blue background
plot(x, y, backgroundcolor=0x88ccff)
# Remove the background
background(nothing)


text:
module Coloring

import DataStructures

include("topological_sort.jl")

# workaround for slow tuples
struct MyPair{T}
    first::T
    second::T
end

# workaround for julia issue #10208
Base.hash(x::MyPair{Int}, h::UInt) = hash(x.first, hash(x.second, h))

# indexed sparse set of integers
mutable struct IndexedSet
    nzidx::Vector{Int}
    empty::BitArray{1}
    nnz::Int
end

IndexedSet(n::Integer) = IndexedSet(zeros(Int, n), trues(n), 0)

function Base.push!(v::IndexedSet, i::Integer)
    if v.empty[i]  # new index
        v.nzidx[v.nnz+=1] = i
        v.empty[i] = false
    end
    return nothing
end

function Base.empty!(v::IndexedSet)
    nzidx = v.nzidx
    empty = v.empty
    for i in 1:v.nnz
        empty[nzidx[i]] = true
    end
    return v.nnz = 0
end

Base.length(v::IndexedSet) = length(v.nzidx)
function Base.resize!(v::IndexedSet, n::Integer)
    if n > length(v)
        @assert v.nnz == 0 # only resize empty vector
        resize!(v.nzidx, n)
        resize!(v.empty, n)
        fill!(v.empty, true)
    end
end

Base.collect(v::IndexedSet) = v.nzidx[1:v.nnz]
function Base.union!(v::IndexedSet, s)
    for x in s
        push!(v, x)
    end
    return nothing
end

# compact storage for an undirected graph
# neighbors of vertex i start at adjlist[offsets[i]]
struct UndirectedGraph
    adjlist::Vector{Int}
    edgeindex::Vector{Int} # corresponding edge number, indexed as adjlist
    offsets::Vector{Int}
    edges::Vector{MyPair{Int}}
end
num_vertices(g::UndirectedGraph) = length(g.offsets) - 1
num_edges(g::UndirectedGraph) = length(g.edges)
num_neighbors(i::Int, g::UndirectedGraph) = g.offsets[i+1] - g.offsets[i]
start_neighbors(i::Int, g::UndirectedGraph) = g.offsets[i]

function gen_adjlist(I, J, nel)
    adjcount = zeros(Int, nel)
    n_edges = 0
    for k in 1:length(I)
        i = I[k]
        j = J[k]
        i == j && continue
        n_edges += 1
        adjcount[i] += 1
        adjcount[j] += 1
    end
    offsets = Array{Int}(undef, nel + 1)
    offsets[1] = 1
    for k in 1:nel
        offsets[k+1] = offsets[k] + adjcount[k]
    end
    fill!(adjcount, 0)

    edges = Array{MyPair{Int}}(undef, n_edges)
    adjlist = Array{Int}(undef, offsets[nel+1] - 1)
    edgeindex = Array{Int}(undef, length(adjlist))
    edge_count = 0

    for k in 1:length(I)
        i = I[k]
        j = J[k]
        i == j && continue
        edge_count += 1
        adjlist[offsets[i]+adjcount[i]] = j
        edgeindex[offsets[i]+adjcount[i]] = edge_count
        adjcount[i] += 1

        adjlist[offsets[j]+adjcount[j]] = i
        edgeindex[offsets[j]+adjcount[j]] = edge_count
        adjcount[j] += 1

        edges[edge_count] = MyPair(i, j)
    end
    @assert edge_count == n_edges

    return UndirectedGraph(adjlist, edgeindex, offsets, edges)
end

export gen_adjlist

struct Edge
    index::Int
    source::Int
    target::Int
end

# convert to lower triangular indices, using Pairs
normalize(i, j) = (j > i) ? (j, i) : (i, j)
function normalize_p(p::MyPair)
    return (p.second > p.first) ? MyPair(p.second, p.first) :
           MyPair(p.first, p.second)
end
normalize_p(i, j) = normalize_p(MyPair(i, j))

macro colored(i)
    return esc(:((color[$i] != 0)))
end

function prevent_cycle(
    v,
    w,
    x,
    e_idx1,
    e_idx2,
    S,
    firstVisitToTree,
    forbiddenColors,
    color,
)
    er = DataStructures.find_root!(S, e_idx2)
    @inbounds first = firstVisitToTree[er]
    p = first.source # but this depends on the order?
    q = first.target
    @inbounds if p != v
        firstVisitToTree[er] = Edge(e_idx1, v, w)
    elseif q != w
        forbiddenColors[color[x]] = v
    end
    return nothing
end

function grow_star(v, w, e_idx, firstNeighbor, color, S)
    @inbounds e = firstNeighbor[color[w]]
    p = e.source
    q = e.target
    @inbounds if p != v
        firstNeighbor[color[w]] = Edge(e_idx, v, w)
    else
        union!(S, e_idx, e.index)
    end
    return nothing
end

function merge_trees(eg, eg1, S)
    e1 = DataStructures.find_root!(S, eg)
    e2 = DataStructures.find_root!(S, eg1)
    if e1 != e2
        union!(S, eg, eg1)
    end
    return nothing
end

# acyclic coloring algorithm of Gebremdehin, Tarafdar, Manne, and Pothen
# "New Acyclic and Star Coloring Algorithms with Application to Computing Hessians"
# SIAM J. Sci. Comput. 2007
function acyclic_coloring(g::UndirectedGraph)
    if num_edges(g) == 0
        return fill(1, num_vertices(g)), 1
    end
    num_colors = 0
    forbiddenColors = Int[]
    firstNeighbor = Edge[]
    firstVisitToTree = fill(Edge(0, 0, 0), num_edges(g))
    color = fill(0, num_vertices(g))
    # disjoint set forest of edges in the graph
    S = DataStructures.IntDisjointSets(num_edges(g))

    @inbounds for v in 1:num_vertices(g)
        n_neighbor = num_neighbors(v, g)
        start_neighbor = start_neighbors(v, g)
        for k in 0:(n_neighbor-1)
            w = g.adjlist[start_neighbor+k]
            @colored(w) || continue
            forbiddenColors[color[w]] = v
        end
        for k in 0:(n_neighbor-1)
            w = g.adjlist[start_neighbor+k]
            e_idx = g.edgeindex[start_neighbor+k]
            @colored(w) || continue
            n_neighbor_w = num_neighbors(w, g)
            start_neighbor_w = start_neighbors(w, g)
            for k2 in 0:(n_neighbor_w-1)
                x = g.adjlist[start_neighbor_w+k2]
                e2_idx = g.edgeindex[start_neighbor_w+k2]
                @colored(x) || continue
                if forbiddenColors[color[x]] != v
                    prevent_cycle(
                        v,
                        w,
                        x,
                        e_idx,
                        e2_idx,
                        S,
                        firstVisitToTree,
                        forbiddenColors,
                        color,
                    )
                end
            end
        end

        # find feasible color
        found = false
        for k in 1:num_colors
            if forbiddenColors[k] != v
                color[v] = k
                found = true
                break
            end
        end
        if !found
            num_colors += 1
            push!(forbiddenColors, 0)
            push!(firstNeighbor, Edge(0, 0, 0))
            color[v] = num_colors
        end

        for k in 0:(n_neighbor-1)
            w = g.adjlist[start_neighbor+k]
            e_idx = g.edgeindex[start_neighbor+k]
            @colored(w) || continue
            grow_star(v, w, e_idx, firstNeighbor, color, S)
        end
        for k in 0:(n_neighbor-1)
            w = g.adjlist[start_neighbor+k]
            e_idx = g.edgeindex[start_neighbor+k]
            @colored(w) || continue
            n_neighbor_w = num_neighbors(w, g)
            start_neighbor_w = start_neighbors(w, g)
            for k2 in 0:(n_neighbor_w-1)
                x = g.adjlist[start_neighbor_w+k2]
                e2_idx = g.edgeindex[start_neighbor_w+k2]
                (@colored(x) && x != v) || continue
                if color[x] == color[v]
                    merge_trees(e_idx, e2_idx, S)
                end
            end
        end
    end

    return color, num_colors
end

struct RecoveryInfo
    vertexmap::Vector{Vector{Int}}
    postorder::Vector{Vector{Int}}
    parents::Vector{Vector{Int}}
    color::Vector{Int}
    num_colors::Int
    nnz::Int # number of off-diagonal nonzeros
    local_indices::Vector{Int} # map back to global indices
end

function RecoveryInfo()
    return RecoveryInfo(
        Vector{Vector{Int}}(undef, 0),
        Vector{Vector{Int}}(undef, 0),
        Vector{Vector{Int}}(undef, 0),
        Vector{Int}(undef, 0),
        0,
        0,
        Vector{Int}(undef, 0),
    )
end

function recovery_preprocess(
    g::UndirectedGraph,
    color,
    num_colors,
    local_indices,
)
    # represent two-color subgraph as:
    # list of vertices (with map to global indices)
    # adjacency list in a single vector (with list of offsets)

    # linear index of pair of colors
    twocolorindex = zeros(Int32, num_colors, num_colors)
    seen_twocolors = 0
    # count of edges in each subgraph
    edge_count = Int[]
    for k in 1:length(g.edges)
        e = g.edges[k]
        u = e.first
        v = e.second
        i = min(color[u], color[v])
        j = max(color[u], color[v])
        if twocolorindex[i, j] == 0
            seen_twocolors += 1
            twocolorindex[i, j] = seen_twocolors
            push!(edge_count, 0)
        end
        idx = twocolorindex[i, j]
        edge_count[idx] += 1
    end
    # edges sorted by twocolor subgraph
    sorted_edges = Array{Vector{MyPair{Int}}}(undef, seen_twocolors)
    for idx in 1:seen_twocolors
        sorted_edges[idx] = MyPair{Int}[]
        sizehint!(sorted_edges[idx], edge_count[idx])
    end

    for i in 1:length(g.edges)
        e = g.edges[i]
        u = e.first
        v = e.second
        i = min(color[u], color[v])
        j = max(color[u], color[v])
        idx = twocolorindex[i, j]
        push!(sorted_edges[idx], MyPair(u, v))
    end

    # list of unique vertices in each twocolor subgraph
    vertexmap = Array{Vector{Int}}(undef, seen_twocolors)

    postorder = Array{Vector{Int}}(undef, seen_twocolors)
    parents = Array{Vector{Int}}(undef, seen_twocolors)

    # temporary lookup map from global index to subgraph index
    revmap = zeros(Int, num_vertices(g))

    adjcount = zeros(Int, num_vertices(g))

    cmap = zeros(Int, 0) # shared storage for DFS
    vertex_stack = Int[]
    index_stack = Int[]

    for idx in 1:seen_twocolors
        my_edges = sorted_edges[idx]
        vlist = Int[]

        # build up the vertex list and adjacency count
        for k in 1:length(my_edges)
            e = my_edges[k]
            u = e.first
            v = e.second
            # seen these vertices yet?
            if revmap[u] == 0
                push!(vlist, u)
                revmap[u] = length(vlist)
                adjcount[u] = 0
            end
            if revmap[v] == 0
                push!(vlist, v)
                revmap[v] = length(vlist)
                adjcount[v] = 0
            end
            adjcount[u] += 1
            adjcount[v] += 1
        end

        # set up offsets for adjlist
        offset = Array{Int}(undef, length(vlist) + 1)
        offset[1] = 1
        for k in 1:length(vlist)
            offset[k+1] = offset[k] + adjcount[vlist[k]]
            adjcount[vlist[k]] = 0
        end
        # adjlist for node u in twocolor idx starts at
        # vec[offset[u]]
        # u has global index vlist[u]
        vec = Array{Int}(undef, offset[length(vlist)+1] - 1)

        # now fill in
        for k in 1:length(my_edges)
            e = my_edges[k]
            u = e.first
            v = e.second

            u_rev = revmap[u] # indices in the subgraph
            v_rev = revmap[v]
            vec[offset[u_rev]+adjcount[u]] = v_rev
            vec[offset[v_rev]+adjcount[v]] = u_rev

            adjcount[u] += 1
            adjcount[v] += 1
        end

        resize!(cmap, length(vlist))
        order, parent =
            reverse_topological_sort_by_dfs(vec, offset, length(vlist), cmap)

        for k in 1:length(vlist)
            # clear for reuse
            revmap[vlist[k]] = 0
        end

        postorder[idx] = order
        parents[idx] = parent
        vertexmap[idx] = vlist
    end

    return RecoveryInfo(
        vertexmap,
        postorder,
        parents,
        color,
        num_colors,
        num_edges(g),
        local_indices,
    )
end

function indirect_recover_structure(rinfo::RecoveryInfo)
    N = length(rinfo.color)
    nnz = rinfo.nnz

    I = zeros(Int, nnz + N)
    J = zeros(Int, nnz + N)

    # diagonal entries
    k = 0
    for i in 1:N
        k += 1
        I[k] = i
        J[k] = i
    end

    for t in 1:length(rinfo.postorder)
        vmap = rinfo.vertexmap[t]
        order = rinfo.postorder[t]
        parent = rinfo.parents[t]

        for z in 1:length(order)
            v = order[z]
            p = parent[v]
            (p == 0) && continue
            i = vmap[v]
            j = vmap[p]
            i, j = normalize(i, j)
            k += 1
            I[k] = i
            J[k] = j
        end
    end

    @assert k == nnz + N

    return I, J
end

# edgelist is nonzeros in hessian, *including* nonzeros on the diagonal
function hessian_color_preprocess(
    edgelist,
    num_total_var,
    seen_idx = IndexedSet(0),
)
    resize!(seen_idx, num_total_var)
    I = Int[]
    J = Int[]
    for (i, j) in edgelist
        push!(seen_idx, i)
        push!(seen_idx, j)
        push!(I, i)
        push!(J, j)
    end
    local_indices = sort!(seen_idx.nzidx[1:seen_idx.nnz])
    empty!(seen_idx)

    global_to_local_idx = seen_idx.nzidx # steal for storage
    for k in 1:length(local_indices)
        global_to_local_idx[local_indices[k]] = k
    end
    # only do the coloring on the local indices
    for k in 1:length(I)
        I[k] = global_to_local_idx[I[k]]
        J[k] = global_to_local_idx[J[k]]
    end

    g = gen_adjlist(I, J, length(local_indices))

    color, num_colors = acyclic_coloring(g)

    @assert length(color) == num_vertices(g)

    rinfo = recovery_preprocess(g, color, num_colors, local_indices)

    I, J = indirect_recover_structure(rinfo)
    #nnz = num_edges(g)
    # convert back to global indices
    for k in 1:length(I)
        I[k] = local_indices[I[k]]
        J[k] = local_indices[J[k]]
    end

    return I, J, rinfo
end

export hessian_color_preprocess

# allocate a seed matrix
function seed_matrix(rinfo::RecoveryInfo)
    return Array{Float64}(undef, length(rinfo.local_indices), rinfo.num_colors)
end

export seed_matrix

function prepare_seed_matrix!(R, rinfo::RecoveryInfo)
    N = length(rinfo.color) # number of local variables
    @assert N == length(rinfo.local_indices)
    @assert size(R, 1) == N
    @assert size(R, 2) == rinfo.num_colors

    fill!(R, 0.0)
    for i in 1:N
        R[i, rinfo.color[i]] = 1
    end
    return
end

export prepare_seed_matrix!

# recover the hessian values from the hessian-matrix solution
# stored_values is a temporary vector of length >= length(rinfo.local_indices)
function recover_from_matmat!(V, R, rinfo::RecoveryInfo, stored_values)
    N = length(rinfo.color) # number of local variables
    nnz = rinfo.nnz

    @assert length(V) == nnz + N
    @assert length(stored_values) >= length(rinfo.local_indices)

    # diagonal entries
    k = 0

    for i in 1:N
        k += 1
        V[k] = R[i, rinfo.color[i]]
    end

    for t in 1:length(rinfo.vertexmap)
        vmap = rinfo.vertexmap[t]
        order = rinfo.postorder[t]
        parent = rinfo.parents[t]

        for z in 1:length(order)
            @inbounds stored_values[z] = 0.0
        end

        @inbounds for z in 1:length(order)
            v = order[z]
            p = parent[v]
            (p == 0) && continue

            i = vmap[v]
            j = vmap[p]

            value = R[i, rinfo.color[j]] - stored_values[v]
            stored_values[p] += value

            k += 1
            V[k] = value
        end
    end

    @assert k == nnz + N

    return
end

export recover_from_matmat!

end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function forward(backend::GPUBackend, state::InnerProductLayerState, inputs::Vector{Blob})
  M = size(state.W, 2)   # target dim
  K = size(state.W, 1)   # source dim
  dtype = eltype(state.W)
  for i = 1:length(inputs)
    input = inputs[i]
    N = get_num(input)   # batch size
    output = state.blobs[i]
    # output = W^T * X
    CuBLAS.gemm(backend.cublas_ctx, CuBLAS.OP_T, CuBLAS.OP_N, M, N, K, convert(dtype, 1),
                state.W.ptr, K, input.ptr, K, convert(dtype, 0), output.ptr, M)
    # output += bias
    CuBLAS.gemm(backend.cublas_ctx, CuBLAS.OP_N, CuBLAS.OP_N, M, N, 1, convert(dtype, 1),
                state.b.ptr, M, state.bias_multipliers[i].ptr, 1, convert(dtype, 1), output.ptr, M)
  end
end

function backward(backend::GPUBackend, state::InnerProductLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
  target_dim = size(state.W, 2)
  source_dim = size(state.W, 1)
  data_type  = eltype(state.W)

  # used in BLAS, at first it is zero, indicating overwriting the data
  # then it becomes one, indicating adding to the data
  zero_and_then_one = convert(data_type, 0)

  for i = 1:length(inputs)
    # ∂f/∂W = input * [∂f/∂o]^T
    input = inputs[i]
    batch_size = get_num(input)
    ∂f_∂o = state.blobs_diff[i]

    if !state.frozen
      CuBLAS.gemm(backend.cublas_ctx, CuBLAS.OP_N, CuBLAS.OP_T, source_dim, target_dim, batch_size,
          one(data_type), input.ptr, source_dim, ∂f_∂o.ptr, target_dim, zero_and_then_one, state.∇W.ptr, source_dim)

      # ∂f/∂b = sum(∂f/∂o, 2)
      CuBLAS.gemm(backend.cublas_ctx, CuBLAS.OP_N, CuBLAS.OP_N, target_dim, 1, batch_size,
          one(data_type), ∂f_∂o.ptr, target_dim, state.bias_multipliers[i].ptr, batch_size, zero_and_then_one, state.∇b.ptr, target_dim)
    end

    zero_and_then_one = convert(data_type, 1)

    # if back propagate down
    if isa(diffs[i], CuTensorBlob)
      # ∂f/∂x = W * [∂f/∂o]
      CuBLAS.gemm(backend.cublas_ctx, CuBLAS.OP_N, CuBLAS.OP_N, source_dim, batch_size, target_dim,
          convert(data_type, 1), state.W.ptr, source_dim, ∂f_∂o.ptr, target_dim, convert(data_type, 0), diffs[i].ptr, source_dim)
    end
  end
end


text:
function load_rawdata(pattern; benchmarkstore = default_benchmarkstore())
    rawdata = map(JSON.parsefile, readdir(pattern, benchmarkstore))
    sort!(rawdata, by = data -> data["metadata"]["nthreads"])
    return rawdata
end

load_task_overheads(pattern = glob"*/task_overheads.json"; kw...) =
    Iterators.map(load_rawdata(pattern; kw...)) do data
        nthreads = data["metadata"]["nthreads"]
        Iterators.map(data["results"]) do result
            (
                nthreads = nthreads,
                spawn_overhead = result["spawn_overhead"],
                sync_overhead = result["sync_overhead"],
            )
        end
    end |>
    Iterators.flatten

load_notify_overheads(pattern = glob"*/notify_overheads.json"; kw...) =
    Iterators.map(load_rawdata(pattern; kw...)) do data
        nthreads = data["metadata"]["nthreads"]
        Iterators.map(data["results"]) do result
            (
                nthreads = nthreads,
                notified = collect(Int, result["notified"]),
                pre = result["pre"],
                post = result["post"],
            )
        end
    end |>
    Iterators.flatten


text:
import WriteVTK: vtk_grid, vtk_save

@doc raw"""
    save_scfres(filename, scfres, Val(:vtk))

The function takes in the VTK filename and the scfres structure and stores into a VTK file.

Parameters
- `save_ψ`: Store the orbitals or not. By default they are not stored.

Grid Values:
- ``\rho`` -> Density in real space
- ``\psi \_k(i)\_band(j)\_real`` -> Real values of Bloch waves in real space
- ``\psi \_k(i)\_band(j)\_imag`` -> Imaginary values of Bloch waves in real space
- ``\rho spin`` -> Real value of ρspin are stored if ρspin in present

MetaData:
- energies, eigenvalues, Fermi level and occupations.
"""
function save_scfres_master(filename::AbstractString, scfres::NamedTuple, ::Val{:vts};
                            save_ψ=false, extra_data=Dict{String,Any}())
    !mpi_master() && error(
        "This function should only be called on MPI master after the k-point data has " *
        "been gathered with `gather_kpts`."
    )

    # Initialzing the VTK Grid
    basis = scfres.basis
    grid = zeros(3, basis.fft_size...)
    for (idcs, r) in zip(CartesianIndices(basis.fft_size), r_vectors_cart(basis))
        grid[:, Tuple(idcs)...] = r
    end
    vtkfile = vtk_grid(filename, grid)

    # Storing the bloch waves
    if save_ψ
        for ik in 1:length(basis.kpoints)
            for iband in 1:size(scfres.ψ[1])[2]
                ψnk_real = G_to_r(basis, basis.kpoints[ik], scfres.ψ[ik][:, iband])
                vtkfile["ψ_k$(ik)_band$(iband)_real"] = real.(ψnk_real)
                vtkfile["ψ_k$(ik)_band$(iband)_imag"] = imag.(ψnk_real)
            end
        end
    end

    # Storing the  density in real space
    vtkfile["ρ"] = scfres.ρ.real

    # Storing ρspin if it is present.
    isnothing(scfres.ρspin)  || (vtkfile["ρspin"] = scfres.ρspin.real)

    # Storing the energy components
    for key in keys(scfres.energies)
        vtkfile["energy_$(key)"] = scfres.energies[key]
    end

    # Storing the Fermi level
    vtkfile["εF"] = scfres.εF

    # Storing the EigenValues as a matrix
    vtkfile["eigenvalues"] = hcat(scfres.eigenvalues...)

    # Storing the Occupation as a matrix
    vtkfile["occupation"] = hcat(scfres.occupation...)

    for (key, value) in pairs(extra_data)
        vtkfile[key] = value
    end

    only(vtk_save(vtkfile))
end


text:
using CSV

input = CSV.read("./src/Day1/input.txt", Tuple, header = false)
inputarray = input[1]
# inputarray = [1721
# 979
# 366
# 299
# 675
# 1456]

for i in 1:(length(inputarray) - 1), j in (i+1):length(inputarray)
    sum = inputarray[i] + inputarray[j]
    if sum == 2020
        prod = inputarray[i] * inputarray[j]
        println(inputarray[i])
        println(inputarray[j])
        println(prod)
    end
end

function sumfinder(inputarray, numentries, target, currentsum, currentvals)
    if length(inputarray) < numentries
        # no solution is possible
        return nothing
    elseif length(inputarray) == numentries
        if sum(inputarray) + currentsum == target
            # solution found
            soln = vcat(currentvals, inputarray)
            return soln
        else
            # solution does not exist
            return nothing
        end
    elseif length(inputarray) == 0
        # solution does not exist
        return nothing
    elseif numentries <= 0
        if currentsum == target
            # solution found
            soln = currentvals
            return soln
        else
            # solution does not exist
            return nothing
        end
    else
        for i in 1:(length(inputarray) - numentries)
            picked = inputarray[i]
            newarray = inputarray[(i+1):end]
            newsum = currentsum + picked
            newvals = vcat(currentvals, picked)
            test = sumfinder(newarray, numentries - 1, target, newsum, newvals)
            if !isnothing(test)
                soln = test
                return soln
            end
        end
        return nothing
    end
end

a = sumfinder(inputarray, 3, 2020, 0, [])
prod(a)


text:
function _contains(convex_sets::COSMO.CompositeConvexSet, t::Type{<:COSMO.AbstractConvexCone})
  for set in convex_sets.sets
    if typeof(set) <: t
      return true
    end
  end
  return false
end

function chordal_decomposition!(ws::COSMO.Workspace{T}) where {T <: AbstractFloat}
  # do nothing if no psd cones present in the problem
  if !_contains(ws.p.C, DecomposableCones)
    ws.ci.decompose = false
    return nothing
  end
  ws.ci = ChordalInfo{T}(ws.p, ws.settings)

  find_sparsity_patterns!(ws)

  if ws.ci.num_decomposable > 0

    # Do transformation similar to clique tree based transformation in SparseCoLo
    if ws.settings.compact_transformation
      augment_clique_based!(ws)
    else
      # find transformation matrix H and new composite convex set
      find_decomposition_matrix!(ws)
      # augment the original system
      augment_system!(ws)
    end
    pre_allocate_variables!(ws)
    ws.ci.decompose = true
  else
    ws.ci.decompose = false
  end
  nothing
end

# analyse PSD cone constraints for chordal sparsity pattern
function find_sparsity_patterns!(ws::COSMO.Workspace)
  row_ranges = get_set_indices(ws.p.C.sets)
  sp_ind = 1
  for (k, C) in enumerate(ws.p.C.sets)
    psd_row_range = row_ranges[k]
    csp = find_aggregate_sparsity(ws.p.A, ws.p.b, psd_row_range, C)
    sp_ind = analyse_sparsity_pattern!(ws.ci, csp, ws.p.C.sets, C, k, psd_row_range, sp_ind, ws.settings.merge_strategy)
  end
end

analyse_sparsity_pattern!(ci::ChordalInfo, csp::Array{Int64, 1}, sets::Vector{AbstractConvexSet}, C::AbstractConvexSet, k::Int64, psd_row_range::UnitRange{Int64}, sp_ind::Int64, merge_strategy::Union{Type{<: AbstractMergeStrategy}, OptionsFactory{<: AbstractMergeStrategy}}) = sp_ind

function analyse_sparsity_pattern!(ci::ChordalInfo, csp::Array{Int64, 1}, sets::Vector{AbstractConvexSet}, C::DecomposableCones{T}, k::Int64, psd_row_range::UnitRange{Int64}, sp_ind::Int64, merge_strategy::Union{Type{<: AbstractMergeStrategy}, OptionsFactory{<: AbstractMergeStrategy}}) where {T <: Real}
  if length(csp) < C.dim
    return _analyse_sparsity_pattern(ci, csp, sets, C, k, psd_row_range, sp_ind, merge_strategy)
  else
   sets[k] = COSMO.DenseEquivalent(C, C.dim)
   return sp_ind
 end
end

function _analyse_sparsity_pattern(ci::ChordalInfo{T}, csp::Array{Int64, 1}, sets::Vector{AbstractConvexSet}, C::Union{PsdCone{T}, PsdConeTriangle{T}}, k::Int64, psd_row_range::UnitRange{Int64}, sp_ind::Int64, merge_strategy::Union{Type{<: AbstractMergeStrategy}, OptionsFactory{<: AbstractMergeStrategy}}) where {T <: AbstractFloat}
  ordering, nz_ind_map = find_graph!(ci, csp, C.sqrt_dim, C)
  sp = COSMO.SparsityPattern(ci.L, C.sqrt_dim, ordering, merge_strategy, psd_row_range, k, nz_ind_map)
  # if after analysis of SparsityPattern & clique merging only one clique remains, don't bother decomposing
  if num_cliques(sp.sntree) == 1
    sets[k] = DenseEquivalent(C, C.dim)
    return sp_ind
  else
    ci.sp_arr[sp_ind] = sp
    push!(ci.psd_cones_ind, k)
    ci.num_decomposable += 1
    return sp_ind + 1
  end
end

DenseEquivalent(C::COSMO.PsdCone{T}, dim::Int64) where {T <: AbstractFloat} = COSMO.DensePsdCone{T}(dim)
DenseEquivalent(C::COSMO.PsdConeTriangle{T}, dim::Int64) where {T <: AbstractFloat} = COSMO.DensePsdConeTriangle{T}(dim)

function nz_rows(a::SparseMatrixCSC{T}, ind::UnitRange{Int64}, DROP_ZEROS_FLAG::Bool) where {T <: AbstractFloat}
  DROP_ZEROS_FLAG && dropzeros!(a)
  active = falses(length(ind))
  for r in a.rowval
    if in(r, ind)
      active[r - ind.start + 1] = true
    end
  end
  return findall(active)
end

function number_of_overlaps_in_rows(A::SparseMatrixCSC{T}) where {T <: AbstractFloat}
  # sum the entries row-wise
  numOverlaps = sum(A, dims = 2)
  ri = findall(x -> x > 1, numOverlaps)
  return ri, numOverlaps[ri]
end

function find_aggregate_sparsity(A::SparseMatrixCSC{T}, b::AbstractVector{T}, ind::UnitRange{Int64}, C::DecomposableCones{T}) where {T <: AbstractFloat}
  AInd = nz_rows(A, ind, false)
  # commonZeros = AInd[find(x->x==0,b[AInd])]
  bInd = findall(x -> x != 0, view(b, ind))
  commonNZeros = union(AInd, bInd)
  return commonNZeros
end
find_aggregate_sparsity(A::SparseMatrixCSC{T}, b::AbstractVector{T}, ind::UnitRange{Int64}, C::AbstractConvexSet{T}) where {T <: AbstractFloat} = Int64[]

"""
    reverse_decomposition!(ws::COSMO.Workspace, settings::COSMO.Settings)

After the problem has beend solved, undo the chordal decomposition.

Depending on the kind of transformation that was used, this involves:

- Reassembling the original matrix S from its blocks
- Reassembling the dual variable MU and performing a positive semidefinite completion.
"""
function reverse_decomposition!(ws::COSMO.Workspace{T}, settings::COSMO.Settings{T}) where {T <: AbstractFloat}

  mO = ws.ci.originalM
  nO = ws.ci.originalN

  vars = Variables{T}(mO, nO, ws.ci.originalC)
  vars.x .= ws.vars.x[1:nO]

  if settings.compact_transformation
    # reassemble the original variables s and μ
    add_sub_blocks!(vars.s, ws.vars.s, vars.μ, ws.vars.μ, ws.ci, ws.p.C, ws.ci.originalC, ws.ci.cone_map)
  else
    H = ws.ci.H
    vars.s  .= SplitVector{T}(H * ws.vars.s[mO + 1:end], ws.ci.originalC)
    fill_dual_variables!(ws, vars)
  end

  ws.p.C = ws.ci.originalC
  # if user requests, perform positive semidefinite completion on entries of μ that were not in the decomposed blocks
  ws.vars = vars
  settings.complete_dual && psd_completion!(ws)

  return nothing
end

function fill_dual_variables!(ws::COSMO.Workspace{T}, vars::COSMO.Variables{T}) where {T <: AbstractFloat}
  mO = ws.ci.originalM
  H = ws.ci.H

  # this performs the operation μ = sum H_k^T *  μ_k causing an addition of (identical valued) overlapping blocks
  vars.μ .= H * ws.vars.μ[mO + 1:end]

  # # to remove the overlaps we take the average of the values for each overlap by dividing by the number of blocks that overlap in a particular entry, i.e. number of 1s in each row of H
  rowInd, nnzs = number_of_overlaps_in_rows(H)

  for iii=1:length(rowInd)
    ri = rowInd[iii]
    vars.μ[ri] .= vars.μ[ri] / nnzs[iii]
  end
  return nothing
end

function add_sub_blocks!(s::SplitVector{T}, s_decomp::SplitVector{T}, μ::AbstractVector{T}, μ_decomp::AbstractVector{T}, ci::ChordalInfo{T}, C::CompositeConvexSet{T}, C0::CompositeConvexSet{T}, cone_map::Dict{Int64, Int64}) where {T <: AbstractFloat}
  sp_arr = ci.sp_arr
  row_start = 1 # the row pointer in the decomposed problem
  row_ranges = get_set_indices(C0.sets) # the row ranges of the same cone (or "parent" cone) in the original problem

  # iterate over all the cones of the decomposed problems and add the entries into the correct positions of the original problem
  for (k, C) in enumerate(C.sets)
    row_range = row_ranges[cone_map[k]]
    row_start = add_blocks!(s, μ, row_start, row_range, sp_arr, s_decomp, μ_decomp, C)
  end
  return nothing
end

function add_blocks!(s::SplitVector{T}, μ::AbstractVector{T}, row_start::Int64, row_range::UnitRange{Int64}, sp_arr::Array{SparsityPattern, 1}, s_decomp::SplitVector{T}, μ_decomp::AbstractVector{T}, C::AbstractConvexSet{T}) where {T <: AbstractFloat}

  @. s.data[row_range] = s_decomp.data[row_start:row_start + C.dim - 1]
  @. μ[row_range] = μ_decomp[row_start:row_start + C.dim - 1]
  return row_start + C.dim
end

function add_blocks!(s::SplitVector{T}, μ::AbstractVector{T}, row_start::Int64, row_range::UnitRange{Int64}, sp_arr::Array{SparsityPattern, 1}, s_decomp::SplitVector{T}, μ_decomp::AbstractVector{T}, C::DecomposableCones{T}) where {T <: AbstractFloat}
  # load the appropriate sparsity_pattern
  sp = sp_arr[C.tree_ind]
  sntree = sp.sntree
  ordering = sp.ordering
  #row_start = sp.row_range.start
  N = length(ordering)
  clique = map(v -> ordering[v], get_clique(sntree, C.clique_ind))
  sort!(clique)
  counter = 0
  for j in clique, i in clique
    if isa(C, PsdCone) || i <= j
      offset = COSMO.vectorized_ind(i, j, N, C) - 1
      s.data[row_range.start + offset] += s_decomp.data[row_start + counter]

      # notice: this overwrites the overlapping entries
      μ[row_range.start + offset] = μ_decomp[row_start + counter]
      counter += 1
    end
  end
  row_start += get_blk_rows(length(clique), C)

  return row_start
end

# The psd entries of μ that correspond to the zeros in s are not constrained by the problem
# however, in order to make the dual psd cone positive semidefinite we have to do a
# positive semidefinite completion routine to choose the values
function psd_completion!(ws::COSMO.Workspace)

  # loop over psd cones
  row_ranges = get_set_indices(ws.p.C.sets)
  sp_ind = 1
  for (kkk, C) in enumerate(ws.p.C.sets)
    sp_ind = complete!(ws.vars.μ, C, ws.ci.sp_arr, sp_ind, row_ranges[kkk])
  end

  return nothing
end

complete!(μ::AbstractVector{T}, ::AbstractConvexSet{T}, sp_arr::Array{SparsityPattern}, sp_ind::Int64, rows::UnitRange{Int64}) where {T <: AbstractFloat} = sp_ind

function complete!(μ::AbstractVector{T}, C::PsdCone{T}, sp_arr::Array{SparsityPattern}, sp_ind::Int64, rows::UnitRange{Int64}) where {T <: AbstractFloat}
  sp = sp_arr[sp_ind]

  μ_view = view(μ, rows)
  # make this y = -μ
  @. μ_view *= -one(T)

  M = reshape(μ_view, C.sqrt_dim, C.sqrt_dim)

  psd_complete!(M, C.sqrt_dim, sp.sntree, sp.ordering)

  @. μ_view *= -one(T)
  return sp_ind + 1
end

function complete!(μ::AbstractVector{T}, C::PsdConeTriangle{T}, sp_arr::Array{SparsityPattern}, sp_ind::Int64, rows::UnitRange{Int64}) where {T <: AbstractFloat}
  sp = sp_arr[sp_ind]

  μ_view = view(μ, rows)

  # I want to psd complete y, which is -μ
  populate_upper_triangle!(C.X, -μ_view, one(T) / sqrt(T(2)))
  psd_complete!(C.X, C.sqrt_dim, sp.sntree, sp.ordering)
  extract_upper_triangle!(C.X, μ_view, sqrt(2))
  @. μ_view *= -one(T)
  return sp_ind + 1
end

# positive semidefinite completion (from Vandenberghe - Chordal Graphs..., p. 362)
# input: A - positive definite completable matrix
function psd_complete!(A::AbstractMatrix{T}, N::Int64, sntree::SuperNodeTree, p::Array{Int64}) where {T <: AbstractFloat}

  # if a clique graph based merge strategy was used for this sparsity pattern, recompute a valid clique tree
  #recompute_clique_tree(sntree.strategy) && clique_tree_from_graph!(sntree, sntree.strategy)

  ip = invperm(p)

  As = Symmetric(A, :U)
  # permutate matrix based on ordering p (p must be a vector type), W is in the order that the cliques are based on
  W = copy(As[p, p])
  W = Matrix(W)
  num_cliques = COSMO.num_cliques(sntree)

  # go through supernode tree in descending order (given a post-ordering). This is ensured in the get_snd, get_sep functions
  for j = (num_cliques - 1):-1:1

    # in order to obtain ν, α the vertex numbers of the supernode are mapped to the new position of the permuted matrix
    # index set of snd(i) sorted using the numerical ordering i,i+1,...i+ni
    ν = get_snd(sntree, j)
    #clique = get_clique(sntree, snd_id)
    # index set containing the elements of col(i) \ snd(i) sorted using numerical ordering σ(i)
    α = get_sep(sntree, j)

    # index set containing the row indizes of the lower-triangular zeros in column i (i: representative index) sorted by σ(i)
    i = ν[1]
    η = collect(i+1:1:N)
    # filter out elements in lower triangular part of column i that are non-zero
    filter!(x -> !in(x, α) && !in(x, ν), η)

    Waa = W[α, α]
    Wαν = view(W, α, ν)
    Wηα = view(W, η, α)

    Y = zeros(length(α), length(ν))
    try
     Y[:, :] = Waa \ Wαν
   catch
    Waa_pinv = pinv(Waa)
    Y[:, :] = Waa_pinv * Wαν
  end

  W[η, ν] =  Wηα * Y
  # symmetry condition
  W[ν, η] = view(W, η, ν)'
end

# invert the permutation
A[:, :] =  W[ip, ip]
end


text:
export  PhyPolyHermite, Cphy,
        degree, PhyPolyH,
        phyhermite_coeffmatrix,
        FamilyPhyPolyHermite, FamilyScaledPhyPolyHermite,
        derivative,
        evaluate!, evaluate,
        vander, vander!

"""
PhyPolyHermite <: Hermite

An immutable structure for physicist Hermite polynomials

## Fields
-  `m` : order of the polynomial
-  `P` : physicist Hermite polynomial of order m
- `scaled` : with rescaling to have unitary norm
"""
struct PhyPolyHermite <: Hermite
    m::Int64
    P::ImmutablePolynomial{Float64}
    scaled::Bool
end

function Base.show(io::IO, P::PhyPolyHermite)
    println(io, string(P.m)*"-th order physicist Hermite polynomial"*string(P.P)*", scaled = "*string(P.scaled))
end
# Hn(x)  = (-1)ⁿ*exp(x²)dⁿ/dxⁿ exp(-x²)
# Hn′(x) = 2n*Hn-1(x)
# Hn″(x) = 2n*(n-1)*Hn-1(x)

# Cphy(m::Int64) =sqrt(gamma(m+1) * 2^m * sqrt(π))
"""
$(TYPEDSIGNATURES)

Normalization constant of the m-th order physicist Hermite polynomial
"""
Cphy(m::Int64) = sqrt(gamma(m+1) * 2^m) # sqrt(π) is already accounted for by sampling

"""
$(TYPEDSIGNATURES)

Normalization constant of the m-th order physicist Hermite polynomial
"""
Cphy(P::PhyPolyHermite) = Cphy(P.m)

"""
$(TYPEDSIGNATURES)

Return the degree of the physicist Hermite function `P`.
"""
degree(P::PhyPolyHermite) = P.m

# Adapted https://people.sc.fsu.edu/~jburkardt/m_src/hermite_polynomial/h_polynomial_coefficients.m

#    Input, integer N, the highest order polynomial to compute.
#    Note that polynomials 0 through N will be computed.
#
#    Output, real C(1:N+1,1:N+1), the coefficients of the Hermite
#    polynomials.

# Recurrence relation:
# Hn+1 = 2*x*Hn - 2*n*Hn-1

"""
$(TYPEDSIGNATURES)

Computes the coefficients of the physicist Hermite polynomials up to order m.
"""
function phyhermite_coeffmatrix(m::Int64)

    if m < 0
        return Float64[]
    end

    coeff = zeros(m+1, m+1)
    coeff[1, 1] = 1.0

    if m == 0
        return coeff
    end

    coeff[2, 2] = 2.0

    if m==1
        return coeff
    end

    for i  = 2:m
        coeff[i+1, 1]      = - 2.0 * (i - 1) * coeff[i-1, 1]
        coeff[i+1, 2:i-1]  =   2.0 *           coeff[i  , 1:i-2] -
                               2.0 * (i - 1) * coeff[i-1, 2:i-1]
        coeff[i+1, i]      =   2.0 *           coeff[i  , i-1]
        coeff[i+1, i+1]    =   2.0 *           coeff[i  , i]
    end
    return coeff
end

const PhyPolyH = phyhermite_coeffmatrix(CstMaxDegree)

function PhyPolyHermite(m::Int64;scaled::Bool= false)
    @assert m>=0 "The order of the polynomial should be >=0"
    if scaled ==false
            return PhyPolyHermite(m, ImmutablePolynomial(view(PhyPolyH,m+1,1:m+1)), scaled)
    else
        C = 1/Cphy(m)
            return PhyPolyHermite(m, ImmutablePolynomial(C*view(PhyPolyH,m+1,1:m+1)), scaled)
    end
end

"""
$(TYPEDSIGNATURES)

Evaluates the physicist Hermite polynomial `P`.
"""
(P::PhyPolyHermite)(x) = P.P(x)

const FamilyPhyPolyHermite = ntuple(i->PhyPolyHermite(i-1), CstMaxDegree+1)
const FamilyScaledPhyPolyHermite = ntuple(i->PhyPolyHermite(i-1; scaled = true), CstMaxDegree+1)

# Compute the k-th derivative of a physicist Hermite polynomial according to
# H_{n}^(k)(x) = 2^{k} n!/(n-k)! H_{n-k}(x)
"""
$(TYPEDSIGNATURES)

Computes the k-th derivative of the physicist Hermite polynomial `P`.
"""
function derivative(P::PhyPolyHermite, k::Int64)
    m = P.m
    @assert k>=0 "This function doesn't compute anti-derivatives of Hermite polynomials"
    if m>=k
        factor = 2^k*exp(loggamma(m+1) - loggamma(m+1-k))
        if P.scaled == false
            return ImmutablePolynomial(factor*PhyPolyH[m-k+1,1:m-k+1])
        else
            C = 1/Cphy(m)
            return ImmutablePolynomial(factor*C*PhyPolyH[m-k+1,1:m-k+1])
        end
    else
        return ImmutablePolynomial((0.0))
    end
end

"""
$(TYPEDSIGNATURES)

Evaluates in-place the physicist Hermite polynomial `P`.
"""
function evaluate!(dV, P::PhyPolyHermite, x)
    m = P.m
    N = size(x,1)
    @assert size(dV) == (N, m+1) "Wrong dimension of the Vander matrix"

    # H_0(x) = 1, H_1(x) = 2*x
    col = view(dV,:,1)
    @avx @. col = 1.0

    if P.scaled
        rmul!(col, 1/Cphy(0))
    end
    if m == 0
        return dV
    end

    col = view(dV,:,2)
    @avx @. col = 2.0*x

    if P.scaled
        rmul!(col, 1/Cphy(1))
    end

    if m == 1
        return dV
    end

    if P.scaled
        @inbounds for i=2:m
            colp1 = view(dV,:,i+1)
            col   = view(dV,:,i)
            colm1 = view(dV,:,i-1)
            @avx @. colp1 = 2.0*(Cphy(i-1)* x * col - Cphy(i-2)*(i-1)*colm1)
            # dV[:,i+1] = 2.0*Cphy(i-1)*x .* dV[:,i] - 2.0*Cphy(i-2)*i*dV[:,i-1]
            rmul!(colp1, 1/Cphy(i))
        end
    else
        @inbounds for i=2:m
            colp1 = view(dV,:,i+1)
            col   = view(dV,:,i)
            colm1 = view(dV,:,i-1)
            @avx @. colp1 = 2.0*(x * col - (i-1)*colm1)
        end
    end
    return dV
end

"""
$(TYPEDSIGNATURES)

Evaluates the physicist Hermite function `P`.
"""
evaluate(P::PhyPolyHermite, x::Array{Float64,1}) = evaluate!(zeros(size(x,1), P.m+1), P, x)

# H_{n}^(k)(x) = 2^{k} n!/(n-k)! H_{n-k}(x)

# vander!(dV::Array{Float64,2}, P::PhyPolyHermite{m}, k=0, x) where {m} = evaluate!(dV::Array{Float64,2}, P::PhyPolyHermite{m}, x)
"""
$(TYPEDSIGNATURES)

Computes in-place the Vandermonde matrix of the k-th derivative of the physicist Hermite polynomial `P`.
"""
function vander!(dV, P::PhyPolyHermite, k::Int64, x)
    m = P.m

    if k==0
        evaluate!(dV, P, x)
    # Derivative
    elseif k==1
        # Use recurrence relation H′_{n+1} = (1 + 1/n)*(2*x*H′_{n} - 2*n*H′_{n-1})
        N = size(x,1)
        @assert size(dV) == (N, m+1) "Wrong dimension of the Vander matrix"

        col = view(dV,:,1)
        @avx @. col = 0.0
        if m==0
            return dV
        end

        col = view(dV,:,2)
        @avx @. col = 2.0

        if m==1
            return dV
        end

        col = view(dV,:,3)
        @avx @. col = 8.0*x

        if m==2
            return dV
        end

        @inbounds for i=3:m
            colp1 = view(dV,:,i+1)
            col   = view(dV,:,i)
            colm1 = view(dV,:,i-1)

            @avx @. colp1 = 2.0*(1.0 + 1.0/(i-1.0))*(x * col - (i-1.0) * colm1)
        end

        if P.scaled
            @inbounds for i=1:m
                colp1 = view(dV,:,i+1)
                factor = 2^k*exp(loggamma(i+1) - loggamma(i+1-k))
                rmul!(colp1, 1/(Cphy(i)*sqrt(factor)))
            end

        end
        return dV
    # Second Derivative
    elseif k==2
        # Use recurrence relation H′_{n+1} = (1 + 1/n)*(2*x*H′_{n} - 2*n*H′_{n-1})
        N = size(x,1)
        @assert size(dV) == (N, m+1) "Wrong dimension of the Vander matrix"

        col = view(dV,:,1)
        @avx @. col = 0.0
        if m==0
            return dV
        end

        col = view(dV,:,2)
        @avx @. col = 0.0

        if m==1
            return dV
        end

        col = view(dV,:,3)
        @avx @. col = 8.0

        if m==2
            return dV
        end

        col = view(dV,:,4)
        @avx @. col = 48.0 * x

        if m==3
            return dV
        end

        @inbounds for i=4:m
            colp1 = view(dV,:,i+1)
            col   = view(dV,:,i)
            colm1 = view(dV,:,i-1)

            @avx @. colp1 = 2.0*((i)/(i-2.0))*(x * col - (i-1.0) * colm1)
        end

        if P.scaled
            @inbounds for i=1:m
                colp1 = view(dV,:,i+1)
                factor = 2^k*exp(loggamma(i+1.0) - loggamma(i+1.0-k))
                @avx rmul!(colp1, 1/(Cphy(i)*sqrt(factor)))
            end

        end
        return dV
    end
end

"""
$(TYPEDSIGNATURES)

Computes the Vandermonde matrix of the k-th derivative of the physicist Hermite polynomial `P`.
"""
vander(P::PhyPolyHermite, k::Int64, x::Array{Float64,1}) = vander!(zeros(size(x,1), P.m+1), P, k, x)


text:
exponentialsample(λ) = -log(rand()) / λ

function denseperiodic(f, dur, fs; magnitude=20.0)
	N_pat = Int(round(dur * fs))
	dn_pat = Int(round(fs/f))
	train = zeros(N_pat)
	idx = 1
	while idx <= N_pat
		train[idx] = magnitude
		idx += dn_pat
	end
	train
end

function sparseperiodic(f, dur)
	numspikes = Int(floor(dur*f))
	dt = 1.0/f
	[dt*(k-1) for k ∈ 1:numspikes]
end

function sparsepoisson(λ, dur)
	ts = Array{Float64, 1}(undef, 0)
	t = 0.0
	while t < dur
		push!(ts, t)
		t += exponentialsample(λ)
	end
	ts
end

function sparserefractorypoisson(λ, dur, deaddur)
	ts = Array{Float64, 1}(undef, 0)
	t = 0.0
	while t < dur
		push!(ts, t)
		sample = exponentialsample(λ)
		while sample < deaddur
			sample = exponentialsample(λ)
		end
		t += sample
	end
	ts
end

function densepoisson(λ, dur, fs; magnitude=20.0)
	ts = sparsepoisson(λ, dur)
	output = zeros(Int64(round(dur*fs))+1)
	for t ∈ ts
		output[Int64(round(t*fs))+1] = magnitude
	end
	output
end

function channeldup(times::Array{T,1}, channels) where T <: Number
	numspikes = length(times)
	numchan = length(channels)
	spikes = Array{Spike,1}(undef,numspikes*numchan)
	count = 1
	for i ∈ 1:numspikes
		for channel ∈ channels
			spikes[count] = Spike(channel, times[i])
			count += 1
		end
	end
	spikes
end

function channelscatter(times::Array{T,1}, available_channels) where T <: Number
	n = length(times)
	spikes = Array{Spike, 1}(undef, n)
	if length(available_channels) > n
		channels = shuffle(available_channels)[1:n]
		for i ∈ 1:n
			spikes[i] = Spike(channels[i], times[i])
		end
	else
		for i ∈ 1:n
			spikes[i] = Spike(rand(available_channels), times[i])
		end
	end
	spikes
end

function stochasticblock(times::Array{T,1}, p, channels) where T <: Number
	n = length(times)
	spikes = Array{Spike,1}(undef,0)
	for time ∈ times
		for channel ∈ channels
			if rand() < p
				push!(spikes, Spike(channel, time))
			end
		end
	end
	spikes
end

function genrandinput(λ, dur, nchannels)
	times = sparsepoisson(λ, dur)
	channelscatter(times, 1:nchannels)
end

function gendupinput(λ, dur, channels)
	times = sparsepoisson(λ, dur)
	channeldup(times, channels)
end

max_stars_repo_name: analogouscircuit/DelayNetTools

max_issues_repo_name: analogouscircuit/DelayNetTools

max_forks_repo_name: analogouscircuit/DelayNetTools


text:
module GPUArrays

using Serialization
using Random
using LinearAlgebra
using Printf

using LinearAlgebra.BLAS
using Base.Cartesian

using AbstractFFTs

using Adapt

# device functionality
include("device/device.jl")
include("device/execution.jl")
## executed on-device
include("device/abstractarray.jl")
include("device/indexing.jl")
include("device/memory.jl")
include("device/synchronization.jl")

# host abstractions
include("host/abstractarray.jl")
include("host/construction.jl")
## integrations and specialized methods
include("host/base.jl")
include("host/indexing.jl")
include("host/broadcast.jl")
include("host/mapreduce.jl")
include("host/linalg.jl")
include("host/math.jl")
include("host/random.jl")
include("host/quirks.jl")

# CPU reference implementation
include("reference.jl")

end # module


text:
"""
Quicksort!
Alex Ellison
@xaellison

Usage:
quicksort!(my_cuarray)

The main quicksort kernel uses dynamic parallelism. Let's call blocksize M. The
first part of the kernel bubble sorts M elements with maximal stride between
lo and hi. If the sublist is <= M elements, stride = 1 and no recursion
happens. Otherwise, we pick element lo + M / 2 * stride as a pivot. This is
an efficient choice for random lists and pre-sorted lists.

Partition is done in stages:
1. In a separate kernel: Merge-sort batches of M values using their comparison
    to pivot as a key. The comparison alternates between < and <= with recursion
    depth. This makes no difference when there are many unique values, but when
    there are many duplicates, this effectively partitions into <, =, and >.
2. Consolidate batches. This runs inside the quicksort kernel.

Information about the partitions of each batch is stored in shared memory if
there are fewer than N batches (list length <= M ^ 2). Otherwise, it is stored in
a global array. Any sub-kernels which also need that array will be executed
synchronously so they don't conflict over that memory.

Naming conventions:
Sublists (ranges of the list being sorted) are denoted by `my_floor` and one of
    `L` and `my_ceil`. `my_floor` is an exclusive lower bound, `my_ceil` is an
    inclusive upperboard, `L` is their difference.

`lo` / `hi` are used in place of `my_floor` / `my_ceil` in `qsort_kernel`
    `lo` / `hi` are indices of  batches in `consolidate_batch_partition`

`b_sums` is "batch sums", the number of values in a batch which are >= pivot
"""

#-------------------------------------------------------------------------------
# Integer arithmetic
"""
Equivalent to Int(ceil(a / b))
"""
function ceil_div(a, b)
    return a ÷ b + 1 & (a % b != 0)
end

"""
Returns smallest power of 2 >= x
"""
function pow2_ceil(x)
    out = 1
    while out < x
        out *= 2
    end
    out
end

function pow2_floor(x)
    out = 1
    while out * 2 <= x
        out *= 2
    end
    out
end

"""
For a batch of size `n` what is the lowest index of the batch `i` is in
"""
function batch_floor(idx, n)
    return idx - (idx - 1) % n
end

"""
For a batch of size `n` what is the highest index of the batch `i` is in
"""
function batch_ceil(idx, n)
    return idx + n - 1 - (idx - 1) % n
end

"""
GPU friendly step function (step at i = 1)
"""
function Θ(i)
    return 1 & (1 <= i)
end

"""
Suppose we are merging two lists of size n, each of which has all falses before
all trues. Together, they will be indexed 1:2n. This is a fast stepwise function
for the destination index of a value at index `x` in the concatenated input,
where `a` is the number of falses in the first half, b = n - a, and false is the
number of falses in the second half.
"""
function step_swap(x, a, b, c)
    return x + Θ(x - a) * b - Θ(x - (a + c)) * (b + c) + Θ(x - (a + b + c)) * c
end

"""
Generalizes `step_swap` for when the floor index is not 1
"""
function batch_step_swap(x, n, a, b, c)
    idx = (x - 1) % n + 1
    return batch_floor(x, n) - 1 + step_swap(idx, a, b, c)
end

@inline function flex_lt(a, b, eq :: Bool)
    if eq
        return a >= b
    else
        return a > b
    end
end

#-------------------------------------------------------------------------------
# Batch partitioning
"""
For thread `idx` with current value `value`, merge two batches of size `n` and
return the new value this thread takes. `sums` and `swap` are shared mem
"""
function merge_swap_shmem(value, idx, n, sums, swap)
    @inbounds begin
    sync_threads()
    b = sums[batch_floor(idx, 2 * n)]
    a = n - b
    d = sums[batch_ceil(idx, 2 * n)]
    c = n - d
    swap[idx] = value
    sync_threads()
    sums[idx] = d + b
    return swap[batch_step_swap(idx, 2 * n, a, b, c)]
    end
end

@inline function batch_partition(values, swap, pivot, sums, i, idx0, my_floor, L, N_t, parity)
    if idx0 - my_floor <= L
        @inbounds swap[i] = values[idx0]
        @inbounds sums[i] = 1 & flex_lt(swap[i], pivot, parity)
    else
        @inbounds sums[i] = 1
    end
    sync_threads()
    val = merge_swap_shmem(swap[i], i, 1, sums, swap)
    temp = 2
    while temp < N_t
        val = merge_swap_shmem(val, i, temp, sums, swap)
        temp *= 2
    end
    sync_threads()

    if idx0 - my_floor <= L
        @inbounds values[idx0] = val
    end
end

"""
Partition `values` such that all < `pivot` come first, then all >= `pivot`.
NOTE: It is very important that this is a stable sort. We effectively allow
    unallocated memory to be our null values, and set them to true, at the top
    of the list. That way, they can be operated on, but at the end they aren't
    transferred out to global mem.
"""
function partition_batches_kernel(values :: CuDeviceArray{T}, pivot, my_floor, L, parity :: Val{P}) where {T, P}
    N_t = blockDim().x
    i = threadIdx().x
    idx0 = my_floor + (blockIdx().x - 1) * blockDim().x + i # for addressing
    swap = @cuDynamicSharedMem(T, N_t, 4 * N_t)
    sums = @cuDynamicSharedMem(Int32, N_t)
    batch_partition(values, swap, pivot, sums, i, idx0, my_floor, L, N_t, P)
    return nothing
end

"""
Partition `values` such that all < `pivot` come first, then all >= `pivot`.
Method stores the number of values >= `pivot` in each batch in `batch_sums`
"""
function partition_batches_kernel(values :: CuDeviceArray{T}, pivot, batch_sums, my_floor, L, parity :: Val{P}) where {T, P}
    N_t = blockDim().x
    i = threadIdx().x
    idx0 = my_floor + (blockIdx().x - 1) * blockDim().x + i # for addressing
    swap = @cuDynamicSharedMem(T, N_t, 4 * N_t)
    sums = @cuDynamicSharedMem(Int32, N_t)
    batch_partition(values, swap, pivot, sums, i, idx0, my_floor, L, N_t, P)

    if i == 1
        # this is the amount by which our kernel "overhangs" the real list
        # and how many unnecessary trues are included in the sum
        overhang = (idx0 - my_floor <= (L - L % N_t) || L % N_t == 0 ? 0 : N_t - (L % N_t))
        @inbounds batch_sums[blockIdx().x] = sums[i] - overhang
    end
    return nothing
end

#-------------------------------------------------------------------------------
# Batch consolidation
"""
Finds the index in `array` of the last value less than `pivot`.
Searches after index `my_floor` up to (inclusive) index `my_ceil`
"""
function seek_last_less_than(array, pivot, my_floor, my_ceil, parity) :: Int32
    low = my_floor + 1
    high = my_ceil
    while low <= high
        mid = (low + high) ÷ 2
        if flex_lt(array[mid], pivot, parity)
            high = mid - 1
        else
            low = mid + 1
        end
    end
    return low - 1
end

"""
Performs binary search on batches to find their partitions. Stores result in
`b_sums`
"""
@inline function init_bsums(vals, pivot, my_floor, L, N_t, N_b, i, b_sums, parity)
    if i <= N_b
        seek_lo = my_floor + (i - 1) * N_t
        seek_hi = my_floor + min(Int32(L), Int32(i * N_t))
        b_sums[i] = seek_hi - seek_last_less_than(vals, pivot, seek_lo, seek_hi, parity)
    end
    sync_threads()
end

"""
In linear time, merge `N_b` consecutive, partitioned batches so that the
`L` length region of `vals` after `lo` is partitioned. Must only run on 1 SM.
`i` - thread id
`b_sums` - number of values >= `pivot` in each batch
"""
@inline function consolidate_batch_partition(vals, pivot, my_floor, L, N_t, N_b, i, b_sums)
    i = threadIdx().x
    sync_threads()
    b = b_sums[1]
    a = N_t - b
    sync_threads()
    for hi in 2:N_b
        n_eff =  (hi != N_b || L % N_t == 0) ? N_t : L % N_t
        sync_threads()
        d = b_sums[hi]
        c = n_eff - d
        to_move = min(b, c)
        sync_threads()
        if i <= n_eff && i <= to_move
            swap = vals[my_floor + a + i]
        end
        sync_threads()
        if i <= n_eff && i <= to_move
            vals[my_floor + a + i] = vals[my_floor + a + b + c - to_move + i]
        end
        sync_threads()
        if i <= n_eff && i <= to_move
            vals[my_floor + a + b + c - to_move + i] = swap
        end
        sync_threads()
        a += c
        b += d
    end
    sync_threads()
    return my_floor + a
end

#-------------------------------------------------------------------------------
# Sorting
"""
Performs bubble sort on `vals` starting after `lo` and going for min(`L`, 1024)
elements spaced by `stride`. Good for sampling pivot values as well as short
sorts.
"""
@inline function bubble_sort(vals, swap, my_floor, L, i, stride, N_t)
    if i <= L
        L = min(Int32(N_t), Int32(L))
        @inbounds begin
        swap[i] = vals[my_floor + i * stride]
        for level in 0:L
            # get left/right neighbor depending on even/odd level
            buddy = i - 1 + 2 * (1 & (i % 2 != level % 2))
            if 1 <= buddy <= L
                buddy_val = swap[buddy]
            end
            sync_threads()
            if 1 <= buddy <= L
                if (i < buddy) != (swap[i] < buddy_val)
                    swap[i] = buddy_val
                end
            end
            sync_threads()
        end
        vals[my_floor + i * stride] = swap[i]
        end
    end
end

"""
Perform quicksort on `vals` for the region with `lo` as an exclusive floor and
`hi` as an inclusive ceiling.
`block_dim` must be a Val{N_t} where `N_t` must equal the number of threads this
    kernel is launched with.

Assumes that RUNTIME_SYNC_DEPTH = 24. Recursion ends when `depth` = 24, or when
`hi` - `lo` <= `N_t`. If `depth` = 24 but `hi` - `lo` > `N_t`, sort may not be
complete, but further recursion might involve invalid memory accesses.
"""
function qsort_kernel(vals :: CuDeviceArray{T}, lo, hi, global_b_sums, parity :: Val{P}, depth :: Int8) where {T, P}
    #N_t = blockDim().x
    b_sums = @cuDynamicSharedMem(Int32, blockDim().x)
    swap = @cuDynamicSharedMem(T, blockDim().x, 4 * blockDim().x)
    i = threadIdx().x
    L = Int32(hi - lo)
    N_b = ceil_div(L, blockDim().x)

    #= step 1 bubble sort. It'll either finish sorting a subproblem or help
    select a pivot value =#
    bubble_sort(vals, swap, lo, L, i, L <= blockDim().x ? Int32(1) : Int32(L ÷ blockDim().x), blockDim().x)
    if L <= blockDim().x || depth >= 24
        return
    end

    @inbounds pivot = vals[lo + (blockDim().x ÷ 2) * (L ÷ blockDim().x)]

    # step 2: use pivot to partition into batches
    if i == 1
        if L <= blockDim().x ^ 2
            # perform batch partition and do not store partition indices
            @cuda blocks=N_b threads=blockDim().x dynamic=true shmem=blockDim().x*(4+sizeof(T)) partition_batches_kernel(vals, pivot, lo, L, parity)
        else
            # store batch partition indices in global array
            @cuda blocks=N_b threads=blockDim().x dynamic=true shmem=blockDim().x*(4+sizeof(T)) partition_batches_kernel(vals, pivot, global_b_sums, lo, L, parity)
        end
    end
    CUDA.device_synchronize()
    sync_threads()
    #= step 3: consolidate the partitioned batches so that the sublist from
    [lo, hi) is partitioned, and the partition is stored in `partition`=#

    partition = -1
    if L <= blockDim().x ^ 2
        #= if small enough, perform binary search to find the partitions in each
        batch, store in shmem, and use that=#
        init_bsums(vals, pivot, lo, L, blockDim().x, N_b, i, b_sums, P)
        partition = consolidate_batch_partition(vals, pivot, lo, L, blockDim().x, N_b, i, b_sums)
    else
        #= if more batches than threads per SM, use the global array written to
        by `partition_batches_kernel`=#
        partition = consolidate_batch_partition(vals, pivot, lo, L, blockDim().x, N_b, i, global_b_sums)
    end

    #= step 4: recursion. If the first subproblem doesn't require global mem,
    run it on a temporary stream. Otherwise run in default stream, so that
    the two subproblems run serially, and don't conflict over the global array=#
    if i == 1 && partition > lo
        if partition - lo <= blockDim().x ^ 2
            s = CuDeviceStream()
            @cuda threads=blockDim().x dynamic=true stream=s shmem=blockDim().x*(4+sizeof(T)) qsort_kernel(vals, lo, Int32(partition), global_b_sums, Val(!P), depth + Int8(1))
            CUDA.unsafe_destroy!(s)
        else
            @cuda threads=blockDim().x dynamic=true shmem=blockDim().x*(4+sizeof(T)) qsort_kernel(vals, lo, Int32(partition), global_b_sums, Val(!P), depth + Int8(1))
        end
    end

    if i == 1 && partition < hi
        if hi - partition <= blockDim().x ^ 2
            s = CuDeviceStream()
            @cuda threads=blockDim().x dynamic=true stream=s shmem=blockDim().x*(4+sizeof(T)) qsort_kernel(vals, Int32(partition), hi, global_b_sums, Val(!P), depth + Int8(1))
            CUDA.unsafe_destroy!(s)
        else
            @cuda threads=blockDim().x dynamic=true shmem=blockDim().x*(4+sizeof(T)) qsort_kernel(vals, Int32(partition), hi, global_b_sums, Val(!P), depth + Int8(1))
        end
    end
    return nothing
end

function quicksort!(c :: CuArray{T}) where T
    CUDA.limit!(CUDA.LIMIT_DEV_RUNTIME_SYNC_DEPTH, 24)
    N = length(c)

    function get_config(kernel)
        get_shmem(threads) = threads * (sizeof(Int32) + max(4, sizeof(T)))
        fun = kernel.fun
        config = launch_configuration(fun, shmem=threads->get_shmem(threads))
        threads = pow2_floor(config.threads)
        @assert threads <= config.threads
        @assert threads >= 32
        return (blocks=1, threads=threads, shmem=get_shmem(threads))
    end

    b_sums = CuArray{Int32}(undef, ceil_div(N, 32))
    CUDA.@sync(@cuda config=get_config qsort_kernel(c, Int32(0), N, b_sums, Val(true), Int8(1)))
end

max_stars_repo_name: StatisticalMice/experiments

max_issues_repo_name: StatisticalMice/experiments

max_forks_repo_name: StatisticalMice/experiments


text:
@testset "plotting" begin

    x = rand(10)
    f = GP(SqExponentialKernel())
    gp = f(x)
    plt1 = sampleplot(gp, 10)
    @test plt1.n == 10

    rec1 = RecipesBase.apply_recipe(Dict{Symbol, Any}(), gp)
    @test rec1[1].args[1] ≈ x atol=1e-5
    @test rec1[1].args[2] ≈ zero(x) atol=1e-5

    rec2 = RecipesBase.apply_recipe(Dict{Symbol, Any}(), f, rand(10))
    @test typeof(rec2[1].args[1]) <: AbstractGPs.FiniteGP

    rec3 = RecipesBase.apply_recipe(Dict{Symbol, Any}(), f, 0:0.01:1)
    @test typeof(rec3[1].args[1]) <: AbstractGPs.FiniteGP

    rec4 = RecipesBase.apply_recipe(Dict{Symbol, Any}(), f, 0, 1)
    @test typeof(rec4[1].args[1]) <: AbstractGPs.FiniteGP
end


text:
function cstrt_bingane_ngon(n)
if mod(log2(n),1)==0 && n>=8
    # initialization
    d = pi/n - asin(sin(2*pi/n)/2);
    # construction
    a = zeros(n-1); b = zeros(n-1);
    for i = 1:Int(n/4)
        a[Int(mod(i*(n/2-1),n))] = sin(2*i*pi/n)*sin(d + (-1)^(i-1)*pi/n)/sin(2*pi/n);
        a[Int(n-mod(i*(n/2-1),n))] = -a[Int(mod(i*(n/2-1),n))];
        b[Int(mod(i*(n/2-1),n))] = (sin(pi/n - d) + cos(2*i*pi/n)*sin(d + (-1)^(i-1)*pi/n))/sin(2*pi/n);
        b[Int(n-mod(i*(n/2-1),n))] = b[Int(mod(i*(n/2-1),n))];
    end
    for i = 1:Int(n/4-1)
        a[Int(mod(i*(n/2-1)+n/2,n))] = a[Int(mod(i*(n/2-1),n))] + (-1)^i*sin(2*i*pi/n);
        a[Int(n-mod(i*(n/2-1)+n/2,n))] = -a[Int(mod(i*(n/2-1)+n/2,n))];
        b[Int(mod(i*(n/2-1)+n/2,n))] = b[Int(mod(i*(n/2-1),n))] + (-1)^i*cos(2*i*pi/n);
        b[Int(n-mod(i*(n/2-1)+n/2,n))] = b[Int(mod(i*(n/2-1)+n/2,n))];
    end
    a[Int(n/2)] = 0; b[Int(n/2)] = 1;
    return(a,b)
end
end

max_stars_repo_name: cbingane/optigon

max_issues_repo_name: cbingane/optigon

max_forks_repo_name: cbingane/optigon


text:
mini = minimum
maxi = maximum

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
using JuliaSets
using Base.Test

# Define function
function R{T <: Real}(a::Complex{T})
    a^2 - 0.75
end

# Run test problem
x = [-2:.01:2]
y = [-2:.01:2]
A = visJuliaSet(R, x, y)

@test A[1, 1] < 10

max_stars_repo_name: Jgmorton/JuliaSets

max_issues_repo_name: Jgmorton/JuliaSets

max_forks_repo_name: Jgmorton/JuliaSets


text:
using MatrixCompletion.Utilities.PrettyPrinter
using Printf

@testset "$(format("PrettyPrinter: Tables Header"))" begin
  let
    header_list = ["Iter", "R(primal)", " R(dual)",  "ℒ(Gaussian)", "ℒ(Bernoulli)", "ℒ(Poisson)", "ℒ(Gamma)", "λ‖diag(Z)‖ᵢ", " μ⟨I, X⟩", " ‖Z₁₂‖ᵢ "]
    row = map(x -> @sprintf("%3.2e", x), rand(10))
    row[1] = "100"
    table_header(header_list)
    add_row(header_list, data = row)
  end
end


text:
# unit test the bGrad

include("../src/DMD_Util.jl");
include("../src/DMD_Type.jl");
include("../src/DMD_Func.jl");

# Generate DMD Synthetic Data
#------------------------------------------------------------------------------
# dimensions
m = 100;    # temporal dimension
n = 100;    # spatial dimension
k = 3;      # number of modes
T = Float64;
# generate data
X, t, at, Bt = simDMD(m, n, k, T; seed=123456);
# specify the loss
function lossFunc(r)
	return sum(abs2, r)
end
#
function lossGrad(r)
	conj!(r)
end
#
# Create object
#------------------------------------------------------------------------------
params = DMDParams(k, X, t, lossFunc, lossGrad);
#
# Initialize with true parameter
#------------------------------------------------------------------------------
copyto!(params.a, at);
copyto!(params.B, Bt);
#
# Evaluate bGrad
#------------------------------------------------------------------------------
Random.randn!(params.B);
#
err = aBFunc(params);
#
if err < 1e-6
	println("aBFunc: OK");
else
	println("aBFunc: Wrong, err: $err");
end
#
id = rand(1:params.n);
#
Random.randn!(params.b[id]);
#
err = abFunc(params, id);
#
if err < 1e-6
	println("abFunc: OK");
else
	println("abFunc: Wrong, err: $err");
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

Main.Core.eval(Main.Core, :(baremodule Inference
using Core.Intrinsics
import Core: print, println, show, write, unsafe_write, STDOUT, STDERR

ccall(:jl_set_istopmod, Void, (Bool,), false)

eval(x) = Core.eval(Inference, x)
eval(m, x) = Core.eval(m, x)

const include = Core.include
# conditional to allow redefining Core.Inference after base exists
isdefined(Main, :Base) || ((::Type{T}){T}(arg) = convert(T, arg)::T)

## Load essential files and libraries
include("essentials.jl")
include("ctypes.jl")
include("generator.jl")
include("reflection.jl")
include("options.jl")

# core operations & types
include("promotion.jl")
include("tuple.jl")
include("traits.jl")
include("range.jl")
include("expr.jl")
include("error.jl")

# core numeric operations & types
include("bool.jl")
include("number.jl")
include("int.jl")
include("operators.jl")
include("pointer.jl")
const checked_add = +
const checked_sub = -

# core array operations
include("array.jl")
include("abstractarray.jl")

include("hashing.jl")
include("nofloat_hashing.jl")

# map-reduce operators
macro simd(forloop)
    esc(forloop)
end
include("reduce.jl")

## core structures
include("bitarray.jl")
include("intset.jl")
include("associative.jl")

# core docsystem
include("docs/core.jl")

# compiler
include("inference.jl")
ccall(:jl_set_typeinf_func, Void, (Any,), typeinf_ext)

end # baremodule Inference
))


text:
###########################################################
##
## MAPS
##
## maps between residue rings of polynomials and in general
##
###########################################################

function domain(M::Map(HeckeMap))
  return M.header.domain
end

function codomain(M::Map(HeckeMap))
  return M.header.codomain
end

function image_function(f::Map(HeckeMap))
  if isdefined(f.header, :image)
    return f.header.image
  else
    return x -> image(f, x)
  end
end

function preimage_function(f::Map(HeckeMap))
  if isdefined(f.header, :preimage)
    return f.header.preimage
  else
    return x -> preimage(f, x)
  end
end

export Map

mutable struct MapCache{D, C, De, Ce}
  lim::Int

  im::Dict{De, Ce}
  imStat::Dict{De, Int}

  pr::Dict{Ce, De}
  prStat::Dict{Ce, Int}

  old_im::Function
  old_pr::Function

  function MapCache(dom::D, cod::C, ::Type{De}, ::Type{Ce}, lim::Int=100) where {D, C, De, Ce}
    r = new{D, C, De, Ce}()
    r.lim = lim
    r.im = Dict{De, Ce}()
    r.imStat = Dict{De, Int}()
    r.pr = Dict{Ce, De}()
    r.prStat = Dict{Ce, Int}()
    return r
  end
end

mutable struct MapHeader{D, C}
  domain::D
  codomain::C
  image::Function
  preimage::Function
  cache::MapCache
  @declare_other

  function MapHeader{D, C}() where {D, C}
    z = new{D, C}()
    return z
  end

  function MapHeader{D, C}(domain::D, codomain::C) where {D, C}
    z = new{D, C}()
    z.domain = domain
    z.codomain = codomain
    return z
  end

  function MapHeader{D, C}(domain::D, codomain::C, image::Function) where {D, C}
    z = new{D, C}()
    z.domain = domain
    z.codomain = codomain
    z.image = image
    return z
  end

  function MapHeader{D, C}(domain::D, codomain::C, image::Function, preimage::Function) where {D, C}
    z = new{D, C}()
    z.domain = domain
    z.codomain = codomain
    z.image = image
    z.preimage = preimage
    return z
  end
end

function MapHeader(domain::D, codomain::C) where {D, C}
  return MapHeader{D, C}(domain, codomain)
end

function MapHeader(domain::D, codomain::C, image::Function) where {D, C}
  return MapHeader{D, C}(domain, codomain, image)
end

function MapHeader(domain::D, codomain::C, image::Function, preimage::Function) where {D, C}
  return MapHeader{D, C}(domain, codomain, image, preimage)
end

# this type represents a -> f(g(a))
#mutable struct CompositeMap{D, C, R} <: Map{D, C, HeckeMap, CompositeMap}
#  header::MapHeader{D, C}
#  f::Map{R, C}
#  g::Map{D, R}
#
#  function CompositeMap{D, C, R}(f::Map, g::Map) where {D, C, R}
#  ##CF should be function CompositeMap(f::Map{R, C}, g::Map{D, R})
#  ## but that seems to not work:
#  # U, m = UnitGroup(ResidueRing(FlintZZ, 2^9));n
#  # q, mq = Hecke.quo(U, [preimage(m, codomain(m)(fmpz(-1)))])
#  # z = Hecke.compose(m, inv(mq))
#  # btw: m*inv(mq) also fails.
#
#
#    domain(f) == codomain(g) || throw("maps not compatible")
#    z = new{D, C, R}()
#    z.f = f
#    z.g = g
#
#    image = function(x)#x::elem_type(domain(z)))
#      parent(x) != domain(g) && error("Element not in domain of map")
#      return f(g(x))::elem_type(codomain(z))
#    end
#
#    if isdefined(f.header, :preimage) && isdefined(g.header, :preimage)
#      _preimage = function(x)#x::elem_type(codomain(z)))
#        return preimage(g, preimage(f, x))::elem_type(domain(z))
#      end
#      z.header = MapHeader(domain(g), codomain(f), image, _preimage)
#    else
#      z.header = MapHeader(domain(g), codomain(f), image)
#    end
#
#    return z
#  end
#end

#function *(f::Map{R, C}, g::Map{D, R}) where {D, C, R}
#  println("======")
#  println(stacktr()[1:2])
#  return CompositeMap{D, C, R}(f, g)
#end

function preimage(f::AbstractAlgebra.Generic.CompositeMap, a)
  return preimage(f.map1, preimage(f.map2, a))
end

preimage_function(f) = a -> preimage(f, a)
image_function(f) = a -> image(f, a)

#function _compose(f, g)
#  return AbstractAlgebra.Generic.compose(g, f)
#end

#function compose(f::Map(HeckeMap){R, C}, g::Map(HeckeMap){D, R}) where {D, C, R}
#  return CompositeMap{D, C, R}(f, g)
#end

mutable struct InverseMap{D, C} <: Map{D, C, HeckeMap, InverseMap}
  header::MapHeader{D, C}
  origin::Map{C, D}

  function InverseMap{D, C}(f::Map{C, D}) where {D, C}
    z = new{D, C}()
    z.header = MapHeader(codomain(f), domain(f), preimage_function(f), image_function(f))
    z.origin = f
    return z
  end
end

function InverseMap(f::Map{C, D}) where {D, C}
  return InverseMap{D, C}(f)
end

###########################################################
# To turn a Function (method) into a map.
###########################################################
mutable struct MapFromFunc{R, T} <: Map{R, T, HeckeMap, MapFromFunc}
  header::Hecke.MapHeader{R, T}
  f::Function
  g::Function

  function MapFromFunc{R, T}(f::Function, D::R, C::T) where {R, T}
    n = new{R, T}()
    n.header = Hecke.MapHeader(D, C, x-> f(x))
    n.f = f
    return n
  end

  function MapFromFunc{R, T}(f::Function, g::Function, D::R, C::T) where {R, T}
    n = new{R, T}()
    n.header = Hecke.MapHeader(D, C, x-> f(x), y->g(y))
    n.f = f
    n.g = g
    return n
  end
end

function Base.show(io::IO, M::MapFromFunc)
  @show_name(io, M)

  io = IOContext(io, :compact => true)
#  println(io, "Map from the $(M.f) julia-function")
  println(io, "Map from")
  show(io, domain(M))
  print(io, " to ")
  show(io, codomain(M))
  print(io, " defined by a julia-function")
  if isdefined(M, :g)
#    println(io, "with inverse by $(M.g)")
    print(io, " with inverse")
  end
end

function MapFromFunc(f::Function, D, C)
  return MapFromFunc{typeof(D), typeof(C)}(f, D, C)
end

function MapFromFunc(f::Function, g::Function, D, C)
  return MapFromFunc{typeof(D), typeof(C)}(f, g, D, C)
end

function Base.inv(M::MapFromFunc)
  if isdefined(M, :g)
     return MapFromFunc(M.g, M.f, codomain(M), domain(M))
  else
     return MapFromFunc(x->preimage(M, x), codomain(M), domain(M))
  end
end

export MapFromFunc


text:
using DynamicalBilliards, Plots
include("../billiard_utils.jl")
include("../plot_utils.jl")
include("../randomfinitewall.jl")
include("../simulation.jl")
include("../vgate.jl")

include("agent.jl")
include("setup.jl")

bt = construct_billiard(poly, sp, ep, clr)
u = construct_controlmap(length(clr))

agt = create_agents(num_agents, length(clr)+1, init_width, init_center, u)

# Loop manually from here if you want to see step by step plots
# for j = 1:1000
    simulate_once!(agt, bt, 0.1)
# end

plot_map(agt, bt)

max_stars_repo_name: gbarsih/Representations-in-Robotics

max_issues_repo_name: gbarsih/Representations-in-Robotics

max_forks_repo_name: gbarsih/Representations-in-Robotics


text:
import DelimitedFiles
import OrderedCollections

function madsmodelrun_adm1(mddata::AbstractDict)
    # function runs the ADM1F model and return selected values from the last result file as predictions

    function get_last_indicator()
        # finds the last indicator file in the current folder
        flist = readdir()
        last_indicator = ""
        for i = 1:length(flist)
            if length(flist[i]) > 9
                if flist[i][1:9] == "indicator"
                    last_indicator = flist[i]
                end
            end
        end
        return last_indicator
    end

    #reads .dat and saves it in mddata
    DelimitedFiles.writedlm("params_tim.dat", Mads.getparamsinit(mddata))

    #removes all previous .out files
    foreach(rm, filter(endswith(".out"), readdir()))

    #runs the ADM1F model with the command line below
    #change the path to the executable below
    command_line=`/Users/elchin/project/ADM1F_WM/build/adm1f -steady -ts_monitor
             -influent_file influent_tim.dat -params_file params_tim.dat
             -ts_type beuler -ts_adapt_type basic -snes_rtol 1.e-5
             -ts_max_snes_failures -1 -Vliq 4981 -t_resx 240`
    Mads.runcmd(command_line; quiet=true, pipe=true)

    #load the outputs from last indicator
    #we are interested only in 7 output values out of 67
    #we return those 7 values as predictions
    last_indicator = get_last_indicator()
    if last_indicator != ""
        o = DelimitedFiles.readdlm(last_indicator)[3:end,:]
        outputs = OrderedCollections.OrderedDict{String, Float64}("o$i"=>o[i] for i=1:67)
        sCOD=0
        for i = 1:7
            sCOD=sCOD+outputs["o$i"]
        end
        Acetate=outputs["o7"]
        Propionate = outputs["o6"]
        Butyrate = outputs["o5"]
        Valerate = outputs["o4"]
        ph = outputs["o27"]
        gaseous = outputs["o41"]*outputs["o44"]*1000
    else
        @error "Output is missing!"
    end
    p=Dict([("o1", sCOD),
            ("o2", Acetate),
            ("o3", Propionate),
            ("o4", Butyrate),
            ("o5", Valerate),
            ("o6", ph),
            ("o7", gaseous)])
    #println("p:",p)
    return p
end