text:
include_ddp()

# Model
include_model("car")
n = model.n
m = model.m

# Time
T = 101
h = 0.05

# Reference trajectory
xx = range(0.0, stop = 1.0, length = T)
yy = 0.1 * sin.(2 * π * xx)
θθ = 0.2 * π * cos.(2 * π * xx)
plot(xx, yy, color = :black)
scatter!(xx, yy, color = :black)

p_ref = [[xx[t]; yy[t]] for t = 1:T]

# Initial conditions, controls, disturbances
x1 = [xx[1]; yy[1]; θθ[1]]#[0.0, 0.0, 0.0]
xT = [xx[T]; yy[T]; θθ[T]]#[1.0, 0.0, 0.0] # goal state
ū = [1.0e-3 * rand(model.m) for t = 1:T-1]
w = [zeros(model.d) for t = 1:T-1]

# Rollout
x̄ = rollout(model, x1, ū, w, h, T)

# Objective
Q = [(t < T ? 0.0 * Diagonal([10.0; 10.0; 0.1])
        : 0.0 * Diagonal([10.0; 10.0; 0.1])) for t = 1:T]
q = [-2.0 * Q[t] * [xx[t]; yy[t]; θθ[t]]  for t = 1:T]

R = [Diagonal(1.0e-3 * ones(model.m)) for t = 1:T-1]
r = [zeros(model.m) for t = 1:T-1]

obj = StageCosts([QuadraticCost(Q[t], q[t],
	t < T ? R[t] : nothing, t < T ? r[t] : nothing) for t = 1:T], T)

function g(obj::StageCosts, x, u, t)
	T = obj.T
    if t < T
		Q = obj.cost[t].Q
		q = obj.cost[t].q
	    R = obj.cost[t].R
		r = obj.cost[t].r
        return x' * Q * x + q' * x + u' * R * u + r' * u
    elseif t == T
		Q = obj.cost[T].Q
		q = obj.cost[T].q
        return x' * Q * x + q' * x
    else
        return 0.0
    end
end

# Constraints
p = [t < T ? 2 * m + 2 : n for t = 1:T]
info_t = Dict(:p => p_ref, :ul => [-10.0; -10.0], :uu => [10.0; 10.0], :inequality => (1:2 * m))
info_T = Dict(:xT => xT)
con_set = [StageConstraint(p[t], t < T ? info_t : info_T) for t = 1:T]

function c!(c, cons::StageConstraints, x, u, t)
	T = cons.T
	p = cons.con[t].p

	if t < T
		ul = cons.con[t].info[:ul]
		uu = cons.con[t].info[:uu]
		c[1:(2 * m)] .= [ul - u; u - uu]

		if t % 10 == 0
			pt = cons.con[t].info[:p]
			c[(2 * m) .+ (1:2)] .= x[1:2] - pt[t]
		end
	elseif t == T
		xT = cons.con[T].info[:xT]
		c .= x - xT
	else
		nothing
	end
end

prob = problem_data(model, obj, con_set, copy(x̄), copy(ū), w, h, T)

# Solve
@time constrained_ddp_solve!(prob,
    max_iter = 1000, max_al_iter = 10,
	con_tol = 1.0e-4,
	ρ_init = 1.0, ρ_scale = 10.0)

x, u = current_trajectory(prob)
x̄, ū = nominal_trajectory(prob)

# Visualize
# using Plots
# plot(π * ones(T),
#     width = 2.0, color = :black, linestyle = :dash)
# plot!(hcat(x...)', width = 2.0, label = "")
#
# plot(hcat([con_set[1].info[:ul] for t = 1:T]...)',
#     width = 2.0, color = :black, label = "")
# plot!(hcat([con_set[1].info[:uu] for t = 1:T]...)',
#     width = 2.0, color = :black, label = "")
# plot!(hcat(u..., u[end])',
#     width = 2.0, linetype = :steppost,
# 	label = "", color = :orange)

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
# open(vis)
visualize!(vis, model, x, Δt = h)

plot(xx, yy, color = :black, width = 2.0)
plot!([xt[1] for xt in x], [xt[2] for xt in x], color = :cyan, width = 1.0)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
using Test
using MOID
using Printf

@testset "Test MOID" begin

   SHOW = true

   rad = pi / 180

   elem = [ [" 1 Ceres"     , 2.7691652, 0.0760091,  73.59764,  80.30553,10.59407],
            ["29 Amphitrite", 2.5541136, 0.0726956,  63.36319, 356.34176, 6.08252],
            ["30 Urania"    , 2.3655722, 0.127581 ,  87.42605, 307.46872, 2.09575],
            ["50 Virginia"  , 2.6487939, 0.2859856, 200.08054, 173.52874, 2.83822],
            ["51 Nemausa"   , 2.3658354, 0.0675594,   2.58053, 175.9785 , 9.97718], ];

   # Expected results
   #Asteroid1  Asteroid2         MOID [AU]
   # 1 Ceres    1 Ceres          0.00000000000001
   # 1 Ceres   29 Amphitrite     0.15677463452737
   # 1 Ceres   30 Urania         0.24521440655832
   # 1 Ceres   50 Virginia       0.08934734026105
   # 1 Ceres   51 Nemausa        0.35972678460706

   result = [ 0.00000000000001, 0.15677463452737, 0.24521440655832, 0.08934734026105, 0.35972678460706 ]

   SHOW && println("\nAsteroid1  Asteroid2         MOID [AU]")

   for i in 1:length(elem)
      moid = moidF(elem[1][2:6]...,elem[i][2:6]...)
      @test moid ≈ result[i] atol=1E-14
      SHOW && @printf("%s   %-15s   %0.14f\n", elem[1][1], elem[i][1], moid)
   end

   SHOW && println()

end

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

ext: jl


text:
# Archimedean copulas
"""
  getV0(θ::Union{Float64, Int}, v::Vector{Float64}, copula::String)

Returns Vector{Float} or Vector{Int} of realisations of axiliary variable V0
used to ganarate data from 1d archimedean copula with parameter Θ, given v:
realisations of 1d variable uniformly distributed on [0,1]

```jldoctest

julia> getV0(2., [0.2, 0.4, 0.6, 0.8], "clayton")
4-element Array{Float64,1}:
 0.0641848
 0.274996
 0.708326
 1.64237
```
"""
function getV0(θ::Union{Float64, Int}, v::Vector{Float64}, copula::String)
  if copula == "clayton"
    return quantile.(Gamma(1/θ, 1), v)
  elseif copula == "amh"
    return 1 .+quantile.(Geometric(1-θ), v)
  elseif copula == "frank"
    return logseriesquantile(1-exp(-θ), v)
  elseif copula == "gumbel"
    return levygen(θ, v)
  end
  throw(AssertionError("$(copula) not supported"))
end

"""
  phi(u::Matrix{Float64}, θ::Union{Float64, Int}, copula::String)

Given a matrix t realisations of n variate data ℜᵗⁿ ∋ u = -log(rand(t,n))./V0
returns it transformed through an inverse generator of archimedean copula. Output
is distributed uniformly on [0,1]ⁿ

```jldoctest

julia> julia> phi([0.2 0.6; 0.4 0.8], 2., "clayton")
2×2 Array{Float64,2}:
 0.845154  0.6742
 0.745356  0.620174
```
"""
function phi(u::Matrix{Float64}, θ::Union{Float64, Int}, copula::String)
  if copula == "clayton"
    return (1 .+ u).^(-1/θ)
  elseif copula == "amh"
    return (1-θ) ./(exp.(u) .-θ)
  elseif copula == "frank"
    return -log.(1 .+exp.(-u) .*(exp(-θ)-1)) ./θ
  elseif copula == "gumbel"
    return exp.(-u.^(1/θ))
  end
  throw(AssertionError("$(copula) not supported"))
end

"""
  copulagen(copula::String, r::Matrix{Float}, θ::Union{Float64, Int})

Auxiliary function used to generate data from archimedean copula (clayton, gumbel, frank or amh)
parametrised by a single parameter θ given a matrix of independent [0,1] distributerd
random vectors.

```jldoctest
  julia> copulagen("clayton", [0.2 0.6 0.9; 0.4 0.5 0.8], 2.)
  2×2 Array{Float64,2}:
   0.675778  0.851993
   0.687482  0.736394
```

"""
function copulagen(copula::String, r::Matrix{T}, θ::Union{Float64, Int}) where T <:AbstractFloat
  u = -log.(r[:,1:end-1])./getV0(θ, r[:,end], copula)
  phi(u, θ, copula)
end

"""

  archcopulagen(t::Int, n::Int, θ::Union{Float64, Int}, copula::String; rev::Bool = false)

Returns: t x n Matrix{Float}, t realisations of n-variate data generated from Archimedean
one parameter copula.

Following copula families are supported: clayton, frank, gumbel and amh --
Ali-Mikhail-Haq.

If rev == true, reverse the copula output i.e. u → 1-u (we call it reversed copula).
It cor == pearson, kendall, uses correlation coeficient as a parameter

```jldoctest
julia> Random.seed!(43);

julia> archcopulagen(10, 2, 1, "clayton")
10×2 Array{Float64,2}:
 0.770331  0.932834
 0.472847  0.0806845
 0.970749  0.653029
 0.622159  0.0518025
 0.402461  0.228549
 0.946375  0.842883
 0.809076  0.129038
 0.747983  0.433829
 0.374341  0.437269
 0.973066  0.910103

 ```
"""
function archcopulagen(t::Int, n::Int, θ::Union{Float64, Int}, copula::String;
                                                              rev::Bool = false,
                                                              cor::String = "")
  cor in ["Spearman", "Kendall", ""] || throw(AssertionError("$(cor) correlation not supported"))
  copula in ["clayton", "amh", "frank", "gumbel"] || throw(AssertionError("$(copula) copula is not supported"))
  if (n == 2)*((copula != "gumbel")*(θ < 0) | (copula == "amh")*(θ in [0,1]))
    return chaincopulagen(t, [θ], copula; rev=rev, cor=cor)
  elseif cor == "Spearman"
    θ = useρ(θ , copula)
  elseif cor == "Kendall"
    θ = useτ(θ , copula)
  else
    testθ(θ, copula)
  end
  u = copulagen(copula, rand(t,n+1), θ)
  rev ? 1 .-u : u
end

"""
  function testθ(θ::Union{Float64, Int}, copula::String)

Tests the parameter θ value for archimedean copula, returns void
"""
function testθ(θ::Union{Float64, Int}, copula::String)
  if copula == "gumbel"
    θ >= 1 || throw(DomainError("generaton not supported for θ < 1"))
  elseif copula == "amh"
    1 > θ > 0 || throw(DomainError("amh multiv. copula supported only for 0 < θ < 1"))
  else
    θ > 0 || throw(DomainError("generaton not supported for θ ≤ 0"))
  end
  Nothing
