text:
# ---
# title: CVI Simple Example
# id: cvi_example
# cover: ../assets/logo.png
# date: 2021-12-6
# author: "[Sasha Petrenko](https://github.com/AP6YC)"
# julia: 1.6
# description: This demo is a simple example of how to use a CVI in batch mode.
# ---

# ## Overview

# This demo is a simple example of how to use CVIs in batch mode.
# Here, we load a simple dataset and run a basic clustering algorithm to prescribe a set of clusters to the features.
# It is a combination of these features and the prescribed labels that are used to compute the criterion value.
# This simple example demonstrates the usage of a single CVI, but it may be substituted for any other CVI in the `ClusterValidityIndices.jl` package.

# ## Clustering

# ### Data Setup

# First, we must load all of our dependencies.
# We will load the `ClusterValidityIndices.jl` along with some data utilities and the Julia `Clustering.jl` package to cluster that data.
using ClusterValidityIndices    # CVI/ICVI
using Clustering                # DBSCAN
using MLDatasets                # Iris dataset
using MLDataUtils               # Shuffling and splitting
using Printf                    # Formatted number printing

# We will download the Iris dataset for its small size and benchmark use for clustering algorithms.
Iris.download(i_accept_the_terms_of_use=true)
features, labels = Iris.features(), Iris.labels()

# Because the MLDatasets package gives us Iris labels as strings, we will use the `MLDataUtils.convertlabel` method with the `MLLabelUtils.LabelEnc.Indices` type to get a list of integers representing each class:
labels = convertlabel(LabelEnc.Indices{Int}, labels)
unique(labels)

# ### Fuzzy C-Means

# Get the Fuzzy C-Means clustering result
results = fuzzy_cmeans(features, 3, 2)

# Because the results are fuzzy weights, find the maximum elements along each sample
indices = argmax(results.weights, dims=2)

# Get those labels as a vector of integers
c_labels = vec([c[2] for c in indices])

# ## CVI Criterion Value Extraction

# Now that we have some data and a clustering algorithm's prescribed labels, we can compute a criterion value using a CVI in batch mode.
# First, we create a CVI object with the default constructor:

## Create a CVI object
my_cvi = CH()

# Finally we can simply get the criterion value in batch by passing all of the data and Fuzzy C-Means labels at once.

## Get the batch criterion value
criterion_value = get_cvi!(my_cvi, features, c_labels)


text:
function Tangent_to_NamedTuple(t)
  return NamedTuple((k => t[k] for k in keys(t)))
end

Tangent_to_NamedTuple(::ZeroTangent) = ZeroTangent()

function rrule(::Type{Op}, x1, x2, x3)
  y = Op(x1, x2, x3)
  function Op_pullback(ȳ)
    x̄1 = x1
    x̄2 = x2
    t = ȳ.params
    x̄3 = Tangent_to_NamedTuple(t)
    return (NoTangent(), x̄1, x̄2, x̄3)
  end
  return y, Op_pullback
end

@non_differentiable Ops.sites(::Any)

ext: jl


text:
# THE BASIC DIFFERENTIATION RULES
# https://www.mathbootcamps.com/power-rule-derivatives/ = power-rule-derivatives
# https://www.dummies.com/education/math/calculus/the-basic-differentiation-rules/ for dummies

# The constant rule
 f(x) = 5 # 5 is a horizontal line with a slope of zero, and thus its derivative is also zero

# The power rule
 f(x) = f^5
# To find its derivative take the 5 power, pull it to the front and reduce the power by 1 thus
  f(x) = 5f^4 # derivative

# Negative power
  f(x) = x^-2
# negative power derivative
 f(x) = 2x^-3 # pull power in front and reduce power by 1 , this ^-2 becomes ^-3

# find the derivative of the function
 y = 2x^4 - 5x^2 + 1
# enclose in brackets
 y' = (2x^4)' - (5x^2)' + (1)'
# factor out coefficients
 y' = 2(x^4)' - 5(x^2)' + (1)'
# apply the power rule
y' = 2(4x^3)' - 5(2x^1)' + (0) # note the constant rule of 1, 1 is a line with a slope of 0 and derivative also 0
# answer is to now simplify
= 8x^3 - 10x

# Find the derivative of the function
y = 5x^3 - 3x^2 + 10x - 8
# enclose in brackets
y' = (5x^3)' - (3x^2)' + (10x)' - (8)'
# factor out the coefficients
y' = 5(x^3)' - 3(x^2)' + 10(x^0)' - (8)' # note - since x is byitself we can write 1x^0 anything > 0 to the power of 0 is 1
# apply the power rule
y' = 5(3x^2)' - 3(2x^1)' + 10(x^0)' - (0)' # note the constant rule of 8, 8 is a line with a slope of 0 and derivative also 0
# now simplify
= 15x^2 - 6x + 10x

# find the derivative of the function
y = 2/x^4 - 1/x^2

# apply rule of negative exponents
y = 2x^-4 - x^-2

# enclose in brackets
y = (2x^-4) - (x^-2)

# factor out the coefficients
y = 2(x^-4) - (x^-2)

# power rule
y = 2(4x^-4-1) - (2x^-2-1) # OR reduce by 1 = 2(4x^-5) - (2x^-3)

# simplify
 -8x^-5 + (2x^-3)
 = 8 / x^5 + 2 / x^3 # only negative terms go to the bottom

 # Derivatives of functions with radicals (square roots and other roots)
 # You can differentiate radical functions by rewriting them as power functions and then using the power rule.
 f(x) = 3 sqrt(x^2)
 # rewrite as
 f(x) = x^2/3

 # find the derivative of the function
 y = 4 sqrt(x) -6 3 sqrt(x^2)
# rewrite function using algebra
y = 4x^1/2 - 6x^2/3
# enclose in brackets
y = (4x^1/2) - (6x^2/3)
# factor out the coefficients
y = 4(x^1/2) - 6(x^2/3)
# power rule
y = 4(1/2x^1/2-1) - 6(2/3x^2/3-1) # note ^2/3 becomes 1/3 as subtract from numerator
y = 2x^-1.5 - 4x^-1/3
# convert to fraction pulling negative exponents to the bottom
y = 2 / x^1.5 - 4 / x^1/3

= 2 / sqrt(x) - 4 / 3 sqrt(x)

# power rule is then used in the chain rule and quotient rule etc...

# https://www.mathbootcamps.com/the-product-rule/ - Product rule

# tomorrow 5.14.2019

# https://www.mathbootcamps.com/the-chain-rule/ - chain rule
# chain rule
f'(x) = ( g(h(x)))' = g'(h(x)) h'(x)

# g = outside function
# h = inside function
""" to find the derivative of f(x) then you can do this by taking the derivative of outside function g and then multiplying it by its derivative of the function h.
In other words you are finding the derivative of f(x) by finding the derivative of its pieces
"""

# Examples of applying the chain rule
""" as we apply the chain rule we will always focus on figuring out what the outside and inside functions are
"""
# find the derivative of the function
f'(x) = (3x +1) ^5
# inside function = 3x + 1
# outside function = ^5

""" you know by the power rule that the derivative of x^5 = 5x^4
So, cover up that 3x+1, and pretend it is an x for a minute. The only deal is, you will have to pay a penalty. Since it was actually not just an x, you will have to multiply by the derivative of the 3x+1
"""

# Step 1 = keep the inside
# Step 2 = take the derivative of the outside
# step 3 = multiply inside by derivative of the outside
f'(x) = 5(3x +1)^4 (3x +1)'
f'(x) = 5(3x +1)^4 (3)

# simplyfy
f'(x) = 15(3x + 4)^4  # when multiply the inside by the outside = 5 * 3 = 15 then that leaves the x and x * 3 = 3x and 1 * 3 = 1

#=============================================
# final = https://medium.com/@pdquant/all-the-backpropagation-derivatives-d5275f727f60
#=============================================


text:
using PSCOPF

using Test
using Dates

@testset "test_modes" begin

    @testset "mode1" begin
        @test PSCOPF.PSCOPF_MODE_1.name == "mode_1"
        @test PSCOPF.get_fo_length(PSCOPF.PSCOPF_MODE_1) == Dates.Minute(60)
    end

    @testset "mode2" begin
        @test PSCOPF.PSCOPF_MODE_2.name == "mode_2"
        @test PSCOPF.get_fo_length(PSCOPF.PSCOPF_MODE_2) == Dates.Minute(60)
    end

    @testset "mode3" begin
        @test PSCOPF.PSCOPF_MODE_3.name == "mode_3"
        @test PSCOPF.get_fo_length(PSCOPF.PSCOPF_MODE_3) == Dates.Minute(60)
    end

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text: include("tux.jl")

max_stars_repo_name: milktrader/Tuxedo

max_issues_repo_name: milktrader/Tuxedo

max_forks_repo_name: milktrader/Tuxedo


text:
# function backward_pass!(p_data::PolicyData, m_data::ModelData)
#     T = m_data.T
#     fx = m_data.dyn_deriv.fx
#     fu = m_data.dyn_deriv.fu
#     gx = m_data.obj_deriv.gx
#     gu = m_data.obj_deriv.gu
#     gxx = m_data.obj_deriv.gxx
#     guu = m_data.obj_deriv.guu
#     gux = m_data.obj_deriv.gux
#
#     # policy
#     K = p_data.K
#     k = p_data.k
#
#     # value function approximation
#     P = p_data.P
#     p = p_data.p
#
#     # state-action value function approximation
#     Qx = p_data.Qx
#     Qu = p_data.Qu
#     Qxx = p_data.Qxx
#     Quu = p_data.Quu
#     Qux = p_data.Qux
#
#     # terminal value function
#     P[T] .= gxx[T]
#     p[T] .=  gx[T]
#
#     for t = T-1:-1:1
#         Qx[t] .= gx[t] + fx[t]' * p[t+1]
#         Qu[t] .= gu[t] + fu[t]' * p[t+1]
#         Qxx[t] .= gxx[t] + fx[t]' * P[t+1] * fx[t]
#         Quu[t] .= guu[t] + fu[t]' * P[t+1] * fu[t]
#         Qux[t] .= gux[t] + fu[t]' * P[t+1] * fx[t]
#
#         K[t] .= -1.0 * Quu[t] \ Qux[t]
#         k[t] .= -1.0 * Quu[t] \ Qu[t]
#
#         P[t] .=  Qxx[t] + K[t]' * Quu[t] * K[t] + K[t]' * Qux[t] + Qux[t]' * K[t]
#         p[t] .=  Qx[t] + K[t]' * Quu[t] * k[t] + K[t]' * Qu[t] + Qux[t]' * k[t]
#     end
# end

function backward_pass!(p_data::PolicyData, m_data::ModelData;
    mode = :nominal)

    T = m_data.T
    fx = m_data.dyn_deriv.fx
    fu = m_data.dyn_deriv.fu
    gx = m_data.obj_deriv.gx
    gu = m_data.obj_deriv.gu
    gxx = m_data.obj_deriv.gxx
    guu = m_data.obj_deriv.guu
    gux = m_data.obj_deriv.gux

    # policy
    if mode == :nominal
        K = p_data.K
        k = p_data.k
    else
        K = p_data.K_cand
        k = p_data.k_cand
    end

    # value function approximation
    P = p_data.P
    p = p_data.p

    # state-action value function approximation
    Qx = p_data.Qx
    Qu = p_data.Qu
    Qxx = p_data.Qxx
    Quu = p_data.Quu
    Qux = p_data.Qux

    # terminal value function
    P[T] .= gxx[T]
    p[T] .=  gx[T]

    for t = T-1:-1:1
        # Qx[t] .= gx[t] + fx[t]' * p[t+1]
        mul!(Qx[t], transpose(fx[t]), p[t+1])
        Qx[t] .+= gx[t]

        # Qu[t] .= gu[t] + fu[t]' * p[t+1]
        mul!(Qu[t], transpose(fu[t]), p[t+1])
        Qu[t] .+= gu[t]

        # Qxx[t] .= gxx[t] + fx[t]' * P[t+1] * fx[t]
        mul!(p_data.xx̂_tmp[t], transpose(fx[t]), P[t+1])
        mul!(Qxx[t], p_data.xx̂_tmp[t], fx[t])
        Qxx[t] .+= gxx[t]

        # Quu[t] .= guu[t] + fu[t]' * P[t+1] * fu[t]
        mul!(p_data.ux̂_tmp[t], transpose(fu[t]), P[t+1])
        mul!(Quu[t], p_data.ux̂_tmp[t], fu[t])
        Quu[t] .+= guu[t]

        # Qux[t] .= gux[t] + fu[t]' * P[t+1] * fx[t]
        mul!(p_data.ux̂_tmp[t], transpose(fu[t]), P[t+1])
        mul!(Qux[t], p_data.ux̂_tmp[t], fx[t])
        Qux[t] .+= gux[t]

        # K[t] .= -1.0 * Quu[t] \ Qux[t]
        # k[t] .= -1.0 * Quu[t] \ Qu[t]
		p_data.uu_tmp[t] .= Quu[t]
        LAPACK.potrf!('U', p_data.uu_tmp[t])
        K[t] .= Qux[t]
        k[t] .= Qu[t]
        LAPACK.potrs!('U', p_data.uu_tmp[t], K[t])
		LAPACK.potrs!('U', p_data.uu_tmp[t], k[t])
		K[t] .*= -1.0
		k[t] .*= -1.0

        # P[t] .=  Qxx[t] + K[t]' * Quu[t] * K[t] + K[t]' * Qux[t] + Qux[t]' * K[t]
        # p[t] .=  Qx[t] + K[t]' * Quu[t] * k[t] + K[t]' * Qu[t] + Qux[t]' * k[t]
		mul!(p_data.ux_tmp[t], Quu[t], K[t])

		mul!(P[t], transpose(K[t]), p_data.ux_tmp[t])
		mul!(P[t], transpose(K[t]), Qux[t], 1.0, 1.0)
		mul!(P[t], transpose(Qux[t]), K[t], 1.0, 1.0)
		P[t] .+= Qxx[t]

		mul!(p[t], transpose(p_data.ux_tmp[t]), k[t])
		mul!(p[t], transpose(K[t]), Qu[t], 1.0, 1.0)
		mul!(p[t], transpose(Qux[t]), k[t], 1.0, 1.0)
		p[t] .+= Qx[t]
    end
end

max_issues_repo_name: thowell/DirectMotionPlanning


text:
include("logger.jl")
const path = "./class_basket/intentional_invalid_magic.class"
function intentionally_invalid()
    bytes = read(path)
    bytes[1:4] = zeros(4)
    write(path, bytes)

    ScriptLog("made class with invalid magic (Should be 0xCAFEBABE, now contains zeros)")
end

max_stars_repo_name: TapVM/Tap

max_issues_repo_name: TapVM/Tap

max_forks_repo_name: TapVM/Tap


text:
abstract type AbstractSeldonianProblem end

struct SeldonianProblem{TF,TG} <: AbstractSeldonianProblem where {TF,TG}
    f::TF
    g::TG
