text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -3.357520876571253 <= q <= 3.7351340364874854)

                     add_NL_constraint(m, :(gelu(0.618782057809276 + -0.4828154761503023*gelu(0.8842828756523886 + -0.6089179668773923*gelu(0.019018219753813348 + 0.8244023518746615*gelu(-0.7862452642307147 + -0.6499989200319587*$(x[1]) + -0.39995929749666725*$(x[2])) + 0.3296833580573324*gelu(0.2947581350949262 + -0.3447038526107806*$(x[1]) + -0.5996590697259889*$(x[2]))) + 0.44138625638024864*gelu(-0.9126485111761089 + 0.9176394232569307*gelu(-0.7862452642307147 + -0.6499989200319587*$(x[1]) + -0.39995929749666725*$(x[2])) + 0.371990627187599*gelu(0.2947581350949262 + -0.3447038526107806*$(x[1]) + -0.5996590697259889*$(x[2])))) + 0.3197591788003944*gelu(0.7251254505916935 + -0.923704860259535*gelu(0.019018219753813348 + 0.8244023518746615*gelu(-0.7862452642307147 + -0.6499989200319587*$(x[1]) + -0.39995929749666725*$(x[2])) + 0.3296833580573324*gelu(0.2947581350949262 + -0.3447038526107806*$(x[1]) + -0.5996590697259889*$(x[2]))) + 0.28767280168552967*gelu(-0.9126485111761089 + 0.9176394232569307*gelu(-0.7862452642307147 + -0.6499989200319587*$(x[1]) + -0.39995929749666725*$(x[2])) + 0.371990627187599*gelu(0.2947581350949262 + -0.3447038526107806*$(x[1]) + -0.5996590697259889*$(x[2]))))) + gelu(-0.0992210488643468 + 0.772244110024173*gelu(0.8842828756523886 + -0.6089179668773923*gelu(0.019018219753813348 + 0.8244023518746615*gelu(-0.7862452642307147 + -0.6499989200319587*$(x[1]) + -0.39995929749666725*$(x[2])) + 0.3296833580573324*gelu(0.2947581350949262 + -0.3447038526107806*$(x[1]) + -0.5996590697259889*$(x[2]))) + 0.44138625638024864*gelu(-0.9126485111761089 + 0.9176394232569307*gelu(-0.7862452642307147 + -0.6499989200319587*$(x[1]) + -0.39995929749666725*$(x[2])) + 0.371990627187599*gelu(0.2947581350949262 + -0.3447038526107806*$(x[1]) + -0.5996590697259889*$(x[2])))) + 0.07492471407403789*gelu(0.7251254505916935 + -0.923704860259535*gelu(0.019018219753813348 + 0.8244023518746615*gelu(-0.7862452642307147 + -0.6499989200319587*$(x[1]) + -0.39995929749666725*$(x[2])) + 0.3296833580573324*gelu(0.2947581350949262 + -0.3447038526107806*$(x[1]) + -0.5996590697259889*$(x[2]))) + 0.28767280168552967*gelu(-0.9126485111761089 + 0.9176394232569307*gelu(-0.7862452642307147 + -0.6499989200319587*$(x[1]) + -0.39995929749666725*$(x[2])) + 0.371990627187599*gelu(0.2947581350949262 + -0.3447038526107806*$(x[1]) + -0.5996590697259889*$(x[2]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# Tv: sparse matrix value type
# Ti: sparse matrix ordinal space
function readadj(Tv::Type, Ti::Type, alpha::Float64, path::String)
  io = open(path, "r")
  header = readline(io)
  edges, order, max_ordinal = map(int64, split(header))
  readline(io)

  I, J, V = Ti[], Ti[], Tv[]
  vertices, sources  = Set{Ti}(), Set{Ti}()

  sinks = 0
  outdegree = 0
  for line in EachLine(io)
    if sinks == outdegree
      separator = search(line, ' ')
      source = int32(line[1:separator - 1])
      outdegree = int32(line[separator:end])
      add!(sources, source)
      add!(vertices, source)

      sinks = 0
    else
      sink = int32(line)
      push!(I, source)
      push!(J, sink)
      push!(V, alpha / outdegree)
      add!(vertices, sink)

      sinks += 1
    end
  end

  close(io)
  M::SparseMatrixCSC = sparse(I, J, V, max_ordinal, max_ordinal)

  absorbing = Ti[]
  for i = keys(setdiff(vertices, sources).dict)
    push!(absorbing, i)
  end

  return M, absorbing, order, max_ordinal
end

# Tv: sparse matrix value type
# Ti: sparse matrix ordinal space
function fastadj(Tv::Type, Ti::Type, alpha::Float64, path::String, bufsize::Int64)
  io = open(path, "r")
  header = readline(io)
  edges, order, max_ordinal = map(int64, split(header))
  readline(io)

  I, J, V = Array(Ti, edges), Array(Ti, edges), Array(Tv, edges)
  buffer = Array(Uint8, bufsize)
  ordinal = 0
  sinks = 0
  outdegree = 0
  edge = 0

  while !eof(io)
    fill!(buffer, 0)
    try
      read(io, buffer)
    end

    for i = 1:length(buffer)
      if buffer[i] == '\n'
        if sinks == 0
          outdegree = ordinal
        else
          sink = ordinal
          edge += 1
          I[edge] = source
          J[edge] = sink
          V[edge] = alpha / outdegree
        end

        sinks += 1
        ordinal = 0
      elseif buffer[i] == ' '
        source = ordinal
        sinks = 0
        ordinal = 0
      else
        ordinal = (ordinal * 10) + (buffer[i] - 48)
      end
    end
  end

  close(io)
  M::SparseMatrixCSC = sparse(I, J, V, max_ordinal, max_ordinal)

  return M, order, max_ordinal
end

function fastadj(Tv::Type, Ti::Type, alpha::Float64, path::String)
  fastadj(Tv, Ti, alpha, path, 8192)
end

function count_edges(path::String)
  lines = countlines(path)
  spaces = countlines(path, ' ') - 1 # broken?
  lines - spaces
end

ext: jl


text:
function _svmreader(files::Vector{<:AbstractString}, ::Type{Index}, ::Type{Value}) where {Index,Value}
  labelmask = r"^(?<val>-?\d)"
  featmask = r"(?<idx>\d+):(?<val>\d+\.\d+)"
  row = Index(1)
  rows, cols, vals, labels = Vector{Index}(), Vector{Index}(), Vector{Value}(), Vector{Value}()
  for file in files
    open(file, "r") do io
      for line in eachline(io)
        m = match(labelmask, line)
        push!(labels, parse(Float64, m[:val]))
        for m in eachmatch(featmask, line)
          push!(rows, row)
          push!(cols, parse(Int, m[:idx]))
          push!(vals, parse(Float64, m[:val]))
        end
        row += one(Index)
      end
    end
  end
  return rows, cols, vals, labels
end

function svmreader(files::Vector{<:AbstractString}, ::Type{Index}, ::Type{Value}) where {Index,Value}
  rows, cols, vals, labels = _svmreader(files, Index, Value)
  return sparse(rows, cols, vals), labels
end

function svmreader(files::Vector{<:AbstractString}, m, n, ::Type{Index}, ::Type{Value}) where {Index, Value}
  rows, cols, vals, labels = _svmreader(files, Index, Value)
  return sparse(rows, cols, vals, m, n), labels
end

svmreader(files::Vector{<:AbstractString}) = svmreader(files, Int, Float64)
svmreader(files::Vector{<:AbstractString}, m, n) = svmreader(files, m, n, Int, Float64)

svmreader(file::AbstractString) = svmreader([file])
svmreader(file::AbstractString, ::Type{Index}, ::Type{Value}) where {Index, Value} = svmreader([file], Index, Value)
svmreader(file::AbstractString, m, n) = svmreader([file], m, n)
svmreader(file::AbstractString, m, n, ::Type{Index}, ::Type{Value}) where {Index, Value} = svmreader([file], m, n, Index, Value)


text:
using BAT
using Distributions
using IntervalSets
using TypedTables
using NamedArrays
using LinearAlgebra
using PrettyTables

likelihood = params -> begin
    r1 = logpdf.(Normal(-10.0, 1.2), params.a)
    r2 = logpdf.(Normal(-5.0, 0.8), params.c[1])
    r3 = logpdf.(Normal(12.0, 1.4), params.c[2])

    return LogDVal(r1+r2+r3)
end

prior = BAT.NamedTupleDist(
    a = Normal(-3, 4.5),
    b = -10..10,
    c = [-20.0..20.0, -10..10]
)

posterior = PosteriorDensity(likelihood, prior);

samples, chains = bat_sample(posterior, MCMCSampling(sampler = MetropolisHastings(), nsteps = 10^5));
#samples = bat_sample(posterior, SobolSampler(nsamples = 10^5)).result;

sd = SampledDensity(posterior, samples, generator=BAT.MCMCSampleGenerator(chains))
display(sd)

# # write summary to txt-file
# io = open("sampled_density.txt", "w")
# show(io, "text/plain", sd)
# close(io)
#
# # ### write output to html-file
# io = open("sampled_density.html", "w")
# show(io, "text/html", sd)
# close(io)


text:
# This file is part of Kpax3. License is MIT.

# parameters to test are
# R = [1; 1; 2; 2; 3; 2]
# S = [2; 3; 1; 2]
# C = [2 3 1 2;
#      2 4 1 2;
#      2 3 1 2]

#=
# compute log(normalization constant) as accurate as possible
include("data/partitions.jl")

function Kpax3.computelogp(b::Int,
                     c::Vector{UInt8},
                     n1s::Matrix{Float64},
                     priorC::Kpax3.AminoAcidPriorCol)
  k = size(n1s, 1)
  logp = 0.0

  if c[1] == UInt8(1)
    logp = priorC.logγ[1]

    for g in 1:k
      logp += Kpax3.logmarglik(n1s[g, b], v[g], priorC.A[1, b], priorC.B[1, b])
    end
  elseif c1[1] == UInt8(2)
    logp = priorC.logγ[2]

    for g in 1:k
      logp += Kpax3.logmarglik(n1s[g, b], v[g], priorC.A[2, b], priorC.B[2, b])
    end
  else
    logp = priorC.logγ[3]

    for g in 1:k
      logp += priorC.logω[k][c[g] - 2] + Kpax3.logmarglik(n1s[g, b], v[g], priorC.A[c[g], b], priorC.B[c[g], b])
    end
  end

  logp
end

function addvalue!(b::Int,
                   c::Vector{UInt8},
                   value::Float64,
                   S::Matrix{Float64})
  if c[1] == UInt8(1)
    S[1, b] += value
  elseif c[1] == UInt8(2)
    S[2, b] += value
  else
    S[3, b] += value
  end

  nothing
end

function Kpax3.computelognormconst(ck,
                             k::Int,
                             lumpp::Float64,
                             data::Matrix{UInt8},
                             po::TestPartition,
                             γ::Vector{Float64},
                             r::Float64,
                             priorR::Kpax3.PriorRowPartition)
  (m, n) = size(data)

  priorC = Kpax3.AminoAcidPriorCol(data, γ, r)

  st = po.index[po.k .== k][1]
  en = any(po.k .== k + 1) ? po.index[po.k .== k + 1][1] - 1 : st

  R = zeros(Int, n)
  v = zeros(Float64, k)
  n1s = zeros(Float64, k, m)

  M = -Inf

  g = 0

  p = 0.0
  logprR = 0.0
  logpost = 0.0
  logp = zeros(Float64, m)

  for l in st:en
    fill!(R, 0)
    fill!(v, 0)
    fill!(n1s, 0.0)
    fill!(logp, 0.0)

    for a in 1:n
      g = po.partition[a, l]

      R[a] = g
      v[g] += 1

      for b in 1:m
        n1s[g, b] += float(data[b, a])
      end
    end

    logprR = Kpax3.logdPriorRow(R, priorR)

    for c1 in ck, c2 in ck, c3 in ck, c4 in ck
      logp[1] += Kpax3.computelogp(1, c1, n1s, priorC)
      logp[2] += Kpax3.computelogp(2, c2, n1s, priorC)
      logp[3] += Kpax3.computelogp(3, c3, n1s, priorC)
      logp[4] += Kpax3.computelogp(4, c4, n1s, priorC)

      logpost = logprR + logp[1] + logp[2] + logp[3] + logp[4]

      if logpost > M
        M = logpost
      end

      p += exp(logpost - lumpp)
    end
  end

  (M, p)
end

function lognormconst(cs,
                      data::Matrix{UInt8},
                      po::TestPartition,
                      γ::Vector{Float64},
                      r::Float64,
                      priorR::Kpax3.PriorRowPartition)
  # log unnormalized maximum posterior probability
  lumpp = -Inf

  println("Computing 'lumpp'...")
  for k in 1:size(data, 2)
    println("k = ", k)
    t1, t2 = Kpax3.computelognormconst(cs[k], k, 0.0, data, po, γ, r, priorR)

    if t1 > lumpp
      lumpp = t1
    end
  end
  println("Done.")

  # now that we know the maximum value, we can compute the logarithm of the
  # normalization constant
  z = 0.0

  println("Computing 'z'...")
  for k in 1:size(data, 2)
    println("k = ", k)
    t1, t2 = Kpax3.computelognormconst(cs[k], k, lumpp, data, po, γ, r, priorR)
    z += t2
  end
  println("Done.")

  (log(z), lumpp)
end

function computeProbs(cs,
                      lz::Float64,
                      lumpp::Float64,
                      data::Matrix{UInt8},
                      po::TestPartition,
                      γ::Vector{Float64},
                      r::Float64,
                      priorR::Kpax3.PriorRowPartition)
  (m, n) = size(data)

  P = zeros(Float64, div(n * (n - 1), 2))
  S = zeros(Float64, 3, m)
  K = zeros(Float64, n)

  u = falses(div(n * (n - 1), 2))

  R = zeros(Int, n)

  logprR = 0.0
  logpost = 0.0
  logp = zeros(Float64, m)

  println("Computing probabilities...")
  for k in 1:(n - 1)
    println("k = ", k)

    priorC = Kpax3.AminoAcidPriorCol(data, γ, r)

    st = po.index[po.k .== k][1]
    en = po.index[po.k .== k + 1][1] - 1

    v = zeros(Float64, k)
    n1s = zeros(Float64, k, m)

    for l in st:en
      fill!(R, 0)
      fill!(v, 0)
      fill!(n1s, 0.0)
      fill!(logp, 0.0)

      for a in 1:n
        g = po.partition[a, l]

        R[a] = g
        v[g] += 1

        for b in 1:m
          n1s[g, b] += float(data[b, a])
        end
      end

      logprR = Kpax3.logdPriorRow(R, priorR)

      idx = 1
      for i in 1:(n - 1)
        for j in (i + 1):n
          u[idx] = (R[i] == R[j])
          idx += 1
        end
      end

      for c1 in cs[k], c2 in cs[k], c3 in cs[k], c4 in cs[k]
        logp[1] += Kpax3.computelogp(1, c1, n1s, priorC)
        logp[2] += Kpax3.computelogp(2, c2, n1s, priorC)
        logp[3] += Kpax3.computelogp(3, c3, n1s, priorC)
        logp[4] += Kpax3.computelogp(4, c4, n1s, priorC)

        tmp = exp(logpost - lumpp)

        P[u] += tmp
        K[k] += tmp

        addvalue!(1, c1, tmp, S)
        addvalue!(2, c2, tmp, S)
        addvalue!(3, c3, tmp, S)
        addvalue!(4, c4, tmp, S)
      end
    end
  end

  # no units are in the same cluster
  k = n
  println("k = ", k)

  priorC = Kpax3.AminoAcidPriorCol(data, γ, r)

  v = ones(Float64, k)
  n1s = zeros(Float64, k, m)

  fill!(R, 0)
  fill!(logp, 0.0)

  for a in 1:n
    R[a] = a

    for b in 1:m
      n1s[a, b] = float(data[b, a])
    end
  end

  logprR = Kpax3.logdPriorRow(R, priorR)

  idx = 1
  for i in 1:(n - 1)
    for j in (i + 1):n
      u[idx] = (R[i] == R[j])
      idx += 1
    end
  end

  for c1 in cs[k], c2 in cs[k], c3 in cs[k], c4 in cs[k]
    logp[1] += Kpax3.computelogp(1, c1, n1s, priorC)
    logp[2] += Kpax3.computelogp(2, c2, n1s, priorC)
    logp[3] += Kpax3.computelogp(3, c3, n1s, priorC)
    logp[4] += Kpax3.computelogp(4, c4, n1s, priorC)

    tmp = exp(logpost - lumpp)

    P[u] += tmp
    K[k] += tmp

    addvalue!(1, c1, tmp, S)
    addvalue!(2, c2, tmp, S)
    addvalue!(3, c3, tmp, S)
    addvalue!(4, c4, tmp, S)
  end
  println("Done.")

  (exp(log(P) - lz), exp(log(S) - lz), exp(log(K) - lz))
end

settings = Kpax3.KSettings(ifile, ofile, alpha=α, theta=θ)

x = Kpax3.AminoAcidData(settings)

priorR = Kpax3.EwensPitman(settings.α, settings.θ)
po = TestPartition(size(x.data, 2))
cs = ((UInt8[1], UInt8[2], UInt8[3], UInt8[4]),
      (UInt8[1; 1], UInt8[2; 2], UInt8[3; 3], UInt8[3; 4], UInt8[4; 3],
       UInt8[4; 4]),
      (UInt8[1; 1; 1], UInt8[2; 2; 2], UInt8[3; 3; 3], UInt8[3; 3; 4],
       UInt8[3; 4; 3], UInt8[4; 3; 3], UInt8[3; 4; 4], UInt8[4; 3; 4],
       UInt8[4; 4; 3], UInt8[4; 4; 4]),
      (UInt8[1; 1; 1; 1], UInt8[2; 2; 2; 2], UInt8[3; 3; 3; 3],
       UInt8[3; 3; 3; 4], UInt8[3; 3; 4; 3], UInt8[3; 4; 3; 3],
       UInt8[4; 3; 3; 3], UInt8[3; 3; 4; 4], UInt8[3; 4; 3; 4],
       UInt8[3; 4; 4; 3], UInt8[4; 3; 3; 4], UInt8[4; 3; 4; 3],
       UInt8[4; 4; 3; 3], UInt8[3; 4; 4; 4], UInt8[4; 3; 4; 4],
       UInt8[4; 4; 3; 4], UInt8[4; 4; 4; 3], UInt8[4; 4; 4; 4]),
      (UInt8[1; 1; 1; 1; 1], UInt8[2; 2; 2; 2; 2], UInt8[3; 3; 3; 3; 3],
       UInt8[3; 3; 3; 3; 4], UInt8[3; 3; 3; 4; 3], UInt8[3; 3; 4; 3; 3],
       UInt8[3; 4; 3; 3; 3], UInt8[4; 3; 3; 3; 3], UInt8[3; 3; 3; 4; 4],
       UInt8[3; 3; 4; 3; 4], UInt8[3; 3; 4; 4; 3], UInt8[3; 4; 3; 3; 4],
       UInt8[3; 4; 3; 4; 3], UInt8[3; 4; 4; 3; 3], UInt8[4; 3; 3; 3; 4],
       UInt8[4; 3; 3; 4; 3], UInt8[4; 3; 4; 3; 3], UInt8[4; 4; 3; 3; 3],
       UInt8[3; 3; 4; 4; 4], UInt8[3; 4; 3; 4; 4], UInt8[3; 4; 4; 3; 4],
       UInt8[3; 4; 4; 4; 3], UInt8[4; 3; 3; 4; 4], UInt8[4; 3; 4; 3; 4],
       UInt8[4; 4; 3; 3; 4], UInt8[4; 3; 4; 4; 3], UInt8[4; 4; 3; 4; 3],
       UInt8[4; 4; 4; 3; 3], UInt8[3; 4; 4; 4; 4], UInt8[4; 3; 4; 4; 4],
       UInt8[4; 4; 3; 4; 4], UInt8[4; 4; 4; 3; 4], UInt8[4; 4; 4; 4; 3],
       UInt8[4; 4; 4; 4; 4]),
      (UInt8[1; 1; 1; 1; 1; 1], UInt8[2; 2; 2; 2; 2; 2],
       UInt8[3; 3; 3; 3; 3; 3], UInt8[3; 3; 3; 3; 3; 4],
       UInt8[3; 3; 3; 3; 4; 3], UInt8[3; 3; 3; 4; 3; 3],
       UInt8[3; 3; 4; 3; 3; 3], UInt8[3; 4; 3; 3; 3; 3],
       UInt8[4; 3; 3; 3; 3; 3], UInt8[3; 3; 3; 3; 4; 4],
       UInt8[3; 3; 3; 4; 3; 4], UInt8[3; 3; 3; 4; 4; 3],
       UInt8[3; 3; 4; 3; 3; 4], UInt8[3; 3; 4; 3; 4; 3],
       UInt8[3; 3; 4; 4; 3; 3], UInt8[3; 4; 3; 3; 3; 4],
       UInt8[3; 4; 3; 3; 4; 3], UInt8[3; 4; 3; 4; 3; 3],
       UInt8[3; 4; 4; 3; 3; 3], UInt8[4; 3; 3; 3; 3; 4],
       UInt8[4; 3; 3; 3; 4; 3], UInt8[4; 3; 3; 4; 3; 3],
       UInt8[4; 3; 4; 3; 3; 3], UInt8[4; 4; 3; 3; 3; 3],
       UInt8[3; 3; 3; 4; 4; 4], UInt8[3; 3; 4; 3; 4; 4],
       UInt8[3; 3; 4; 4; 3; 4], UInt8[3; 3; 4; 4; 4; 3],
       UInt8[3; 4; 3; 3; 4; 4], UInt8[3; 4; 3; 4; 3; 4],
       UInt8[3; 4; 4; 3; 3; 4], UInt8[3; 4; 3; 4; 4; 3],
       UInt8[3; 4; 4; 3; 4; 3], UInt8[3; 4; 4; 4; 3; 3],
       UInt8[4; 3; 3; 3; 4; 4], UInt8[4; 3; 3; 4; 3; 4],
       UInt8[4; 3; 4; 3; 3; 4], UInt8[4; 4; 3; 3; 3; 4],
       UInt8[4; 3; 3; 4; 4; 3], UInt8[4; 3; 4; 3; 4; 3],
       UInt8[4; 4; 3; 3; 4; 3], UInt8[4; 3; 4; 4; 3; 3],
       UInt8[4; 4; 3; 4; 3; 3], UInt8[4; 4; 4; 3; 3; 3],
       UInt8[3; 3; 4; 4; 4; 4], UInt8[3; 4; 3; 4; 4; 4],
       UInt8[3; 4; 4; 3; 4; 4], UInt8[3; 4; 4; 4; 3; 4],
       UInt8[3; 4; 4; 4; 4; 3], UInt8[4; 3; 3; 4; 4; 4],
       UInt8[4; 3; 4; 3; 4; 4], UInt8[4; 3; 4; 4; 3; 4],
       UInt8[4; 3; 4; 4; 4; 3], UInt8[4; 4; 3; 3; 4; 4],
       UInt8[4; 4; 3; 4; 3; 4], UInt8[4; 4; 3; 4; 4; 3],
       UInt8[4; 4; 4; 3; 3; 4], UInt8[4; 4; 4; 3; 4; 3],
       UInt8[4; 4; 4; 4; 3; 3], UInt8[3; 4; 4; 4; 4; 4],
       UInt8[4; 3; 4; 4; 4; 4], UInt8[4; 4; 3; 4; 4; 4],
       UInt8[4; 4; 4; 3; 4; 4], UInt8[4; 4; 4; 4; 3; 4],
       UInt8[4; 4; 4; 4; 4; 3], UInt8[4; 4; 4; 4; 4; 4]))

# lmpp = lumpp - lc = lumpp - lumpp - lz = - lz
# lc + lmpp = lumpp + lz - lz = lumpp
lz, lumpp = lognormconst(cs, x.data, po, settings.γ, settings.r, priorR)
probs = computeProbs(cs, lz, lumpp, x.data, po, settings.γ, settings.r, priorR)

@printf("%.100f\n", lz)
@printf("%.100f\n", lumpp)
@printf("%.100f\n", lumpp + lz)

for i in 1:6
  @printf("%.100f\n", probs[3][i])
end

for i in 1:15
  @printf("%.100f\n", probs[1][i])
end

for i in 1:12
  @printf("%.100f\n", probs[2][i])
end

# EwensPitmanPAUT
# lz = 2.82383600433172521348978989408351480960845947265625
# lumpp = -20.077619086471660381221226998604834079742431640625
# lc = lumpp + lz = -17.253783082139936055909856804646551609039306640625

# EwensPitmanPAZT
# lz = 3.718733440312136817595956017612479627132415771484375
# lumpp = -20.701675744062715267546082031913101673126220703125
# lc = lumpp + lz = -16.982942303750579782217755564488470554351806640625

# EwensPitmanZAPT
# lz = 4.09068229126453619670655825757421553134918212890625
# lumpp = -21.0913924952027400649967603385448455810546875
# lc = lumpp + lz = -17.00071020393820475646862178109586238861083984375

# EwensPitmanNAPT
# lz = 5.15291171975398309967886234517209231853485107421875
# lumpp = -21.99795555234826593959951424039900302886962890625
# lc = lumpp + lz = -16.845043832594281951742232195101678371429443359375
=#

function test_mcmc_algorithm()
  ifile = "data/mcmc_6.fasta"
  ofile = "../build/mcmc_6"

  partition = "data/mcmc_6.csv"

  x = Kpax3.AminoAcidData(Kpax3.KSettings(ifile, ofile))

  n = 6

  # EwensPitmanZAPT
  α = 0.0
  θ = 1.0

  settings = Kpax3.KSettings(ifile, ofile, maxclust=1, maxunit=1, alpha=α, theta=θ, op=[0.6; 0.2; 0.2])

  trueProbK = [0.0572011678382854799052026351091626565903425216674804687500000;
               0.3147244396041969372035396190767642110586166381835937500000000;
               0.3980266432259659814540952993411337956786155700683593750000000;
               0.1895011450641402861450046657409984618425369262695312500000000;
               0.0378534754499949971373595758450392168015241622924804687500000;
               0.0026931288174161633307279739568684817641042172908782958984375]

  tmp = [0.5647226958512603367523752240231260657310485839843750000;
         0.3859601712857648747601047034549992531538009643554687500;
         0.3859601712857648747601047034549992531538009643554687500;
         0.3877266490089570361021742428420111536979675292968750000;
         0.2495816768360757109679326504192431457340717315673828125;
         0.3859601712857648747601047034549992531538009643554687500;
         0.3859601712857648747601047034549992531538009643554687500;
         0.3877266490089570361021742428420111536979675292968750000;
         0.2495816768360757109679326504192431457340717315673828125;
         0.5647226958512601147077702989918179810047149658203125000;
         0.2495816768360775983470745131853618659079074859619140625;
         0.3877266490089560369014520802011247724294662475585937500;
         0.2495816768360764881240498880288214422762393951416015625;
         0.3877266490089572026356279366154922172427177429199218750;
         0.3775410735825202035442771375528536736965179443359375000]

  trueProbP = ones(Float64, n, n)
  idx = 0
  for j in 1:(n - 1), i in (j + 1):n
    idx += 1
    trueProbP[i, j] = trueProbP[j, i] = tmp[idx]
  end

  trueProbS = permutedims(
    reshape([0.6466162109955354564405638484458904713392257690429688;
             0.6466162109955360115520761610241606831550598144531250;
             0.6854223966792249989055108017055317759513854980468750;
             0.6810176138774772791606437749578617513179779052734375;
             0.3217860669917309568432983724051155149936676025390625;
             0.3217860669917303462206348285690182819962501525878906;
             0.2886555906155557904568809135525953024625778198242188;
             0.2863666389656327471158192565781064331531524658203125;
             0.0315977220123564370157787095649837283417582511901856;
             0.0315977220123591293066134255695942556485533714294434;
             0.0259220127048742574049633446975349215790629386901856;
             0.0326157471566144094299311007034702925011515617370606],
            (4, 3)))

  Kpax3.kpax3mcmc(x, partition, settings)

  (k, estimK) = Kpax3.readposteriork(settings.ofile)
  (id, estimP) = Kpax3.readposteriorP(settings.ofile)
  (site, aa, freq, estimS) = Kpax3.readposteriorC(settings.ofile)

  @test isapprox(estimK, trueProbK, rtol=0.02)
  @test isapprox(estimP, trueProbP, rtol=0.02)
  @test isapprox(estimS, trueProbS, rtol=0.02)

  # EwensPitmanPAUT
  α = 0.5
  θ = -0.25

  settings = Kpax3.KSettings(ifile, ofile, maxclust=1, maxunit=1, alpha=α, theta=θ, op=[0.6; 0.2; 0.2])

  trueProbK = [0.20304357310850917883726651780307292938232421875000000000;
               0.21850248716493331224697271863988135010004043579101562500;
               0.23654721412654236556427633786370279267430305480957031250;
               0.19408622743422676570901330705964937806129455566406250000;
               0.11197185302918269411698304338642628863453865051269531250;
               0.03584864513660558638097342054606997407972812652587890625]

  tmp = [0.570577011488977547948309165803948417305946350097656250;
         0.458274033588720208776123854477191343903541564941406250;
         0.458274033588720208776123854477191343903541564941406250;
         0.441580861997908746818808367606834508478641510009765625;
         0.356836131052323146661819919245317578315734863281250000;
         0.458274033588720208776123854477191343903541564941406250;
         0.458274033588720208776123854477191343903541564941406250;
         0.441580861997908746818808367606834508478641510009765625
         0.356836131052323146661819919245317578315734863281250000;
         0.570577011488982877018827366555342450737953186035156250;
         0.356836131052323146661819919245317578315734863281250000;
         0.441580861997908746818808367606834508478641510009765625;
         0.356836131052323146661819919245317578315734863281250000;
         0.441580861997908746818808367606834508478641510009765625;
         0.425947160536098878846900106509565375745296478271484375]

  trueProbP = ones(Float64, n, n)
  idx = 0
  for j in 1:(n - 1), i in (j + 1):n
    idx += 1
    trueProbP[i, j] = trueProbP[j, i] = tmp[idx]
  end

  trueProbS = permutedims(
    reshape([0.6628554623025003644798403001914266496896743774414063;
             0.6628554623024998093683279876131564378738403320312500;
             0.6958954347146397712009502356522716581821441650390625;
             0.6907786000474358534262364628375507891178131103515625;
             0.3161542446605755674049476056097773835062980651855469;
             0.3161542446605746792265279054845450446009635925292969;
             0.2858092016103817578631662854604655876755714416503907;
             0.2834957069790527972585891802737023681402206420898438;
             0.0209902930365564316383952814248914364725351333618164;
             0.0209902930365640644216895793761068489402532577514648;
             0.0182953636746360850939829845174244837835431098937988;
             0.0257256929731574761344159441023293766193091869354248],
            (4, 3)))

  Kpax3.kpax3mcmc(x, partition, settings)

  (k, estimK) = Kpax3.readposteriork(settings.ofile)
  (id, estimP) = Kpax3.readposteriorP(settings.ofile)
  (site, aa, freq, estimS) = Kpax3.readposteriorC(settings.ofile)

  @test isapprox(estimK, trueProbK, rtol=0.02)
  @test isapprox(estimP, trueProbP, rtol=0.02)
  @test isapprox(estimS, trueProbS, rtol=0.02)

  # EwensPitmanPAZT
  α = 0.5
  θ = 0.0

  settings = Kpax3.KSettings(ifile, ofile, maxclust=1, maxunit=1, alpha=α, theta=θ, op=[0.6; 0.2; 0.2])

  trueProbK = [0.08297365650265563219445397180606960318982601165771484375;
               0.17858186829001693185503540917125064879655838012695312500;
               0.25777307456537906782401137206761632114648818969726562500;
               0.25380237308441272459091919699858408421277999877929687500;
               0.16734077362357668850023628692724741995334625244140625000;
               0.05952825393395902442428280210151569917798042297363281250]

  tmp = [0.4594673179563065201769234136008890345692634582519531250;
         0.3358129241361126893217203814856475219130516052246093750;
         0.3358129241361126893217203814856475219130516052246093750;
         0.3205116334957199830668628237617667764425277709960937500;
         0.2303208685969153901584860477669280953705310821533203125;
         0.3358129241361126893217203814856475219130516052246093750;
         0.3358129241361126893217203814856475219130516052246093750;
         0.3205116334957199830668628237617667764425277709960937500;
         0.2303208685969170277374473698728252202272415161132812500;
         0.4594673179563261156133080476138275116682052612304687500;
         0.2303208685969153068917592008801875635981559753417968750;
         0.3205116334957179846654184984799940139055252075195312500;
         0.2303208685969174440710816043065278790891170501708984375;
         0.3205116334957179846654184984799940139055252075195312500;
         0.3050339698727991821769478519854601472616195678710937500]

  trueProbP = ones(Float64, n, n)
  idx = 0
  for j in 1:(n - 1), i in (j + 1):n
    idx += 1
    trueProbP[i, j] = trueProbP[j, i] = tmp[idx]
  end

  trueProbS = permutedims(
    reshape([0.6487705348032424268467366346158087253570556640625000;
             0.6487705348032450913819957349915057420730590820312500;
             0.6895045716381241751591346655914094299077987670898438;
             0.6820507372914499377714037109399214386940002441406250;
             0.3276120832340845279162522274418734014034271240234375;
             0.3276120832340849720054620775044895708560943603515625;
             0.2910537694838980504741243748867418617010116577148438;
             0.2876286629947675255181138709303922951221466064453125;
             0.0236173819628183144503363877220181166194379329681397;
             0.0236173819628183456753589553045458160340785980224609;
             0.0194416588781203998614888206475370679982006549835205;
             0.0303205997137780264294448784312407951802015304565430],
            (4, 3)))

  Kpax3.kpax3mcmc(x, partition, settings)

  (k, estimK) = Kpax3.readposteriork(settings.ofile)
  (id, estimP) = Kpax3.readposteriorP(settings.ofile)
  (site, aa, freq, estimS) = Kpax3.readposteriorC(settings.ofile)

  @test isapprox(estimK, trueProbK, rtol=0.02)
  @test isapprox(estimP, trueProbP, rtol=0.02)
  @test isapprox(estimS, trueProbS, rtol=0.02)

  # EwensPitmanNAPT
  α = -1
  θ = 4

  settings = Kpax3.KSettings(ifile, ofile, maxclust=1, maxunit=1, alpha=α, theta=θ, op=[0.6; 0.2; 0.2])

  trueProbK = [0.013987232861628222380101504995764116756618022918701171875;
               0.277957780373217555602849415663513354957103729248046875000;
               0.542335567963948439462740225280867889523506164550781250000;
               0.165719418801206025415595490812847856432199478149414062500]

  tmp = [0.5514954001719989395979837354389019310474395751953125000;
         0.3318939556633083709513698522641789168119430541992187500;
         0.3318939556633083709513698522641789168119430541992187500;
         0.3588148423035991685381418392353225499391555786132812500;
         0.1920174063238060946368079839885467663407325744628906250;
         0.3318939556633083709513698522641789168119430541992187500;
         0.3318939556633083709513698522641789168119430541992187500;
         0.3588148423035966150251852013752795755863189697265625000;
         0.1920174063238067885261983747113845311105251312255859375;
         0.5514954001719960530181197100318968296051025390625000000;
         0.1920174063238073713932863029185682535171508789062500000;
         0.3588148423035962819582778138283174484968185424804687500;
         0.1920174063238073713932863029185682535171508789062500000;
         0.3588148423035962819582778138283174484968185424804687500;
         0.3581303270326594012651355569687439128756523132324218750]

  trueProbP = ones(Float64, n, n)
  idx = 0
  for j in 1:(n - 1), i in (j + 1):n
    idx += 1
    trueProbP[i, j] = trueProbP[j, i] = tmp[idx]
  end

  trueProbS = permutedims(
    reshape([0.6350275618390684995162587256345432251691818237304688;
             0.6350275618390672782709316379623487591743469238281250;
             0.6821394899702427627730116910242941230535507202148438;
             0.6769128361810061722536602246691472828388214111328125;
             0.3271348717192160138189649387641111388802528381347656;
             0.3271348717192145705290329260606085881590843200683594;
             0.2895518564454515608730389431002549827098846435546875;
             0.2864517226506904656879726189799839630722999572753906;
             0.0378375664417737456179935406908043660223484039306641;
             0.0378375664417758481028464245810027932748198509216309;
             0.0283086535842254523320793424545627203769981861114502;
             0.0366354411681910074882750905089778825640678405761719],
            (4, 3)))

  Kpax3.kpax3mcmc(x, partition, settings)

  (k, estimK) = Kpax3.readposteriork(settings.ofile)
  (id, estimP) = Kpax3.readposteriorP(settings.ofile)
  (site, aa, freq, estimS) = Kpax3.readposteriorC(settings.ofile)

  @test isapprox(estimK, trueProbK, rtol=0.02)
  @test isapprox(estimP, trueProbP, rtol=0.02)
  @test isapprox(estimS, trueProbS, rtol=0.02)

  nothing
end

test_mcmc_algorithm()


text:
global const __REFDIC__ = Dict(
"functional" => 523,
"process" => 178,
"want" => 522,
"Methfessel-Paxton" => 333,
"units)" => 232,
"drho" => 534,
"Pseudo" => 148,
"isym" => 19,
"mixing" => 248,
"(lmaxx)" => 451,
"spent" => 42,
"type:" => 419,
"occupation" => 14,
"et" => 276,
"__QEelROUTINES__" => 50,
"cutoff" => 156,
"General" => 336,
"=" => 3,
"renormalized" => 397,
"relaxed" => 496,
"Dal" => 235,
"Reading" => 395,
"To" => 127,
"accuracy" => 67,
"force" => 27,
"processes" => 394,
"Forces" => 210,
"(-TS)" => 193,
"star" => 417,
"follows:" => 80,
"iter" => 47,
"Any" => 521,
"plain" => 332,
"k(" => 6,
"those" => 495,
"quantum" => 393,
"WALL" => 37,
"Convergence" => 126,
"*" => 392,
")" => 1,
"band" => 307,
"code," => 346,
"atomic" => 98,
"Self-consistent" => 103,
"failed," => 543,
"there" => 470,
"calculation" => 99,
"core" => 134,
"contribution" => 86,
"time" => 30,
"shift" => 82,
"star:" => 416,
"and" => 175,
"__UNITkbar__" => 43,
"operations" => 469,
"AE" => 249,
"radius" => 254,
"stress" => 132,
"inside" => 415,
"transformation" => 414,
"(criteria:" => 494,
"Final" => 493,
"MB" => 52,
"dvanqq" => 533,
"cor," => 324,
"bravais-lattice" => 231,
"curvature" => 542,
"End" => 97,
"pseudopotential" => 299,
"width" => 230,
"Quantum" => 391,
"routines" => 250,
"__KWUPF__" => 345,
"smearing" => 192,
"cite" => 390,
"volume" => 153,
"bfgs:" => 436,
"threshold" => 182,
"angular" => 450,
"a(" => 114,
"augmentation" => 323,
"c" => 485,
"__POINTGROUP__" => 468,
"dynamical" => 120,
"axes" => 161,
"mode" => 69,
"potential" => 174,
"site" => 160,
")=" => 75,
"__GRPSYMBOL__" => 90,
"avg" => 59,
"or" => 266,
"deg" => 23,
"numbers" => 13,
"division:" => 275,
"coordinates" => 438,
"Giannozzi" => 274,
"addusddens" => 532,
"if" => 476,
"first" => 467,
"enthalpy" => 166,
"be" => 116,
"with" => 44,
"acting" => 209,
"The" => 70,
"cholesky" => 560,
"static" => 298,
"Called" => 155,
"simulation" => 389,
"negative" => 133,
"(ntypx)" => 449,
"MPI" => 388,
"at" => 183,
"Ultrasoft" => 269,
"program" => 303,
"#:" => 89,
");" => 273,
"mass" => 229,
"irreducible" => 413,
"displacements:" => 412,
"s" => 228,
"renormalised" => 196,
"q" => 123,
"kpt" => 88,
"Representation" => 87,
"Computing" => 83,
"iterations" => 49,
"the" => 71,
"RAM" => 121,
"matrix" => 281,
"preceding" => 492,
"G" => 264,
"__HHMMSS__" => 284,
"secs" => 29,
"XC" => 520,
"procs)" => 387,
"__UNITSTRESS__" => 208,
"extrapolated" => 277,
"read" => 147,
"!" => 191,
"__KWPW__" => 237,
"problem:" => 386,
"data" => 236,
"Error" => 559,
"cdiaghg" => 558,
"Total" => 131,
"INITIALIZATION:" => 526,
"SLA" => 477,
"Sum" => 263,
"PAW" => 110,
"standard" => 448,
"Max" => 163,
"achieved" => 107,
"npool" => 385,
"table:" => 466,
"NUM" => 519,
"analytical," => 481,
"__KWCPU__" => 36,
"GB" => 327,
"structure" => 457,
"(alat" => 227,
"PBX" => 502,
"__UNITVOLb__" => 245,
"-" => 17,
"eigenvalues" => 311,
"stopping" => 435,
"Norm-conserving" => 546,
"__DATEa__" => 283,
"__KWfhiupfx__" => 453,
"uphill" => 554,
"superposition" => 340,
"__au__" => 226,
"__UNITTWOPIALAT__" => 146,
"h,s,v(r/c" => 297,
"__KWSCF__" => 130,
"separable" => 480,
"sum" => 169,
"wfc" => 553,
"__MILLER__" => 22,
"Generated" => 145,
"space" => 384,
"different" => 447,
"Dynamical" => 262,
"Starting" => 301,
"addusdynmat" => 531,
"by" => 102,
"axis)" => 207,
"free" => 339,
"ethr" => 58,
"Atomic" => 411,
"have" => 566,
"Diagonalizing" => 246,
"__UNITDENSITY__" => 244,
"electrons" => 296,
"Alpha" => 410,
"scf" => 54,
"parameter" => 225,
"Writing" => 409,
"terminated" => 424,
"you" => 518,
"Ewald" => 408,
"large:" => 318,
"__PRESSUREEQS__" => 206,
"MATRIX:" => 525,
"group" => 383,
"rho" => 181,
"definition" => 517,
"autoval" => 91,
"run" => 423,
"buffer" => 295,
"PWs)" => 12,
"sticks:" => 261,
"coefficients" => 234,
"processors" => 382,
"steps" => 243,
"distributed-memory" => 381,
">" => 125,
"not" => 539,
"Program" => 380,
"(cartesian" => 205,
"final" => 312,
"__QEgenROUTINES__" => 122,
"for" => 85,
"scalapack" => 379,
"please" => 378,
"character" => 465,
"symmetry" => 464,
"previous" => 434,
"arising" => 377,
"axes:" => 144,
"Band" => 456,
"BFGS" => 431,
"up" => 41,
"given" => 516,
"(MPI)," => 376,
"index" => 224,
"terms:" => 190,
"units" => 113,
"reset" => 433,
"(size" => 375,
"presentations" => 374,
"URL" => 373,
"Crystallographic" => 294,
"trust" => 253,
"alat)" => 223,
"are" => 72,
"ecut=" => 62,
"now" => 40,
"k-points" => 446,
"types" => 222,
"done" => 124,
"converted" => 300,
"one" => 372,
"code" => 186,
"with:" => 143,
"init/wfcrot:" => 293,
"__KWPHONON__" => 515,
"nstep" => 486,
"classes" => 463,
"this" => 329,
"proc/nbgrp/npool/nimage" => 371,
"reason" => 565,
"Fermi" => 73,
"Results" => 491,
"bands" => 11,
"MaX" => 483,
"eigenvalue" => 370,
"__VERSIONa__" => 162,
"*ecutwfc" => 564,
"distributed" => 369,
"celldm(" => 81,
"Ultrasoft," => 498,
"coefficients," => 344,
"resetting" => 398,
"(crystal)" => 242,
"Dense" => 292,
"Corso" => 233,
"density" => 118,
"Message" => 428,
"file" => 239,
"cell:" => 407,
"JOB" => 422,
"on:" => 421,
"calls)" => 39,
"inversion," => 291,
"dense" => 167,
"iterative" => 368,
"PBE" => 325,
"Projector" => 322,
"deallocated" => 418,
"was" => 247,
"Smooth" => 527,
"Waiting" => 445,
"open-source" => 367,
"cpu" => 28,
"Theta=" => 541,
"G-vecs:" => 260,
"will" => 328,
"Zval" => 142,
"modes" => 101,
"beta=" => 61,
"SLA-PZ-NOGX-NOGC" => 549,
"dimensions:" => 265,
"crystal" => 221,
"total" => 18,
"Calculation" => 95,
"publications" => 366,
"Estimated" => 149,
"__DURATION__" => 25,
"#" => 16,
"species" => 238,
"newdq" => 530,
"__QEforceKW__" => 84,
"really" => 514,
"__UNITFORCEb__" => 204,
"old" => 305,
"+" => 158,
"identity" => 164,
"as" => 79,
"to" => 35,
"Parallelization" => 259,
"augmented-wave" => 321,
"file:" => 141,
"output" => 406,
"Optimization" => 430,
"K-points" => 365,
"(ethr)" => 317,
"using" => 173,
"atom" => 26,
"q-points" => 484,
"step:" => 426,
"name" => 462,
"format" => 343,
"wk" => 5,
"__QEDynRAMfor__" => 63,
"solution" => 364,
"in" => 76,
"recalculated" => 337,
"overlap" => 57,
"There" => 405,
"(" => 2,
"addusdbec" => 529,
"-q+G" => 280,
"verify" => 513,
"__BOHR__" => 425,
"PBESOL" => 503,
"dimensions" => 444,
"max" => 290,
"Check:" => 268,
"condition" => 540,
"DFT" => 512,
"starts" => 363,
"used" => 170,
"form" => 479,
"unit-cell" => 152,
"smearing," => 220,
"__REPSYMBOL__" => 128,
"estimated" => 66,
"term" => 203,
"NOGX" => 538,
"zero" => 552,
"matrices" => 511,
"sub-group:" => 362,
"materials;" => 361,
"functions" => 140,
"Norm-conserving," => 309,
"," => 51,
"(Cartesian" => 202,
"smooth" => 165,
"Shape" => 320,
"discarded" => 510,
"List" => 404,
"version" => 360,
"Exchange-correlation" => 197,
"states=" => 289,
"cut-off" => 279,
"bfgs" => 172,
"reciprocal" => 219,
"atoms/cell" => 218,
"convergence" => 154,
"__KWPBC__" => 501,
"__RELPATH__" => 347,
"cycles" => 241,
"some" => 475,
"kinetic-energy" => 217,
"PSX" => 544,
"Subspace" => 359,
"found" => 288,
"nodes" => 278,
"grid" => 139,
"__FULLPATH__" => 119,
"mp" => 310,
"__UNITCMINV__" => 94,
"freq" => 93,
"Harris-Foulkes" => 65,
"no" => 563,
"drhodvus" => 499,
"atoms" => 171,
"lowered" => 316,
"charge=" => 267,
"SLA-PW-PBX-PBC" => 550,
"diagonalization" => 53,
"correction" => 129,
"__UNITVOLc__" => 240,
"momentum" => 443,
"charge" => 115,
"PZ" => 471,
"may" => 490,
"Hamann" => 482,
"positions" => 159,
"already" => 432,
"suite" => 358,
"Parallel" => 282,
"thresh=" => 46,
"problems" => 557,
"(npk)" => 442,
"__CHEM__" => 38,
"pseudized" => 177,
"Atoms" => 403,
"A" => 489,
"Current" => 441,
"all-electron" => 420,
"pressure" => 201,
"algorithm" => 357,
"b(" => 112,
"points=" => 216,
"Irreps" => 78,
")," => 4,
"Using" => 138,
"Cartesian" => 215,
"routine" => 427,
"->" => 109,
"interpolate" => 545,
"randomized" => 348,
"Checking" => 474,
"R" => 356,
"Q(r)" => 176,
"svn" => 472,
"from" => 96,
"computing" => 556,
"type" => 32,
"setup:" => 562,
"IMPORTANT:" => 509,
"bohr" => 252,
"info" => 258,
"(alat=" => 306,
"following" => 189,
"__Ry__" => 15,
"paw" => 335,
"history" => 302,
"<" => 55,
"sub-group" => 355,
"lattice" => 214,
"WARNING:" => 536,
"adddvscf" => 528,
"wavefunction(s)" => 396,
"xc" => 151,
"G-vectors" => 184,
"input" => 272,
"down" => 180,
"enforced" => 508,
"SLA-PW" => 547,
"hartree" => 150,
"of" => 31,
"running" => 354,
"More" => 353,
"Please," => 507,
"__SYMBOLtypeA__" => 24,
"cell" => 437,
"element:" => 461,
"grid:" => 195,
"(up," => 179,
"rotation" => 21,
"xq(" => 454,
"pseudopotentials" => 440,
"Davidson" => 56,
"is" => 33,
"one-electron" => 188,
"estimate" => 64,
"Geometry" => 429,
"sends" => 402,
"PSQ" => 330,
"has" => 104,
":" => 7,
"__EV__" => 9,
"k" => 10,
"one-center" => 334,
"norm" => 551,
"wfcs" => 168,
"[THz]" => 92,
"(with" => 401,
"random" => 331,
"This" => 211,
"__KWFHIPP__" => 452,
"ESPRESSO" => 352,
"on" => 111,
"are:" => 439,
"__QESTOPPING__" => 555,
"s(" => 8,
"representations" => 400,
"tau(" => 77,
"FFT" => 198,
"Min" => 257,
"Structure" => 455,
"axes," => 200,
"__UNITVOLa__" => 213,
"dynmatcc" => 535,
"each" => 460,
"PW" => 256,
"starting" => 341,
"radial" => 137,
"DYNAMICAL" => 524,
"been" => 106,
"directory:" => 506,
"Kohn-Sham" => 287,
"__THREETUPLES__" => 194,
"self-consistent" => 105,
"forces" => 199,
"axis" => 20,
"converged" => 478,
"__UNITFORCEa__" => 488,
"rinner" => 342,
"per" => 157,
"ewald" => 187,
"BESSEL" => 500,
"CASE:" => 304,
"NEW-OLD" => 251,
"__ALAT__" => 212,
"__CHKSUM__" => 136,
"points," => 135,
"new" => 100,
"can" => 473,
"Number" => 399,
"valence" => 286,
"|" => 45,
"Initial" => 338,
"number" => 74,
"small" => 458,
"energy" => 34,
"what" => 505,
"LDA" => 548,
"NOGC" => 537,
"__QEphROUTINES__" => 108,
"&" => 351,
"iteration" => 60,
"part" => 350,
"details" => 349,
"correction," => 255,
"__QEstressKW__" => 48,
"class" => 459,
"__KWPWSCF__" => 313,
"too" => 315,
"__ATOMORBIT__" => 308,
"Title:" => 285,
"beta" => 117,
"ecutrho>" => 561,
"further" => 504,
"charge:" => 319,
"Threshold" => 314,
"inversion" => 185,
"unit" => 326,
"differ" => 487,
"__URL__" => 271,
"l(" => 68,
"Matter" => 270,
"Begin" => 497
) ;

max_stars_repo_name: algorithmx/SmartParser

max_issues_repo_name: algorithmx/SmartParser

max_forks_repo_name: algorithmx/SmartParser


text:
module CorrelationSpectroscopy

using ToeplitzMatrices: Toeplitz

export NodaHilbert, N
export synchronous, asynchronous, disrelation

# TODO: support 2D correlation spectra for distinct spectra.
include("nodahilbert.jl")
include("synchronous.jl")
include("asynchronous.jl")
include("disrelation.jl")

end


text:
# Wrap netcdf.h header file using Clang.jl
# This produces two files, common.jl and netcdf.jl
#
# common.jl is filtered manually, leaving out constants not starting with the NC_ prefix
# or otherwise unneccesary or invalid lines
# the manually selected constant and wrapped types are then copied into src/netcdf_c.jl
#
# from netcdf.jl first ::Cint and ::nc_type in the type signature are loosened to ::Integer
# then Cstring is replaced with Ptr{UInt8}
# then it is copied entirely to the end of src/netcdf_c.jl
# then all ccalls that return an error code are wrapped in the `check` function

using Clang

workdir = dirname(@__FILE__)
srcdir = joinpath(workdir, "../src")

headerpath = ["/home/martijn/src/netcdf-4.4.1/include/netcdf.h"]
includedirs = ["/home/martijn/src/netcdf-4.4.1/include"]

context = wrap_c.init(
    headers = headerpath,
    output_dir = workdir,
    common_file = "common.jl",
    clang_includes = includedirs,
    header_library = x -> "libnetcdf",
)

run(context)


text:
# pour le nnls voir optim.jl + ajouter les gros poids pour full et empty comme dans R !!

function cost_function(
    forest              :: EnsembleCausal{S},
    set_u               :: Set{Int},
    base_causal_effect  :: Vector{Float64}
    ) where {S}

    n = length(forest.Y)
    som = 0
    projected_causal_effect = prf_causal_effect(forest, set_u)

    for i in 1:n
        som += (base_causal_effect[i] - projected_causal_effect[i])^2
    end

    som /= n * var(base_causal_effect)
    return 1 - som

end

function weight(set_u :: Set{Int}, n_imp :: Int)
    card_u = length(set_u)
    combi = binomial(n_imp, card_u)
    res = (n_imp - 1) / (combi * card_u * (n_imp - card_u))
    if res > 1000     #TODO histoire des poids infinis
        return 1000
    else
        return res
    end
end

function to_binary_vector(set_u :: Set{Int}, n_imp :: Int)
    res = [0 for i=1:n_imp]
    for j in set_u
        res[j] = 1
    end
    return res
end

function to_binary_vector_comp(set_u :: Set{Int}, n_imp :: Int)
    res = [1 for i=1:n_imp]
    for j in set_u
        res[j] = 0
    end
    return res
end

#TODO faire estimate Shapley en regardant bien le code R de shaff + attention besoin des probas de U !!

function estimate_Shapley(
    forest    :: EnsembleCausal{S},
    dict_sets :: Dict{Set{Int}, Int},
    K         :: Int,
    proba_set :: Dict{Set{Int}, Float64}
    ) where {S}

    n_imp = size(forest.X, 2)

    w_vect = cat([sqrt((weight(U, n_imp)*dict_sets[U])/(K*proba_set[U])) for U in keys(dict_sets)],
                 [sqrt((weight(U, n_imp)*dict_sets[U])/(K*proba_set[U])) for U in keys(dict_sets)],
                 10000, 10000; dims = 1)
    w_matrix = diagm(w_vect)

    z_vect = cat([to_binary_vector(U, n_imp) for U in keys(dict_sets)],
                 [to_binary_vector_comp(U, n_imp) for U in keys(dict_sets)],
                 [[0 for i=1:n_imp]], [[1 for i=1:n_imp]]; dims = 1)
    z_matrix = vec(z_vect)

    cause = apply_forest_oob(forest)

    cost_matrix = cat([transpose(cost_function(forest, set_u, cause)) for set_u in keys(dict_sets)],
                      [transpose(cost_function(forest, setdiff(Set(1:n_imp), set_u), cause)) for set_u in keys(dict_sets)],
                      0, 1; dims = 1)

    A_matrix = w_matrix * z_matrix
    A_matrix = reduce(hcat, A_matrix)'

    b_matrix = w_matrix * cost_matrix

    beta = nonneg_lsq(A_matrix, b_matrix)

    return beta

end

function shaff(forest :: EnsembleCausal{S}, K :: Int) where {S}
    freq = get_occurence_frequencies(forest)
    samp = sample_U(freq, K)
    dict_samp = StatsBase.countmap(samp)
    return estimate_Shapley(forest, dict_samp, K, freq)
end


text:
"""
Functions to calculate different line broadenings.
"""

# For now
const mass_H = 1.008 * m_u
const mass_He = 4.003 * m_u
const abund_He = 10^10.99 / 10^12  # From RH

#=----------------------------------------------------------------------------
                    Linear Stark broadening: C_2 / r^2
----------------------------------------------------------------------------=#
"""
    γ_linear_stark(
        electron_density::NumberDensity{T},
        n_upper::Integer,
        n_lower::Integer
    ) where T <: AbstractFloat

Calculate linear Stark broadening according to the approximation of
[Sutton (1978)](https://ui.adsabs.harvard.edu/abs/1978JQSRT..20..333S/), eq (24),
so that it can be added into a Voigt profile and avoid the computation of a Holtsmark
profile. Valid up to electron densities of 1e19 m^-3 in the chromosphere.

# Arguments
- `electron_density`: electron density per volume
- `n_upper`: principal quantum number of upper level
- `n_lower`: principal quantum number of lower level
"""
function γ_linear_stark(
    electron_density::NumberDensity{T},
    n_upper::Integer,
    n_lower::Integer
) where T <: AbstractFloat
    @assert n_upper > n_lower
    @assert n_upper > 1
    if n_upper - n_lower == 1
        a1 = convert(T, 0.642)
    else
        a1 = convert(T, 1)
    end
    γβ = convert(T, 0.425)
    power = convert(T, 2/3)
    zcoeff = convert(T, 6e-5)u"s^-1"  # should be m^2/s, dropping m to allow ne=ustrip(...)
    ne = ustrip(electron_density |> u"m^-3")
    # Factor of 2 to convert from half half-width to half-width
    return 2 * γβ * a1 * zcoeff * (n_upper^2 - n_lower^2) * ne^power
end

#=----------------------------------------------------------------------------
                    Quadratic Stark broadening: C_4 / r^4
----------------------------------------------------------------------------=#
"""
    c4_traving(line::AtomicLine)

Calculate the \$C_4\$ interaction constant (quadratic Stark effect) using the
recipe of Traving (1960), "Uber die Theorie der Druckverbreiterung von Spektrallinien",
p 93.
"""
function c4_traving(line::AtomicLine)
    n_eff_u = n_eff(line.χ∞, line.χj, line.Z)
    n_eff_l = n_eff(line.χ∞, line.χi, line.Z)
    C4 = (e^2 * inv_4πε0 * a_0^3 * 2 * π / (h * 18 * line.Z^4) *
        ((n_eff_u * (5 * n_eff_u^2 + 1))^2 - (n_eff_l * (5 * n_eff_l^2 + 1))^2))
    return C4 |> u"m^4 / s"
end

"""
    const_quadratic_stark(line::AtomicLine;
                          mean_atomic_weight::Unitful.Mass=28 * m_u,
                          scaling::Real=1)

Calculate height-independent constant to use in `γ_quadratic_stark`, using the recipe
from RH, which is based on the following estimate:

\$\$
\\gamma = 11.37 \\cdot vrel^{1/3} * C_4^{2/3} * (n_e + n_{ion}),
\$\$

Using the estimate for \$C_4\$ from Traving (1960), "Uber die Theorie der
Druckverbreiterung von Spektrallinien", p 93., and \$n_{ion}\\approx n_e\$
(following Gray).
"""
function const_quadratic_stark(line::AtomicLine;
                               mean_atomic_weight::Unitful.Mass=28 * m_u,
                               scaling::Real=1)
    C = ustrip(8 * k_B / (π * line.atom_weight) |> u"J/(K * kg)")
    Cm = ((1 + line.atom_weight / m_e)^(1/6) +
          (1 + line.atom_weight / mean_atomic_weight)^(1/6))
    C4 = ustrip(c4_traving(line) |> u"m^4 / s")
    cStark23 = 11.37u"m^3 / s" * (scaling * C4)^(2/3)
    return C^(1/6) * cStark23 * Cm
end

"""
    γ_quadratic_stark(
        electron_density::NumberDensity,
        temperature::Unitful.Temperature;
        stark_constant::Unitful.VolumeFlow=1.0u"m^3 / s",
    )

Compute quadratic Stark broadening for a given `electron_density`.  If `temperature` is
nonzero, then it will apply the standard recipe of RH (using \$C_4\$ from Traving 1960).
The `stark_constant` can be obtained either from atomic data sources, or, if using the RH
recipe, using the function `const_quadratic_stark`.
"""
function γ_quadratic_stark(
    electron_density::NumberDensity,
    temperature::Unitful.Temperature;
    stark_constant::Unitful.VolumeFlow=1.0u"m^3 / s",
)
    if temperature > 0u"K"
        t_factor = ustrip(temperature |> u"K")^(1/6)
    else
        t_factor = 1.0
    end
    return stark_constant * t_factor * electron_density
end

"""
    γ_quadratic_stark_gray(
        electron_density::NumberDensity,
        temperature::Unitful.Temperature,
        c4::Quantity{<:AbstractFloat, Unitful.𝐋^4 / Unitful.𝐓},
    )

Compute quadratic Stark broadening using the recipe of Gray (2005), page 244, eq 11.27.
The interaction constant `c4` should be provided, either from atomic data or from the
estimate of Traving (1960) using `c4_traving`.
"""
function γ_quadratic_stark_gray(
    electron_density::NumberDensity,
    temperature::Unitful.Temperature,
    c4::Quantity{<:AbstractFloat, Unitful.𝐋^4 / Unitful.𝐓},
)
    # Formula assumes CGS
    ne_term = ustrip(electron_density * k_B |> u"erg / (K * cm^3)")
    c4_term = ustrip(c4 |> u"cm^4 / s")
    t_term = ustrip(temperature |> u"K")
    log10γ = 19 + (2/3) * log10(c4_term) + log10(ne_term) + (1/6) * log10(t_term)
    return (10^log10γ) * u"s^-1"
end

#=----------------------------------------------------------------------------
                 van der Waals broadening, C_6 / r^6
----------------------------------------------------------------------------=#
"""
    function const_unsold(line::AtomicLine; H_scaling::Real=1, He_scaling::Real=1)

Compute atmosphere-independent constant for γ_unsold, to be used in function `γ_unsold`.
Based on expressions from RH broad.c, which uses formula in Mihalas (1978),
pp 282, 286-287, eq. (9-50) for v_rel, table 9-1 and eq. (9-76) for the interaction
coefficient C6.  Input is an `AtomicLine`, which contains the necessary energies and
element weight. The van der Waals broadening can be scaled for both H and He perturbers
using `H_scaling` and `He_scaling`.
"""
function const_unsold(line::AtomicLine; H_scaling::Real=1, He_scaling::Real=1)
    Δr = (Ry^2 * (1 / (line.χ∞ - line.χj)^2 - 1 / (line.χ∞ - line.χi)^2)) |> u"J/J"
    C6 = ustrip((2.5 * e^2 * αp * inv_4πε0^2 * 2 * π *
                 (line.Z * a_0)^2 / h * Δr) |> u"C^2 * m^6 / (F * J * s)")
    v_rel_const = ustrip(8 * k_B / (π * line.atom_weight) |> u"J/(K * kg)")
    v_rel_H = v_rel_const * (1 + line.atom_weight / mass_H)
    v_rel_He = v_rel_const * (1 + line.atom_weight / mass_He)
    return 8.08 * (H_scaling * v_rel_H^0.3 + He_scaling * abund_He * v_rel_He^0.3) * C6^0.4
end

"""
    function γ_unsold(unsold_const::AbstractFloat, temperature::Unitful.Temperature,
                      h_neutral_density::NumberDensity)

Compute van der Waals broadening in Lindholm theory using Unsöld's approximation
for the interaction coefficient \$C_6\$. Based on Mihalas (1978), pp 282, 286-287.
Takes the atmosphere-indepenent `unsold_const` from `γ_unsold_const`, temperature,
and populations of neutral hydrogen, and returns broadening in units of s^-1.
"""
function γ_unsold(
    unsold_const::AbstractFloat,
    temperature::Unitful.Temperature,
    h_neutral_density::NumberDensity
)
    return (unsold_const * ustrip(temperature |> u"K")^0.3 *
            ustrip(h_neutral_density |> u"m^-3") * u"s^-1")
end

"""
    function const_barklem(atomic_weight::Unitful.Mass, α::Real, σ::Real)

Compute the atmosphere-independent constant used to calculate broadening from collisions
with neutral hydrogen using the recipes of Barklem/O'Mara/Anstee, in the function
`γ_barklem`. The calculation performed here follows eq (3) of
[Anstee & O'Mara (1995)](https://ui.adsabs.harvard.edu/abs/1995MNRAS.276..859A).

# Arguments
- `atomic_weight::Unitful.Mass`: atomic weight of element
- `α::Real`: velocity exponent from Barklem/O'Mara/Anstee tables
- `σ::Real`: line broadening cross section in atomic units (a_0^2) for a collision
   velocity of 10 km/s, from Barklem/O'Mara/Anstee tables.

# Returns
- `Unitful.VolumeFlow`: line broadening width per neutral hydrogen atom. Needs
   to be multiplied by temperature ^ ((1 - α)/2) to give proper temperature dependence.
"""
function const_barklem(atomic_weight::Unitful.Mass, α::Real, σ::Real)
    α < 0 && error("α must be non-negative")
    σ < 0 && error("σ must be non-negative")
    μ = m_u / (1 / Ar_H + 1 / (atomic_weight / m_u))
    # Using 1 K to keep units right for later multiplication by correct temperature
    v_bar = sqrt(8 * k_B * u"K"/ (π * μ)) |> u"m/s"
    v_ratio = (1e4u"m/s" / v_bar) |> u"m/m"
    # Squared Bohr radius is to convert from atomic units to m^2, factor of 2 from HW to FW
    return (a_0^2 * 2 * (4 / π)^(α / 2) * gamma((4 - α) / 2) * v_bar * σ *
            v_ratio^α) |> u"m^3 / s"
end

"""
    function γ_barklem(
        α::AbstractFloat,
        barklem_const::Unitful.VolumeFlow,
        temperature::Unitful.Temperature,
        h_neutral_density::NumberDensity,
    )

Compute van der Waals broadening from collisions with neutral hydrogen atoms following the
theory from Barklem/O'Mara/Anstee.

# Arguments
- `α::AbstractFloat`: velocity exponent from Barklem/O'Mara/Anstee tables
- `barklem_const::Unitful.VolumeFlow`: atmosphere-independent constant computed from
   `const_barklem()`.
- `temperature::Unitful.Temperature`
- `h_neutral_density::NumberDensity`: number density of neutral hydrogen atoms

# Returns
- `γ::Unitful.Frequency`: broadening in units of s^-1.

# Notes
This computes broadening only from hydrogen atoms. For collisions with helium atoms,
it is recommended to add van der Waals broadening using Unsöld's approximation (see example).

# Examples
Example for Ca II 854.2 nm line:
```
julia> Ca8542 = AtomicLine(25414.400u"cm^-1", 13710.880u"cm^-1", 95785.470u"cm^-1",
4, 6, 7.242e-02, 40.08 * m_u, 20);

julia> temp = 6000u"K";

julia> h_density = 1e23u"m^-3";

julia> bconst = const_barklem(Ca8542.atom_weight, 0.275, 291)
7.495208174533257e-16 m³ s⁻¹

julia> γ = γ_barklem(0.275, bconst, temp, h_density)
1.3596876505340942e11 s⁻¹
```

Now adding van der Waals broadening for helium as well:
```
julia> uconst = const_unsold(Ca8542; H_scaling=0, He_scaling=1)
2.482484115415461e-16

julia> γ = γ_barklem(0.275, bconst, temp, h_density) + γ_unsold(uconst, temp, h_density)
1.3630631018876224e11 s⁻¹
```

# References
- [Anstee & O'Mara (1995)](https://ui.adsabs.harvard.edu/abs/1995MNRAS.276..859A)
- [Barklem & O'Omara (1997)](https://ui.adsabs.harvard.edu/abs/1997MNRAS.290..102B)
- [Barklem, O'Mara & Ross (1998)](https://ui.adsabs.harvard.edu/abs/1998MNRAS.296.1057B)
- [Barklem & O'Mara (1998)](https://ui.adsabs.harvard.edu/abs/1998MNRAS.300..863B)
"""
function γ_barklem(
    α::AbstractFloat,
    barklem_const::Unitful.VolumeFlow,
    temperature::Unitful.Temperature,
    h_neutral_density::NumberDensity,
)
    return barklem_const * ustrip(temperature |> u"K")^((1 - α)/2) * h_neutral_density
end

"""
    function const_deridder_rensbergen(
        atomic_mass_pert::Unitful.Mass,
        atomic_mass_rad::Unitful.Mass
        α::Real,
        β::Real
    )

Compute the atmosphere-independent constant used to calculate broadening using
the recipes of [Deridder & Rensbergen (1976)](https://ui.adsabs.harvard.edu/abs/1976A%26AS...23..147D),
in the function `γ_deridder_rensbergen`.

# Arguments
- `atomic_weight_pert::Unitful.Mass`: atomic mass of perturbing element (hydrogen or helium)
- `atomic_weight_rad::Unitful.Mass`: atomic mass of line-producing species
- `α::Real`: α parameter as taken from the tables of Deridder & Rensbergen (1976),
   in units of 10^-8 cm^3/s (ignoring the dimensions of the temperature exponent)
- `β::Real`: β parameter as taken from the tables of Deridder & Rensbergen (1976),
   dimensionless.

# Returns
- `Unitful.VolumeFlow`: line broadening width per perturber atom. Needs
   to be multiplied by temperature ^ β to give proper temperature dependence.
"""
function const_deridder_rensbergen(
    atomic_mass_pert::Unitful.Mass,
    atomic_mass_rad::Unitful.Mass,
    α::Real,
    β::Real,
)
    α < 0 && error("α must be non-negative")
    α = (α * 1e-8u"cm^3/s") |> u"m^3/s"  # Convert from paper's 1e-8 units to SI
    mass_corr = (1 + atomic_mass_pert / atomic_mass_rad) ^ β
    return α * mass_corr
end

"""
    function γ_deridder_rensbergen(
        β::Real,
        deridder_const::Unitful.VolumeFlow,
        temperature::Unitful.Temperature,
        h_neutral_density::NumberDensity,
    )

Compute van der Waals broadening from collisions with neutral hydrogen or helium
atoms following [Deridder & Rensbergen (1976)](https://ui.adsabs.harvard.edu/abs/1976A%26AS...23..147D).

# Arguments
- `β::Real`: β parameter as taken from the tables of Deridder & Rensbergen (1976),
   dimensionless.
- `deridder_const::Unitful.VolumeFlow`: atmosphere-independent constant computed from
   `const_deridder_rensbergen()`.
- `temperature::Unitful.Temperature`
- `perturber_density::NumberDensity`: number density of perturber atoms, either
   neutral hydrogen or helium.

# Returns
- `γ::Unitful.Frequency`: broadening in units of s^-1.
"""
function γ_deridder_rensbergen(
    β::Real,
    deridder_const::Unitful.VolumeFlow,
    temperature::Unitful.Temperature,
    perturber_density::NumberDensity,
)
    return deridder_const * ustrip(temperature |> u"K")^β * perturber_density
end

#=----------------------------------------------------------------------------
                    Radiation utilities
----------------------------------------------------------------------------=#
"""
    function calc_Aji(λ0::Unitful.Length, g_ratio::Real, f_value::AbstractFloat)

Compute the spontaneous deexcitation rate \$A_{ul}\$ (natural broadening)
for a bound-bound transition, using the SI expression:

\$\$
A_{ul} = \\frac{2\\pi e^2}{\\varepsilon_0 m_e c} \\frac{g_l}{g_u} \\frac{f_{lu}}{\\lambda^2}
\$\$

for a given rest wavelength `λ0`, ration between statistical weights of lower and
upper levels (`g_ratio` = gl / gu), and `f_value` .
"""
function calc_Aji(λ0::Unitful.Length, g_ratio::Real, f_value::AbstractFloat)
    (2π * e^2 / (ε_0 * m_e * c_0) * g_ratio * f_value / λ0^2) |> u"s^-1"
end

"""
    function calc_Bji(λ0::Unitful.Length, Aji::Unitful.Frequency)

Compute the induced deexcitation rate \$B_{ul}\$ for a bound-bound transition,
using the SI expression *per wavelength*:

\$\$
B_{ul} = \\frac{\\lambda^5}{2 h c} A_{ul}
\$\$

for a given rest wavelength `λ0`, and spontaneous deexcitation rate `Aji.`
"""
calc_Bji(λ0::Unitful.Length, Aji::Unitful.Frequency) = (λ0^5 * Aji / (2h * c_0^2)) |> u"m^3 / J"

"""
Compute damping parameter.
"""
function damping(γ::Unitful.Frequency, λ::Unitful.Length, ΔλD::Unitful.Length)
    return (γ * λ^2 / (4 * π * c_0 * ΔλD)) |> u"m/m"
end


text:
macro primitive(ex)
    @capture(shortdef(ex), (logpdf(argname_::name_, args__) = body_) | (logpdf(argname_::name_, args__) where {Ts__} = body_)) || error("PrimitiveError: defining new contextual primitive calls requires a logpdf definition for the call.")
    argname = gensym(argname)
    expr = quote
        $ex

        @inline function (ctx::Jaynes.SimulateContext)(call::typeof(rand), addr::Jaynes.Address, $argname::$name, args...)
            Jaynes.visit!(ctx.visited, addr)
            s = $argname(args...)
            Jaynes.add_choice!(ctx, addr, Jaynes.ChoiceSite(logpdf($argname, args..., s), s))
            return s
        end

        function (ctx::Jaynes.GenerateContext)(call::typeof(rand), addr::T, $argname::$name, args...) where {T <: Jaynes.Address, K}
            Jaynes.visit!(ctx.visited, addr)
            if Jaynes.has_top(ctx.select, addr)
                s = Jaynes.get_top(ctx.select, addr)
                score = logpdf($argname, args..., s)
                Jaynes.add_choice!(tr, addr, Jaynes.ChoiceSite(score, s))
                increment!(ctx, score)
            else
                s = $argname(args...)
                score = logpdf($argname, args..., s)
                Jaynes.add_choice!(ctx.tr, addr, Jaynes.ChoiceSite(score, s))
            end
            return s
        end

        @inline function (ctx::Jaynes.ProposeContext)(call::typeof(rand), addr::T, $argname::$name, args...) where {T <: Jaynes.Address, K}
            s = $argname(args...)
            score = logpdf($argname, args..., s)
            Jaynes.add_choice!(ctx.tr, addr, Jaynes.ChoiceSite(score, s))
            increment!(ctx, score)
            return s
        end

        @inline function (ctx::Jaynes.RegenerateContext)(call::typeof(rand), addr::T, $argname::$name, args...) where {T <: Jaynes.Address, K}
            visit!(ctx.visited, addr)
            in_prev_chm = has_top(ctx.prev, addr)
            in_sel = has_top(ctx.select, addr)
            if in_prev_chm
                prev = get_top(ctx.prev, addr)
                if in_sel
                    ret = $argname(args...)
                    add_choice!(ctx.discard, addr, prev)
                else
                    ret = prev.val
                end
            end
            score = logpdf($argname, args..., ret)
            if in_prev_chm && !in_sel
                increment!(ctx, score - prev.score)
            end
            add_choice!(ctx.tr, addr, ChoiceSite(score, ret))
            return ret
        end

        @inline function (ctx::Jaynes.UpdateContext)(call::typeof(rand),
                                                     addr::T,
                                                     $argname::$name,
                                                     args...) where {T <: Jaynes.Address, K}
            in_prev_chm = Jaynes.has_top(ctx.prev, addr)
            in_prev_chm && begin
                prev = Jaynes.get_top(ctx.prev, addr)
                prev_ret = prev.val
                prev_score = prev.score
            end
            in_selection = Jaynes.has_top(ctx.select, addr)
            if in_selection
                ret = Jaynes.get_top(ctx.select, addr)
                in_prev_chm && begin
                    Jaynes.add_choice!(ctx.discard, addr, prev)
                end
                Jaynes.visit!(ctx.visited, addr)
            elseif in_prev_chm
                ret = prev_ret
            else
                ret = $argname(args...)
            end
            score = logpdf($argname, args..., ret)
            if in_prev_chm
                increment!(ctx, score - prev_score)
            elseif in_selection
                increment!(ctx, score)
            end
            Jaynes.add_choice!(ctx.tr, addr, ChoiceSite(score, ret))
            return ret
        end

        @inline function (ctx::Jaynes.ScoreContext)(call::typeof(rand), addr::T, $argname::$name, args...) where {T <: Jaynes.Address, K}
            Jaynes.has_top(ctx.select, addr) || error("ScoreError: constrained selection must provide constraints for all possible addresses in trace. Missing at address $addr.")
            val = Jaynes.get_top(ctx.select, addr)
            increment!(ctx, logpdf(d, val))
            return val

        end
    end
    expr = MacroTools.prewalk(unblock ∘ rmlines, expr)
    esc(expr)
end

macro primitive(inject, ex)
    @capture(inject, nm_ -> injbody_) || error("InjectionError: defining an injection requires a valid lambda expression.")
    @capture(shortdef(ex), (logpdf(argname_::name_, args__) = body_) | (logpdf(argname_::name_, args__) where {Ts__} = body_)) || error("PrimitiveError: defining new contextual primitive calls requires a logpdf definition for the call.")
    argname = gensym(argname)
    inject_name = gensym("inject")
    expr = quote
        $ex
        $inject_name = $nm -> $injbody
        @inline function (ctx::Jaynes.SimulateContext)(call::typeof(rand), addr::Jaynes.Address, $argname::$name, args...)
            Jaynes.visit!(ctx.visited, addr)
            s = $argname(args...)
            s = $inject_name(s)
            Jaynes.add_choice!(ctx, addr, Jaynes.ChoiceSite(logpdf($argname, args..., s), s))
            return s
        end

        function (ctx::Jaynes.GenerateContext)(call::typeof(rand), addr::T, $argname::$name, args...) where {T <: Jaynes.Address, K}
            Jaynes.visit!(ctx.visited, addr)
            if Jaynes.has_top(ctx.select, addr)
                s = Jaynes.get_top(ctx.select, addr)
                s = $inject_name(s)
                score = logpdf($argname, args..., s)
                Jaynes.add_choice!(tr, addr, Jaynes.ChoiceSite(score, s))
                increment!(ctx, score)
            else
                s = $argname(args...)
                s = $inject_name(s)
                score = logpdf($argname, args..., s)
                Jaynes.add_choice!(ctx.tr, addr, Jaynes.ChoiceSite(score, s))
            end
            return s
        end

        @inline function (ctx::Jaynes.ProposeContext)(call::typeof(rand), addr::T, $argname::$name, args...) where {T <: Jaynes.Address, K}
            s = $argname(args...)
            s = $inject_name(s)
            score = logpdf($argname, args..., s)
            Jaynes.add_choice!(ctx.tr, addr, Jaynes.ChoiceSite(score, s))
            increment!(ctx, score)
            return s
        end

        @inline function (ctx::Jaynes.RegenerateContext)(call::typeof(rand), addr::T, $argname::$name, args...) where {T <: Jaynes.Address, K}
            visit!(ctx.visited, addr)
            in_prev_chm = has_top(ctx.prev, addr)
            in_sel = has_top(ctx.select, addr)
            if in_prev_chm
                prev = get_top(ctx.prev, addr)
                if in_sel
                    ret = $argname(args...)
                    ret = $inject_name(ret)
                    add_choice!(ctx.discard, addr, prev)
                else
                    ret = prev.val
                end
            end
            score = logpdf($argname, args..., ret)
            if in_prev_chm && !in_sel
                increment!(ctx, score - prev.score)
            end
            add_choice!(ctx.tr, addr, ChoiceSite(score, ret))
            return ret
        end

        @inline function (ctx::Jaynes.UpdateContext)(call::typeof(rand),
                                                     addr::T,
                                                     $argname::$name,
                                                     args...) where {T <: Jaynes.Address, K}
            in_prev_chm = Jaynes.has_top(ctx.prev, addr)
            in_prev_chm && begin
                prev = Jaynes.get_top(ctx.prev, addr)
                prev_ret = prev.val
                prev_score = prev.score
            end
            in_selection = Jaynes.has_top(ctx.select, addr)
            if in_selection
                ret = Jaynes.get_top(ctx.select, addr)
                in_prev_chm && begin
                    Jaynes.add_choice!(ctx.discard, addr, prev)
                end
                Jaynes.visit!(ctx.visited, addr)
            elseif in_prev_chm
                ret = prev_ret
            else
                ret = $argname(args...)
                ret = $inject_name(ret)
            end
            score = logpdf($argname, args..., ret)
            if in_prev_chm
                increment!(ctx, score - prev_score)
            elseif in_selection
                increment!(ctx, score)
            end
            Jaynes.add_choice!(ctx.tr, addr, ChoiceSite(score, ret))
            return ret
        end

        @inline function (ctx::Jaynes.ScoreContext)(call::typeof(rand), addr::T, $argname::$name, args...) where {T <: Jaynes.Address, K}
            Jaynes.has_top(ctx.select, addr) || error("ScoreError: constrained selection must provide constraints for all possible addresses in trace. Missing at address $addr.")
            val = Jaynes.get_top(ctx.select, addr)
            increment!(ctx, logpdf(d, val))
            return val

        end
    end
    expr = MacroTools.prewalk(unblock ∘ rmlines, expr)
    esc(expr)
end


text:
include("Commandline.jl")
include("Params.jl")
include("Solution.jl")
include("LocalSearch.jl")

function main()

    c = Commandline()

    println("Input File: \"",c.instance_path,"\"")
    println("Running code with seed=", c.seed)
    # Initialization of the problem data from the commandline
    params = Params(c.instance_path, c.output_path, c.seed, c.cpu_time)

    println("----- STARTING METAHEURISTIC OPTIMIZATION")

    # TODO: Run algorithm
    #p = [4  13  14   7  16  26  25  17   1  15  20  18  12  19   3  8  21   9   5   6  10  24   2  22  11  23]
    p = [i for i = 1:params.datasetSize]
    sol = Solution(params, p)
    printSolution(sol)
    println()
    #------------------------------------------------

    println("Random solution:")
    Rsol = Solution(params)
    printSolution(Rsol)
    println()
    #------------------------------------------------

    # Local Search
    println("Initial Solution:")
    printSolution(sol)
    println()
    #------------------------------------------------

    bestSol = descentHeuristic(sol)

    println("Best solution after local search (descent heuristic):")
    printSolution(bestSol)
    println()
    #------------------------------------------------

    bestSol = multistartLS(params, 1000)
    println("Multistart Local Search solution:")
    printSolution(bestSol)
    println()
    #------------------------------------------------

    #cost = sum(A[i,j]*B[p[i],p[j]] for i = 1:n, j = 1:n)
    #println("Cost=$cost")

    println("----- METAHEURISTIC OPTIMIZATION COMPLETED IN ")

    # TODO: Print final solution and export results
    # Printing the solution and exporting statistics (also export results into a file)

    println("----- END OF ALGORITHM")
end

main()

max_stars_repo_name: BSAraujo/QAP-Metaheuristics

max_issues_repo_name: BSAraujo/QAP-Metaheuristics

max_forks_repo_name: BSAraujo/QAP-Metaheuristics


text:
@testset "36.valid-sudoku.jl" begin
    board = [
        ['5', '3', '.', '.', '7', '.', '.', '.', '.'],
        ['6', '.', '.', '1', '9', '5', '.', '.', '.'],
        ['.', '9', '8', '.', '.', '.', '.', '6', '.'],
        ['8', '.', '.', '.', '6', '.', '.', '.', '3'],
        ['4', '.', '.', '8', '.', '3', '.', '.', '1'],
        ['7', '.', '.', '.', '2', '.', '.', '.', '6'],
        ['.', '6', '.', '.', '.', '.', '2', '8', '.'],
        ['.', '.', '.', '4', '1', '9', '.', '.', '5'],
        ['.', '.', '.', '.', '8', '.', '.', '7', '9'],
    ]
    @test is_valid_sudoku(board) == true

    board = [
        ['.', '.', '4', '.', '.', '.', '6', '3', '.'],
        ['.', '.', '.', '.', '.', '.', '.', '.', '.'],
        ['5', '.', '.', '.', '.', '.', '.', '9', '.'],
        ['.', '.', '.', '5', '6', '.', '.', '.', '.'],
        ['4', '.', '3', '.', '.', '.', '.', '.', '1'],
        ['.', '.', '.', '7', '.', '.', '.', '.', '.'],
        ['.', '.', '.', '5', '.', '.', '.', '.', '.'],
        ['.', '.', '.', '.', '.', '.', '.', '.', '.'],
        ['.', '.', '.', '.', '.', '.', '.', '.', '.'],
    ]
    @test is_valid_sudoku(board) == false

    board = [
        ['8', '3', '.', '.', '7', '.', '.', '.', '.'],
        ['6', '.', '.', '1', '9', '5', '.', '.', '.'],
        ['.', '9', '8', '.', '.', '.', '.', '6', '.'],
        ['8', '.', '.', '.', '6', '.', '.', '.', '3'],
        ['4', '.', '.', '8', '.', '3', '.', '.', '1'],
        ['7', '.', '.', '.', '2', '.', '.', '.', '6'],
        ['.', '6', '.', '.', '.', '.', '2', '8', '.'],
        ['.', '.', '.', '4', '1', '9', '.', '.', '5'],
        ['.', '.', '.', '.', '8', '.', '.', '7', '9'],
    ]
    @test is_valid_sudoku(board) == false
end


text:
#Haario et al. 2001 An  adaptive  Metropolis  algorithm

#proposal function
function mcmc_proposal_func_ad(par,sig)
    return par.+rand(MvNormal(sig))
end

function mcmc_update_Σ_n_μ(x,Σ,μ,ϵ,s_d,T)

    μ_new = T*μ/(T+1)+x/(T+1)
    #Σ_new = Hermitian((T-1)/T*Σ+s_d/T*(T*μ*transpose(μ)-(T+1)*μ_new*transpose(μ_new)+x*transpose(x)+UniformScaling(ϵ)))
    Σ_new = (T-1)/T*Σ+s_d/T*(T*(μ*transpose(μ))-(T+1)*(μ_new*transpose(μ_new))+(x*transpose(x))+UniformScaling(ϵ))

    return Σ_new,μ_new
end

#n runs of the Metropolis algorithm
function mcmc_run_ad!(x_current,P_current,p,sig,n::Integer,n_chains::Integer)
    n_accept = zeros(n_chains)

    Threads.@threads for j = 1:n_chains
    #for j = 1:n_chains
        for i = 1:n
            x_new = mcmc_proposal_func_ad(x_current[j],sig[j]) #ceate new candidate
            p_new = p(x_new) #Probability of candidate
            A = min(1,p_new/P_current[j]) #calcualte acceptance ratio
            u = rand() #Probability of accepting candidate
            if u<=A
                x_current[j],P_current[j]= x_new,p_new
                n_accept[j] += 1
            end
        end
    end
    return n_accept
end

#n runs of the Metropolis algorithm
function mcmc_run_ad!(x_current,P_current,p,Σ,μ,ϵ,s_d,T,n::Integer,n_chains::Integer)
    n_accept = zeros(n_chains)

    Threads.@threads for j = 1:n_chains
    #for j = 1:n_chains
        T_local = T
        for i = 1:n
            x_new = mcmc_proposal_func_ad(x_current[j],Σ[j]) #ceate new candidate
            p_new = p(x_new) #Probability of candidate
            A = min(1,p_new/P_current[j]) #calcualte acceptance ratio
            u = rand() #Probability of accepting candidate
            if u<=A
                x_current[j],P_current[j]= x_new,p_new
                n_accept[j] += 1
            end
            Σ[j],μ[j] = mcmc_update_Σ_n_μ(reshape(x_current[j],(length(x_current[j]),1)),Σ[j],μ[j],ϵ,s_d,T_local)
            T_local += 1
        end
    end
    T += n
    return n_accept
end

function mcmc_create_Σ_n_μ!(x_current,P_current,p,sig,n::Integer,n_chains::Integer)

    μ = Array{Array{Float64,2},1}(undef,n_chains) #Sample mean
    Σ = Array{Array{Float64,2},1}(undef,n_chains) #Sample Covariance matrix

    Threads.@threads for j = 1:n_chains
    #for j = 1:n_chains

        samples = Array{Float64,2}(undef,n,length(x_current[j]))

        for i = 1:n
            x_new = mcmc_proposal_func_ad(x_current[j],sig[j]) #ceate new candidate
            p_new = p(x_new) #Probability of candidate
            A = min(1,p_new/P_current[j]) #calcualte acceptance ratio
            u = rand() #Probability of accepting candidate
            if u<=A
                x_current[j],P_current[j]= x_new,p_new
            end
            samples[i,:] = x_current[j]
        end
        μ[j] = transpose(mean(samples, dims=1))
        Σ[j] = cov(samples)
    end
    return Σ,μ
end

function metropolis_mcmc_ad!(x_current,P_current,sig,p,s_d;
    n_samples::Integer=100,burn_in_init::Integer=100,burn_in::Integer=100,
    sample_freq::Integer=5,init_sample::Integer=300,ϵ::Float64=1.0e-9,every_n_sample_accept_check::Integer=10)

    n_chains = length(x_current)

    #Run initial burn-in
    println("Start: Init Burn-in")
    mcmc_run_ad!(x_current,P_current,p,sig,burn_in,n_chains)
    println("Done: Init Burn-in")

    #Create sample Covariance matrix matrix and mean
    println("Start: Create sample Σ and μ")
    Σ,μ = mcmc_create_Σ_n_μ!(x_current,P_current,p,sig,init_sample,n_chains)
    T = init_sample #Total number of iteration after initial burn-in
    println("Done: Create sample Σ and μ")

    #Run burn-in
    println("Start: Burn-in")
    mcmc_run_ad!(x_current,P_current,p,Σ,μ,ϵ,s_d,T,burn_in,n_chains)
    println("Done: Burn-in")

    #Create container for samples
    samples = Array{typeof(x_current[1]),2}(undef,n_samples,n_chains)
    accept_rate = zeros(n_chains)
    accept_rate_monitor = zeros(n_chains)

    println("Start: Sampling")
    #Run samples
    for i = 1:n_samples
        println("---Sample: ",i,"/",n_samples)
        #Run mcmc
        n_accept = mcmc_run_ad!(x_current,P_current,p,Σ,μ,ϵ,s_d,T,sample_freq,n_chains)
        accept_rate = accept_rate.+n_accept
        accept_rate_monitor = accept_rate_monitor.+n_accept
        #Save sample
        for j = 1:n_chains
            samples[i,j] = x_current[j]
        end

        if mod(i,every_n_sample_accept_check)==0
            println("Acceptance rate (short term): ",accept_rate_monitor/(every_n_sample_accept_check*sample_freq))
            println("Acceptance rate (long term): ",accept_rate/(i*sample_freq))
            accept_rate_monitor = zeros(n_chains)
        end
    end
    println("Done: Sampling")
    accept_rate /= n_samples*sample_freq

    return samples,accept_rate
end


text:
mutable struct DataBase
    T::Dict{String, Dict{String, Any}}
    Av::Dict{Int, Dict{String, Any}}
    update_id::Int
    K::Dict{DateTime, Dict{String, Any}}
    groups::Dict{Int, Dict{String, String}}
end


text:
using Distributed
@everywhere using StochasticDiffEq, DiffEqDevTools, Test
@everywhere using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems
@everywhere importsdeproblems()
using DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_linear, prob_sde_2Dlinear, prob_sde_additive
Random.seed!(100)
dts = 1./2.^(10:-1:2) #14->7 good plot

prob = prob_sde_linear
sim  = test_convergence(dts,prob,EM(),trajectories=Int(1e4),
                        weak_timeseries_errors=true,weak_dense_errors=true)
@test abs(sim.𝒪est[:weak_final]-1) < 0.3
@test abs(sim.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3
@test abs(sim.𝒪est[:weak_L2]-1) < 0.3
@test abs(sim.𝒪est[:weak_L∞]-1) < 0.3
sim  = test_convergence(dts,prob,SimplifiedEM(),trajectories=Int(1e5),
                        weak_timeseries_errors=true,weak_dense_errors=true)
@test abs(sim.𝒪est[:weak_final]-1) < 0.3
@test abs(sim.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim.𝒪est[:weak_l∞]-1) < 0.35
sim2 = test_convergence(dts,prob,RKMil(),trajectories=Int(1e4),
                        weak_timeseries_errors=true,dense_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,RKMil_General(),trajectories=Int(1e4),
                        weak_timeseries_errors=true,dense_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,SROCK1(),trajectories=Int(1e4),
                        weak_timeseries_errors=true,dense_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,SROCK2(),trajectories=Int(1e5),
                        weak_timeseries_errors=true,dense_errors=true)
@test abs(sim2.𝒪est[:weak_final]-2) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-2) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-2) < 0.3
sim2 = test_convergence(dts,prob,SROCKEM(strong_order_1=false),trajectories=Int(1e4),
                        weak_timeseries_errors=true,dense_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,SROCKEM(),trajectories=Int(1e4),
                        weak_timeseries_errors=true,dense_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,SKSROCK(),trajectories=Int(1e4),
                        weak_timeseries_errors=true,dense_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,SROCKC2(),trajectories=Int(1e5),
                        weak_timeseries_errors=true,dense_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-2) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-2) < 0.3

#omitting tests for incomplete methods
# sim2 = test_convergence(dts,prob,TangXiaoSROCK2(version_num=1),trajectories=Int(1e5),
#                           weak_timeseries_errors=true,dense_errors=true)
# @test abs(sim2.𝒪est[:weak_final]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l2]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l∞]-2) < 0.3
# sim2 = test_convergence(dts,prob,TangXiaoSROCK2(version_num=2),trajectories=Int(1e5),
#                           weak_timeseries_errors=true,dense_errors=true)
# @test abs(sim2.𝒪est[:weak_final]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l2]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l∞]-2) < 0.3
# sim2 = test_convergence(dts,prob,TangXiaoSROCK2(version_num=3),trajectories=Int(1e5),
#                           weak_timeseries_errors=true,dense_errors=true)
# @test abs(sim2.𝒪est[:weak_final]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l2]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l∞]-2) < 0.3
# sim2 = test_convergence(dts,prob,TangXiaoSROCK2(version_num=4),trajectories=Int(1e5),
#                           weak_timeseries_errors=true,dense_errors=true)
# @test abs(sim2.𝒪est[:weak_final]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l2]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l∞]-2) < 0.3
# sim2 = test_convergence(dts,prob,TangXiaoSROCK2(version_num=5),trajectories=Int(1e5),
#                           weak_timeseries_errors=true,dense_errors=true)
# @test abs(sim2.𝒪est[:weak_final]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l2]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l∞]-2) < 0.3

sim2 = test_convergence(dts,prob,WangLi3SMil_A(),trajectories=Int(1e4),
                        weak_timeseries_errors=true,dense_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,WangLi3SMil_B(),trajectories=Int(1e4),
                        weak_timeseries_errors=true,dense_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,WangLi3SMil_C(),trajectories=Int(1e4),
                        weak_timeseries_errors=true,dense_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,WangLi3SMil_D(),trajectories=Int(1e4),
                        weak_timeseries_errors=true,dense_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,WangLi3SMil_E(),trajectories=Int(1e4),
                        weak_timeseries_errors=true,dense_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,WangLi3SMil_F(),trajectories=Int(1e4),
                        weak_timeseries_errors=true,dense_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim3 = test_convergence(dts,prob,SRI(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim3.𝒪est[:weak_final]-2) < 0.3
@test abs(sim3.𝒪est[:weak_l2]-2) < 0.3
@test abs(sim3.𝒪est[:weak_l∞]-2) < 0.3
sim4 = test_convergence(dts,prob,SRIW1(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim4.𝒪est[:weak_final]-2) < 0.3
@test abs(sim4.𝒪est[:weak_l2]-2) < 0.3
@test abs(sim4.𝒪est[:weak_l∞]-2) < 0.3

prob = prob_sde_2Dlinear
sim  = test_convergence(dts,prob,EM(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim.𝒪est[:weak_final]-1) < 0.3
@test abs(sim.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3
sim  = test_convergence(dts,prob,SimplifiedEM(),trajectories=Int(1e5),
                        weak_timeseries_errors=true)
@test abs(sim.𝒪est[:weak_final]-1) < 0.3
@test abs(sim.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim.𝒪est[:weak_l∞]-1) < 0.35
sim2 = test_convergence(dts,prob,RKMil(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,RKMil_General(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,SROCK1(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,SROCK2(),trajectories=Int(1e5),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-2) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-2) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-2) < 0.3
sim2 = test_convergence(dts,prob,SROCKEM(strong_order_1=false),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,SROCKEM(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,SKSROCK(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,SROCKC2(),trajectories=Int(1e5),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-2) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-2) < 0.3

#omitting tests for incomplete methods
# sim2 = test_convergence(dts,prob,TangXiaoSROCK2(version_num=1),trajectories=Int(1e5),
#                    weak_timeseries_errors=true)
# @test abs(sim2.𝒪est[:weak_final]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l2]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l∞]-2) < 0.3
# sim2 = test_convergence(dts,prob,TangXiaoSROCK2(version_num=2),trajectories=Int(1e5),
#                    weak_timeseries_errors=true)
# @test abs(sim2.𝒪est[:weak_final]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l2]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l∞]-2) < 0.3
# sim2 = test_convergence(dts,prob,TangXiaoSROCK2(version_num=3),trajectories=Int(1e5),
#                    weak_timeseries_errors=true)
# @test abs(sim2.𝒪est[:weak_final]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l2]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l∞]-2) < 0.3
# sim2 = test_convergence(dts,prob,TangXiaoSROCK2(version_num=4),trajectories=Int(1e5),
#                    weak_timeseries_errors=true)
# @test abs(sim2.𝒪est[:weak_final]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l2]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l∞]-2) < 0.3
# sim2 = test_convergence(dts,prob,TangXiaoSROCK2(version_num=5),trajectories=Int(1e5),
#                    weak_timeseries_errors=true)
# @test abs(sim2.𝒪est[:weak_final]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l2]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l∞]-2) < 0.3

sim2 = test_convergence(dts,prob,WangLi3SMil_A(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,WangLi3SMil_B(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,WangLi3SMil_C(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,WangLi3SMil_D(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,WangLi3SMil_E(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,WangLi3SMil_F(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim3 = test_convergence(dts,prob,SRI(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim3.𝒪est[:weak_final]-2) < 0.3
@test abs(sim3.𝒪est[:weak_l2]-2) < 0.3
@test abs(sim3.𝒪est[:weak_l∞]-2) < 0.3
sim4 = test_convergence(dts,prob,SRIW1(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim4.𝒪est[:weak_final]-2) < 0.3
@test abs(sim4.𝒪est[:weak_l2]-2) < 0.3
@test abs(sim4.𝒪est[:weak_l∞]-2) < 0.35

prob = prob_sde_additive
sim  = test_convergence(dts,prob,EM(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim.𝒪est[:weak_final]-1) < 0.3
@test abs(sim.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3
sim  = test_convergence(dts,prob,SimplifiedEM(),trajectories=Int(1e5),
                        weak_timeseries_errors=true)
@test abs(sim.𝒪est[:weak_final]-1) < 0.3
@test abs(sim.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim.𝒪est[:weak_l∞]-1) < 0.35
sim2 = test_convergence(dts,prob,RKMil(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,RKMil_General(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,SROCK1(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,SROCK2(),trajectories=Int(1e5),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-2) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-2) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-2) < 0.3
sim2 = test_convergence(dts,prob,SROCKEM(strong_order_1=false),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,SROCKEM(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,SKSROCK(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,SROCKC2(),trajectories=Int(1e5),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-2) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-2) < 0.3

#omitting tests for incomplete methods
# sim2 = test_convergence(dts,prob,TangXiaoSROCK2(version_num=1),trajectories=Int(1e5),
#                    weak_timeseries_errors=true)
# @test abs(sim2.𝒪est[:weak_final]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l2]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l∞]-2) < 0.3
# sim2 = test_convergence(dts,prob,TangXiaoSROCK2(version_num=2),trajectories=Int(1e5),
#                    weak_timeseries_errors=true)
# @test abs(sim2.𝒪est[:weak_final]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l2]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l∞]-2) < 0.3
# sim2 = test_convergence(dts,prob,TangXiaoSROCK2(version_num=3),trajectories=Int(1e5),
#                    weak_timeseries_errors=true)
# @test abs(sim2.𝒪est[:weak_final]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l2]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l∞]-2) < 0.3
# sim2 = test_convergence(dts,prob,TangXiaoSROCK2(version_num=4),trajectories=Int(1e5),
#                    weak_timeseries_errors=true)
# @test abs(sim2.𝒪est[:weak_final]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l2]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l∞]-2) < 0.3
# sim2 = test_convergence(dts,prob,TangXiaoSROCK2(version_num=5),trajectories=Int(1e5),
#                    weak_timeseries_errors=true)
# @test abs(sim2.𝒪est[:weak_final]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l2]-2) < 0.3
# @test abs(sim2.𝒪est[:weak_l∞]-2) < 0.3

sim2 = test_convergence(dts,prob,WangLi3SMil_A(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,WangLi3SMil_B(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,WangLi3SMil_C(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,WangLi3SMil_D(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,WangLi3SMil_E(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim2 = test_convergence(dts,prob,WangLi3SMil_F(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim2.𝒪est[:weak_final]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l2]-1) < 0.3
@test abs(sim2.𝒪est[:weak_l∞]-1) < 0.3
sim3 = test_convergence(dts,prob,SRI(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim3.𝒪est[:weak_final]-2) < 0.3
@test abs(sim3.𝒪est[:weak_l2]-2) < 0.3
@test abs(sim3.𝒪est[:weak_l∞]-2) < 0.3
sim4 = test_convergence(dts,prob,SRIW1(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim4.𝒪est[:weak_final]-2) < 0.3
@test abs(sim4.𝒪est[:weak_l2]-2) < 0.3
@test abs(sim4.𝒪est[:weak_l∞]-2) < 0.35
sim5 = test_convergence(dts,prob,SRA(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim5.𝒪est[:weak_final]-2) < 0.3
@test abs(sim5.𝒪est[:weak_l2]-2) < 0.3
@test abs(sim5.𝒪est[:weak_l∞]-2) < 0.3
sim6 = test_convergence(dts,prob,SRA1(),trajectories=Int(1e4),
                        weak_timeseries_errors=true)
@test abs(sim6.𝒪est[:weak_final]-2) < 0.3
@test abs(sim6.𝒪est[:weak_l2]-2) < 0.3
@test abs(sim6.𝒪est[:weak_l∞]-2) < 0.3


text:
abstract type Generator <: Injection end
const Generators = Array{<: Generator, 1}

abstract type HydroGen <: Generator end
abstract type RenewableGen <: Generator end
abstract type ThermalGen <: Generator end

function IS.get_limits(valid_range::Union{NamedTuple{(:min,:max)}, NamedTuple{(:max,:min)}},
                       unused::T) where T <: Generator
    # Gets min and max value defined for a field,
    # e.g. "valid_range": {"min":-1.571, "max":1.571}.
    return (min = valid_range.min, max = valid_range.max, zero = 0.0)
end


text:
function lower_varname(var::Variable, idv, order)
    order == 0 && return var
    name = Symbol(var.name, :ˍ, string(idv.name)^order)
    return Variable{vartype(var)}(name)
end

function flatten_differential(O::Operation)
    @assert is_derivative(O) "invalid differential: $O"
    is_derivative(O.args[1]) || return (O.args[1], O.op.x, 1)
    (x, t, order) = flatten_differential(O.args[1])
    isequal(t, O.op.x) || throw(ArgumentError("non-matching differentials on lhs: $t, $(O.op.x)"))
    return (x, t, order + 1)
end

"""
$(TYPEDSIGNATURES)

Takes a Nth order ODESystem and returns a new ODESystem written in first order
form by defining new variables which represent the N-1 derivatives.
"""
function ode_order_lowering(sys::ODESystem)
    eqs_lowered, new_vars = ode_order_lowering(equations(sys), sys.iv, states(sys))
    return ODESystem(eqs_lowered, sys.iv, new_vars, sys.ps)
end

function ode_order_lowering(eqs, iv, states)
    var_order = OrderedDict{Variable,Int}()
    D = Differential(iv())
    diff_eqs = Equation[]
    diff_vars = Variable[]
    alge_eqs = Equation[]
    alge_vars = Variable[]

    for (i, (eq, ss)) ∈ enumerate(zip(eqs, states))
        if isequal(eq.lhs, Constant(0))
            push!(alge_vars, ss)
            push!(alge_eqs, eq)
        else
            var, maxorder = var_from_nested_derivative(eq.lhs)
            # only save to the dict when we need to lower the order to save memory
            maxorder > get(var_order, var, 1) && (var_order[var] = maxorder)
            var′ = lower_varname(var, iv, maxorder - 1)
            rhs′ = rename_lower_order(eq.rhs)
            push!(diff_vars, var′)
            push!(diff_eqs, D(var′(iv())) ~ rhs′)
        end
    end

    for (var, order) ∈ var_order
        for o in (order-1):-1:1
            lvar = lower_varname(var, iv, o-1)
            rvar = lower_varname(var, iv, o)
            push!(diff_vars, lvar)

            rhs = rvar(iv())
            eq = Differential(iv())(lvar(iv())) ~ rhs
            push!(diff_eqs, eq)
        end
    end

    # we want to order the equations and variables to be `(diff, alge)`
    return (vcat(diff_eqs, alge_eqs), vcat(diff_vars, alge_vars))
end

function rename_lower_order(O::Expression)
    isa(O, Operation) || return O
    if is_derivative(O)
        (x, t, order) = flatten_differential(O)
        return lower_varname(x.op, t.op, order)(x.args...)
    end
    return Operation(O.op, rename_lower_order.(O.args))
end