end

"""
  useρ(ρ::Float64, copula::String)

Tests the available pearson correlation for archimedean copula, returns Float,
corresponding copula parameter θ.
```jldoctest

julia> useρ(0.75, "gumbel")
2.294053859606698
```

"""
function useρ(ρ::Float64, copula::String)
  0 < ρ < 1 || throw(DomainError("correlation coeficiant must fulfill 0 < ρ < 1"))
  if copula == "amh"
    0 < ρ < 0.5 || throw(DomainError("correlation coeficiant must fulfill 0 < ρ < 0.5"))
  end
  ρ2θ(ρ, copula)
end

"""
  useτ(ρ::Float64, copula::String)

Tests the available kendall's correlation for archimedean copula, returns Float,
corresponding copula parameter θ.

```jldoctest

julia> useτ(0.5, "clayton")
2.0
```
"""
function useτ(τ::Float64, copula::String)
  0 < τ < 1 || throw(DomainError("correlation coeficiant must fulfill 0 < τ < 1"))
  if copula == "amh"
    0 < τ < 1/3 || throw(DomainError("correlation coeficiant must fulfill 0 < τ < 1/3"))
  end
  τ2θ(τ, copula)
end


text:
using Statistics
using PyPlot
using DataFrames
using Distributed

# Set-up for Experimental protocol (section 2.4)
#
# Sgm stands for Bayes Factor Surprise
# Ssh stands for Shannon Surprise

mutable struct SurpriseExperiment{S}
    sigma::Float64
    changeprobability::Float64
    nsteps::Int
    seed::Array{Int, 1}
    specs::S
    learnertype::Any
    learnerparam::Any
end
function SurpriseExperiment(;sigma = 0.5, changeprobability = 0.05,
                            nsteps = 1000, seed = [3, 4],
                            learnertype = GNas12,
                            Npar = 1, etaleak = 0.1, m = 0.1
                            )

    specs = SpecsGaussian(; mu_0 = 0., sigma_0 = 1.,
                    sigmas = [sigma],
                    changeprobabilities = [changeprobability],
                    nseeds = length(seed), seedenvpolicy = seed,
                    seedlearner = seed .* 99,
                    nsteps = nsteps)

    if learnertype == GParticleFilter
        learnerparam = (Npar, sigma, changeprobability)
    elseif learnertype == GLeakyIntegrator
        learnerparam = etaleak
    elseif learnertype == GVarSmile
        learnerparam = (m, sigma, changeprobability)
    elseif learnertype == GMPN
        learnerparam = (Npar, sigma, changeprobability)
    elseif learnertype == GSOR
        learnerparam = (Npar, sigma, changeprobability)
    elseif learnertype == GNas12
        learnerparam = (sigma, changeprobability)
    end
    SurpriseExperiment(sigma, changeprobability, nsteps, seed, specs,
                    learnertype, learnerparam)
end

function runSingleExperiment(exp::SurpriseExperiment)
    callback = run!(exp.specs.envtype, exp.learnertype, exp.learnerparam,
                    sigma = exp.sigma[1],
                    changeprobability = exp.changeprobability[1],
                    nsteps = exp.nsteps,
                    iseedenvpolicy = exp.specs.seedenvpolicy[1],
                    iseedlearner = exp.specs.seedlearner[1]
                    )

end


text:
using GNSS
using Base.Test
using Lint

# code cleanliness
@test isempty(lintpkg("GNSS", returnMsgs=true))

# try to load type-b file
@test_throws ErrorException readsp3(joinpath(dirname(@__FILE__), "data/igs11720.sp3"))

# load type-c file
@test size(readsp3(joinpath(dirname(@__FILE__), "data/igs18000.sp3"))) == (96, 10, 32)


text:
tests = ["basic"]

println("Running tests ...")

for t in tests
    fn = "test_$t.jl"
    println("* $fn ...")
    include(fn)
end


text:
using HermiteNormalForm
using LinearAlgebra
using Test

@testset "HermiteNormalForm.jl" begin
    for m in 1:20, n in 1:20
        A = rand(-big(3):3, m, n)
        F = hnf(A)
        if F !== nothing
            H, U = F
            @test ishnf(H)
            @test abs(det(U.//1)) == 1
            @test LinearAlgebra.issuccess(lu(Rational.(A)', check=false))
        else
            @test !LinearAlgebra.issuccess(lu(Rational.(A)', check=false))
        end
    end
end


text:
#=
Structure ConcaveProperties stores the following properties. They are needed for function VHACD for decomposing a 3D surface
into a set of near convex solids.

outputDirectory   ... Output directory
outputBasename    ... Base name for output file (a number and file extension will be added automatically)
showLog           ... show log file upon completion
resolution        ... Maximum number of voxels generated during the voxelization stage (default=100,000, range=10,000-16,000,000)
depth             ... Maximum number of clipping stages. During each split stage, solids with a concavity higher than the user defined threshold are clipped according the 'best' clipping plane (default=20, range=1-32)
concavity         ... Maximum allowed concavity (default=0.0025, range=0.0-1.0)
planeDownsampling ... Controls the granularity of the search for the 'best' clipping plane (default=4, range=1-16)
convexhullDownsampling ... Controls the precision of the convex-hull generation process during the clipping plane selection stage (default=4, range=1-16)
alpha    ... Controls the bias toward clipping along symmetry planes (default=0.05, range=0.0-1.0)
beta     ... Controls the bias toward clipping along revolution axes (default=0.05, range=0.0-1.0)
gamma    ... Controls the maximum allowed concavity during the merge stage (default=0.00125, range=0.0-1.0)
delta    ... Controls the bias toward maximaxing local concavity (default=0.05, range=0.0-1.0)
pca      ... Enable/disable normalizing the mesh before applying the convex decomposition (default=false)
mode     ... false: voxel-based approximate convex decomposition, true: tetrahedron-based approximate convex decomposition
maxNumVerticesPerCH     ... Controls the maximum number of triangles per convex-hull (default=64, range=4-1024)
minVolumePerCH          ... Controls the adaptive sampling of the generated convex-hulls (default=0.0001, range=0.0-0.01)
convexhullApproximation ... Enable/disable approximation when computing convex-hulls
=#

struct ConcaveProperties
   outputDirectory::String
   outputBasename::String
   showLog::Bool
   resolution::Int
   depth::Int
   concavity::Float64
   planeDownsampling::Int
   convexhullDownsampling::Int
   alpha::Float64
   beta::Float64
   gamma::Float64
   delta::Float64
   pca::Bool
   mode::Bool
   convexhullApproximation::Bool
   maxNumVerticesPerCH::Int
   minVolumePerCH::Float64

   function ConcaveProperties(filename::String; showLog=false, resolution=100000, depth=20,
               concavity=0.0025, planeDownsampling=4, convexhullDownsampling=4,
               alpha=0.05, beta=0.05, gamma=0.00125, delta=0.05, pca=false, mode=false, convexhullApproximation=true, maxNumVerticesPerCH=64, minVolumePerCH=0.0001)
     @assert(isfile(filename))
     @assert(filename[end-3:end] == ".obj")
     @assert(isinteger(resolution))
     @assert(resolution>=10000 && resolution<=16000000)
     @assert(isinteger(depth))
     @assert(depth>=1 && depth<=32)
     @assert(concavity>=0.0 && concavity<=1.0)
     @assert(isinteger(planeDownsampling))
     @assert(planeDownsampling>=1 && planeDownsampling<=16)
     @assert(isinteger(convexhullDownsampling))
     @assert(convexhullDownsampling>=1 && convexhullDownsampling<=16)
     @assert(alpha>=0.0 && alpha<=1.0)
     @assert(beta>=0.0 && beta<=1.0)
     @assert(gamma>=0.0 && gamma<=1.0)
     @assert(delta>=0.0 && delta<=1.0)
     @assert(isinteger(maxNumVerticesPerCH))
     @assert(maxNumVerticesPerCH>=4 && maxNumVerticesPerCH<=1024)
     @assert(minVolumePerCH>=0.0 && minVolumePerCH<=0.01)

     outputDirectory = joinpath(dirname(filename),"convexSolids_" * basename(filename))
     outputBasename = rsplit(basename(filename), ".")[1]

     new(outputDirectory, outputBasename, showLog, resolution, depth, concavity,
         planeDownsampling, convexhullDownsampling, alpha, beta, gamma, delta, pca,
         mode, convexhullApproximation, maxNumVerticesPerCH, minVolumePerCH)
end; end


text:
using Documenter, Kroki

makedocs(
  authors = "Joris Kraak",
  sitename = "Kroki.jl",
  pages = ["Home" => "index.md", "API" => "api.md"],
)

deploydocs(
  devbranch = "development",
  devurl = "latest",
  push_preview = true,
  repo = "github.com/bauglir/Kroki.jl.git",
)


text:
# Smooth ReLU
t_ss = 25.0
m_ss = tan(deg2rad(10.0)) # 10 degree slope
x_off_ss = 0.5
slope_smooth(x) = m_ss * 1.0 / t_ss .* log.(1.0 .+ exp.(t_ss .* (x[1] .- x_off_ss)))

# x = range(-1.0, stop = 4.0, length = 1000)
# plot(x, ss.(x), aspect_ratio = :equal)

slope_smooth_2D_lc = environment_2D(slope_smooth)

slope1_2D_lc = environment_2D(x -> 0.5 * x[1])


text:
module Time

import Calendar
using Calendar

import Base.string

include("business_calendars.jl")

include("day_count.jl")

end


text:
export MPIMeasurementProtocol, MPIMeasurementProtocolParams
"""
Parameters for the MPIMeasurementProtocol
"""
Base.@kwdef mutable struct MPIMeasurementProtocolParams <: ProtocolParams
  "Foreground frames to measure. Overwrites sequence frames"
  fgFrames::Int64 = 1
  "Background frames to measure. Overwrites sequence frames"
  bgFrames::Int64 = 1
  "If set the tx amplitude and phase will be set with control steps"
  controlTx::Bool = false
  "If unset no background measurement will be taken"
  measureBackground::Bool = true
  "Sequence to measure"
  sequence::Union{Sequence, Nothing} = nothing
  "Remember background measurement"
  rememberBGMeas::Bool = false
end
function MPIMeasurementProtocolParams(dict::Dict, scanner::MPIScanner)
  sequence = nothing
  if haskey(dict, "sequence")
    sequence = Sequence(scanner, dict["sequence"])
    dict["sequence"] = sequence
    delete!(dict, "sequence")
  end
  params = params_from_dict(MPIMeasurementProtocolParams, dict)
  params.sequence = sequence
  return params
end
MPIMeasurementProtocolParams(dict::Dict) = params_from_dict(MPIMeasurementProtocolParams, dict)

Base.@kwdef mutable struct MPIMeasurementProtocol <: Protocol
  name::AbstractString
  description::AbstractString
  scanner::MPIScanner
  params::MPIMeasurementProtocolParams
  biChannel::Union{BidirectionalChannel{ProtocolEvent}, Nothing} = nothing
  executeTask::Union{Task, Nothing} = nothing

  bgMeas::Array{Float32, 4} = zeros(Float32,0,0,0,0)
  done::Bool = false
  cancelled::Bool = false
  finishAcknowledged::Bool = false
  measuring::Bool = false
  txCont::Union{TxDAQController, Nothing} = nothing
end

function requiredDevices(protocol::MPIMeasurementProtocol)
  result = [AbstractDAQ]
  if protocol.params.controlTx
    push!(result, TxDAQController)
  end
  return result
end

function _init(protocol::MPIMeasurementProtocol)
  if isnothing(protocol.params.sequence)
    throw(IllegalStateException("Protocol requires a sequence"))
  end
  protocol.done = false
  protocol.cancelled = false
  protocol.finishAcknowledged = false
  if protocol.params.controlTx
    controllers = getDevices(protocol.scanner, TxDAQController)
    if length(controllers) > 1
      throw(IllegalStateException("Cannot unambiguously find a TxDAQController as the scanner has $(length(controllers)) of them"))
    end
    protocol.txCont = controllers[1]
    protocol.txCont.currTx = nothing
  else
    protocol.txCont = nothing
  end
  protocol.bgMeas = zeros(Float32,0,0,0,0)
end

function timeEstimate(protocol::MPIMeasurementProtocol)
  est = "Unknown"
  if !isnothing(protocol.params.sequence)
    params = protocol.params
    seq = params.sequence
    totalFrames = (params.fgFrames + params.bgFrames) * acqNumFrameAverages(seq)
    samplesPerFrame = rxNumSamplingPoints(seq) * acqNumAverages(seq) * acqNumPeriodsPerFrame(seq)
    totalTime = (samplesPerFrame * totalFrames) / (125e6/(txBaseFrequency(seq)/rxSamplingRate(seq)))
    time = totalTime * 1u"s"
    est = string(time)
    @show est
  end
  return est
end

function enterExecute(protocol::MPIMeasurementProtocol)
  protocol.done = false
  protocol.cancelled = false
  protocol.finishAcknowledged = false
end

function _execute(protocol::MPIMeasurementProtocol)
  @info "Measurement protocol started"

  performMeasurement(protocol)

  put!(protocol.biChannel, FinishedNotificationEvent())
  while !(protocol.finishAcknowledged)
    handleEvents(protocol)
    protocol.cancelled && throw(CancelException())
    sleep(0.05)
  end

  @info "Protocol finished."
  close(protocol.biChannel)
end

function performMeasurement(protocol::MPIMeasurementProtocol)
  if (length(protocol.bgMeas) == 0 || !protocol.params.rememberBGMeas) && protocol.params.measureBackground
    askChoices(protocol, "Press continue when background measurement can be taken", ["Continue"])
    acqNumFrames(protocol.params.sequence, protocol.params.bgFrames)
    measurement(protocol)
    protocol.bgMeas = protocol.scanner.seqMeasState.buffer
    askChoices(protocol, "Press continue when foreground measurement can be taken", ["Continue"])
  end

  acqNumFrames(protocol.params.sequence, protocol.params.fgFrames)
  measurement(protocol)
end

function measurement(protocol::MPIMeasurementProtocol)
  # Start async measurement
  protocol.measuring = true
  measState = asyncMeasurement(protocol)
  producer = measState.producer
  consumer = measState.consumer

  # Handle events
  while !istaskdone(consumer)
    handleEvents(protocol)
    protocol.cancelled && throw(CancelException())
    sleep(0.05)
  end
  protocol.measuring = false

  # Check tasks
  ex = nothing
  if Base.istaskfailed(producer)
    @error "Producer failed"
    stack = Base.catch_stack(producer)[1]
    @error stack[1]
    @error stacktrace(stack[2])
    ex = stack[1]
  end
  if Base.istaskfailed(consumer)
    @error "Consumer failed"
    stack = Base.catch_stack(consumer)[1]
    @error stack[1]
    @error stacktrace(stack[2])
    if isnothing(ex)
      ex = stack[1]
    end
  end
  if !isnothing(ex)
    throw(ErrorException("Measurement failed, see logged exceptions and stacktraces"))
  end

end

function asyncMeasurement(protocol::MPIMeasurementProtocol)
  scanner = protocol.scanner
  sequence = protocol.params.sequence
  prepareAsyncMeasurement(scanner, sequence)
  if protocol.params.controlTx
    controlTx(protocol.txCont, sequence, protocol.txCont.currTx)
  end
  scanner.seqMeasState.producer = @tspawnat scanner.generalParams.producerThreadID asyncProducer(scanner.seqMeasState.channel, scanner, sequence, prepTx = !protocol.params.controlTx)
  bind(scanner.seqMeasState.channel, scanner.seqMeasState.producer)
  scanner.seqMeasState.consumer = @tspawnat scanner.generalParams.consumerThreadID asyncConsumer(scanner.seqMeasState.channel, scanner)
  return scanner.seqMeasState
end

function cleanup(protocol::MPIMeasurementProtocol)
  # NOP
end

function stop(protocol::MPIMeasurementProtocol)
  put!(protocol.biChannel, OperationNotSupportedEvent(StopEvent()))
end

function resume(protocol::MPIMeasurementProtocol)
   put!(protocol.biChannel, OperationNotSupportedEvent(ResumeEvent()))
end

function cancel(protocol::MPIMeasurementProtocol)
  protocol.cancelled = true
  #put!(protocol.biChannel, OperationNotSupportedEvent(CancelEvent()))
  # TODO stopTx and reconnect for pipeline and so on
end

function handleEvent(protocol::MPIMeasurementProtocol, event::DataQueryEvent)
  data = nothing
  if event.message == "CURRFRAME"
    data = max(protocol.scanner.seqMeasState.nextFrame - 1, 0)
  elseif startswith(event.message, "FRAME")
    frame = tryparse(Int64, split(event.message, ":")[2])
    if !isnothing(frame) && frame > 0 && frame <= protocol.scanner.seqMeasState.numFrames
        data = protocol.scanner.seqMeasState.buffer[:, :, :, frame:frame]
    end
  elseif event.message == "BUFFER"
    data = protocol.scanner.seqMeasState.buffer
  else
    put!(protocol.biChannel, UnknownDataQueryEvent(event))
    return
  end
  put!(protocol.biChannel, DataAnswerEvent(data, event))
end

function handleEvent(protocol::MPIMeasurementProtocol, event::ProgressQueryEvent)
  reply = nothing
  if length(protocol.bgMeas) > 0 && !protocol.measuring
    reply = ProgressEvent(0, 1, "No bg meas", event)
  elseif !isnothing(protocol.scanner.seqMeasState)
    framesTotal = protocol.scanner.seqMeasState.numFrames
    framesDone = min(protocol.scanner.seqMeasState.nextFrame - 1, framesTotal)
    reply = ProgressEvent(framesDone, framesTotal, "Frames", event)
  else
    reply = ProgressEvent(0, 0, "N/A", event)
  end
  put!(protocol.biChannel, reply)
end

handleEvent(protocol::MPIMeasurementProtocol, event::FinishedAckEvent) = protocol.finishAcknowledged = true

function handleEvent(protocol::MPIMeasurementProtocol, event::DatasetStoreStorageRequestEvent)
  store = event.datastore
  scanner = protocol.scanner
  params = event.params
  data = protocol.scanner.seqMeasState.buffer
  bgdata = nothing
  if length(protocol.bgMeas) > 0
    bgdata = protocol.bgMeas
  end
  filename = saveasMDF(store, scanner, protocol.params.sequence, data, params, bgdata = bgdata)
  @show filename
  put!(protocol.biChannel, StorageSuccessEvent(filename))
end


text:
# A module that collects all necessary helper functions
module JuTools
export compute_accuracy, scale_data, scale_data!, shuffle_data, shuffle_data!, split_data
export data_generate_linear_2d

import Statistics
import Random

"""
Compute accuracy given prediction and real values
"""
function compute_accuracy(y_pred::Array, y_real::Array)::AbstractFloat
    @assert ndims(y_pred) == ndims(y_real) == 1
    @assert size(y_pred) == size(y_real)
    acc = Statistics.mean(y_pred .== y_real)
    return acc
end

"""
Scale a 2d data, to unit variance
"""
function scale_data(X::Array{T} where T<:Number)::Array
    @assert ndims(X) == 2
    u = Statistics.mean(X, dims=1)
    s = Statistics.std(X, dims=1)
    res = (X .- u) ./ s
    return res
end

"""
Scale a 2d data in place, to unit variance
"""
function scale_data!(X::Array{T} where T<:Number)
    @assert ndims(X) == 2
    u = Statistics.mean(X, dims=1)
    s = Statistics.std(X, dims=1)
    X .= (X .- u) ./ s
    return nothing
end

"""
Shuffle X and y data\\
`X` is 2 dim, `y` is 1 dim\\
Returns shuffled `X` and `y` in Tuple
"""
function shuffle_data(X::Array{T} where T<:Number, y::Array{T} where T<:Number)::Tuple{Array,Array}
    @assert ndims(X) == 2
    @assert ndims(y) == 1
    @assert size(X)[1] == size(y)[1]
    data = hcat(X, y)
    data .= data[Random.randperm(size(data)[1]), :]
    X_shuffled = data[:, 1:end-1]
    y_shuffled = data[:, end]
    return (X_shuffled, y_shuffled)
end

"""
Shuffle X and y data in place\\
`X` is 2 dim, `y` is 1 dim\\
Returns shuffled `X` and `y` in Tuple
"""
function shuffle_data!(X::Array{T} where T<:Number, y::Array{T} where T<:Number)
    @assert ndims(X) == 2
    @assert ndims(y) == 1
    @assert size(X)[1] == size(y)[1]
    data = hcat(X, y)
    data .= data[Random.randperm(size(data)[1]), :]
    X .= data[:, 1:end-1]
    y .= data[:, end]
    return nothing
end

"""
Split X and y data to training set and test set\\
`X` is 2 dim, `y` is 1 dim\\
`ratio` should be in range (0, 1), defines the test data percentage\\
Returns `X_train`, `X_test`, `y_train`, `y_test`
"""
function split_data(X::Array{T} where T<:Number, y::Array{T} where T<:Number;
        shuffle::Bool=true, ratio::AbstractFloat=0.3)::Tuple{Array,Array,Array,Array}
    @assert ndims(X) == 2
    @assert ndims(y) == 1
    @assert size(X)[1] == size(y)[1]
    @assert 0.0 < ratio < 1.0
    if shuffle
        X, y = shuffle_data(X, y)
    end
    m = size(X)[1]
    n = floor(Int, float(m) * ratio)
    # make sure at least one value in the arrays
    n = min(n, m-1)
    n = max(1, n)
    X_test = X[1:n, :]
    y_test = y[1:n]
    X_train = X[n+1:end, :]
    y_train = y[n+1:end]
    return (X_train, X_test, y_train, y_test)
end

# functions for data generation

"""
Generate 2 dimensional data for testing\\
`linear_func` should define a linear function\\
e.g. x1 -> x1 * 2 + 1\\
Returns a Tuple for X and y
"""
function data_generate_linear_2d(;linear_func::Function=k -> k * 2 + 1, data_size::Integer=1000,
        range_min::AbstractFloat=0.0, range_max::AbstractFloat=100.0, range_step::AbstractFloat=0.1,
        random_scale::AbstractFloat=20.0)::Tuple{Array, Array}
    @assert data_size > 0
    @assert range_max > range_min
    @assert range_step > 0
    @assert range_step < (range_max - range_min)
    X_data = Random.rand(range_min:range_step:range_max, (data_size, 2))
    Y_data = Array{Int64}(undef, data_size)
    for i = 1:size(X_data)[1]
        if X_data[i, 2] > linear_func(X_data[i, 1]) + Random.randn() * random_scale
            Y_data[i] = 0
        else
            Y_data[i] = 1
        end
    end
    shuffle_data!(X_data, Y_data)
    return (X_data, Y_data)
end

"""
Generate 2 dimensional data for testing\\
`pos1` and `pos2` define the cluster center of 2 targets\\
`radius1` and `radius2` define the cluster radius\\
`random_scale` defines the scale of randomness\\
Returns a Tuple for X and y, where 2 targets are balanced
"""
function data_generate_cluster_2d(;pos1::Tuple{AbstractFloat, AbstractFloat}=(10.0, 10.0), pos2::Tuple{AbstractFloat, AbstractFloat}=(30.0, 30.0),
        radius1::AbstractFloat=5.0, radius2::AbstractFloat=5.0, random_scale::AbstractFloat=10.0, data_size::Integer=1000)::Tuple{Array, Array}
    @assert radius1 > 0
    @assert radius2 > 0
    @assert data_size >= 2
    data_size_1 = div(data_size, 2)
    data_size_2 = data_size - data_size_1
    step_scale = min(radius1 / 10.0, radius2 / 10.0)
    step_scale = 10.0 ^ (trunc(Int32, step_scale) - 2)
    X_data_1 = Random.rand(-radius1:step_scale:radius1, (data_size_1, 2))
    X_data_2 = Random.rand(-radius2:step_scale:radius2, (data_size_2, 2))
    X_data_1[:, 1] .= X_data_1[:, 1] .+ pos1[1]
    X_data_1[:, 2] .= X_data_1[:, 2] .+ pos1[2]
    X_data_2[:, 1] .= X_data_2[:, 1] .+ pos2[1]
    X_data_2[:, 2] .= X_data_2[:, 2] .+ pos2[2]
    X_data_1 .= X_data_1 .+ (Random.randn(size(X_data_1)) .* random_scale)
    X_data_2 .= X_data_2 .+ (Random.randn(size(X_data_2)) .* random_scale)
    Y_data_1 = ones(data_size_1)
    Y_data_2 = zeros(data_size_2)
    X_data = cat(X_data_1, X_data_2, dims=1)
    Y_data = cat(Y_data_1, Y_data_2, dims=1)
    shuffle_data!(X_data, Y_data)
    return (X_data, Y_data)
end

end

max_stars_repo_name: teamclouday/JuliaLearn

max_issues_repo_name: teamclouday/JuliaLearn

max_forks_repo_name: teamclouday/JuliaLearn


text:
using ExcelFiles
using ExcelReaders
using IteratorInterfaceExtensions
using TableTraits
using TableTraitsUtils
using Dates
using DataValues
using DataFrames
using Test

@testset "ExcelFiles" begin

    filename = normpath(dirname(pathof(ExcelReaders)), "..", "test", "TestData.xlsx")

    efile = load(filename, "Sheet1")

    @test sprint((stream, data) -> show(stream, "text/html", data), efile) == "

Some Float64s	Some Strings	Some Bools	Mixed column	Mixed with NA	Float64 with NA	String with NA	Bool with NA	Some dates	Dates with NA	Some errors	Errors with NA	Column with NULL and then mixed

1.0	"A"	true	2.0	9.0	3.0	"FF"	#NA	2015-03-03T00:00:00	1965-04-03T00:00:00	#DIV/0!	#DIV/0!	#NA

1.5	"BB"	false	"EEEEE"	"III"	#NA	#NA	true	2015-02-04T10:14:00	1950-08-09T18:40:00	#N/A	#N/A	3.4

2.0	"CCC"	false	false	#NA	3.5	"GGG"	#NA	1988-04-09T00:00:00	19:00:00	#REF!	#NAME?	"HKEJW"

2.5	"DDDD"	true	1.5	true	4.0	"HHHH"	false	15:02:00	#NA	#NAME?	#NA	#NA

"

    @test sprint((stream, data) -> show(stream, "application/vnd.dataresource+json", data), efile) == "{\"schema\":{\"fields\":[{\"name\":\"Some Float64s\",\"type\":\"number\"},{\"name\":\"Some Strings\",\"type\":\"string\"},{\"name\":\"Some Bools\",\"type\":\"boolean\"},{\"name\":\"Mixed column\",\"type\":\"string\"},{\"name\":\"Mixed with NA\",\"type\":\"string\"},{\"name\":\"Float64 with NA\",\"type\":\"number\"},{\"name\":\"String with NA\",\"type\":\"string\"},{\"name\":\"Bool with NA\",\"type\":\"boolean\"},{\"name\":\"Some dates\",\"type\":\"string\"},{\"name\":\"Dates with NA\",\"type\":\"string\"},{\"name\":\"Some errors\",\"type\":\"string\"},{\"name\":\"Errors with NA\",\"type\":\"string\"},{\"name\":\"Column with NULL and then mixed\",\"type\":\"string\"}]},\"data\":[{\"Some Float64s\":1.0,\"Some Strings\":\"A\",\"Some Bools\":true,\"Mixed column\":2.0,\"Mixed with NA\":9.0,\"Float64 with NA\":3.0,\"String with NA\":\"FF\",\"Bool with NA\":null,\"Some dates\":\"2015-03-03T00:00:00\",\"Dates with NA\":\"1965-04-03T00:00:00\",\"Some errors\":{\"errorcode\":7},\"Errors with NA\":{\"errorcode\":7},\"Column with NULL and then mixed\":null},{\"Some Float64s\":1.5,\"Some Strings\":\"BB\",\"Some Bools\":false,\"Mixed column\":\"EEEEE\",\"Mixed with NA\":\"III\",\"Float64 with NA\":null,\"String with NA\":null,\"Bool with NA\":true,\"Some dates\":\"2015-02-04T10:14:00\",\"Dates with NA\":\"1950-08-09T18:40:00\",\"Some errors\":{\"errorcode\":42},\"Errors with NA\":{\"errorcode\":42},\"Column with NULL and then mixed\":3.4},{\"Some Float64s\":2.0,\"Some Strings\":\"CCC\",\"Some Bools\":false,\"Mixed column\":false,\"Mixed with NA\":null,\"Float64 with NA\":3.5,\"String with NA\":\"GGG\",\"Bool with NA\":null,\"Some dates\":\"1988-04-09T00:00:00\",\"Dates with NA\":\"19:00:00\",\"Some errors\":{\"errorcode\":23},\"Errors with NA\":{\"errorcode\":29},\"Column with NULL and then mixed\":\"HKEJW\"},{\"Some Float64s\":2.5,\"Some Strings\":\"DDDD\",\"Some Bools\":true,\"Mixed column\":1.5,\"Mixed with NA\":true,\"Float64 with NA\":4.0,\"String with NA\":\"HHHH\",\"Bool with NA\":false,\"Some dates\":\"15:02:00\",\"Dates with NA\":null,\"Some errors\":{\"errorcode\":29},\"Errors with NA\":null,\"Column with NULL and then mixed\":null}]}"

    @test sprint(show, efile) == "4x13 Excel file\nSome Float64s │ Some Strings │ Some Bools │ Mixed column │ Mixed with NA\n──────────────┼──────────────┼────────────┼──────────────┼──────────────\n1.0           │ A            │ true       │ 2.0          │ 9.0          \n1.5           │ BB           │ false      │ \"EEEEE\"      │ \"III\"        \n2.0           │ CCC          │ false      │ false        │ #NA          \n2.5           │ DDDD         │ true       │ 1.5          │ true         \n... with 8 more columns: Float64 with NA, String with NA, Bool with NA, Some dates, Dates with NA, Some errors, Errors with NA, Column with NULL and then mixed"

    @test TableTraits.isiterabletable(efile) == true
    @test IteratorInterfaceExtensions.isiterable(efile) == true
    @test showable("text/html", efile) == true
    @test showable("application/vnd.dataresource+json", efile) == true

    @test isiterable(efile) == true

    full_dfs = [create_columns_from_iterabletable(load(filename, "Sheet1!C3:O7")), create_columns_from_iterabletable(load(filename, "Sheet1"))]
    for (df, names) in full_dfs
        @test length(df) == 13
        @test length(df[1]) == 4

        @test df[1] == [1., 1.5, 2., 2.5]
        @test df[2] == ["A", "BB", "CCC", "DDDD"]
        @test df[3] == [true, false, false, true]
        @test df[4] == [2, "EEEEE", false, 1.5]
        @test df[5] == [9., "III", NA, true]
        @test df[6] == [3., NA, 3.5, 4]
        @test df[7] == ["FF", NA, "GGG", "HHHH"]
        @test df[8] == [NA, true, NA, false]
        @test df[9] == [Date(2015, 3, 3), DateTime(2015, 2, 4, 10, 14), Date(1988, 4, 9), Dates.Time(15, 2, 0)]
        @test df[10] == [Date(1965, 4, 3), DateTime(1950, 8, 9, 18, 40), Dates.Time(19, 0, 0), NA]
        @test eltype(df[11]) == ExcelReaders.ExcelErrorCell
        @test df[12][1][] isa ExcelReaders.ExcelErrorCell
        @test df[12][2][] isa ExcelReaders.ExcelErrorCell
        @test df[12][3][] isa ExcelReaders.ExcelErrorCell
        @test df[12][4] == NA
        @test df[13] == [NA, 3.4, "HKEJW", NA]
    end

    df, names = create_columns_from_iterabletable(load(filename, "Sheet1!C4:O7", header=false))
    @test names == [:x1,:x2,:x3,:x4,:x5,:x6,:x7,:x8,:x9,:x10,:x11,:x12,:x13]
    @test length(df[1]) == 4
    @test length(df) == 13
    @test df[1] == [1., 1.5, 2., 2.5]
    @test df[2] == ["A", "BB", "CCC", "DDDD"]
    @test df[3] == [true, false, false, true]
    @test df[4] == [2, "EEEEE", false, 1.5]
    @test df[5] == [9., "III", NA, true]
    @test df[6] == [3, NA, 3.5, 4]
    @test df[7] == ["FF", NA, "GGG", "HHHH"]
    @test df[8] == [NA, true, NA, false]
    @test df[9] == [Date(2015, 3, 3), DateTime(2015, 2, 4, 10, 14), DateTime(1988, 4, 9), Dates.Time(15, 2, 0)]
    @test df[10] == [Date(1965, 4, 3), DateTime(1950, 8, 9, 18, 40), Dates.Time(19, 0, 0), NA]
    @test isa(df[11][1], ExcelReaders.ExcelErrorCell)
    @test isa(df[11][2], ExcelReaders.ExcelErrorCell)
    @test isa(df[11][3], ExcelReaders.ExcelErrorCell)
    @test isa(df[11][4], ExcelReaders.ExcelErrorCell)
    @test isa(df[12][1][], ExcelReaders.ExcelErrorCell)
    @test isa(df[12][2][], ExcelReaders.ExcelErrorCell)
    @test isa(df[12][3][], ExcelReaders.ExcelErrorCell)
    @test DataValues.isna(df[12][4])
    @test df[13] == [NA, 3.4, "HKEJW", NA]

    good_colnames = [:c1, :c2, :c3, :c4, :c5, :c6, :c7, :c8, :c9, :c10, :c11, :c12, :c13]

    df, names = create_columns_from_iterabletable(load(filename, "Sheet1!C4:O7", header=false, colnames=good_colnames))
    @test names == good_colnames
    @test length(df[1]) == 4
    @test length(df) == 13
    @test df[1] == [1., 1.5, 2., 2.5]
    @test df[2] == ["A", "BB", "CCC", "DDDD"]
    @test df[3] == [true, false, false, true]
    @test df[4] == [2, "EEEEE", false, 1.5]
    @test df[5] == [9., "III", NA, true]
    @test df[6] == [3, NA, 3.5, 4]
    @test df[7] == ["FF", NA, "GGG", "HHHH"]
    @test df[8] == [NA, true, NA, false]
    @test df[9] == [Date(2015, 3, 3), DateTime(2015, 2, 4, 10, 14), DateTime(1988, 4, 9), Dates.Time(15, 2, 0)]
    @test df[10] == [Date(1965, 4, 3), DateTime(1950, 8, 9, 18, 40), Dates.Time(19, 0, 0), NA]
    @test isa(df[11][1], ExcelReaders.ExcelErrorCell)
    @test isa(df[11][2], ExcelReaders.ExcelErrorCell)
    @test isa(df[11][3], ExcelReaders.ExcelErrorCell)
    @test isa(df[11][4], ExcelReaders.ExcelErrorCell)
    @test isa(df[12][1][], ExcelReaders.ExcelErrorCell)
    @test isa(df[12][2][], ExcelReaders.ExcelErrorCell)
    @test isa(df[12][3][], ExcelReaders.ExcelErrorCell)
    @test DataValues.isna(df[12][4])
    @test df[13] == [NA, 3.4, "HKEJW", NA]

# Test for saving DataFrame to XLSX
    input = (Day = ["Nov. 27","Nov. 28","Nov. 29"], Highest = [78,79,75]) |> DataFrame
    file = save("file.xlsx", input)
    output = load("file.xlsx", "Sheet1") |> DataFrame
    @test input == output
    rm("file.xlsx")

# Test for saving DataFrame to XLSX with sheetname keyword
    input = (Day = ["Nov. 27","Nov. 28","Nov. 29"], Highest = [78,79,75]) |> DataFrame
    file = save("file.xlsx", input, sheetname="SheetName")
    output = load("file.xlsx", "SheetName") |> DataFrame
    @test input == output
    rm("file.xlsx")

    df, names = create_columns_from_iterabletable(load(filename, "Sheet1", colnames=good_colnames))
    @test names == good_colnames
    @test length(df[1]) == 4
    @test length(df) == 13
    @test df[1] == [1., 1.5, 2., 2.5]
    @test df[2] == ["A", "BB", "CCC", "DDDD"]
    @test df[3] == [true, false, false, true]
    @test df[4] == [2, "EEEEE", false, 1.5]
    @test df[5] == [9., "III", NA, true]
    @test df[6] == [3, NA, 3.5, 4]
    @test df[7] == ["FF", NA, "GGG", "HHHH"]
    @test df[8] == [NA, true, NA, false]
    @test df[9] == [Date(2015, 3, 3), DateTime(2015, 2, 4, 10, 14), DateTime(1988, 4, 9), Dates.Time(15, 2, 0)]
    @test df[10] == [Date(1965, 4, 3), DateTime(1950, 8, 9, 18, 40), Dates.Time(19, 0, 0), NA]
    @test isa(df[11][1], ExcelReaders.ExcelErrorCell)
    @test isa(df[11][2], ExcelReaders.ExcelErrorCell)
    @test isa(df[11][3], ExcelReaders.ExcelErrorCell)
    @test isa(df[11][4], ExcelReaders.ExcelErrorCell)
    @test isa(df[12][1][], ExcelReaders.ExcelErrorCell)
    @test isa(df[12][2][], ExcelReaders.ExcelErrorCell)
    @test isa(df[12][3][], ExcelReaders.ExcelErrorCell)
    @test DataValues.isna(df[12][4])
    @test df[13] == [NA, 3.4, "HKEJW", NA]

# Too few colnames
    @test_throws ErrorException create_columns_from_iterabletable(load(filename, "Sheet1!C4:O7", header=true, colnames=[:c1, :c2, :c3, :c4]))

# Test for constructing DataFrame with empty header cell
    data, names = create_columns_from_iterabletable(load(filename, "Sheet2!C5:E7"))
    @test names == [:Col1, :x1, :Col3]

end


text:
using Lazy

cfrac(n::Real) = let m = floor(n)
    if (n-m)==0
        @lazy Integer(m)
    else
        @lazy Integer(m):cfrac(1/(n-m))
    end
end

restore_cfrac(cf::List) = begin
    ans = Inf
    for x in reverse(cf)
        ans = x + 1/ans
    end
    ans
end

max_stars_repo_name: sclereid/collection

max_issues_repo_name: sclereid/collection

max_forks_repo_name: sclereid/collection


text:
using EcologicalNetworks
using EcologicalNetworksPlots
using Plots
using CSV
using DataFrames

##################################
##################################
df = CSV.read("WoodEtal_Append1_v2.csv")

# Separate based on WebScale
scale = [df[df.WebScale .== u,:] for u in unique(df.WebScale)]

# select for edge interaction columns only
edgelist_A = scale[3][:, [:PredTSN, :PreyTSN]]

# remove the "san" from the dataframe itself
edgelist_A.PredTSN = replace.(edgelist_A.PredTSN, "san" .=> "")

### Do the same for PreyTSN
# remove the "san" from the dataframe itself
edgelist_A.PreyTSN = replace.(edgelist_A.PreyTSN, "san" .=> "")
edgelist_A

# convert "string" to "integer"
edgelist_A.PredTSN = parse.(Int64, edgelist_A.PredTSN)
edgelist_A.PreyTSN = parse.(Int64, edgelist_A.PreyTSN)

#############################################
##Change PredTSN to newID starting from 1
# Sort by ID
sort!(edgelist_A, :PredTSN)
# Get unique IDs
unique_ids = unique(edgelist_A.PredTSN)
#create new IDs
newid_pred = collect(1:length(unique_ids))
# Create Dict where id => newid
id_dict = Dict(unique_ids[i] => newid_pred[i] for i in 1:length(unique_ids))
# Add new ids in dataframe
edgelist_A.newPredID = [id_dict[id] for id in edgelist_A.PredTSN]

### Do same for prey
## Change PreyTSN to new ID starting from end of newPredID
# Sort by ID
sort!(edgelist_A, :PreyTSN)
# Get unique IDs
unique_ids_prey = unique(edgelist_A.PreyTSN)
# Create new IDs
newid_prey = collect(length(id_dict)+1:length(id_dict)+length(unique_ids_prey))
# Create Dict
id_dict2 = Dict(unique_ids_prey[i] => newid_prey[i] for i in 1:length(unique_ids_prey))
# Add new ids in dataframe
edgelist_A.newPreyID = [id_dict2[id] for id in edgelist_A.PreyTSN]

A = edgelist_A[:, [:newPredID, :newPreyID]]

#sort
sort!(A, :newPredID)

#get unique values
unique!(A)

#add undirected interactions
y = vcat(A.newPredID, A.newPreyID)
x = vcat(A.newPreyID, A.newPredID)
A = DataFrame(hcat(y, x))

rename!(A, :x1 => :newPredID)
rename!(A, :x2 => :newPreyID)

sort!(A, :newPredID)

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

##Create Adjacency matrix
#create matrix of zeros with n = colPred, m = colPrey
M = zeros(Int16, length(unique(A.newPredID)), length(unique(A.newPreyID)))

## Fill adjacency matrix with ones where there are interactions
for i in 1:length(A.newPredID)
    M[A[i, 1], A[i, 2]] = 1
end
M

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

# Convert ones and zeros in boolean
M_bool = convert(Array{Bool}, M .== 1)

# Convert matrix in bipartite network
B = BipartiteNetwork(M_bool)

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

#create figure
graph_A = heatmap(B, c=ColorGradient([:white, :green]), legend=false, title="Matrice d'Intéractions de l'Archipel", xaxis="ProieID", yaxis="PrédateurID")
png("figures/undirected_unipartite_int_mat.png")

max_stars_repo_name: kiristern/reseaux

max_issues_repo_name: kiristern/reseaux

max_forks_repo_name: kiristern/reseaux


text:
include("iPredict.jl/src/iPredict.jl")
using .iPredict
using DataFrames
using Plots
using DelimitedFiles
using Plots.Measures

function make_hist(path, titleArg)
    Arr = readdlm(path)

    if titleArg == "3 Isoforms" && path == "figures/figure5_data/Weibull_H1_24_scnorm_counts.txt_3_isoforms_array.txt"
        #Untransformed plot
        p1 = histogram(Arr,
        title = titleArg,
        xlims = (0, 4),
        ylims = (0, 5),
        ylabel = "Probability Density",
        normalize=:pdf,
        bins = 100,
        axis = (font(150)), titlefontsize=150, yticks = 0:5:5,
        margin = 35mm,
        left_margin = 50mm,
        right_margin = 50mm)

    elseif titleArg == "Real"
        p1 = histogram(Arr,
        title = titleArg,
        xlims = (0, 4),
        ylims = (0, 5),
        normalize=:pdf,
        xlabel = "Mean No. Isoforms\nper Gene per Cell",
        bins = 100,
        axis = (font(150)), titlefontsize=150, yticks = 0:5:5,
        margin = 35mm,
        left_margin = 50mm,
        right_margin = 50mm)

    else
        #Untransformed plot
        p1 = histogram(Arr,
        title = titleArg,
        xlims = (0, 4),
        ylims = (0, 5),
        normalize=:pdf,
        bins = 100,
        axis = (font(150)), titlefontsize=150, yticks = 0:5:5,
        margin = 35mm,left_margin = 50mm,
        right_margin = 50mm)
    end
    #vline!([mean(Arr)], color = :black, linewidth = 8)
    return p1
end
#for countsMatrix in scnorm_counts
# "Weibull", "Random", "UniformObserved", "CellVariable"
Weibull_1 = make_hist("figures/figure5_data/Weibull_H9_96_scnorm_counts.txt_1_isoforms_array.txt", "1 Isoform")
Weibull_2 = make_hist("figures/figure5_data/Weibull_H9_96_scnorm_counts.txt_2_isoforms_array.txt", "2 Isoforms")
Weibull_3 = make_hist("figures/figure5_data/Weibull_H9_96_scnorm_counts.txt_3_isoforms_array.txt", "3 Isoforms")
Weibull_4 = make_hist("figures/figure5_data/Weibull_H9_96_scnorm_counts.txt_4_isoforms_array.txt", "4 Isoforms")
Weibull_real = make_hist("figures/figure5_data/Weibull_H9_96_scnorm_counts.txt_real_isoforms_array.txt", "Real")

Random_1 = make_hist("figures/figure5_data/Random_H9_96_scnorm_counts.txt_1_isoforms_array.txt", "1 Isoform")
Random_2 = make_hist("figures/figure5_data/Random_H9_96_scnorm_counts.txt_2_isoforms_array.txt", "2 Isoforms")
Random_3 = make_hist("figures/figure5_data/Random_H9_96_scnorm_counts.txt_3_isoforms_array.txt", "3 Isoforms")
Random_4 = make_hist("figures/figure5_data/Random_H9_96_scnorm_counts.txt_4_isoforms_array.txt", "4 Isoforms")
Random_real = make_hist("figures/figure5_data/Random_H9_96_scnorm_counts.txt_real_isoforms_array.txt", "Real")

UniformObserved_1 = make_hist("figures/figure5_data/UniformObserved_H9_96_scnorm_counts.txt_1_isoforms_array.txt", "1 Isoform")
UniformObserved_2 = make_hist("figures/figure5_data/UniformObserved_H9_96_scnorm_counts.txt_2_isoforms_array.txt", "2 Isoforms")
UniformObserved_3 = make_hist("figures/figure5_data/UniformObserved_H9_96_scnorm_counts.txt_3_isoforms_array.txt", "3 Isoforms")
UniformObserved_4 = make_hist("figures/figure5_data/UniformObserved_H9_96_scnorm_counts.txt_4_isoforms_array.txt", "4 Isoforms")
UniformObserved_real = make_hist("figures/figure5_data/UniformObserved_H9_96_scnorm_counts.txt_real_isoforms_array.txt", "Real")

CellVariable_1 = make_hist("figures/figure5_data/CellVariable_H9_96_scnorm_counts.txt_1_isoforms_array.txt", "1 Isoform")
CellVariable_2 = make_hist("figures/figure5_data/CellVariable_H9_96_scnorm_counts.txt_2_isoforms_array.txt", "2 Isoforms")
CellVariable_3 = make_hist("figures/figure5_data/CellVariable_H9_96_scnorm_counts.txt_3_isoforms_array.txt", "3 Isoforms")
CellVariable_4 = make_hist("figures/figure5_data/CellVariable_H9_96_scnorm_counts.txt_4_isoforms_array.txt", "4 Isoforms")
CellVariable_real = make_hist("figures/figure5_data/CellVariable_H9_96_scnorm_counts.txt_real_isoforms_array.txt", "Real")

p3 = plot( Weibull_1, Random_1, UniformObserved_1, CellVariable_1,
       Weibull_2, Random_2,UniformObserved_2, CellVariable_2,
       Weibull_3, Random_3,UniformObserved_3, CellVariable_3,
       Weibull_4, Random_4,UniformObserved_4, CellVariable_4,
       Weibull_real, Random_real,UniformObserved_real, CellVariable_real,
       layout = (5,4), legend = false, size = (750 *8, 750*8), axis = (font(90)), titlefontsize=90,
       margin = 45mm)

savefig(p3, "figures/supplementary_figs/figure4_sup.png")

###############################################################
using JLD2
using DelimitedFiles
using HypothesisTests
using DataFrames

Weibull_1 = readdlm("figures/figure5_data/Weibull_H9_96_scnorm_counts.txt_1_isoforms_array.txt")
Weibull_2 = readdlm("figures/figure5_data/Weibull_H9_96_scnorm_counts.txt_2_isoforms_array.txt")
Weibull_3 = readdlm("figures/figure5_data/Weibull_H9_96_scnorm_counts.txt_3_isoforms_array.txt")
Weibull_4 = readdlm("figures/figure5_data/Weibull_H9_96_scnorm_counts.txt_4_isoforms_array.txt")
Weibull_real = readdlm("figures/figure5_data/Weibull_H9_96_scnorm_counts.txt_real_isoforms_array.txt")

Random_1 = readdlm("figures/figure5_data/Random_H9_96_scnorm_counts.txt_1_isoforms_array.txt")
Random_2 = readdlm("figures/figure5_data/Random_H9_96_scnorm_counts.txt_2_isoforms_array.txt")
Random_3 = readdlm("figures/figure5_data/Random_H9_96_scnorm_counts.txt_3_isoforms_array.txt")
Random_4 = readdlm("figures/figure5_data/Random_H9_96_scnorm_counts.txt_4_isoforms_array.txt")
Random_real = readdlm("figures/figure5_data/Random_H9_96_scnorm_counts.txt_real_isoforms_array.txt")

UniformObserved_1 = readdlm("figures/figure5_data/UniformObserved_H9_96_scnorm_counts.txt_1_isoforms_array.txt")
UniformObserved_2 = readdlm("figures/figure5_data/UniformObserved_H9_96_scnorm_counts.txt_2_isoforms_array.txt")
UniformObserved_3 = readdlm("figures/figure5_data/UniformObserved_H9_96_scnorm_counts.txt_3_isoforms_array.txt")
UniformObserved_4 = readdlm("figures/figure5_data/UniformObserved_H9_96_scnorm_counts.txt_4_isoforms_array.txt")
UniformObserved_real = readdlm("figures/figure5_data/UniformObserved_H9_96_scnorm_counts.txt_real_isoforms_array.txt")

CellVariable_1 = readdlm("figures/figure5_data/CellVariable_H9_96_scnorm_counts.txt_1_isoforms_array.txt")
CellVariable_2 = readdlm("figures/figure5_data/CellVariable_H9_96_scnorm_counts.txt_2_isoforms_array.txt")
CellVariable_3 = readdlm("figures/figure5_data/CellVariable_H9_96_scnorm_counts.txt_3_isoforms_array.txt")
CellVariable_4 = readdlm("figures/figure5_data/CellVariable_H9_96_scnorm_counts.txt_4_isoforms_array.txt")
CellVariable_real = readdlm("figures/figure5_data/CellVariable_H9_96_scnorm_counts.txt_real_isoforms_array.txt")

df = DataFrame(IsoformsSimulated=Int64[], Comparison = String[], pValue = Float64[])

append!(df, DataFrame(IsoformsSimulated = [1], Comparison = ["All"], pValue=[pvalue(KSampleADTest(vec(Weibull_1), vec(Random_1), vec(UniformObserved_1), vec(CellVariable_1)))]))
append!(df, DataFrame(IsoformsSimulated = [2], Comparison = "All", pValue=[pvalue(KSampleADTest(vec(Weibull_2), vec(Random_2), vec(UniformObserved_2), vec(CellVariable_2)))]))
append!(df, DataFrame(IsoformsSimulated = [3], Comparison = "All", pValue=[pvalue(KSampleADTest(vec(Weibull_3), vec(Random_3), vec(UniformObserved_3), vec(CellVariable_3)))]))
append!(df, DataFrame(IsoformsSimulated = [4], Comparison = "All", pValue=[pvalue(KSampleADTest(vec(Weibull_4), vec(Random_4), vec(UniformObserved_4), vec(CellVariable_4)))]))

append!(df, DataFrame(IsoformsSimulated = [1], Comparison = "CellVarUnif", pValue=[pvalue(KSampleADTest(vec(UniformObserved_1), vec(CellVariable_1)))]))
append!(df, DataFrame(IsoformsSimulated = [2], Comparison = "CellVarUnif", pValue=[pvalue(KSampleADTest( vec(UniformObserved_2), vec(CellVariable_2)))]))
append!(df, DataFrame(IsoformsSimulated = [3], Comparison = "CellVarUnif", pValue=[pvalue(KSampleADTest(vec(UniformObserved_3), vec(CellVariable_3)))]))
append!(df, DataFrame(IsoformsSimulated = [4], Comparison = "CellVarUnif", pValue=[pvalue(KSampleADTest(vec(UniformObserved_4), vec(CellVariable_4)))]))

jldopen("figures/figure4_data/Supfig4_AD_test_results.jld2", "w") do file
    file["df"] = df
end
#####################
#overlap

function make_hist(path, titleArg)

     Arr = readdlm(path)
     p1 = histogram(Arr, title = titleArg, normalize=:pdf, bins = 100,
     ylims = (0,4), color = palette(:default)[2],
     xlims = (0,1),
     axis = (font(150)), titlefontsize=150, yticks = 0:2:4,
     xticks = 0:0.5:1,
     margin = 35mm,left_margin = 50mm,
     right_margin = 50mm)
     return p1

end

Weibull_1 = make_hist("figures/figure5_data/Weibull_H9_96_scnorm_counts.txt_1_isoforms_overlap_array.txt", "1 Isoform")
Weibull_2 = make_hist("figures/figure5_data/Weibull_H9_96_scnorm_counts.txt_2_isoforms_overlap_array.txt", "2 Isoforms")
Weibull_3 = make_hist("figures/figure5_data/Weibull_H9_96_scnorm_counts.txt_3_isoforms_overlap_array.txt", "3 Isoforms")
Weibull_4 = make_hist("figures/figure5_data/Weibull_H9_96_scnorm_counts.txt_4_isoforms_overlap_array.txt", "4 Isoforms")

Random_1 = make_hist("figures/figure5_data/Random_H9_96_scnorm_counts.txt_1_isoforms_overlap_array.txt", "1 Isoform")
Random_2 = make_hist("figures/figure5_data/Random_H9_96_scnorm_counts.txt_2_isoforms_overlap_array.txt", "2 Isoforms")
Random_3 = make_hist("figures/figure5_data/Random_H9_96_scnorm_counts.txt_3_isoforms_overlap_array.txt", "3 Isoforms")
Random_4 = make_hist("figures/figure5_data/Random_H9_96_scnorm_counts.txt_4_isoforms_overlap_array.txt", "4 Isoforms")

UniformObserved_1 = make_hist("figures/figure5_data/UniformObserved_H9_96_scnorm_counts.txt_1_isoforms_overlap_array.txt", "1 Isoform")
UniformObserved_2 = make_hist("figures/figure5_data/UniformObserved_H9_96_scnorm_counts.txt_2_isoforms_overlap_array.txt", "2 Isoforms")
UniformObserved_3 = make_hist("figures/figure5_data/UniformObserved_H9_96_scnorm_counts.txt_3_isoforms_overlap_array.txt", "3 Isoforms")
UniformObserved_4 = make_hist("figures/figure5_data/UniformObserved_H9_96_scnorm_counts.txt_4_isoforms_overlap_array.txt", "4 Isoforms")

CellVariable_1 = make_hist("figures/figure5_data/CellVariable_H9_96_scnorm_counts.txt_1_isoforms_overlap_array.txt", "1 Isoform")
CellVariable_2 = make_hist("figures/figure5_data/CellVariable_H9_96_scnorm_counts.txt_2_isoforms_overlap_array.txt", "2 Isoforms")
CellVariable_3 = make_hist("figures/figure5_data/CellVariable_H9_96_scnorm_counts.txt_3_isoforms_overlap_array.txt", "3 Isoforms")
CellVariable_4 = make_hist("figures/figure5_data/CellVariable_H9_96_scnorm_counts.txt_4_isoforms_overlap_array.txt", "4 Isoforms")

p3 = plot( Weibull_1, Random_1, UniformObserved_1, CellVariable_1,
       Weibull_2, Random_2,UniformObserved_2, CellVariable_2,
       Weibull_3, Random_3,UniformObserved_3, CellVariable_3,
       Weibull_4, Random_4,UniformObserved_4, CellVariable_4,
       layout = (4,4), legend = false, size = (750 *8, 750*8), axis = (font(90)), titlefontsize=90,
       margin = 45mm)

savefig(p3, "figures/supplementary_figs/figure4_sup_overlap.png")

max_stars_repo_name: jenni-westoby/Obstacles

max_issues_repo_name: jenni-westoby/Obstacles

max_forks_repo_name: jenni-westoby/Obstacles


text:
###
# Reflection tooling
##

using Base.Meta
using .Compiler: widenconst, argextype, Const, MethodMatchInfo,
    UnionSplitApplyCallInfo, UnionSplitInfo, ConstCallInfo,
    MethodResultPure, ApplyCallInfo

const sptypes_from_meth_instance = Core.Compiler.sptypes_from_meth_instance
const may_invoke_generator = Base.may_invoke_generator
function code_for_method(method, metharg, methsp, world, preexisting=false)
    @static if VERSION ≥ v"1.8.0-DEV.369"
        # https://github.com/JuliaLang/julia/pull/41920
        Core.Compiler.specialize_method(method, metharg, methsp; preexisting)
    else
        Core.Compiler.specialize_method(method, metharg, methsp, preexisting)
    end
end

transform(::Val, callsite) = callsite
function transform(::Val{:CuFunction}, callsite, callexpr, CI, mi, slottypes; params=current_params())
    sptypes = sptypes_from_meth_instance(mi)
    tt = argextype(callexpr.args[4], CI, sptypes, slottypes)
    ft = argextype(callexpr.args[3], CI, sptypes, slottypes)
    isa(tt, Const) || return callsite
    return Callsite(callsite.id, CuCallInfo(callinfo(Tuple{widenconst(ft), tt.val.parameters...}, Nothing, params)), callsite.head)
end

function find_callsites(interp::CthulhuInterpreter, CI::Union{Core.CodeInfo, IRCode},
                        stmt_info::Union{Vector, Nothing}, mi::Core.MethodInstance,
                        slottypes::Vector{Any}, optimize::Bool=true;
                        params=current_params())
    sptypes = sptypes_from_meth_instance(mi)
    callsites = Callsite[]

    for (id, c) in enumerate(isa(CI, IRCode) ? CI.stmts.inst : CI.code)
        callsite = nothing
        isa(c, Expr) || continue
        if stmt_info !== nothing && is_call_expr(c, optimize)
            info = stmt_info[id]
            if info !== nothing
                rt = argextype(SSAValue(id), CI, sptypes, slottypes)
                # in unoptimized IR, there may be `slot = rhs` expressions, which `argextype` doesn't handle
                # so extract rhs for such an case
                args = c.args
                if !optimize
                    args = (ignorelhs(c)::Expr).args
                end
                types = mapany(function (@nospecialize(arg),)
                                   t = argextype(arg, CI, sptypes, slottypes)
                                   return widenconst(ignorelimited(t))
                               end, args)
                was_return_type = false
                if isa(info, Core.Compiler.ReturnTypeCallInfo)
                    info = info.info
                    was_return_type = true
                end

                callinfos = process_info(interp, info, types, rt, optimize)
                if !was_return_type && isempty(callinfos)
                    continue
                end
                callsite = let
                    if length(callinfos) == 1
                        callinfo = callinfos[1]
                    else
                        callinfo = MultiCallInfo(Compiler.argtypes_to_type(types), rt, callinfos)
                    end
                    if was_return_type
                        callinfo = ReturnTypeCallInfo(callinfo)
                    end
                    Callsite(id, callinfo, c.head)
                end
            end
        end

        if callsite === nothing && c isa Expr
            c = ignorelhs(c)
            (c isa Expr) || continue
            if c.head === :invoke
                rt = argextype(SSAValue(id), CI, sptypes, slottypes)
                at = argextype(c.args[2], CI, sptypes, slottypes)
                callsite = Callsite(id, MICallInfo(c.args[1], rt), c.head)
            elseif c.head === :foreigncall
                # special handling of jl_new_task
                length(c.args) > 0 || continue
                if c.args[1] isa QuoteNode
                    cfunc = c.args[1].value
                    if cfunc === :jl_new_task
                        func = c.args[6]
                        ftype = widenconst(argextype(func, CI, sptypes, slottypes))
                        sig = Tuple{ftype}
                        callsite = Callsite(id, TaskCallInfo(callinfo(sig, nothing; params)), c.head)
                    end
                end
            end
        end

        if callsite !== nothing
            info = callsite.info
            if info isa MICallInfo
                mi = get_mi(info)
                meth = mi.def
                if isa(meth, Method) && nameof(meth.module) === :CUDAnative && meth.name === :cufunction
                    callsite = transform(Val(:CuFunction), callsite, c, CI, mi, slottypes; params)
                end
            end

            push!(callsites, callsite)
        end
    end
    return callsites
end

function process_info(interp, @nospecialize(info), types, @nospecialize(rt), optimize::Bool)
    is_cached(@nospecialize(key)) = haskey(optimize ? interp.opt : interp.unopt, key)
    process_recursive(@nospecialize(newinfo)) = process_info(interp, newinfo, types, rt, optimize)

    if isa(info, MethodMatchInfo)
        if info.results === missing
            return []
        end

        matches = info.results.matches
        return mapany(matches) do match::Core.MethodMatch
            mi = Core.Compiler.specialize_method(match)
            mici = MICallInfo(mi, rt)
            return is_cached(mi) ? mici : UncachedCallInfo(mici)
        end
    elseif isa(info, MethodResultPure)
        return Any[PureCallInfo(types, rt)]
    elseif isa(info, UnionSplitInfo)
        return mapreduce(process_recursive, vcat, info.matches; init=[])::Vector{Any}
    elseif isa(info, UnionSplitApplyCallInfo)
        return mapreduce(process_recursive, vcat, info.infos; init=[])::Vector{Any}
    elseif isa(info, ApplyCallInfo)
        # XXX: This could probably use its own info. For now,
        # we ignore any implicit iterate calls.
        return process_recursive(info.call)
    elseif isa(info, ConstCallInfo)
        infos = process_recursive(info.call)
        @assert length(infos) == length(info.results)
        return mapany(enumerate(info.results)) do (i, result)
            if isnothing(result)
                infos[i]
            else
                linfo = result.linfo
                mici = MICallInfo(linfo, rt)
                ConstPropCallInfo(is_cached(optimize ? linfo : result) ? mici : UncachedCallInfo(mici), result)
            end
        end
    elseif (@static isdefined(Compiler, :InvokeCallInfo) && true) && isa(info, Compiler.InvokeCallInfo)
        return Any[InvokeCallInfo(Core.Compiler.specialize_method(info.match), rt)]
    elseif (@static isdefined(Compiler, :OpaqueClosureCallInfo) && true) && isa(info, Compiler.OpaqueClosureCallInfo)
        return Any[OCCallInfo(Core.Compiler.specialize_method(info.match), rt)]
    elseif (@static isdefined(Compiler, :OpaqueClosureCreateInfo) && true) && isa(info, Compiler.OpaqueClosureCreateInfo)
        # TODO: Add ability to descend into OCs at creation site
        return []
    elseif info == false
        return []
    else
        @show CI
        @show c
        @show info
        error()
    end
end

ignorelhs(@nospecialize(x)) = isexpr(x, :(=)) ? last(x.args) : x
function is_call_expr(x::Expr, optimize::Bool)
    optimize && isexpr(x, :invoke) && return true
    return isexpr(ignorelhs(x), :call)
end

function dce!(ci, mi)
    if VERSION >= v"1.7.0-DEV.705"
        argtypes = Core.Compiler.matching_cache_argtypes(mi, nothing, false)[1]
    else
        argtypes = Core.Compiler.matching_cache_argtypes(mi, nothing)[1]
    end
    ir = Compiler.inflate_ir(ci, sptypes_from_meth_instance(mi),
                             argtypes)
    ir = dce!(ir, mi)
    Core.Compiler.replace_code_newstyle!(ci, ir, length(argtypes)-1)
end

function dce!(ir::IRCode, mi)
    compact = Core.Compiler.IncrementalCompact(ir, true)
    # Just run through the iterator without any processing
    Core.Compiler.foreach(x -> nothing, compact)
    ir = Core.Compiler.finish(compact)
end

function preprocess_ci!(ci::CodeInfo, mi, optimize, config::CthulhuConfig)
    if optimize && config.dead_code_elimination
        # if the optimizer hasn't run, the IR hasn't been converted
        # to SSA form yet and dce is not legal
        dce!(ci, mi)
        dce!(ci, mi)
    end
    return ci
end

function preprocess_ci!(ir::IRCode, mi, optimize, config::CthulhuConfig)
    if optimize && config.dead_code_elimination
        # if the optimizer hasn't run, the IR hasn't been converted
        # to SSA form yet and dce is not legal
        ir = dce!(ir, mi)
        ir = dce!(ir, mi)
    end
    return ir
end

const CompilerParams = Core.Compiler.NativeInterpreter
current_params() = CompilerParams()

function callinfo(sig, rt, max_methods=-1; params=current_params())
    methds = Base._methods_by_ftype(sig, max_methods, params.world)
    methds isa Bool && return FailedCallInfo(sig, rt)
    length(methds) < 1 && return FailedCallInfo(sig, rt)
    callinfos = CallInfo[]
    for x in methds
        meth = x[3]
        atypes = x[1]
        sparams = x[2]
        if isdefined(meth, :generator) && !may_invoke_generator(meth, atypes, sparams)
            push!(callinfos, GeneratedCallInfo(sig, rt))
        else
            mi = code_for_method(meth, atypes, sparams, params.world)
            if mi !== nothing
                push!(callinfos, MICallInfo(mi, rt))
            else
                push!(callinfos, FailedCallInfo(sig, rt))
            end
        end
    end

    @assert length(callinfos) != 0
    length(callinfos) == 1 && return first(callinfos)
    return MultiCallInfo(sig, rt, callinfos)
end

function find_caller_of(callee::MethodInstance, caller::MethodInstance)
    interp = CthulhuInterpreter()
    do_typeinf!(interp, caller)
    params = current_params()
    locs = Tuple{Core.LineInfoNode,Int}[]
    for optimize in (true, false)
        (CI, rt, infos, slottypes) = lookup(interp, caller, optimize)
        CI = preprocess_ci!(CI, caller, optimize, CONFIG)
        callsites = find_callsites(interp, CI, infos, caller, slottypes, optimize; params=params)
        callsites = filter(cs->is_callsite(cs, callee), callsites)
        foreach(cs -> add_sourceline!(locs, CI, cs.id), callsites)
    end
    # Consolidate by method, but preserve the order
    prlookup = Dict{Tuple{Symbol,Symbol},Int}()
    ulocs = Pair{Tuple{Symbol,Symbol,Int},Vector{Int}}[]
    if !isempty(locs)
        for (loc, depth) in locs
            idx = get(prlookup, (loc.method, loc.file), nothing)
            if idx === nothing
                push!(ulocs, (loc.method, loc.file, depth) => Int[])
                prlookup[(loc.method, loc.file)] = idx = length(ulocs)
            end
            lines = ulocs[idx][2]
            line = loc.line
            if line ∉ lines
                push!(lines, line)
            end
        end
    end
    return ulocs
end

function add_sourceline!(locs, CI, stmtidx::Int)
    if isa(CI, IRCode)
        stack = Base.IRShow.compute_loc_stack(CI.linetable, CI.stmts.line[stmtidx])
        for (i, idx) in enumerate(stack)
            line = CI.linetable[idx]
            line.line == 0 && continue
            push!(locs, (CI.linetable[idx], i-1))
        end
    else
        push!(locs, (CI.linetable[CI.codelocs[stmtidx]], 0))
    end
    return locs
end


text:
# reduction.jl: Array->Scalar and Array->Vector reductions.

# The entry format is (cudaname, julianame, merge, item, init)
# ai is the accumulator, xi is the array element
reduction_ops = [
("sum","sum","ai+xi","xi","0"),
("prod","prod","ai*xi","xi","1"),
("maximum","maximum","(ai>xi?ai:xi)","xi","(-INFINITY)"),
("minimum","minimum","(ai<xi?ai:xi)","xi","INFINITY"),
("sumabs","sumabs","ai+xi","(xi<0?-xi:xi)","0"),
("sumabs2","sumabs2","ai+xi","(xi*xi)","0"),
("countnz","countnz","ai+xi","(xi!=0)","0"),
]

function reduction_op(f, j=f, o...)
    J=Symbol(j)
    if isdefined(Base, J); eval(Expr(:import,:Base,J)); end
    for S in (32,64)
        T = Symbol("Float$S")
        F20 = "$(f)_$(S)_20"
        F21 = "$(f)_$(S)_21"
        @eval begin
            # Array->Scalar reduction:
            function $J(x::KnetArray{$T})
                ccall(($F20,$libknet8),$T,(Cint,Ptr{$T}),length(x),x)
            end
            # Array->Vector reduction:
            function $J(x::KnetArray{$T}, region)
                rdims = Base.reduced_dims(size(x), region)
                vdims = count(x->x>1,rdims)
                if vdims == 0   # falls back to Array->Scalar reduction
                    return fill!(similar(x,rdims), $J(x))
                elseif vdims == 1
                    i0 = 0
                    ysize = ntuple(ndims(x)) do i
                        if in(i,region)
                            1
                        elseif i0>0
                            error("Bad region $region")
                        else
                            i0=i
                            size(x,i)
                        end
                    end
                    y = similar(x, ysize)
                    nx = length(x); ny = length(y); sy = stride(x,i0)
                    ccall(($F21,$libknet8),Void,(Cint,Ptr{$T},Cint,Cint,Ptr{$T}),nx,x,sy,ny,y)
                    return y
                else
                    error("Only scalar and vector reductions supported: $((size(x),region))")
                end
            end
        end
    end
end

for f in reduction_ops
    if !isa(f,Tuple); f=(f,); end
    reduction_op(f...)
end

@zerograd countnz(a,d...)

# Norm primitives:

import Base.LinAlg: norm, vecnorm

norm(x::KnetVector, p::Real=2) = vecnorm(x, p)

function vecnorm{T}(x::KnetArray{T}, p::Real=2)
    if length(x) == 0
        zero(T)
    elseif p == 2
        sqrt(sumabs2(x))
    elseif p == 1
        sumabs(x)
    elseif p == Inf
        maximum(abs(x))
    elseif p == 0
        countnz(x)
    elseif p == -Inf
        minimum(abs(x))
    else
        sum(abs(x).^p)^(1/p)
    end
end

"""

    logsumexp(x,[dims])

Compute `log(sum(exp(x),dims))` in a numerically stable manner.

`dims` is an optional argument, if not specified the summation is over
the whole `x`, otherwise the summation is performed over the given
dimensions.  In particular if `x` is a matrix, `dims=1` sums columns
of `x` and `dims=2` sums rows of `x`.

"""
function logsumexp(x,d...)
    xmax = maximum(x,d...)
    xmax + log(sum(exp(x .- xmax),d...))
end

@primitive logsumexp(x,d...),dy,y  (dy .* exp(x .- y))

# # The xentloss interface is no good because of double normalization.

# """
# xentloss(x, p [,dims])

# Compute cross-entropy loss for unnormalized log probability estimates
# x and normalized probabilities p normalizing over the given
# dimensions.  By default normalization is over the whole array, dims=1
# normalizes over the columns, dims=2 normalizes over the rows of a 2-D
# array.

# """
# function xentloss(x,p,d...)
#     x = x .- maximum(x,d...)
#     z = log(sum(exp(x),d...))
#     return sum(p .* (x .- z)) / length(z)
# end

# function xentback(x,p,d...)
#     x = x .- maximum(x,d...)
#     x = exp(x)
#     z = sum(x,d...)
#     return x./z - p
# end

# @primitive xentloss(x,p,d...),dy,y  (dy.*xentback(x,p,d...))


text:
include("Density.jl")

function calc_design_matrix(umbrella_centers, data, sigma)
    design_matrix = zeros(Float64, length(umbrella_centers) * length(data[1]), length(umbrella_centers))
    for i = 1:length(umbrella_centers)
        for j = 1:length(data[i])
            for k = 1:length(umbrella_centers)
                diff = get_gaussian(data[i][j], umbrella_centers[k], sigma) - get_gaussian(umbrella_centers[i], umbrella_centers[k], sigma)
                design_matrix[j + length(data[i]) * (i - 1), k] = diff
            end
        end
    end

    return design_matrix
end


text:
using Test

using Beamforming
using Seis

@testset "Types" begin
    n = 5
    sx = -2:0.1:3
    sy = -1:0.1:2
    x = 1:n
    y = 0:(n - 1)
    lon = rand()
    lat = rand()
    freqs = 1:0.2:2
    @testset "$T" for T in (Float32, Float64)
        power = zeros(T, length(sx), length(sy))
        @testset "==, isequal, hash" begin
            @testset "BeamformGrid" begin
                fk = Beamforming.BeamformGrid{T, Trace{T, Vector{T}, Seis.Geographic{T}}}(
                    lon, lat, x, y, sx, sy, power, 1, [1], [2], [1], [120], [],
                    :linear, nothing, -10, 20, :normal, [], nothing
                    )
                fk′ = deepcopy(fk)
                @test fk == fk′
                @test fk !== fk′
                @test isequal(fk, fk′)
                @test hash(fk) == hash(fk′)
                @test unique([fk, fk′]) == [fk]
                # Modify
                fk′.power[1] = 1
                @test fk != fk′
                @test !isequal(fk, fk′)
                @test unique([fk, fk′]) == [fk, fk′]
            end

            @testset "$ARF" for ARF in (Beamforming.ArrayResponse, Beamforming.ArrayResponseCrosscorrelation)
                arf = ARF{T}(freqs, x, y, sx, sy, power)
                arf′ = deepcopy(arf)
                @test arf == arf′
                @test arf !== arf′
                @test isequal(arf, arf′)
                @test hash(arf) == hash(arf′)
                @test unique([arf, arf′]) == [arf]
                # Modify
                arf′.power[1] = 1
                @test arf != arf′
                @test !isequal(arf, arf′)
                @test unique([arf, arf′]) == [arf, arf′]
            end

            @testset "VespaGrid" begin
                times = freqs
                slows = sx
                grid = zeros(T, length(times), length(slows))
                vg = Beamforming.VespaGrid{T, Trace{T, Vector{T}, Seis.Geographic{T}}}(
                    lon, lat, x, y, times, slows,
                    grid, 1, [1], [1], [], :nthroot, 3, :plane, false)
                vg′ = deepcopy(vg)
                @test vg == vg′
                @test vg !== vg′
                @test isequal(vg, vg′)
                @test hash(vg) == hash(vg′)
                @test unique([vg, vg′]) == [vg]
                # Modify
                vg.power[1] = 1
                @test vg != vg′
                @test !isequal(vg, vg′)
                @test unique([vg, vg′]) == [vg, vg′]
            end
        end
    end
end

max_stars_repo_name: anowacki/Beamforming

max_issues_repo_name: anowacki/Beamforming

max_forks_repo_name: anowacki/Beamforming