end

function solve(prob::SeldonianProblem, candidate_search, initialsolution, Dtrain, Dsafety)
    candidate = candidate_search(prob.f, Dtrain, initialsolution)

    safetycheck = safetytest(prob.g, Dsafety, candidate)
    if safetycheck
        return candidate
    else
        return :NSF
    end
end

function safetytest(g, Dsafety, candidate)
    return g(Dsafety, candidate) ≤ 0
end

function safetytest(g::AbstractArray, Dsafety, candidate)
    return all(safetytest.(g, (Dsafety,), (candidate,)))
end


text:
# Definition of tree and leaves
mutable struct Leaf
  parent :: Union{Nothing,Leaf}
  depth  :: Int
  left   :: Union{Int,Leaf}
  right  :: Union{Int,Leaf}
end

Leaf() = Leaf(nothing, 0, 0, 0)

# Set depth of a leaf
function setDepth(l::Leaf,d::Int)
  l.depth = d
  (typeof(l.left) !== Int) && setDepth(l.left,d+1)
  (typeof(l.right) !== Int) && setDepth(l.right,d+1)
end

# Given a string, separate two leaves
function splitLeaves(s::String) :: Tuple{String,String}
  begin k = 0; i = 2; end
  while k != 0 || i == 2
    (s[i] == '[') && (k += 1)
    (s[i] == ']') && (k -= 1)
    i += 1
  end
  return s[2:i-1], s[i+1:end-1]
end

# Create a new leaf from a string
function newLeaf(parent::Union{Leaf,Nothing}, depth::Int, s::String) :: Union{Int,Leaf}
  (length(s) == 1) && return parse(Int,s) # Base case -> Integer

  # General case -> Leaf
  left, right = splitLeaves(s)
  res :: Leaf = Leaf(parent,depth,0,0)
  res.left    = newLeaf(res,depth+1,left)
  res.right   = newLeaf(res,depth+1,right)
  return res
end

newTree(s::String) = newLeaf(nothing,0,s)

# Magnitude of a leaf
magnitude(l::Leaf) :: Int = 3*magnitude(l.left) + 2*magnitude(l.right)
magnitude(n::Int)  :: Int = n

# Get left-most sub-leaf
function leftMost(l::Leaf) :: Leaf
  l2 = l
  while typeof(l2.left) == Leaf
    l2 = l2.left
  end
  return l2
end

# Get right-most sub-leaf
function rightMost(l::Leaf) :: Leaf
  l2 = l
  while typeof(l2.right) === Leaf
    l2 = l2.right
  end
  return l2
end

# Get left neighbor of a leaf
function leftNbor(l::Leaf) :: Union{Leaf,Nothing}
  l2 = l
  while l2.parent !== nothing && l2 === l2.parent.left
    l2 = l2.parent
  end
  if (l2.parent === nothing)
    return nothing
  else
    (typeof(l2.parent.left) === Int) ? (return l2.parent) : (return rightMost(l2.parent.left))
  end
end

# Get right neighbor of a leaf
function rightNbor(l::Leaf) :: Union{Leaf,Nothing}
  l2 = l
  while l2.parent !== nothing && l2 === l2.parent.right
    l2 = l2.parent
  end
  if (l2.parent === nothing)
    return nothing
  else
    (typeof(l2.parent.right) === Int) ? (return l2.parent) : (return leftMost(l2.parent.right))
  end
end

# Check if any leaf can be exploded recursively (left ones first).
# Returns true if some leaf has been exploded.
function explode(l::Leaf) :: Bool
  if (typeof(l.left) !== Int) || (typeof(l.right) !== Int)
    return ((typeof(l.left) !== Int) &&  explode(l.left)) || ((typeof(l.right) !== Int) &&  explode(l.right))
  elseif l.depth >= 4
    left = leftNbor(l)
    right = rightNbor(l)

    if (left  !== nothing)
      (typeof(left.right) === Int) ? (left.right += l.left) : (left.left += l.left)
    end
    if (right !== nothing)
      (typeof(right.left) === Int) ?  (right.left += l.right) : (right.right += l.right)
    end

    (l.parent.left === l) ? (l.parent.left = 0) : (l.parent.right = 0)
    return true
  end
  return false
end

# Check if any leaf can be split recursively (left ones first).
# Returns true if some leaf has been split.
function split(l::Leaf) :: Bool
  if (typeof(l.left) == Int) && (l.left > 9)
    l.left = Leaf(l,l.depth+1,l.left÷2,l.left÷2+l.left%2)
    return true
  elseif (typeof(l.left) == Leaf) && split(l.left)
    return true
  elseif (typeof(l.right) == Int) && (l.right > 9)
    l.right = Leaf(l,l.depth+1,l.right÷2,l.right÷2+l.right%2)
    return true
  elseif (typeof(l.right) == Leaf) && split(l.right)
    return true
  end
  return false
end

# Reduce a leave until none of its sub-leaves can be exploded or split
function reduce(l::Leaf)
  keepGoing = true
  while keepGoing
    keepGoing = explode(l) || split(l)
  end
end

# Add two leaves together
function add(l::Leaf, r::Leaf) :: Leaf
  res = Leaf(nothing,0,l,r)
  l.parent = res
  r.parent = res
  setDepth(res,0)
  reduce(res)
  return res
end

# Display Leaves like in the inputs
printLeaf(l::Int) = print(l)

function printLeaf(l::Leaf)
  (l.parent === nothing) && println("")
  print("[")
  printLeaf(l.left)
  print(",")
  printLeaf(l.right)
  print("]")
  (l.parent === nothing) && println("")
end

function readInput() :: Array{String}
  f = open("./inputs/day18.txt")
  return map( x -> String(x), eachline(f))
end

function problem1(data::Array{String}) :: Int
  t = newTree(data[1])
  for line in data[2:end]
    printLeaf(t)
    tnew = newTree(line)
    t = add(t,tnew)
  end
  printLeaf(t)
  return magnitude(t)
end

function problem2(data::Array{String}) :: Int
  maxMag = 0
  for num1 in data
    for num2 in data
      if (num1 !== num2)
        t = add(newTree(num1),newTree(num2))
        m = magnitude(t)
        (m > maxMag) && (maxMag = m)
      end
    end
  end
  return maxMag
end

println("Reading Data... ")
data = readInput()

println("Problem 1: ")
res = problem1(data)
println(res)

println("Problem 2: ")
res = problem2(data)
println(res)


