text:
import HTTP
import Markdown: htmlesc
import UUIDs: UUID

# Serve everything from `/frontend`, and create HTTP endpoints to open notebooks.

"Attempts to find the MIME pair corresponding to the extension of a filename. Defaults to `text/plain`."
function mime_fromfilename(filename)
    # This bad boy is from: https://developer.mozilla.org/en-US/docs/Web/HTTP/Basics_of_HTTP/MIME_types/Common_types
    mimepairs = Dict(".aac" => "audio/aac", ".bin" => "application/octet-stream", ".bmp" => "image/bmp", ".css" => "text/css", ".csv" => "text/csv", ".eot" => "application/vnd.ms-fontobject", ".gz" => "application/gzip", ".gif" => "image/gif", ".htm" => "text/html", ".html" => "text/html", ".ico" => "image/vnd.microsoft.icon", ".jpeg" => "image/jpeg", ".jpg" => "image/jpeg", ".js" => "text/javascript", ".json" => "application/json", ".jsonld" => "application/ld+json", ".mjs" => "text/javascript", ".mp3" => "audio/mpeg", ".mp4" => "video/mp4", ".mpeg" => "video/mpeg", ".oga" => "audio/ogg", ".ogv" => "video/ogg", ".ogx" => "application/ogg", ".opus" => "audio/opus", ".otf" => "font/otf", ".png" => "image/png", ".pdf" => "application/pdf", ".rtf" => "application/rtf", ".sh" => "application/x-sh", ".svg" => "image/svg+xml", ".tar" => "application/x-tar", ".tif" => "image/tiff", ".tiff" => "image/tiff", ".ttf" => "font/ttf", ".txt" => "text/plain", ".wav" => "audio/wav", ".weba" => "audio/webm", ".webm" => "video/webm", ".webp" => "image/webp", ".woff" => "font/woff", ".woff2" => "font/woff2", ".xhtml" => "application/xhtml+xml", ".xml" => "application/xml", ".xul" => "application/vnd.mozilla.xul+xml", ".zip" => "application/zip")
    file_extension = getkey(mimepairs, '.' * split(filename, '.')[end], ".txt")
    MIME(mimepairs[file_extension])
end

function asset_response(path)
    if !isfile(path) && !endswith(path, ".html")
        return asset_response(path * ".html")
    end
    try
        @assert isfile(path)
        response = HTTP.Response(200, read(path, String))
        push!(response.headers, "Content-Type" => string(mime_fromfilename(path)))
        push!(response.headers, "Access-Control-Allow-Origin" => "*")
        response
    catch e
        HTTP.Response(404, "Not found!")
    end
end

function error_response(status_code::Integer, title, advice, body="")
    template = read(project_relative_path("frontend", "error.jl.html"), String)

    body_title = body == "" ? "" : "Error message:"
    filled_in = replace(replace(replace(replace(template, "\$TITLE" => title), "\$ADVICE" => advice), "\$BODYTITLE" => body_title), "\$BODY" => htmlesc(body))

    response = HTTP.Response(status_code, filled_in)
    push!(response.headers, "Content-Type" => string(mime_fromfilename(".html")))
    response
end

function notebook_redirect_response(notebook; home_url="./")
    response = HTTP.Response(302, "")
    push!(response.headers, "Location" => home_url * "edit?id=" * string(notebook.notebook_id))
    return response
end

function http_router_for(session::ServerSession)
    router = HTTP.Router()

    function create_serve_onefile(path)
        return request::HTTP.Request -> asset_response(normpath(path))
    end

    HTTP.@register(router, "GET", "/", create_serve_onefile(project_relative_path("frontend", "index.html")))
    HTTP.@register(router, "GET", "/edit", create_serve_onefile(project_relative_path("frontend", "editor.html")))

    HTTP.@register(router, "GET", "/ping", r -> HTTP.Response(200, "OK!"))
    HTTP.@register(router, "GET", "/websocket_url_please", r -> HTTP.Response(200, string(session.secret)))
    HTTP.@register(router, "GET", "/possible_binder_token_please", r -> HTTP.Response(200, session.binder_token === nothing ? "" : session.binder_token))
    HTTP.@register(router, "GET", "/favicon.ico", create_serve_onefile(project_relative_path("frontend", "img", "favicon.ico")))

    function try_launch_notebook_response(action::Function, path_or_url::AbstractString; title="", advice="", home_url="./")
        try
            nb = action(session, path_or_url)
            notebook_redirect_response(nb; home_url=home_url)
        catch e
            if e isa SessionActions.NotebookIsRunningException
                notebook_redirect_response(e.notebook; home_url=home_url)
            else
                error_response(500, title, advice, sprint(showerror, e, stacktrace(catch_backtrace())))
            end
        end
    end

    function serve_newfile(req::HTTP.Request)
        return notebook_redirect_response(SessionActions.new(session))
    end
    HTTP.@register(router, "GET", "/new", serve_newfile)

    function serve_openfile(req::HTTP.Request)
        uri = HTTP.URI(req.target)
        try
            query = HTTP.queryparams(uri)

            if session.options.security.require_token_for_open_links && (UUID(get(query, "secret", string(uuid1()))) != session.secret)
                error_response(405, "Functionality disabled", "This Pluto server does not allow the requested action. If you are running the server yourself, have a look at the security keyword argument to Pluto.run. Please report this error if you did not expect it!")
            else
                if haskey(query, "path")
                    path = tamepath(query["path"])
                    if isfile(path)
                        return try_launch_notebook_response(SessionActions.open, path, title="Failed to load notebook", advice="The file $(htmlesc(path)) could not be loaded. Please report this error!")
                    else
                        return error_response(404, "Can't find a file here", "Please check whether $(htmlesc(path)) exists.")
                    end
                elseif haskey(query, "url")
                    url = query["url"]
                    return try_launch_notebook_response(SessionActions.open_url, url, title="Failed to load notebook", advice="The notebook from $(htmlesc(url)) could not be loaded. Please report this error!")
                else
                    error("Empty request")
                end
            end
        catch e
            return error_response(400, "Bad query", "Please report this error!", sprint(showerror, e, stacktrace(catch_backtrace())))
        end
    end

    HTTP.@register(router, "GET", "/open", serve_openfile)

    function serve_sample(req::HTTP.Request)
        uri = HTTP.URI(req.target)
        sample_path = HTTP.URIs.unescapeuri(split(uri.path, "sample/")[2])
        sample_path_without_dotjl = "sample " * sample_path[1:end - 3]

        path = numbered_until_new(joinpath(tempdir(), sample_path_without_dotjl))
        readwrite(project_relative_path("sample", sample_path), path)

        return try_launch_notebook_response(SessionActions.open, path, home_url="../", title="Failed to load sample", advice="Please report this error!")
    end
    HTTP.@register(router, "GET", "/sample/*", serve_sample)

    function serve_notebookfile(req::HTTP.Request)
        uri = HTTP.URI(req.target)
        try
            query = HTTP.queryparams(uri)
            id = UUID(query["id"])
            notebook = session.notebooks[id]

            response = HTTP.Response(200, sprint(save_notebook, notebook))
            push!(response.headers, "Content-Type" => "text/plain; charset=utf-8")
            push!(response.headers, "Content-Disposition" => "inline; filename=\"$(basename(notebook.path))\"")
            response
        catch e
            return error_response(400, "Bad query", "Please report this error!", sprint(showerror, e, stacktrace(catch_backtrace())))
        end
    end
    HTTP.@register(router, "GET", "/notebookfile", serve_notebookfile)

    function serve_asset(req::HTTP.Request)
        reqURI = req.target |> HTTP.URIs.unescapeuri |> HTTP.URI

        filepath = project_relative_path("frontend", relpath(reqURI.path, "/"))
        asset_response(filepath)
    end
    HTTP.@register(router, "GET", "/*", serve_asset)

    return router
end


text:
################################################################################
#                              DocStringExtension                              #
################################################################################

struct DocThemer <: DocStringExtensions.Abbreviation end

const ATTRIBUTES = DocThemer()

function DocStringExtensions.format(::DocThemer, buf, doc)
	local binding = doc.data[:binding] |> Docs.resolve
	help_attributes(buf, binding; extended=true)
end

ext: jl


text:
abstract type ArrayBackend end

abstract type BackendArrayFunction end
abstract type BackendValue <: ArrayExpr end
abstract type BackendArray <: BackendValue end
abstract type BackendRandomState <: BackendValue end

default_array_backend() = JuliaBackend()

compile_array_function(afunc::ArrayFunction) =
    compile_array_function(default_array_backend(), afunc)
random_state(; seed=nothing) = random_state(default_array_backend(), seed=seed)
empty_array(dtype, shape...) = empty_array(default_array_backend(), dtype, shape...)
to_backend(arr::AbstractArray) = to_backend(default_array_backend(), arr)
from_backend(arr::BackendArray) = error("from_backend() not implemented for $(typeof(arr))")

expr_type(arr::BackendArray) = ArrayType([Nullable{Int64}(x) for x in size(arr)], eltype(arr))

#to_backend(backend, serialized_rs::SerializedRandomState)
#from_backend(rs::BackendSideRandomState)


text:
import FunctionWrappers.FunctionWrapper
using BenchmarkTools
using MittagLeffler
using DSP

#############
# Motivation:
# test best possible functionwrappers and tuple implementation
# with realistic cost functions

#########################
# function wrapper set-up

struct DummyModelFW
  params::Vector{Symbol}
  Ga::FunctionWrapper{Vector{Float64}, Tuple{Vector{Float64},Vector{Float64}}}
end

# reduced version of what is in definitions.jl
function DummyModelFW(p::Array{Symbol}, G::Expr)
    unpack_expr = Meta.parse(string(join(string.(p), ","), ",=params"))
    @eval return(DummyModelFW($p,
        ((ta, params) -> begin $unpack_expr; [$G for t in ta]; end) |> FunctionWrapper{Vector{Float64}, Tuple{Vector{Float64},Vector{Float64}}})
    )
end

# simper modulus, SLS_zener.G
exmoresimple = quote kᵧ + kᵦ*exp(-t*kᵦ/η) end
# complicated modulus, FractSLS_Zener.G
exmorecomplex = quote kᵦ*mittleff(a, -(kᵦ/cₐ)*t^a) + kᵧ end

fwstruct_simple = DummyModelFW([:kᵦ, :kᵧ, :η], exmoresimple)
fwstruct_complex = DummyModelFW([:cₐ, :a, :kᵦ, :kᵧ], exmorecomplex)

##############
# tuple set-up

struct DummyModelTuples{T<:Function}
  free_params::Tuple
  G::T
end

# simper modulus, SLS_zener.G
function moresimple(t; kᵦ, kᵧ, η)
  kᵧ + kᵦ*exp(-t*kᵦ/η)
end
# complicated modulus, FractSLS_Zener.G
function morecomplex(t; cₐ, a, kᵦ, kᵧ)
  kᵦ*mittleff(a, -(kᵦ/cₐ)*t^a) + kᵧ
end

params_syms_simple = (:kᵦ, :kᵧ, :η)
tupstruct_simple = DummyModelTuples(params_syms_simple, moresimple)

params_syms_complex = (:cₐ, :a, :kᵦ, :kᵧ)
tupstruct_complex = DummyModelTuples(params_syms_complex, morecomplex)

######################
# cost function set-up

# reduced & modified versions of what is in base.jl

function cost_fw(params_numbers, modulus, time_series, dt, prescribed_dot, measured)

    modulus_numbers = modulus(time_series, params_numbers)

    convolved = conv(modulus_numbers, prescribed_dot)
    convolved = convolved[1:length(time_series)]*dt

    cost = sum((measured - convolved).^2)
    return cost
end

function cost_tuple(params_numbers, params_symbols, modulus, time_series, dt, prescribed_dot, measured)

    kwargz = NamedTuple{params_symbols}(params_numbers)
    modulus_numbers = modulus.(time_series; kwargz...)

    convolved = conv(modulus_numbers, prescribed_dot)
    convolved = convolved[1:length(time_series)]*dt

    cost = sum((measured - convolved).^2)
    return cost
end

###########
# benchmark

dt = 0.1
ta = Float64.(0.0:dt:100) # representative number of samples, 1000
prescribed = 2*ta .+ 1 # representative loading derivative
measured = 5*ta .- 1 # again, exact form doesn't matter

params_numbers_simple = [5.0, 10.0, 10.0]
params_numbers_complex = [5.0, 0.5, 10.0, 10.0]

@btime cost_fw($params_numbers_simple, $fwstruct_simple.Ga, $ta, $dt, $prescribed, $measured)
@btime cost_tuple($params_numbers_simple, $params_syms_simple, $tupstruct_simple.G, $ta, $dt, $prescribed, $measured)
# 113.005 μs (106 allocations: 126.23 KiB)
# 120.539 μs (128 allocations: 126.92 KiB)
# tuples add ~6% overhead for simple modulus

@btime cost_fw($params_numbers_complex, $fwstruct_complex.Ga, $ta, $dt, $prescribed, $measured)
@btime cost_tuple($params_numbers_complex, $params_syms_complex, $tupstruct_complex.G, $ta, $dt, $prescribed, $measured)
# 203.411 μs (106 allocations: 126.23 KiB)
# 210.341 μs (129 allocations: 127.05 KiB)
# named tuples add ~3% overhead complex modulus

max_stars_repo_name: JuliaRheology/RheoBenchAndTest

max_issues_repo_name: JuliaRheology/RheoBenchAndTest

max_forks_repo_name: JuliaRheology/RheoBenchAndTest


text:
# ──────────────────────────────────────────────────────────────
#  Various KKT Solvers
# ──────────────────────────────────────────────────────────────

"""
Solves the 3x3 system
```
┌             ┐ ┌    ┐   ┌   ┐
│ Q   G'  -A' │ │ y' │ = │ y │
│ G           │ │ w' │   │ w │
│ A       FᵀF │ │ v' │   │ v │
└             ┘ └    ┘   └   ┘
```
by the double QR method described in CVXOPT
http://www.seas.ucla.edu/~vandenbe/publications/coneprog.pdf
section 10.2
"""
function kktsolver_qr(Q, A, G, cone_dims)

  n = size(Q,1) # Number of variables
  m = size(A,1) # Number of inequality constraints
  p = size(G,1) # Number of equality constraints

  (Q0,R1) = qr(full(G'), thin = false)
  Q1      = Q0[:,1:p]
  Q2      = Q0[:,p+1:end]

  function solve3x3gen(F, F⁻ᵀ)

    F⁻ᵀ     = full(inv(F))'
    Atil    = F⁻ᵀ*full(A)
    QpAᵀA   = Q + Atil'Atil
    L = qrfact(Q2'*(QpAᵀA)*Q2)

    function solve3x3(bx, by, bz)

      Q1ᵀx = R1'\by
      Q2ᵀx = L \ ( Q2'*(bx + Atil'*(F⁻ᵀ*bz)) -
                   Q2'*((QpAᵀA)*(Q1*(Q1ᵀx))) )
      y    = R1 \ ( Q1'*(bx + Atil'*(F⁻ᵀ*bz))    -
                    Q1'*(QpAᵀA*(Q1*Q1ᵀx)) -
                    Q1'*(QpAᵀA*(Q2*Q2ᵀx)) )
      x    = Q0'\[Q1ᵀx; Q2ᵀx]
      Fz   = ( F⁻ᵀ*bz -
               Atil*(Q1*(Q1ᵀx)) - Atil*(Q2*(Q2ᵀx)) )
      z    = inv(F)*Fz

      return (x,y,z)

    end

    return solve3x3

  end

end

function lift(F::Block)

  d = zeros(0)

  IA, JA, VA = Int[], Int[], Float64[]
  IB, JB, VB = Int[], Int[], Float64[]
  ID, JD, VD = Int[], Int[], Float64[]

  n = block_idx(F)[end][end]
  Ir = 0   # Index of top right coordinate for expansion

  for (In,Blk) = zip(block_idx(F), F.Blocks)

    if isa(Blk, SymWoodbury)

      for i = 1:length(Blk.A.diag)
        push!(IA,In[i]); push!(JA,In[i]); push!(VA,Blk.A.diag[i])
      end

      for i = 1:size(Blk.B,1), j = 1:size(Blk.B,2)
        push!(IB,In[i]); push!(JB,Ir+j); push!(VB,Blk.B[i,j])
      end

      for i = 1:size(Blk.D,1), j = 1:size(Blk.D,2)
        invD = inv(Blk.D)
        if Blk.D[i,j] != 0
          push!(ID,Ir+i); push!(JD,Ir+j); push!(VD,-invD[i,j])
        end
      end

      Ir = Ir + size(Blk.B,2)

    end

    if isa(Blk, Diag)

      for i = 1:length(Blk.diag)
        push!(IA, In[i]); push!(JA, In[i]); push!(VA, Blk.diag[i])
      end

    end

  end
  return (sparse(IA,JA,VA), sparse(IB,JB,VB,n,Ir), sparse(ID,JD,VD));

end

"""
Estimates for the number of nonzeros of lift(F)
"""
function count_lift(cone_dims)
  n = 0
  for (btype, k) = cone_dims
    if btype == "Q"; n = n + k + 2*(2*k) + 4;  end
    if btype == "R"; n = n + k; end
    if btype == "S"; n = n + k^2; end
  end
  return n
end

"""
Estimates the number of nonzeros of F
"""
function count_dense(cone_dims)
  n = 0
  for (btype, k) = cone_dims
    if btype == "Q"; n = n + k^2;  end
    if btype == "R"; n = n + k; end
    if btype == "S"; n = n + k^2; end
  end
  return n
end

"""
Creates a matrix with the same sparsity structure as F
"""
function placeholder(cone_dims)
  num_cones = length(cone_dims)
  B = Block(num_cones);
  for i = 1:num_cones
    (ctype, k) = cone_dims[i]
    if ctype == "R"; B[i] = ConicIP.Diag(2*rand(k)); end
    if ctype == "Q"; B[i] = SymWoodbury(ConicIP.Diag(3*rand(k)), rand(k), 1.); end
    if ctype == "S"; B[i] = ConicIP.VecCongurance(ConicIP.mat(rand(k)) + I); end
  end
  return B
end

"""
Checks if two sparse matrices have the same sparse structure
"""
function identical_sparse_structure(A::SparseMatrixCSC,B::SparseMatrixCSC)
  if length(A.nzval) != length(B.nzval)
    return false
  end
  if ( all(i -> (A.rowval[i] == B.rowval[i]), 1:length(A.rowval)) &&
       all(i -> (A.colptr[i] == B.colptr[i]), 1:length(A.colptr)) )
    return true
  end
  return false
end

"""
Solves the 3x3 system
```
┌             ┐ ┌    ┐   ┌   ┐
│ Q   G'  -A' │ │ y' │ = │ y │
│ G           │ │ w' │   │ w │
│ A       FᵀF │ │ v' │   │ v │
└             ┘ └    ┘   └   ┘
```

By lifting the large diagonal plus rank 3 blocks of FᵀF

Intelligently chooses between solve3x3gen_sparse_lift and
solve3x3gen_sparse_dense by approximating the number of non-zeros in
both and choosing the form with more sparsity. The former is better
for large second order cones, while the latter is better if the
constraints are the product of many small cones.
"""
function kktsolver_sparse(Q, A, G, cone_dims)

  n = size(Q,1) # Number of variables
  m = size(A,1) # Number of inequality constraints
  p = size(G,1) # Number of equality constraints

  Q = sparse(Q) # TODO: remove the need for this ?
  A = sparse(A)
  G = sparse(G)

  Fp = placeholder(cone_dims)

  if count_lift(cone_dims) < count_dense(cone_dims)

    # precompute sparse structure and analyze it
    (FᵀFA, FᵀFB, invFᵀFD) = lift(Fp'Fp); r = size(invFᵀFD,1)
    Z = [ Q             G'            -A'            spzeros(n,r)
          G             spzeros(p,p)   spzeros(p,m)  spzeros(p,r)
          A             spzeros(m,p)   FᵀFA          FᵀFB
          spzeros(r,n)  spzeros(r,p)   FᵀFB'         invFᵀFD        ]
    Zᶠ  = lufact(Z)
    Z.nzval[:] = 1:length(Z.nzval)
    I₁₁ = round( Int, Z[n+p+1:n+p+m,n+p+1:n+p+m].nzval )
    I₁₂ = round( Int, Z[n+p+1:n+p+m,n+p+m+1:end].nzval )
    I₂₁ = round( Int, Z[n+p+m+1:end,n+p+1:n+p+m].nzval )
    I₂₂ = round( Int, Z[n+p+m+1:end,n+p+m+1:end].nzval )

    function solve3x3gen_lift(F, F⁻ᵀ)

      (FᵀFA, FᵀFB, invFᵀFD) = lift(F'F); r = size(invFᵀFD,1)
      # In the first iteration, FᵀF is the identity. This
      # detects that
      if r == 0
        Z₀ = [ Q        G'             -A'
               G        spzeros(p,p)   spzeros(p,m)
               A        spzeros(m,p)   FᵀFA         ]
        Z₀ᶠ = lufact(Z₀)
        function solve3x3I(Δy, Δw, Δv)
          z = Z₀ᶠ\[Δy; Δw; Δv]
          return (z[1:n,:], z[n+1:n+p,:], z[n+p+1:end,:])
        end
        return solve3x3I
      else
        # If the sparsity structure is the same, you can reuse
        # the symbolic factorization.
        Zᶠ.nzval[I₁₁] = FᵀFA.nzval
        Zᶠ.nzval[I₁₂] = FᵀFB.nzval
        Zᶠ.nzval[I₂₁] = FᵀFB'.nzval
        Zᶠ.nzval[I₂₂] = invFᵀFD.nzval
        Zᶠ.numeric = C_NULL; SparseArrays.UMFPACK.umfpack_numeric!(Zᶠ)
        function solve3x3lift(Δy, Δw, Δv)
          z = Zᶠ\[Δy; Δw; Δv; zeros(r,1)]
          return (z[1:n,:], z[n+1:n+p,:], z[(n+p+1):(n+m+p),:])
        end
        return solve3x3lift
      end
    end

    return solve3x3gen_lift

  else

    # Compute sparse structure and analyze it
    FᵀFp = sparse(Fp'Fp)
    Z = [ Q        G'             -A'
          G        spzeros(p,p)   spzeros(p,m)
          A        spzeros(m,p)   FᵀFp         ]
    Zᶠ = lufact(Z)
    Z.nzval[:] = 1:length(Z.nzval)
    Iᵤ = round(Int,Z[end-m+1:end,end-m+1:end].nzval)

    function solve3x3gen_nolift(F, F⁻ᵀ)

      FᵀF = sparse(F'F)

      # Check if FᵀF has the same sparsity structure as FᵀF₀
      # If not, construct a new matrix
      if !identical_sparse_structure(FᵀF, FᵀFp)
        Z₀ = [ Q        G'             -A'
               G        spzeros(p,p)   spzeros(p,m)
               A        spzeros(m,p)   FᵀF          ]
        Z₀ᶠ = lufact(Z₀)
        function solve3x3I(Δy, Δw, Δv)
          z = Z₀ᶠ\[Δy; Δw; Δv]
          return (z[1:n,:], z[n+1:n+p,:], z[n+p+1:end,:])
        end
        return solve3x3I
      else
        # If the sparsity structure is the same, you can reuse
        # the symbolic factorization.
        Zᶠ.nzval[Iᵤ] = FᵀF.nzval
        Zᶠ.numeric = C_NULL
        SparseArrays.UMFPACK.umfpack_numeric!(Zᶠ)
        function solve3x3_nolift(Δy, Δw, Δv)
          z = Zᶠ\[Δy; Δw; Δv]
          return (z[1:n,:], z[n+1:n+p,:], z[n+p+1:end,:])
        end
        return solve3x3_nolift
      end

    end

    return solve3x3gen_nolift

  end

end

"""
Solves the 2x2 system
```
┌                   ┐ ┌    ┐   ┌   ┐
│ Q + A'F⁻¹F⁻ᵀA  G' │ │ y' │ = │ y │
│ G                 │ │ w' │   │ w │
└                   ┘ └    ┘   └   ┘
```
"""
function kktsolver_2x2(Q, A, G, cone_dims)

  n = size(Q,1) # Number of variables
  m = size(A,1) # Number of inequality constraints
  p = size(G,1) # Number of equality constraints

  function solve2x2gen(F, F⁻ᵀ)

    F⁻ᵀ = sparse(F⁻ᵀ)
    AᵀF⁻¹F⁻ᵀA = A'*(F⁻ᵀ'*(F⁻ᵀ*A))

    Z = [ Q + AᵀF⁻¹F⁻ᵀA   G'
          G               spzeros(p,p) ]

    Z = lufact(Z)

    function solve2x2(Δy, Δw)

      z = Z\[Δy; Δw]
      return (z[1:n,:], z[n+1:end,:])

    end

    return solve2x2

  end

  return solve2x2gen

end

"""
Wrapper around solve2xegen to solve 3x3 systems by pivoting
on the third component.
"""
function pivotgen(kktsolver_2x2,Q,A,G,cone_dims)

  solve2x2gen = kktsolver_2x2(Q,A,G,cone_dims)

  function solve3x3gen(F, F⁻ᵀ)

    solve2x2 = solve2x2gen(F, F⁻ᵀ)

    function solve3x3(y, w, v)

      t1 = F⁻ᵀ*(F⁻ᵀ*v)
      (Δy, Δw) = solve2x2(y + A'*t1, w)
      axpy!(-1, F⁻ᵀ*(F⁻ᵀ*(A*Δy)), t1)  # Δv = F⁻²*(v - A*Δy)

      return(Δy, Δw, t1)

    end

  end

  return solve3x3gen

end

pivot(kktsolver_2x2) = (Q,A,G,cone_dims) -> pivotgen(kktsolver_2x2,Q,A,G,cone_dims)


text:
""" Scheduling and evaluation of undirected wiring diagrams.

In category-theoretic terms, this module is about evaluating arbitrary
composites of morphisms in hypergraph categories.
"""
module ScheduleUndirectedWiringDiagrams
export AbstractNestedUWD, AbstractScheduledUWD, NestedUWD, ScheduledUWD,
  SchedulingAlgorithm, SequentialSchedule,
  eval_schedule, to_nested_diagram, schedule

import Base: schedule
using DataStructures: IntDisjointSets, union!, in_same_set

using ...Present, ...CategoricalAlgebra.CSets, ...CategoricalAlgebra.FinSets
using ..UndirectedWiringDiagrams, ..WiringDiagramAlgebras
using ..UndirectedWiringDiagrams: TheoryUWD, flat

# Data types
############

""" Abstract type for C-sets that contain a scheduled UWDS.

This type includes [`ScheduledUWD`](@ref) and [`NestedUWD`](@ref).
"""
@abstract_acset_type HasScheduledUWD <: HasUWD

@present TheoryScheduledUWD <: TheoryUWD begin
  Composite::Ob

  parent::Hom(Composite, Composite)
  box_parent::Hom(Box, Composite)

  # Po-category axiom enforcing that composites form a rooted forest.
  # parent <= id(Composite)
end

@abstract_acset_type AbstractScheduledUWD <: HasScheduledUWD

""" Scheduled undirected wiring diagram.

A scheduled UWD consists of a UWD together with a set of *composite* nodes
forming a rooted tree, or in general a rooted forest, whose leaves are the
diagram's boxes.

See also: [`NestedUWD`](@ref).
"""
@acset_type ScheduledUWD(TheoryScheduledUWD,
  index=[:box, :junction, :outer_junction, :parent, :box_parent]) <: AbstractScheduledUWD

ncomposites(x::HasScheduledUWD) = nparts(x, :Composite)
composites(x::HasScheduledUWD) = parts(x, :Composite)
parent(x::HasScheduledUWD, args...) = subpart(x, args..., :parent)
children(x::HasScheduledUWD, c::Int) =
  filter(c′ -> c′ != c, incident(x, c, :parent))
box_parent(x::HasScheduledUWD, args...) = subpart(x, args..., :box_parent)
box_children(x::HasScheduledUWD, args...) = incident(x, args..., :box_parent)

""" Abstract type for C-sets that contained a nested UWD.
"""
@abstract_acset_type HasNestedUWD <: HasScheduledUWD

@present TheoryNestedUWD <: TheoryScheduledUWD begin
  CompositePort::Ob

  composite::Hom(CompositePort, Composite)
  composite_junction::Hom(CompositePort, Junction)
end

@abstract_acset_type AbstractNestedUWD <: HasNestedUWD

""" Nested undirected wiring diagram.

A nested UWD is a scheduled UWD whose composite nodes have been given ports,
making explicit the intermediate boxes in the composition.

Nested UWDs are very similar but not quite identical to Robin Milner's
[bigraphs](https://doi.org/10.1017/CBO9780511626661).

See also: [`ScheduledUWD`](@ref).
"""
@acset_type NestedUWD(TheoryNestedUWD,
  index=[:box, :junction, :outer_junction,
         :composite, :composite_junction, :parent, :box_parent]) <: AbstractNestedUWD

composite_ports(x::HasNestedUWD, args...) = incident(x, args..., :composite)
composite_junction(x::HasNestedUWD, args...) =
  subpart(x, args..., :composite_junction)
composite_ports_with_junction(x::HasNestedUWD, args...) =
  incident(x, args..., :composite_junction)

# Evaluation
############

""" Evaluate a scheduled UWD given a set of generators for the boxes.

The optional first argument `f` should be a callable with the same signature as
[`oapply`](@ref) for UWD algebras, which by default is `oapply` itself.

```julia
f(composite::UndirectedWiringDiagram, values::AbstractVector)
```

Nested UWDs are used as an auxiliary data structure during the evaluation. If
the schedule will be evaluated multiple times, you should explicitly convert it
to a nested UWD using [`to_nested_diagram`](@ref) and then call `eval_schedule`
on the resulting object instead.
"""
eval_schedule(schedule, generators::AbstractVector) =
  eval_schedule(oapply, schedule, generators)

eval_schedule(f, schedule::AbstractScheduledUWD, generators::AbstractVector) =
  eval_schedule(f, to_nested_diagram(schedule), generators)

function eval_schedule(f, d::AbstractNestedUWD, generators::AbstractVector)
  # Evaluate `f` after constructing UWD for the composite.
  function do_eval(values, juncs, outer_junc)
    # Create diagram, adding boxes and ports.
    composite = UndirectedWiringDiagram(length(outer_junc))
    for junc in juncs
      add_box!(composite, length(junc))
    end

    # Set junctions, normalizing junction IDs to consecutive numbers `1:n`.
    jmap = Dict{Int,Int}()
    mapj!(j) = get!(jmap, j) do; add_junction!(composite) end
    for (port, j) in enumerate(Iterators.flatten(juncs))
      set_subpart!(composite, port, :junction, mapj!(j))
    end
    for (port, j) in enumerate(outer_junc)
      set_subpart!(composite, port, :outer_junction, mapj!(j))
    end

    f(composite, values)
  end

  # Mutually recursively evaluate children of composite `c`.
  function eval_children(c::Int)
    bs, cs = box_children(d, c), children(d, c)
    values = isempty(cs) ? generators[bs] : # XXX: Avoid Any[].
      [ generators[bs]; map(eval_composite, cs) ]
    juncs = [ [junction(d, ports(d, b)) for b in bs];
              [composite_junction(d, composite_ports(d, c′)) for c′ in cs] ]
    (values, juncs)
  end

  # Mutually recursively evaluate composite `c`.
  function eval_composite(c::Int)
    values, juncs = eval_children(c)
    outer_junc = composite_junction(d, composite_ports(d, c))
    do_eval(values, juncs, outer_junc)
  end

  # Evaluate diagram starting at the root, assumed to be unique. The root
  # composite must be handled specially due to outer ports.
  root = only(filter(c -> parent(d, c) == c, composites(d)))
  values, juncs = eval_children(root)
  do_eval(values, juncs, junction(d, outer=true))
end

""" Convert a scheduled UWD to a nested UWD.
"""
function to_nested_diagram(s::AbstractScheduledUWD)
  d = NestedUWD()
  copy_parts!(d, s)

  n = nboxes(s)
  sets = IntDisjointSets(n + ncomposites(s))

  function add_composite_ports!(c::Int)
    # Recursively add ports to all child composites.
    foreach(add_composite_ports!, children(d, c))

    # Get all junctions incident to any child box or child composite.
    js = unique!(sort!(flat([
      [ junction(d, ports(d, b)) for b in box_children(d, c) ];
      [ composite_junction(d, composite_ports(d, c′)) for c′ in children(d, c) ]
    ])))

    # Filter for "outgoing" junctions, namely those incident to any outer port
    # or to any port of a box that is not a descendant of this composite.
    c_rep = n+c
    for b in box_children(d, c); union!(sets, c_rep, b) end
    for c′ in children(d, c); union!(sets, c_rep, n+c′) end
    js = filter!(js) do j
      !(all(in_same_set(sets, c_rep, box(d, port))
            for port in ports_with_junction(d, j)) &&
        isempty(ports_with_junction(d, j, outer=true)))
    end

    # Add port for each outgoing junction.
    add_parts!(d, :CompositePort, length(js), composite=c, composite_junction=js)
  end

  # Add ports to each composite, starting at roots.
  roots = filter(c -> parent(d, c) == c, composites(d))
  foreach(add_composite_ports!, roots)
  return d
end

# Scheduling
############

abstract type SchedulingAlgorithm end

""" Schedule diagram as a sequential chain of binary composites.

This is the simplest possible scheduling algorithm, the equivalent of `foldl`
for undirected wiring diagrams. Unless otherwise specified, the boxes are folded
according to the order of their IDs, which is arbitrary.
"""
struct SequentialSchedule <: SchedulingAlgorithm end

""" Schedule an undirected wiring diagram.

By default, a simple sequential schedule is used.
"""
function schedule(d::AbstractUWD;
                  alg::SchedulingAlgorithm=SequentialSchedule(), kw...)
  schedule(d, alg; kw...)
end

function schedule(d::AbstractUWD, ::SequentialSchedule;
                  order::Union{AbstractVector{Int},Nothing}=nothing)
  if isnothing(order)
    order = boxes(d)
  end
  nb = nboxes(d)
  @assert length(order) == nb
  nc = max(1, nb-1)

  schedule = ScheduledUWD()
  copy_parts!(schedule, d)
  add_parts!(schedule, :Composite, nc, parent=[2:nc; nc])
  set_subpart!(schedule, order[1:min(2,nb)], :box_parent, 1)
  set_subpart!(schedule, order[3:nb], :box_parent, 2:nc)
  schedule
end

end


text:
#! /usr/bin/env julia

using CSV
using DataFrames
using MAT
using Statistics

@info "Read $(ARGS[1])"

file = matopen(ARGS[1]); # "goodPoints.mat"

good_points = read(file, "goodPoints");
good_bases = vec(read(file, "goodBases"));

close(file)

@info "Read $(ARGS[2])"

gene_per_base = CSV.read(ARGS[2], header=0); # "cell_barcode_names.csv"
gene_per_base = Dict("$(row[2])" => row[3] for row in eachrow(gene_per_base));

df_spatial = DataFrame(good_points, [:x, :y, :z]);
df_spatial[!, :gene] = [gene_per_base[b] for b in good_bases];

CSV.write("molecules.csv", df_spatial)

if length(ARGS) > 2
    import NPZ
    import Images

    @info "Read $(ARGS[3])"
    Images.save("segmentation.tiff", UInt16.(NPZ.npzread(ARGS[3])["labels"]))
end

max_stars_repo_name: YinuoJin/Baysor

max_issues_repo_name: YinuoJin/Baysor

max_forks_repo_name: YinuoJin/Baysor


text:
module LakeBiquads

export EQBiquads, BiquadArray, AllPass, HighPass, LowPass, Peak, HighShelf, LowShelf, Bessel

using DelimitedFiles, DSP, NLsolve

include("Peak.jl")
include("Shelf.jl")
include("LowHighPass.jl")
include("Allpass.jl")

"""
	BiquadArray(x::SecondOrderSections)

Construct 2D Array of Biquad Parameters from `x`.

The returned Array for 2 Biquads is of the structure:
```
2×5 Array{Float64,2}:
 b0  b1  b2  a1  a2
 b0  b1  b2  a1  a2
```

Note that the Gain of the SecondOrderSections is spread evenly over all Biquads.
"""
function BiquadArray(x::SecondOrderSections)
	biquads = x.biquads
	n = length(biquads)
	gain = x.g^(1/n)
	out = zeros(n, 5)
	for i in 1:n
		out[i, :] = [gain*biquads[i].b0 gain*biquads[i].b1 gain*biquads[i].b2 biquads[i].a1 biquads[i].a2]
	end
	return out
end

"""
	EQBiquads(filename)

Parse the given EQ file as Lake EQs and create a new file "nameBiquads.format"
with 96kHz Biquads in space seperated format.

The EQ file allows specification of the Allpass, Peak, Highshelf, Lowshelf,
Lowpass and Highpass Filters in the following format:

```
AP f order BWoct
P f dBGain BWoct
HS f dBGain BWoct
LS f dBGain BWoct
[LP/HP][BT/BS/LR][6/12/18/24/30/36/42/48] f
```

# Example

	EQBiquads("myEQ.txt")

would create a new file named "myEQBiquads.txt" containing Biquads.

# EQ File Examples

```
LPBT30 1000  #Butterworth Lowpass with 30 dB/oct fall-off at 1000 Hz
P 1k 10 1    #Peak at 1 kHz with 10 dB Gain and a Bandwidth of 1 octave
AP 3e3 1 2   #Allpass at 3 kHz of order 1 with a Bandwidth of 2 octaves (order may only be 1 or 2)
```
"""
function EQBiquads(filename)
	EQs = readdlm(filename, ' ')
	n = size(EQs, 1)
	BiquadParams = Core.Array{Float64,2}(undef, 0,5)
	for i in 1:n
		#try
			EQtype = EQs[i,1]
			f = EQs[i,2]
			if f isa SubString
				f = replace(f, "k"=>".")
				f = f * "e3"
				f = parse(Float64, f)
			end
			if size(EQs,2)>=3
				dBGain = EQs[i,3]
				BWoct = EQs[i,4]
			end
			if EQtype == "P"
				BiquadParams = [BiquadParams; Peak(f, dBGain, BWoct)|>BiquadArray]
			elseif EQtype[1:2] in ["HP", "LP"]
				type = Symbol(EQtype[3:4])
				dBperOct = parse(Int64, EQtype[5:end])
				if EQtype[1:2] == "HP"
					bf = HighPass(type, f, dBperOct)
				else
					bf = LowPass(type, f, dBperOct)
				end
				BiquadParams = [BiquadParams; BiquadArray(bf)]
			elseif EQtype == "HS"
				BiquadParams = [BiquadParams; HighShelf(f, dBGain, BWoct)|>BiquadArray]
			elseif EQtype == "LS"
				BiquadParams = [BiquadParams; LowShelf(f, dBGain, BWoct)|>BiquadArray]
			elseif EQtype == "AP"
				BiquadParams = [BiquadParams; AllPass(f, dBGain, BWoct)|>BiquadArray] #dBGain is used here as specifying the order of the AllPass, because it controls the amount of added phase as Gain in a Peak EQ
			else
				error(string("Unknown Filter type. Execution stopped and no changes were written to filesystem. "))
			end
		# catch err
		# 	println(string("Error while parsing line ", i, ":"))
		# 	println(stacktrace())
		# 	throw(err)
		# end
	end
	#filename processing
	sep = split(filename, ".")
	newFilenamePrevious = join(sep[1:end-1], '.')
	newFilename = string(newFilenamePrevious, "Biquads.", sep[end])

	writedlm(newFilename, BiquadParams, ' ')
	#post processing to add DOS style line breaks
	str = read(newFilename, String)
	str = replace(str, "\n"=>"\r\n")
	write(newFilename, str)
end

end # module


text:
### A Pluto.jl notebook ###
# v0.17.4

using Markdown
using InteractiveUtils

# ╔═╡ 36f67d36-7d13-11ec-31df-d3b214c29ec9
using Documenter

# ╔═╡ 36f67d22-7d13-11ec-3a19-d9fc26ca1480
html"""

Happy holiday! Remember to take care of yourself and your loved ones!

"""

# ╔═╡ 868b3d5a-ba27-4b1e-b9ca-402ceef79748
using ../src/colortransfer

# ╔═╡ fbdbcae0-ab65-4a95-abbe-5a0ff9a23cf4
makedocs(sitename="STMO ZOO",
    format = Documenter.HTML(),
    modules=[Example], # add your module
    pages=Any[
        "Example"=> "man/example.md",  # add the page to your documentation
    ])

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"

[compat]
Documenter = "~0.27.12"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.0"
manifest_format = "2.0"

[[deps.ANSIColoredPrinters]]
git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c"
uuid = "a4c015fc-c6ff-483c-b24f-f7ea428134e9"
version = "0.0.1"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[deps.Documenter]]
deps = ["ANSIColoredPrinters", "Base64", "Dates", "DocStringExtensions", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"]
git-tree-sha1 = "75c6cf9d99e0efc79b724f5566726ad3ad010a01"
uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
version = "0.27.12"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "92f91ba9e5941fc781fecf5494ac1da87bdac775"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.0"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
"""

# ╔═╡ Cell order:
# ╟─36f67d22-7d13-11ec-3a19-d9fc26ca1480
# ╠═36f67d36-7d13-11ec-31df-d3b214c29ec9
# ╠═868b3d5a-ba27-4b1e-b9ca-402ceef79748
# ╠═fbdbcae0-ab65-4a95-abbe-5a0ff9a23cf4
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
#= Copyright (c) 2016 Facebook

   This program is free software; you can redistribute it and/or
   modify it under the terms of version 2 of the GNU General Public
   License as published by the Free Software Foundation.

   Ported to Julia by Julian P Samaroo, 2021
=#

ccall(:jl_exit_on_sigint, Cvoid, (Cint,), 0)

if length(ARGS) > 0
    const out_path = ARGS[1]
    const root_pids = parse.(Cint, ARGS[2:end])
else
    const out_path = nothing
    const root_pids = [getpid()]
end
_all_pids = ()
for root_pid in root_pids
    global _all_pids = (_all_pids..., parse.(Int, readdir("/proc/$root_pid/task"))...)
end
const all_pids = _all_pids::NTuple

using BPFnative
import BPFnative: API, RT, Host, VMLinux
using CBinding

const MINBLOCK_US = 1
const TASK_COMM_LEN = 32

struct Key
    waker::NTuple{TASK_COMM_LEN,UInt8}
    target::NTuple{TASK_COMM_LEN,UInt8}
    wret::UInt32
    tret::UInt32
    wpid::UInt32
    tpid::UInt32
end

const counts = RT.RTMap(;name = "counts",
                   maptype = API.BPF_MAP_TYPE_HASH,
                   keytype = Key,
                   valuetype = UInt64,
                   maxentries = 10000)

const start = RT.RTMap(;name = "start",
                  maptype = API.BPF_MAP_TYPE_HASH,
                  keytype = UInt32,
                  valuetype = UInt64,
                  maxentries = 10000)

struct WokeBy
    name::NTuple{TASK_COMM_LEN,UInt8}
    ret::UInt32
    pid::UInt32
end

const wokeby = RT.RTMap(;name = "wokeby",
                   maptype = API.BPF_MAP_TYPE_HASH,
                   keytype = UInt32,
                   valuetype = WokeBy,
                   maxentries = 10000)

const stackmap = RT.RTMap(;name = "stackmap",
                     maptype = API.BPF_MAP_TYPE_STACK_TRACE,
                     keytype = UInt32,
                     valuetype = NTuple{API.PERF_MAX_STACK_DEPTH,UInt64},
                     maxentries = 10000)

# Syscall-specific data
#const

function real_pid(ctx)
    level = unsafe_load(API.@elemptr(task.nsproxy.pid_ns_for_children.level))
    return unsafe_load(API.@elemptr(task.thread_pid.numbers[level].nr))
end

# kprobe/try_to_wake_up
function waker(ctx)
    task_ptr = API.get_param(ctx, Val(1))
    task_ptr = reinterpret(API.pointertype(API.task_struct), task_ptr)
    woke = RT.ZeroInitRef(WokeBy)

    if (pid = RT.safe_load(task_ptr.pid)) === nothing
        return 0
    end
    #=
    if !(pid in all_pids)
        return 0
    end
    =#

    RT.get_current_comm(RT.SizedBuffer(API.@elemptr(woke.name), TASK_COMM_LEN))
    sid = RT.get_stackid(ctx, stackmap, RT.BPF_F_FAST_STACK_CMP)
    if sid > 0
        unsafe_store!(API.@elemptr(woke.ret), unsafe_trunc(UInt32, sid))
        unsafe_store!(API.@elemptr(woke.pid), unsafe_trunc(UInt32, pid))
    end

    wokeby[pid] = woke[]
    return 0
end

@inline function update_counts(ctx, pid::UInt32, delta::UInt64)
    woke = RT.ZeroInitRef(WokeBy)
    key = RT.ZeroInitRef(Key)

    RT.get_current_comm(RT.SizedBuffer(API.@elemptr(key.target), TASK_COMM_LEN))
    tsid = RT.get_stackid(ctx, stackmap, RT.BPF_F_FAST_STACK_CMP)
    unsafe_store!(API.@elemptr(key.tret), reinterpret(UInt32, tsid))
    unsafe_store!(API.@elemptr(key.wret), 0)
    unsafe_store!(API.@elemptr(key.tpid), pid)
    unsafe_store!(API.@elemptr(key.wpid), 0)

    woke = wokeby[pid]
    if woke !== nothing
        unsafe_store!(API.@elemptr(key.wret), woke.ret)
        #RT.memcpy!(API.@elemptr(key.waker), woke.name, TASK_COMM_LEN)
        unsafe_store!(API.@elemptr(key.waker), woke.name)
        unsafe_store!(API.@elemptr(key.wpid), woke.pid)
        delete!(wokeby, pid)
    end

    GC.@preserve key begin
    _key = Base.unsafe_convert(Ptr{Key}, key)
    val = counts[_key]
    if val == C_NULL
        counts[_key] = 0
        #RT.map_update_elem(counts, _key, _val, UInt64(1)) # TODO: BPF_NOEXIST
        val = counts[_key]
        if val == C_NULL
            return 0
        end
    end
    counts[_key] = unsafe_load(val) + delta
    end
    return 0
end

# taken from /sys/kernel/debug/tracing/events/sched/sched_switch/format
struct sched_switch_args
    pad::UInt64
    prev_comm::NTuple{16,UInt8}
    prev_pid::Cint
    prev_prio::Cint
    prev_state::UInt64
    next_comm::NTuple{16,UInt8}
    next_pid::Cint
    next_prio::Cint
end
# tracepoint/sched/sched_switch
function oncpu(#=struct sched_switch_args *=#ctx::Ptr{sched_switch_args})
    #= record previous thread sleep time =#
    pid = unsafe_load(API.@elemptr(ctx.prev_pid))

    ts = RT.ktime_get_ns()
    start[pid] = ts
    #=
    if pid in all_pids
    end
    =#

    #= calculate current thread's delta time =#
    pid = RT.get_current_pid_tgid()[1]
    tsp = start[pid]
    if (tsp === nothing)
        #= missed start or filtered =#
        return 0
    end

    delta = RT.ktime_get_ns() - tsp
    delete!(start, pid)
    delta = delta ÷ 1000
    if (delta < MINBLOCK_US)
        return 0
    end

    return update_counts(ctx, pid, delta)
end

@info "KProbe"
kp = KProbe(waker, "try_to_wake_up"; license="GPL")
API.load(kp)

@info "Tracepoint"
tp = Tracepoint(oncpu, Tuple{Ptr{sched_switch_args}}, "sched", "sched_switch"; license="GPL", merge_with=[kp.obj])
API.load(tp)

host_counts = Host.hostmap(tp.obj, counts)
host_wokeby = Host.hostmap(tp.obj, wokeby)
host_stacks = Host.hostmap(tp.obj, stackmap)

println("Press enter to exit...")
try
    readline()
catch err
    if !(err isa InterruptException)
        rethrow(err)
    end
end

frames = Tuple{Symbol,UInt64,NTuple{(API.PERF_MAX_STACK_DEPTH*2),UInt64}}[]
proc_names = Dict{Symbol,UInt64}()
for key in keys(host_counts)
    if key.wret != 0 && key.tret != 0
        waker_str = String([key.waker...])
        target_str = String([key.target...])
        waker_sym = Symbol(replace(waker_str, '\0'=>""))
        target_sym = Symbol(replace(target_str, '\0'=>""))
        try
            #waker_addr = get!(()->rand(UInt64), proc_names, waker_sym)
            #target_addr = get!(()->rand(UInt64), proc_names, target_sym)
            waker_stack = host_stacks[key.wret]
            target_stack = host_stacks[key.tret]
            stacks = (reverse(waker_stack)..., target_stack...)
            push!(frames, (target_sym, host_counts[key], stacks))

            #= TODO: Log some extra info
            if occursin("julia", waker_str) || occursin("julia", target_str)
                println("[$waker_str ($(key.wpid)) | $target_str ($(key.tpid)): $(host_counts[key])]")
                println("====WAKER STACK====")
                print(Host.stack_to_string(waker_stack; reverse=true))
                println("===================")
                print(Host.stack_to_string(target_stack))
                println("====TARGET STACK====")
                println()
                if in(key.wpid, all_pids) || in(key.tpid, all_pids)
                    exit()
                end
            end
            =#
        catch err
            @error exception=err
        end
    end
end

API.unload(kp)
API.unload(tp)

@show length(frames)
data, lidict = Host.stack_to_frames(frames)

if out_path !== nothing
    using Serialization
    open(out_path, "w") do io
        serialize(io, (data, lidict))
    end
else
    using PProf, Profile

    pprof(data, lidict)
end


text:
# a rule within a decision tree for a single attribute value
mutable struct OneNode
    atr_value :: String     # symbol name of an attribute value (categorical data)
    prediction :: String    # predicted class for `atr_value`
    err_rate :: Real        # fraction [0 .. 1[ of how many instances are wrongly classified by this rule
    err_count :: Integer    # number of wrongly classified instances
end

# make all nodes for attribute `column`
function get_nodes(column::AbstractVector, observations::AbstractVector, uniq_obs :: AbstractVector)
    vals = unique(column)                         # possible attribute values
    nodes = Vector{OneNode}(undef, length(vals))  # for the results

    freq_table = Dict{String, Dict{String, Int}}()    # Dict of Dicts for counting frequencies

    # for each attribute value we need a Dict to count frequencies of each possible observation
    for v in vals
        atr_freq = Dict(cat => 0 for cat in uniq_obs)
        freq_table[v] = atr_freq
    end

    # let's count the observations for each attribute value
    for i in eachindex(column)
        freq_table[column[i]][observations[i]] += 1
    end

    # build the nodes of the tree (one for each attribute value)
    for i in eachindex(vals)
        maxfreq, cat = findmax(freq_table[vals[i]]) # find the most frequent category (and the number of it's occurences)
        sumfreq = sum(values(freq_table[vals[i]]))  # number or all occurences
        errfreq = sumfreq - maxfreq                 # wrongly classified occurences
        nodes[i] = OneNode(vals[i], cat, errfreq / sumfreq, errfreq)
    end
    return(nodes)
end

# print a 'OneNode' nicely formatted
function Base.show(io::IO, ::MIME"text/plain", node::OneNode)
    print(io, "OneNode: $(node.atr_value) -> $(node.prediction) (Err $(node.err_count), $(node.err_rate))")
end

# use the same formatting in compact mode
Base.show(io::IO, node::OneNode) = show(io, "text/plain", node)

max_stars_repo_name: roland-KA/OneRule

max_issues_repo_name: roland-KA/OneRule

max_forks_repo_name: roland-KA/OneRule


text:
# Some economics here
using LinearAlgebra

f(x) = x.^3

x = [2, 4]

@show (norm(f(x)))


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct DedicatedHostUpdate <: SwaggerModel
    tags::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: tags

    function DedicatedHostUpdate(;tags=nothing)
        o = new()
        validate_property(DedicatedHostUpdate, Symbol("tags"), tags)
        setfield!(o, Symbol("tags"), tags)
        o
    end
end # type DedicatedHostUpdate

const _property_map_DedicatedHostUpdate = Dict{Symbol,Symbol}(Symbol("tags")=>Symbol("tags"))
const _property_types_DedicatedHostUpdate = Dict{Symbol,String}(Symbol("tags")=>"Dict{String, String}")
Base.propertynames(::Type{ DedicatedHostUpdate }) = collect(keys(_property_map_DedicatedHostUpdate))
Swagger.property_type(::Type{ DedicatedHostUpdate }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_DedicatedHostUpdate[name]))}
Swagger.field_name(::Type{ DedicatedHostUpdate }, property_name::Symbol) =  _property_map_DedicatedHostUpdate[property_name]

function check_required(o::DedicatedHostUpdate)
    true
end

function validate_property(::Type{ DedicatedHostUpdate }, name::Symbol, val)
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -0.27492067550218746 <= q <= 2.084201684402027)

                     add_NL_constraint(m, :(sigmoid(-0.0416569430659357 + 0.13675627334278007*sigmoid(0.1752604319892206 + -0.8557289632032381*sigmoid(-0.028885055816546323 + 0.15091441453709065*sigmoid(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])) + -0.8760277243880852*sigmoid(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))) + 0.21149394336428795*sigmoid(-0.29077418744555983 + -0.11119024401630329*sigmoid(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])) + 0.4817613487326513*sigmoid(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3])))) + 0.703085236566185*sigmoid(-0.4914467550189556 + 0.5078124686195609*sigmoid(-0.028885055816546323 + 0.15091441453709065*sigmoid(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])) + -0.8760277243880852*sigmoid(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))) + 0.7523099505151869*sigmoid(-0.29077418744555983 + -0.11119024401630329*sigmoid(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])) + 0.4817613487326513*sigmoid(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))))) + sigmoid(0.6991595157542867 + 0.5132010992296534*sigmoid(0.1752604319892206 + -0.8557289632032381*sigmoid(-0.028885055816546323 + 0.15091441453709065*sigmoid(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])) + -0.8760277243880852*sigmoid(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))) + 0.21149394336428795*sigmoid(-0.29077418744555983 + -0.11119024401630329*sigmoid(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])) + 0.4817613487326513*sigmoid(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3])))) + 0.6682334830395109*sigmoid(-0.4914467550189556 + 0.5078124686195609*sigmoid(-0.028885055816546323 + 0.15091441453709065*sigmoid(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])) + -0.8760277243880852*sigmoid(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))) + 0.7523099505151869*sigmoid(-0.29077418744555983 + -0.11119024401630329*sigmoid(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])) + 0.4817613487326513*sigmoid(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module ReferenceFrameRotations

export DCM, EulerAngleAxis, EulerAngles, Quaternion

import Base: +, -, *, /, conj, copy, getindex, inv, imag, real, show
import Base: zeros
import LinearAlgebra: norm

using LinearAlgebra
using StaticArrays

# Re-export `I` from LinearAlgebra.
export I

################################################################################
#                                    Types
################################################################################

"""
The Direction Cosine Matrix of type `T` is a `SMatrix{3,3,T}`, which is a 3x3
static matrix of type `T`.

"""
DCM{T} = SMatrix{3,3,T}

"""
    struct EulerAngles{T}

The definition of Euler Angles, which is composed of three angles `a1`, `a2`,
and `a3` together with a rotation sequence `rot_seq`. The latter is provided by
a symbol with three characters, each one indicating the rotation axis of the
corresponding angle (for example, `:ZYX`). The valid values for `rot_seq` are:

* `:XYX`, `:XYZ`, `:XZX`, `:XZY`, `:YXY`, `:YXZ`, `:YZX`, `:YZY`, `:ZXY`,
  `:ZXZ`, `:ZYX`, and `ZYZ`.

"""
struct EulerAngles{T}
    a1::T
    a2::T
    a3::T
    rot_seq::Symbol
end

"""
    struct EulerAngleAxis{T}

The definition of Euler Angle and Axis to represent a 3D rotation.

# Fields

* `a`: The Euler angle [rad].
* `v`: The unitary vector aligned with the Euler axis.

"""
struct EulerAngleAxis{T}
    a::T
    v::Vector{T}
end

"""
### struct Quaternion{T}

The definition of the quaternion. It has four values of type `T`. The
quaternion representation is:

    q0 + q1.i + q2.j + q3.k

"""
struct Quaternion{T}
    q0::T
    q1::T
    q2::T
    q3::T
end

################################################################################
#                                 Deprecations
################################################################################

@deprecate eye(::Type{Quaternion{T}}) where T Quaternion{T}(I)
@deprecate eye(::Type{Quaternion}) where T Quaternion{Float64}(I)
@deprecate eye(q::Quaternion{T}) where T Quaternion(I,q)

################################################################################
#                                   Includes
################################################################################

include("compose_rotations.jl")
include("DCM.jl")
include("euler_angle_axis.jl")
include("euler_angles.jl")
include("inv_rotations.jl")
include("quaternion.jl")

end # module


text:
function get_position_aim(direction, distance, verticle_distance, horizontal_distance, aim)
  if direction == "forward"
    horizontal_distance += distance
    verticle_distance += distance * aim
  elseif direction == "up"
    aim -= distance
  else
    aim += distance
  end
  return verticle_distance, horizontal_distance, aim
end

verticle_distance = 0
horizontal_distance = 0
aim = 0

instructions = readlines("input.txt", keep=true)
function final_position_aim(instructions, get_position_aim, verticle_distance, horizontal_distance, aim)
  for instruction in instructions
    direction, distance = split(instruction)
    distance = parse(Int64, distance)
    verticle_distance, horizontal_distance, aim = get_position_aim(direction, distance, verticle_distance, horizontal_distance, aim)
  end
  return verticle_distance*horizontal_distance
end

print("The final position inclusive of aim is: ", final_position_aim(instructions, get_position_aim, verticle_distance, horizontal_distance, aim))


text:
module ols

export tols, summary

# Author: Adam Savitzky
# Email: asavitzky@forio.com
# Github: github.com/adambom

# Ported from the Python implemented by Vincent Nijs
# http://www.scipy.org/Cookbook/OLS?action=AttachFile&do=get&target=ols.0.2.py

# Julia type for multiple (multivariate) regression using OLS
# For least squared regression on linear equations of multiple independent variables
# y = a1 * x1 + a2 * x2 + ... an * xn
# Y = AX + E

# Input
## y = dependent variable
## y_varnm = string with the variable label for y
## x = independent variables, note that a constant is added by default
## x_varnm = list of variable labels for the independent variables

# Usage
## Instantiate a new ols type
### reg = ols(y, x, "y", ["x1", "x2", "x3"])
### Coefficients: reg.b
### R-Squared: reg.R2
### F-Statistic: reg.F
### Summary: summary(reg)

type tols
    y::Array{Float64}
    x::Array{Float64}
    y_varnm::String
    x_varnm::Array{String, 1}
    inv_xx::Array{Float64}
    b::Array{Float64, 1}
    nobs::Int
    ncoef::Int
    df_e::Int
    df_r::Int
    er::Array
    sse::Float64
    se::Array{Float64, 1}
    t::Array{Float64}
    #p::Array
    R2::Float64
    R2adj::Float64
    F::Float64
    #Fpv::Float64

    function tols(y, x, y_varnm, x_varnm)
        x = hcat(ones(size(x, 1)), x)
        xT = transpose(x)

        inv_xx = inv(xT * x)
        xy = xT * y
        b = inv_xx * xy           # estimate coefficients

        nobs = size(y, 1)         # number of observations
        ncoef = size(x, 2)        # number of coefficients
        df_e = nobs - ncoef       # degrees of freedom, error
        df_r = ncoef - 1          # degrees of freedom, regression

        er = y - x * b            # residuals
        sse = e^2/df_e            # SSE
        se = sqrt(diag(sse * inv_xx)) # coef. standard errors
        t = b / se                # coef. t-statistics
        # p = (1 - cdf(abs(t), df_e)) * 2 # coef. p-values

        R2 = 1 - var(er) / var(y) # model R-squared
        R2adj = 1 - (1 - R2) * ((nobs - 1) / (nobs - ncoef)) # adjusted R-square

        F = (R2 / df_r) / ((1 - R2) / df_e) # model F-statistic
        # Fpv = 1 - cdf(F, df_r, df_e) # F-statistic p-value

        new(y, x, y_varnm, x_varnm, inv_xx, b, nobs, ncoef, df_e, df_r, er, sse, se, t, R2, R2adj, F)
    end
end

function dw(self::tols)
    # Calculates the Durbin-Waston statistic
    de = self.er - 1.
    result = dot(de, de) / dot(self.er, self.er)
    return result
end

function ll(self::tols)
    # Calculate model log-likelihood and two information criteria

    # Model log-likelihood, AIC, and BIC criterion values
    loglike = -(self.nobs / 2) * (1 + log(2pi)) - (self.nobs / 2) * log(dot(self.er, self.er) / self.nobs)
    aic = -2loglike / self.nobs + (2 * self.ncoef / self.nobs)
    bic = -2loglike / self.nobs + (self.ncoef * log(self.nobs)) / self.nobs

    return loglike, aic, bic
end

function summary(self::tols)
    # print model output to screen

    t = time()

    # extra stats
    loglike, aic, bic = ll(self)
    #JB, JBpv, skew, kurtosis = self.JB()
    #omni, omnipv = self.omni()

    println("==============================================================================")
    #println("Dependent Variable: " + self.y_varnm)
    println("Method: Least Squares")
    println("Time: $t")
    println("No. obs:               $(self.nobs)")
    println("No. variables:         $(self.ncoef)")
    println("==============================================================================")
    println("variable             coefficient             std. Error              t-statistic")
    println("==============================================================================")
    for i in 1:length(self.x_varnm)
        println("$(self.x_varnm[i])         $(self.b[i])        $(self.se[i])        $(self.t[i])")
    end
    println("========================================================================================")
    println("Model Stats                                    Residual Stats")
    println("========================================================================================")
    println("R-Squared            $(self.R2)        Durbin-Watson Stat: $(dw(self))")
    println("Adjusted R-Squared   $(self.R2adj)        Omnimbus Stat:      ?")
    println("F-Statistic          $(self.F)        Prob(Omnibus stat): ?")
    println("Log-Likelihood       $loglike       Prob(JB):           ?")
    println("AIC Criterion        $aic        Skew:               ?")
    println("BIC Criterion        $bic         Kurtosis:           ?")
    println("========================================================================================")
end

function linreg{T<:Number}(X::StridedVecOrMat{T}, y::Vector{T})
    hcat(ones(T, size(X,1)), X)\y
end

end

max_issues_repo_path: Module 3/Chapter05/ols.jl

max_forks_repo_path: Module 3/Chapter05/ols.jl


text:
# These are the additions to the types_and_consts.jl file not from clang.jl

function Base.convert(::Type{Matrix}, J::DlsMat)
    _dlsmat = unsafe_load(J)
    # own is false as memory is allocated by sundials
    unsafe_wrap(Array, _dlsmat.data, (_dlsmat.M, _dlsmat.N), own = false)
end

CVRhsFn_wrapper(fp::CVRhsFn) = fp
CVRhsFn_wrapper(f) = @cfunction($f, Cint, (realtype, N_Vector, N_Vector, Ptr{Cvoid})).ptr

ARKRhsFn_wrapper(fp::ARKRhsFn) = fp
ARKRhsFn_wrapper(f) = @cfunction($f, Cint, (realtype, N_Vector, N_Vector, Ptr{Cvoid})).ptr

CVRootFn_wrapper(fp::CVRootFn) = fp
CVRootFn_wrapper(f) =
    @cfunction($f, Cint, (realtype, N_Vector, Ptr{realtype}, Ptr{Cvoid})).ptr

CVQuadRhsFn_wrapper(fp::CVQuadRhsFn) = fp
CVQuadRhsFn_wrapper(f) =
    @cfunction($f, Cint, (realtype, N_Vector, N_Vector, Ptr{Cvoid})).ptr

IDAResFn_wrapper(fp::IDAResFn) = fp
IDAResFn_wrapper(f) =
    @cfunction($f, Cint, (realtype, N_Vector, N_Vector, N_Vector, Ptr{Cvoid})).ptr

IDARootFn_wrapper(fp::IDARootFn) = fp
IDARootFn_wrapper(f) =
    @cfunction($f, Cint, (realtype, N_Vector, N_Vector, Ptr{realtype}, Ptr{Cvoid})).ptr

KINSysFn_wrapper(fp::KINSysFn) = fp
KINSysFn_wrapper(f) = @cfunction($f, Cint, (N_Vector, N_Vector, Ptr{Cvoid})).ptr

function Base.convert(::Type{Matrix}, J::SUNMatrix)
    _sunmat = unsafe_load(J)
    _mat = convert(SUNMatrixContent_Dense, _sunmat.content)
    mat = unsafe_load(_mat)
    # own is false as memory is allocated by sundials
    unsafe_wrap(Array, mat.data, (mat.M, mat.N), own = false)
end

function Base.convert(::Type{SparseArrays.SparseMatrixCSC}, J::SUNMatrix)
    _sunmat = unsafe_load(J)
    _mat = convert(SUNMatrixContent_Sparse, _sunmat.content)
    mat = unsafe_load(_mat)
    # own is false as memory is allocated by sundials
    # TODO: Get rid of allocation for 1-based index change
    rowval = unsafe_wrap(Array, mat.indexvals, (mat.NNZ), own = false)
    colptr = unsafe_wrap(Array, mat.indexptrs, (mat.NP + 1), own = false)
    colptr .+= 1
    m = mat.M
    n = mat.N
    nzval = unsafe_wrap(Array, mat.data, (mat.NNZ), own = false)
    SparseArrays.SparseMatrixCSC(m, n, colptr, rowval, nzval)
end

abstract type SundialsMatrix end
struct DenseMatrix <: SundialsMatrix end
struct BandMatrix <: SundialsMatrix end
struct SparseMatrix <: SundialsMatrix end

abstract type SundialsLinearSolver end
struct Dense <: SundialsLinearSolver end
struct Band <: SundialsLinearSolver end
struct SPGMR <: SundialsLinearSolver end
struct SPFGMR <: SundialsLinearSolver end
struct SPBCGS <: SundialsLinearSolver end
struct PCG <: SundialsLinearSolver end
struct PTFQMR <: SundialsLinearSolver end
struct KLU <: SundialsLinearSolver end
struct LapackBand <: SundialsLinearSolver end
struct LapackDense <: SundialsLinearSolver end

abstract type SundialsNonLinearSolver end
struct Newton <: SundialsNonLinearSolver end
struct FixedPoint <: SundialsNonLinearSolver end

abstract type StiffnessChoice end
struct Explicit <: StiffnessChoice end
struct Implicit <: StiffnessChoice end


text:
Try.map(f, ok::Ok) = Ok(f(ok.value))
Try.map(_, err::Err) = err

function Try.map(f, result)
    if Try.isok(result)
        Ok(f(Try.unwrap(result)))
    else
        # Need to "forget" about any Ok value type information:
        Err(Try.unwrap_err(result))
    end
end

Try.map(f, some::Some) = Some(f(something(some)))
Try.map(_, ::Nothing) = nothing

struct TryMap{F} <: Function
    f::F
end
TryMap(::Type{T}) where {T} = TryMap{Type{T}}(T)

(f::TryMap)(x) = Try.map(f.f, x)

Try.map(f::F) where {F} = TryMap(f)


text:
function main(r::Robot)
    path = go_to_west_south_corner_and_return_path!(r; go_around_barriers = true)

    side = East
    while !isborder(r,North)
        go_to_border_and_return_path!(r, side; markers = true, go_around_barriers = true)
        side = inverse_side(side)
        go!(r, North)
    end

    go_to_border_and_return_path!(r, side; markers = true, go_around_barriers = true)

    go_to_west_south_corner_and_return_path!(r, go_around_barriers = true)
    go_by_path!(r, path)
end

max_stars_repo_name: Savko-fokin/mirea-progs

max_issues_repo_name: Savko-fokin/mirea-progs

max_forks_repo_name: Savko-fokin/mirea-progs


text:
function tensorproduct(x::Array{Float64,2},porder::Int)::Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,2},Array{Float64,2}}
# #TENSORPRODUCT Vandermonde matrices for tensor product polynomials in [0,1]^d
# #
# #   [F,FX,FY,FZ] = TENSORPRODUCT(X,P)
# #
# #      X:         Coordinates of the points wherethe polynomials
# #                 are to be evaluated (npoints,dim)
# #      PORDER:    Maximum order of the polynomials consider. That
# #                 is all polynomials of complete degree up to p,
# #                 npoly = (PORDER+1)*(PORDER+1)*(PORDER+1)
# #      F:         Vandermonde matrix (npoints,npoly)
# #      FX:        Vandermonde matrix for the derivative of the Koornwinder
# #                 polynomials w.r.t. x (npoints,npoly)
# #      FY:        Vandermonde matrix for the derivative of the Koornwinder
# #                 polynomials w.r.t. y (npoints,npoly)
# #      FZ:        Vandermonde matrix for the derivative of the Koornwinder
# #                 polynomials w.r.t. z (npoints,npoly)
# #

n,dim = size(x);

if dim== 1 # 1D
        f,fx = legendrepolynomial(x,porder);  # Legendre basis
        fy=[0.0 0.0];
        fz=[0.0 0.0];
elseif dim == 2 # 2D
        g1,gx = legendrepolynomial(x[:,1],porder); # Legendre basis in x direction
        g2,gy = legendrepolynomial(x[:,2],porder); # Legendre basis in y direction
        f  = zeros(Float64,n,(porder+1)*(porder+1));
        fx = 0*f;
        fy = 0*f;
        fz=[0.0 0.0];
        # perform tensor product to obtain the shape functions and their
        # derivatives on the unit square
        for ii=1:n
            tm = kron(g2[ii,:],g1[ii,:]);
            f[ii,:] =  kron(g2[ii,:],g1[ii,:]);
            fx[ii,:] = kron(g2[ii,:],gx[ii,:]);
            fy[ii,:] = kron(gy[ii,:],g1[ii,:]);
        end
elseif dim==3
        g1,gx=legendrepolynomial(x[:,1],porder); # Legendre basis in x direction
        g2,gy=legendrepolynomial(x[:,2],porder); # Legendre basis in y direction
        g3,gz=legendrepolynomial(x[:,3],porder); # Legendre basis in z direction
        f  = zeros(Float64,n,(porder+1)*(porder+1)*(porder+1));
        fx = 0*f;
        fy = 0*f;
        fz = 0*f;
        # perform tensor product to obtain the shape functions and their
        # derivatives on the unit cube
        for ii=1:n
            f[ii,:] =  kron(g3[ii,:],kron(g2[ii,:],g1[ii,:]));
            fx[ii,:] = kron(g3[ii,:],kron(g2[ii,:],gx[ii,:]));
            fy[ii,:] = kron(g3[ii,:],kron(gy[ii,:],g1[ii,:]));
            fz[ii,:] = kron(gz[ii,:],kron(g2[ii,:],g1[ii,:]));
        end
else
        error("Only can handle dim=1, dim=2 or dim=3");
end

return f, fx, fy, fz;

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
module PointerBasedTrees
using FMMTrees

mutable struct Node{T}
    data::T
    num_children::Int
    next_sibling::Int
    parent::Int
    first_child::Int
end

data(n::Node) = n.data

abstract type APBTree end
struct PointerBasedTree{N<:Node} <: APBTree
    nodes::Vector{N}
    root::Int
end

FMMTrees.root(tree::PointerBasedTree) = tree.root

getnode(tree::PointerBasedTree, node_idx) = tree.nodes[node_idx]
nextsibling(tree::PointerBasedTree, node_idx) = getnode(tree, node_idx).next_sibling
parent(tree::PointerBasedTree, node_idx) = getnode(tree, node_idx).parent
firstchild(tree::PointerBasedTree, node_idx) = getnode(tree, node_idx).first_child

struct ChildView{T<:APBTree,N}
    tree::T
    node_idx::N
end

Base.iterate(itr::ChildView) = iterate(itr, firstchild(itr.tree, itr.node_idx))
function Base.iterate(itr::ChildView, state)
    # Check if state is a valid node pointer
    state < 1 && return nothing
    # If yes, check that the parent of the corresponding node equals
    # the original parent. If not, we are looking at a distant *cousin*
    sibling_par = parent(itr.tree, state)
    sibling_par != itr.node_idx && return nothing
    sibling_idx = nextsibling(itr.tree, state)
    return (state, sibling_idx)
end

Base.IteratorSize(cv::ChildView) = Base.SizeUnknown()

FMMTrees.children(tree::APBTree, node=FMMTrees.root(tree)) = ChildView(tree, node)
FMMTrees.haschildren(tree::APBTree, node) = (firstchild(tree,node) >= 1)
FMMTrees.data(tree::PointerBasedTree, node=FMMTrees.root(tree)) = data(getnode(tree, node))

"""
    insert!(tree, parent, data)

Insert a node carrying `data` as the first child of `parent`
"""
function FMMTrees.insert!(tree::PointerBasedTree, parent_idx, data)
    next = firstchild(tree, parent_idx)
    FMMTrees.insert!(tree, data, parent=parent_idx, next=next, prev=0)
end

setfirstchild!(node::Node, child) = Node(node.data, node.num_children, node.next_sibling, node.parent, child)
setfirstchild!(tree::PointerBasedTree, node, child) = tree.nodes[node] = setfirstchild!(getnode(tree, node), child)

setnextsibling!(node::Node, next) = Node(node.data, node.num_children, next, node.parent, node.first_child)
setnextsibling!(tree::PointerBasedTree, node, next) = tree.nodes[node] = setnextsibling!(getnode(tree, node), next)

"""
    insert!(tree, data, before=bef, parent=par, prev=prev)

Insert a new node carrying data such that the next sibling of this new node
will be `bef`. The parent node and the previous node to `bef` are required
to update all links in the data structure. If `bef` is the first child of `par`,
the previous node `prev` should be pasesed as `0`. Similarly, if the new node
should be inserted as the last child, `bef` should be passed as `0`. In the
special case of a parent node without any existing children, both `bef` and
`prev` should be set to `0`.

It may seem redundant to pass also `prev`, but this is a necessity of the single
linked list implementation of `PointerBasedTree`.
"""
function FMMTrees.insert!(tree::PointerBasedTree, data; parent, next, prev)
    push!(tree.nodes, Node(data, 0, next, parent, 0))
    @assert !(parent < 1)

    fs = firstchild(tree, parent)
    if fs < 1 || fs == next
        # getnode(tree, parent).first_child = length(tree.nodes)
        setfirstchild!(tree, parent, length(tree.nodes))
    end
    if !(prev < 1)
        # getnode(tree, prev).next_sibling = length(tree.nodes)
        setnextsibling!(tree, prev, length(tree.nodes))
    end
    return length(tree.nodes)
end

end


text:
abstract type AdversarialMDP{S,A} <: MDP{S,A} end
abstract type AdversarialPOMDP{S,A,O} <: POMDP{S,A,O} end

struct AdditiveAdversarialMDP{S,A}  <: AdversarialMDP{S,A}
    mdp::MDP{S,A}
    x_distribution
end

POMDPs.initialstate(mdp::AdditiveAdversarialMDP) = initialstate(mdp.mdp)

POMDPs.actions(mdp::AdditiveAdversarialMDP) = actions(mdp.mdp)
POMDPs.actionindex(mdp::AdditiveAdversarialMDP, a) = actionindex(mdp.mdp, a)

disturbances(mdp::AdditiveAdversarialMDP) = mdp.x_distribution.vs
disturbanceindex(mdp::AdditiveAdversarialMDP, x) = findfirst(mdp.x_distribution.vs .== x)

POMDPs.isterminal(mdp::AdditiveAdversarialMDP, s) = isterminal(mdp.mdp, s)
POMDPs.discount(mdp::AdditiveAdversarialMDP) = discount(mdp.mdp)

render(mdp::AdditiveAdversarialMDP, s, a = nothing; kwargs...) = render(mdp.mdp, s, a; kwargs...)
render(mdp::AdditiveAdversarialMDP; kwargs...) = render(mdp.mdp; kwargs...)

POMDPs.gen(mdp::AdditiveAdversarialMDP, s, a, x, rng::AbstractRNG = Random.GLOBAL_RNG; kwargs...) = gen(mdp.mdp, s, a .+ x, rng; kwargs...)

POMDPs.gen(mdp::AdditiveAdversarialMDP, s, a, rng::AbstractRNG = Random.GLOBAL_RNG; kwargs...) = gen(mdp, s, a, typeof(a)(rand(mdp.x_distribution)), rng; kwargs...)

isfailure(mdp::AdditiveAdversarialMDP, s) = isfailure(mdp.mdp, s)

struct AdditiveAdversarialPOMDP{S,A,O}  <: AdversarialPOMDP{S,A, O}
    pomdp::POMDP{S,A,O}
    x_distribution
end

POMDPs.initialstate(pomdp::AdditiveAdversarialPOMDP) = initialstate(pomdp.pomdp)
POMDPs.initialobs(pomdp::AdditiveAdversarialPOMDP, s) = initialobs(pomdp.pomdp, s)

POMDPs.isterminal(pomdp::AdditiveAdversarialPOMDP, s) = isterminal(pomdp.pomdp, s)
POMDPs.discount(pomdp::AdditiveAdversarialPOMDP) = discount(pomdp.pomdp)

render(pomdp::AdditiveAdversarialPOMDP, s, a = nothing; kwargs...) = render(pomdp.pomdp, s, a; kwargs...)
render(pomdp::AdditiveAdversarialPOMDP; kwargs...) = render(pomdp.pomdp; kwargs...)

POMDPs.gen(pomdp::AdditiveAdversarialPOMDP, s, a, x, rng::AbstractRNG = Random.GLOBAL_RNG; kwargs...) = gen(pomdp.pomdp, s, a .+ x, rng; kwargs...)

POMDPs.gen(pomdp::AdditiveAdversarialPOMDP, s, a, rng::AbstractRNG = Random.GLOBAL_RNG; kwargs...) = gen(pomdp, s, a, typeof(a)(rand(pomdp.x_distribution)), rng; kwargs...)

isfailure(pomdp::AdditiveAdversarialPOMDP, s) = isfailure(pomdp.pomdp, s)

max_forks_repo_name: ancorso/POMDPGym


text:
verify = "verify" in ARGS

struct Benchmark
    name::String
    benchmark::Union{Number, String}
    verify::Union{Number, String}
end

const BENCHMARKS = [
   #Benchmark("binarytrees", 21, 10),
  #  Benchmark("fannkuchredux", 12, 7),
   # Benchmark("fasta", 25000000, 1000),
   # Benchmark("knucleotide", 25000000, "knucleotide-input.txt"),
  #  Benchmark("mandelbrot", 16000, 200),
   # Benchmark("nbody", 50000000, 1000),
  #  Benchmark("pidigits", 10000, 27),
 #   Benchmark("regexredux", 5000000, "regexredux-input.txt"),
 #   Benchmark("revcomp", "0 < revcomp-input.txt"),
    Benchmark("spectralnorm", 5500, 100),
]

verify && println("VERIFYING!")
for benchmark in BENCHMARKS
    dir = benchmark.name
    _arg = verify ? benchmark.verify : benchmark.benchmark
    arg = _arg isa String ? "0 < $(_arg)" : string(_arg)
    println("Running $dir")
    bdir = joinpath(@__DIR__, dir)
    for file in readdir(bdir)
        endswith(file, ".jl") || continue
        println("    $file:")
        cmd = `$(Base.julia_cmd()) $(joinpath(bdir, file)) $(arg)`
        if verify
            bench_output = read(cmd, String)
            correct_output = read(joinpath(bdir, string(dir, "-output.txt")), String)
            if bench_output != correct_output
                @show bench_output
                @show correct_output
                error()
            end
        else
            time = @elapsed run(cmd)
            println("Time: $time")
        end
    end
end


text:
"""
    getNodeAges(net)

vector of node ages in pre-order, as in `nodes_changed`,
which is assumed to have been calculated before.
"""
function getNodeAges(net::HybridNetwork)
    x = Vector{Float64}(undef, length(net.nodes_changed))
    for i in reverse(1:length(net.nodes_changed)) # post-order
        n = net.nodes_changed[i]
        if n.leaf
            x[i] = 0.0
            continue
        end
        for e in n.edge
            if getParent(e) == n # n parent of e
                childnode = getChild(e)
                childIndex = getIndex(childnode, net.nodes_changed)
                x[i] = x[childIndex] + e.length
                break # use 1st child, ignores all others
            end
        end
    end
    return x
end

"""
    pairwiseTaxonDistanceMatrix(net; keepInternal=false,
                                checkPreorder=true, nodeAges=[])
    pairwiseTaxonDistanceMatrix!(M, net, nodeAges)

Return the matrix `M` of pairwise distances between nodes in the network:
- between all nodes (internal and leaves) if `keepInternal=true`,
  in which case the nodes are listed in `M` in the
  order in which they appear in `net.nodes_changed`
- between taxa only otherwise, in which case the nodes are listed
  in `M` in the order in which they appear in `tipLabels(net)`
  (i.e. same order as in `net.leaf`)

The second form modifies `M` in place, assuming all nodes.

The distance between the root and a given hybrid node (to take an example)
is the weighted average of path lengths from the root to that node,
where each path is weighted by the product of γs of all edges on that path.
This distance measures the average genetic distance across the genome,
if branch lengths are in substitutions/site.

optional arguments:

- `checkPreorder`: if true, `net.nodes_changed` is updated to get a
  topological ordering of nodes.
- `nodeAges`: if not provided, i.e. empty vector, the network is *not* modified.
  If provided and non-empty, `nodeAges` should list node ages in the
  pre-order in which nodes are listed in `nodes_changed` (including leaves),
  and **edge lengths** in `net` **are modified** accordingly.

Providing node ages hence makes the network time consistent: such that
all paths from the root to a given hybrid node have the same length.
If node ages are not provided, the network need not be time consistent.
"""
function pairwiseTaxonDistanceMatrix(net::HybridNetwork;
            keepInternal=false::Bool, checkPreorder=true::Bool,
            nodeAges=Float64[]::Vector{Float64})
    net.isRooted || error("net needs to be rooted for preorder recursion")
    if(checkPreorder)
        preorder!(net)
    end
    nnodes = net.numNodes
    M = zeros(Float64,nnodes,nnodes)
    if length(nodeAges)>0
        length(nodeAges) == net.numNodes ||
          error("there should be $(net.numNodes) node ages")
    end
    pairwiseTaxonDistanceMatrix!(M,net,nodeAges)
    if !keepInternal
        M = getTipSubmatrix(M, net)
    end
    return M
end
"""
    getTipSubmatrix(M, net; indexation=:both)

Extract submatrix of M, with rows and/or columns corresponding to
tips in the network, ordered like in `net.leaf`.
In M, rows and/or columns are assumed ordered as in `net.nodes_changed`.

indexation: one of `:rows`, `:cols` or `:both`: are nodes numbers indexed
in the matrix by rows, by columns, or both? Subsetting is taken accordingly.

"""
function getTipSubmatrix(M::Matrix, net::HybridNetwork; indexation=:both)
    nodenames = [n.name for n in net.nodes_changed]
    tipind = Int[]
    for l in tipLabels(net)
        push!(tipind, findfirst(isequal(l), nodenames))
    end
    if indexation == :both
        return M[tipind, tipind]
    elseif indexation == :rows
        return M[tipind, :]
    elseif indexation == :cols
        return M[:, tipind]
    else
        error("indexation must be one of :both :rows or :cols")
    end
end

function pairwiseTaxonDistanceMatrix!(M::Matrix{Float64},net::HybridNetwork,nodeAges)
    recursionPreOrder!(net.nodes_changed, M, # updates M in place
            updateRootSharedPathMatrix!, # does nothing
            updateTreePairwiseTaxonDistanceMatrix!,
            updateHybridPairwiseTaxonDistanceMatrix!,
            nodeAges)
end

function updateTreePairwiseTaxonDistanceMatrix!(V::Matrix,
            i::Int,parentIndex::Int,edge::Edge,
            params)
    nodeAges = params # assumed pre-ordered, as in nodes_changed
    if length(nodeAges)>0
        edge.length= nodeAges[parentIndex] - nodeAges[i]
    end
    for j in 1:(i-1)
        V[i,j] = V[parentIndex,j]+edge.length
        V[j,i] = V[i,j]
    end
    V[i,i] = 0.0
end

function updateHybridPairwiseTaxonDistanceMatrix!(V::Matrix,
        i::Int, parentIndex1::Int, parentIndex2::Int,
        edge1::Edge, edge2::Edge,
        params)
    nodeAges = params # should be pre-ordered
    if length(nodeAges)>0
        edge1.length= nodeAges[parentIndex1] - nodeAges[i]
        edge2.length= nodeAges[parentIndex2] - nodeAges[i]
    end
    for j in 1:(i-1)
        V[i,j] = edge1.gamma*(edge1.length+V[parentIndex1,j]) +
                 edge2.gamma*(edge2.length+V[parentIndex2,j])
        V[j,i] = V[i,j]
    end
end

"""
    pairwiseTaxonDistanceGrad(net; checkEdgeNumber=true, nodeAges=[])

3-dim array: gradient of pairwise distances between all nodes.
(internal and leaves); gradient with respect to edge lengths
if `nodeAges` is empty; with respect to node ages otherwise.
Assume correct `net.nodes_changed` (preorder).
This gradient depends on the network's topology and γ's only,
not on branch lengths or node ages (distances are linear in either).

WARNING: edge numbers need to range between 1 and #edges.
"""
function pairwiseTaxonDistanceGrad(net::HybridNetwork;
        checkEdgeNumber=true::Bool, nodeAges=Float64[]::Vector{Float64})
    if checkEdgeNumber
      sort([e.number for e in net.edge]) == collect(1:net.numEdges) ||
        error("edge numbers must range between 1 and #edges")
    end
    n = (length(nodeAges)==0 ? net.numEdges : net.numNodes)
    M = zeros(Float64, net.numNodes, net.numNodes, n)
    recursionPreOrder!(net.nodes_changed, M,
            updateRootSharedPathMatrix!,  # does nothing
            updateTreePairwiseTaxonDistanceGrad!,
            updateHybridPairwiseTaxonDistanceGrad!,
            nodeAges) # nodeAges assumed pre-ordered, like nodes_changed
    return M
end
function updateTreePairwiseTaxonDistanceGrad!(V::Array{Float64,3}, i::Int,
            parentIndex::Int, edge::Edge, params)
    nodeAges = params # assumed pre-ordered
    for j in 1:(i-1)
        if length(nodeAges) == 0 # d/d(edge length)
            V[i,j,edge.number] = 1.0
        else                     # d/d(node age)
            V[i,j,parentIndex] = 1.0
            V[i,j,i] = -1.0
        end
        for k in 1:size(V)[3] # brute force...
            V[i,j,k] += V[parentIndex,j,k]
            V[j,i,k] = V[i,j,k]
        end
    end
    # V[i,i,k] initialized to 0.0 already
end
function updateHybridPairwiseTaxonDistanceGrad!(V::Array{Float64,3},i::Int,
            parentIndex1::Int, parentIndex2::Int,
            edge1::Edge, edge2::Edge, params)
    nodeAges = params
    for j in 1:(i-1)
        if length(nodeAges) == 0 # d/d(edge length)
            V[i,j,edge1.number] = edge1.gamma
            V[i,j,edge2.number] = edge2.gamma
        else                     # d/d(node age)
            V[i,j,parentIndex1] = edge1.gamma
            V[i,j,parentIndex2] = edge2.gamma
            V[i,j,i] = - 1.0 # because γ1+γ2 = 1
        end
        for k in 1:size(V)[3]
          V[i,j,k] += edge1.gamma*V[parentIndex1,j,k] + edge2.gamma*V[parentIndex2,j,k]
          V[j,i,k] = V[i,j,k]
        end
    end
end

"""
    calibrateFromPairwiseDistances!(net, distances::Matrix{Float64},
        taxon_names::Vector{<:AbstractString})

Calibrate the network to match (as best as possible) input
pairwise distances between taxa, such as observed from sequence data.
`taxon_names` should provide the list of taxa, in the same order
in which they they are considered in the `distances` matrix.
The optimization criterion is the sum of squares between the
observed distances, and the distances from the network
(weighted average of tree distances, weighted by γ's).
The network's edge lengths are modified.

Warning: for many networks, mutiple calibrations can fit the pairwise
distance data equally well (lack of identifiability).
This function will output *one* of these equally good calibrations.

optional arguments (default):
- checkPreorder (true)
- forceMinorLength0 (false) to force minor hybrid edges to have a length of 0
- NLoptMethod (`:LD_MMA`) for the optimization algorithm.
  Other options include `:LN_COBYLA` (derivative-free); see NLopt package.
- tolerance values to control when the optimization is stopped:
  ftolRel (1e-12), ftolAbs (1e-10) on the criterion, and
  xtolRel (1e-10), xtolAbs (1e-10) on branch lengths / divergence times.
- verbose (false)
"""
function calibrateFromPairwiseDistances!(net::HybridNetwork,
      D::Array{Float64,2}, taxNames::Vector{<:AbstractString};
      checkPreorder=true::Bool, forceMinorLength0=false::Bool, verbose=false::Bool,
      ultrametric=true::Bool, NLoptMethod=:LD_MMA::Symbol,
      ftolRel=fRelBL::Float64, ftolAbs=fAbsBL::Float64,
      xtolRel=xRelBL::Float64, xtolAbs=xAbsBL::Float64)

    checkPreorder && preorder!(net)
    # fixit: remove root node if of degree 2, and if ultrametric=false
    defaultedgelength = median(D)/(length(net.edge)/2)
    for e in net.edge
        if e.length == -1.0 e.length=defaultedgelength; end
        # get smarter starting values: NJ? fast dating?
    end
    if ultrametric # get all node ages in pre-order
        na = getNodeAges(net)
    else na = Float64[]; end
    # get number and indices of edge/nodes to be optimized
    if forceMinorLength0 && !ultrametric
        parind = filter(i -> net.edge[i].isMajor, 1:net.numEdges)
        nparams = length(parind) # edges to be optimized: indices
        par = [e.length for e in net.edge][parind] # and lengths
        for i in 1:net.numEdges
            if !net.edge[i].isMajor net.edge[i].length=0.0; end
        end
    elseif !forceMinorLength0 && !ultrametric
        nparams = length(net.edge) # number of parameters to optimize
        par = [e.length for e in net.edge]
        parind = 1:nparams
    elseif !forceMinorLength0 && ultrametric
        nparams = net.numNodes - net.numTaxa # internal nodes to be optimized
        parind = filter(i -> !net.nodes_changed[i].leaf, 1:net.numNodes)
        par = na[parind]
    else # forceMinorLength0 && ultrametric
        nparams = net.numNodes - net.numTaxa - net.numHybrids
        parind = filter(i -> !(net.nodes_changed[i].leaf || net.nodes_changed[i].hybrid),
                        1:net.numNodes)
        par = na[parind]
        hybInd = filter(i -> net.nodes_changed[i].hybrid, 1:net.numNodes)
        hybParentInd = Int[] # index in nodes_changed of minor parent
        hybGParentI = Int[] # index in 1:nparams of minor (grand-)parent in param list
        for i in hybInd
            n = net.nodes_changed[i]
            p = getMinorParent(n)
            pi = findfirst(n -> n===p, net.nodes_changed)
            push!(hybParentInd, pi)
            pii = findfirst(isequal(pi), parind)
            while pii===nothing # in case minor parent of n is also hybrid node
                p = getMinorParent(p)
                pi = findfirst(n -> n===p, net.nodes_changed)
                pii = findfirst(isequal(pi), parind)
            end
            push!(hybGParentI, pii)
        end
    end
    # initialize M=dist b/w all nodes, G=gradient (constant)
    M = pairwiseTaxonDistanceMatrix(net, keepInternal=true,
            checkPreorder=false, nodeAges=na)
    if !ultrametric && sort([e.number for e in net.edge]) != collect(1:net.numEdges)
        for i in 1:net.numEdges # renumber edges, needed for G
            net.edge[i].number = i
        end
    end # G assumes edges numbered 1:#edges, if optim edge lengths
    G = pairwiseTaxonDistanceGrad(net, checkEdgeNumber=false, nodeAges=na) .* 2
    # match order of leaves in input matrix, versus pre-order
    nodenames = [n.name for n in net.nodes_changed] # pre-ordered
    ntax = length(taxNames)
    tipind = Int[] # pre-order index for leaf #i in dna distances
    for l in taxNames
        i = findfirst(isequal(l), nodenames)
        i !== nothing || error("taxon $l not found in network")
        push!(tipind, i)
    end
    # contraints: to force a parent to be older than its child
    if ultrametric
      numConstraints = length(parind) -1 + net.numHybrids
      # calculate indices in param list of child & (grand-)parent once
      chii = Int[] # non-root internal node, can be repeated: once per constraint
      anii = Int[] # closest ancestor in param list
      ci = 1 # index of constraint
      for i in 2:length(net.nodes_changed) # 1=root, can skip
        n = net.nodes_changed[i]
        if n.leaf continue; end # node ages already bounded by 0
        if n.hybrid && forceMinorLength0          # get index in param list of
          nii = hybGParentI[findfirst(isequal(i), hybInd)] # minor grand-parent (same age)
        else
          nii = findfirst(isequal(i), parind)
        end
        for e in n.edge
          if getChild(e) == n # n child of e
            p = getParent(e)  # parent of n
            if forceMinorLength0 && n.hybrid && !e.isMajor
                continue; end # p and n at same age already
            pi = findfirst(isequal(p.number), [no.number for no in net.nodes_changed])
            if forceMinorLength0 && p.hybrid
              pii = hybGParentI[findfirst(isequal(pi), hybInd)]
            else
              pii = findfirst(isequal(pi), parind)
            end
            push!(chii, nii)
            push!(anii, pii)
            verbose && println("node $(net.nodes_changed[parind[nii]].number) constrained by age of parent $(net.nodes_changed[parind[pii]].number)")
            ci += 1
          end
        end
      end
      length(chii) == numConstraints ||
        error("incorrect number of node age constraints: $numConstraints")
      function ageConstraints(result, nodeage, grad)
        if length(grad) > 0 # grad: nparams x nConstraints: ∂cj/∂xi = grad[i,j]
            fill!(grad, 0.0)
        end
        for j in 1:numConstraints
          nii = chii[j]; pii = anii[j]
          result[j] = nodeage[nii] - nodeage[pii] # jth constraint: cj ≤ 0
          if length(grad) > 0 # nparams x nConstraints
            grad[nii, j] =  1.
            grad[pii, j] = -1.
          end
        end
      end
    end
    opt = NLopt.Opt(NLoptMethod,nparams) # :LD_MMA to use gradient
    # :LN_COBYLA for (non)linear constraits, :LN_BOBYQA for bound constraints
    NLopt.maxeval!(opt,1000) # max iterations
    NLopt.ftol_rel!(opt,ftolRel)
    NLopt.ftol_abs!(opt,ftolAbs)
    NLopt.xtol_rel!(opt,xtolRel)
    NLopt.xtol_abs!(opt,xtolAbs)
    # NLopt.maxtime!(opt, t::Real)
    NLopt.lower_bounds!(opt, zeros(nparams))
    if ultrametric
      NLopt.inequality_constraint!(opt,ageConstraints,fill(0.0,numConstraints))
    end
    counter = [0]
    function obj(x::Vector{Float64}, grad::Vector{Float64})
        verbose && println("mismatch objective, BL = $(x)")
        counter[1] += 1
        # update edge lengths or node ages
        if ultrametric # update na using x, in place
            for i in 1:nparams # na=0 at leaves already
                na[parind[i]] = x[i]
            end
            if forceMinorLength0 # set hybrid age to minor parent age
                for i in 1:net.numHybrids
                    na[hybInd[i]] = na[hybParentInd[i]] # pre-order important
                end
            end
        else # not ultrametric: optimize branch lengths
            for i in 1:nparams # update network
                net.edge[parind[i]].length = x[i]
            end
        end
        # update distances in M, in place
        pairwiseTaxonDistanceMatrix!(M,net,na)
        ss = 0.0 # sum of squares between M and observed distances
        for i in 2:ntax; for j in 1:(i-1)
            ss += (M[tipind[i],tipind[j]]-D[i,j])^2
        end; end
        if length(grad) > 0 # sum_ij 2 * dM_ij/dx_t * (M_ij-D_ij)
          for t in 1:nparams grad[t] = 0.0; end;
          for i in 2:ntax; for j in 1:(i-1);
            for t in 1:nparams
              G[tipind[i],tipind[j],parind[t]] != 0.0 || continue
              grad[t] += G[tipind[i],tipind[j],parind[t]] *
                        (M[tipind[i],tipind[j]]-D[i,j])
            end
            if ultrametric && forceMinorLength0
              for i in 1:net.numHybrids # na[hybrid] set to na[minor parent]
                grad[hybGParentI[i]] += G[tipind[i],tipind[j],hybInd[i]] *
                        (M[tipind[i],tipind[j]]-D[i,j])
              end
            end
          end; end
        end
        return ss
    end
    NLopt.min_objective!(opt,obj)
    fmin, xmin, ret = NLopt.optimize(opt,par) # optimization here!
    verbose && println("got $(round(fmin, digits=5)) at $(round.(xmin, digits=5)) after $(counter[1]) iterations (return code $(ret))")
    return fmin,xmin,ret
end


text:
module Get_Sim_Params
    using Constants
    import YAML

    path = ARGS[2] #"./runs/run01/run01.yml"
    data = YAML.load(open(path))

    export geom, nz, nregions, L, smoothing
    export coeff_CFL, maxTime, smoothing, dr_smooth, tm_quiet
    export dt_print, efield_switch, electron_switch
    export heattransfer_switch, flimit, electron_heat_conduction
    export friction_switch, Lab_min, check_frequency, smax, ion_viscosity
    export bremsstrahlung, restart, nedges
    export nspec, Ai, Zi, r0, den0, temp0, vel0
    export rmin, neqi, eps_visc_max, eps_compress
    export smooth_vars, smooth_type
    export boundary, phi

    key = "grid"
    const geom = data[key]["geom"]
    const nz = data[key]["nz"]
    const nregions = data[key]["nregions"]
    const L = data[key]["L"]
    const boundary = data[key]["boundary"]

    key = "output"
    const dt_print = data[key]["dt_print"]

    key = "numerics"
    const restart = data[key]["restart"]
    const tm_quiet = data[key]["tm_quiet"]
    const coeff_CFL = data[key]["coeff_CFL"]
    const maxTime = data[key]["maxTime"]
    const smoothing = data[key]["smoothing"]
    const dr_smooth = data[key]["dr_smooth"]
    const eps_visc_max = data[key]["eps_visc_max"]
    const eps_compress = data[key]["eps_compress"]
    const nedges = data[key]["nedges"]
    const smooth_vars = data[key]["smooth_vars"]
    const phi = data[key]["phi"]

    key = "physics"
    subkey = "electric_field"
    const efield_switch = data[key][subkey]["efield_switch"]
    const electron_switch = data[key][subkey]["electron_switch"]
    subkey = "heat_conduction"
    const heattransfer_switch = data[key][subkey]["heattransfer_switch"]
    const flimit = data[key][subkey]["flimit"]
    const electron_heat_conduction = data[key][subkey]["electron_heat_conduction"]
    subkey = "collisions"
    const friction_switch = data[key][subkey]["friction_switch"]
    const Lab_min = data[key][subkey]["Lab_min"]
    const check_frequency = data[key][subkey]["check_frequency"]
    const smax = data[key][subkey]["smax"]
    const ion_viscosity = data[key][subkey]["ion_viscosity"]
    subkey = "radiation"
    const bremsstrahlung = data[key][subkey]["bremsstrahlung"]

    key = "species"
    const nspec = data[key]["nspec"]
    ai = zeros(Float64, nspec)
    zi = zeros(Float64, nspec)
    r = zeros(Float64, nspec, nregions + 1)
    den = zeros(Float64, nspec, nregions)
    temp = zeros(Float64, nspec, nregions)
    vel = zeros(Float64, nspec, nregions)
    smooth_type = Array{String}(nedges, nspec)
    for ispec = 1:nspec
        ai[ispec] = data[key]["species"*string(ispec)]["Ai"]
        zi[ispec] = data[key]["species"*string(ispec)]["Zi"]
        r_species = data[key]["species"*string(ispec)]["r"]
        for iedge = 1:nregions + 1
            if isa(r_species[iedge], Float64)
                r[ispec, iedge] = r_species[iedge]
            else
                r[ispec, iedge] = L
            end
        end
        den_species = data[key]["species"*string(ispec)]["den"]
        temp_species = qe_C * data[key]["species"*string(ispec)]["temp"]
        vel_species = data[key]["species"*string(ispec)]["vel"]
        for ireg = 1:nregions
            den[ispec, ireg] = den_species[ireg]
            temp[ispec, ireg] = temp_species[ireg]
            vel[ispec, ireg] = vel_species[ireg]
        end
        smooth_type_species = data[key]["species"*string(ispec)]["smooth_type"]
        for iedge = 1:nedges
            smooth_type[iedge, ispec] = smooth_type_species
        end

    end
    const Ai = ai
    const Zi = zi
    const r0 = r
    const den0 = den
    const temp0 = temp
    const vel0 = vel

    if geom=="slab"
        const rmin = 0.
    else #spherical
        const rmin = 2.5e-6
    end

    const neqi = 3*nspec
end

max_stars_repo_name: cbellei/JUPL

max_issues_repo_name: cbellei/JUPL

max_forks_repo_name: cbellei/JUPL


text:
mutable struct Historied
    num::Number
    history::Vector{Number}
    Historied(n) = new(n, Vector{Number}())
end

assign(y::Historied, z) = (push!(y.history, y.num); y.num = z; y)

x = Historied(1)

assign(x, 3)
assign(x, 5)
assign(x, 4)

println("Past history of variable x: $(x.history). Current value is $(x.num)")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode