text:
################################################################################
# This set of functions provides various general utilities.
################################################################################
#
# Required external packages.
#
# using Distributions
# using LinearAlgebra
# using Statistics

export blanks, empties, repeated_string, select_set_element, random_category
export standarize!, sample_mean_std, sample_stats, print_sample_stats
export simes_fdr, fast_regress, fast_score_test
export glm_regress, glm_score_statistic

"""
Creates an array of blank strings.
"""
function blanks(n::Int)
#
  blank = Array{String}(undef, n)
  for i = 1:n
    blank[i] = ""
  end
  return blank
end # function blanks

"""
Creates an array of empty ordered sets.
"""
function empties(n::Int)
#
  empty = Array{BitSet}(undef, n)
  for i = 1:n
    empty[i] = BitSet()
  end
  return empty
end # function empties

"""
Identifies a repeated string in an array of strings.
"""
function repeated_string(name::Array{String, 1})
#
  sorted = sort(name)
  for i = 2:length(name)
    if sorted[i - 1] == sorted[i]
      return (true, sorted[i])
    end
  end
  return (false, "")
end # function repeated_string

"""
Selects element number k of the integer set S.
"""
function select_set_element(set_s::BitSet, k::Int)
#
  j = 1
  for s in set_s
    if j == k
      return s
    else
      j = j + 1
    end
  end
  return nothing
end # function select_set_element

"""
Returns a random category (numbered 1, 2, and so forth)
according to a vector of category frequencies.
"""
function random_category(frequency::Vector{T} where T <: Real)
#
  u = rand(1)
  for i = 1:length(frequency)
    if u[1] <= frequency[i]
      return i
    else
      u[1] = u[1] - frequency[i]
    end
  end
  return length(frequency)
end # function random_category

"""
Standardizes the vector x to have mean 0.0 and variance 1.0.
Missing values are ignored.
"""
function standarize!(x::Vector)
#
  avg = mean(skipmissing(x))
  stdev = std(skipmissing(x))
  if stdev > zero(eltype(stdev))
    x .= (x .- avg) ./ stdev
  else
    x .= x .- avg
  end
end # function normalize!

"""
Computes the sample mean and standard deviation.
Missing values are ignored.
"""
function sample_mean_std(x::Vector)
#
  return (mean(skipmissing(x)), std(skipmissing(x)))
end # function sample_mean_std

"""
Computes sample statistics for the data vector x.
Missing values are ignored. The output consists of
the number of values present, the number of values missing,
the minimum, the 0.25 quantile, the median, the 0.75 quantile,
the maximum, the mean, the standard deviation, the skewness,
and the excess kurtosis.
"""
function sample_stats(x::Vector)
#
  y = collect(skipmissing(x))
  (p, n) = (length(y), length(x))
  return p, n - p, minimum(y), quantile(y, 0.25), median(y),
    quantile(y, 0.75), maximum(y), mean(y), std(y),
    skewness(y), kurtosis(y)
end # function sample_stats

"""
Prints the sample statistics output by the function sample_stats.
"""
function print_sample_stats(s, io::IO = STDOUT)
#
  println(io, "Summary Stats:")
  println(io, "Values Present:     ", s[1])
  println(io, "Values Absent:      ", s[2])
  println(io, "Minimum:            ", s[3])
  println(io, "1st Quartile:       ", s[4])
  println(io, "Median:             ", s[5])
  println(io, "3rd Quartile:       ", s[6])
  println(io, "Maximum:            ", s[7])
  println(io, "Mean:               ", s[8])
  println(io, "Standard Deviation: ", s[9])
  println(io, "Skewness:           ", s[10])
  println(io, "Excess Kurtosis:    ", s[11])
end # function print_sample_stats

"""
Performs the Simes false discovery rate (FDR) procedure discussed
by Benjamini and Hochberg. All p-values at or below the threshold
are declared significant for the given FDR rate and number of tests.
"""
function simes_fdr(pvalue::Vector{T} where T <: Real,
                   fdr::Vector{T} where T <: Real, tests::Int)
#
  n = length(fdr)
  threshold = zeros(typeof(pvalue[1]), n)
  number_passing = zeros(Int, n)
  perm = sortperm(pvalue)
  k = 1 # previous value of i
  for j = 1:n
    i = 0
    for m = k:length(pvalue)
      i = m
      if tests * pvalue[perm[m]] > m * fdr[j]; break; end
    end
    if i == 1; continue; end
    k = i - 1
    threshold[j] = pvalue[perm[i - 1]]
    number_passing[j] = i - 1
  end
  return (number_passing, threshold)
end # function simes_fdr

"""
Performs either linear, logistic, or Poisson regression with a canonical
link function. X is the design matrix, and y is the response vector.
The parameter estimates and loglikelihood are returned.
"""
function fast_regress(X::Matrix{Float64}, y::Vector{Float64}, model::AbstractString)

  if model != "linear" && model != "logistic" && model != "poisson"
    throw(ArgumentError(
      "The only model choices are linear, logistic, and poisson.\n \n"))
  end
  #
  # Create the score vector, information matrix, estimate, a work
  # vector z, and the loglikelihood.
  #
  (n, p) = size(X)
  @assert n == length(y)
  score = zeros(p)
  information = zeros(p, p)
  estimate = zeros(p)
  z = zeros(n)
  #
  # Handle linear regression separately.
  #
  if model == "linear"
    BLAS.gemv!('N', 1.0, X, estimate, 0.0, z) # z = X * estimate
    BLAS.axpy!(-1.0, y, z) # z = z - y
    score = BLAS.gemv('T', -1.0, X, z) # score = - X' * (z - y)
    information = BLAS.gemm('T', 'N', X, X) # information = X' * X
    estimate = information \ score
    BLAS.gemv!('N', 1.0, X, estimate, 0.0, z) # z = X * estimate
    obj = - 0.5 * n * log(sum(abs2, y - z) / n) - 0.5 * n
    return (estimate, obj)
  end
  #
  # Prepare for logistic and Poisson regression by estimating the
  # intercept.
  #
  if model == "logistic"
    estimate[1] = log(mean(y) / (1.0 - mean(y)))
  elseif model == "poisson"
    estimate[1] = log(mean(y))
  else
    throw(ArgumentError(
      "The only model choices are linear, logistic, and poisson.\n \n"))
  end
  #
  # Initialize the loglikelihood and the convergence criterion.
  #
  v = zeros(p)
  obj = 0.0
  old_obj = 0.0
  epsilon = 1e-6
  #
  #  Estimate parameters by the scoring algorithm.
  #
  for iteration = 1:10
    #
    # Initialize the score and information.
    #
    fill!(score, 0.0)
    fill!(information, 0.0)
    #
    # Compute the score, information, and loglikelihood (obj).
    #
    BLAS.gemv!('N', 1.0, X, estimate, 0.0, z) # z = X * estimate
    clamp!(z, -20.0, 20.0)
    if model == "logistic"
      z = exp.(-z)
      z = 1.0 ./ (1.0 .+ z)
      w = z .* (1.0 .- z)
      BLAS.axpy!(-1.0, y, z) # z = z - y
      score = BLAS.gemv('T', -1.0, X, z) # score = - X' * (z - y)
      w = sqrt.(w)
      lmul!(Diagonal(w), X) # X = diag(w) * X
      information = BLAS.gemm('T', 'N', X, X) # information = X' * W * X
      w = 1.0 ./ w
      lmul!(Diagonal(w), X) # X = diag(w) * X
    elseif model == "poisson"
      z = exp.(z)
      w = copy(z)
      BLAS.axpy!(-1.0, y, z) # z = z - y
      score = BLAS.gemv('T', -1.0, X, z) # score = - X' * (z - y)
      w = sqrt.(w)
      lmul!(Diagonal(w), X) # X = diag(w) * X
      information = BLAS.gemm('T', 'N', X, X) # information = X' * W * X
      w = 1.0 ./ w
      lmul!(Diagonal(w), X) # X = diag(w) * X
    end
    #
    # Compute the scoring increment.
    #
    increment = information \ score
    #
    # Step halve to produce an increase in the loglikelihood.
    #
    steps = -1
    for step_halve = 0:3
      steps = steps + 1
      obj = 0.0
      estimate = estimate + increment
      BLAS.gemv!('N', 1.0, X, estimate, 0.0, z) # z = X * estimate
      clamp!(z, -20.0, 20.0)
      #
      # Compute the loglikelihood under the appropriate model.
      #
      if model == "logistic"
        z = exp.(-z)
        z = 1.0 ./ (1.0 .+ z)
        for i = 1:n
          if y[i] > 0.0
            obj = obj + log(z[i])
          else
            obj = obj + log(1.0 - z[i])
          end
        end
      elseif model == "poisson"
        for i = 1:n
          q = exp(z[i])
          obj = obj + y[i] * z[i] - q
        end
      end
      #
      # Check for an increase in the loglikelihood.
      #
      if old_obj < obj
        break
      else
        estimate = estimate - increment
        increment = 0.5 * increment
      end
    end
    #
    # Check for convergence.
    #
    if iteration > 1 && abs(obj - old_obj) < epsilon * (abs(old_obj) + 1.0)
      return (estimate, obj)
    else
      old_obj = obj
    end
  end
  return (estimate, obj)
end # function fast_regress

"""
Performs a score test for either linear, logistic, or Poisson regression
with a canonical link function. X is the design matrix, y is the
response vector, and estimate is the MLE under the null hypothesis.
"""
function fast_score_test(X::Matrix{Float64}, y::Vector{Float64},
  estimate::Vector{Float64}, model::AbstractString)

  if model != "linear" && model != "logistic" && model != "poisson"
    throw(ArgumentError(
      "The only model choices are linear, logistic, and poisson.\n \n"))
  end
  #
  # Initialize the score vector and information matrix.
  #
  (n, p) = size(X)
  @assert n == length(y)
  @assert p == length(estimate)
  score = zeros(p)
  information = zeros(p, p)
  z = zeros(n)
  #
  # Handle each model separately.
  #
  BLAS.gemv!('N', 1.0, X, estimate, 0.0, z) # z = X * estimate
  if model == "linear"
    BLAS.axpy!(-1.0, y, z) # z = z - y
    var_inv = n / norm(z)^2 # var = residual sum of squares / n
    score = BLAS.gemv('T', -var_inv, X, z) # score = - X' * (z - y) / var
    information = BLAS.gemm('T', 'N', var_inv, X, X) # informat = X' * X / var
  elseif model == "logistic"
    clamp!(z, -20.0, 20.0)
    z = exp.(-z)
    z = 1.0 ./ (1.0 .+ z)
    w = z .* (1.0 .- z)
    BLAS.axpy!(-1.0, y, z) # z = z - y
    score = BLAS.gemv('T', -1.0, X, z) # score = - X' * (z - y)
    w = sqrt.(w)
    lmul!(Diagonal(w), X) # X = diag(w) * X
    information = BLAS.gemm('T', 'N', X, X) # information = X' * W * X
    w = 1.0 ./ w
    lmul!(Diagonal(w), X) # X = diag(w) * X
  elseif model == "poisson"
    clamp!(z, -20.0, 20.0)
    z = exp.(z)
    w = copy(z)
    BLAS.axpy!(-1.0, y, z) # z = z - y
    score = BLAS.gemv('T', -1.0, X, z) # score = - X' * (z - y)
    w = sqrt.(w)
    lmul!(Diagonal(w), X) # X = diag(w) * X
    information = BLAS.gemm('T', 'N', X, X) # information = X' * W * X
    w = 1.0 ./ w
    lmul!(Diagonal(w), X) # X = diag(w) * X
  end
  #
  # Compute the score statistic from the score and information.
  #
  z = information \ score
  score_test = dot(score, z)
  return score_test
end # function fast_score_test

"""
Performs generalized linear regression.
X is the design matrix, y is the response vector,
meanf is the value and derivative of the inverse link,
and varf is the variance function of the mean.
"""
function glm_regress(X::Matrix, y::Vector, meanf::Function, varf::Function,
  loglikelihood::Function)
#
  (n, p) = size(X)
  @assert n == length(y)
  (score, inform, beta) = (zeros(p), zeros(p, p), zeros(p))
  (x, z) = (zeros(p), zeros(n))
  ybar = mean(y)
  for iteration = 1:20 # find the intercept by Newton's method
    g = meanf(beta[1])
    beta[1] = beta[1] - clamp((g[1] - ybar) / g[2], -1.0, 1.0)
    if abs(g[1] - ybar) < 1e-10
      break
    end
  end
  (obj, old_obj, c, v) = (0.0, 0.0, 0.0, 0.0)
  epsilon = 1e-8
  for iteration = 1:100 # scoring algorithm
    fill!(score, 0.0)
    fill!(inform, 0.0)
    mul!(z, X, beta) # z = X * beta
    for i = 1:n
      f = meanf(z[i])
      v = varf(f[1])
      c = ((y[i] - f[1]) / v) * f[2]
      copyto!(x, X[i, :])
      BLAS.axpy!(c, x, score) # score = score + c * x
      c = f[2]^2 / v
      BLAS.ger!(c, x, x, inform) # inform = inform + c * x * x'
    end
    increment = inform \ score
    beta = beta + increment
    steps = -1
    fill!(score, 0.0)
    for step_halve = 0:3 # step halving
      obj = 0.0
      mul!(z, X, beta) # z = X * beta
      steps = steps + 1
      for i = 1:n
        f = meanf(z[i])
        v = varf(f[1])
        c = ((y[i] - f[1]) / v) * f[2]
        copyto!(x, X[i, :])
        BLAS.axpy!(c, x, score) # score = score + c * x
        obj = obj + loglikelihood(y[i], f[1])
      end
      if obj > old_obj
        break
      else
        beta = beta - increment
        increment = 0.5 * increment
      end
    end
    println(iteration," ",old_obj," ",obj," ",steps)
    if iteration > 1 && abs(obj - old_obj) < epsilon * (abs(old_obj) + 1.0)
      return (beta, obj)
    else
      old_obj = obj
    end
  end
  return (beta, obj)
end # function glm_regress

"""
Computes the score statistic for a generalized linear regression.
X is the design matrix, y is the response vector,
beta is the vector of regression coefficients,
meanf is the value and derivative of the inverse link,
and varf is the variance function of the mean.
"""
function glm_score_statistic(X::Matrix, y::Vector, beta::Vector,
  meanf::Function, varf::Function)
#
  (n, p) = size(X)
  @assert n == length(y)
  @assert p == length(beta)
  (score, inform) = (zeros(p), zeros(p, p))
  (x, z) = (zeros(p), zeros(n))
  mul!(z, X, beta) # z = X * beta
  for i = 1:n
    f = meanf(z[i])
    v = varf(f[1])
    c = ((y[i] - f[1]) / v) * f[2]
    copyto!(x, X[i, :])
    BLAS.axpy!(c, x, score) # score = score + c * x
    c = f[2]^2 / v
    BLAS.ger!(c, x, x, inform) # inform = inform + c * x * x'
  end
  x = inform \ score
  score_statistic = dot(score, x)
  return score_statistic
end # function glm_score_statistic

# Testing Routines
#
# n = 10
# a = blanks(n)
# println(a,"  ",typeof(a))
# b = empties(n)
# println(b,"  ",typeof(b))
# c = repeated_string(["av", "bc", "qrb", "bc"])
# println(c)
# integerset = BitSet([7, 5, 4, 10])
# d = select_set_element(integerset, 3)
# println(integerset,"  ",d)
# frequency = rand(n)
# frequency = frequency / sum(frequency)
# @time for i = 1:1000
#   s = random_category(frequency)
#   if s < 1 || s > n
#     println("category error")
#   end
# end
# n = 100
# x = Vector{Union{Missing, Float64}}(undef, n);
# x[1:n - 1] = randn(n - 1)
# println(typeof(x), size(x))
# standarize!(x)
# println(sample_mean_std(x))
# s = sample_stats(x)
# print_sample_stats(s)
# pvalue = 0.1*rand(n)
# fdr = collect(0.1:0.1:0.5)
# tests = 1000
# (number_passing, threshold) = simes_fdr(pvalue, fdr, tests)
# function GaussMean(u)
#   return [u, one(eltype(u))]
# end
# function GaussVar(mu)
#   return one(eltype(mu))
# end
# function GaussLoglikelihood(y, mu)
#   - (y - mu)^2
# end
# X = [68., 49, 60, 68, 97, 82, 59, 50, 73, 39, 71, 95, 61, 72, 87,
#   40, 66, 58, 58, 77];
# y = [75., 63, 57, 88, 88, 79, 82, 73, 90, 62, 70, 96, 76, 75, 85,
#   40, 74, 70, 75, 72];
# X = [ones(size(X, 1)) X];
# # Jennrich section 1.4 least squares problem
# (beta, obj) = glm(X, y, GaussMean, GaussVar, GaussLoglikelihood)
# println(beta," ",obj)
# println(" ")
# function LogisticMean(u)
#   p = exp(u)
#   p = p / (1 + p)
#   return [p, p * (1 - p)]
# end
# function LogisticVar(mu)
#   return mu * (1 - mu)
# end
# function LogisticLoglikelihood(y, mu)
#   y * log(mu) + (1 - y) * log(1 - mu)
# end
# X = [0.50, 0.75, 1.00, 1.25, 1.50, 1.75, 1.75, 2.00, 2.25, 2.50, 2.75,
#   3.00, 3.25, 3.50, 4.00, 4.25, 4.50, 4.75, 5.00, 5.50];
# y = [0., 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1];
# X = [ones(size(X,1)) X];
# # Wikipedia logistic regression problem
# (beta, obj) = glm(X, y, LogisticMean, LogisticVar, LogisticLoglikelihood)
# println(beta," ",obj)
# function PoissonMean(u)
#   p = exp(u)
#   return [p, p]
# end
# function PoissonVar(mu)
#   return mu
# end
# function PoissonLoglikelihood(y, mu)
#   y * log(mu) - mu
# end
# println(" ")
# X = [1.0, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14];
# X = reshape(X, 14, 1);
# y = [0., 1 ,2 ,3, 1, 4, 9, 18, 23, 31, 20, 25, 37, 45];
# # AIDs Poisson regression problem
# (beta, obj) = glm(X, y, PoissonMean, PoissonVar, PoissonLoglikelihood)
# println(beta," ",obj)
# beta = [zero(1); beta];
# X = [ones(size(X,1)) X];
# score_statistic = glm_score_statistic(X, y, beta, PoissonMean, PoissonVar)
# println("score_statistic = ",score_statistic)

max_stars_repo_name: OpenMendel/MendelBase

max_issues_repo_name: OpenMendel/MendelBase

max_forks_repo_name: OpenMendel/MendelBase


text:
import Base: setindex!, getindex

"""
    HierarchicalImage{T<:Number} <: AbstractArray{T,2}

Image type which dynamically allocated memory for pixels when their value is set, the value of unset pixels is assumed to be zero.

This is used for the detector image of [`OpticalSystem`](@ref)s which can typically be very high resolution, but often have a large proportion of the image blank.
"""
struct HierarchicalImage{T<:Number} <: AbstractArray{T,2} #allow for complex values. May want to expand this even further to allow more complicated detector image types
    toplevel::Array{Array{T,2},2}
    secondlevelrows::Int64
    secondlevelcols::Int64
    apparentrows::Int64
    apparentcols::Int64

    function HierarchicalImage{T}(height, width) where {T<:Number}
        heightfirst, heightsecond = nearestsqrts(height)
        widthfirst, widthsecond = nearestsqrts(width)
        return HierarchicalImage{T}(height, width, heightfirst, widthfirst, heightsecond, widthsecond)
    end
    HierarchicalImage{T}(apparentheight, apparentwidth, height, width, secondlevelheight, secondlevelwidth) where {T<:Number} = new{T}(Array{Array{T,2},2}(undef, height, width), secondlevelheight, secondlevelwidth, apparentheight, apparentwidth)
end
export HierarchicalImage

function nearestsqrts(num)
    rt1 = Int64(ceil(sqrt(num)))
    rt2 = Int64(floor(sqrt(num)))
    if rt1 * rt2 >= num
        return (rt1, rt2)
    else
        return (rt1, rt2 + 1)
    end
end

toplevel(a::HierarchicalImage) = a.toplevel
secondlevelrows(a::HierarchicalImage) = a.secondlevelrows
secondlevelcols(a::HierarchicalImage) = a.secondlevelcols
apparentrows(a::HierarchicalImage) = a.apparentrows
apparentcols(a::HierarchicalImage) = a.apparentcols

topindex(a::HierarchicalImage, index1, index2)::Int = linearindex(size(toplevel(a), 1), 1 + (index1 - 1) ÷ secondlevelrows(a), 1 + (index2 - 1) ÷ secondlevelcols(a))
lowerindex(a::HierarchicalImage, index1, index2)::Int = linearindex(secondlevelrows(a), 1 + mod(index1 - 1, secondlevelrows(a)), 1 + mod(index2 - 1, secondlevelcols(a)))

linearindex(rows::Int, index1::Int, index2::Int)::Int = index1 + rows * (index2 - 1)

function outofbounds(a::HierarchicalImage, indices::Tuple{Int,Int})
    if any(indices .> size(a))
        throw(BoundsError(a, indices))
    end
end

Base.eltype(::HierarchicalImage{T}) where {T<:Number} = T
Base.length(a::HierarchicalImage{T}) where {T<:Number} = reduce(*, size(a))
Base.size(a::HierarchicalImage{T}) where {T<:Number} = (apparentrows(a), apparentcols(a))

function Base.getindex(a::HierarchicalImage{T}, indices::Vararg{Int,2}) where {T<:Number}
    outofbounds(a, indices)

    top = topindex(a, indices[1], indices[2])
    bott = lowerindex(a, indices[1], indices[2])

    return getindexdirect(a, top, bott)
end

function getindexdirect(a::HierarchicalImage{T}, topidx::Int, botidx::Int) where {T<:Real}
    if !isassigned(toplevel(a), topidx) #by default any part of the image that was not written to has the value zero
        return zero(T)
    else
        return toplevel(a)[topidx][botidx]
    end
end

function Base.setindex!(a::HierarchicalImage, v::T, indices::Vararg{Int,2}) where {T<:Number} #if use HierarchicalImage{T} get error "setindex! not defined". No idea why this happens
    outofbounds(a, indices)

    top = topindex(a, indices[1], indices[2])
    bott = lowerindex(a, indices[1], indices[2])

    setindexdirect!(a, v, top, bott)
end

function setindexdirect!(a::HierarchicalImage{T}, v::T, topidx::Int, botidx::Int) where {T<:Real}
    if !isassigned(toplevel(a), LinearIndices(toplevel(a))[topidx]) #by default any part of the image that has not been not written to is set to zero
        toplevel(a)[topidx] = fill(zero(T), secondlevelrows(a), secondlevelcols(a))
    end
    toplevel(a)[topidx][botidx] = v
end

"""
    reset!(a::HierarchicalImage{T})

Resets the pixels in the image to `zero(T)`.
Do this rather than `image .= zero(T)` because that will cause every pixel to be accessed, and therefore allocated.
For large images this can cause huge memory traffic.
"""
function reset!(a::HierarchicalImage{T}) where {T}
    for index in 1:length(toplevel(a))
        if isassigned(toplevel(a), index)
            temp = toplevel(a)[index]
            for k in 1:length(temp)
                temp[k] = zero(T)
            end
        end
    end
end

"""
    sum!(a::HierarchicalImage{T}, b::HierarchicalImage{T})

Add the contents of `b` to `a` in an efficient way.
"""
function sum!(a::HierarchicalImage{T}, b::HierarchicalImage{T}) where {T}
    @assert size(a) == size(b)
    for index in 1:length(toplevel(b))
        if isassigned(toplevel(b), index)
            temp = toplevel(b)[index]
            for k in 1:length(temp)
                setindexdirect!(a, getindexdirect(a, index, k) + temp[k], index, k)
            end
        end
    end
end

export sum!, reset!

ext: jl


text:
# Benchmark Optim.jl on logistic regression problems

using Optim, LineSearches
using LogisticOptTools
using Plots

const LOT = LogisticOptTools

# Choose dataset to load
DATASET = "covtype.libsvm.binary.bz2"
# Choose penalty with given coefficient λ
PENALTY = LOT.L2Penalty(1e-2)

if !isfile(DATASET)
    println("Download dataset from $LIBSVM_URL")
    download("$LIBSVM_URL/$DATASET", DATASET)
end

# Some plots utilities
"Display evolution of gradients along iteration."
function display_grad(results::Optim.OptimizationResults,
                      label::String,
                      p::Plots.Plot)
    exectime = [t.metadata["time"] for t in results.trace]
    exectime .-= exectime[1]
    grad = [t.g_norm for t in results.trace]
    plot!(p, exectime, log10.(grad), marker=true, label=label, markersize=1.)
end

"Display evolution of functions along iteration."
function display_func(results::Optim.OptimizationResults,
                      f♯,
                      label::String,
                      p::Plots.Plot)
    exectime = [t.metadata["time"] for t in results.trace]
    exectime .-= exectime[1]
    values = [t.value for t in results.trace]
    values .-= f♯
    plot!(p, exectime, log10.(values), marker=true, label=label, markersize=1.)
end

"Benchmark algorithms from Optim.jl"
function benchmark(algos, options, logreg)
    p = LOT.nfeatures(logreg)
    x = zeros(p)
    g = zeros(p)
    f, g!, _, _ = LOT.generate_callbacks(logreg)
    # Activate precompilation
    f(x)
    g!(g, x)

    results = []
    for algo in algos
        res = Optim.optimize(f, g!, x, algo, options)
        push!(results, res)
    end
    return results
end

bfgs1 = BFGS(alphaguess=InitialHagerZhang(α0=1.0),
             linesearch=LineSearches.HagerZhang())
bfgs2 = BFGS(alphaguess=InitialStatic(alpha=1.0),
             linesearch=LineSearches.MoreThuente())
cg = ConjugateGradient()

options = Optim.Options(iterations=100, store_trace=true, show_trace=true,
                        g_tol=1e-5, allow_f_increases=true)

# Load data
dataset = LOT.LogitData(DATASET, scale_data=true)
# Build logistic regression
logreg = LOT.LogisticRegressor(dataset.X, dataset.y, penalty=PENALTY, fit_intercept=true)
results = benchmark([bfgs1, bfgs2, cg], options, logreg)

# Plot evolution of gradient
pg = plot(markershape=:rect)
display_grad.(results, ["Optim BFGS HZ", "Optim BFGS MT", "Optim CG"], Ref(pg))
xlabel!(pg, "Time (s)")
ylabel!(pg, "Gradient (log scale)")

# Plot evolution of function
pf = plot(markershape=:rect)
# Get optimum
f♯ = minimum([res.minimum for res in results])
display_func.(results, f♯, ["Optim BFGS HZ", "Optim BFGS MT", "Optim CG"], Ref(pf))
xlabel!(pf, "Time (s)")
ylabel!(pf, "log(f - fopt) ")


text:
# using Documenter, SOM
#
# makedocs(modules=[SOM],
#          doctest=true)
#
# deploydocs(deps   = Deps.pip("mkdocs", "python-markdown-math"),
#            repo   = "github.com/andreasdominik/SOM.jl.git",
#            julia  = "0.6",
#            osname = "linux")
#
using Documenter, SOM

makedocs(modules=[SOM],
         clean = false,
         format = :html,
         assets = ["assets/favicon.ico"],
         sitename = "SOM.jl",
         authors = "Andreas Dominik",
         pages = [
                  "Introduction" => "index.md",
                  "Tutorials" => Any[
                            "Quick Start" => "tutorials/firstTutorial.md"
                                     ],
                   "API" => Any[
                            "Types" => "api/types.md",
                            "Training" => "api/soms.md",
                            "Visualisation" => "api/visualisations.md"
                               ],
                    "License" => "LICENSE.md"
                  ],
                  # Use clean URLs, unless built as a "local" build
          html_prettyurls = !("local" in ARGS),
          html_canonical = "https://LiScI-Lab.github.io/SOM.jl/stable/",
         )

deploydocs(repo   = "github.com/LiScI-Lab/SOM.jl.git",
           julia  = "1.0",
           osname = "linux",
           target = "build",
           deps = nothing,
           make = nothing)


text: #Written By Johnathan Bizzano


text:
const _JULIA_ORDERINGS =
    (:unordered, :monotonic, :acquire, :release, :acquire_release, :sequentially_consistent)

@inline function llvm_ordering_from_juila(julia_ordering::Symbol)
    julia_ordering in _JULIA_ORDERINGS || error("unknown ordering: ", julia_ordering)
    return getfield(UnsafeAtomics, julia_ordering)
end


text:
module actuator_examples
using FinEtools
using FinEtoolsHeatDiff
using FinEtools.MeshExportModule
using UnicodePlots

function actuator2()
    # MEMS actuator.   Thermal analysis.
    x0 =  0.0*phun("micro*m");
    x1 = x0+5.0*phun("micro*m");
    x2 = x1+10.0*phun("micro*m");
    x3 = x2+10.0*phun("micro*m");
    x4 = x3+10.0*phun("micro*m");
    y0 = 0.0*phun("micro*m");
    y4 = 250.0*phun("micro*m");
    y3 = y4-10.0*phun("micro*m");
    y2 = y3-10.0*phun("micro*m");
    y1 = y2-10.0*phun("micro*m");
    t = 2.0*phun("micro*m");
    h = 0.1*phun("micro*m");
    z0 = 0.0*phun("micro*m");
    z3 = 2*t+h;
    z2 = z3-t;
    z1 = z2-h;
    m1 = 2*2;
    m2 = 2*2;
    m3 = 2*2;
    m4 = 3*2;
    n1 = 20*2;
    n2 = 4*2;
    n3 = 2*2;
    n4 = 2*2;
    n5 = 7*2;
    p1 = 1*2;
    p2 = 1*2;
    p3 = 1*2;
    kappa = 157*[1 0 0; 0 1 0; 0 0 1]*phun("W/m/K"); # W/m/K, conductivity matrix
    DV = 5*phun("V"); # voltage drop in volt
    l  = 2*(y1+y2)/2+2*(x1+x2)/2; # length of the conductor
    resistivity  =  1.1e-5*phun("Ohm*m"); # Ohm m
    Q = DV^2/resistivity/l^2; # rate of Joule heating, W/m^3
    T_substrate = 293; # substrate temperature in degrees Kelvin
    tol = 1e6*eps(h)

    fens,fes =  H8hexahedron([x1 y0 z0; x2 y1 z1],m2,n1,p1);
    fens1,fes1  =  H8hexahedron([x1 y1 z0;x2 y2 z1],m2,n2,p1);
    fens,fes1,fes2  =  mergemeshes(fens1, fes1, fens, fes, tol);
    fes =  cat(fes1,fes2);
    fens1,fes1  =  H8hexahedron([x0 y1 z0;x1 y2 z1],m1,n2,p1);
    fens,fes1,fes2  =  mergemeshes(fens1, fes1, fens, fes, tol);
    fes =  cat(fes1,fes2);
    fens1,fes1  =  H8hexahedron([x0 y1 z1;x1 y2 z2], m1,n2,p2);
    fens,fes1,fes2  =  mergemeshes(fens1, fes1, fens, fes, tol);
    fes =  cat(fes1,fes2);
    fens1,fes1  =  H8hexahedron([x0 y1 z2;x1 y2 z3],m1,n2,p3);
    fens,fes1,fes2  =  mergemeshes(fens1, fes1, fens, fes, tol);
    fes =  cat(fes1,fes2);
    fens1,fes1  =  H8hexahedron([x0 y2 z2;x1 y3 z3],m1,n3,p3);
    fens,fes1,fes2  =  mergemeshes(fens1, fes1, fens, fes, tol);
    fes =  cat(fes1,fes2);
    fens1,fes1  =  H8hexahedron([x0 y3 z2;x1 y4 z3], m1,n4,p3);
    fens,fes1,fes2  =  mergemeshes(fens1, fes1, fens, fes, tol);
    fes =  cat(fes1,fes2);
    fens1,fes1  =  H8hexahedron([x1 y3 z2;x3 y4 z3],m4,n4,p3);
    fens,fes1,fes2  =  mergemeshes(fens1, fes1, fens, fes, tol);
    fes =  cat(fes1,fes2);
    fens1,fes1  =  H8hexahedron([x3 y3 z2;x4 y4 z3],m3,n4,p3);
    fens,fes1,fes2  =  mergemeshes(fens1, fes1, fens, fes, tol);
    fes =  cat(fes1,fes2);
    fens1,fes1  =  H8hexahedron([x3 y0 z2;x4 y3 z3], m3,n5,p3);
    fens,fes1,fes2  =  mergemeshes(fens1, fes1, fens, fes, tol);
    fes =  cat(fes1,fes2);

    fens,fes  =  H8toH20(fens,fes);

    hotmater = MatHeatDiff(kappa)
    coldmater = MatHeatDiff(kappa)
    cl =  selectelem(fens, fes, box=[x0,x2,y0,y2,z0,z1],inflate = t/100);

    hotfemm  =  FEMMHeatDiff(IntegDomain(subset(fes,cl), GaussRule(3, 3)), hotmater)
    coldfemm  = FEMMHeatDiff(IntegDomain(subset(fes,setdiff(collect(1:count(fes)), cl)),  GaussRule(3, 3)), coldmater)
    geom = NodalField(fens.xyz)
    Temp = NodalField(zeros(size(fens.xyz,1),1))
    fenids = selectnode(fens, box=[x0,x4,y0,y0,z0,z3], inflate=t/1000) ; # fixed temperature on substrate
    setebc!(Temp, fenids, true, 1, T_substrate)
    applyebc!(Temp)
    numberdofs!(Temp)

    K = conductivity(hotfemm, geom, Temp) + conductivity(coldfemm, geom, Temp)
    fi = ForceIntensity(FFlt[Q]);
    F = distribloads(hotfemm, geom, Temp, fi, 3) + nzebcloadsconductivity(hotfemm, geom, Temp) + nzebcloadsconductivity(coldfemm, geom, Temp)

    U = K\F
    scattersysvec!(Temp,U[:])

    nList = selectnode(fens, box=[x1,x1,y0,y1,z1,z1], inflate=t/1000)
    y_i = geom.values[nList, 2]
    T_i = Temp.values[nList, 1]
    ixi = sortperm(y_i)
    nList = selectnode(fens, box=[x3,x3,y0,y3,z2,z2], inflate=t/1000)
    y_o = geom.values[nList, 2]
    T_o = Temp.values[nList, 1]
    ixo = sortperm(y_o)
    # @pgf a = Axis({
    #     xlabel = "Distance",
    #     ylabel = "Temperature",
    #     title = "Transient pressure",
    #     legend_pos  = "south east"
    # },
    # Plot(Table([:x => y_o[ixo], :y => T_o[ixo]])), LegendEntry("cold leg"),
    # Plot(Table([:x => y_i[ixi], :y => T_i[ixi]])), LegendEntry("hot leg"))
    a = lineplot(y_o[ixo], T_o[ixo], name = "cold leg", ylim = [0, 1400])
    a = lineplot!(a, y_i[ixi], T_i[ixi], color = :red, name = "hot leg")
    display(a)

    File =  "a.vtk"
    MeshExportModule.vtkexportmesh(File, fens, fes; scalars =[("Temperature", Temp.values)])
    @async run(`"paraview.exe" $File`)
    true
end # actuator2

function allrun()
    println("#####################################################")
    println("# actuator2 ")
    actuator2()
    return true
end # function allrun

end # module actuator_examples

max_stars_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff

max_issues_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff

max_forks_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff


text:
"""
ДАНО: Где-то на неограниченном со всех сторон поле и без внутренних перегородок имеется единственный маркер. Робот - в произвольной клетке поля.
РЕЗУЛЬТАТ: Робот - в клетке с тем маркером.
"""

function detective!(r)
    steps=1
    side=Nord
    while ismarker(r)==false
        for a in 1:2
            detective!(r,side,steps)
            side=next_dir(side)
        end
        steps+=1
    end
end

function detective!(r,side,steps)
    for a in 1:steps
        if ismarker(r) == true
            return nothing
        end
        move!(r,side)
    end
end

next_dir(side::HorizonSide)=HorizonSide(mod(Int(side)+1,4))

max_stars_repo_name: alexgruzberg/mirea-progs

max_issues_repo_name: alexgruzberg/mirea-progs

max_forks_repo_name: alexgruzberg/mirea-progs


text:
# There are two files in .gz with .data extension, need to pick the right one

dataset_path = joinpath(dirname(@__FILE__), "cpu-small.data.orig")
file_path = joinpath(dirname(@__FILE__), "cpu-small.data.gz-contents",
                     "cpu_small.data")

write("$dataset_path.custom", readstring(file_path))

max_stars_repo_name: CoventryResearch/uci-data

max_issues_repo_name: CoventryResearch/uci-data

max_forks_repo_name: CoventryResearch/uci-data


text:
using LinearAlgebra
using RecipesBase

# polarization ellipse
@userplot PolEllipse

"""
    polellipse(S)
    polellipse!(S)

Plot the polarizatio ellipse for the given Stokes parameters. The polarization ellipse draws the path of the eletric field over one wavelength. The input vector should be a Stokes vector `[I, Q, U, V]`. This is a user recipe and can be composed with other plot recipes.
"""
polellipse

@recipe function f(pe::PolEllipse, N=1000)
    aspect_ratio --> 1
    xguide --> "x"
    yguide --> "y"
    alpha -->  range(0.2, 1, length=N)

    # extract elliptical params (a, b, θ) from Stokes vector
    stokes = first(pe.args)
    L = hypot(stokes[begin + 1], stokes[begin + 2])
    Ip = hypot(L, stokes[begin + 3])
    θ = 0.5 * atan(stokes[begin + 2], stokes[begin + 1])
    a = sqrt(0.5 * (Ip + L))
    b = sqrt(0.5 * (Ip - L))
    # draw ellipse
    t = range(0, 2π, length=N)
    x = @. a * cos(t) * cos(θ) - b * sin(t) * sin(θ)
    y = @. a * cos(t) * sin(θ) + b * sin(t) * cos(θ)
    return x, y
end


text:
export iauSxpv
"""
Multiply a pv-vector by a scalar.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Given:
   s       double          scalar
   pv      double[2][3]    pv-vector

Returned:
   spv     double[2][3]    s * pv

Note:
   It is permissible for pv and spv to be the same array

Called:
   iauS2xpv     multiply pv-vector by two scalars

This revision:  2013 August 7

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauSxpv(double s, double pv[2][3], double spv[2][3])

function iauSxpv(s::Real, pv::AbstractMatrix{<:Real})
   spv = zeros(Float64, 3, 2)

   ccall((:iauSxpv, libsofa_c), Cvoid,
       (Cdouble, Ptr{Cdouble}, Ptr{Cdouble}),
       convert(Float64, s),
       convert(Array{Float64, 2}, pv'),
       spv)

   return SMatrix{2,3}(spv')
end


text:
# start these lines in the REPL:
using Sockets
errormonitor(@async begin                            # 1
  server = listen(8080)
  while true
    conn = accept(server)
    @async while isopen(conn)
      write(conn, readline(conn, keep=true))         # 2
    end
  end
end)
# Task (runnable) @0x0000000008d40da0

client = connect(8080)                               # 3
# TCPSocket(Base.Libc.WindowsRawSocket(0x000000000000039c) open, 0 bytes waiting)

errormonitor(@async while isopen(client)             # 4
  write(stdout, readline(client, keep=true))
end)
# Task (runnable) @0x0000000008d41370

println(client, "Hello World from the Echo Server")  # 5
# Hello World from the Echo Server

close(client)                                         # 6


text:
module AperturePhotometry

include("calibration.jl")
include("photometry.jl")
include("stacking.jl")
include("load_images.jl")

end # module


text:
export Reshape

"""
`Reshape(A::AbstractOperator, dim_out...)`

Shorthand constructor:

`reshape(A, idx...)`

Reshape the codomain dimensions of an `AbstractOperator`.

```julia
julia> A = Reshape(DFT(10),2,5)
¶ℱ  ℝ^10 -> ℂ^(2, 5)

julia> R = reshape(Conv((19,),randn(10)),7,2,2)
¶★  ℝ^19 -> ℝ^(7, 2, 2)

```

"""
struct Reshape{N,L<:AbstractOperator} <: AbstractOperator
	A::L
	dim_out::NTuple{N,Int}

	function Reshape(A::L,dim_out::NTuple{N,Int}) where {N, L<:AbstractOperator}
		if prod(size(A,1)) != prod(dim_out)
			throw(DimensionMismatch("new dimensions $(dim_out) must be consistent with AbstractOperator codomain size $(size(A,1))"))
		end
		new{N,L}(A,dim_out)
	end
end

# Constructors

Reshape(A::L, dim_out::Vararg{Int,N}) where {N,L<:AbstractOperator} =
Reshape(A, dim_out)

# Mappings

function mul!(y::C, R::Reshape{N,L}, b::D) where {N,L,C,D}
	y_res = reshape(y,size(R.A,1))
	mul!(y_res, R.A, b)
end

function mul!(y::D, A::AdjointOperator{Reshape{N,L}}, b::C) where {N,L,C,D}
    R = A.A
	b_res = reshape(b,size(R.A,1))
	mul!(y, R.A', b_res)
end

# Properties

size(R::Reshape) = (R.dim_out, size(R.A,2))

  domainType(  R::Reshape) =   domainType(R.A)
codomainType(  R::Reshape) = codomainType(R.A)

is_linear(      R::Reshape) = is_linear(R.A)
is_null(        R::Reshape) = is_null(R.A)
is_eye(        R::Reshape)  = is_eye(R.A)
is_diagonal(    R::Reshape) = is_diagonal(R.A)
is_AcA_diagonal(R::Reshape) = is_AcA_diagonal(R.A)
is_AAc_diagonal(R::Reshape) = is_AAc_diagonal(R.A)
is_orthogonal(  R::Reshape) = is_orthogonal(  R.A)
is_invertible(  R::Reshape) = is_invertible(R.A)
is_full_row_rank(  R::Reshape)    = is_full_row_rank(     R.A)
is_full_column_rank(  R::Reshape) = is_full_column_rank(  R.A)

fun_name(R::Reshape) = "¶"*fun_name(R.A)
remove_displacement(R::Reshape) = Reshape(remove_displacement(R.A), R.dim_out)

function permute(T::Reshape{N,L}, p::AbstractVector{Int}) where {N,L}
    A = AbstractOperators.permute(T.A,p)
    return Reshape(A,T.dim_out)
end


text:
function ConfigurableLowThrustManeuver(arg0::ThrustDirectionAndAttitudeProvider, arg1::AbstractDetector, arg2::AbstractDetector, arg3::jdouble, arg4::jdouble)
    return ConfigurableLowThrustManeuver((ThrustDirectionAndAttitudeProvider, AbstractDetector, AbstractDetector, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4)
end

function acceleration(obj::Maneuver, arg0::FieldSpacecraftState, arg1::Vector{RealFieldElement})
    return jcall(obj, "acceleration", FieldVector3D, (FieldSpacecraftState, Vector{RealFieldElement}), arg0, arg1)
end

function acceleration(obj::Maneuver, arg0::SpacecraftState, arg1::Vector{jdouble})
    return jcall(obj, "acceleration", Vector3D, (SpacecraftState, Vector{jdouble}), arg0, arg1)
end

function add_contribution(obj::Maneuver, arg0::FieldSpacecraftState, arg1::FieldTimeDerivativesEquations)
    return jcall(obj, "addContribution", void, (FieldSpacecraftState, FieldTimeDerivativesEquations), arg0, arg1)
end

function add_contribution(obj::Maneuver, arg0::SpacecraftState, arg1::TimeDerivativesEquations)
    return jcall(obj, "addContribution", void, (SpacecraftState, TimeDerivativesEquations), arg0, arg1)
end

function depends_on_position_only(obj::Maneuver)
    return jcall(obj, "dependsOnPositionOnly", jboolean, ())
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_attitude_override(obj::Maneuver)
    return jcall(obj, "getAttitudeOverride", AttitudeProvider, ())
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_events_detectors(obj::Maneuver)
    return jcall(obj, "getEventsDetectors", Stream, ())
end

function get_field_events_detectors(obj::Maneuver, arg0::Field)
    return jcall(obj, "getFieldEventsDetectors", Stream, (Field,), arg0)
end

function get_isp(obj::ConfigurableLowThrustManeuver)
    return jcall(obj, "getISP", jdouble, ())
end

function get_maneuver_triggers(obj::Maneuver)
    return jcall(obj, "getManeuverTriggers", ManeuverTriggers, ())
end

function get_name(obj::Maneuver)
    return jcall(obj, "getName", JString, ())
end

function get_parameter_driver(obj::AbstractForceModel, arg0::JString)
    return jcall(obj, "getParameterDriver", ParameterDriver, (JString,), arg0)
end

function get_parameters(obj::ForceModel)
    return jcall(obj, "getParameters", Vector{jdouble}, ())
end

function get_parameters(obj::ForceModel, arg0::Field)
    return jcall(obj, "getParameters", Vector{RealFieldElement}, (Field,), arg0)
end

function get_parameters_drivers(obj::Maneuver)
    return jcall(obj, "getParametersDrivers", Vector{ParameterDriver}, ())
end

function get_propulsion_model(obj::Maneuver)
    return jcall(obj, "getPropulsionModel", PropulsionModel, ())
end

function get_thrust(obj::ConfigurableLowThrustManeuver)
    return jcall(obj, "getThrust", jdouble, ())
end

function get_thrust_direction_provider(obj::ConfigurableLowThrustManeuver)
    return jcall(obj, "getThrustDirectionProvider", ThrustDirectionAndAttitudeProvider, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function init(obj::Maneuver, arg0::SpacecraftState, arg1::AbsoluteDate)
    return jcall(obj, "init", void, (SpacecraftState, AbsoluteDate), arg0, arg1)
end

function is_supported(obj::AbstractForceModel, arg0::JString)
    return jcall(obj, "isSupported", jboolean, (JString,), arg0)
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
export BrickletMultiTouchIdentity
struct BrickletMultiTouchIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletMultiTouch
"""
Capacitive touch sensor for 12 electrodes
"""
mutable struct BrickletMultiTouch <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletMultiTouch(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_multi_touch")
        deviceInternal = package.BrickletMultiTouch(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_touch_state
"""
    $(SIGNATURES)

Returns the current touch state. The state is given as a bitfield.

Bits 0 to 11 represent the 12 electrodes and bit 12 represents
the proximity.

If an electrode is touched, the corresponding bit is *true*. If
a hand or similar is in proximity to the electrodes, bit 12 is
*true*.

Example: The state 4103 = 0x1007 = 0b1000000000111 means that
electrodes 0, 1 and 2 are touched and that something is in the
proximity of the electrodes.

The proximity is activated with a distance of 1-2cm. An electrode
is already counted as touched if a finger is nearly touching the
electrode. This means that you can put a piece of paper or foil
or similar on top of a electrode to build a touch panel with
a professional look.
"""
function get_touch_state(device::BrickletMultiTouch)
    return device.deviceInternal.get_touch_state()
end

export recalibrate
"""
    $(SIGNATURES)

Recalibrates the electrodes. Call this function whenever you changed
or moved you electrodes.
"""
function recalibrate(device::BrickletMultiTouch)
    device.deviceInternal.recalibrate()
end

export set_electrode_config
"""
    $(SIGNATURES)

Enables/disables electrodes with a bitfield (see :func:`Get Touch State`).

*True* enables the electrode, *false* disables the electrode. A
disabled electrode will always return *false* as its state. If you
don't need all electrodes you can disable the electrodes that are
not needed.

It is recommended that you disable the proximity bit (bit 12) if
the proximity feature is not needed. This will reduce the amount of
traffic that is produced by the :cb:`Touch State` callback.

Disabling electrodes will also reduce power consumption.

Default: 8191 = 0x1FFF = 0b1111111111111 (all electrodes and proximity feature enabled)
"""
function set_electrode_config(device::BrickletMultiTouch, enabled_electrodes)
    device.deviceInternal.set_electrode_config(enabled_electrodes)
end

export get_electrode_config
"""
    $(SIGNATURES)

Returns the electrode configuration, as set by :func:`Set Electrode Config`.
"""
function get_electrode_config(device::BrickletMultiTouch)
    return device.deviceInternal.get_electrode_config()
end

export set_electrode_sensitivity
"""
    $(SIGNATURES)

Sets the sensitivity of the electrodes. An electrode with a high sensitivity
will register a touch earlier then an electrode with a low sensitivity.

If you build a big electrode you might need to decrease the sensitivity, since
the area that can be charged will get bigger. If you want to be able to
activate an electrode from further away you need to increase the sensitivity.

After a new sensitivity is set, you likely want to call :func:`Recalibrate`
to calibrate the electrodes with the newly defined sensitivity.
"""
function set_electrode_sensitivity(device::BrickletMultiTouch, sensitivity)
    device.deviceInternal.set_electrode_sensitivity(sensitivity)
end

export get_electrode_sensitivity
"""
    $(SIGNATURES)

Returns the current sensitivity, as set by :func:`Set Electrode Sensitivity`.
"""
function get_electrode_sensitivity(device::BrickletMultiTouch)
    return device.deviceInternal.get_electrode_sensitivity()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletMultiTouch)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletMultiTouch, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
"""
    get_notes(csv::String) -> Vector{String}

Extract the values of the `Note` column from the CSV file.

# Arguments
- `csv::String`: $(CSV_ARGUMENT)

# Returns
- `Vector{String}`: the notes

# Throws
- Exceptions from: [`get_all`](@ref)

# Example
```jldoctest; output = false
using PDFHighlights
HEADER = PDFHighlights.Internal.Constants.HEADER

_file, io = mktemp()
println(io, HEADER, '\\n', ",,,Journal,")
flush(io)
file = _file * ".csv"
mv(_file, file)

get_notes(file) == ["Journal"]

# output

true
```
"""
function get_notes(csv::String)::Vector{String}
    _, _, _, notes, _ = get_all(csv)
    return notes
end


text:
using CSV
using DataFrames
using Plots
using JLD2
using DiscretisedFluidQueues
include("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/plot_utils.jl")
include("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/error_metrics.jl")

@halfwidth_plot_defaults()

plotlyjs()
function p()
ks_data = CSV.read((@__DIR__)*"/../data2/meshs_ks_func_count_4.csv",DataFrame)

plot()
linestyles_vec = [:solid,:dash,:dashdot, :dot]
markerstyles_vec = [:cross,:diamond,:circle,:dot]
for (c,col) in enumerate(names(ks_data))
    plot!(log10.(1:2:21),ks_data[:,col],
        label=col,linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],
        linewidth=2)
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="KS error - CDF")
plot!(legend=(0.225,0.55))
plot!()
@add_lines!(ks_data,("Unif","QBDRAP"),@__DIR__)
savefig((@__DIR__)*"/meshs_ks_error_formatted.pdf")
end
p()

function p()
l2_pdf_data = CSV.read((@__DIR__)*"/../data2/meshs_l2_pdf_func_count_4.csv",DataFrame)

plot()
linestyles_vec = [:solid,:dash,:dashdot, :dot]
markerstyles_vec = [:cross,:diamond,:circle,:cross]
for (c,col) in enumerate(names(l2_pdf_data))
    plot!(log10.(1:2:21),l2_pdf_data[:,col],
        label=col,linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],
        linewidth=2)
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="L² error - PDF")
plot!(legend=(0.225,0.55))
plot!()
@add_lines!(l2_pdf_data,("Unif","QBDRAP"),@__DIR__)
savefig((@__DIR__)*"/meshs_l2_pdf_error_formatted.pdf")
end
p()

function p()
@fullwidth_plot_defaults()
dir = @__DIR__()
os = 1:2:7
plot(layout=(4,length(os)))
linestyles_vec = [:solid,:solid,:solid,:solid]
markerstyles_vec = [:cross,:diamond,:circle]
xs = range(2+eps(),7-eps(),length=200)
for (c,o) in enumerate(os)
    model = BoundedFluidQueue(fill(0.0,2,2),[1.0; 0.0],10.0)

    mesh = DGMesh(0.0:10.0,o)
    dq = DiscretisedFluidQueue(model,mesh)
    jldopen(dir*"/../../../coeffs_matrix.jld2") do f
        coeffs = f["coeff_matrix"][4,c,1]
        dt = SFMDistribution(coeffs,dq)
        plot!(xs,pdf(dt).(xs,1),subplot=c,label=(c==1 ? "DG" : false),
            linestyle=linestyles_vec[1],ylims = (-0.1,1.1),  yticks=false,
            # marker=markerstyles_vec[1],
            linewidth=2,xticks=false)
    end
    plot!(subplot=c, yticks=(c==1 ? (0.1:0.4:1.2) : false),
        ylabel=(c==1 ? "DG" : false), grid=false,
        title = string("Dim. ",o))

    mesh = DGMesh(0.0:(1.0./o):10.0,1)
    dq = DiscretisedFluidQueue(model,mesh)
    jldopen(dir*"/../../../coeffs_matrix.jld2") do f
        coeffs = f["coeff_matrix"][4,c,2]
        dt = SFMDistribution(coeffs,dq)
        plot!(xs,pdf(dt).(xs,1),subplot=c+4,label=(c==1 ? "Unif." : false),
            linestyle=linestyles_vec[2],ylims = (-0.1,1.1),  yticks=false,
            # marker=markerstyles_vec[2],
            linewidth=2,xticks=false)
    end
    (c==1)&&plot!(subplot=c+4, yticks=(c==1 ? (0.1:0.4:1.2) : false),
        ylabel=(c==1 ? "Unif." : false), grid=false)

    mesh = DGMesh(0.0:10.0,o)
    dq = DiscretisedFluidQueue(model,mesh)
    jldopen(dir*"/../../../coeffs_matrix.jld2") do f
        coeffs = f["coeff_matrix"][4,c,3]
        dt = SFMDistribution(coeffs,dq)
        plot!(xs,pdf(dt).(xs,1),subplot=c+8,label=(c==1 ? "DG (limiter)" : false),
            linestyle=linestyles_vec[3],ylims = (-0.1,1.1),  yticks=false,
            # marker=markerstyles_vec[3],
            linewidth=2,xticks=false)
    end
    (c==1)&&plot!(subplot=c+8, yticks=(c==1 ? (0.1:0.4:1.2) : false),
        ylims = (-0.1,1.1),
        ylabel=(c==1 ? "DG-limit" : false), grid=false)

    mesh = FRAPMesh(0.0:10.0,o)
    dq = DiscretisedFluidQueue(model,mesh)
    jldopen(dir*"/../../../coeffs_matrix.jld2") do f
        coeffs = f["coeff_matrix"][4,c,4]
        dt = SFMDistribution(coeffs,dq)
        plot!(xs,pdf(dt).(xs,1),subplot=c+12,label=(c==1 ? "QBD-RAP" : false),
            linestyle=linestyles_vec[3],ylims = (-0.1,1.1), yticks=false,
            # marker=markerstyles_vec[3],
            linewidth=2,xticks=0:2:10, grid=false)
    end
    (c==1)&&plot!(subplot=c+12, yticks=(c==1 ? (0.1:0.4:1.2) : false),
        ylims = (-0.1,1.1),
        ylabel=(c==1 ? "QBD-RAP" : false), grid=false)
    # plot!(xs,x->2*x*(x<0.5) + Float64(x>=0.5),subplot=c,label=(c==1 ? "truth" : false),
    #     linestyle=:dot,
    #     # marker=markerstyles_vec[3],
    #     color=:grey,
    #     linewidth=2)

    # plot!(subplot=c, yticks=(c==1 ? (0.1:0.3:1.2) : false),
    #     ylims = (-0.1,1.1), xticks = 0.0:3:10,
    #     ylabel=(c==1 ? "probability" : false), grid=false,
    #     title = string("Dim. ",o),
    #     legend=:outertopright)
    plot!(legend=false)
end
plot!()
savefig((@__DIR__)*"/pdfs_formatted.pdf")
end
p()

function p()
L1_cell_probs_errors = DataFrame(DG = Float64[], Unif = Float64[], DG_limiter = Float64[], QBDRAP = Float64[])
truth = zeros(10)
truth[5] = 1.0
truth_pm = zeros(3)
for (c_o,o) in enumerate(1:2:21)
    row = zeros(4)
    for (c_m,m) in enumerate(names(L1_cell_probs_errors))
        jldopen((@__DIR__)*"/../../../coeffs_matrix.jld2") do f
            model = BoundedFluidQueue(fill(0.0,2,2),[1.0; 0.0],10.0)

            mesh = (c_m<4) ? DGMesh(0.0:10.0,o) : FRAPMesh(0.0:10.0,o)
            dq = DiscretisedFluidQueue(model,mesh)

            coeffs = f["coeff_matrix"][4,c_o,c_m]
            dt = SFMDistribution(coeffs,dq)

            x_vals = range(0.5,9.5;length=10)

            approx = sum(cell_probs(dt).(x_vals,(1:2)'),dims=2)
            approx_pm = [coeffs[1:N₋(model)];coeffs[end-N₊(model)+1:end]]

            row[c_m] = log10(sum(abs.(approx-truth)) + sum(abs.(approx_pm-truth_pm)))
        end
    end
    push!(L1_cell_probs_errors,row)
end
@halfwidth_plot_defaults()
file = (@__DIR__)*"/../data2/meshs_l1_cell_"*"func_count_4"
CSV.write(file*".csv",L1_cell_probs_errors)
q = plot()
linestyles_vec = [:solid,:dash,:dashdot, :dot]
markerstyles_vec = [:cross,:diamond,:circle,:dot]
for (c,reconstruction) in enumerate(names(L1_cell_probs_errors))
    plot!(log10.(1:2:21),L1_cell_probs_errors[:,reconstruction],label=reconstruction,
        linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],
        linewidth=2)
end
plot!(xlabel="Dimension", ylabel="Error",
    title="Error between cell probabilities",
    legend=(0.225,0.55))
error_ticks!(plot!())
@add_lines!(L1_cell_probs_errors,("Unif","QBDRAP"),@__DIR__)
savefig((@__DIR__)*"/L1_cell_probs.pdf")
end
p()

# os = 1:2:7
# plot(layout=(1,length(os)))
# linestyles_vec = [:solid,:dash,:dashdot]
# markerstyles_vec = [:cross,:diamond,:circle]
# xs = range(eps(),1-eps(),length=200)
# for (c,o) in enumerate(os)
#     model = BoundedFluidQueue(fill(0.0,1,1),[1.0],1.0)

#     mesh = DGMesh(0.0:1.0,o)
#     dq = DiscretisedFluidQueue(model,mesh)
#     d0 = SFMDistribution((x,i)->2*Float64(x<=0.5),dq,TrapezoidRule;fun_evals=10_001)
#     plot!(xs,pdf(d0).(xs,1),subplot=c,label=(c==1 ? "DG" : false),
#         linestyle=linestyles_vec[1],
#         # marker=markerstyles_vec[1],
#         linewidth=2)

#     mesh = DGMesh(0.0:(1.0./o):1.0,1)
#     dq = DiscretisedFluidQueue(model,mesh)
#     d0 = SFMDistribution((x,i)->2*Float64(x<=0.5),dq,TrapezoidRule;fun_evals=10_001)
#     plot!(xs,pdf(d0).(xs,1),subplot=c,label=(c==1 ? "Unif." : false),
#         linestyle=linestyles_vec[2],
#         # marker=markerstyles_vec[2],
#         linewidth=2)

#     mesh = FRAPMesh(0.0:1.0,o)
#     dq = DiscretisedFluidQueue(model,mesh)
#     d0 = SFMDistribution((x,i)->2*Float64(x<=0.5),dq,TrapezoidRule;fun_evals=10_001)
#     plot!(xs,pdf(d0).(xs,1),subplot=c,label=(c==1 ? "QBD-RAP" : false),
#         linestyle=linestyles_vec[3],
#         # marker=markerstyles_vec[3],
#         linewidth=2)

#     plot!(xs,x->2*Float64(x<=0.5),subplot=c,label=(c==1 ? "truth" : false),
#         linestyle=:dot,
#         # marker=markerstyles_vec[3],
#         color=:grey,
#         linewidth=2)

#     plot!(subplot=c, yticks=(c==1 ? (-2:2:6) : false),
#         ylims = (-3,7.2), xticks = 0.25:0.5:0.75,
#         ylabel=(c==1 ? "probability" : false), grid=false,
#         title = string("Dim. ",o),
#         legend=(0.1,0.8))
# end
# plot!()
# savefig((@__DIR__)*"/pdfs_formatted.pdf")


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    SingleBridgeOptimizer{BT<:AbstractBridge, OT<:MOI.ModelLike} <:
    AbstractBridgeOptimizer

The `SingleBridgeOptimizer` bridges any constrained variables supported by the
bridge `BT`. This is in contrast with the [`MathOptInterface.Bridges.LazyBridgeOptimizer`](@ref)
which only bridges the constrained variables that are unsupported by the
internal model, even if they are supported by one of its bridges.

!!! note
    Two bridge optimizers using variable bridges cannot be used together as both
    of them assume that the underlying model only returns variable indices with
    nonnegative values.
"""
mutable struct SingleBridgeOptimizer{BT<:AbstractBridge,OT<:MOI.ModelLike} <:
               MOIB.AbstractBridgeOptimizer
    model::OT
    map::Map # index of bridged variable -> variable bridge
    var_to_name::Dict{MOI.VariableIndex,String}
    name_to_var::Union{Dict{String,MOI.VariableIndex},Nothing}
    con_to_name::Dict{MOI.ConstraintIndex,String}
    name_to_con::Union{Dict{String,MOI.ConstraintIndex},Nothing}
end

function SingleBridgeOptimizer{BT}(model::OT) where {BT,OT<:MOI.ModelLike}
    return SingleBridgeOptimizer{BT,OT}(
        model,
        Map(),
        Dict{MOI.VariableIndex,String}(),
        nothing,
        Dict{MOI.ConstraintIndex,String}(),
        nothing,
    )
end

function bridges(::MOI.Bridges.AbstractBridgeOptimizer)
    return EmptyMap()
end

bridges(bridge::SingleBridgeOptimizer) = bridge.map

MOIB.supports_constraint_bridges(::SingleBridgeOptimizer) = false

function MOIB.supports_bridging_constrained_variable(
    ::SingleBridgeOptimizer{BT},
    S::Type{<:MOI.AbstractSet},
) where {BT}
    return supports_constrained_variable(BT, S)
end

function MOIB.is_variable_bridged(
    ::SingleBridgeOptimizer,
    ::Type{<:MOI.AbstractSet},
)
    return true
end

function MOIB.is_bridged(b::SingleBridgeOptimizer, S::Type{<:MOI.AbstractSet})
    return MOIB.supports_bridging_constrained_variable(b, S)
end

function MOIB.is_bridged(
    ::SingleBridgeOptimizer,
    ::Type{<:MOI.AbstractFunction},
    ::Type{<:MOI.AbstractSet},
)
    return false
end

function MOIB.is_bridged(
    ::SingleBridgeOptimizer,
    ::Type{<:MOI.AbstractScalarFunction},
)
    return false
end

function MOIB.bridge_type(
    ::SingleBridgeOptimizer{BT},
    ::Type{<:MOI.AbstractSet},
) where {BT}
    return BT
end

MOIB.bridging_cost(::SingleBridgeOptimizer, args...) = 1.0

MOIB.recursive_model(b::SingleBridgeOptimizer) = b.model

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


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

#Initialise variables
GeneIsoformRelationships = "data/gencode_human_gene_isoform_relationships.txt"
NumSimulations = 100
NumIsoformsToSimulate = 4
filteringThreshold = 4
alpha = nothing
beta = nothing
scnorm_counts = ["H1_24_scnorm_counts.txt",
"H1_96_scnorm_counts.txt",
"H9_24_scnorm_counts.txt",
"H9_96_scnorm_counts.txt"]
error = collect(0:0.1:0.5)

df = DataFrame(counts=String[], pFP = Float64[], pFN = Float64[],
NumIsos = Int64[], mean = Float64[], std = Float64[], median = Float64[])

for countsMatrix in scnorm_counts

    for i in 1:length(error) # i = pFN
        for j in 1:length(error) # j = pFP

            #output directories
            output_both = "figures/figure3/" * string(error[i]) * "_" * string(error[j]) * "_" * countsMatrix
            output_both_data = "figures/figure3_data/" * string(error[i]) * "_" * string(error[j]) * "_" * countsMatrix

            #simulate
            out = globalPredict("data/" * countsMatrix, GeneIsoformRelationships,
            NumSimulations, output_both, NumIsoformsToSimulate, filteringThreshold,
            false, alpha, beta, error[i], error[j], "none", "Weibull", output_both_data)

            newRow = DataFrame(counts = repeat([countsMatrix], 4),
            pFP = repeat([error[j]], 4), pFN = repeat([error[i]], 4),
            NumIsos = [1,2,3,4], mean = out[1], std = out[2], median = out[3])
            append!(df, newRow)

        end

    end
    # globalPredict("data/" * countsMatrix, GeneIsoformRelationships,
    # NumSimulations, "figures/figure3/random_" * countsMatrix , NumIsoformsToSimulate, filteringThreshold,
    # false, alpha, beta, 0.04, 0.01, "random", "Weibull", "none")
    # globalPredict("data/" * countsMatrix, GeneIsoformRelationships,
    # NumSimulations, "figures/figure3/IsoformDependence_" * countsMatrix , NumIsoformsToSimulate, filteringThreshold,
    # false, alpha, beta, 0.04, 0.01, "IsoformDependence", "Weibull", "none")
end

jldopen("figures/figure3_data/means_df.jld2", "w") do file
    file["df"] = df
end

p1 = Plots.plot(df.error, df.mean, markershapes =[:circle, :star5, :x, :square], group = df.NumIsos)
savefig(p1, "figures/figure3/mean_scnorm_err.png")
p1 = Plots.plot(df.error, df.median, markershapes =[:circle, :star5, :x, :square], group = df.NumIsos)
savefig(p1, "figures/figure3/median_scnorm_err.png")
p1 = Plots.plot(df.error, df.std, markershapes =[:circle, :star5, :x, :square], group = df.NumIsos)
savefig(p1, "figures/figure3/std_scnorm_err.png")

using DataFrames
using Plots, Plots.Measures
using JLD2
using DelimitedFiles

error = collect(0:0.1:0.5)

function make_hist(path)
    Arr = readdlm(path)

    #Untransformed plot
    p1 = histogram(Arr,
    xlims = (0, 4),
    ylims = (0, 5),
    normalize=:pdf,
    bins = 100,
    axis = (font(80)), titlefontsize=80, yticks = 0:5:10, xticks = 0:1:4,
    margin = 25mm)

    return p1
end

plot_arr = []

for i in 1:length(error) # i = pFN
    for j in 1:length(error) # j = pFP

        push!(plot_arr, make_hist("figures/figure3_data/" *
        string(error[i]) * "_" * string(error[j]) * "_" *
        "H1_24_scnorm_counts.txt_1_isoforms_array.txt"))
    end
end

tmp = plot(plot_arr..., layout =(6,6),
size = (720 * 8, 720 * 8), legend = false, margin = 25mm)

savefig(tmp, "figures/figure3/figure3A.png")

#make a plot for FP=0, FN=0
Arr = readdlm("figures/figure3_data/0.0_0.0_H1_24_scnorm_counts.txt_1_isoforms_array.txt")

#Untransformed plot
p_intro = histogram(Arr,
title = "1 Isoform",
xlims = (0, 4),
ylims = (0, 5),
ylabel = "Probability Density",
xlabel = "Mean No. Isoforms per Gene per Cell",
normalize=:pdf,
bins = 100,
axis = (font(16)), titlefontsize=16, yticks = 0:5,
legend = false)

savefig(p_intro, "figures/figure3/zero_errors_plot.png")

c = jldopen("figures/figure3_data/means_df.jld2", "r")
df = c["df"]

both = df[df.pFP .== df.pFN, :]
pFP = df[df.pFN .== 0.00, :]
pFN = df[df.pFP .== 0.00, :]

both_arr =[]
pFN_arr =[]
pFP_arr =[]

for i in 1:4
    tmp = both[both.NumIsos .== i, :]
    both_plot = Plots.plot(tmp.pFP, tmp.mean, group = tmp.counts, title = string(i))
    push!(both_arr, both_plot)

    tmp = pFN[pFN.NumIsos .== i, :]
    pFN_plot = Plots.plot(tmp.pFN, tmp.mean, group = tmp.counts, title = string(i))
    push!(pFN_arr, pFN_plot)

    tmp = pFP[pFP.NumIsos .== i, :]
    pFP_plot = Plots.plot(tmp.pFP, tmp.mean, group = tmp.counts, title = string(i))
    push!(pFP_arr, pFP_plot)
end

p1 = plot(both_arr[1], legend = :none, xlabel = "pFP and pFN", ylabel = "mean",
       size = (180 * 8, 180*8),titlefontsize = 80, title = "1 Isoform",
       linestyle=:solid, linealpha=0.5, linewidth=2*8, axis = font(80), margin = 15mm)
savefig(p1, "figures/figure3/both_means.png")
p1 = plot(pFN_arr[1], legend = :none, xlabel = "pFN", ylabel = "mean",
       size = (180 * 8, 180*8),titlefontsize = 80, title = "1 Isoform",
       linestyle=:solid, linealpha=0.5, linewidth=2*8, axis = font(80), margin = 10mm)
savefig(p1, "figures/figure3/pFN_means.png")
p1 = plot(pFP_arr[1], legend = :none, xlabel = "pFP", ylabel = "mean",
       size = (180 * 8, 180*8),titlefontsize = 80, title = "1 Isoform",
       linestyle=:solid, linealpha=0.5, linewidth=2*8, axis = font(80), margin = 10mm)
savefig(p1, "figures/figure3/pFP_means.png")

#get a legend
#get a legend
       p1 = plot(both_arr[1], legend = :inside,
              labels = ["H1 4 million reads", "H1 1 million reads", "H9 4 million reads",
              "H9 1 million reads"], size = (180 * 8, 180*8),titlefontsize = 50, legendfontsize=80, grid=false,
                  xlims=(-1,-0.5), showaxis=false,linestyle=:solid, linealpha=0.5, linewidth=2*8, xlabel = "", ylabel = "", title="")
savefig(p1, "figures/figure3/means_legend.png")

#########################
# overlap
#########################

function make_hist(path)
    Arr = readdlm(path)

    #Untransformed plot
    p1 = histogram(Arr,
    xlims = (0, 1),
    ylims = (0, 10),
    normalize=:pdf,
    bins = 50,
    color = palette(:default)[2],
    axis = (font(80)), titlefontsize=80, yticks = 0:10:10, xticks = 0:1:1,
    margin = 25mm)

    return p1
end

error = collect(0:0.1:0.5)
plot_arr = []

for i in 1:length(error) # i = pFN
    for j in 1:length(error) # j = pFP

        push!(plot_arr, make_hist("figures/figure3_data/" *
        string(error[i]) * "_" * string(error[j]) * "_" *
        "H1_24_scnorm_counts.txt_1_isoforms_overlap_array.txt"))
    end
end

tmp = plot(plot_arr..., layout =(6,6),
size = (720 * 8, 720 * 8), legend = false, margin = 25mm)

savefig(tmp, "figures/figure3/figure3A_overlap.png")

max_stars_repo_name: jenni-westoby/Obstacles

max_issues_repo_name: jenni-westoby/Obstacles

max_forks_repo_name: jenni-westoby/Obstacles


text:
pages = [
        "NeuralPDE.jl: Scientific Machine Learning (SciML) for Partial Differential Equations" => "index.md",
        "Physics-Informed Neural Network Tutorials" => Any[
            "pinn/poisson.md",
            "pinn/wave.md",
            "pinn/2D.md",
            "pinn/system.md",
            "pinn/3rd.md",
            "pinn/low_level.md",
            "pinn/ks.md",
            "pinn/fp.md",
            "pinn/parm_estim.md",
            "pinn/heterogeneous.md",
            "pinn/integro_diff.md",
            "pinn/debugging.md",
            "pinn/neural_adapter.md",
        ],
        "Specialized Neural PDE Tutorials" => Any[
            "examples/kolmogorovbackwards.md",
            "examples/optimal_stopping_american.md",
        ],
        "Specialized Neural ODE Tutorials" => Any[
            "examples/ode.md",
            "examples/nnrode_example.md",
        ],
        "API Documentation" => Any[
            "solvers/pinns.md",
            "solvers/kolmogorovbackwards_solver.md",
            "solvers/optimal_stopping.md",#TODO
            "solvers/ode.md",
            "solvers/nnrode.md",#TODO
        ]
    ]


text:
mutable struct LDAPConnection
    handle::Ptr{Cvoid}
    url::String
end

# typedef struct ldap_url_desc
struct LDAPURLDesc
    lud_next::Ptr{Cvoid}
    lud_scheme::Cstring
    lud_host::Cstring
    lud_port::Cint
    lud_dn::Cstring
    lud_attrs::Ptr{Ptr{UInt8}}
    lud_scope::Cint
    lud_filter::Cstring
    lud_exts::Ptr{Ptr{UInt8}}
    lud_crit_exts::Cint
end

struct URL
    scheme::String
    host::String
    port::Int
    dn::Union{Nothing, String}
    attrs::Vector{String}
    scope::Int
    filter::Union{Nothing, String}
    exts::Vector{String}
    crit_exts::Int
end

Base.:(==)(ua::URL, ub::URL) = ua.scheme == ub.scheme && ua.host == ub.host && ua.port == ub.port && ua.dn == ub.dn && ua.attrs == ub.attrs && ua.scope == ub.scope && ua.filter == ub.filter && ua.exts == ub.exts && ua.crit_exts == ub.crit_exts

abstract type AuthenticationResult end

struct AuthOk <: AuthenticationResult
    uri::String
    who::String
end

struct AuthErr <: AuthenticationResult
    uri::String
    who::String
    err_code::Cint
end

mutable struct MessageChain
    handle::Ptr{Cvoid}
    ldap::LDAPConnection

    function MessageChain(handle::Ptr{Cvoid}, ldap::LDAPConnection)
        chain = new(handle, ldap)
        finalizer(destroy!, chain)
        return chain
    end
end

abstract type AbstractMessage end

mutable struct Message{msgtype} <: AbstractMessage
    handle::Ptr{Cvoid}
    id::Int
    chain::MessageChain

    function Message(handle::Ptr{Cvoid}, chain::MessageChain)
        id = ldap_msgid(handle)
        msgtype = LDAPMessageType(ldap_msgtype(handle))
        return new{msgtype}(handle, id, chain)
    end
end

const Entry = Message{LDAP_RES_SEARCH_ENTRY}
const Reference = Message{LDAP_RES_SEARCH_REFERENCE}
const Result = Message{LDAP_RES_SEARCH_RESULT}

struct MessageIterator{T<:AbstractMessage}
    chain::MessageChain
end

#=
typedef struct berval {
    ber_len_t   bv_len; # long
    char        *bv_val;
} BerValue;
=#
struct Berval
    len::Clong
    val::Ptr{UInt8}
end

#=
typedef struct ldapcontrol {
    char *          ldctl_oid;          /* numericoid of control */
    struct berval   ldctl_value;        /* encoded value of control */
    char            ldctl_iscritical;   /* criticality */
} LDAPControl;
=#
struct LDAPControl
    oid::Cstring
    value::Berval
    iscritical::Cchar
end

struct Control
    oid::String
    value::String
    iscritical::Char
end

struct ParsedResult
    errcodep::Cint
    matched::String
    error_message::String
    referrals::Vector{String}
    controls::Vector{Control}
end

mutable struct AttributesIterator
    entry::Entry
    first_ber_element::Ptr{Cvoid}
    first_attribute_name::Union{Nothing, String}

    function AttributesIterator(entry::Entry, first_ber_element::Ptr{Cvoid}, attr::String)
        result = new(entry, first_ber_element, attr)
        if first_ber_element != C_NULL
            finalizer(destroy!, result)
        end
        return result
    end
end


text:
using ConstrainedDynamics

# Parameters
ex = [1.;0.;0.]

l1 = 1.0
l2 = 1.0# sqrt(2)/2
x, y = .1, .1

vert11 = [0.;0.;l1 / 2]
vert12 = -vert11

vert21 = [0.;0.;l2 / 2]
vert22 = -vert21

verts = [
    [vert11]
    [vert12]
    [vert21]
    [vert22]
]

# Initial orientation
phi1 = pi / 5
q1 = UnitQuaternion(RotX(phi1))

# Links
origin = Origin{Float64}()

link1 = Box(x, y, l1, l1, color = RGBA(1., 1., 0.))
link2 = Box(x, y, l2, l2, color = RGBA(1., 1., 0.))
link3 = Box(x, y, l1, l1, color = RGBA(1., 0., 0.))
link4 = Box(x, y, l2, l2, color = RGBA(1., 0., 0.))

link5 = deepcopy(link1)
link6 = deepcopy(link2)
link7 = deepcopy(link3)
link8 = deepcopy(link4)

# Constraints
function fourbar(links, vertices, axis)
    j1 = EqualityConstraint(Revolute(links[1], links[2], axis; p2=vertices[1]))
    j2 = EqualityConstraint(Revolute(links[2], links[3], axis; p1=vertices[2], p2=vertices[3]), Cylindrical(links[2], links[4], axis; p1=vertices[1], p2=vertices[1]))
    j3 = EqualityConstraint(Revolute(links[4], links[5], axis; p1=vertices[2], p2=vertices[3]))
    j4 = EqualityConstraint(Revolute(links[3], links[5], axis; p1=vertices[4], p2=vertices[4]))

    return j1, j2, j3, j4
end

function initfourbar!(mechanism, links, vertices, Δq1, Δq2)
    setPosition!(links[1], links[2], p2 = vertices[1], Δq = Δq1)
    setPosition!(links[2], links[3], p1 = vertices[2], p2 = vertices[3], Δq = inv(Δq2))
    setPosition!(links[2], links[4], p1 = vertices[1], p2 = vertices[1], Δq = inv(Δq2))
    setPosition!(links[4], links[5], p1 = vertices[2], p2 = vertices[3], Δq = Δq2)

    return
end

links = [link1; link2; link3; link4; link5; link6; link7; link8]
constraints = [fourbar([origin;links[1:4]], verts, ex)...,fourbar([links[4];links[5:8]], verts, ex)...]

mech = Mechanism(origin, links, constraints, InequalityConstraint{Float64}[])

initfourbar!(mech,[origin;links[1:4]],verts,q1,q1)
initfourbar!(mech,[links[4];links[5:8]],verts,q1 / q1,q1)


text:
using Documenter, Fermi

Basics = "Basics" => "index.md"
Methods = "Methods" => ["hartreefock.md", "mp.md", "cc.md"]
Index = "Index" => "indice.md"

PAGES = [
    Basics,
    Methods,
    Index
]

makedocs(
    sitename="Fermi.jl",
    authors = "Gustavo Aroeira",
    format = Documenter.HTML(
        sidebar_sitename = false
    ),
    pages = PAGES
)

deploydocs(
    repo = "github.com/FermiQC/Fermi.jl.git",
)

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
using ContactStructures

plotlyjs()
plot_along_xy_grid(-2:0.5:2,-2:0.5:2,0,standard_structure;  markersize = 2.5,ε = 0.25)
plot!(xlims=(-3, 3), ylims=(-3,3), zlims=(-3,3),legend=false)


text:
function solve()
    sx, sy, gx, gy = [parse(Int, x) for x in split(readline())]
    sx - sy*(sx-gx)/(sy+gy)
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
# cfitsio.jl - C-style interface to CFITSIO functions:
#
# - Function names closely mirror the C interface (e.g., `fits_open_file()`).
# - Functions operate on `FITSFile`, a thin wrapper for `fitsfile` C struct
#   (`FITSFile` has concept of "current HDU", as in CFITSIO).
# - Note that the wrapper functions *do* check the return status from CFITSIO
#   and throw an error with the appropriate message.
#
#
# Syntax compatibility notes:
#
# - `convert(Vector{T}, x)` can be changed to `Vector{T}(x)` once
#   v0.3 is no longer supported, or can be changed to
#   `@compat Vector{T}(x)` once syntax support is in Compat.
# - `convert(Cint, x)` can be changed to `Cint(x)` once v0.3 is not supported
#   or `@compat Cint(x)` once syntax support is in Compat.
# - Once v0.3 is no longer supported, the type `Nothing` should be changed to
#   `Void` (this is the type of `nothing`).
#
# The following table gives the correspondances between CFITSIO "types",
# the BITPIX keyword and Julia types.
#
#     -------------------------------------------------
#     CODE  CFISTIO         Julia     Comments
#     -------------------------------------------------
#           int             Cint
#           long            Clong
#           LONGLONG        Int64     64-bit integer
#     -------------------------------------------------
#     -------- FITS BITPIX ----------------------------
#        8  BYTE_IMG        Uint8
#       16  SHORT_IMG       Int16
#       32  LONG_IMG        Int32
#       64  LONGLONG_IMG    Int64
#      -32  FLOAT_IMG       Float32
#      -64  DOUBLE_IMG      Float64
#      -------- cfitsio "aliases" ---------------------
#       10  SBYTE_IMG       Int8     written as: BITPIX = 8, BSCALE = 1,
#                                                BZERO = -128
#       20  USHORT_IMG      Uint16   written as: BITPIX = 16, BSCALE = 1,
#                                                BZERO = 32768
#       40  LONG_IMG        Uint32   written as: BITPIX = 32, BSCALE = 1,
#                                                BZERO = 2147483648
#     -------------------------------------------------
#     -------- FITS TABLE DATA TYPES ------------------
#        1  TBIT
#       11  TBYTE           Cuchar = Uint8
#       12  TSBYTE          Cchar = Int8
#       14  TLOGICAL        Bool
#       16  TSTRING         ASCIIString
#       20  TUSHORT         Cushort
#       21  TSHORT          Cshort
#       30  TUINT           Cuint
#       31  TINT            Cint
#       40  TULONG          Culong
#       41  TLONG           Clong
#       42  TFLOAT          Cfloat
#       81  TLONGLONG       Int64
#       82  TDOUBLE         Cdouble
#       83  TCOMPLEX        Complex{Cfloat}
#      163  TDBLCOMPLEX     Complex{Cdouble}
#     -------------------------------------------------
#

const bitpix_to_type = Dict{Cint, DataType}()
for (T, code) in ((Uint8,     8), # BYTE_IMG
                  (Int16,    16), # SHORT_IMG
                  (Int32,    32), # LONG_IMG
                  (Int64,    64), # LONGLONG_IMG
                  (Float32, -32), # FLOAT_IMG
                  (Float64, -64), # DOUBLE_IMG
                  (Int8,     10), # SBYTE_IMG
                  (Uint16,   20), # USHORT_IMG
                  (Uint32,   40)) # ULONG_IMG
    local value = convert(Cint, code)
    @eval begin
        bitpix_to_type[$value] = $T
        _cfitsio_bitpix(::Type{$T}) = $value
    end
end

for (T, code) in ((Uint8,       11),
                  (Int8,        12),
                  (Bool,        14),
                  (ASCIIString, 16),
                  (Cushort,     20),
                  (Cshort,      21),
                  (Cuint,       30),
                  (Cint,        31),
                  (Culong,      40),
                  (Clong,       41),
                  (Float32,     42),
                  (Int64,       81),
                  (Float64,     82),
                  (Complex64,   83),
                  (Complex128, 163))
    @eval _cfitsio_datatype(::Type{$T}) = convert(Cint, $code)
end

type FITSFile
    ptr::Ptr{Void}

    function FITSFile(ptr::Ptr{Void})
        f = new(ptr)
        finalizer(f, fits_close_file)
        f
    end
end

# -----------------------------------------------------------------------------
# errors and assertions

function fits_assert_open(f::FITSFile)
    if f.ptr == C_NULL
        error("attempt to access closed FITS file")
    end
end

function fits_get_errstatus(status::Cint)
    msg = Array(Uint8, 31)
    ccall((:ffgerr,libcfitsio), Void, (Cint,Ptr{Uint8}), status, msg)
    bytestring(pointer(msg))
end

function fits_assert_ok(status::Cint)
    if status != 0
        error(fits_get_errstatus(status))
    end
end

# -----------------------------------------------------------------------------
# file access & info functions

function fits_create_file(filename::String)
    ptr = Array(Ptr{Void}, 1)
    status = Cint[0]
    ccall((:ffinit,libcfitsio), Cint, (Ptr{Ptr{Void}},Ptr{Uint8},Ptr{Cint}),
          ptr, bytestring(filename), status)
    fits_assert_ok(status[1])
    FITSFile(ptr[1])
end

fits_clobber_file(filename::String) = fits_create_file("!"*filename)

for (a,b) in ((:fits_open_data, "ffdopn"),
              (:fits_open_file, "ffopen"),
              (:fits_open_image,"ffiopn"),
              (:fits_open_table,"fftopn"))
    @eval begin
        function ($a)(filename::String, mode::Integer=0)
            ptr = Array(Ptr{Void}, 1)
            status = Cint[0]
            ccall(($b,libcfitsio), Cint,
                  (Ptr{Ptr{Void}},Ptr{Uint8},Cint,Ptr{Cint}),
                  ptr, bytestring(filename), mode, status)
            fits_assert_ok(status[1])
            FITSFile(ptr[1])
        end
    end
end

for (a,b) in ((:fits_close_file, "ffclos"),
              (:fits_delete_file,"ffdelt"))
    @eval begin
        function ($a)(f::FITSFile)

            # fits_close_file() is called during garbage collection, but file
            # may already be closed by user, so we need to check if it is open.
            if f.ptr != C_NULL
                status = Cint[0]
                ccall(($b,libcfitsio), Cint,
                      (Ptr{Void},Ptr{Cint}),
                      f.ptr, status)
                fits_assert_ok(status[1])
                f.ptr = C_NULL
            end
        end
    end
end

close(f::FITSFile) = fits_close_file(f)

function fits_file_name(f::FITSFile)
    value = Array(Uint8, 1025)
    status = Cint[0]
    ccall((:ffflnm,libcfitsio), Cint,
          (Ptr{Void},Ptr{Uint8},Ptr{Cint}),
          f.ptr, value, status)
    fits_assert_ok(status[1])
    bytestring(pointer(value))
end

function fits_file_mode(f::FITSFile)
    result = Cint[0]
    status = Cint[0]
    ccall(("ffflmd", libcfitsio), Cint, (Ptr{Void}, Ptr{Cint}, Ptr{Cint}),
          f.ptr, result, status)
    fits_assert_ok(status[1])
    result[1]
end

# -----------------------------------------------------------------------------
# header access functions

function fits_get_hdrspace(f::FITSFile)
    keysexist = Cint[0]
    morekeys = Cint[0]
    status = Cint[0]
    ccall((:ffghsp,libcfitsio), Cint,
        (Ptr{Void},Ptr{Cint},Ptr{Cint},Ptr{Cint}),
        f.ptr, keysexist, morekeys, status)
    fits_assert_ok(status[1])
    (keysexist[1], morekeys[1])
end

function fits_read_keyword(f::FITSFile, keyname::ASCIIString)
    value = Array(Uint8, 71)
    comment = Array(Uint8, 71)
    status = Cint[0]
    ccall((:ffgkey,libcfitsio), Cint,
        (Ptr{Void},Ptr{Uint8},Ptr{Uint8},Ptr{Uint8},Ptr{Cint}),
        f.ptr, bytestring(keyname), value, comment, status)
    fits_assert_ok(status[1])
    bytestring(pointer(value)), bytestring(pointer(comment))
end

function fits_read_record(f::FITSFile, keynum::Integer)
    card = Array(Uint8, 81)
    status = Cint[0]
    ccall((:ffgrec,libcfitsio), Cint,
        (Ptr{Void},Cint,Ptr{Uint8},Ptr{Cint}),
        f.ptr, keynum, card, status)
    fits_assert_ok(status[1])
    bytestring(pointer(card))
end

function fits_read_keyn(f::FITSFile, keynum::Integer)
    keyname = Array(Uint8, 9)
    value = Array(Uint8, 71)
    comment = Array(Uint8, 71)
    status = Cint[0]
    ccall((:ffgkyn,libcfitsio), Cint,
        (Ptr{Void},Cint,Ptr{Uint8},Ptr{Uint8},Ptr{Uint8},Ptr{Cint}),
        f.ptr, keynum, keyname, value, comment, status)
    fits_assert_ok(status[1])
    (bytestring(pointer(keyname)), bytestring(pointer(value)),
     bytestring(pointer(comment)))
end

function fits_write_key(f::FITSFile, keyname::ASCIIString,
                        value::Union(FloatingPoint,ASCIIString),
                        comment::ASCIIString)
    cvalue = isa(value,ASCIIString) ?  bytestring(value) :
             isa(value,Bool) ? Cint[value] : [value]
    status = Cint[0]
    ccall((:ffpky,libcfitsio), Cint,
        (Ptr{Void},Cint,Ptr{Uint8},Ptr{Uint8},Ptr{Uint8},Ptr{Cint}),
        f.ptr, _cfitsio_datatype(typeof(value)), bytestring(keyname),
        cvalue, bytestring(comment), status)
    fits_assert_ok(status[1])
end

function fits_write_comment(f::FITSFile, comment::ASCIIString)
    status = Cint[0]
    ccall((:ffpcom, libcfitsio), Cint, (Ptr{Void}, Ptr{Uint8}, Ptr{Cint}),
          f.ptr, bytestring(comment), status)
    fits_assert_ok(status[1])
end

function fits_write_history(f::FITSFile, history::ASCIIString)
    status = Cint[0]
    ccall((:ffphis, libcfitsio), Cint, (Ptr{Void}, Ptr{Uint8}, Ptr{Cint}),
          f.ptr, bytestring(history), status)
    fits_assert_ok(status[1])
end

# update key: if already present, update it, otherwise add it.
for (a,T,S) in (("ffukys", :ASCIIString, :(Ptr{Uint8})),
                ("ffukyl", :Bool,        :Cint),
                ("ffukyj", :Integer,     :Int64))
    @eval begin
        function fits_update_key(f::FITSFile, key::ASCIIString, value::$T,
                                 comment::Union(ASCIIString, Ptr{Void})=C_NULL)
            status = Cint[0]
            ccall(($a, libcfitsio), Cint,
                  (Ptr{Void}, Ptr{Uint8}, $S, Ptr{Uint8}, Ptr{Cint}),
                  f.ptr, key, value, comment, status)
            fits_assert_ok(status[1])
        end
    end
end

function fits_update_key(f::FITSFile, key::ASCIIString, value::FloatingPoint,
                         comment::Union(ASCIIString, Ptr{Void})=C_NULL)
    status = Cint[0]
    ccall(("ffukyd", libcfitsio), Cint,
          (Ptr{Void}, Ptr{Uint8}, Cdouble, Cint, Ptr{Uint8}, Ptr{Cint}),
          f.ptr, key, value, -15, comment, status)
    fits_assert_ok(status[1])
end

function fits_update_key(f::FITSFile, key::ASCIIString, value::Nothing,
                         comment::Union(ASCIIString, Ptr{Void})=C_NULL)
    status = Cint[0]
    ccall(("ffukyu", libcfitsio), Cint,
          (Ptr{Void}, Ptr{Uint8}, Ptr{Uint8}, Ptr{Cint}),
          f.ptr, key, comment, status)
    fits_assert_ok(status[1])
end

function fits_write_record(f::FITSFile, card::ASCIIString)
    status = Cint[0]
    ccall((:ffprec,libcfitsio), Cint,
        (Ptr{Void},Ptr{Uint8},Ptr{Cint}),
        f.ptr, bytestring(card), status)
    fits_assert_ok(status[1])
end

function fits_delete_record(f::FITSFile, keynum::Integer)
    status = Cint[0]
    ccall((:ffdrec,libcfitsio), Cint,
        (Ptr{Void},Cint,Ptr{Cint}),
        f.ptr, keynum, status)
    fits_assert_ok(status[1])
end

function fits_delete_key(f::FITSFile, keyname::ASCIIString)
    status = Cint[0]
    ccall((:ffdkey,libcfitsio), Cint,
        (Ptr{Void},Ptr{Uint8},Ptr{Cint}),
        f.ptr, bytestring(keyname), status)
    fits_assert_ok(status[1])
end

function fits_hdr2str(f::FITSFile, nocomments::Bool=false)
    status = Cint[0]
    header = Array(Ptr{Uint8}, 1)
    nkeys = Cint[0]
    ccall((:ffhdr2str, libcfitsio), Cint,
          (Ptr{Void}, Cint, Ptr{Ptr{Uint8}}, Cint,
           Ptr{Ptr{Uint8}}, Ptr{Cint}, Ptr{Cint}),
          f.ptr, nocomments, &C_NULL, 0, header, nkeys, status)
    result = bytestring(header[1])

    # free header pointer allocated by cfitsio (result is a copy)
    ccall((:fffree, libcfitsio), Ptr{Cint},
          (Ptr{Uint8}, Ptr{Cint}),
          header[1], status)
    fits_assert_ok(status[1])
    result
end

# -----------------------------------------------------------------------------
# HDU info functions and moving the current HDU

function hdu_int_to_type(hdu_type_int)
    if hdu_type_int == 0
        return :image_hdu
    elseif hdu_type_int == 1
        return :ascii_table
    elseif hdu_type_int == 2
        return :binary_table
    end

    :unknown
end

for (a,b) in ((:fits_movabs_hdu,"ffmahd"),
              (:fits_movrel_hdu,"ffmrhd"))
    @eval begin
        function ($a)(f::FITSFile, hduNum::Integer)
            hdu_type = Cint[0]
            status = Cint[0]
            ccall(($b,libcfitsio), Cint,
                  (Ptr{Void}, Cint, Ptr{Cint}, Ptr{Cint}),
                  f.ptr, hduNum, hdu_type, status)
            fits_assert_ok(status[1])
            hdu_int_to_type(hdu_type[1])
        end
    end
end

function fits_movnam_hdu(f::FITSFile, extname::ASCIIString, extver::Integer=0,
                         hdu_type::Integer=-1)
    status = Cint[0]
    ccall((:ffmnhd,libcfitsio), Cint,
          (Ptr{Void}, Cint, Ptr{Uint8}, Cint, Ptr{Cint}),
          f.ptr, hdu_type, bytestring(extname), extver, status)
    fits_assert_ok(status[1])
end

function fits_get_hdu_num(f::FITSFile)
    hdunum = Cint[0]
    ccall((:ffghdn,libcfitsio), Cint,
          (Ptr{Void}, Ptr{Cint}),
          f.ptr, hdunum)
    hdunum[1]
end

function fits_get_hdu_type(f::FITSFile)
    hdutype = Cint[0]
    status = Cint[0]
    ccall((:ffghdt, libcfitsio), Cint,
          (Ptr{Void}, Ptr{Cint}, Ptr{Cint}),
          f.ptr, hdutype, status)
    fits_assert_ok(status[1])
    hdu_int_to_type(hdutype[1])
end

# -----------------------------------------------------------------------------
# image HDU functions

for (a, b) in ((:fits_get_img_type,      "ffgidt"),
               (:fits_get_img_equivtype, "ffgiet"),
               (:fits_get_img_dim,       "ffgidm"))
    @eval function ($a)(f::FITSFile)
        result = Cint[0]
        status = Cint[0]
        ccall(($b, libcfitsio), Cint, (Ptr{Void}, Ptr{Cint}, Ptr{Cint}),
              f.ptr, result, status)
        fits_assert_ok(status[1])
        result[1]
    end
end

function fits_get_img_size(f::FITSFile)
    ndim = fits_get_img_dim(f)
    naxes = Array(Clong, ndim)
    status = Cint[0]
    ccall((:ffgisz, libcfitsio), Cint,
        (Ptr{Void}, Cint, Ptr{Clong}, Ptr{Cint}),
        f.ptr, ndim, naxes, status)
    fits_assert_ok(status[1])
    naxes
end

function fits_create_img{S<:Integer}(f::FITSFile, t::Type, naxes::Vector{S})
    status = Cint[0]
    ccall((:ffcrimll, libcfitsio), Cint,
          (Ptr{Void}, Cint, Cint, Ptr{Int64}, Ptr{Cint}),
          f.ptr, _cfitsio_bitpix(t), length(naxes),
          convert(Vector{Int64}, naxes), status)
    fits_assert_ok(status[1])
end

function fits_write_pix{S<:Integer,T}(f::FITSFile, fpixel::Vector{S},
                                      nelements::Integer, data::Array{T})
    status = Cint[0]
    ccall((:ffppxll, libcfitsio), Cint,
          (Ptr{Void}, Cint, Ptr{Int64}, Int64, Ptr{Void}, Ptr{Cint}),
          f.ptr, _cfitsio_datatype(T), convert(Vector{Int64}, fpixel),
          nelements, data, status)
    fits_assert_ok(status[1])
end

function fits_write_pix(f::FITSFile, data::Array)
    fits_write_pix(f, ones(Int64, length(size(data))), length(data), data)
end

function fits_read_pix{S<:Integer,T}(f::FITSFile, fpixel::Vector{S},
                                     nelements::Int, nullval::T,
                                     data::Array{T})
    anynull = Cint[0]
    status = Cint[0]
    ccall((:ffgpxv, libcfitsio), Cint,
          (Ptr{Void}, Cint, Ptr{Clong}, Int64, Ptr{Void}, Ptr{Void},
           Ptr{Cint}, Ptr{Cint}),
          f.ptr, _cfitsio_datatype(T), convert(Vector{Int64}, fpixel),
          nelements, &nullval, data, anynull, status)
    fits_assert_ok(status[1])
    anynull[1]
end

function fits_read_pix{S<:Integer,T}(f::FITSFile, fpixel::Vector{S},
                                     nelements::Int, data::Array{T})
    anynull = Cint[0]
    status = Cint[0]
    ccall((:ffgpxv, libcfitsio), Cint,
          (Ptr{Void}, Cint, Ptr{Clong}, Int64, Ptr{Void}, Ptr{Void},
           Ptr{Cint}, Ptr{Cint}),
          f.ptr, _cfitsio_datatype(T), convert(Vector{Int64}, fpixel),
          nelements, C_NULL, data, anynull, status)
    fits_assert_ok(status[1])
    anynull[1]
end

function fits_read_pix(f::FITSFile, data::Array)
    fits_read_pix(f, ones(Int64,length(size(data))), length(data), data)
end

function fits_read_subset{S1<:Integer,S2<:Integer,S3<:Integer,T}(
             f::FITSFile, fpixel::Vector{S1}, lpixel::Vector{S2},
             inc::Vector{S3}, data::Array{T})
    anynull = Cint[0]
    status = Cint[0]
    ccall((:ffgsv, libcfitsio), Cint,
          (Ptr{Void}, Cint, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Ptr{Void},
           Ptr{Void}, Ptr{Cint}, Ptr{Cint}),
          f.ptr, _cfitsio_datatype(T),
          convert(Vector{Clong}, fpixel),
          convert(Vector{Clong}, lpixel),
          convert(Vector{Clong}, inc),
          C_NULL, data, anynull, status)
    fits_assert_ok(status[1])
    anynull[1]
end

function fits_copy_image_section(fin::FITSFile, fout::FITSFile,
                                 section::ASCIIString)
    status = Cint[0]
    ccall((:fits_copy_image_section,libcfitsio), Cint,
          (Ptr{Void}, Ptr{Void}, Ptr{Uint8}, Ptr{Cint}),
          fin.ptr, fout.ptr, bytestring(section), status)
    fits_assert_ok(status[1])
end

# -----------------------------------------------------------------------------
# ASCII/binary table HDU functions

# The three fields are: ttype, tform, tunit (CFITSIO's terminology)
typealias ColumnDef @compat Tuple{ASCIIString, ASCIIString, ASCIIString}

for (a,b) in ((:fits_create_binary_tbl, 2),
              (:fits_create_ascii_tbl,  1))
    @eval begin
        function ($a)(f::FITSFile, numrows::Integer,
                      coldefs::Array{ColumnDef}, extname::ASCIIString)

            ntype = length(coldefs)
            ttype = map((x) -> pointer(x[1].data), coldefs)
            tform = map((x) -> pointer(x[2].data), coldefs)
            tunit = map((x) -> pointer(x[3].data), coldefs)
            status = Cint[0]

            ccall(("ffcrtb", libcfitsio), Cint,
                  (Ptr{Void}, Cint, Int64, Cint,
                   Ptr{Ptr{Uint8}}, Ptr{Ptr{Uint8}},
                   Ptr{Ptr{Uint8}}, Ptr{Uint8}, Ptr{Cint}),
                  f.ptr, $b, numrows, ntype,
                  ttype, tform, tunit, bytestring(extname),
                  status)
            fits_assert_ok(status[1])
        end
    end
end

for (a,b,T) in ((:fits_get_num_cols,  "ffgncl",  :Cint),
                (:fits_get_num_hdus,  "ffthdu",  :Cint),
                (:fits_get_num_rows,  "ffgnrw",  :Clong),
                (:fits_get_num_rowsll,"ffgnrwll",:Int64),
                (:fits_get_rowsize,   "ffgrsz",  :Clong))
    @eval begin
        function ($a)(f::FITSFile)
            result = $T[0]
            status = Cint[0]
            ccall(($b,libcfitsio), Cint,
                  (Ptr{Void}, Ptr{$T}, Ptr{Cint}),
                  f.ptr, result, status)
            fits_assert_ok(status[1])
            result[1]
        end
    end
end

# The function `fits_read_tdim()` returns the dimensions of a table column in a
# binary table. Normally this information is given by the TDIMn keyword, but if
# this keyword is not present then this routine returns `[r]` with `r` equals
# to the repeat count in the TFORM keyword.
let fn, T, ffgtdm, ffgtcl, ffeqty
    if  promote_type(Int, Clong) == Clong
        T = Clong
        ffgtdm = "ffgtdm"
        ffgtcl = "ffgtcl"
        ffeqty = "ffeqty"
    else
        T = Int64
        ffgtdm = "ffgtdmll"
        ffgtcl = "ffgtclll"
        ffeqty = "ffeqtyll"
    end
    @eval begin
        function fits_get_coltype(ff::FITSFile, colnum::Integer)
            typecode = Cint[0]
            repcnt = $T[0]
            width = $T[0]
            status = Cint[0]
            ccall(($ffgtcl,libcfitsio), Cint,
                  (Ptr{Void}, Cint, Ptr{Cint}, Ptr{$T}, Ptr{$T}, Ptr{Cint}),
                  ff.ptr, colnum, typecode, repcnt, width, status)
            fits_assert_ok(status[1])
            return @compat Int(typecode[1]), Int(repcnt[1]), Int(width[1])
        end
        function fits_get_eqcoltype(ff::FITSFile, colnum::Integer)
            typecode = Cint[0]
            repcnt = $T[0]
            width = $T[0]
            status = Cint[0]
            ccall(($ffeqty,libcfitsio), Cint,
                  (Ptr{Void}, Cint, Ptr{Cint}, Ptr{$T}, Ptr{$T}, Ptr{Cint}),
                  ff.ptr, colnum, typecode, repcnt, width, status)
            fits_assert_ok(status[1])
            return @compat Int(typecode[1]), Int(repcnt[1]), Int(width[1])
        end
        function fits_read_tdim(ff::FITSFile, colnum::Integer)
            naxes = Array($T, 99) # 99 is the maximum allowed number of axes
            naxis = Cint[0]
            status = Cint[0]
            ccall(($ffgtdm,libcfitsio), Cint,
                  (Ptr{Void}, Cint, Cint, Ptr{Cint}, Ptr{$T}, Ptr{Cint}),
                  ff.ptr, colnum, length(naxes), naxis, naxes, status)
            fits_assert_ok(status[1])
            return naxes[1:naxis[1]]
        end
    end
end

function fits_read_col(f::FITSFile,
                       colnum::Integer,
                       firstrow::Integer,
                       firstelem::Integer,
                       data::Array{ASCIIString})

    # Make sure there is enough room for each string
    _, repcount, width = fits_get_coltype(f, colnum)
    for i in 1:length(data)
        # We need to call `repeat' N times in order to allocate N
        # strings
        data[i] = repeat(" ", repcount)
    end

    # Call the CFITSIO function
    anynull = Cint[0]
    status = Cint[0]
    ccall((:ffgcvs, libcfitsio), Cint,
          (Ptr{Void}, Cint, Int64, Int64, Int64,
           Ptr{Uint8}, Ptr{Ptr{Uint8}}, Ptr{Cint}, Ptr{Cint}),
          f.ptr, colnum, firstrow, firstelem, length(data),
          "", data, anynull, status)
    fits_assert_ok(status[1])

    # Truncate the strings to the first NULL character (if present)
    for idx in 1:length(data)
        zeropos = search(data[idx], '\0')
        if zeropos >= 1
            data[idx] = (data[idx])[1:(zeropos-1)]
        end
    end
end

function fits_read_col{T}(f::FITSFile,
                          colnum::Integer,
                          firstrow::Integer,
                          firstelem::Integer,
                          data::Array{T})
    anynull = Cint[0]
    status = Cint[0]
    ccall((:ffgcv,libcfitsio), Cint,
          (Ptr{Void}, Cint, Cint, Int64, Int64, Int64,
           Ptr{T}, Ptr{T}, Ptr{Cint}, Ptr{Cint}),
          f.ptr, _cfitsio_datatype(T), colnum,
          firstrow, firstelem, length(data),
          T[0], data, anynull, status)
    fits_assert_ok(status[1])
end

function fits_write_col(f::FITSFile,
                        colnum::Integer,
                        firstrow::Integer,
                        firstelem::Integer,
                        data::Array{ASCIIString})
    status = Cint[0]
    ccall((:ffpcls, libcfitsio), Cint,
          (Ptr{Void}, Cint, Int64, Int64, Int64,
           Ptr{Ptr{Uint8}}, Ptr{Cint}),
          f.ptr, colnum, firstrow, firstelem, length(data),
          data, status)
    fits_assert_ok(status[1])
end

function fits_write_col{T}(f::FITSFile,
                           colnum::Integer,
                           firstrow::Integer,
                           firstelem::Integer,
                           data::Array{T})
    status = Cint[0]
    ccall((:ffpcl, libcfitsio), Cint,
          (Ptr{Void}, Cint, Cint, Int64, Int64, Int64,
           Ptr{T}, Ptr{Cint}),
          f.ptr, _cfitsio_datatype(T), colnum,
          firstrow, firstelem, length(data),
          data, status)
    fits_assert_ok(status[1])
end

for (a,b) in ((:fits_insert_rows, "ffirow"),
              (:fits_delete_rows, "ffdrow"))
    @eval begin
        function ($a)(f::FITSFile, firstrow::Integer, nrows::Integer)
            status = Cint[0]
            ccall(($b,libcfitsio), Cint,
                  (Ptr{Void}, Int64, Int64, Ptr{Cint}),
                  f.ptr, firstrow, nrows, status)
            fits_assert_ok(status[1])
        end
    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libudfread_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libudfread")
JLLWrappers.@generate_main_file("libudfread", UUID("037e6697-03b9-52b7-b841-7aee0d773eb5"))
end  # module libudfread_jll


text:
function findExtrema(signal::AbstractArray)
    maxima  = [x.I[1] for x in Images.findlocalmaxima(signal, 1)]
    maxvals = signal[maxima]

    minima  = [x.I[1] for x in Images.findlocalminima(signal, 1)]
    minvals = signal[minima]

    ## make sure at length extrema > 1
    if isempty(maxvals)
        mm,tt = findmax(tempy)
        push!(maxima,tt)
        push!(maxvals,mm)
    end

    if isempty(minvals)
        mm,tt = findmin(tempy)
        push!(minima,tt)
        push!(minvals,mm)
    end

    return maxima,minima,maxvals,minvals
end

max_stars_repo_name: alexmorley/EmpiricalModeDecomposition

max_issues_repo_name: alexmorley/EmpiricalModeDecomposition

max_forks_repo_name: alexmorley/EmpiricalModeDecomposition