text:
@testset "LSQR" begin

	@testset "small test" begin
		A   = [1 1; 0 -1; 1 0.]
		r1  = A[:,1]   # compatible
		r2  = randn(3) # incompatible
		r3  = [1;1;-1.] # in the nullpace of A'

		x1 = lsqr(A,r1,out=2,atol=1e-8,doBidiag=true,storeInterm=true)
		@test norm(x1[1][:,end]-[1;0])<1e-10
		@test x1[2]==1
		x2 = lsqr(A,r2,out=1,atol=1e-8)
		@test x2[2]==2
		x3 = lsqr(A,r3,out=2)
		@test all(x3[1].==0)
	end

	@testset "sparse random" begin
		A   = sprandn(100,10,.2)
		rhs = randn(100)
		xgt = full(A)\rhs
		Af = (flag,x,a=0.0,v=0.0) -> ((flag=='F') ? A*x +a*v : A'*x+a*v)

		x1  = lsqr(LinearOperator(A),rhs,atol=1e-10,btol=1e-10,condlim=1e5,maxIter=10,out=1)
		x2  = lsqr(A,rhs,atol=1e-10,btol=1e-10,condlim=1e5,maxIter=10)
		x3  = lsqr(Af,rhs,atol=1e-10,btol=1e-10,condlim=1e5,maxIter=10)

		@test norm(xgt-x1[1])/norm(xgt) < 1e-6
		@test norm(xgt-x2[1])/norm(xgt) < 1e-6
		@test norm(xgt-x3[1])/norm(xgt) < 1e-6
		@test norm(x1[1]-x2[1])/norm(x1[1]) < 1e-12
		@test norm(x2[1]-x3[1])/norm(x3[1]) < 1e-12
	end
	@testset "hilbert matrix" begin
		A   = matrixdepot("hilb",3,3)
		rhs = [1;0;0.]
		xgt = A\rhs

		x1  = lsqr(LinearOperator(A),rhs,atol=1e-10,btol=1e-10,condlim=1e5,maxIter=10)
		x2  = lsqr(A,rhs,atol=1e-10,btol=1e-10,condlim=1e5,maxIter=10)

		@test norm(xgt-x1[1])/norm(xgt) < 1e-6
		@test norm(xgt-x2[1])/norm(xgt) < 1e-6
		@test norm(x1[1]-x2[1])/norm(x1[1]) < 1e-12
		conA = norm(vec(A))*norm(vec(pinv(A)))
		@test norm(conA-x1[3][end,end])/conA < 1e-12
		@test norm(conA-x2[3][end,end])/conA < 1e-12
	end

end


text:
module TestWiringDiagramExpressions

using Test

using Catlab.Theories, Catlab.Graphs, Catlab.WiringDiagrams
using Catlab.Syntax: syntax_module
using Catlab.WiringDiagrams.WiringDiagramExpressions: find_parallel,
  find_series, transitive_reduction!

# Expression -> Diagram
#######################

# Directed wiring diagrams
#-------------------------

A, B, C = Ob(FreeSymmetricMonoidalCategory, :A, :B, :C)
I = munit(FreeSymmetricMonoidalCategory.Ob)
f, g = Hom(:f,A,B), Hom(:g,B,C)

# Functorality of conversion.
fd, gd = to_wiring_diagram(f), to_wiring_diagram(g)
@test to_wiring_diagram(compose(f,g)) == compose(fd,gd)
@test to_wiring_diagram(otimes(f,g)) == otimes(fd,gd)
@test to_wiring_diagram(I) == munit(Ports)

# Undirected wiring diagrams
#---------------------------

const to_uwd = to_undirected_wiring_diagram
const HDOb = HypergraphDiagramOb{Symbol,Symbol}
const HDHom = HypergraphDiagramHom{Symbol,Symbol}

A, B, C = Ob(FreeHypergraphCategory, :A, :B, :C)
I = munit(FreeHypergraphCategory.Ob)
f, g = Hom(:f,A,B), Hom(:g,B,C)

# Generators.
@test to_uwd(A) == HDOb([:A])
@test to_uwd(f) == HDHom(singleton_diagram(HDOb([:A,:B]), name=:f),
                         dom=BitArray([1,0]))

# Functorality of conversion.
fd, gd = to_uwd(f), to_uwd(g)
@test to_uwd(compose(f,g)) == compose(fd,gd)
@test to_uwd(otimes(f,g)) == otimes(fd,gd)
@test to_uwd(I) == munit(HypergraphDiagramOb{Symbol,Symbol})

# Diagram -> Expression
#######################

roundtrip(A::ObExpr) = to_ob_expr(syntax_module(f), to_wiring_diagram(A))
roundtrip(f::HomExpr) = to_hom_expr(syntax_module(f), to_wiring_diagram(f))

# Symmetric monoidal category
#----------------------------

A, B, C, D = Ob(FreeSymmetricMonoidalCategory, :A, :B, :C, :D)
I = munit(FreeSymmetricMonoidalCategory.Ob)
f, g, h, k = Hom(:f,A,B), Hom(:g,B,C), Hom(:h,C,D), Hom(:k,D,C)

# Base case.
@test roundtrip(f) == f

# Series reduction.
@test roundtrip(compose(f,g)) == compose(f,g)
@test roundtrip(compose(f,g,h)) == compose(f,g,h)

# Parallel reduction.
@test roundtrip(otimes(f,g)) == otimes(f,g)
@test roundtrip(otimes(f,g,h)) == otimes(f,g,h)

# Series-parallel reduction.
@test roundtrip(otimes(compose(f,g),compose(h,k))) ==
  # roundtrip(compose(otimes(f,h),otimes(g,k))) ==
  otimes(compose(f,g),compose(h,k))
@test roundtrip(otimes(compose(f,g),h)) == otimes(compose(f,g),h)
@test roundtrip(otimes(f,compose(g,h))) == otimes(f,compose(g,h))

m = Hom(:m, otimes(B,B), otimes(C,C))
expr = compose(otimes(f,f),m,otimes(h,h))
@test roundtrip(expr) == expr

m, n = Hom(:m, A, otimes(A,B)), Hom(:n, otimes(B,C), C)
@test roundtrip(compose(m,otimes(f,g),n)) == compose(m,otimes(f,g),n)

m, n = Hom(:m, A, otimes(B,B)), Hom(:n, otimes(B,B), C)
expr = compose(otimes(m,f),otimes(g,n))
@test roundtrip(expr) == expr

m, n = Hom(:m, A, otimes(B,B)), Hom(:n, otimes(C,C), D)
l = Hom(:l, otimes(B,B), otimes(C,C))
expr = compose(otimes(m,m),otimes(g,l,g),otimes(n,n))
@test roundtrip(expr) == expr

# Transitive reduction (necessarily with series and/or parallel reduction).
@test roundtrip(otimes(f,id(C))) == otimes(f,id(C))
@test roundtrip(otimes(id(A),g)) == otimes(id(A),g)

m = Hom(:m, otimes(B,B), otimes(C,C))
expr = compose(otimes(f,id(B)),m,otimes(id(C),h))
@test roundtrip(expr) == expr

m = Hom(:m, otimes(B,B,B), C)
@test roundtrip(compose(otimes(f,id(otimes(B,B))),m)) ==
  compose(otimes(f,id(B),id(B)),m)

m, n = Hom(:m, A, otimes(A,B)), Hom(:n, otimes(B,B), C)
expr = compose(m,otimes(f,id(B)),n)
@test roundtrip(expr) == expr

m = Hom(:m, otimes(A,A), A)
expr = compose(otimes(m,id(A),id(A)), otimes(m,id(A)), m)
@test roundtrip(expr) == expr

m, n = Hom(:m, A, otimes(B,B)), Hom(:n, otimes(B,B), C)
expr = compose(otimes(m,id(B)), otimes(id(B),n))
@test roundtrip(expr) == expr
b = Hom(:b, B, B)
expr = compose(otimes(m,id(B)), otimes(id(B),b,id(B)), otimes(id(B),n))
@test roundtrip(expr) == expr

# Identities.
@test roundtrip(id(A)) == id(A)
@test roundtrip(otimes(id(A),id(B))) == id(otimes(A,B))

# Monoidal units.
@test roundtrip(Hom(:m,A,I)) == Hom(:m,A,I)
@test roundtrip(Hom(:m,I,B)) == Hom(:m,I,B)
@test roundtrip(Hom(:m,I,I)) == Hom(:m,I,I)
@test roundtrip(Hom(:m,I,B)⋅Hom(:n,B,I)) == Hom(:m,I,B)⋅Hom(:n,B,I)
@test roundtrip(Hom(:m,A,I)⋅Hom(:n,I,B)) == Hom(:m,A,I)⊗Hom(:n,I,B)
@test roundtrip(Hom(:m,I,I)⊗Hom(:n,I,I)) == Hom(:m,I,I)⊗Hom(:n,I,I)

# Monoidal units with transitive reduction.
@test roundtrip(Hom(:m,A,I)⊗id(B)) == Hom(:m,A,I)⊗id(B)
@test roundtrip(id(A)⊗Hom(:m,I,B)) == id(A)⊗Hom(:m,I,B)
@test roundtrip(id(A)⊗Hom(:m,I,I)) == id(A)⊗Hom(:m,I,I)

expr = compose(id(A)⊗Hom(:m,I,A), Hom(:l, A⊗A, B⊗B), Hom(:n,B,I)⊗id(B))
@test roundtrip(expr) == expr

# Braidings.
@test roundtrip(braid(A,B)) == braid(A,B)
@test roundtrip(otimes(id(A),braid(B,C))) == otimes(id(A),braid(B,C))

# Diagonals and codiagonals
#--------------------------

A, B, C = Ob(FreeBiproductCategory, :A, :B, :C)
f, g = Hom(:f,A,B), Hom(:g,B,C)

# Diagonals.
@test roundtrip(mcopy(A)) == mcopy(A)
@test roundtrip(delete(A)) == delete(A)
@test roundtrip(mcopy(A)⋅(mcopy(A)⊗id(A))) == mcopy(A)⋅(mcopy(A)⊗id(A))
@test roundtrip(mcopy(A)⋅(id(A)⊗mcopy(A))) == mcopy(A)⋅(mcopy(A) ⊗ id(A))
@test roundtrip(mcopy(A⊗B)) == (mcopy(A)⊗mcopy(B))⋅(id(A)⊗braid(A,B)⊗id(B))
@test roundtrip(delete(A⊗B)) == delete(A)⊗delete(B)

# Codiagonals.
@test roundtrip(mmerge(A)) == mmerge(A)
@test roundtrip(create(A)) == create(A)
@test roundtrip((mmerge(A)⊗id(A)⋅mmerge(A))) == (mmerge(A)⊗id(A))⋅mmerge(A)
@test roundtrip((id(A)⊗mmerge(A))⋅mmerge(A)) == (mmerge(A)⊗id(A))⋅mmerge(A)

# Compound morphisms.
@test roundtrip(compose(f,mcopy(B))) == compose(f,mcopy(B))
@test roundtrip(compose(mcopy(A),otimes(f,f))) == compose(mcopy(A),otimes(f,f))
@test roundtrip(compose(f,delete(B))) == compose(f,delete(B))
@test roundtrip(braid(A,B)⋅(mcopy(B)⊗mcopy(A))) == braid(A,B)⋅(mcopy(B)⊗mcopy(A))

# Dagger category
#----------------

A, B, C = Ob(FreeDaggerSymmetricMonoidalCategory, :A, :B, :C)
f, g = Hom(:f,A,B), Hom(:g,B,C)

@test roundtrip(dagger(f)) == dagger(f)
@test roundtrip(dagger(compose(f,g))) == compose(dagger(g),dagger(f))

# Compact closed category
#------------------------

A, B, C = Ob(FreeCompactClosedCategory, :A, :B, :C)
f, g = Hom(:f,A,B), Hom(:g,B,C)

# Duals.
@test roundtrip(dual(A)) == dual(A)
@test roundtrip(dual(otimes(A,B))) == otimes(dual(B),dual(A))

# Units and counits.
@test roundtrip(dunit(A)) == dunit(A)
@test roundtrip(dcounit(A)) == dcounit(A)
# TODO: The following expressions are correct, but it would be nice to return
# them in the commented simplified form.
@test roundtrip(dunit(otimes(A,B))) ==
  #dunit(B) ⋅ (id(dual(B)) ⊗ dunit(A) ⊗ id(B))
  (dunit(A) ⊗ dunit(B)) ⋅
    (id(dual(A))⊗braid(A,dual(B))⊗id(B)) ⋅ (braid(dual(A),dual(B))⊗id(A⊗B))
@test roundtrip(dcounit(otimes(A,B))) ==
  #(id(A) ⊗ dcounit(B) ⊗ id(dual(A))) ⋅ dcounit(A)
  (id(A⊗B)⊗braid(dual(B),dual(A))) ⋅ (id(A)⊗braid(B,dual(A))⊗id(dual(B))) ⋅
    (dcounit(A) ⊗ dcounit(B))

# Adjoint mates.
@test roundtrip(mate(f)) == mate(f)
@test roundtrip(mate(compose(f,g))) == compose(mate(g),mate(f))

# Graph operations
##################

# Parallel compositions in graph.
graph = Graph(6)
add_edges!(graph, [1,2,3,3,4,5], [2,3,4,5,6,6])
@test find_parallel(graph) == Dict((3 => 6) => [4,5])

# Series compositions in graph.
graph = Graph(13)
add_edges!(graph, [5,8,9,11,12], [6,9,10,12,13])
@test Set(find_series(graph)) == Set([[5,6],[8,9,10],[11,12,13]])

end


text:
module Dummy_AI

using ClientAI

# the DummyAI class
type DummyAI <: AbstractAI
  config::Config
  bots::Vector{AbstractBot}
end

# This function is called when the game starts and gets information about own and enemy bots.
function on_start(ai::DummyAI, own_bots::Vector, enemy_bots::Vector)
end

# This function is called at the beginning of every round, and gets information about own bots and
# the events that happend in the last round. It is called before events are processed by the
# on_event handlers.
function init_round(ai::DummyAI, bots::Vector{AbstractBot}, events::Vector{AbstractEvent}, round_id::Integer)
  ai.bots = bots
end

# this function is called after event handling and has to return an array of actions that the bot wants to
# perform.
function decide(ai::DummyAI)
  bots = filter_valid(ai.bots)
  actions = AbstractAction[]
  for b in bots
    push!(actions, MoveAction(botid(b), rand(collect(move_area(b, ai.config)))))
  end
  return actions
end

# ignore all events
ClientAI.on_event(ai::DummyAI, event::AbstractEvent) = false

# this function is called from main to create the AI
function create(team_id, config::Config)
	info("create ai $team_id")
	return DummyAI(config, AbstractBot[])
end

end

max_stars_repo_name: Condimentum/spacetyckiting

max_issues_repo_name: Condimentum/spacetyckiting

max_forks_repo_name: Condimentum/spacetyckiting


text:
using Statistics
using Printf
using Flux
using Plots

# Flux uses Float32 by default so let's just use that!
const FT = Float32

function generate_training_data(f; N, domain)
    x = range(domain[1], domain[2], length=N)
    return @. FT(x), FT(f(x))
end

function generate_neural_network(; layers, nodes=1, activation=identity)
    return Chain([Dense(l == 1 ? 1 : nodes, l == layers ? 1 : nodes, activation) for l in 1:layers]...)
end

function learn_scalar_function!(NN, f; N, domain, epochs, optimizers=[ADAM()])
    xs, ys = generate_training_data(f; N, domain)
    training_data = zip(xs, ys) |> collect

    loss(x, y) = Flux.mse(NN([x]), y)
    function callback()
        μ_loss = mean(loss(d...) for d in training_data)
        @info @sprintf("mean squared error loss = %.16e", μ_loss)
        return μ_loss
    end

    for opt in optimizers, e in 1:epochs
        @info "Training with $(typeof(opt))(η=$(opt.eta)) [epoch $e/$epochs]..."
        Flux.train!(loss, Flux.params(NN), training_data, opt, cb=callback)
    end

    return NN
end

function test_learned_mapping(f, NN; N, domain)
    xs, ys_truth = generate_training_data(f; N, domain)
    ys_NN = [NN([x])[1] for x in xs]

    kwargs = (linewidth=3, linealpha=0.8, xlims=extrema(xs),
              grid=false, legend=:bottomright, framestyle=:box,
              foreground_color_legend=nothing, background_color_legend=nothing)

    fit_plot = plot(xs, ys_truth, label="truth", title="Fit", xlabel="x", ylabel="y"; kwargs...)
    plot!(fit_plot, xs, ys_NN, label="neural network"; kwargs...)

    loss = @. abs(ys_NN - ys_truth) / abs(ys_truth)
    ylims = (10^floor(log10(minimum(loss))), 10^ceil(log10(maximum(loss))))
    loss_plot = plot(xs, loss, label="", title="extrapolation accuracy", xlabel="x", ylabel="|y′ - y| / y", yaxis=:log, ylims=ylims; kwargs...)

    return plot(fit_plot, loss_plot, layout=(2, 1))
end

linear(x) = x
NN = generate_neural_network(layers=1)
learn_scalar_function!(NN, linear, N=10, domain=(-1, 1), epochs=20, optimizers=[Descent()])
test_learned_mapping(linear, NN, N=100, domain=(-100, 100))

quadratic(x) = x^2
NN = generate_neural_network(layers=20, nodes=20)
learn_scalar_function!(NN, quadratic, N=10, domain=(0, 1), epochs=10, optimizers=[ADAM()])
test_learned_mapping(quadratic, NN, N=100, domain=(0, 2))


text:
# This file was generated, do not modify it. # hide
using RDatasets, DataFrames

boston = dataset("MASS", "Boston");


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -8.186501529651974 <= q <= 9.622483160840265)

                     add_NL_constraint(m, :(1/(1 + exp(-(-0.9342689023247241 + -0.06643065035242302*1/(1 + exp(-(-0.6967397920955185 + 0.4612305646899797*1/(1 + exp(-(-0.7651411836742423 + 0.736339264438953*$(x[1]) + 0.06323249947984699*$(x[2]) + 0.6326865446578034*$(x[3]) + 0.8935006364653741*$(x[4])))) + 0.5917923558177942*1/(1 + exp(-(0.8959549878207418 + 0.6412415470262633*$(x[1]) + 0.6383904188579423*$(x[2]) + 0.9706157764666923*$(x[3]) + 0.8108132764672078*$(x[4]))))))) + 0.13022424784864484*1/(1 + exp(-(-0.8092681495291991 + -0.6577545463201151*1/(1 + exp(-(-0.7651411836742423 + 0.736339264438953*$(x[1]) + 0.06323249947984699*$(x[2]) + 0.6326865446578034*$(x[3]) + 0.8935006364653741*$(x[4])))) + 0.816846389391447*1/(1 + exp(-(0.8959549878207418 + 0.6412415470262633*$(x[1]) + 0.6383904188579423*$(x[2]) + 0.9706157764666923*$(x[3]) + 0.8108132764672078*$(x[4])))))))))) + 1/(1 + exp(-(0.9502905692637253 + 0.2731932353707758*1/(1 + exp(-(-0.6967397920955185 + 0.4612305646899797*1/(1 + exp(-(-0.7651411836742423 + 0.736339264438953*$(x[1]) + 0.06323249947984699*$(x[2]) + 0.6326865446578034*$(x[3]) + 0.8935006364653741*$(x[4])))) + 0.5917923558177942*1/(1 + exp(-(0.8959549878207418 + 0.6412415470262633*$(x[1]) + 0.6383904188579423*$(x[2]) + 0.9706157764666923*$(x[3]) + 0.8108132764672078*$(x[4]))))))) + 0.7525923835832526*1/(1 + exp(-(-0.8092681495291991 + -0.6577545463201151*1/(1 + exp(-(-0.7651411836742423 + 0.736339264438953*$(x[1]) + 0.06323249947984699*$(x[2]) + 0.6326865446578034*$(x[3]) + 0.8935006364653741*$(x[4])))) + 0.816846389391447*1/(1 + exp(-(0.8959549878207418 + 0.6412415470262633*$(x[1]) + 0.6383904188579423*$(x[2]) + 0.9706157764666923*$(x[3]) + 0.8108132764672078*$(x[4])))))))))) - $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:
"""
$(TYPEDEF)

A lorenz attractor with with two inputs and parameters `σ`, `ρ` and `β`.
"""
struct Lorenz3D{ΔT} <: ControlSystem{ΔT, 3, 2}
    σ::Float64;
    ρ::Float64;
    β::Float64;
end

function dx(cs::Lorenz3D, x::SVector{3}, u::SVector{2}, t::AbstractFloat)
    dx1 = cs.σ*(x[2] - x[1])
    dx2 = x[1]*(cs.ρ - x[3]) - x[2] + u[2]
    dx3 = x[1]*x[2] - cs.β*x[3] + u[1]

    return @SVector[dx1, dx2, dx3]
end


text:
@testset "Locally free class group of group algebras" begin
  G = small_group(8, 4)
  A = AlgGrp(FlintQQ, G)
  O = Order(A, basis(A))
  C = locally_free_class_group(O)
  @test C.snf == fmpz[ 2 ]

  A = AlgAss(A)[1]
  O = Order(A, basis(A))
  C = locally_free_class_group(O)
  @test C.snf == fmpz[ 2 ]

  G = small_group(10, 2)
  A = AlgGrp(FlintQQ, G)
  O = Order(A, basis(A))
  C = locally_free_class_group(O)
  @test C.snf == fmpz[]

  A = AlgAss(A)[1]
  O = Order(A, basis(A))
  C = locally_free_class_group(O)
  @test C.snf == fmpz[]

  G = small_group(12, 3)
  A = AlgAss(AlgGrp(FlintQQ, G))[1]
  O = Order(A, basis(A))
  C = locally_free_class_group(O)
  @test C.snf == fmpz[]
end

@testset "Locally free class group of matrix algebras" begin
  Qx, x = FlintQQ["x"]
  f = x^2 + 47
  K, a = number_field(f, "a")
  A = AlgAss(MatrixAlgebra(K, 2))
  A, _, _ = Hecke.restrict_scalars(A, FlintQQ)
  O = MaximalOrder(A)
  C = Hecke.locally_free_class_group(O) # == class_group(K)
  @test C.snf == fmpz[ 5 ]

  f = x^2 + 26
  K, a = number_field(f, "a")
  A = AlgAss(MatrixAlgebra(K, 2))
  A, _, _ = Hecke.restrict_scalars(A, FlintQQ)
  O = MaximalOrder(A)
  C = Hecke.locally_free_class_group(O) # == class_group(K)
  @test C.snf == fmpz[ 6 ]
end

@testset "Discrete logarithm of locally free class group" begin
  Qx, x = FlintQQ["x"]
  f = x^12 - x^11 + x^9 - x^8 + x^6 - x^4 + x^3 - x + 1
  K, a = number_field(f, "a") # Gal(K/Q) = C2 x C6 (aka 12T2 aka small_group(12, 5))
  OK = maximal_order(K)
  A, KtoA, AtoK = Hecke._as_group_algebra(K, a)
  basisOK = [ KtoA(b.elem_in_nf) for b in basis(OK) ]
  d = lcm([ denominator(b) for b in basisOK ])
  ZG = Order(A, basis(A))
  I = Hecke.ideal_from_lattice_gens(A, ZG, [ d*b for b in basisOK ])
  S, mS = locally_free_class_group_with_disc_log(ZG)
  @test S.snf == fmpz[ 2, 2 ]
  @test iszero(mS(I))

  # Check whether one can also call it with AlgAss
  B, BtoA = AlgAss(A)
  basisOK2 = [ BtoA\b for b in basisOK ]
  d2 = lcm([ denominator(b) for b in basisOK2 ])
  ZG = Order(B, basis(B))
  I = Hecke.ideal_from_lattice_gens(B, ZG, [ d*b for b in basisOK2 ])
  S, mS = locally_free_class_group_with_disc_log(ZG, check = false)
  @test S.snf == fmpz[ 2, 2 ]
  @test iszero(mS(I))

  f = x^8 - 3*x^7 + 22*x^6 - 60*x^5 + 201*x^4 - 450*x^3 + 1528*x^2 - 3069*x + 4561
  K, a = number_field(f, "a") # Gal(K/Q) = Q8 (aka 8T5 aka small_group(8, 4))
  OK = maximal_order(K)
  A, KtoA, AtoK = Hecke._as_group_algebra(K, a)
  basisOK = [ KtoA(b.elem_in_nf) for b in basis(OK) ]
  d = lcm([ denominator(b) for b in basisOK ])
  ZG = Order(A, basis(A))
  I = Hecke.ideal_from_lattice_gens(A, ZG, [ d*b for b in basisOK ])
  S, mS = locally_free_class_group_with_disc_log(ZG)
  @test S.snf == fmpz[ 2 ]
  @test mS(I) == S[1]
  @test iszero(mS(I^2))
end


text:
baremodule MouseButton
    const LEFT = 0
    const RIGHT = 1
    const MIDDLE = 2
    const XBUTTON1 = 3
    const XBUTTON2 = 4
end

function is_mouse_pressed(button::Int)
    return Bool(ccall((:sfMouse_isButtonPressed, libcsfml_window), Int32, (Int32,), button))
end

function get_mousepos(window::RenderWindow)
    return ccall((:sfMouse_getPosition, libcsfml_window), Vector2i, (Ptr{Void},), window.ptr)
end

function set_mousepos(pos::Vector2, window::RenderWindow)
    pos = to_vec2i(pos)
    ccall((:sfMouse_setPosition, libcsfml_window), Void, (Vector2i, Ptr{Void},), pos, window.ptr)
end


text:
export UniformEnergyDistribution

"Initial uniform energy distribution."
struct UniformEnergyDistribution{N} <: InitialEnergyDistribution{N}
    "Specific energy."
    u::Float

    UniformEnergyDistribution{N}(u::Number) where {N} = @check(u > 0) && new(u)
end

Base.show(io::IO, ::Type{UniformEnergyDistribution{N}}) where {N} =
    print(io, "UniformEnergyDistribution")

function (distribution::UniformEnergyDistribution{N})(
    positions::AbstractVector{SVectorF{N}},
    mass_densities::AbstractVector{Float},
    eos::EquationOfState,
) where {N}
    @checkeq(length(positions), length(mass_densities))
    fill(distribution.u, length(positions))
end

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
using Pkg
Pkg.activate(normpath(joinpath(@__DIR__, "..")))

using ATools
using NReco

using ArgParse
using Configurations
using DataFrames
using Glob
using HDF5
using Logging
using Printf

logger = SimpleLogger(stdout, Logging.Warn)
old_logger = global_logger(logger)

function makenema(args)

	lorf    = args["loralgo"]
	detconf = args["detconf"]
	dr      = args["dir"]
	outd    = args["odir"]
	outf    = args["ofile"]
	file_i  = args["filei"]
	file_l  = args["filel"]

	lor_algo = NReco.lor_kmeans
	if lorf == "lor_qmax"
		lor_algo = NReco.lor_maxq
	end

	if isdir(outd) == false
		mkdir(outd)
	end
	output = string(outd,"/", outf)

	## Patch here so we can read actually select in a directory which
	## is already partially processed. Needs to be improved in general.
	files  = sort(glob("*.h5", dr), by=x->parse(Int64, split(x, "-")[end][1:end-3]))

	if detconf != "default"
		dconf = from_toml(NReco.DetConf, detconf)
	else
		dconf = NReco.DetConf()
	end

	println("makenema configuration")
	println("detector configuration", dconf)
	println(" lor_algo  = $lorf")

	println("number of files in data dir = $(length(files))")
	println("reading = $(file_l - file_i + 1) files")
	println("output file  = $output")

	h5open(output, "w") do h5out
		pars_grp    = create_group(h5out, "selected_events")
		pars_dtype  = ATools.generate_hdf5_datatype(ATools.EventParameters)
		conf_grp    = create_group(h5out, "configuration")
		conf_dtype  = ATools.generate_hdf5_datatype(ATools.SimConfiguration)

		nevt, evts  = NReco.nemareco(files, dconf, file_i, file_l, lor_algo)

		rmin, rmax  = extrema(getfield.(evts, :r1))
		# Maybe combine the DetConf and SimConfiguration types?
		config_pars = ATools.SimConfiguration(SiPMPDE=dconf.pde,
						ElecSTD=dconf.sigma_tof, SiPMThr=dconf.ecut,
						SumQmin=dconf.qmin, SumQmax=dconf.qmax,
						Ntof=dconf.ntof, NEvent=nevt, Rmin=rmin,
						Rmax=rmax)#TODO: Calibration option implementation

		## Save to file.
		pars_dspace = dataspace(evts)
		pars_dset   = create_dataset(pars_grp, "EventParameters", pars_dtype, pars_dspace)
		write_dataset(pars_dset, pars_dtype, evts)
		conf_dspace = dataspace([config_pars])
		conf_dset   = create_dataset(conf_grp, "RunConfiguration", conf_dtype, conf_dspace)
		write_dataset(conf_dset, conf_dtype, [config_pars])#Must be a better way
	end
end

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--dir", "-d"
            help = "directory with nema simulations"
            arg_type = String
		"--odir", "-o"
            help = "output directory"
            arg_type = String
		"--ofile", "-x"
            help = "output file"
            arg_type = String
            default = "evtdf.h5"
		"--filei", "-i"
	        help = "number of initial file in glob list"
	        default  = 1
			arg_type = Int
		"--filel", "-l"
		    help = "number of last file in glob list"
		    default  = 1
			arg_type = Int
		"--loralgo", "-g"
			help = "algorithm to use for LOR reconstruction "
			arg_type = String
			default = "lor_kmeans"
		"--detconf", "-c"
			help = "Detector configuration"
			arg_type = String
			default = "default"
    end

    return parse_args(s)
end

function main()
	parsed_args = parse_commandline()
	println("Running makenema with arguments", parsed_args)
	makenema(parsed_args)
end

@time main()

max_stars_repo_name: andLaing/NReco

max_issues_repo_name: andLaing/NReco

max_forks_repo_name: andLaing/NReco


text:
"""
Functions to calculate atomic collisional rates.
"""

# Dimension for CE and CI collisional tables, in SI units m^3 K^-1/2 s-1
@derived_dimension CEI_dimension Unitful.𝐋^3 / (Unitful.𝚯^(1/2)  * Unitful.𝐓)

const Ω_c0 = Ry / sqrt(m_e) * π * a_0^2 * sqrt(8 / (π * k_B))

"""
    coll_CE(
        rate_interpolant::Interpolations.AbstractInterpolation{<:CEI_dimension, 1},
        g_ratio::Real,
        electron_density::NumberDensity,
        temperature::Unitful.Temperature
    )

Calculate collisional de-excitation by electrons of a bound-bound transition, using
the CE expression from RH / MULTI.

# Arguments
- `rate_interpolant`: a generic interpolant that takes as argument a value of
temperature, and returns in units similar to m^3 K^-1/2 s^-1
- `g_ratio`: ratio g_l / g_u, statistical weights of lower and upper level
- `electron_density`: electron density
- `temperature`: gas temperature

# Returns
- `coll_deexc`: collisional de-excitations per second (from upper level to lower level)
"""
function coll_CE(
    rate_interpolant::Interpolations.AbstractInterpolation{<:CEI_dimension, 1},
    g_ratio::Real,
    electron_density::NumberDensity,
    temperature::Unitful.Temperature
)
    return (rate_interpolant(temperature) * g_ratio *
            electron_density * sqrt(temperature)) |> u"s^-1"
end

"""
    coll_CI(
        rate_interpolant::Interpolations.AbstractInterpolation{<:CEI_dimension, 1},
        dE::Unitful.Energy,
        electron_density::NumberDensity,
        temperature::Unitful.Temperature
    )

Calculate collisional ionisation by electrons of a given level, using
the CI expression from RH / MULTI.

# Arguments
- `rate_interpolant`: a generic interpolant that takes as argument a value of
temperature, and returns in units similar to m^3 K^-1/2 s^-1
- `dE`: energy difference between continuum and level, dE = E_cont - E_level.
- `electron_density`: electron density
- `temperature`: gas temperature

# Returns
- `coll_ion`: collisional ionisations per second (from level to continuum)
"""
function coll_CI(
    rate_interpolant::Interpolations.AbstractInterpolation{<:CEI_dimension, 1},
    dE::Unitful.Energy,
    electron_density::NumberDensity,
    temperature::Unitful.Temperature
)
    @assert dE > 0u"J" "dE should be positive"
    return (rate_interpolant(temperature) * exp(-dE / (k_B * temperature)) *
            electron_density * sqrt(temperature)) |> u"s^-1"
end

"""
    coll_Ω(
        rate_interpolant::Interpolations.AbstractInterpolation{<:CEI_dimension, 1},
        g_u::Integer,
        electron_density::NumberDensity,
        temperature::Unitful.Temperature
    )

Calculate collisional de-excitation by electrons of a bound-bound transition, using
the OMEGA expression from RH (OHMEGA from MULTI), from a tabulated dimensionless Ω
(collision strength).

# Arguments
- `rate_interpolant`: a generic interpolant that takes as argument a value of
temperature, and returns a dimensionless Ω
- `g_u`: statistical weight of upper level
- `electron_density`: electron density
- `temperature`: gas temperature

# Returns
- `coll_deexc`: collisional de-excitations per second (from upper level to lower level)
"""
function coll_Ω(
    rate_interpolant::Interpolations.AbstractInterpolation{<:Real, 1},
    g_u::Integer,
    electron_density::NumberDensity,
    temperature::Unitful.Temperature
)
    return (rate_interpolant(temperature) * Ω_c0 * electron_density /
            (g_u * sqrt(temperature))) |> u"s^-1"
end

#=----------------------------------------------------------------------------
Collisional rates for hydrogen from Johnson (1972), based on make_h.c from RH
----------------------------------------------------------------------------=#
const johnson_c0 = sqrt(8 * k_B / (π * m_e)) * 2 * π * a_0^2
const johnson_c1 = 32 / (3 * sqrt(3) * π)

"""
    coll_exc_hydrogen_johnson(
        n_l::Integer,
        n_u::Integer,
        electron_density::NumberDensity,
        temperature::Unitful.Temperature,
    )

Calculates the rates for collisional excitation for a hydrogen line, resulting
from collisions with electrons, using the recipe of
[Johnson (1972)](https://ui.adsabs.harvard.edu/abs/1972ApJ...174..227J/abstract)

# Arguments
- `n_l`: lower level (principal quantum number) of transition
- `n_u`: upper level (principal quantum number) of transition
- `electron_density`: electron density
- `temperature`: gas temperature

# Returns
- `coll_exc`: collisional excitations per second from n_l -> n_u
"""
function coll_exc_hydrogen_johnson(
    n_l::Integer,
    n_u::Integer,
    electron_density::NumberDensity,
    temperature::Unitful.Temperature,
)
    @assert n_l > 0 "n_l must be > 0"
    @assert n_u > n_l "n_u must be > n_l"
    rn = _rn(n_l)
    bn = _bn(n_l)
    x = 1 - (n_l / n_u)^2
    rnn = rn * x
    fnn = johnson_c1 * n_l / (n_u * x)^3 * (g0(n_l) + (g1(n_l) + g2(n_l) / x) / x)
    Ann = 2 * n_l^2 / x * fnn
    Bnn = 4 * n_l * (n_l / n_u)^3 * (1 + 4 / (3 * x) + bn / x^2) / x^2
    y = x * Ryh  / (n_l^2 * k_B * temperature)
    z = rnn + y
    coll_exc = (johnson_c0 * (n_l * y)^2 / x * sqrt(temperature) * electron_density *
                (Ann * ((1 / y + 0.5) * expint(1, y) - (1 / z + 0.5) * expint(1, z)) +
                 (Bnn - Ann * log(2 * n_l^2 / x)) * (expint(2, y) / y - expint(2, z) / z)))
    return coll_exc |> u"s^-1"
end

"""
    coll_ion_hydrogen_johnson(
        n::Integer,
        electron_density::NumberDensity,
        temperature::Unitful.Temperature,
    )

Calculates the rates for collisional ionisation for a hydrogen level, resulting
from collisions with electrons, using the recipe of
[Johnson (1972)](https://ui.adsabs.harvard.edu/abs/1972ApJ...174..227J/abstract).

# Arguments
- `n`: principal quantum number of hydrogen level
- `electron_density`: electron density
- `temperature`: gas temperature

# Returns
- `coll_ion`: collisional ionisations per second from level n
"""
function coll_ion_hydrogen_johnson(
    n::Integer,
    electron_density::NumberDensity,
    temperature::Unitful.Temperature,
)
    @assert n > 0 "n must be > 0"
    An = johnson_c1 * n * (g0(n) / 3 + g1(n) / 4 + g2(n) / 5)
    Bn = 2 * n^2 / 3 * (5 + _bn(n))
    yn = Ryh / (n^2 * k_B * temperature)
    zn = _rn(n) + yn
    coll_ion = (johnson_c0 * (n * yn)^2 * sqrt(temperature) * electron_density *
                (An * (expint(1, yn) / yn - expint(1, zn) / zn) +
                 (Bn - An * log(2 * n^2)) * (ξ(yn) - ξ(zn))))
    return coll_ion |> u"s^-1"
end

"""
    CI_RH_hydrogen(n::Integer, temperature::Unitful.Temperature)

Calculate CI coefficients for collisional ionisation for a hydrogen level,
in the SI units of RH, using the recipe of
[Johnson (1972)](https://ui.adsabs.harvard.edu/abs/1972ApJ...174..227J/abstract).

# Arguments
- `n`: principal quantum number of hydrogen level
- `temperature`: gas temperature
"""
function CI_RH_hydrogen(n::Integer, temperature::Unitful.Temperature)
    tmp = coll_ion_hydrogen_johnson(n, 1.0u"m^-3", temperature)
    yn = Ryh / (n^2 * k_B * temperature)
    return tmp * exp(yn) / sqrt(temperature) / 1.0u"m^-3"
end

"""
    CE_RH_hydrogen(n_l::Integer, n_u::Integer, temperature::Unitful.Temperature)

Calculate CE coefficients for collisional deexcitation for a hydrogen transition,
in the SI units of RH, using the recipe of
[Johnson (1972)](https://ui.adsabs.harvard.edu/abs/1972ApJ...174..227J/abstract).

# Arguments
- `n_l`: lower level (principal quantum number) of transition
- `n_u`: upper level (principal quantum number) of transition
- `temperature`: gas temperature
"""
function CE_RH_hydrogen(n_l::Integer, n_u::Integer, temperature::Unitful.Temperature)
    tmp = coll_exc_hydrogen_johnson(n_l, n_u, 1.0u"m^-3", temperature)
    x = 1 - (n_l / n_u)^2
    y = x * Ryh  / (n_l^2 * k_B * temperature)
    return tmp * exp(y) / sqrt(temperature) / 1.0u"m^-3"
end

#=----------------------------------------------------------------------------
                        Utility functions from Johnson (1972)
----------------------------------------------------------------------------=#
# bn and rn expressions, from eqs (24) and (32) of Johnson (1972)
_bn(n::Integer) = (n == 1) ? -0.603 : (4.0+ (-18.63 + (36.24 - 28.09 / n) / n) / n) / n
_rn(n::Integer) = (n == 1) ? 0.45 : 1.94 * n^-1.57

# ξ(t), eq (42)
ξ(t::AbstractFloat) = expint(0, t) - 2 * expint(1, t) + expint(2, t)

# Gaunt factor coefficients, Table 1 of Johnson (1972)
function g0(n::Integer)
    if n == 1
        return 1.1330f0
    elseif n == 2
        return 1.0785f0
    else
        return 0.9935f0 + (0.2328f0 - 0.1296f0 / n) / n
    end
end

function g1(n::Integer)
    if n == 1
        return -0.4059f0
    elseif n == 2
        return -0.2319f0
    else
        return -(0.6282 - (0.5598 - 0.5299 / n) / n) / n
    end
end

function g2(n::Integer)
    if n == 1
        return 0.07014f0
    elseif n == 2
        return 0.02947f0
    else
        return (0.3887f0 - (1.181f0 - 1.4700f0 / n) / n) / n^2
    end
end

#=----------------------------------------------------------------------------
        Collisional rates for hydrogen from Przybilla & Butler (2004)
----------------------------------------------------------------------------=#

const PB04_temp = [0.25, 0.5, 0.75, 1, 1.5, 2, 2.5, 3,
                    4, 5, 6, 8, 10, 15, 20, 25] * 10000u"K"
const PB04_Ω =
    [[6.40f-1 6.98f-1 7.57f-1 8.09f-1 8.97f-1 9.78f-1 1.06f+0 1.15f+0   #= 1->2
=#    1.32f+0 1.51f+0 1.68f+0 2.02f+0 2.33f+0 2.97f+0 3.50f+0 3.95f+0];
     [2.20f-1 2.40f-1 2.50f-1 2.61f-1 2.88f-1 3.22f-1 3.59f-1 3.96f-1   #= 1->3
=#    4.64f-1 5.26f-1 5.79f-1 6.70f-1 7.43f-1 8.80f-1 9.79f-1 1.06f+0];
     [9.93f-2 1.02f-1 1.10f-1 1.22f-1 1.51f-1 1.80f-1 2.06f-1 2.28f-1   #= 1->4
=#    2.66f-1 2.95f-1 3.18f-1 3.55f-1 3.83f-1 4.30f-1 4.63f-1 4.88f-1];
     [4.92f-2 5.84f-2 7.17f-2 8.58f-2 1.12f-1 1.33f-1 1.50f-1 1.64f-1   #= 1->5
=#    1.85f-1 2.01f-1 2.12f-1 2.29f-1 2.39f-1 2.59f-1 2.71f-1 2.81f-1];
     [2.97f-2 4.66f-2 6.28f-2 7.68f-2 9.82f-2 1.14f-1 1.25f-1 1.33f-1   #= 1->6
=#    1.45f-1 1.53f-1 1.58f-1 1.65f-1 1.70f-1 1.77f-1 1.82f-1 1.85f-1];
     [5.03f-2 6.72f-2 7.86f-2 8.74f-2 1.00f-1 1.10f-1 1.16f-1 1.21f-1   #= 1->7
=#    1.27f-1 1.31f-1 1.34f-1 1.36f-1 1.37f-1 1.39f-1 1.39f-1 1.40f-1];
     [2.35f+1 2.78f+1 3.09f+1 3.38f+1 4.01f+1 4.71f+1 5.45f+1 6.20f+1   #= 2->3
=#    7.71f+1 9.14f+1 1.05f+2 1.29f+2 1.51f+2 1.93f+2 2.26f+2 2.52f+2];
     [1.07f+1 1.15f+1 1.23f+1 1.34f+1 1.62f+1 1.90f+1 2.18f+1 2.44f+1   #= 2->4
=#    2.89f+1 3.27f+1 3.60f+1 4.14f+1 4.56f+1 5.31f+1 5.83f+1 6.23f+1];
     [5.22f+0 5.90f+0 6.96f+0 8.15f+0 1.04f+1 1.23f+1 1.39f+1 1.52f+1   #= 2->5
=#    1.74f+1 1.90f+1 2.03f+1 2.23f+1 2.37f+1 2.61f+1 2.78f+1 2.89f+1];
     [2.91f+0 4.53f+0 6.06f+0 7.32f+0 9.17f+0 1.05f+1 1.14f+1 1.21f+1   #= 2->6
=#    1.31f+1 1.38f+1 1.44f+1 1.51f+1 1.56f+1 1.63f+1 1.68f+1 1.71f+1];
     [5.25f+0 7.26f+0 8.47f+0 9.27f+0 1.03f+1 1.08f+1 1.12f+1 1.14f+1   #= 2->7
=#    1.17f+1 1.18f+1 1.19f+1 1.19f+1 1.20f+1 1.19f+1 1.19f+1 1.19f+1];
     [1.50f+2 1.90f+2 2.28f+2 2.70f+2 3.64f+2 4.66f+2 5.70f+2 6.72f+2   #= 3->4
=#    8.66f+2 1.04f+3 1.19f+3 1.46f+3 1.67f+3 2.08f+3 2.39f+3 2.62f+3];
     [7.89f+1 9.01f+1 1.07f+2 1.26f+2 1.66f+2 2.03f+2 2.37f+2 2.68f+2   #= 3->5
=#    3.19f+2 3.62f+2 3.98f+2 4.53f+2 4.95f+2 5.68f+2 6.16f+2 6.51f+2];
     [4.13f+1 6.11f+1 8.21f+1 1.01f+2 1.31f+2 1.54f+2 1.72f+2 1.86f+2   #= 3->6
=#    2.08f+2 2.24f+2 2.36f+2 2.53f+2 2.65f+2 2.83f+2 2.94f+2 3.02f+2];
     [7.60f+1 1.07f+2 1.25f+2 1.37f+2 1.52f+2 1.61f+2 1.68f+2 1.72f+2   #= 3->7
=#    1.78f+2 1.81f+2 1.83f+2 1.85f+2 1.86f+2 1.87f+2 1.86f+2 1.87f+2];
     [5.90f+2 8.17f+2 1.07f+3 1.35f+3 1.93f+3 2.47f+3 2.96f+3 3.40f+3   #= 4->5
=#    4.14f+3 4.75f+3 5.25f+3 6.08f+3 6.76f+3 8.08f+3 9.13f+3 1.00f+4];
     [2.94f+2 4.21f+2 5.78f+2 7.36f+2 1.02f+3 1.26f+3 1.46f+3 1.64f+3   #= 4->6
=#    1.92f+3 2.15f+3 2.33f+3 2.61f+3 2.81f+3 3.15f+3 3.36f+3 3.51f+3];
     [4.79f+2 7.06f+2 8.56f+2 9.66f+2 1.11f+3 1.21f+3 1.29f+3 1.34f+3   #= 4->7
=#    1.41f+3 1.46f+3 1.50f+3 1.55f+3 1.57f+3 1.61f+3 1.62f+3 1.63f+3];
     [1.93f+3 2.91f+3 4.00f+3 5.04f+3 6.81f+3 8.20f+3 9.29f+3 1.02f+4   #= 5->6
=#    1.15f+4 1.26f+4 1.34f+4 1.49f+4 1.63f+4 1.97f+4 2.27f+4 2.54f+4];
     [1.95f+3 3.24f+3 4.20f+3 4.95f+3 6.02f+3 6.76f+3 7.29f+3 7.70f+3   #= 5->7
=#    8.26f+3 8.63f+3 8.88f+3 9.21f+3 9.43f+3 9.78f+3 1.00f+4 1.02f+4];
     [6.81f+1 1.17f+4 1.50f+4 1.73f+4 2.03f+4 2.21f+4 2.33f+4 2.41f+4   #= 6->7
=#    2.52f+4 2.60f+4 2.69f+4 2.90f+4 3.17f+4 3.94f+4 4.73f+4 5.50f+4]]
# Index of transition data in matrix form: PB04_index[n_l, n_u] = i,
# where i is row index in PB04_Ω
const PB04_index = [0  1  2  3  4  5  6
                    0  0  7  8  9 10 11
                    0  0  0 12 13 14 15
                    0  0  0  0 16 17 18
                    0  0  0  0  0 19 20
                    0  0  0  0  0  0 21]
const PB04_interp = [LinearInterpolation(PB04_temp, PB04_Ω[i, :],
                                         extrapolation_bc=Line()) for i in 1:16]

"""
    coll_deexc_hydrogen_PB04(
        n_l::Integer,
        n_u::Integer,
        g_u::Integer,
        electron_density::NumberDensity,
        temperature::Unitful.Temperature,
    )

Calculate collisional de-excitation by electrons of a bound-bound hydrogen transition,
using the Ω data from Table 3 of
[Przybilla & Butler (2004)](https://ui.adsabs.harvard.edu/abs/2004ApJ...610L..61P/abstract).
Data available only up to n=7.

# Arguments
- `n_l`: lower level (principal quantum number) of transition
- `n_u`: upper level (principal quantum number) of transition
- `g_u`: statistical weight of upper level
- `electron_density`: electron density
- `temperature`: gas temperature

# Returns
- `coll_deexc`: collisional de-excitations per second (from upper level to lower level)
"""
function coll_deexc_hydrogen_PB04(
    n_l::Integer,
    n_u::Integer,
    g_u::Integer,
    electron_density::NumberDensity,
    temperature::Unitful.Temperature,
)
    @assert 7 > n_l > 0 "Must have 7 > n_l > 0"
    @assert 8 > n_u > n_l "Must have 8 > n_u > n_l"
    return coll_Ω(PB04_interp[PB04_index[n_l, n_u]], g_u, electron_density, temperature)
end


text:
function FirstOrderConverter(arg0::SecondOrderODE)
    return FirstOrderConverter((SecondOrderODE,), arg0)
end

function compute_derivatives(obj::FirstOrderConverter, arg0::jdouble, arg1::Vector{jdouble})
    return jcall(obj, "computeDerivatives", Vector{jdouble}, (jdouble, Vector{jdouble}), arg0, arg1)
end

function get_dimension(obj::FirstOrderConverter)
    return jcall(obj, "getDimension", jint, ())
end

function init(obj::OrdinaryDifferentialEquation, arg0::jdouble, arg1::Vector{jdouble}, arg2::jdouble)
    return jcall(obj, "init", void, (jdouble, Vector{jdouble}, jdouble), arg0, arg1, arg2)
end


text:
"""
	$(SIGNATURES)

Bar graph. Simple wrapper for `barplot` that applies default settings.

# Arguments

- yerror: errorbars around `dataV`.
"""
function bar_graph(groupLabelV, dataV;
    fig = Figure(), pos = (1,1), yerror = nothing, kwargs...)

    @assert size(dataV) == size(groupLabelV);
    args = merge(bar_defaults(), kwargs);
    n = length(groupLabelV);
    ax = make_axis(fig, pos; xticks = (1:n, string.(groupLabelV)), args...);
    barplot!(ax, dataV; args...);
    add_errorbars!(dataV, yerror);
    return fig, ax
end

"""
	$(SIGNATURES)

Plots bar graph into existing axis or FigurePosition.
"""
function bar_graph!(ax :: Axis, dataV; yerror = nothing, kwargs...)
    args = merge(bar_defaults(), kwargs);
    barplot!(ax, dataV;  args...);
    add_errorbars!(dataV, yerror);
end

add_errorbars!(dataV, yerror) =
    errorbars!(1 : length(dataV), dataV, yerror; whiskerwidth = 10);
add_errorbars!(dataV, yerror :: Nothing) = nothing;

bar_defaults() = Dict([
    :strokewidth => 0,
    :color => main_color()
]);

function test_bar_graph()
    n = 10;
    dataV = LinRange(4, 5, n);
    groupLabelV = ["G$j" for j = 1 : n];
    fig, ax = bar_graph(groupLabelV, dataV;
        xlabel = "x label", ylabel = "y label",
        xticks = (1:n, groupLabelV));
    return fig, ax
end

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


text:
# script that wants to numerically test chatteris diaconis (misspelled with 99% prob)
# for the estimates of beta, fitness,ergm (many names..) in the DirBin1 case
using Utilities,AReg,StaticNets,JLD,MLBase,StatsBase#,DynNets using PyCall; pygui(:qt); using PyPlot

Nsample = 100

# load my data
halfPeriod = false
## Load dataj
fold_Path =  "/home/Domenico/Dropbox/Dynamic_Networks/data/emid_data/juliaFiles/"
loadFilePartialName = "Weekly_eMid_Data_from_"
halfPeriod? periodEndStr =  "2012_03_12_to_2015_02_27.jld": periodEndStr =  "2009_06_22_to_2015_02_27.jld"
@load(fold_Path*loadFilePartialName*periodEndStr, AeMidWeekly_T,banksIDs,inactiveBanks, YeMidWeekly_T,weekInd,datesONeMid ,degsIO_T,strIO_T)

matY_T = YeMidWeekly_T

function logLikelihoodGammaVectorObs(μ_vec::Array{<:Real,1},Y_vec::Array{<:Real,1};α::T where T<:Real = 1.0)
    #check that the vector contains only non zero elements
    #prod(Y_vec .>0)?():error()
    ratioVec = Y_vec./μ_vec
    if α!=1
        error("Need to add the term counting the likelihood of zero observations")
        loglike =  α * sum(log.(ratioVec)) - sum(ratioVec) - sum(log.(Y_vec)) - log(gamma(α))
    else
        loglike = - sum(log.(μ_vec) ) - sum(ratioVec)
    end
    return loglike
end

function logLikeGammaDirMatNodePar(parR::Array{<:Real,1},parC::Array{<:Real,1},Y_mat::Array{<:Real,2};α::T where T<:Real = 1.0)

    # create vectors of nonzero observations needed to compute the likelihood of
    # gamma distributions for a matrix with parameters μ_i + μ_j
    A_mat = Y_mat.>0
    indVec =  find(A_mat)
    indR,indC = ind2sub(A_mat,indVec)
    Nnnz = sum(A_mat)
    Y_vec = zeros(Real,Nnnz)
    μ_vec= zeros(Real,Nnnz)
    for n=1:Nnnz
         i,j = indR[n],indC[n]
        #println((i,j))
        Y_vec[n] = Y_mat[i,j]
        μ_vec[n] = parR[i] + parC[j]
    end
    return logLikelihoodGammaVectorObs(μ_vec,Y_vec;α = α )
end

Y_mat = YeMidWeekly_T[:,:,20]
A_mat = Y_mat.>0
N = length(Y_mat[1,:])
logLikeGammaDirMatNodePar(ones(N),ones(N),Y_mat)

using Optim
tol = eps()*10
opt = Optim.Options(  g_tol = tol,
                 x_tol = tol,
                 f_tol = tol,
                 iterations = 5000,
                 show_trace = true,#false,#
                 show_every=1)

algo = Newton(; alphaguess = LineSearches.InitialHagerZhang(),
                 linesearch = LineSearches.BackTracking())

algo = NewtonTrustRegion(; initial_delta = 0.5,
                delta_hat = 1.0,
                eta = 0.1,
                rho_lower = 0.25,
                rho_upper = 0.75)# Newton(;alphaguess = Optim.LineSearches.InitialStatic(),

function objfun(vUnPar::Array{<:Real,1})
    reParR = exp.(vUnPar[1:N])
    reParC = exp.(vUnPar[1+N:end])
    return  - logLikeGammaDirMatNodePar(reParR,reParC,Y_mat)
end

parRC_0 = ones(2N)
ADobjfun = TwiceDifferentiable(objfun, parRC_0; autodiff = :forward);

objfun(parRC_0)

optim_out = optimize(ADobjfun,parRC_0,algo,opt)

Optim.minimizer(optim_out)

#


text:
#! /usr/bin/julia

# Rosetta Code, Permutation test

delta_mean{T<:Real}(a::Array{T,1}, b::Array{T,1}) = mean(a) - mean(b)

function bifurcate{T<:Integer}(a::AbstractVector, sel::Array{T,1})
    x = a[sel]
    asel = trues(length(a))
    asel[sel] = false
    y = a[asel]
    return (x, y)
end

function perm_sig_test{T<:Real}(treat::Array{T,1}, control::Array{T,1})
    base_effect = delta_mean(treat, control)
    pool = [treat, control]
    tlen = length(treat)
    plen = length(pool)
    better = 0
    worse = 0
    for s in combinations(1:plen, tlen)
        (t, c) = bifurcate(pool, s)
        if base_effect < delta_mean(t, c)
            better += 1
        else
            worse += 1
        end
    end
    return (better, worse)
end

treat = [85, 88, 75, 66, 25, 29, 83, 39, 97]
control = [68, 41, 10, 49, 16, 65, 32, 92, 28, 98]

(better, worse) = perm_sig_test(treat, control)

tot = better + worse

println("Permutation test using the following data:")
println("Treated:  ", treat)
println("Control:  ", control)

println()
println("There are ", tot, " different permuted groups of these data.")
print(@sprintf("%8d, %5.2f%% ", better, 100better/tot))
println("showed better than actual results.")
print(@sprintf("%8d, %5.2f%% ", worse, 100worse/tot))
println("showed equalivalent or worse results.")

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
module FFTransforms

using LinearAlgebra
using FFTW
using XFields

const module_dir  = joinpath(@__DIR__, "..") |> normpath

C64 = Complex{Float64}
C32 = Complex{Float32}
F64 = Float64
F32 = Float32

FFTR = Union{F32,F64}
FFTC = Union{C32,C64}
FFTN = Union{FFTR, FFTC}

Plan{T,d,G} = Union{
	FFTW.cFFTWPlan{T,-1,false,d,G},
 	FFTW.rFFTWPlan{T,-1,false,d,G},
	FFTW.cFFTWPlan{T,1,false,d,G},
 	FFTW.rFFTWPlan{T,1,false,d,G}
}

# ft::𝕎{Tf,d,...} <: Transform{Tf,d}
# =========================================
# Adds a lightweight layer between a container for concrete plan
# and the inputs to the planning methods.
# This allows one to easily store 𝕎 as an field in an array type
# wrapper.

struct 𝕎{Tf<:FFTN, d, Tsf<:Number, Tp<:Real} <: Transform{Tf,d}
	sz::NTuple{d,Int}
	region::NTuple{d,Bool}
	scale::Tsf
	period::NTuple{d,Tp}
	function 𝕎{Tf,d}(sz::NTuple{d,Int}, rg::NTuple{d,Bool}, sc::Tsf, pd::NTuple{d,Tp}) where {Tf<:FFTN,d,Tsf,Tp}
		new{Tf,d,Tsf,Tp}(sz,rg,sc,pd)
	end
end

@inline XFields.size_in(w::𝕎) = w.sz

XFields.size_out(w::𝕎{Tf}) where {Tf<:FFTC} = w.sz

function XFields.size_out(w::𝕎{Tf,d})::NTuple{d,Int} where {Tf<:FFTR,d}
    ir = findfirst(w.region)
    return map(w.sz, tuple(1:d...)) do nᵢ, i
        i==ir ? nᵢ÷2+1 : nᵢ
    end
end

@inline XFields.eltype_in(w::𝕎{Tf,d}) where {Tf,d}  = Tf

@inline XFields.eltype_out(w::𝕎{Tf,d}) where {Tf,d} = Complex{real(Tf)}

include("plan_fft.jl")

function XFields.plan(w::𝕎{Tf,d,Tsf}) where {d,Tf<:FFTR,Tsf}
	Ti   = Complex{Tf}
	Tsi  = promote_type(Tf, Tsf)
	plan(Tf,SizeInt{w.sz},RegionBool{w.region},w.scale)::FFTplan{Tf,d,Ti,Tsf,Tsi}
end

function XFields.plan(w::𝕎{Tf,d,Tsf}) where {d,Tf<:FFTC,Tsf}
	Ti  = Tf
	Tsi = promote_type(real(Tf), Tsf)
	plan(Tf,SizeInt{w.sz},RegionBool{w.region},w.scale)::FFTplan{Tf,d,Ti,Tsf,Tsi}
end

export size_in, size_out, eltype_in, eltype_out, plan, FFTplan, AdjointFFTplan

## Extra grid information available
# =====================================
# TODO: incorperate the above into a file for include so this
# fits better in the code

include("grid.jl")

## rand_in, rand_out (incomplete), dot_in, dot_out (incomplete)
# =====================================

include("methods.jl")

## Extra convienent constructors
# =====================================
# The following structs allow lazy construction of an fft plan.
# Mixing 𝕀 and 𝕎 with ⊗ creates another 𝕎 (or 𝕀) untill
# passed to the plan method

include("constructors.jl")

export 𝕀, 𝕎, 𝕌, 𝔽, ⊗

end # Module


text:
using Pkg
Pkg.add(["Pluto", "PlutoUI", "Plots"]) # "Images", "ImageIO", "ImageMagick"])
Pkg.precompile()


text:
using HScss

#cd(Base.Filesystem.homedir() * "/.julia/packages/HScss/3YW9k")
s = open(Base.Filesystem.homedir() * "/.julia/packages/HScss/3YW9k/src/test/testPage.Scss") do file
    String(read(file))
end
println("\nParse file contents")
vars, classes, body = hscss(s)

println("\ntranslate data structure to HScss and print")
printScss(vars, classes, body)

println("\nget head element")
print( hscss_get(body, "head"))

println("\nget all divs")
print( hscss_getall(body, "div") )

# This will be a future work... maybe
#compile(s)


text:
using Test: @test, @testset
using Distributions: Uniform, Normal
using Statistics: mean
using LinearAlgebra: norm
using PERK: perk, GaussianKernel, EuclideanKernel, GaussianRFF

function test_perk_1()

    rng = StableRNG(0)
    f = x -> exp(-30 / x)
    xtrue = 100
    y = f(xtrue)
    T = 200
    xDists = Uniform(10, 500)
    noiseDist = Normal(0, 0.01)
    signalModels = f
    λ = 2.0^-1.5
    kernel = GaussianKernel([λ * mean(y)])
    ρ = 2.0^-20
    xhat = perk(rng, y, T, xDists, noiseDist, signalModels, kernel, ρ)

    error_rel = abs(xhat - xtrue) / xtrue
    return isapprox(error_rel, 0.042677934398487306, atol = 1e-7)

end

function test_perk_2()

    rng = StableRNG(0)
    f = x -> exp(-30 / x)
    xtrue = 20:490
    N = 1
    T = 200
    xDists = Uniform(10, 500)
    noiseDist = Normal(0, 0.01)
    signalModels = f
    λ = 2.0^-1.5
    ρ = 2.0^-20
    error_rel = zeros(length(xtrue))
    for i = 1:length(xtrue)
        y = fill(f(xtrue[i]), 1, N)
        kernel = GaussianKernel([λ * mean(y)])
        xhat = perk(rng, y, T, xDists, noiseDist, signalModels, kernel, ρ)
        error_rel[i] = abs(xhat[] - xtrue[i]) / xtrue[i]
    end

    error_rel_avg = sum(error_rel) / length(error_rel)
    return isapprox(error_rel_avg, 0.043934535569840415, atol = 1e-7)

end

function test_perk_3()

    rng = StableRNG(0)
    f = x -> exp(-30 / x)
    xtrue = 100
    N = 1
    y = fill(f(xtrue), 1, N)
    T = 200
    xDists = [Uniform(10, 500)]
    noiseDist = Normal(0, 0.01)
    signalModels = f
    λ = 2.0^-1.5
    H = 100
    kernel = GaussianRFF(H, [λ * mean(y)])
    ρ = 2.0^-20
    xhat = perk(rng, y, T, xDists, noiseDist, signalModels, kernel, ρ)

    error_rel = abs(xhat[] - xtrue) / xtrue
    return isapprox(error_rel, 0.07054474887124002, atol = 1e-7)

end

function test_perk_4()

    rng = StableRNG(0)
    f = x -> exp(-30 / x)
    xtrue = 20:490
    T = 200
    xDists = [Uniform(10, 500)]
    noiseDist = Normal(0, 0.01)
    signalModels = f
    λ = 2.0^-1.5
    ρ = 2.0^-20
    H = 100
    error_rel = zeros(length(xtrue))
    for i = 1:length(xtrue)
        y = f(xtrue[i])
        kernel = GaussianRFF(H, [λ * mean(y)])
        xhat = perk(rng, y, T, xDists, noiseDist, signalModels, kernel, ρ)
        error_rel[i] = abs(xhat[] - xtrue[i]) / xtrue[i]
    end

    error_rel_avg = sum(error_rel) / length(error_rel)
    return isapprox(error_rel_avg, 0.05827088471817421, atol = 1e-7)

end

function test_perk_5()

    rng = StableRNG(0)
    f = (x, ν) -> exp(-ν / x)
    xtrue = 100
    ν = 30
    N = 1
    y = fill(f(xtrue, ν), 1, N)
    T = 200
    xDists = [Uniform(10, 500)]
    νDists = Uniform(30, 30.000000001)
    noiseDist = Normal(0, 0.01)
    signalModels = [f]
    λ = 2.0^-1.5
    kernel = GaussianKernel([λ * mean(y); λ * mean(ν)])
    ρ = 2.0^-20
    xhat = perk(rng, y, ν, T, xDists, νDists, noiseDist, signalModels, kernel, ρ)

    error_rel = abs(xhat[] - xtrue) / xtrue
    return isapprox(error_rel, 0.0057175117436099755, atol = 1e-7)

end

function test_perk_6()

    rng = StableRNG(0)
    f = (x, ν) -> [x + ν, exp(-ν / x)]
    xtrue = 100
    ν = 30
    D = 2
    N = 5
    y = f(xtrue, ν) .+ zeros(D, N)
    T = 200
    H = 100
    kernel = GaussianRFF(H, [vec(mean(y, dims = 2)); ν])
    xDists = Uniform(10, 500)
    νDists = Uniform(30, 30.000000001)
    noiseDist = Normal(0, 0.01)
    signalModels = f
    ρ = 0.001
    xhat = perk(rng, y, fill(ν, 1, N), T, xDists, νDists, noiseDist,
        signalModels, kernel, ρ)

    error_rel = norm(xhat .- xtrue) / (sqrt(N) * xtrue)
    return isapprox(error_rel, 0.0009325470666789215, atol = 1e-7)

end

function test_perk_7()

    rng = StableRNG(0)
    f = (x, ν) -> ν * log(x)
    xtrue = 10
    ν = ones(1, 1)
    y = f(xtrue, ν[])
    T = 200
    kernel = GaussianKernel([1, 1])
    xDists = Uniform(0, 20)
    νDists = [Uniform(0, 2)]
    noiseDist = Normal(0, 0.01)
    signalModels = [f]
    ρ = 0.0001
    xhat = perk(rng, y, ν, T, xDists, νDists, noiseDist, signalModels, kernel, ρ)

    error_rel = abs(xhat[] - xtrue) / xtrue
    return isapprox(error_rel, 0.1083022379521612, atol = 1e-7)

end

function test_perk_8()

    rng = StableRNG(0)
    f = (x, ν) -> ν * log(x)
    xtrue = 10
    ν = 1
    y = f(xtrue, ν)
    T = 200
    kernel = GaussianKernel([1, 1])
    xDists = Uniform(0, 20)
    νDists = Uniform(0, 2)
    noiseDist = Normal(0, 0.01)
    signalModels = [f]
    ρ = 0.0001
    xhat = perk(rng, y, ν, T, xDists, νDists, noiseDist, signalModels, kernel, ρ)

    error_rel = abs(xhat[] - xtrue) / xtrue
    return isapprox(error_rel, 0.1083022379521612, atol = 1e-7)

end

function test_perk_9()

    rng = StableRNG(0)
    f = x -> x^2
    xtrue = 4
    y = fill(f(xtrue), 1, 1)
    T = 200
    kernel = EuclideanKernel()
    xDists = Uniform(1, 10)
    noiseDist = Normal(0, 0.01)
    signalModels = [f]
    ρ = 0.001
    xhat = perk(rng, y, T, xDists, noiseDist, signalModels, kernel, ρ)

    error_rel = abs(xhat[] - xtrue) / xtrue
    return isapprox(error_rel, 0.09566274665119057, atol = 1e-6)

end

function test_perk_10()

    rng = StableRNG(0)
    f = (x, ν) -> x * ν
    xtrue = 5.5
    N = 100
    ν = fill(2, N)
    y = f.(xtrue, ν) .+ 0.01 .* randn(rng, N)
    T = 200
    kernel = GaussianKernel([mean(y), mean(ν)])
    xDists = [Uniform(1, 10)]
    νDists = [Uniform(0, 3)]
    noiseDist = Normal(0, 0.01)
    signalModels = f
    ρ = 0.01
    xhat = perk(rng, y, ν, T, xDists, νDists, noiseDist, signalModels, kernel, ρ)

    error_rel = norm(xhat .- xtrue) / (sqrt(N) * xtrue)
    return isapprox(error_rel, 0.060384227201893494, atol = 1e-7)

end

function test_perk_11()

    rng = StableRNG(0)
    f = (x, ν) -> x * ν
    xtrue = 5.5
    N = 100
    ν = fill(2, N)
    y = f.(xtrue, ν) .+ 0.01 .* randn(rng, N)
    T = 200
    kernel = GaussianKernel([mean(y), mean(ν)])
    xDists = [Uniform(1, 10)]
    νDists = Uniform(0, 3)
    noiseDist = Normal(0, 0.01)
    signalModels = f
    ρ = 0.01
    xhat = perk(rng, y, ν, T, xDists, νDists, noiseDist, signalModels, kernel, ρ)

    error_rel = norm(xhat .- xtrue) / (sqrt(N) * xtrue)
    return isapprox(error_rel, 0.060384227201893494, atol = 1e-7)

end

function test_perk_12()

    rng = StableRNG(0)
    f = (x, ν) -> x * ν
    xtrue = 5.5
    N = 100
    ν = fill(2, N)
    y = f.(xtrue, ν) .+ 0.01 .* randn(rng, N)
    T = 200
    kernel = GaussianKernel([mean(y), mean(ν)])
    xDists = Uniform(1, 10)
    νDists = [Uniform(0, 3)]
    noiseDist = Normal(0, 0.01)
    signalModels = f
    ρ = 0.01
    xhat = perk(rng, y, ν, T, xDists, νDists, noiseDist, signalModels, kernel, ρ)

    error_rel = norm(xhat .- xtrue) / (sqrt(N) * xtrue)
    return isapprox(error_rel, 0.06038422720189351, atol = 1e-7)

end

@testset "PERK" begin

    @test test_perk_1()
    @test test_perk_2()
    @test test_perk_3()
    @test test_perk_4()
    @test test_perk_5()
    @test test_perk_6()
    @test test_perk_7()
    @test test_perk_8()
    @test test_perk_9()
    @test test_perk_10()
    @test test_perk_11()
    @test test_perk_12()

end


text:
using MySQL_Octo

HOST = "127.0.0.1"

if true
    adapter = SQL.connect(host=HOST, user="root", pass="", reconnect=true)
    SQL.execute(adapter, """
    DROP DATABASE IF EXISTS mysqltest;
    GRANT USAGE ON *.* TO 'test'@'127.0.0.1' IDENTIFIED BY 'test';
    DROP USER 'test'@'127.0.0.1';
    """)
    SQL.execute(adapter, """
    CREATE DATABASE mysqltest;
    CREATE USER test@$HOST IDENTIFIED BY 'test';
    GRANT ALL ON mysqltest.* TO 'test'@$HOST;
    FLUSH PRIVILEGES;
    """)
    SQL.disconnect(adapter)
end

adapter = SQL.connect(host=HOST, user="test", pass="test", db="mysqltest", reconnect=true)

SQL.execute(adapter, """
CREATE TABLE users (
    id INT NOT NULL AUTO_INCREMENT,
    name VARCHAR(20),
    PRIMARY KEY (id)
);
CREATE TABLE comments (
    id INT NOT NULL AUTO_INCREMENT,
    user_id INT,
    body VARCHAR(300),
    FOREIGN KEY (user_id) REFERENCES users(id),
    PRIMARY KEY (id)
); """)

SQL.execute(adapter, "INSERT INTO users (name) VALUES ('John');")


text:
struct BanditExperience{T,TA} <: Any where {T,TA}
    action::TA
    blogp::T
    reward::T
end

struct Trajectory{T,TS,TA} <: Any where {T,TS,TA}
    states::Array{TS,1}
    actions::Array{TA,1}
    blogps::Array{T,1}
    rewards::Array{T,1}

    function Trajectory(::Type{T}, ::Type{TS}, ::Type{TA}) where {T,TS,TA}
        new{T,TS,TA}(Array{TS,1}(), Array{TA,1}(), Array{T,1}(), Array{T,1}())
    end
end

function length(τ::Trajectory)
    return length(τ.rewards)
end

function push!(τ::Trajectory{T,TS,TA}, state::TS, action::TA, blogp::T, reward::T) where {T,TS,TA}
    push!(τ.states, state)
    push!(τ.actions, action)
    push!(τ.blogps, blogp)
    push!(τ.rewards, reward)
end


text:
# GFF3 Reader

struct Reader{S <: TranscodingStream} <: BioGenerics.IO.AbstractReader
    io::S
end

"""
    GFF.Reader(input::IO)

Create a data reader of the GFF3 file format.
"""
function Reader(input::IO)
    if input isa TranscodingStream
        Reader(input)
    else
        stream = TranscodingStreams.NoopStream(input)
        Reader(stream)
    end
end

function Base.eltype(::Type{<:Reader})
    return Record
end

function Base.close(reader::Reader)
    close(reader.io)
end

function Base.open(::Type{<:Reader}, input::AbstractString)
    if input[end-2:end] == ".gz"
        return Reader(GzipDecompressorStream(open(input)))
    else
        return Reader(TranscodingStreams.NoopStream(open(input)))
    end
end

function BioGenerics.IO.stream(reader::Reader)
    reader.io
end

function Base.iterate(reader::Reader, nextone::Record = Record())
    if BioGenerics.IO.tryread!(reader, nextone) === nothing
        return nothing
    end
    return nextone, Record()
end

function BioGenerics.IO.tryread!(reader::Reader, output)
    parsechromosome!(reader.io, output)
end

function parsechromosome!(input, record::Record{G}) where G <: AbstractGene
    record.genedata[!, :source] = Union{Missing, String}[]
    record.genedata[!, :score] = Union{Missing, Float64}[]
    record.genedata[!, :phase] = Union{Missing, Int}[]
    iobuffer = IOBuffer()
    isheader = true
    isfooter = false
    qualifier = String("")
    content = String("")
    header = IOBuffer()
    currentfasta = ""
    linecount = 0
    while !eof(input)
        linecount += 1
        line = readline(input)
        if length(line) == 0
            continue
        end
        # Store the header
        if line[1] == '#' && line != "##FASTA"
            isheader = true
            println(header, line)
        else
            isheader = false
        end
        if line == "##FASTA"
            isfooter = true
        elseif !isfooter && !isheader
            (seqid, source, feature, sstart, send, score, strand, phase, attributes) = split(line, '\t')
            if isempty(record.name)
                record.header = String(take!(header))
                record.name = seqid
            elseif seqid != record.name
                ### The following contig lacks a header
                TranscodingStreams.unread(input, UInt8.(c for c in line))
                return record
            end
            locus = Locus(parse(Int, sstart):parse(Int, send), strand[1])
            addgene!(record, Symbol(feature), locus)
            if source != "." && ismissing(record.genes[end].source)
                record.genes[end].source = source
            end
            if score != "." && ismissing(record.genes[end].score)
                record.genes[end].score = parse(Float64, score)
            end
            if phase != "."
                record.genes[end].phase = parse(Int, phase)
            end
            if attributes != "."
                for attribute in split(attributes, ';')
                    isempty(attribute) && continue
                    qualifier, values = split(attribute, '=')
                    values = split(values, ',')
                    for value in values
                        if occursin(r"^\d+$", value)
                            value = parse(Int, value)
                        else
                            value = String(value)
                        end
                        pushproperty!(record.genes[end], Symbol(qualifier), value)
                    end
                end
            end
        elseif isfooter
            if line[1] == '>'
                currentfasta = line[2:end]
                seq = String(take!(iobuffer))
                record.sequence = LongDNASeq(seq)
            else
                if line[1] ∉ ['A', 'T', 'G', 'C']
                    println(line)
                end
                print(iobuffer, line)
            end
        end
    end
    if !isempty(currentfasta)
        record.sequence = LongDNASeq(String(take!(iobuffer)))
    end
    record.header = headerstring = String(take!(header))
    if isempty(record.name)
        return nothing
    else
        return record
    end
end


text:
# This is a simple and naive tool for stripping metadata information (dbg, tbaa,
# etc.) from LLVM assembly code as produced by Julia's LLVM code generator.
# Takes as input the name of an assembly (`.ll`) file and prints the
# transformed code to STDOUT.

llvm_ir = readall(ARGS[1])
# removes metadata from function signatures
llvm_ir = replace(llvm_ir, r"\) \#.*\{", "\) \{")
# removes calls to `llvm.dbg.declare` and `llvm.dbg.value`
llvm_ir = replace(llvm_ir, r".*call void .*llvm\.dbg.*\n", "")
# removes metadata which is passed to instructions
llvm_ir = replace(llvm_ir, r", \!dbg.*", "")
# removes tbaa nodes
llvm_ir = replace(llvm_ir, r", \!tbaa.*", "")
print(llvm_ir)


text:
using BayesNets
using Discretizers
using JLD

basepath = "./results/ByIntersection/"
doingPeach = false
testnums = ["111","011"]#"000","001","010","011", "100"]
intersections = collect(1:9)
test_intersections = collect(1:9)
doSubtest = false #2500 features only

function makePaths(basepath, testnums, intersections)
    paths = []
    for testnum in testnums
        for intersection in intersections
            push!(paths, string(basepath, testnum, "/", intersection))
        end
    end
    return paths
end

paths = makePaths(basepath, testnums, intersections)

typeToEdges = Dict(
    "V"=>[0.0,0.5,20,40,60],
    "A"=>[-5000,-20,-5,-0.5,0.5,5,20,5000],
    "dist" => [0.0, 5, 25, 100, 500, 10000],
    "hdwy" => [0.0, 5, 25, 100, 250, 10000],)
idToType = Dict(
    3=>"V",
    4=>"A",
    5=> "yaw",
    #6 is indicator
    7=> "hdwy",
    8=> "dist",
)

#returns allFeatures, allTargets dictionaries wheren "testnum/intersection" are the keys
function loadAllFeaturesTargets(paths)
    allFeatures = Dict()
    allTargets = Dict()
    allFids = Dict()
    nFids = 0
    numFeatures = 0
    for path in paths #path is result/ByIntersection/yyy/x
        i = path[end-4:end]
        f = open("$path/featuresAndTargets")
        allFeatures[i] = readdlm(f)
        allTargets[i] = allFeatures[i][:,end]  #nextMove is the last thing only
        allFeatures[i] = allFeatures[i][:,1:end-3]

        if doSubtest
            allFeatures[i] = allFeatures[i][1:2500,:]
            allTargets[i] = allTargets[i][1:2500,:]
        end
        close(f)

        if numFeatures == 0
            numFeatures = length(allFeatures[i][1,:])
        end
    end
    return allFeatures, allTargets, numFeatures
end

function makeDiscretizedThings(trainlines, testlines, trainactuals, testactuals, numFeatures,
                                           indexToID, numYawBins, moveIndex)
  allLines = vcat(trainlines, testlines)
  allactuals = vcat(trainactuals, testactuals)
  discAllLines = zeros(Int64, (length(trainlines[:,1])+length(testlines[:,1]), numFeatures+1))
  discTrainlines = zeros(Int64, (length(trainlines[:,1]), numFeatures+1))
  discTestlines = zeros(Int64, (length(testlines[:,1]), numFeatures+1))
  for k in 1:numFeatures
    id = indexToID[k]
    if haskey(idToType, id) #if not, already discretized
      F_type = idToType[id]
      if haskey(typeToEdges, F_type)
        Disc = LinearDiscretizer(typeToEdges[F_type])
      else
        Disc = LinearDiscretizer(binedges(DiscretizeUniformWidth(numYawBins),allLines[:,k]))
      end
      DiscretizedAll = encode(Disc, allLines[:,k])
      DiscretizedTrain = encode(Disc, trainlines[:,k])
      DiscretizedTest = encode(Disc, testlines[:,k])
    else
      #println("index ", k, " has no id, aka already disc")
      DiscretizedAll = allLines[:,k]
      DiscretizedTrain = trainlines[:,k]
      DiscretizedTest = testlines[:,k]
    end
    discTrainlines[:,k] = encode(CategoricalDiscretizer(DiscretizedAll), DiscretizedTrain)
    discTestlines[:,k] = encode(CategoricalDiscretizer(DiscretizedAll), DiscretizedTest)
    discAllLines[:,k] = encode(CategoricalDiscretizer(DiscretizedAll), DiscretizedAll)
  end
  moveDiscretizer = CategoricalDiscretizer(trainactuals)
  discTrainlines[:,moveIndex] = encode(moveDiscretizer,trainactuals)
  discTestlines[:,moveIndex] = encode(moveDiscretizer,testactuals)
  discAllLines[:,moveIndex] = encode(moveDiscretizer,allactuals)

  alldata = convert(DataFrame, discAllLines) #used to find num bins
  traindata = convert(DataFrame, discTrainlines)
  testdata = convert(DataFrame, discTestlines)
  return alldata, traindata, testdata, moveDiscretizer
end

function makeIndexToID(numFeatures, testnum, laneTypeEncodedLen)
    useLaneType = (testnum[1] == '1')
    useHist = (testnum[2] == '1')
    useTraffic = (testnum[3] == '1')
    println(testnum, useLaneType, useHist, useTraffic)
    indexToID = zeros(Int64, (1, numFeatures+1+laneTypeEncodedLen))   #+1 for when no lane type, shift needs buffer
    coreLen = 8
    if useLaneType
        coreLen += laneTypeEncodedLen
    end
    histLen = 9 #indicator is first
    numHist = 4
    totHist = histLen * numHist
    traffLen = 8
    #base = lane, lane, v, a, yaw, 1{hdwy}, hdwy, dist
    #hist is the same, 4 times
    #traffic is dx, dy, v, a, yaw, 1{hdwy}, hdwy
    for index in 1:numFeatures
        ID = index
        if index <= coreLen
          if useLaneType
            if index > 4
              ID -= laneTypeEncodedLen
            else #laneTypeEncoding features are 0 ID
              ID = 0
            end
          end
          ID += 1 #because all others get shifted, limits amount of confusion i think
        elseif useHist == false
          ID = ((index - coreLen) % traffLen)
          if ID == 0
            ID = traffLen
          elseif ID == 2 || ID == 3
            ID = 9 #distance + 1
          end
        elseif useHist && index <= coreLen + totHist
            ID = ((index - coreLen) % histLen)
            if ID == 0
                ID = histLen
            end
        else  #use hist and in traffic section
          ID = ((index - (coreLen + totHist)) % traffLen)
          if ID == 0
            ID = traffLen
          elseif ID == 2 || ID == 3
            ID = 9 #distance + 1
          end
        end
        ID -= 1 #shift for indicator at front of each traff
        indexToID[index] = ID
    end
    #for i in 1:numFeatures
    #    id = indexToID[i]
    #    println("index: ",i, "id: ", id)
    #    if haskey(idToType, id)
    #        println("key:", idToType[id])
    #    end
    #end
    return indexToID
end

# fit function --- modified from https://github.com/sisl/BayesNets.jl/blob/master/src/DiscreteBayesNet/greedy_hill_climbing.jl
function Distributions.fit(::Type{DiscreteBayesNet}, data::DataFrame, params::GreedyHillClimbing, indexOfOut::Int;
    ncategories::Vector{Int} = map!(i->infer_number_of_instantiations(data[i]), Array(Int, ncol(data)), 1:ncol(data)),
    )

    n = ncol(data)
    parent_list = map!(i->Int[], Array(Vector{Int}, n), 1:n)
    datamat = convert(Matrix{Int}, data)'
    score_components = bayesian_score_components(parent_list, ncategories, datamat, params.prior, params.cache)

    while true
      best_diff = 0.0
      best_parent_list = parent_list
      i = indexOfOut
      parent_list[i] = collect(1:n)
      break
      #for i in 1:n #this is added for full fitting
        # 1) add an edge (j->i)
        if length(parent_list[i]) < params.max_n_parents
            for j in deleteat!(collect(1:n), parent_list[i])
                    if adding_edge_preserves_acyclicity(parent_list, j, i)
                        new_parents = sort!(push!(copy(parent_list[i]), j))
                        new_component_score = bayesian_score_component(i, new_parents, ncategories, datamat, params.prior, params.cache)
                    if new_component_score - score_components[i] > best_diff
                        best_diff = new_component_score - score_components[i]
                        best_parent_list = deepcopy(parent_list)
                        best_parent_list[i] = new_parents
                    end
                end
             end
        end

        # 2) remove an edge
        for (idx, j) in enumerate(parent_list[i])

            new_parents = deleteat!(copy(parent_list[i]), idx)
            new_component_score = bayesian_score_component(i, new_parents, ncategories, datamat, params.prior, params.cache)
            if new_component_score - score_components[i] > best_diff
                best_diff = new_component_score - score_components[i]
                best_parent_list = deepcopy(parent_list)
                best_parent_list[i] = new_parents
            end
            # 3) flip an edge
            new_parent_list = deepcopy(parent_list) # TODO: make this more efficient
            deleteat!(new_parent_list[i], idx)

            if adding_edge_preserves_acyclicity(new_parent_list, i, j)
                sort!(push!(new_parent_list[j], i))
                new_diff = bayesian_score_component(i, new_parent_list[i], ncategories, datamat, params.prior,
                                                       params.cache) - score_components[i]
                new_diff += bayesian_score_component(j, new_parent_list[j], ncategories, datamat, params.prior,
                                                       params.cache) - score_components[j]
                if new_diff > best_diff
                    best_diff = new_diff
                    best_parent_list = new_parent_list
                end
            end
        end

        if best_diff > 0.0
            parent_list = best_parent_list
            score_components = bayesian_score_components(parent_list, ncategories, datamat, params.prior, params.cache)
        else
            break
        end
      #end   #this is for the "for i in 1:n"
    end
    # construct the BayesNet
    cpds = Array(DiscreteCPD, n)
    varnames = names(data)
    for j in 1:n
        name = varnames[j]
        parents = varnames[parent_list[j]]
        cpds[j] = Distributions.fit(DiscreteCPD, data, name, parents, params.prior,
                      parental_ncategories=ncategories[parent_list[j]],
                      target_ncategories=ncategories[j])
    end
    #println(cpds)
    println("Parent List: ", parent_list)
    BayesNet(cpds)
end

#big for loop
println("Loading features/targets for BN")
allFeatures, allTargets, numFeatures = loadAllFeaturesTargets(paths)
println("Done loading features/targets for BN")
println("NumFeatures == ", numFeatures)
for inter in test_intersections
  for testnum in testnums
      i = string(testnum, "/", inter)
      scorename = "Ypred_BN"
      scorePath = string("$basepath",testnum, "/TestOn$inter","/$scorename")
      f = open(scorePath, "w")
      #truncate score file
      close(f)
      moveCats = [1,2,3]
      numNextMoves = length(moveCats)
      useLaneType = false
      if testnum[1] == "1"
        useLaneType = true
      end
      laneTypeIndex = 3
      laneTypeEncodedLen = length([0,0,0,0])
      numYawBins = 7
      #truncate the log file
      numFeatures = length(allFeatures[i][1,:])
      indexToID = makeIndexToID(numFeatures, testnum, laneTypeEncodedLen)

      println("index to ID dict: \n ",indexToID)

      moveIndex = numFeatures+1

      println("Done making discretizer dictionaries")
      println("Starting feature/target formatting for CV #: $i")
      testlines = allFeatures[i]
      testactuals = allTargets[i]

      trainlines = zeros(Float64, (1,numFeatures))
      trainactuals = zeros(Float64, (1,numFeatures))
      first = true
      for other_inter in intersections
        if other_inter != inter  #the train lines are all that are not test
          j = string(testnum, "/", other_inter)
          if first == true
            trainlines = allFeatures[j]
            trainactuals = allTargets[j]
            first = false
          else
            trainlines = vcat(trainlines, allFeatures[j])
            trainactuals = vcat(trainactuals, allTargets[j])
          end
        end
      end
      println("Done feature/target formatting for CV #: $i")

      println("Starting to discretize for CV #: $i")
      alldata, traindata, testdata, moveDiscretizer = makeDiscretizedThings(trainlines, testlines, trainactuals, testactuals,
                                                           numFeatures, indexToID, numYawBins, moveIndex)

      println("Done discretizing for CV #: $i")
      println("Starting to fit BN for CV #: $i")
      max_parents = numFeatures
      #max_parents = min(numFeatures, max_parents)
      println("Max parents: ", max_parents)

      params = GreedyHillClimbing(ScoreComponentCache(traindata), max_n_parents=max_parents, prior=UniformPrior())
      println("Done fitting params, starting to fit BN")

      num_bins_all = map!(i->infer_number_of_instantiations(alldata[i]), Array(Int, ncol(alldata)), 1:ncol(alldata))
      bnDis = @time fit(DiscreteBayesNet, traindata, params, moveIndex; ncategories=num_bins_all)
      println("Done fitting BN for CV #: $i")
      save(string("$basepath$i","/BN_model.jld"), "bnDis", bnDis)
      score = 0
      numNaN = 0
      nrows = length(testlines[:,1])
      moveCats = unique(collect(testdata[:,moveIndex]))
      numNextMoves = 3#length(moveCats)
      numRight = 0
      for index in 1:nrows
          actualNextMove = convert(Int64,decode(moveDiscretizer,testdata[index,moveIndex]))
          p_dists = zeros(Float64, (1, numNextMoves))
          for move in moveCats
            featureline = testdata[index,:]
            featureline[moveIndex] = encode(moveDiscretizer, move)
            p_dists[move] = pdf(bnDis, featureline)
          end
          for x in p_dists
            if !(x > 1) && !(x < 1) && !(x == 1)
              p_dists = [1.0/numNextMoves for i in 1:numNextMoves]
              numNaN += 1
              break
            end
          end
          p_dists = p_dists ./ sum(p_dists)
          p_right = p_dists[actualNextMove]
          score += (1 - p_right)
          if p_right == maximum(p_dists)
              numRight += 1
          end
          scoref = open(scorePath, "a")
          writedlm(scoref, p_dists)
          close(scoref)
      end
      println("Score:", score)
      println("Num Right:", numRight)
      println("Num Predicitons:", nrows)
      println("Num NaN:", numNaN)
      scoref = open("$basepath$i/$scorename", "a")
      #writedlm(scoref, score)
      #writedlm(scoref, numNaN)
      close(scoref)
      println("Done scoring BN for CV #: $i")
    end #end over testnums
    println("Done with intersection", inter)
end  #end over intersections

max_stars_repo_name: zhangzhangzf/IntentionPrediction

max_issues_repo_name: zhangzhangzf/IntentionPrediction

max_forks_repo_name: zhangzhangzf/IntentionPrediction


text:
# Autogenerated wrapper script for Sundials_jll for x86_64-apple-darwin-libgfortran3
export libsundials_arkode, libsundials_cvode, libsundials_cvodes, libsundials_ida, libsundials_idas, libsundials_kinsol, libsundials_nvecmanyvector, libsundials_nvecserial, libsundials_sunlinsolband, libsundials_sunlinsoldense, libsundials_sunlinsolklu, libsundials_sunlinsollapackband, libsundials_sunlinsollapackdense, libsundials_sunlinsolpcg, libsundials_sunlinsolspbcgs, libsundials_sunlinsolspfgmr, libsundials_sunlinsolspgmr, libsundials_sunlinsolsptfqmr, libsundials_sunmatrixband, libsundials_sunmatrixdense, libsundials_sunmatrixsparse, libsundials_sunnonlinsolfixedpoint, libsundials_sunnonlinsolnewton

using CompilerSupportLibraries_jll
using OpenBLAS_jll
using SuiteSparse_jll
JLLWrappers.@generate_wrapper_header("Sundials")
JLLWrappers.@declare_library_product(libsundials_arkode, "@rpath/libsundials_arkode.4.dylib")
JLLWrappers.@declare_library_product(libsundials_cvode, "@rpath/libsundials_cvode.5.dylib")
JLLWrappers.@declare_library_product(libsundials_cvodes, "@rpath/libsundials_cvodes.5.dylib")
JLLWrappers.@declare_library_product(libsundials_ida, "@rpath/libsundials_ida.5.dylib")
JLLWrappers.@declare_library_product(libsundials_idas, "@rpath/libsundials_idas.4.dylib")
JLLWrappers.@declare_library_product(libsundials_kinsol, "@rpath/libsundials_kinsol.5.dylib")
JLLWrappers.@declare_library_product(libsundials_nvecmanyvector, "@rpath/libsundials_nvecmanyvector.5.dylib")
JLLWrappers.@declare_library_product(libsundials_nvecserial, "@rpath/libsundials_nvecserial.5.dylib")
JLLWrappers.@declare_library_product(libsundials_sunlinsolband, "@rpath/libsundials_sunlinsolband.3.dylib")
JLLWrappers.@declare_library_product(libsundials_sunlinsoldense, "@rpath/libsundials_sunlinsoldense.3.dylib")
JLLWrappers.@declare_library_product(libsundials_sunlinsolklu, "@rpath/libsundials_sunlinsolklu.3.dylib")
JLLWrappers.@declare_library_product(libsundials_sunlinsollapackband, "@rpath/libsundials_sunlinsollapackband.3.dylib")
JLLWrappers.@declare_library_product(libsundials_sunlinsollapackdense, "@rpath/libsundials_sunlinsollapackdense.3.dylib")
JLLWrappers.@declare_library_product(libsundials_sunlinsolpcg, "@rpath/libsundials_sunlinsolpcg.3.dylib")
JLLWrappers.@declare_library_product(libsundials_sunlinsolspbcgs, "@rpath/libsundials_sunlinsolspbcgs.3.dylib")
JLLWrappers.@declare_library_product(libsundials_sunlinsolspfgmr, "@rpath/libsundials_sunlinsolspfgmr.3.dylib")
JLLWrappers.@declare_library_product(libsundials_sunlinsolspgmr, "@rpath/libsundials_sunlinsolspgmr.3.dylib")
JLLWrappers.@declare_library_product(libsundials_sunlinsolsptfqmr, "@rpath/libsundials_sunlinsolsptfqmr.3.dylib")
JLLWrappers.@declare_library_product(libsundials_sunmatrixband, "@rpath/libsundials_sunmatrixband.3.dylib")
JLLWrappers.@declare_library_product(libsundials_sunmatrixdense, "@rpath/libsundials_sunmatrixdense.3.dylib")
JLLWrappers.@declare_library_product(libsundials_sunmatrixsparse, "@rpath/libsundials_sunmatrixsparse.3.dylib")
JLLWrappers.@declare_library_product(libsundials_sunnonlinsolfixedpoint, "@rpath/libsundials_sunnonlinsolfixedpoint.2.dylib")
JLLWrappers.@declare_library_product(libsundials_sunnonlinsolnewton, "@rpath/libsundials_sunnonlinsolnewton.2.dylib")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll, OpenBLAS_jll, SuiteSparse_jll)
    JLLWrappers.@init_library_product(
        libsundials_arkode,
        "lib/libsundials_arkode.4.2.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_cvode,
        "lib/libsundials_cvode.5.2.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_cvodes,
        "lib/libsundials_cvodes.5.2.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_ida,
        "lib/libsundials_ida.5.2.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_idas,
        "lib/libsundials_idas.4.2.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_kinsol,
        "lib/libsundials_kinsol.5.2.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_nvecmanyvector,
        "lib/libsundials_nvecmanyvector.5.2.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_nvecserial,
        "lib/libsundials_nvecserial.5.2.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsolband,
        "lib/libsundials_sunlinsolband.3.2.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsoldense,
        "lib/libsundials_sunlinsoldense.3.2.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsolklu,
        "lib/libsundials_sunlinsolklu.3.2.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsollapackband,
        "lib/libsundials_sunlinsollapackband.3.2.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsollapackdense,
        "lib/libsundials_sunlinsollapackdense.3.2.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsolpcg,
        "lib/libsundials_sunlinsolpcg.3.2.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsolspbcgs,
        "lib/libsundials_sunlinsolspbcgs.3.2.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsolspfgmr,
        "lib/libsundials_sunlinsolspfgmr.3.2.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsolspgmr,
        "lib/libsundials_sunlinsolspgmr.3.2.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunlinsolsptfqmr,
        "lib/libsundials_sunlinsolsptfqmr.3.2.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunmatrixband,
        "lib/libsundials_sunmatrixband.3.2.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunmatrixdense,
        "lib/libsundials_sunmatrixdense.3.2.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunmatrixsparse,
        "lib/libsundials_sunmatrixsparse.3.2.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunnonlinsolfixedpoint,
        "lib/libsundials_sunnonlinsolfixedpoint.2.2.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsundials_sunnonlinsolnewton,
        "lib/libsundials_sunnonlinsolnewton.2.2.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function update_zeta!(state::State, data::Data, prior::Prior)
  for i in samplenames
    update_zeta!(i, state, data, prior)
  end
end

function update_zeta!(i::Char, state::State, data::Data, prior::Prior)
  yi = ref_yfinite(data, i)
  zetai = ref_zeta(state, i)
  lami = ref_lambda(state, i)
  vi = ref_v(state, i)

  psi = state.psi[lami]
  omega = state.omega[lami]
  mu = state.mu[lami]

  vnew = 1 ./ (vi + (psi .^ 2) .* vi ./ omega)
  mnew = vnew .* vi .* psi .* (yi - mu) ./ omega

  for n in eachindex(zetai)
    zetai[n] = rand(truncated(Normal(mnew[n], sqrt(vnew[n])), 0, Inf))
  end
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation