text:
abstract type DiscreteSource end

struct PulseTrain <: DiscreteSource
    amplitude
    period
    duty
end

struct Step <: DiscreteSource
    amplitude
    steptime
end

struct ContinuousSource
    func::Function
end

Source = Union{DiscreteSource, ContinuousSource}

max_stars_repo_name: jonniedie/ComposableDiffEqs

max_issues_repo_name: jonniedie/ComposableDiffEqs

max_forks_repo_name: jonniedie/ComposableDiffEqs


text:
println("**********Testing Clustering**********")

facts("Random init provides k initial centroids") do
  a = FeatureVector(Dict(["hello"=>4, "word"=>0]))
  b = FeatureVector(Dict(["hello"=>0, "word"=>5]))
  c = FeatureVector(Dict(["hello"=>8, "word"=>4]))
  d = FeatureVector(Dict(["hello"=>4, "word"=>8]))
  e = FeatureVector(Dict(["kim"=>4, "julia"=>8]))
  f = FeatureVector(Dict(["jello"=>4, "kim"=>0]))
  dict = Dict([1=>a,2=>b,3=>c,4=>d,5=>e,6=>f])
  features = collect(Base.values(dict))
  centroids = random_init(features,2)

  @fact Base.length(centroids) => 2
  @fact centroids[1] != centroids[2] => true
end

facts("Max-min init provides k initial distinct centroids") do
  a = FeatureVector(Dict(["hello"=>4, "word"=>0]))
  b = FeatureVector(Dict(["hello"=>0, "word"=>5]))
  c = FeatureVector(Dict(["hello"=>8, "word"=>4]))
  d = FeatureVector(Dict(["hello"=>4, "word"=>8]))
  e = FeatureVector(Dict(["kim"=>4, "julia"=>8]))
  f = FeatureVector(Dict(["jello"=>4, "kim"=>0]))
  dict = Dict([1=>a,2=>b,3=>c,4=>d,5=>e,6=>f])
  features = collect(Base.values(dict))
  centroids = max_min_init(features,3,dist_cos)

  @fact Base.length(centroids) => 3
  @fact centroids[1] != centroids[2] => true
  @fact centroids[1] != centroids[3] => true
  @fact centroids[2] != centroids[3] => true
end

facts("Max-min init provides a random centroid and its furthest centroid") do
  a = FeatureVector(Dict(["hello"=>4, "word"=>0]))
  b = FeatureVector(Dict(["hello"=>0, "word"=>5]))
  c = FeatureVector(Dict(["hello"=>8, "word"=>4]))
  d = FeatureVector(Dict(["hello"=>4, "word"=>8]))
  e = FeatureVector(Dict(["kim"=>4, "julia"=>8]))
  f = FeatureVector(Dict(["jello"=>4, "kim"=>0]))
  dict = Dict([1=>a,2=>b,3=>c,4=>d,5=>e,6=>f])
  features = collect(Base.values(dict))
  centroids = max_min_init(features,2,dist_cos)

  for fv in features
    @fact dist_cos(centroids[1], centroids[2]) >= dist_cos(centroids[1],fv) => true
  end
end

facts("kmeans returns correct amount of clusters") do
  a = FeatureVector(Dict(["hello"=>4, "word"=>0]))
  b = FeatureVector(Dict(["hello"=>0, "word"=>5]))
  c = FeatureVector(Dict(["kim"=>4, "julia"=>8]))
  dict = Dict([1=>a,2=>b,3=>c])
  clusters = kmeans(dict,[],3)

  @fact Base.length(clusters) => 3
end

facts("kmeans returns no empty clusters") do
  a = FeatureVector(Dict(["hello"=>4, "word"=>0]))
  b = FeatureVector(Dict(["hello"=>0, "word"=>5]))
  c = FeatureVector(Dict(["kim"=>4, "julia"=>8]))
  dict = Dict([1=>a,2=>b,3=>c])
  clusters = kmeans(dict,[],3)

  for cluster in clusters
    @fact cluster.vectors != Cluster().vectors => true
  end
end

facts("kmeans returns the correct clustering of features") do
  a = FeatureVector(Dict(["hello"=>4, "word"=>3]))
  b = FeatureVector(Dict(["part"=>3, "science"=>5]))
  c = FeatureVector(Dict(["kim"=>4, "julia"=>8]))
  aa = FeatureVector(Dict(["hello"=>7, "word"=>5]))
  bb = FeatureVector(Dict(["part"=>8, "science"=>9]))
  cc = FeatureVector(Dict(["kim"=>2, "julia"=>3]))
  aaa = FeatureVector(Dict(["hello"=>1, "word"=>3]))
  bbb = FeatureVector(Dict(["part"=>3, "science"=>1]))
  ccc = FeatureVector(Dict(["kim"=>15, "julia"=>12]))
  dict = Dict([1=>a,2=>b,3=>c,4=>aa,5=>bb,6=>cc,7=>aaa,8=>bbb,9=>ccc])
  clusters = kmeans(dict,[],3)

  # what the clusters should be
  clust1_values = [a,aa,aaa]
  clust2_values = [b,bb,bbb]
  clust3_values = [c,cc,ccc]
  array = [clust1_values,clust2_values,clust3_values]

  # getting each cluster's FVs and making sure they are all there.
  for cluster in clusters
    fvs = collect(Base.values(cluster.vectors))
    for value in fvs
      @fact value in array => true
    end
  end
end

facts("kmeans lets you put in array of pre-defined centroids") do
  a = FeatureVector(Dict(["hello"=>4, "word"=>3]))
  b = FeatureVector(Dict(["part"=>3, "science"=>5]))
  c = FeatureVector(Dict(["kim"=>4, "julia"=>8]))
  aa = FeatureVector(Dict(["hello"=>7, "word"=>5]))
  bb = FeatureVector(Dict(["part"=>8, "science"=>9]))
  cc = FeatureVector(Dict(["kim"=>2, "julia"=>3]))
  aaa = FeatureVector(Dict(["hello"=>1, "word"=>3]))
  bbb = FeatureVector(Dict(["part"=>3, "science"=>1]))
  ccc = FeatureVector(Dict(["kim"=>15, "julia"=>12]))
  dict = Dict([1=>a,2=>b,3=>c,4=>aa,5=>bb,6=>cc,7=>aaa,8=>bbb,9=>ccc])
  array = [a,b,c]
  clusters = kmeans(dict,array)

  # what the clusters should be
  clust1_values = [a,aa,aaa]
  clust2_values = [b,bb,bbb]
  clust3_values = [c,cc,ccc]
  array = [clust1_values,clust2_values,clust3_values]

  # getting each cluster's FVs and making sure they are all there.
  for cluster in clusters
    fvs = collect(Base.values(cluster.vectors))
    for value in fvs
      @fact value in array => true
    end
  end
end

ext: jl

max_stars_repo_name: drawkward/SLU-TextMining

max_issues_repo_name: drawkward/SLU-TextMining

max_forks_repo_name: drawkward/SLU-TextMining


text:
# Estimates the basic Nerlove Cobb-Douglas model by MLE
# this is intended to show that MLE with normality is the
# ML estimator (for the b"s, not for sig^2)
using Econometrics, DelimitedFiles
data = readdlm("../Data/nerlove.data")
data = data[:,2:end]
data = log.(data)
n = size(data,1)
y = data[:,1]
x = data[:,2:end]
x = [ones(n,1) x]
names = ["constant", "output", "labor", "fuel", "capital", "sig"]
theta = [zeros(size(x,2)); 1.0] # start values for estimation
model = theta -> normal(theta, y, x)
thetahat, objvalue, V, converged = mleresults(model, theta, "estimate Nerlove model by MLE", names)
nothing

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
using Graphs
n = 10
k = 2
beta = 0.3
svt.S
g = watts_strogatz_graph(n, k, beta)
lap = Symmetric(laplacian_matrix(g))
laps = Symmetric(laplacian_matrix_sparse(g))
lapsp = laplacian_matrix_sparse(g)
function algConnectivity(L::Symmetric;vecs=false)
    # Calculates the algebraic connectivity of graph with laplacian L
    # Input: L a symmetric laplacian matrix
    # Output: value of algebraic connectivity λ_2
    e_vals = eigfact(L,2:2)
    λ = e_vals.values
    ν = e_vals.vectors
    if vecs
        return λ, ν
    else
        return λ
    end
end
function algConnectivity(L::SparseMatrixCSC;vecs=false)
    # Calculates the algebraic connectivity of graph with laplacian L
    # Input: L a symmetric laplacian matrix
    s = svds(L)
    # Output: value of algebraic connectivity λ_2
    λ = s[2][end-1]
    ν = s[1][:,end-1]
    if vecs
        return λ, ν
    else
        return λ
    end
end
@time connectivity = algConnectivity(lapsp,vecs=false)
println(size(lap))
println(typeof(lap))
plot(g)
issym(laps)
laps.data[3,1] = 90
laps.uplo
s = svds(lapsp)
s[1]
sv = svd(full(lapsp))
svt = svdfact(full(lapsp))


text:
# ------------------------------------------------------------------------------
# Jacobi smoother example
# ------------------------------------------------------------------------------

using LFAToolkit
using LinearAlgebra

# setup
mesh = Mesh2D(1.0, 1.0)
p = 3

# diffusion operator
diffusion = GalleryOperator("diffusion", p + 1, p + 1, mesh)

# Jacobi smoother
jacobi = Jacobi(diffusion)

# compute operator symbols
A = computesymbols(jacobi, [1.0], [π, π])
eigenvalues = real(eigvals(A))

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


text:
#---------------------------------------------------------
# # [SPECTrecon overview](@id 1-overview)
#---------------------------------------------------------

# This page explains the Julia package
# [`SPECTrecon`](https://github.com/JeffFessler/SPECTrecon.jl).

# ### Setup

# Packages needed here.

using SPECTrecon
using MIRTjim: jim, prompt
using Plots: scatter, plot!, default; default(markerstrokecolor=:auto)
using Plots # @animate, gif

# The following line is helpful when running this example.jl file as a script;
# this way it will prompt user to hit a key after each figure is displayed.

isinteractive() ? jim(:prompt, true) : prompt(:draw);

# ### Overview

# To perform SPECT image reconstruction,
# one must have a model for the imaging system
# encapsulated in a forward projector and back projector.

# Mathematically, we write the forward projection process in SPECT
# as "y = A * x" where A is a "system matrix"
# that models the physics of the imaging system
# (including depth-dependent collimator/detector response
# and attenuation)
# and "x" is the current guess of the emission image.

# However, in code we usually cannot literally store "A"
# as dense matrix because it is too large.
# A typical size in SPECT is that
# the image `x` is
# `nx × ny × nz = 128 × 128 × 100`
# and the array of projection views `y` is
# `nx × nz × nview = 128 × 100 × 120`.
# So the system matrix `A` has `1536000 × 1638400` elements
# which is far to many to store,
# even accounting for some sparsity.

# Instead, we write functions called forward projectors
# that calculate `A * x` "on the fly".

# Similarly, the operation `A' * y`
# is called "back projection",
# where `A'` denotes the transpose or "adjoint" of `A`.

# ### Example

# To illustrate forward and back projection,
# it is easiest to start with a simulation example
# using a digital phantom.
# The fancy way would be to use a 3D phantom from
# [ImagePhantoms](https://github.com/JuliaImageRecon/ImagePhantoms.jl),
# but instead we just use two simple cubes.

nx,ny,nz = 128,128,80
T = Float32
xtrue = zeros(T, nx,ny,nz)
xtrue[(1nx÷4):(2nx÷3), 1ny÷5:(3ny÷5), 2nz÷6:(3nz÷6)] .= 1
xtrue[(2nx÷5):(3nx÷5), 1ny÷5:(2ny÷5), 4nz÷6:(5nz÷6)] .= 2

average(x) = sum(x) / length(x)
function mid3(x::AbstractArray{T,3}) where {T}
    (nx,ny,nz) = size(x)
    xy = x[:,:,ceil(Int, nz÷2)]
    xz = x[:,ceil(Int,end/2),:]
    zy = x[ceil(Int, nx÷2),:,:]'
    return [xy xz; zy fill(average(x), nz, nz)]
end
jim(mid3(xtrue), "Middle slices of x")

# ### PSF

# Create a synthetic depth-dependent PSF for a single view
nx_psf = 11
psf1 = psf_gauss( ; nx, nx_psf, fwhm_end = 6)
jim(psf1, "PSF for each of $nx planes")

# In general the PSF can vary from view to view
# due to non-circular detector orbits.
# For simplicity, here we illustrate the case
# where the PSF is the same for every view.

nview = 60
psfs = repeat(psf1, 1, 1, 1, nview)

# Plan the PSF modeling (see `3-psf.jl`)

plan = plan_psf(nx, nz, nx_psf)

# ### Basic SPECT forward projection

# Here is a simple illustration
# of a SPECT forward projection operation.
# (This is a memory inefficient way to do it!)

dy = 4 # transaxial pixel size in mm
mumap = zeros(T, size(xtrue)) # μ-map just zero for illustration here
views = project(xtrue, mumap, psfs, dy)

# Display the calculated (i.e., simulated) projection views

jim(views[:,:,1:4:end], "Every 4th of $nview projection views")

# ### Basic SPECT back projection

# This illustrates an "unfiltered backprojection"
# that leads to a very blurry image
# (again, with a simple memory inefficient usage).

# First, back-project two "rays"
# to illustrate the depth-dependent PSF.
tmp = zeros(T, size(views))
tmp[nx÷2, nz÷2, nview÷5] = 1
tmp[nx÷2, nz÷2, 1] = 1
tmp = backproject(tmp, mumap, psfs, dy)
jim(mid3(tmp), "Back-projection of two rays")

# Now back-project all the views of the phantom

back = backproject(views, mumap, psfs, dy)
jim(mid3(back), "Back-projection of ytrue")

# ### Memory efficiency

# For iterative reconstruction,
# one must do forward and back-projection repeatedly.
# It is more efficient to pre-allocate work arrays
# for those operations,
# instead of repeatedly making system calls.

# Here we illustrate the memory efficient versions
# that are recommended for iterative SPECT reconstruction.

# First construction the SPECT plan

#viewangle = (0:(nview-1)) * 2π # default
plan = SPECTplan(mumap, psfs, dy; T)

# Mutating version of forward projection:

tmp = Array{T}(undef, nx, nz, nview)
project!(tmp, xtrue, plan)
@assert tmp == views

# Mutating version of back-projection:

tmp = Array{T}(undef, nx, ny, nz)
backproject!(tmp, views, plan)
@assert tmp == back

# ### Using `LinearMapAA`

# Calling `project!` and `backproject!` repeatedly
# leads to application-specific code.
# More general code uses the fact that SPECT projection and back-projection
# are linear operations,
# so we use `LinearMapAA` to define a "system matrix" for these operations.

using LinearMapsAA: LinearMapAA
forw! = (y,x) -> project!(y, x, plan)
back! = (x,y) -> backproject!(x, y, plan)
idim = (nx,ny,nz)
odim = (nx,nz,nview)
A = LinearMapAA(forw!, back!, (prod(odim),prod(idim)); T, odim, idim)

# Simple forward and back-projection:
@assert A * xtrue == views
@assert A' * views == back

# Mutating version:
using LinearAlgebra: mul!
tmp = Array{T}(undef, nx, nz, nview)
mul!(tmp, A, xtrue)
@assert tmp == views
tmp = Array{T}(undef, nx, ny, nz)
mul!(tmp, A', views)
@assert tmp == back

# ### Units

# The pixel dimensions `deltas` can (and should!) be values with units.

# Here is an example ... (todo)
#using UnitfulRecipes
#using Unitful: mm

# ### Projection view animation

anim = @animate for i in 1:nview
    ymax = maximum(views)
    jim(views[:,:,i],
        "SPECT projection view $i of $nview",
        clim = (0, ymax),
    )
end
gif(anim, "views.gif", fps = 8)


text:
# redundant
struct QMDPState{B,S}
    isstate::Bool
    b::Union{Nothing,B}
    s::Union{Nothing,S}

    function QMDPState{B,S}(isstate::Bool, x) where {B,S}
        if isstate
            return new(isstate, nothing, x)
        else
            return new(isstate, x, nothing)
        end
    end
end

struct QMDPWrapper{M,B,S,A} <: MDP{QMDPState{B,S}, A}
    mdp::M
end

function QMDPWrapper(mdp::MDP, B::Type)
    S = statetype(mdp)
    A = actiontype(mdp)
    return QMDPWrapper{typeof(mdp), B, S, A}(mdp)
end

state(m::QMDPWrapper{<:Any,B,S,<:Any}, s::S) where {B,S} = QMDPState{B,S}(true, s)
state(m::QMDPWrapper{<:Any,B,S,<:Any}, b::B) where {B,S} = QMDPState{B,S}(false, b)

function POMDPs.isterminal(m::QMDPWrapper, s::QMDPState)
    if s.isstate
        return isterminal(m, s.s)
    else
        return isterminal(m, s.b)
    end
end

function POMDPs.gen(m::QMDPWrapper, s::QMDPState, a, rng::AbstractRNG)
    if s.isstate
        sp, r = @gen(:sp,:r)(m.mdp, s.s, a, rng)
    else
        sp, r = @gen(:sp,:r)(m.mdp, rand(rng, s.b), a, rng)
    end
    return (sp=statetype(m)(true, sp), r=r)
end

POMDPs.discount(m::QMDPWrapper) = discount(m.mdp)
function POMDPs.actions(m::QMDPWrapper, s::QMDPState)
    if s.isstate
        return actions(m.mdp, s.s)
    else
        # Gallium.@enter actions(m.mdp, s.b)
        return actions(m.mdp, s.b)
    end
end

struct GenQMDPSolver <: Solver
    solver
end

struct GenQMDPPolicy{P<:Policy, Q<:QMDPWrapper} <: Policy
    policy::P
    qmdp::Q
end

POMDPs.solve(sol::GenQMDPSolver, qmdp::QMDPWrapper) = GenQMDPPolicy(solve(sol.solver, qmdp), qmdp)

function POMDPModelTools.action_info(p::GenQMDPPolicy, b)
    # XXX if this ever makes it into the toolbox, need to get the mdp some other way
    s = QMDPState{typeof(b), statetype(p.qmdp.mdp)}(false, b)
    return action_info(p.policy, s)
end

POMDPs.action(p::GenQMDPPolicy, b) = first(action_info(p, b))


text:
f = open("input02.txt", "r")
    l = readlines(f)
close(f)

horizontal = 0
deep = 0
for a in l
    s = split(a, " ")
    if s[1] == "forward"
        global horizontal += parse(Int, s[2])
    elseif s[1] == "down"
        global deep += parse(Int, s[2])
    elseif s[1] == "up"
        global deep += -parse(Int, s[2])
    else
        println("Unknown entry in input")
    end
end
product = horizontal * deep
println("Part 1")
println("Horizontal : $horizontal")
println("Depth : $deep")
println("Product : $product")

horizontal = 0
deep = 0
aim = 0
for a in l
    s = split(a, " ")
    if s[1] == "down"
        global aim += parse(Int, s[2])
    elseif s[1] == "up"
        global aim += -parse(Int, s[2])
    elseif s[1] == "forward"
        global horizontal += parse(Int, s[2])
        global deep += (parse(Int, s[2]) * aim)
    else
        println("Unknown entry in input")
    end
end
product = horizontal * deep
println("Part 2")
println("Horizontal : $horizontal")
println("Depth : $deep")
println("Product : $product")

max_stars_repo_name: jung-benjamin/advent-of-code

max_issues_repo_name: jung-benjamin/advent-of-code

max_forks_repo_name: jung-benjamin/advent-of-code


text:
type GlmResp{T<:FP} <: ModResp               # response in a glm model
    y::Vector{T}                # response
    d::UnivariateDistribution
    l::Link
    devresid::Vector{T}         # (squared) deviance residuals
    eta::Vector{T}              # linear predictor
    mu::Vector{T}               # mean response
    mueta::Vector{T}            # derivative of mu w.r.t. eta
    offset::Vector{T}           # offset added to linear predictor (usually 0)
    var::Vector{T}              # (unweighted) variance at current mu
    wts::Vector{T}              # prior weights
    wrkresid::Vector{T}         # working residuals
    function GlmResp(y::Vector{T}, d::Distribution, l::Link, eta::Vector{T},
                     mu::Vector{T}, off::Vector{T}, wts::Vector{T})
        if isa(d, Binomial)
            for yy in y; 0. <= yy <= 1. || error("$yy in y is not in [0,1]"); end
        else
            insupport(d, y) || error("y must be in the support of d")
        end
        n = length(y)
        length(eta) == length(mu) == length(wts) == n || error("mismatched sizes")
        lo = length(off); lo == 0 || lo == n || error("offset must have length $n or length 0")
        res = new(y,d,l,Array(T,n),eta,mu,Array(T,n),off,Array(T,n),wts,Array(T,n))
        updatemu!(res, eta)
        res
    end
end

deviance(r::GlmResp) = sum(r.devresid)

devresid!(r::GlmResp) = devresid!(r.d, r.devresid, r.y, r.mu, r.wts)

linkfun!(r::GlmResp) = linkfun!(r.l, r.eta, r.mu)

linkinv!(r::GlmResp) = linkinv!(r.l, r.mu, r.eta)

mueta!(r::GlmResp) = mueta!(r.l, r.mueta, r.eta)

function updatemu!{T<:FP}(r::GlmResp{T}, linPr::Vector{T})
    n = length(linPr)
    length(r.offset) == n ? map!(Add(), r.eta, linPr, r.offset) : copy!(r.eta, linPr)
    linkinv!(r); mueta!(r); var!(r); wrkresid!(r); devresid!(r)
    sum(r.devresid)
end

updatemu!{T<:FP}(r::GlmResp{T}, linPr) = updatemu!(r, convert(Vector{T},vec(linPr)))

var!(r::GlmResp) = var!(r.d, r.var, r.mu)

wrkresid!(r::GlmResp) = map1!(Divide(), map!(Subtract(), r.wrkresid, r.y, r.mu), r.mueta)

function wrkresp(r::GlmResp)
    if length(r.offset) > 0 return map1!(Add(), map(Subtract(), r.eta, r.offset), r.wrkresid) end
    map(Add(), r.eta, r.wrkresid)
end

function wrkwt(r::GlmResp)
    length(r.wts) == 0 && return [r.mueta[i] * r.mueta[i]/r.var[i] for i in 1:length(r.var)]
    [r.wts[i] * r.mueta[i] * r.mueta[i]/r.var[i] for i in 1:length(r.var)]
end

type GlmMod <: LinPredModel
    fr::ModelFrame
    rr::GlmResp
    pp::LinPred
    ff::Formula
    fit::Bool
end

function coeftable(mm::GlmMod)
    cc = coef(mm)
    se = stderr(mm)
    zz = cc ./ se
    DataFrame({cc, se, zz, 2.0 * ccdf(Normal(), abs(zz))},
              ["Estimate","Std.Error","z value", "Pr(>|z|)"])
end

function confint(obj::GlmMod, level::Real)
    cft = coeftable(obj)
    hcat(coef(obj),coef(obj)) + cft["Std.Error"] *
    quantile(Normal(), (1. - level)/2.) * [1. -1.]
end
confint(obj::GlmMod) = confint(obj, 0.95)

deviance(m::GlmMod)  = deviance(m.rr)

function fit(m::GlmMod; verbose::Bool=false, maxIter::Integer=30, minStepFac::Real=0.001, convTol::Real=1.e-6)
    m.fit && return m
    maxIter >= 1 || error("maxIter must be positive")
    0 < minStepFac < 1 || error("minStepFac must be in (0, 1)")

    cvg = false; p = m.pp; r = m.rr
    scratch = similar(p.X.m)
    devold = updatemu!(r, linpred(delbeta!(p, wrkresp(r), wrkwt(r), scratch)))
    installbeta!(p)
    for i=1:maxIter
        f = 1.0
        dev = updatemu!(r, linpred(delbeta!(p, r.wrkresid, wrkwt(r), scratch)))
        while dev > devold
            f /= 2.; f > minStepFac || error("step-halving failed at beta0 = $beta0")
            dev = updatemu!(r, linpred(p, f))
        end
        installbeta!(p, f)
        crit = (devold - dev)/dev
        verbose && println("$i: $dev, $crit")
        if crit < convTol; cvg = true; break end
        devold = dev
    end
    cvg || error("failure to converge in $maxIter iterations")
    m.fit = true
    m
end

function glm(f::Formula, df::AbstractDataFrame, d::UnivariateDistribution, l::Link; dofit::Bool=true, wts=Float64[], offset=Float64[])
    mf = ModelFrame(f, df)
    mm = ModelMatrix(mf)
    y = model_response(mf); T = eltype(y);
    if T <: Integer
        y = float64(y)
        T = Float64
    end
    n = length(y); lw = length(wts)
    lw == 0 || lw == n || error("length(wts) = $lw should be 0 or $n")
    w = lw == 0 ? ones(T,n) : (T <: Float64 ? copy(wts) : convert(Vector{T}, wts))
    mu = mustart(d, y, w)
    off = T <: Float64 ? copy(offset) : convert(Vector{T}, offset)
    rr = GlmResp{T}(y, d, l, linkfun!(l, similar(mu), mu), mu, off, w)
    res = GlmMod(mf, rr, DensePredChol(mm), f, false)
    dofit ? fit(res) : res
end

glm(e::Expr, df::AbstractDataFrame, d::UnivariateDistribution, l::Link) = glm(Formula(e),df,d,l)
glm(e::Expr, df::AbstractDataFrame, d::UnivariateDistribution) = glm(Formula(e),df,d,canonicallink(d))
glm(f::Formula, df::AbstractDataFrame, d::UnivariateDistribution) = glm(f, df, d, canonicallink(d))
glm(s::String, df::AbstractDataFrame, d::UnivariateDistribution) = glm(Formula(parse(s)[1]), df, d)

## scale(m) -> estimate, s, of the scale parameter
## scale(m,true) -> estimate, s^2, of the squared scale parameter
scale(x::GlmMod,sqr::Bool=false) = 1. # generalize this - only appropriate for Bernoulli and Poisson


text:
#Correlations

#basis for three-photon correlation
basis = calculate_basis(L,LMAX,N,K)
#two-photon correlation
println("Calculate two photon correlation.")
c2 = twoPhotons(intensity, basis, K)
#dummy reference c3
c3ref = zeros(N,2*N,K,K,K)
println("calculate energy.")
#energy calculation including three-photon correlation
energy(intensity, basis, c3ref)


text:
using Turing
using ArviZ
import PyPlot as plt
using PrettyTables
using LinearAlgebra
## ========================================================================
# data
J = 8
y = [28.0, 8.0, -3.0, 7.0, -1.0, 1.0, 18.0, 12.0]
σ = [15.0, 10.0, 16.0, 11.0, 9.0, 11.0, 10.0, 18.0]
schools = [
    "Choate",
    "Deerfield",
    "Phillips Andover",
    "Phillips Exeter",
    "Hotchkiss",
    "Lawrenceville",
    "St. Paul's",
    "Mt. Hermon",
];

nwarmup, nsamples, nchains = 10_000, 1000, 4;

## ============================================================================
# model

Turing.@model function turing_model(y, σ, J=length(y))
    μ ~ Normal(10, 15)
    τ ~ InverseGamma2(5, 1)
    θ ~ filldist(Normal(μ, τ), J)
    for i in 1:J
        y[i] ~ Normal(θ[i], σ[i])
    end
end

param_mod = turing_model(y, σ)
sampler = NUTS(nwarmup, 0.8)

rng = Random.MersenneTwister()#16653
turing_chns = sample(
    rng, param_mod, sampler, MCMCThreads(), nsamples, nchains; progress=false
);

plot(turing_chns)
meanplot(turing_chns)

## ============================================================================
# autocorr
plot_autocorr(turing_chns; var_names=["μ", "τ"]);
plt.gcf()

# remove autocorr
turing_chns = turing_chns[1:20:end]
plot_autocorr(turing_chns; var_names=["μ", "τ"]);
plt.gcf()

##
idata = from_mcmcchains(
    turing_chns;
    #coords=Dict("school" => schools),
    #dims=Dict("y" => ["school"], "σ" => ["school"], "θ" => ["school"]),
    library="Turing",
)
idata.posterior

plot_trace(idata);
plt.gcf()

summarystats(idata) |> pretty_table
## ============================================================================
# analyse divergences
plot_pair(
    idata;
    #coords=Dict("school" => ["Choate", "Deerfield", "Phillips Andover"]),
    divergences=true,
);
plt.gcf()
## ============================================================================
# extract additional info
prior = sample(param_mod, Prior(), nsamples; progress=false)

# Instantiate the predictive model
param_mod_predict = turing_model(similar(y, Missing), σ)
# and then sample!
prior_predictive = predict(param_mod_predict, prior)
posterior_predictive = predict(param_mod_predict, turing_chns)

loglikelihoods = Turing.pointwise_loglikelihoods(
    param_mod, MCMCChains.get_sections(turing_chns, :parameters)
)

## ============================================================================
# add to ArviZ
# Ensure the ordering of the loglikelihoods matches the ordering of `posterior_predictive`
ynames = string.(keys(posterior_predictive))
loglikelihoods_vals = getindex.(Ref(loglikelihoods), ynames)
# Reshape into `(nchains, nsamples, size(y)...)`
loglikelihoods_arr = permutedims(cat(loglikelihoods_vals...; dims=3), (2, 1, 3))

idata = from_mcmcchains(
    turing_chns;
    posterior_predictive=posterior_predictive,
    log_likelihood=Dict("y" => loglikelihoods_arr),
    prior=prior,
    prior_predictive=prior_predictive,
    observed_data=Dict("y" => y),
    coords=Dict("school" => schools),
    #dims=Dict("y" => ["school"], "σ" => ["school"], "θ" => ["school"]),
    library="Turing",
)

loo(idata)

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
module CudaUtilsModule

import CUDAdrv
import CUDAdrv: CuModule, CuModuleFile, unsafe_unload!,
                CuFunction, cudacall, CuDevice, CuArray, CuDim

const ptxdict = Dict()
const mdlist = Array{CuModule}(0)

function mdinit(devlist, ptxfile)
  global ptxdict
  global mdlist

  for dev in devlist

    CuDevice(dev)
    md = CuModuleFile(ptxfile)

    # Utils functions
    ptxdict["setup_kernel"] = CuFunction(md, "setup_kernel")
    ptxdict["perturb"]      = CuFunction(md, "perturb")
    # ptxdict[(dev, "veccost")]      = CuFunction(md, "veccost")
    ptxdict["veccost2"]     = CuFunction(md, "veccost2")
    ptxdict["vec_add"]      = CuFunction(md, "vec_add")

    # ICM functions
    # ptxdict[(dev, "condition_icm")]  = CuFunction(md, "condition_icm");
    # ptxdict[(dev, "condition_icm2")] = CuFunction(md, "condition_icm2");
    ptxdict["condition_icm3"] = CuFunction(md, "condition_icm3");

    push!(mdlist, md)
  end
end

mdclose() = (empty!(mdlist); empty!(ptxdict))

function finit( )
  mdclose()
end

function init(devlist, ptxfile)
  mdinit(devlist, ptxfile)
end

# # Accumulates binaries to unaries for icm conditioning
# function condition_icm(
#   nblocks::Integer, nthreads::Integer,
#   d_ub::CudaArray{Cfloat, 2},     # in/out. h-by-n. where we accumulate
#   d_bb::CudaArray{Cfloat, 2},     # in. binaries
#   d_codek::CudaArray{Cuchar, 1},  # in. n-long. indices into d_bb
#   n::Cint,                        # in. size( d_ub, 2 )
#   h::Cint)                        # in. size( d_ub, 1 )
#
#   fun = ptxdict[(device(), "condition_icm")];
#
#   cudacall( fun, nblocks, nthreads,
#     (d_ub, d_bb, d_codek, n, h));
#
#   return nothing
# end

# function condition_icm2(
#   nblocks::Integer, nthreads::CUDArt.CudaDims,
#   d_ub::CudaArray{Cfloat, 2},     # in/out. h-by-n. where we accumulate
#   d_bb::CudaArray{Cfloat, 2},     # in. binaries
#   d_codek::CudaArray{Cuchar, 1},  # in. n-long. indices into d_bb
#   n::Cint,                        # in. size( d_ub, 2 )
#   h::Cint)                        # in. size( d_ub, 1 )
#
#   fun = ptxdict[(device(), "condition_icm2")];
#
#   cudacall( fun, nblocks, nthreads,
#     (d_ub, d_bb, d_codek, n, h));
#
#   return nothing
# end

function condition_icm3(
  nblocks::Integer, nthreads::CuDim,
  # d_ub::CuArray{Cfloat, 2},     # in/out. h-by-n. where we accumulate
  # d_bbs::CuArray{Cfloat, 2},     # in. binaries
  # d_codek::CuArray{Cuchar, 2},  # in. n-long. indices into d_bb
  d_ub::CUDAdrv.Mem.Buffer,     # in/out. h-by-n. where we accumulate
  d_bbs::CUDAdrv.Mem.Buffer,     # in. binaries
  d_codek::CUDAdrv.Mem.Buffer,  # in. n-long. indices into d_bb
  conditioning::Cint,
  m::Cint,
  n::Cint)                        # in. size( d_ub, 2 )

  fun = ptxdict["condition_icm3"]
  cudacall(fun, nblocks, nthreads,
    (Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cuchar}, Cint, Cint, Cint),
    d_ub, d_bbs, d_codek, conditioning, m, n)

  return nothing
end

function vec_add(
  nblocks::CuDim, nthreads::CuDim,
  d_matrix::CUDAdrv.Mem.Buffer,
  d_vec::CUDAdrv.Mem.Buffer,
  n::Cint,
  h::Cint)

  fun = ptxdict["vec_add"]
  cudacall(fun, nblocks, nthreads,
    (Ptr{Cfloat}, Ptr{Cfloat}, Cint, Cint),
    d_matrix, d_vec, n, h)
end

# function veccost(
#   nblocks::CUDArt.CudaDims, nthreads::CUDArt.CudaDims,
#   d_rx::CudaArray{Cfloat},
#   d_codebooks::CudaArray{Cfloat},
#   d_codes::CudaArray{Cuchar},
#   d_veccost::CudaArray{Cfloat}, # out.
#   m::Cint,
#   n::Cint)
#
#   fun = ptxdict[( device(), "veccost" )];
#   cudacall( fun, nblocks, nthreads,
#     ( d_rx, d_codebooks, d_codes, d_veccost, m, n ));
# end

function veccost2(
  nblocks::CuDim, nthreads::CuDim,
  d_rx::CUDAdrv.Mem.Buffer,
  d_codebooks::CUDAdrv.Mem.Buffer,
  d_codes::CUDAdrv.Mem.Buffer,
  d_veccost::CUDAdrv.Mem.Buffer, # out.
  d::Cint,
  m::Cint,
  n::Cint)

  fun = ptxdict["veccost2"]
  cudacall(fun, nblocks, nthreads,
    (Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cuchar}, Ptr{Cfloat}, Cint, Cint, Cint),
    d_rx, d_codebooks, d_codes, d_veccost, d, m, n,
    shmem=Int(d*sizeof(Cfloat)))
end

function perturb(
  nblocks::CuDim, nthreads::CuDim,
  state::CUDAdrv.Mem.Buffer,
  codes::CUDAdrv.Mem.Buffer,
  n::Cint,
  m::Cint,
  k::Cint)

  fun = ptxdict["perturb"]
  cudacall( fun, nblocks, nthreads,
    (Ptr{Void}, Ptr{Cuchar}, Cint, Cint, Cint),
    state, codes, n, m, k)
end

function setup_kernel(
  nblocks::CuDim, nthreads::CuDim,
  n::Cint,
  state::CUDAdrv.Mem.Buffer)

  fun = ptxdict["setup_kernel"]
  cudacall(fun, nblocks, nthreads,
    (Cint, Ptr{Void}),
    n, state)
end

end #cudaUtilsModule

max_issues_repo_name: xiaosu-zhu/local-search-quantization

max_forks_repo_name: xiaosu-zhu/local-search-quantization


text:
#
# provides various loss functions
#

# 0-1 loss
loss01(y,z) = float(y!=z)
# define truncated log loss
clip(z) = max(min(z,0.99),0.01)
logloss(y,z) = -1.0*(y*log(10,z)+(1.0-y)*log(10,1-z))
sqloss(y,z) = (y-z)^2
# note that this is in [0,2]
cliplogloss(y,z) = logloss(y,clip(z))
clipsqloss(y,z) = min((y-z)^2,1)

function empirical_loss(labels::Array{Float64,1},predictions::Array{Float64,1},lossfct)
    return sum(map(lossfct,labels,predictions))/length(labels)
end


text:
using Gadfly

set_default_plot_size(8inch, 3inch)

x = [0, 1, 1, 2, 2, 3, 3, 2, 2, 1, 1, 0, 4, 5, 5, 4]
y = [0, 0, 1, 1, 0, 0, 3, 3, 2, 2, 3, 3, 0, 0, 3, 3]
group = reduce(vcat, fill.(["H", "I"], [12,4]))

p1 = plot(x=x, y=y, group=group, Geom.polygon(preserve_order=false, fill=true))
p2 = plot(x=x, y=y, group=group, Geom.polygon(preserve_order=true, fill=false),
    Theme(line_width=2mm), linestyle=[:dash], color=[colorant"orange"] )

hstack(p1, p2)


text:
# map θ, S, p to sigma 1 surfaces. Test for one ocean column.
# This is a wrapper routine to read files on poseidon.
# ggebbie, 7-Apr-2021

using Revise
using ECCOtour, Dierckx, PyPlot

pz = 0:500:4000 # pressure levels
nz = length(pz)
p₀ = 1000 # reference pressure = 1000 dbar
θz = range(20,stop=10,length=nz)
Sz = range(36,stop=35,length=nz)
sig1grida = 24:0.05:31
sig1gridb = 31.02:0.02:33
sig1grid = vcat(sig1grida,sig1gridb)
sig1grid = sig1grid[1:3:end]
splorder = 3

@time θσ,Sσ,pσ=TSP2sigma1column(collect(θz),collect(Sz),collect(pz),sig1grid,splorder)
@time θσ,Sσ,pσ=TSP2sigma1column(θz,Sz,pz,sig1grid,splorder)

figure(100)
plot(θz,-pz)
plot(θσ,-pσ,"o")

figure(101)
plot(Sz,-pz)
plot(Sσ,-pσ,"o")


text:
# AbstractFFTs interface

const DUALFORW = true
const DUALBACK = false
const INPLACE = true

# ---------------------------------------------------------------------------- #
# DualPlan type (for AbstractFFTs interface)
# ---------------------------------------------------------------------------- #
mutable struct DualPlan{T,forward,inplace} <: Plan{T}
    region
    msize
    pinv::ScaledPlan{T}
    DualPlan{T,forward,inplace}(region,msize) where {T,forward,inplace} = new{T,forward,inplace}(region,msize)
end

# ---------------------------------------------------------------------------- #
# AbstractFFTs interface functions
# ---------------------------------------------------------------------------- #
AbstractFFTs.plan_fft(  Z::Array{Complex{D}}, region=1:ndims(Z); kwargs...) where D <: Dual = DualPlan{Complex{D},DUALFORW,!INPLACE}(region,size(Z))
AbstractFFTs.plan_fft!( Z::Array{Complex{D}}, region=1:ndims(Z); kwargs...) where D <: Dual = DualPlan{Complex{D},DUALFORW,INPLACE}(region,size(Z))
AbstractFFTs.plan_bfft( Z::Array{Complex{D}}, region=1:ndims(Z); kwargs...) where D <: Dual = DualPlan{Complex{D},DUALBACK,!INPLACE}(region,size(Z))
AbstractFFTs.plan_bfft!(Z::Array{Complex{D}}, region=1:ndims(Z); kwargs...) where D <: Dual = DualPlan{Complex{D},DUALBACK,INPLACE}(region,size(Z))

AbstractFFTs.plan_inv(p::DualPlan{Complex{D},forward,inplace}) where {D<:Dual,forward,inplace} = ScaledPlan(DualPlan{Complex{D},!forward,inplace}(p.region,p.msize),
           normalization(basetype(D), p.msize, p.region))

LinearAlgebra.mul!(Y::Array{Complex{D}}, p::DualPlan{Complex{D},DUALFORW,inplace}, X::Array{Complex{D}}) where {D<:Dual,inplace} = (copy!(Y,X); dualfft!(Y,p.region); return Y)
LinearAlgebra.mul!(Y::Array{Complex{D}}, p::DualPlan{Complex{D},DUALBACK,inplace}, X::Array{Complex{D}}) where {D<:Dual,inplace} = (copy!(Y,X); dualbfft!(Y,p.region); return Y)

Base.:*(p::DualPlan{Complex{D},DUALFORW,!INPLACE}, X::Array{Complex{D}}) where D <: Dual = dualfft(X,p.region)
Base.:*(p::DualPlan{Complex{D},DUALFORW,INPLACE}, X::Array{Complex{D}}) where D <: Dual = (dualfft!(X,p.region); return X)
Base.:*(p::DualPlan{Complex{D},DUALBACK,!INPLACE}, X::Array{Complex{D}}) where D <: Dual = dualbfft(X,p.region)
Base.:*(p::DualPlan{Complex{D},DUALBACK,INPLACE}, X::Array{Complex{D}}) where D <: Dual = (dualbfft!(X,p.region); return X)


text:
"""
    GammaCountExpert(m, s)

PMF:

```math
P(X = k) = G(m k, s T) - G(m (k+1), s T), \\quad \\text{ for } k = 0,1,2, \\ldots, n.
```

with

```math
G(m k, s T) = \\frac{1}{\\Gamma(mk)}  \\int^{sT}_{0} u^{mk - 1} e^{-u} du
```

See also: [Gamma Count Distribution](https://arxiv.org/abs/1312.2423) (Arxiv)

"""
struct GammaCountExpert{T<:Real} <: NonZIDiscreteExpert
    m::T
    s::T
    GammaCountExpert{T}(m, s) where {T<:Real} = new{T}(m, s)
end

function GammaCountExpert(m::T, s::T; check_args=true) where {T <: Real}
    check_args && @check_args(GammaCountExpert, m > zero(m) && s > zero(s))
    return GammaCountExpert{T}(m, s)
end

## Outer constructors
GammaCountExpert(m::Real, s::Real) = GammaCountExpert(promote(m, s)...)
GammaCountExpert(m::Integer, s::Integer) = GammaCountExpert(float(m), float(s))
GammaCountExpert() = GammaCountExpert(2.0, 1.0)

## Conversion
function convert(::Type{GammaCountExpert{T}}, m::S, s::S) where {T <: Real, S <: Real}
    GammaCountExpert(T(m), T(s))
end
function convert(::Type{GammaCountExpert{T}}, d::GammaCountExpert{S}) where {T <: Real, S <: Real}
    GammaCountExpert(T(d.m), T(d.s), check_args=false)
end
copy(d::GammaCountExpert) = GammaCountExpert(d.m, d.s, check_args=false)

## Loglikelihood of Expoert
logpdf(d::GammaCountExpert, x...) = isinf(x...) ? -Inf : LRMoE.logpdf.(LRMoE.GammaCount(d.m, d.s), x...)
pdf(d::GammaCountExpert, x...) = isinf(x...) ? 0.0 : LRMoE.pdf.(LRMoE.GammaCount(d.m, d.s), x...)
logcdf(d::GammaCountExpert, x...) = isinf(x...) ? 0.0 : LRMoE.logcdf.(LRMoE.GammaCount(d.m, d.s), x...)
cdf(d::GammaCountExpert, x...) = isinf(x...) ? 1.0 : LRMoE.cdf.(LRMoE.GammaCount(d.m, d.s), x...)

## expert_ll, etc
expert_ll_exact(d::GammaCountExpert, x::Real) = LRMoE.logpdf(LRMoE.GammaCount(d.m, d.s), x)
function expert_ll(d::GammaCountExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    # d_exp = LRMoE.GammaCount(d.m, d.s)
    expert_ll = (yl == yu) ? logpdf.(d, yl) : logcdf.(d, yu) + log1mexp.(logcdf.(d, ceil.(yl) .- 1) - logcdf.(d, yu))
    return expert_ll
end
function expert_tn(d::GammaCountExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    # d_exp = LRMoE.GammaCount(d.m, d.s)
    expert_tn = (tl == tu) ? logpdf.(d, tl) : logcdf.(d, tu) + log1mexp.(logcdf.(d, ceil.(tl) .- 1) - logcdf.(d, tu))
    return expert_tn
end
function expert_tn_bar(d::GammaCountExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    # d_exp = LRMoE.GammaCount(d.m, d.s)
    expert_tn_bar = (tl == tu) ? log1mexp.(logpdf.(d, tl)) : log1mexp.(logcdf.(d, tu) + log1mexp.(logcdf.(d, ceil.(tl) .- 1) - logcdf.(d, tu)))
    return expert_tn_bar
end

exposurize_expert(d::GammaCountExpert; exposure = 1) = GammaCountExpert(d.m, d.s/exposure)

## Parameters
params(d::GammaCountExpert) = (d.m, d.s)
function params_init(y, d::GammaCountExpert)

    # function init_obj(logparams, y)
    #     n = length(y)
    #     m_tmp = exp(logparams[1])
    #     s_tmp = exp(logparams[2])

    #     return -1 * ( sum(logpdf.(GammaCountExpert(m_tmp, s_tmp), y) ))
    # end

    # logparams_init = Optim.minimizer( Optim.optimize(x -> init_obj(x, y),  [log(2.0), 0.0] ))

    # m_init = exp(logparams_init[1])
    # s_init = exp(logparams_init[2])

    μ, σ2 = mean(y), var(y)
    s_init = σ2 / μ
    m_init = μ/(s_init)

    try
        GammaCountExpert(m_init, s_init)
    catch;
        GammaCountExpert()
    end
end

## Simululation
sim_expert(d::GammaCountExpert) = Distributions.rand(LRMoE.GammaCount(d.m, d.s), 1)[1]

## penalty
penalty_init(d::GammaCountExpert) = [2.0 10.0 2.0 10.0]
no_penalty_init(d::GammaCountExpert) = [1.0 Inf 1.0 Inf]
penalize(d::GammaCountExpert, p) = (p[1]-1)*log(d.m) - d.m/p[2] + (p[3]-1)*log(d.s) - d.s/p[4]

## statistics
mean(d::GammaCountExpert) = mean(LRMoE.GammaCount(d.m, d.s))
var(d::GammaCountExpert) = var(LRMoE.GammaCount(d.m, d.s))
quantile(d::GammaCountExpert, p) = quantile(LRMoE.GammaCount(d.m, d.s), p)

## Misc functions for E-Step

function _sum_dens_series(m_new, s_new, d::GammaCountExpert, yl, yu, exposure)
    upper_finite = isinf(yu) ? quantile(d, 1-1e-8) : yu
    series = yl:(max(yl, min(yu, upper_finite+1)))
    return sum(logpdf.(GammaCountExpert(m_new, s_new/exposure), series) .* pdf.(d, series))[1]
end

function _int_obs_dens_raw(m_new, s_new, d::GammaCountExpert, yl, yu, exposure)
    return _sum_dens_series(m_new, s_new, d, yl, yu, exposure)
end

function _int_lat_dens_raw(m_new, s_new, d::GammaCountExpert, tl, tu, exposure)
    return (tl==0 ? 0.0 : _sum_dens_series.(m_new, s_new, d, 0, ceil(tl)-1), exposure) +
        (isinf(tu) ? 0.0 : _sum_dens_series.(m_new, s_new, d, floor(tu)+1, Inf, exposure))
end

function _gammacount_optim_params(lognew,
                        d_old,
                        tl, yl, yu, tu,
                        # expert_ll_pos, expert_tn_pos, expert_tn_bar_pos,
                        exposure,
                        z_e_obs, z_e_lat, k_e; # ,
                        # Y_e_obs, Y_e_lat;
                        penalty = true, pen_pararms_jk = [1.0 Inf 1.0 Inf])
    # Optimization in log scale for unconstrained computation
    m_tmp = exp(lognew[1])
    s_tmp = exp(lognew[2])

    # Further E-Step
    # yl_yu_unique = unique_bounds(yl, yu)
    # int_obs_dens_tmp = _int_obs_dens_raw.(m_tmp, s_tmp, d_old, yl_yu_unique[:,1], yl_yu_unique[:,2])
    # densY_e_obs = exp.(-expert_ll_pos) .* int_obs_dens_tmp[match_unique_bounds(hcat(vec(yl), vec(yu)), yl_yu_unique)]
    # nan2num(densY_e_obs, 0.0) # get rid of NaN

    # tl_tu_unique = unique_bounds(tl, tu)
    # int_lat_dens_tmp = _int_lat_dens_raw.(m_tmp, s_tmp, d_old, tl_tu_unique[:,1], tl_tu_unique[:,2])
    # densY_e_lat = exp.(-expert_tn_bar_pos) .* int_lat_dens_tmp[match_unique_bounds(hcat(vec(tl), vec(tu)), tl_tu_unique)]
    # nan2num(densY_e_lat, 0.0) # get rid of NaN

    densY_e_obs = fill(0.0, length(yl))
    densY_e_lat = fill(0.0, length(yl))
    d_tmp = LRMoE.GammaCountExpert(m_tmp, s_tmp)

    for i in 1:length(yl)
        d_expo = exposurize_expert(d_tmp, exposure = exposure[i])
        # expert_ll_pos = expert_ll(d_expo, tl[i], yl[i], yu[i], tu[i])
        # expert_tn_bar_pos = expert_tn_bar(d_expo, tl[i], yl[i], yu[i], tu[i])
        densY_e_obs[i] = expert_ll(d_expo, tl[i], yl[i], yu[i], tu[i])
        # _int_obs_dens_raw(m_tmp, s_tmp, d_expo, yl[i], yu[i], exposure[i]) # exp(-expert_ll_pos) * _int_obs_dens_raw(m_tmp, s_tmp, d_expo, yl[i], yu[i], exposure[i])
        densY_e_lat[i] = expert_tn_bar(d_expo, tl[i], yl[i], yu[i], tu[i])
        # _int_lat_dens_raw(m_tmp, s_tmp, d_expo, yl[i], yu[i], exposure[i]) # exp(-expert_tn_bar_pos) * _int_lat_dens_raw(m_tmp, s_tmp, d_expo, yl[i], yu[i], exposure[i])
    end

    nan2num(densY_e_obs, 0.0) # get rid of NaN
    nan2num(densY_e_lat, 0.0) # get rid of NaN
    inf2num(densY_e_obs, 0.0) # get rid of inf
    inf2num(densY_e_lat, 0.0) # get rid of inf

    # term_zkz = z_e_obs .+ (z_e_lat .* k_e)
    term_zkz_Y = (z_e_obs .* densY_e_obs) .+ (z_e_lat .* k_e .* densY_e_lat)

    # sum_term_zkz = sum(term_zkz)[1]
    sum_term_zkzy = sum(term_zkz_Y)[1]

    obj = sum_term_zkzy
    p = penalty ? (pen_pararms_jk[1]-1)*log(m_tmp) - m_tmp/pen_pararms_jk[2] + (pen_pararms_jk[3]-1)*log(s_tmp) - s_tmp/pen_pararms_jk[4] : 0.0
    return (obj + p) * (-1.0)
end

## EM: M-Step
function EM_M_expert(d::GammaCountExpert,
                    tl, yl, yu, tu,
                    exposure,
                    z_e_obs, z_e_lat, k_e;
                    penalty = true, pen_pararms_jk = [1.0 Inf 1.0 Inf])

    # Update parameters
    logparams_new = Optim.minimizer( Optim.optimize(x -> _gammacount_optim_params(x, d,
                                                tl, yl, yu, tu,
                                                exposure,
                                                z_e_obs, z_e_lat, k_e,
                                                # Y_e_obs, Y_e_lat,
                                                penalty = penalty, pen_pararms_jk = pen_pararms_jk),
                                                # [log(d.m)-2.0, log(d.s)-2.0],
                                                # [log(d.m)+2.0, log(d.s)+2.0],
                                                [log(d.m), log(d.s)] ))
    # println("$logparams_new")
    m_new = exp(logparams_new[1])
    s_new = exp(logparams_new[2])

    # println("$m_new, $s_new")

    # Deal with zero mass
    if ccdf.(Gamma((0+1)*s_new, 1), m_new*s_new) > 0.999999 || isnan(ccdf.(Gamma((0+1)*s_new, 1), m_new*s_new))
        m_new, s_new = d.m, d.s
    end

    return GammaCountExpert(m_new, s_new)

end

## EM: M-Step, exact observations
function _gammacount_optim_params(lognew,
                        ye, exposure,
                        z_e_obs;
                        penalty = true, pen_pararms_jk = [1.0 Inf 1.0 Inf])
    # Optimization in log scale for unconstrained computation
    m_tmp = exp(lognew[1])
    s_tmp = exp(lognew[2])

    # Further E-Step
    densY_e_obs = fill(NaN, length(exposure))
    for i in 1:length(exposure)
        densY_e_obs[i] = expert_ll_exact(exposurize_expert(LRMoE.GammaCountExpert(m_tmp, s_tmp), exposure = exposure[i]), ye[i])
    end
    # densY_e_obs = logpdf.(GammaCountExpert(m_tmp, s_tmp), ye)
    nan2num(densY_e_obs, 0.0) # get rid of NaN

    term_zkz_Y = z_e_obs .* densY_e_obs

    sum_term_zkzy = sum(term_zkz_Y)[1]

    obj = sum_term_zkzy
    p = penalty ? (pen_pararms_jk[1]-1)*log(m_tmp) - m_tmp/pen_pararms_jk[2] + (pen_pararms_jk[3]-1)*log(s_tmp) - s_tmp/pen_pararms_jk[4] : 0.0
    return (obj + p) * (-1.0)
end

function EM_M_expert_exact(d::GammaCountExpert,
                    ye, exposure,
                    z_e_obs;
                    penalty = true, pen_pararms_jk = [1.0 Inf 1.0 Inf])

    # Update parameters
    logparams_new = Optim.minimizer( Optim.optimize(x -> _gammacount_optim_params(x,
                                                ye, exposure,
                                                z_e_obs,
                                                penalty = penalty, pen_pararms_jk = pen_pararms_jk),
                                                [log(d.m), log(d.s)] ))

    # println("$logparams_new")
    m_new = exp(logparams_new[1])
    s_new = exp(logparams_new[2])

    # Deal with zero mass
    if ccdf.(Gamma((0+1)*s_new, 1), m_new*s_new) > 0.999999 || isnan(ccdf.(Gamma((0+1)*s_new, 1), m_new*s_new))
        m_new, s_new = d.m, d.s
    end

    # println("$m_new, $s_new")
    return GammaCountExpert(m_new, s_new)

end

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
### functions to define useful regions of interest and load gene models

### differeing regions of interest
### coords are 1-based inclusive
function hk1_locus_coords()
    chrom = "chr10"
    mutations = 71108610:71108717

    genelocus = (71029740 - 4000):(71161640+4000)

    roi = 71107650:71110242
    criticalregion = 71108645:71108686
    criticalregion_4k = (71108645 - 2000):(71108686 + 2000)
    roi_200k = (first(roi) - 200_000):(last(roi) + 200_000)
    Dict("mutations"         => (chrom, mutations),
         "roi"               => (chrom, roi),
         "roi_200k"          => (chrom, roi_200k),
         "genelocus"         => (chrom, genelocus),
         "criticalregion"    => (chrom, criticalregion),
         "criticalregion_4k" => (chrom, criticalregion_4k),
         "widecontactview"   => (chrom, 70257650:71960242),
         "ins_contactview"   => ("chr11", 1330974:3032599),
         "sst_contactview"   => ("chr3", 186536699:188238191),
         "gcg_contactview"   => ("chr2", 162149392:163858914))
end

### load gtf file
function loadgtf(;file="gencode.v36lift37.annotation.gois.gtf.gz", projdir=getprojectdir(), fp = joinpath(projdir, "data", file))

    gtf = CSV.read(fp, DataFrame, comment="#", header=[:chrom, :source, :feature, :start, :stop, :dot, :strand, :frame, :description]);

    gtf = @subset(gtf, :feature .== "exon")
    dd = descriptiondict.(gtf.description)
    gtf[!, :GeneID] = getindex.(dd, "gene_id")
    gtf[!, :GeneName] = getindex.(dd, "gene_name")
    gtf[!, :Gene] = string.(gtf.GeneID, "|", gtf.GeneName);
    gtf[!, :TranscriptID] = getindex.(dd, "transcript_id")
    gtf[!, :TranscriptType] = getindex.(dd, "transcript_type")

    @subset(gtf, :TranscriptType .== "protein_coding")
end

### turn gtf file description field into a dictionary
function descriptiondict(des)
    fields = split(des, ";", keepempty=false)
    rc = [split(f, keepempty=false) for f in fields]
    Dict(id => replace(rec, "\"" => "") for (id, rec) in rc)
end

transcriptcoords(gtf) = combine(groupby(gtf, [:GeneName, :TranscriptID, :TranscriptType, :strand]), :start => minimum => :start, :stop => maximum => :stop, nrow => :count, [:start, :stop] => ((a, b) -> [collect(zip(a, b))]) => :coords)


text:
const RSDefaultParameters = ParamsDict(
    :PrecisionRatio    => 0.40, # 40% of the diameter is used as the initial step length
    :PrecisionTreshold => 1e-6  # They use 1e-6 in the paper.
)

"""
The variants of the memetic search algorithms RS and RIS.
However, we have modified them since they did not give very good performance when
implemented as described in the papers below. Possibly, the papers are not
unambigous and I have misinterpreted something from them...

The "Resampling Search" (RS) memetic algorithm is described in:

    F. Caraffini, F. Neri, M. Gongora and B. N. Passow, "Re-sampling Search: A
    Seriously Simple Memetic Approach with a High Performance", 2013.

and its close sibling "Resampling Inheritance Search" (RIS) is described in:

    F. Caraffini, F. Neri, B. N. Passow and G. Iacca, "Re-sampled Inheritance
    Search: High Performance Despite the Simplicity", 2013.
"""
mutable struct ResamplingMemeticSearcher{E<:Evaluator} <: SteppingOptimizer
    name::String
    params::Parameters
    evaluator::E
    resampling_func::Function

    precisions      # Cache the starting precision values so we need not calc them for each step
    diameters       # Cache the diameters...

    elite           # Current elite (best) candidate
    elite_fitness   # Fitness of current elite

    # Constructor for RS:
    function ResamplingMemeticSearcher(evaluator::E, parameters::Parameters,
            resampling_function::Function, name::String) where {E<:Evaluator}
        params = chain(RSDefaultParameters, parameters)
        elite = rand_individual(search_space(evaluator))
        diams = diameters(search_space(evaluator))
        new{E}(name, params, evaluator, resampling_function,
               params[:PrecisionRatio] * diams, diams,
               elite, fitness(elite, evaluator))
    end
end

name(rs::ResamplingMemeticSearcher) = rs.name

const RISDefaultParameters = ParamsDict(
    :InheritanceRatio => 0.30   # On average, 30% of positions are inherited when resampling in RIS
)

ResamplingMemeticSearcher(problem::OptimizationProblem,
        params::Parameters = EMPTY_PARAMS,
        resampling_function = random_resample,
        name = "Resampling Memetic Search (RS)") =
    ResamplingMemeticSearcher(ProblemEvaluator(problem), params, resampling_function, name)

# Constructor for the RIS:
ResamplingInheritanceMemeticSearcher(problem::OptimizationProblem, parameters::Parameters = EMPTY_PARAMS) =
    ResamplingMemeticSearcher(problem,
            chain(RSDefaultParameters, RISDefaultParameters, parameters),
            random_resample_with_inheritance,
            "Resampling Inheritance Memetic Search (RIS)")

resampling_memetic_searcher(problem::OptimizationProblem, params::Parameters) =
    ResamplingMemeticSearcher(problem, params)

resampling_inheritance_memetic_searcher(problem::OptimizationProblem, params::Parameters) =
    ResamplingInheritanceMemeticSearcher(problem, params)

# For Resampling Search (RS) the resample is purely random.
random_resample(rms::ResamplingMemeticSearcher) = rand_individual(search_space(rms.evaluator))

# For Resampling Inheritance Search (RIS) the resample has an inheritance component.
function random_resample_with_inheritance(rms::ResamplingMemeticSearcher)
    xt = random_resample(rms)
    n = numdims(rms.evaluator)
    i = rand(1:n)
    Cr = 0.5^(1/(rms.params[:InheritanceRatio]*n)) # See equation 3 in the RIS paper
    k = 1

    while rand() <= Cr && k < n
        xt[i] = rms.elite[i]
        i = 1 + mod(i, n)
        k += 1
    end

    return xt
end

function step!(rms::ResamplingMemeticSearcher)
    # First randomly sample two candidates and select the best one. It seems
    # RS and RIS might be doing this in two different ways but use the RS way for
    # now.
    trial, fitness = best_of(rms.resampling_func(rms), rms.resampling_func(rms), rms.evaluator)

    # Update elite if new trial is better. This is how they write it in the RIS paper
    # but in the RS paper it seems they always update the elite. Unclear! To me it
    # seems we should always update since we have already done a local search from
    # the current elite so it has a "head start" compared to new sampled points
    # which have not yet gone through local refinement. Since the evaluator/archive
    # keeps the best candidates anyway there is no risk for us in always overwriting the elite...
    set_as_elite_if_better(rms, trial, fitness)

    # Then run the local search on the elite one until step length too small.
    trial, fitness = local_search(rms, trial, fitness)
    set_as_elite_if_better(rms, trial, fitness)

    return trial, fitness
end

function set_as_elite_if_better(rms::ResamplingMemeticSearcher, candidate, fitness)
    if is_better(fitness, rms.elite_fitness, rms.evaluator)
        rms.elite = candidate
        rms.elite_fitness = fitness
        return true
    else
        return false
    end
end

stop_due_to_low_precision(rms::ResamplingMemeticSearcher, precisions) =
    norm(precisions ./ rms.diameters) < rms.params[:PrecisionTreshold]

function local_search(rms::ResamplingMemeticSearcher, xstart, fitness)
    ps = copy(rms.precisions)
    xt = copy(xstart)
    tfitness = copy(fitness)

    searchSpace = search_space(rms.evaluator)
    ssmins, ssmaxs = mins(searchSpace), maxs(searchSpace)
    n = numdims(rms.evaluator)

    indices = collect(1:n)

    while !stop_due_to_low_precision(rms, ps)
        xs = copy(xt)
        # We randomize the order that each decision var is changed. This is not done in the orig papers.
        shuffle!(indices)

        for j in 1:n
            i = indices[j]

            # This is how it is written in orig papers. To me it seems better to
            # take the step in a random direction; why prioritize one direction?
            xs[i] = xt[i] - ps[i]

            # rand bound from target if below min
            if xs[i] < ssmins[i]
                xs[i] = ssmins[i] + rand() * (xt[i] - ssmins[i])
            end

            if is_better(xs, tfitness, rms.evaluator)
                xt[i] = xs[i]
                tfitness = last_fitness(rms.evaluator)
            else
                xs[i] = xt[i] + ps[i]/2

                # rand bound from target if above max
                if xs[i] > ssmaxs[i]
                    xs[i] = xt[i] + rand() * (ssmaxs[i] - xt[i])
                end

                if is_better(xs, tfitness, rms.evaluator)
                    xt[i] = xs[i]
                    tfitness = last_fitness(rms.evaluator)
                end
            end
        end

        if is_better(tfitness, fitness, rms.evaluator)
            fitness = tfitness
            xstart = xt
        else
            ps = ps / 2
        end
    end

    return xstart, fitness
end


text:
### 17 OCT 2016: MUST RECODE THIS FUNCTION
#"""
#    fit_logistic(x, y, mask_n, λ)
#
#If called with an `Int` vector `mask_n`, then `fit_logistic()` will refit logistic effect sizes while masking components `y[i]` where `mask_n[i] = 0`.
#"""
#function fit_logistic{T <: Float}(
#    x        :: DenseMatrix{T},
#    y        :: DenseVector{T},
#    mask_n   :: DenseVector{Int},
#    λ   :: T;
#    n        :: Int    = length(y),
#    p        :: Int    = size(x,2),
#    mn       :: Int    = sum(mask_n),
#    d2b      :: DenseMatrix{T} = zeros(T, p,p),
#    x2       :: DenseMatrix{T} = zeros(T, n,p),
#    b        :: DenseVector{T} = zeros(T, p),
#    b0       :: DenseVector{T} = zeros(T, p),
#    ntb      :: DenseVector{T} = zeros(T, p),
#    db       :: DenseVector{T} = zeros(T, p),
#    xβ       :: DenseVector{T} = zeros(T, n),
#    lxβ      :: DenseVector{T} = zeros(T, n),
#    l2xb     :: DenseVector{T} = zeros(T, n),
#    tol      :: Float = convert(T, 1e-8),
#    max_iter :: Int   = 50,
#    quiet    :: Bool  = true,
#)
#
#    # if b is not warm-started, then ensure that it is not entirely zero
#    if all(b .== 0)
#        b[1] = logit(mean(y[mask_n .== 1]))
#    end
#
#    # initialize intermediate arrays for calculations
#    BLAS.gemv!('N', one(T), x, b, zero(T), xβ)
#    log2xb!(lxβ, l2xb, xβ, n=n)
#    mask!(lxβ, mask_n, 0, zero(T), n=n)
##    lxβ[mask_n .== 0] = zero(T)
#    copy!(b0,b)
#    fill!(db, zero(T))
#
#    # track objective
#    old_obj = Inf
#    new_obj = logistic_loglik(xβ,y,b,mask_n,0,λ,p, n=n)
#
#    # output progress to console
#    quiet || println("Iter\tHalves\tObjective")
#
#    i = 0
#    bktrk = 0
#    # enter loop for Newton's method
#    for i = 1:max_iter
#
#        # db = (x'*(lxβ - y)) / n + λ*b
#        BLAS.axpy!(n, -one(T), sdata(y), 1, sdata(lxβ), 1)
#        mask!(lxβ, mask_n, 0, zero(T), n=n)
##        lxβ[mask_n .== 0] = zero(T)
#        BLAS.gemv!('T', one(T), sdata(x), sdata(lxβ), zero(T), sdata(db))
##        BLAS.scal!(p, 1/n, sdata(db), 1)
#        BLAS.scal!(p, 1/mn, sdata(db), 1)
#        BLAS.axpy!(p, λ, sdata(b), 1, sdata(db), 1)
#
#        # d2b = (x'*diagm(l2xb)*x)/n + λ*I
#        # note that log2xb!() already performs division by n on l2xb
#        copy!(x2,x)
#        mask!(l2xb, mask_n, 0, zero(T), n=n)
##        l2xb[mask_n .== 0] = zero(T)
#        BLAS.scal!(p, n/mn, sdata(l2xb), 1) # rescale to number of unmasked samples
#        scale!(sdata(l2xb), sdata(x2))
#        BLAS.gemm!('T', 'N', one(T), sdata(x), sdata(x2), zero(T), sdata(d2b))
#        d2b += λ*I
#
#        # b = b0 - ntb = b0 - inv(d2b)*db
#        #   = b0 - inv[ x' diagm(pi) diagm(1 - pi) x + λ*I] [x' (pi - y) + λ*b]
#        ntb = d2b\db
#        copy!(b,b0)
#        BLAS.axpy!(p,-one(T),ntb,1,b,1)
#
#        # compute objective
#        new_obj = logistic_loglik(xβ,y,b,mask_n,0,λ,p, n=n)
#
#        # backtrack
#        j = 0
#        while (new_obj > old_obj + tol) && (j < 50)
#
#            # increment iterator
#            j += 1
#
#            # b = b0 - 0.5*ntb
#            copy!(b,b0)
##            BLAS.axpy!(p,-(0.5^j),ntb,1,b,1)
#            BLAS.axpy!(p,-one(T) / (2^j),ntb,1,b,1)
#
#            # recalculate objective
#            new_obj = logistic_loglik(xβ,y,b,mask_n,0,λ,p, n=n)
#
#        end
#
#        # accumulate total backtracking steps
#        bktrk += j
#
#        # track distance between iterates
#        dist = euclidean(sdata(b),sdata(b0)) / (norm(sdata(b0),2) + one(T))
#
#        # track progress
#        quiet || println(i, "\t", j, "\t", dist)
#
#        # check for convergence
#        # if converged, then return b
#        dist < tol && return b, i, div(bktrk,i)
#
#        # unconverged at this point, so update intermediate arrays
#        BLAS.gemv!('N', one(T), sdata(x), sdata(b), zero(T), sdata(xβ))
#        log2xb!(lxβ, l2xb, xβ, n=n)
#
#        # save previous beta
#        copy!(b0, b)
#        old_obj = new_obj
#    end
#
#    warn("fit_logistic failed to converge in $(max_iter) iterations, exiting...")
#    return b, i, div(bktrk,max_iter)
#end

### MENTION ADDITION OF "ACTIVE" VECTOR OUTPUT? ###
"""
If used with a `BEDFile` object `x`, then the additional optional arguments are:

- `pids`, a vector of process IDs. Defaults to `procs()`.
- `means`, a vector of SNP means. Defaults to `mean(T, x, shared=true, pids=procs()`.
- `invstds`, a vector of SNP precisions. Defaults to `invstd(x, means, shared=true, pids=procs()`.
"""
function L0_log{T <: Float, V <: SharedVector}(
    x        :: BEDFile{T},
    y        :: V,
    k        :: Int;
    pids     :: Vector{Int} = procs(x),
    mask_n   :: DenseVector{Int} = ones(Int, n),
    lambda   :: T    = convert(T, sqrt(log(p)/n)),
    mu       :: T    = one(T),
    tol      :: T    = convert(T, 1e-6),
    tolG     :: T    = convert(T, 1e-3),
    tolrefit :: T    = convert(T, 1e-6),
    max_iter :: Int  = 100,
    max_step :: Int  = 100,
    refit    :: Bool = true,
    quiet    :: Bool = true,
)

    # start timer
    tic()

    # check arguments
    n        == size(x,1) || throw(ArgumentError("Length n = $n of response vector y does not match number of rows = $(size(x,1)) in x"))
    k        >  p         && throw(ArgumentError("Value of argument k = $k exceeds number of predictors p = $p"))
    lambda   <  zero(T)   && throw(ArgumentError("Value of argument lambda = $lambda must be nonnegative"))
    mu       >  zero(T)   || throw(ArgumentError("Value of argument mu must be positive"))
    tol      >  eps(T)    || throw(ArgumentError("Value of argument tol must exceed machine precision"))
    tolG     >  eps(T)    || throw(ArgumentError("Value of argument tolG must exceed machine precision"))
    tolrefit >  eps(T)    || throw(ArgumentError("Value of argument tolrefit must exceed machine precision"))
    max_iter >= 0         || throw(ArgumentError("Value of max_iter must be nonnegative\n"))
    max_step >= 0         || throw(ArgumentError("Value of max_step must be nonnegative\n"))

    # enforce all 0/1 in mask_n
    sum((mask_n .== 1) $ (mask_n .== 0)) == n || throw(ArgumentError("Argument mask_n can only contain 1s and 0s"))
    mn = sum(mask_n),

    # initialize return values
    iter      = 0               # number of iterations of L0_reg
    exec_time = zero(T)         # compute time *within* L0_reg
    loss      = convert(T, Inf) # loss function value

    # initialize algorithm parameters
    num_df    = min(p,3*k)      # largest 3*k active components of gradient
    short_df  = min(p,2*k)      # largest 2*k active components of gradient
    converged = false           # is algorithm converged?
    normdf    = convert(T, Inf) # norm of active portion of gradient
    loss0     = convert(T, Inf) # previous loss function value
    bktrk     = 0               # track loss and number of backtracking steps
    nt_iter   = 0               # number of Newton iterations
    lt        = length(active)  # size of current active set?

    # formatted output to monitor algorithm progress
    quiet || print_header_log()

    # main GraSP iterations
    for iter = 1:max_iter

        # notify and break if maximum iterations are reached.
        if iter >= max_iter || stuck

            # warn about hitting maximum iterations
            iter >= max_iter && print_maxiter(max_iter, loss)

            # warn about cycling
            stuck && warn("L0_log appears to be cycling after $iter iterations, aborting early...\nCurrent loss: $loss\n")

            # send elements below tol to zero
            threshold!(v.b, tol)

            # if requested, apply final refit without regularization
            if refit
                copy!(v.idxs0, v.idxs)
                update_indices!(v.idxs, v.b)
                !isequal(v.idxs, v.idxs0) && copy!(v.xk, view(x, :, v.idxs))
                nt_iter, bktrk = fit_logistic!(v, y, zero(T), tol=tolrefit, max_iter=max_step, quiet=true)
                v.b[v.idxs] = v.bk
            end

            # stop timer
            exec_time = toq()

            # these are output variables for function
            # wrap them into a Dict and return
            return IHTLogResults(exec_time, iter, loss, copy(v.b), copy(v.active))
        end

        # save previous loss, iterate
        loss00 = loss0
        loss0 = loss
        copy!(v.b0, v.b)

        # size of current active set?
        if iter > 1
            active = union(v.dfidxs[1:short_df], v.bidxs[1:k])
        end
        lt = length(active)

        # update x*b
        # no need to compute anything if b = 0
        if all(v.b .== zero(T))
            fill!(v.xb,zero(T))
        else
            fill!(x.idxs2,false)
            v.idxs2[active] = true
            A_mul_B!(v.xb, x, v.b, v.idxs2, sum(v.idxs2), mask_n, pids=pids)
        end

        # recompute active loss = (-dot(y,Xb) + sum(log(1.0 + exp(Xb)))) / n + 0.5*lambda*sumabs2(b[active])
        # special case: b = 0 --> Xb = 0 --> loss = n*log(1 + exp(0))/n + 0.5*lambda*norm(0)
        if all(v.xb .== 0)
            loss = log(2)
        else
            #loss = logistic_loglik(v.xb, y, v.b, v.bidxs,mask_n,0,lambda,k, n=n)
            loss = logistic_loglik(v.xb, y, v.b, v.active, mask_n, 0, lambda, k)
        end

        # guard against numerical instabilities in loss function
        iter < 2 && check_finiteness(loss)

        # recompute active gradient df[active] = (x[:,active]'*(logistic(Xb) - y)) / n + lambda*b[active]
        # arrange calculations differently if active set is entire support 1, 2, ..., p
        if lt == p
            logistic_grad!(v.df, v.lxb, x, y, v.b, v.xb, mask_n, lambda, pids=pids)
        else
            logistic_grad!(v.df, v.lxb, x, y, v.b, v.xb, v.active, mask_n, lt, lambda)
        end

        # need norm of largest 3*k components of gradient
        normdf = df_norm(v.df, v.dfidxs, 1, num_df)

        # guard against numerical instabilities in gradient
        check_finiteness(normdf)

        # identify 3*k dominant directions in gradient
        selectperm!(dfidxs, df, 1:num_df, by=abs, rev=true, initialized=true)

        # clean b and fill, b[active] = b0[active] - mu*df[active]
        # note that sparsity level is size(active) which is one of [k, k+1, ..., 3*k]
        #update_x!(v.b, v.b0, v.df, v.active, mu, k=lt)
        fill!(v.b, zero(T))
        v.b[v.active] = v.b0[v.active] .- mu.*v.df[v.active]

        # now apply hard threshold on model to enforce original desired sparsity k
        project_k!(v.b,k)
#        selectperm!(v.bidxs, v.b, 1:k, by=abs, rev=true, initialized=true)
        v.bk = v.b[view(v.bidxs, 1:k)]
        fill!(v.b, 0)
        v.b[view(v.bidxs, 1:k)] = v.bk

        # check for convergence
        converged_obj  = abs(loss - loss0) < tol
        converged_grad = normdf < tolG
        converged      = converged_obj || converged_grad
        stuck          = abs(loss - loss00) < tol

        # output algorithm progress
        quiet || @printf("%d\t%d\t%d\t%3.7f\t%3.7f\n",iter,nt_iter,bktrk,loss,normdf)

        # check for convergence
        # if converged and in feasible set, then algorithm converged before maximum iteration
        # perform final computations and output return variables
        if converged

            # send elements below tol to zero
            threshold!(v.b, tol)

            # if requested, apply final refit without regularization
            # note that regularization is NOT set exactly to zero!
            # this aids numerical stability
            if refit
                copy!(v.idxs0, v.idxs)
                update_indices!(v.idxs, v.b)
                decompress_genotypes!(v.xk, x, v.idxs, mask_n)
                try
                    nt_iter, bktrk = fit_logistic!(v, y, 2*eps(T), tol=tolrefit, max_iter=max_step, quiet=true)
                    v.b[v.idxs] = v.bk
                catch e
                    nt_iter = 0
                    bktrk   = 0
                end
            end

            # stop time
            exec_time = toq()

            # announce convergence
            !quiet && print_log_convergence(iter, loss, exec_time, normdf)

            # these are output variables for function
            # wrap them into a Dict and return
            return IHTLogResults(exec_time, iter, loss, copy(v.b), copy(v.active))

        end # end convergence check

        # unconverged at this point
        # refit nonzeroes in b?
        if refit

            # update boolean vector of nonzeroes
            copy!(v.idxs0, v.idxs)
            update_indices!(v.idxs, v.b)

            # update active set of x, if necessary
            (iter == 1 || !isequal(v.idxs, v.idxs0)) && decompress_genotypes!(v.xk, x, v.idxs, mask_n)

            # attempt refit but guard against possible instabilities or singularities
            # these are more frequent if lambda is small
            # if refitting destabilizes, then leave b alone
            try
                nt_iter, bktrk = fit_logistic!(v, y, lambda, tol=tolrefit, max_iter=max_step, quiet=true)
                v.b[v.idxs] = v.bk
            catch e
#                warn("in refitting, caught error: ", e)
#                warn("skipping refit")
                bktrk   = 0
                nt_iter = 0
            end # end try-catch for refit
        end # end if-else for refit
    end # end main GraSP iterations

    # return null result
    return IHTLogResults(zero(T), 0, zero(T), zeros(T, 1), zeros(Int, 1))
end # end L0_log

"""
    iht_path_log(x::BEDFile, y, path)

If used with a `BEDFile` object `x`, then the additional optional arguments are:

- `pids`, a vector of process IDs. Defaults to `procs()`.
"""
function iht_path_log{T <: Float}(
    x        :: BEDFile,
    y        :: SharedVector{T},
    path     :: DenseVector{Int};
    pids     :: Vector{Int}    = procs(x),
    mask_n   :: BitArray{1}    = trues(y),
    lambdas  :: DenseVector{T} = default_lambda(x, y, length(path)),
    mu       :: T    = one(T),
    tol      :: T    = convert(T, 1e-6),
    tolG     :: T    = convert(T, 1e-3),
    tolrefit :: T    = convert(T, 1e-6),
    max_iter :: Int  = 100,
    max_step :: Int  = 100,
    refit    :: Bool = true,
    quiet    :: Bool = true,
)

    # dimensions of problems?
    n,p = size(x)

    # how many models will we compute?
    nmodels = length(path)

    # preallocate temporary arrays
    v = IHTLogVariables(x, y, 1)

    # betas will be a sparse matrix with models
    betas = spzeros(T, p, nmodels)

    # compute the path
    for i = 1:nmodels

        # model size?
        q = path[i]
        quiet || println("Current model size: $q")

        # update temporary variables with new threshold size q
        update_variables!(v, x, q)

        # store projection of beta onto largest k nonzeroes in magnitude
        project_k!(v.b, q)

        # current regularization parameter?
        λ = lambdas[i]

        # now compute current model
        output = L0_log(x, y, q, v=v, lambda=λ, mu=mu, tol=tol, tolG=tolG, tolrefit=tolrefit, max_iter=max_iter, max_step=max_step, refit=refit, quiet=quiet)

        # extract and save model
        betas[:,i] = sparsevec(sdata(output.beta))

        # run garbage compiler to clear temporary arrays for next step of path
        #gc()
    end

    # return a sparsified copy of the models
    return betas
end

function one_fold_log{T <: Float}(
    x         :: BEDFile{T},
    y         :: DenseVector{T},
    path      :: DenseVector{Int},
    folds     :: DenseVector{Int},
    fold      :: Int;
    lambdas   :: DenseVector{T} = default_lambda(x, y, length(path)),
    pids      :: Vector{Int}    = procs(x),
    criterion :: String = "deviance",
    tol       :: T      = convert(T, 1e-6),
    tolG      :: T      = convert(T, 1e-3),
    tolrefit  :: T      = convert(T, 1e-6),
    max_iter  :: Int    = 100,
    max_step  :: Int    = 100,
    refit     :: Bool   = true,
    quiet     :: Bool   = true,
)
    # problem dimensions?
    n,p = size(x)
    nmodels = length(path)

    # ensure that crossvalidation criterion is valid
    criterion in ["deviance", "class"] || throw(ArgumentError("Argument criterion must be 'deviance' or 'class'"))

    # make vector of indices for folds
    test_idx  = folds .== fold
    test_size = sum(test_idx)

    # train_idx is the vector that indexes the TRAINING set
    train_idx = !test_idx

    # allocate the arrays for the training set
#    train_idx = convert(Vector{Int}, train_idx)
#    test_idx  = convert(Vector{Int}, train_idx)

    # compute the regularization path on the training set
    betas = iht_path_log(x, y, path, tol=tol, tolG=tolG, tolrefit=tolrefit, max_iter=max_iter, max_step=max_step, refit=refit, quiet=quiet, lambdas=lambdas, mask_n=train_idx)

    # allocate an index vector for b
    idx = falses(p)

    # allocate temporary arrays for the test set
    xb  = SharedArray(T, n, pids=pids) :: SharedVector{T}
    b   = SharedArray(T, p, pids=pids) :: SharedVector{T}
    lxb = zeros(T, (n, nmodels))
    xβ  = zeros(T, (n, nmodels))

    # compute the mean out-of-sample error for the TEST set
    # do this for every computed model in regularization path
    for i = 1:size(betas,2)

        # which model are we analyzing?
        m = path[i]

        # pull ith model in dense vector format
        b2 = full(vec(betas[:,i]))

        # copy it into SharedArray b
        copy!(b,b2)

        # indices stores Boolean indexes of nonzeroes in b
        update_indices!(idx, b)

        # compute estimated responses
        A_mul_B!(xb, x, b, idx, m, test_idx, pids=pids)
        setindex!(xβ, xb, :, m)
    end

    # mask data from training set
    fill!(view(xb, train_idx, :), zero(T))
    logistic!(lxb, xb)
    fill!(view(lxb, train_idx, :), zero(T))
    ytest = copy(y)
    fill!(view(ytest, train_idx, :), zero(T))

    # compute out-of-sample error as misclassification (L1 norm) averaged over size of test set
    errors = zeros(T, nmodels)
    if criterion == "class" # compute misclassification error
        errors .= vec(sumabs(lxb .- ytest, 1)) ./ sum(!train_idx)
    else # else -> criterion == "deviance" = 2*(neg logistic loss with no reg)
        errors .= - 2 ./ test_size .* ( vec( sum(xβ .* ytest, 1) .- sum(log(one(T) .+ exp(xβ)), 1) ) )
    end

    return errors :: Vector{T}
end

"""
    pfold_log(xfile, xtfile, x2file, yfile, meanfile, precfile, path, folds, q[, pids=procs()])

This function is the parallel execution kernel in `cv_log()`. It is not meant to be called outside of `cv_iht()`.
It will distribute `q` crossvalidation folds across the processes supplied by the optional argument `pids` and call `one_fold_log()` for each fold.
Each fold will compute a regularization path given by `path`.
`pfold()` collects the vectors of MSEs returned by calling `one_fold_log()` for each process, reduces them, and returns their average across all folds.
"""
function pfold_log(
    T        :: Type,
    xfile    :: String,
    xtfile   :: String,
    x2file   :: String,
    yfile    :: String,
    meanfile :: String,
    precfile :: String;
    n        :: Int = begin
                  # find n from the corresponding FAM file, then make folds
                  famfile = xfile[1:(endof(xfile)-3)] * "fam"
                  countlines(famfile)
                  end,
    p        :: Int = begin
                  # find p from the corresponding BIM file, then make path
                  bimfile = xfile[1:(endof(xfile)-3)] * "bim"
                  countlines(bimfile)
                end,
    path      :: DenseVector{Int} = collect(1:min(20,p)),
    folds     :: DenseVector{Int} = cv_get_folds(n, q)
    q         :: Int = cv_get_num_folds(3,5),
    pids      :: Vector{Int} = procs(),
    lambdas   :: DenseVector{Float} = default_lambda(n, p, q),
    criterion :: String = "deviance",
    tol       :: T      = convert(T, 1e-6),
    tolG      :: T      = convert(T, 1e-3),
    tolrefit  :: T      = convert(T, 1e-6),
    max_iter  :: Int    = 100,
    max_step  :: Int    = 100,
    quiet     :: Bool   = true,
    refit     :: Bool   = true,
    header    :: Bool   = false
)
    # ensure that crossvalidation criterion is valid
    criterion in ["deviance", "class"] || throw(ArgumentError("Argument criterion must be 'deviance' or 'class'"))
    T <: Float || throw(ArgumentError("Argument T must be either Float32 or T"))

    # how many CPU processes can pfold use?
    np = length(pids)

    # report on CPU processes
    quiet || println("pfold: np = ", np)
    quiet || println("pids = ", pids)

    # set up function to share state (indices of folds)
    i = 1
    nextidx() = (idx=i; i+=1; idx)

    # preallocate cell array for results
    results = SharedArray(T, (length(path),q), pids=pids) :: SharedMatrix{T}

    # master process will distribute tasks to workers
    # master synchronizes results at end before returning
    @sync begin

        # loop over all workers
        for worker in pids

            # exclude process that launched pfold, unless only one process is available
            if worker != myid() || np == 1

                # asynchronously distribute tasks
                @async begin
                    while true

                        # grab next fold
                        current_fold = nextidx()

                        # if current fold exceeds total number of folds then exit loop
                        current_fold > q && break

                        # report distribution of fold to worker and device
                        quiet || print_with_color(:blue, "Computing fold $current_fold on worker $worker.\n\n")

                        # launch job on worker
                        # worker loads data from file paths and then computes the errors in one fold
                        r = remotecall_fetch(worker) do
                            processes = [worker]
                            x = BEDFile(T, xfile, xtfile, x2file, meanfile, precfile, pids=processes, header=header)
                            n = x.geno.n
                            y = SharedArray(abspath(yfile), T, (n,), pids=processes) :: SharedVector{T}
                            one_fold_log(x, y, path, folds, current_fold, max_iter=max_iter, max_step=max_step, quiet=quiet, pids=processes, tol=tol, tolG=tolG, tolrefit=tolrefit, refit=refit, criterion=criterion, lambdas=lambdas)

                        end # end remotecall_fetch()
                        setindex!(results, r, :, current_fold)
                    end # end while
                end # end @async
            end # end if
        end # end for
    end # end @sync

    # return reduction (row-wise sum) over results
    return (vec(sum(results, 2) ./ q)) :: Vector{T}
end

# default type for pfold_log is Float64
pfold_log(xfile::String, xtfile::String, x2file::String, yfile::String, meanfile::String, precfile::String; n::Int=begin famfile=xfile[1:(endof(xfile)-3)]*"fam"; countlines(famfile) end, p::Int=begin bimfile=xfile[1:(endof(xfile)-3)]*"bim"; countlines(bimfile) end, q::Int = cv_get_num_folds(3,5), path::DenseVector{Int} = collect(1:min(20,p)), folds::DenseVector{Int} = cv_get_folds(n,q), lambdas::DenseVector{Float64} = default_lambda(n,p,q), pids::Vector{Int}=procs(), criterion::String="deviance", tol::Float64=1e-6, tolG::Float64=1e-3, tolrefit::Float64=1e-6, max_iter::Int=100, max_step::Int=100, quiet::Bool=true, refit::Bool=true, header::Bool=false) = pfold_log(Float64, xfile, xtfile, x2file, yfile, meanfile, precfile, path, folds, q, n=n, p=p, pids=pids, lambdas=lambdas, criterion=criterion, tol=tol, tolG=tolG, tolrefit=tolrefit, max_iter=max_iter, max_step=max_step, quiet=quiet, refit=refit, header=header)

"""
    cv_log(xfile,xtfile,x2file,yfile,meanfile,precfile,path,kernfile,folds,q [, pids=procs()])

This variant of `cv_log()` performs `q` crossvalidation with a `BEDFile` object loaded by `xfile`, `xtfile`, and `x2file`,
with column means stored in `meanfile` and column precisions stored in `precfile`.
The continuous response is stored in `yfile` with data particioned by the `Int` vector `folds`.
The folds are distributed across the processes given by `pids`.
"""
function cv_log(
    T         :: Type,
    xfile     :: String,
    xtfile    :: String,
    x2file    :: String,
    yfile     :: String,
    meanfile  :: String,
    precfile  :: String;
    n         :: Int = begin
                  # find n from the corresponding FAM file, then make folds
                  famfile = xfile[1:(endof(xfile)-3)] * "fam"
                  countlines(famfile)
                 end,
    p         :: Int = begin
                  # find p from the corresponding BIM file, then make path
                  bimfile = xfile[1:(endof(xfile)-3)] * "bim"
                  countlines(bimfile)
                 end,
    q         :: Int = cv_get_num_folds(3,5),
    path      :: DenseVector{Int} = collect(1:min(20,p)),
    folds     :: DenseVector{Int} = cv_get_folds(n,q),
    pids      :: DenseVector{Int} = procs(),
    lambdas   :: DenseVector      = default_lambda(n,p,q), # type unstable?
    criterion :: String = "deviance",
    tol       :: Float  = convert(T, 1e-6),
    tolG      :: Float  = convert(T, 1e-3),
    tolrefit  :: Float  = convert(T, 1e-6),
    max_iter  :: Int    = 100,
    max_step  :: Int    = 100,
    quiet     :: Bool   = true,
    refit     :: Bool   = true,
    header    :: Bool   = false
)
    # ensure that crossvalidation criterion is valid
    criterion in ["deviance", "class"] || throw(ArgumentError("Argument criterion must be 'deviance' or 'class'"))
    T <: Float || throw(ArgumentError("Argument T must be either Float32 or T"))

    # how many elements are in the path?
    nmodels = length(path)

    # want to compute a path for each fold
    # the folds are computed asynchronously
    # only use the worker processes
    errors = pfold_log(T, xfile, xtfile, x2file, yfile, meanfile, precfile, n=n, p=p, path=path, folds=folds, q=q, max_iter=max_iter, max_step=max_step, quiet=quiet, pids=pids, header=header, tolG=tolG, tolrefit=tolrefit, refit=refit, criterion=criterion, lambdas=lambdas)

    # what is the best model size?
    k = path[indmin(errors)] :: Int

    # print results
    !quiet && print_cv_results(errors, path, k)

    # recompute ideal model
    # load data on *all* processes
    x = BEDFile(T, xfile, xtfile, x2file, header=header, pids=pids)
    n = x.geno.n
    y = SharedArray(abspath(yfile), T, (n,), pids=pids) :: SharedVector{T}

    # get lambda value for best model
    λ = lambdas[path .== k][1]

    # use L0_reg to extract model
    output = L0_log(x, y, k, tol=tol, max_iter=max_iter, max_step=max_step, quiet=quiet, tolG=tolG, tolrefit=tolrefit, refit=refit, lambda=λ, pids=pids)

    # which components of beta are nonzero?
    bidx = find(output.beta)

    return IHTCrossvalidationResults(errors, sdata(path), output.beta[bidx], bidx, k)

end

# default XV function type is Float64
cv_log(xfile::String, xtfile::String, x2file::String, yfile::String, meanfile::String, precfile::String; n::Int = begin famfile=xfile[1:(endof(xfile)-3)]*"fam"; countlines(famfile) end, p::Int = begin bimfile=xfile[1:(endof(xfile)-3)]*"bim"; countlines(bimfile) end, q::Int=cv_get_num_folds(3,5), path::DenseVector{Int} = collect(1:min(20,p)), folds::DenseVector{Int} = cv_get_folds(n,q), pids::Vector{Int} = procs(), lambdas::DenseVector{Float64} = default_lambda(n,p,q),  criterion::String="deviance", tol::Float64=1e-6, tolG::Float64=1e-3, tolrefit::Float64=1e-6, max_iter::Int=100, max_step::Int=100, quiet::Bool=true, refit::Bool=true, header::Bool=false)=cv_log(Float64, xfile, xtfile, x2file, yfile, meanfile, precfile, path=path, folds=folds, q=q, pids=pids, lambdas=lambdas, criterion=criterion, tol=tol, tolG=tolG, tolrefit=tolrefit, max_iter=max_iter, max_step=max_step, quiet=quiet, refit=refit, header=header)


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -17.39516631142554 <= q <= 12.050938789121952)

                     add_NL_constraint(m, :(sigmoid(0.1892165831624366 + -0.6797303891528506*sigmoid(0.23382152176968107 + 0.970568112692761*$(x[1]) + 0.9147961292868287*$(x[2]) + 0.25099364821370385*$(x[3])) + -0.5563706493225857*sigmoid(0.1824369939040409 + -0.636682724457545*$(x[1]) + 0.684450108697257*$(x[2]) + -0.3896097552087272*$(x[3])) + -0.5881166064998418*sigmoid(-0.7017137443246066 + -0.579761469455264*$(x[1]) + -0.38090529589576505*$(x[2]) + -0.30474275998045197*$(x[3])) + -0.5324749908347566*sigmoid(0.495407339331571 + -0.11045663155169017*$(x[1]) + -0.19556503402268088*$(x[2]) + 0.8751710578699647*$(x[3]))) + sigmoid(-0.7087227741366906 + -0.3339546766337098*sigmoid(0.23382152176968107 + 0.970568112692761*$(x[1]) + 0.9147961292868287*$(x[2]) + 0.25099364821370385*$(x[3])) + -0.8987724666577948*sigmoid(0.1824369939040409 + -0.636682724457545*$(x[1]) + 0.684450108697257*$(x[2]) + -0.3896097552087272*$(x[3])) + 0.3907505378456695*sigmoid(-0.7017137443246066 + -0.579761469455264*$(x[1]) + -0.38090529589576505*$(x[2]) + -0.30474275998045197*$(x[3])) + -0.7983638213649877*sigmoid(0.495407339331571 + -0.11045663155169017*$(x[1]) + -0.19556503402268088*$(x[2]) + 0.8751710578699647*$(x[3]))) + sigmoid(0.16171979073463483 + 0.31171990946083294*sigmoid(0.23382152176968107 + 0.970568112692761*$(x[1]) + 0.9147961292868287*$(x[2]) + 0.25099364821370385*$(x[3])) + 0.2757921127232579*sigmoid(0.1824369939040409 + -0.636682724457545*$(x[1]) + 0.684450108697257*$(x[2]) + -0.3896097552087272*$(x[3])) + -0.732639957443523*sigmoid(-0.7017137443246066 + -0.579761469455264*$(x[1]) + -0.38090529589576505*$(x[2]) + -0.30474275998045197*$(x[3])) + -0.947011017876227*sigmoid(0.495407339331571 + -0.11045663155169017*$(x[1]) + -0.19556503402268088*$(x[2]) + 0.8751710578699647*$(x[3]))) + sigmoid(-0.17479526548317725 + -0.7121643152938684*sigmoid(0.23382152176968107 + 0.970568112692761*$(x[1]) + 0.9147961292868287*$(x[2]) + 0.25099364821370385*$(x[3])) + -0.20573796119512977*sigmoid(0.1824369939040409 + -0.636682724457545*$(x[1]) + 0.684450108697257*$(x[2]) + -0.3896097552087272*$(x[3])) + 0.920351043144791*sigmoid(-0.7017137443246066 + -0.579761469455264*$(x[1]) + -0.38090529589576505*$(x[2]) + -0.30474275998045197*$(x[3])) + -0.8770524767395065*sigmoid(0.495407339331571 + -0.11045663155169017*$(x[1]) + -0.19556503402268088*$(x[2]) + 0.8751710578699647*$(x[3]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# Tests for Profit Reverse DDF DEA Models
@testset "ProfitReverseDDFDEAModel" begin

    # Test using Book data
    X = [2; 4; 8; 12; 6; 14; 14; 9.412];
    Y = [1; 5; 8; 9; 3; 7; 9; 2.353];
    W = [1; 1; 1; 1; 1; 1; 1; 1];
    P = [2; 2; 2; 2; 2; 2; 2; 2];

    # Profit Reverse DDF :ERG
    profitrddferg = deaprofitrddf(X, Y, W, P, :ERG)

    @test typeof(profitrddferg) == ProfitReverseDDFDEAModel

    @test nobs(profitrddferg) == 8
    @test ninputs(profitrddferg) == 1
    @test noutputs(profitrddferg) == 1
    @test ismonetary(profitrddferg) == false

    @test profitrddferg.measure == :ERG
    @test efficiency(profitrddferg) == efficiency(profitrddferg, :Economic)
    @test efficiency(profitrddferg, :Economic)   ≈ [4.0; 0.5; 0.0; 1/6; 0.8; 4/7; 2/7; 0.949449] atol = 1e-5
    @test efficiency(profitrddferg, :Technical)  ≈ [0.0; 0.0; 0.0; 0.0; 0.6; 11/21; 1/7; 0.8] atol = 1e-5
    @test efficiency(profitrddferg, :Allocative) ≈ [4.0; 0.5; 0.0; 1/6; 0.2; 1/21; 1/7; 0.149449] atol = 1e-5
    @test normfactor(profitrddferg) ≈ [2.0; 4.0; 8.0; 12.0; 10.0; 14.0; 14.0; 13.3825]

    # With monetary = true
    profitrddfergmon = deaprofitrddf(X, Y, W, P, :ERG, monetary = true)

    @test ismonetary(profitrddfergmon) == true

    @test efficiency(profitrddfergmon) ≈ efficiency(profitrddferg) .* normfactor(profitrddferg)
    @test efficiency(profitrddfergmon, :Technical) ≈ efficiency(profitrddferg, :Technical) .* normfactor(profitrddferg)
    @test efficiency(profitrddfergmon, :Allocative) ≈ efficiency(profitrddferg, :Allocative) .* normfactor(profitrddferg)

    # Print
    show(IOBuffer(), profitrddferg)

    # Profit Reverse DDF :MDDF
    profitrddfmddf = deaprofitrddf(X, Y, W, P, :MDDF, Gx = :Observed, Gy = :Observed)

    @test typeof(profitrddfmddf) == ProfitReverseDDFDEAModel

    @test profitrddfmddf.measure == :MDDF
    @test efficiency(profitrddfmddf) == efficiency(profitrddfmddf, :Economic)
    @test efficiency(profitrddfmddf, :Economic)   ≈ [4.0; 0.5; 0.0; 1/6; 4/3; 4/7; 2/7; 2.699958] atol = 1e-5
    @test efficiency(profitrddfmddf, :Technical)  ≈ [0.0; 0.0; 0.0; 0.0; 7/6; 4/7; 1/7; 2.549936] atol = 1e-5
    @test efficiency(profitrddfmddf, :Allocative) ≈ [4.0; 0.5; 0.0; 1/6; 1/6; 0.0; 1/7; 0.150021] atol = 1e-5
    @test normfactor(profitrddfmddf) ≈ [2.0; 4.0; 8.0; 12.0; 6.0; 14.0; 14.0; 4.706] atol = 1e-5

    # Profit Reverse DDF :MDDF with Custom Direction
    profitrddfmddfcurstom = deaprofitrddf(X, Y, W, P, :MDDF, Gx = X, Gy = Y)
    @test efficiency(profitrddfmddfcurstom) == efficiency(profitrddfmddf)

    # Test errors
    @test_throws DimensionMismatch deaprofitrddf([1; 2 ; 3], [4; 5], [1; 1; 1], [1; 1], :ERG) #  Different number of observations
    @test_throws DimensionMismatch deaprofitrddf([1; 2; 3], [4; 5; 6], [1; 2; 3; 4], [1; 1; 1], :ERG) # Different number of observation in input prices
    @test_throws DimensionMismatch deaprofitrddf([1; 2; 3], [4; 5; 6], [1; 1; 1], [1; 2; 3; 4], :ERG) # Different number of observation in output prices
    @test_throws DimensionMismatch deaprofitrddf([1 1; 2 2; 3 3 ], [4; 5; 6], [1 1 1; 2 2 2; 3 3 3], [1; 1; 1], :ERG) # Different number of input prices and inputs
    @test_throws DimensionMismatch deaprofitrddf([1; 2; 3 ], [4; 5; 6], [1; 1; 1], [1 1 1; 2 2 2; 3 3 3], :ERG) # Different number of output prices and outputs
    @test_throws ArgumentError deaprofitrddf([1; 2; 3], [4; 5; 6], [1; 2; 3], [4; 5; 6], :Error) # Invalid associated efficiency measure
    @test_throws ArgumentError deaprofitrddf([1; 2; 3], [4; 5; 6], [1; 2; 3], [4; 5; 6], :MDDF, Gx = :Error, Gy = :Observed) # Invalid Gx direction
    @test_throws ArgumentError deaprofitrddf([1; 2; 3], [4; 5; 6], [1; 2; 3], [4; 5; 6], :MDDF, Gx = :Observed, Gy = :Error) # Invalid Gy direction

end


text:
using Test
using NeXLSpectrum

include("window.jl")
include("detefficiency.jl")
include("spectrum.jl")
include("emsa.jl")
include("brukerspx.jl")
include("filterfit.jl")
include("hyperspectrum.jl")
include("line.jl")
include("multikev.jl")
include("standardize.jl")
include("multidet.jl")


text:
# Autogenerated wrapper script for CompilerSupportLibraries_jll for armv7l-linux-gnueabihf-libgfortran5
export libatomic, libgcc_s, libgfortran, libgomp, libstdcxx

JLLWrappers.@generate_wrapper_header("CompilerSupportLibraries")
JLLWrappers.@declare_library_product(libatomic, "libatomic.so.1")
JLLWrappers.@declare_library_product(libgcc_s, "libgcc_s.so.1")
JLLWrappers.@declare_library_product(libgfortran, "libgfortran.so.5")
JLLWrappers.@declare_library_product(libgomp, "libgomp.so.1")
JLLWrappers.@declare_library_product(libstdcxx, "libstdc++.so.6")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libatomic,
        "lib/libatomic.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgcc_s,
        "lib/libgcc_s.so.1",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgfortran,
        "lib/libgfortran.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgomp,
        "lib/libgomp.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libstdcxx,
        "lib/libstdc++.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@testset "learning transition small functions" begin
    for x0 ∈ (0.25, 0.5, 0.75,)
        vspace = range(-2.0, stop=2.0, length=5)
        @test Calculus.derivative(_ρ, x0) ≈ ShaleDrillingModel._dρdθρ(x0)
        @test Calculus.derivative(_ρ, x0) ≈ ShaleDrillingModel._dρdσ(x0)

        for u in vspace
            for v in vspace
                @test Calculus.derivative((z) -> _ψ1(u,v,_ρ(z)), x0) ≈ _dψ1dθρ(u,v,_ρ(x0), x0)
                @test Calculus.derivative((z) -> _ψ1(u,v,z), x0) ≈ ShaleDrillingModel._dψ1dρ(u,v,x0)
            end
        end
    end
end

@testset "Frechet derivative of ψ transition matrix" begin

    dpi = similar(tmpv.Πψtmp)
    fdpi = similar(tmpv.Πψtmp)
    h = peturb(σv)
    σ1 = σv - h
    σ2 = σv + h
    hh = σ2 - σ1

    ShaleDrillingModel._βΠψ!(dpi, ψspace, σ1, β)
    ShaleDrillingModel._βΠψ!(fdpi, ψspace, σ2, β)
    fdpi .-= dpi
    fdpi ./= hh
    ShaleDrillingModel._βΠψdθρ!(dpi, ψspace, σv, β)

    @views maxv, idx = findmax(abs.(dpi .- fdpi))
    sub = CartesianIndices(fdpi)[idx]
    @show "worst value is $maxv at $sub for β dΠ/dσ"

    @test 0.0 < maxv < 1e-5
    @test fdpi ≈ dpi

end


text:
# target a univariate measure μ(dx) ∝ (x-\mu)^2/2\sigma^2 + δ_0 (x) dx
using ZigZagBoomerang
const ZZB = ZigZagBoomerang
function freezing_time(x, θ)
    if θ*x > 0
        return Inf
    else
        return -x/θ
    end
end

# k determines the weight on the Dirac measure. The smaller k, the higher the wegiht
function ss_pdmp(∇ϕ, x, θ, T, c, k, Flow::ZZB.ContinuousDynamics; adapt=false, factor=2.0)
    t = zero(T)
    Ξ = [(t, x, θ)]
    num = acc = 0
    a, b = ZZB.ab(x, θ, c, Flow)
    t_ref = t + ZZB.waiting_time_ref(Flow)
    t′ =  t + poisson_time(a, b, rand())
    tˣ = t + freezing_time(x, θ)
    while t < T
        if  tˣ < min(t_ref, t′)
            t, x, θ = ZZB.move_forward(tˣ - t, t, x, θ, Flow) # go to 0
            @assert -0.0001 < x < 0.0001 #check
            #t, x , θ = tˣ, 0.0, θ #go to 0
            push!(Ξ, (t, x, 0.0))
            t = t - log(rand())/k #wait exponential time
            push!(Ξ, (t, x, θ))
            tˣ  = Inf
        elseif t_ref < t′
            t, x, θ = ZZB.move_forward(t_ref - t, t, x, θ, Flow)
            θ = sqrt(Flow.Σ)*randn()
            t_ref = t +  ZZB.waiting_time_ref(Flow)
            tˣ = t + freezing_time(x, θ)
            push!(Ξ, (t, x, θ))
        else
            τ = t′ - t
            t, x, θ = ZZB.move_forward(τ, t, x, θ, Flow)
            l, lb = ZZB.λ(∇ϕ, x, θ, Flow), ZZB.λ_bar(τ, a, b)
            num += 1
            if rand()*lb < l
                acc += 1
                if l >= lb
                    !adapt && error("Tuning parameter `c` too small.")
                    c *= factor
                end
                θ = -θ  # In multi dimensions the change of velocity is different:
                        # reflection symmetric on the normal vector of the contour
                tˣ = t + freezing_time(x, θ)
                push!(Ξ, (t, x, θ))
            end
        end
        a, b = ZZB.ab(x, θ, c, Flow)
        t′ = t + poisson_time(a, b, rand())
    end
    return Ξ, acc/num
end

ϕ(x) = x^2/2 #not used
∇ϕ(x) = x

# Example: ZigZag
x0, θ0 = -1.0, 1.0
T = 300.0
c = 1.0
k = 1.0
out1, acc = ss_pdmp(∇ϕ, x0, θ0, T, c, k, ZigZag1d())

# Boomerang
function freezing_time(x, θ)
    if θ*x > 0
        return Inf
    else
        return -x/θ
    end
end
out2, acc = ss_pdmp(∇ϕ, x0, θ0, T, c, k, Boomerang1d(0.1))
@show acc
using Makie, CairoMakie
p1 = lines(eventtime.(out1), eventposition.(out1))

## TEST probability of x being 0 euqal to total time divided by time spent on 0
function prob_of_zero(out1)
    k = 0
    freeze =false
    for event in out1
        if freeze == true
            k += event[1]
        end
        if event[3] == 0.0
            freeze = true
            k -= event[1]
        else
            freeze = false
        end
    end
    k/out1[end][1]
end

const μ = 1.0
ϕ(x) = (x-μ)^2/2 #not used
∇ϕ(x) = x - μ
# Example: ZigZag
x0, θ0 = -1.0, 1.0
T = 2000.0
c = 10.01
k = 0.1
@show acc
out1, acc = ss_pdmp(∇ϕ, x0, θ0, T, c, k, ZigZag1d())
using Makie, CairoMakie
p1 = lines(eventtime.(out1), eventposition.(out1), linewidth = 0.3)
prob_of_zero(out1)

f0 = 1/(sqrt(2*pi))*exp(-1/2*μ^2)
1/(k/f0 + 1)


text:
using CooperationPolarization
using BenchmarkTools
using Profile
using PyPlot: plt, savefig
using Dates

const b = 200. # benefit to cooperating
const c = 0.2 # cost to cooperating

const N = 40 # population size
const M = 1 # number of sets
const K = 1 # number of issues to care about
const P = 2 # number of parties

const β = 0.001    # selection strength
const u = 0.001    # strategy mutation rate
const v = 0.001    # set mutation rate
const p = 1.       # party bias (0 = low, 1 = high)
const ϵ = 1.       # damping factor (0 = no bias in mutation, 1 = max bias (= p) in mutation)
# note: α = 0 automatically

# simulation parameters
const verbose       = false
const verbose_track = false
const generations   = 10000000
const testtype      = "test"

# plotting parameters
const gen_skip  = 1 # max( div(generations, 10000), 1)
const gen_start = max( 1, generations-1000000 )
const gen_end   = generations

const start_time = Dates.format(now(),"yymmdd-HHMMSS.sss")

println("parameters defined")

sets    = random_sets(N, M, K, P)
game    = Game(b, c, β, u, v, p, ϵ)
pop     = Population(sets, game, verbose)
tracker = Tracker(pop, generations, verbose_track)

println("Payoff matrix: \t $(game.A)")
println("Payoffs:\t $(pop.payoffs)")
println("Verbose?:\t $(pop.verbose)")
println("Initial:\t strategies_changed = $(pop.strategies_changed), sets_changed = $(pop.sets_changed)")

@time for g::Int64 in 1:generations
    if generations > 10 && (g-1) % (generations ÷ 10) == 0  # set up to print 100 generations
        println("\tinitiating generation $(tracker.pop.generation)")
    end
    # update population
    evolve!(tracker.pop, 1)  # this is the same as evolve!(pop, 1)
    track!(tracker, g)

    # check pop and tracker.pop are the same
    # println("\t$(pop.sets.h == tracker.pop.sets.h)") # just to check

    # check whether the updating step is working properly
    # if true, there should be a mix of 'true' and 'false' for sets_chagned and strategies_changed
    # println("\tsets_changed = $(pop.sets_changed), \tstrat_changed = $(pop.strategies_changed)")

    if g == generations
        println("replicate complete, ran $(generations) generations:\n")
        println("\tnum_imitating:\t\t\t$(tracker.pop.num_imitating)")
        println("\tnum_not_imitating:\t\t$(tracker.pop.num_not_imitating)")
        println("\tnum_strategies_changed:  \t$(tracker.pop.num_strategies_changed)")
        println("\t\tnum_random_strategies: \t$(tracker.pop.num_random_strategies)")
        println("\tnum_sets_changed: \t\t$(tracker.pop.num_sets_changed)")
        println("\t\tnum_random_opinions: \t$(tracker.pop.num_random_opinions)")
        println("\t\tnum_biased_opinions: \t$(tracker.pop.num_biased_opinions)")
        println("\tnum_both_changed: \t\t$(tracker.pop.num_both_changed)\n")
    end
end

# plotting
fig = plotter( tracker, gen_start, gen_skip, gen_end )

plt.suptitle("Parameters: b=$b, c=$c, N=$N, M=$M, K=$K, β=$β, u=$u, v=$v, p=$p, ϵ=$ϵ", fontsize=12)
fig.tight_layout(rect=[0, 0.03, 1, 0.97])

display(fig)

# save file
# savefig("plots/$(testtype)_$(filename).pdf")