text:
function calculateSplitValue(a::MaxAbsValueSplit, fname::Symbol,number_of_char_features::Int, labellist::Array{T}, sumnumerator::Array{Float64,1}, sumdenominator::Array{Float64,1}, sumweight::Array{Float64,1}, countlistfloat::Array{Float64,1}, minweight::Float64, subs::DTSubsets) where T <: Unsigned
# here randomweight==0
# for subsets, exhaustive search with flipping members (gray code) or "increasing" subset search ({1}, {1,2}, {1,2,3}, .... {1,2,3, ....., n-1,2})
# all input lists (labellist,sumnumerator,sumdenominator,sumweight,countlistfloat) need to be sorted in the same manner
# convention, in the beginning everything is on the right side
    elementsInLeftChildBV = BitVector(undef, length(labellist));fill!(elementsInLeftChildBV, false) # indicates which classes belong to right child
    chosen_subset_bitarray = BitVector(undef, 0)
    numtot = sum(sumnumerator)
    denomtot = sum(sumdenominator)
    weighttot = sum(sumweight)
    weighttot_minw = weighttot - minweight

    sumnl = sumdl = sumwl = 0.0
    val = -Inf
    chosen_sumwl = NaN
    for i in subs
  # @show "vv=$(i)" #should be the element which switches
  # i is an index, it indicates which element of labellist will flip sides for the next calculation
        if elementsInLeftChildBV[i]
  # update elementsInLeftChildBV, i.e. toggle the value of the ith component, $=XOR
            elementsInLeftChildBV[i] = elementsInLeftChildBV[i] $ true # updating needs to occur here before we copy the array (in case it is a better split than what we have seen so far)
    # move class from left to right side
            sumnl -= sumnumerator[i]
	sumdl -= sumdenominator[i]
	sumwl -= sumweight[i]
        else
  # update elementsInLeftChildBV, i.e. toggle the value of the ith component, $=XOR
            elementsInLeftChildBV[i] = elementsInLeftChildBV[i] $ true
      # move class from right to left side
            sumnl += sumnumerator[i]
	sumdl += sumdenominator[i]
	sumwl += sumweight[i]
        end
        if (sumwl > minweight) && (weighttot_minw > sumwl) # do we have enough exposure? is the split valid?
	  # vold=valnew
	  valnew = max(abs(sumnl / sumdl), abs((numtot - sumnl) / (denomtot - sumdl)))
      # if vold==valnew;@show i,valnew,sumnl,sumdl;end;
	  if valnew > val
                val = valnew
                chosen_subset_bitarray = copy(elementsInLeftChildBV)
                chosen_sumwl = sumwl
            end
        end

    end
    if isfinite(val)
  # warning: if the labellist is not sorted 1:n we need to be careful here!
        chosen_subset = labellist[chosen_subset_bitarray]
    else
        chosen_subset = Array{UInt}(undef, 0)
    end
    return val, chosen_subset, chosen_sumwl, weighttot - chosen_sumwl
end

function calculateSplitValue(a::MaxAbsValueSplit, fname, number_of_char_features::Int, labellist::Array{T,1}, sumnumerator::Array{Float64,1}, sumdenominator::Array{Float64,1}, sumweight::Array{Float64,1}, countlistfloat::Array{Float64,1}, minweight::Float64, subs::DTSubsets, feature_column_id::Int) where T <: Unsigned
# here randomweight>0
# for subsets, exhaustive search with flipping members (gray code) or "increasing" subset search ({1}, {1,2}, {1,2,3}, .... {1,2,3, ....., n-1,2})
# all input lists (labellist,sumnumerator,sumdenominator,sumweight,countlistfloat) need to be sorted in the same manner
# convention, in the beginning everything is on the right side
    elementsInLeftChildBV = BitVector(undef, length(labellist));fill!(elementsInLeftChildBV, false) # indicates which classes belong to right child
    numtot = sum(sumnumerator)
    denomtot = sum(sumdenominator)
    weighttot = sum(sumweight)
    weighttot_minw = weighttot - minweight
    sumnl = sumdl = sumwl = 0.0
    this_splitlist = Array{Splitdef}(undef, 0)

    for i in subs
    # i is an index, it indicates which element of labellist will flip sides for the next calculation
        if elementsInLeftChildBV[i]
  # update elementsInLeftChildBV, i.e. toggle the value of the ith component, $=XOR
            elementsInLeftChildBV[i] = elementsInLeftChildBV[i] $ true # updating needs to occur here before we copy the array (in case it is a better split than what we have seen so far)
    # move class from left to right side
            sumnl -= sumnumerator[i]
	sumdl -= sumdenominator[i]
	sumwl -= sumweight[i]
        else
  # update elementsInLeftChildBV, i.e. toggle the value of the ith component, $=XOR
            elementsInLeftChildBV[i] = elementsInLeftChildBV[i] $ true
      # move class from right to left side
            sumnl += sumnumerator[i]
	sumdl += sumdenominator[i]
	sumwl += sumweight[i]
        end
        if (sumwl > minweight) && (weighttot_minw > sumwl) # do we have enough exposure? is the split valid?
            valnew = max(abs(sumnl / sumdl), abs((numtot - sumnl) / (denomtot - sumdl)))
	  # @show i,valnew,sumnl,sumdl
	  push!(this_splitlist, Splitdef(feature_column_id, copy(labellist[elementsInLeftChildBV]), valnew, sumwl, weighttot - sumwl))
        end
    end
    return this_splitlist
end

# Functions for arbitrary splits for RankOpt
function calculateSplitValue(a::MaxAbsValueSplit, number_of_char_features::Int, labels::Array{Float64,1}, labels_orig::Array{Float64,1}, labels_new::Array{Float64,1}, features::Array{Float64,1}, thresh::Float64, minweight::Float64, prem_buffer::Int, moderationfactor::Float64)
    suml = sumr = f = l = 0.0
    countl = countr = 0
    for i = 1:size(features, 1)
        l = labels[i]
        f = features[i]
        if f < thresh
            suml += l
            countl += 1
        else
            sumr += l
            countr += 1
        end
    end
    return max(abs(sumr / Float64(countr)), abs(suml / Float64(countl))), Float64(countl), Float64(countr)
end

function calculateSplitValue(a::MaxAbsValueSplit, number_of_char_features::Int, labels::Array{Float64,1}, labels_orig::Array{Float64,1}, labels_new::Array{Float64,1}, features::PooledArray, subset::Array{UInt8,1}, minweight::Float64, prem_buffer::Int, moderationfactor::Float64)
    suml = sumr = 0.0
    countl = countr = 0
    for i in 1:length(features.pda)
        if features.pda.refs[i] in subset
            suml += labels[i]
            countl += 1
        else
            sumr += labels[i]
            countr += 1
        end
    end
    return max(abs(sumr / Float64(countr)), abs(suml / Float64(countl))), Float64(countl), Float64(countr)
end


text:
using Test
using AES

mutable struct Testcase
  key::String
  iv::String
  input::String
  output::String
  encrypt::Bool
end

function loadtests(filename::String)
  println("Loading tests from $filename")
  lines = readlines(filename)

  testcases = Testcase[]
  encrypt = true
  testcase = Testcase("", "", "", "", true)
  i = 0
  for line in lines
    i += 1
    if line == ""
      continue

    elseif occursin(r"^#", line)
      continue

    elseif line == "[ENCRYPT]"
      encrypt = true

    elseif line == "[DECRYPT]"
      encrypt = false

    elseif occursin(r"^KEY",line)
      x = match(r"^KEY = ([0-9A-Fa-f]+)", line)
      testcase.key = x.captures[1]

    elseif occursin(r"^IV",line)
      x = match(r"^IV = ([0-9A-Fa-f]+)", line)
      testcase.iv = x.captures[1]

    elseif occursin(r"^PLAINTEXT",line)
      x = match(r"^PLAINTEXT = ([0-9A-Fa-f]+)", line)
      if encrypt
        testcase.input = x.captures[1]
      else
        testcase.output = x.captures[1]
      end

    elseif occursin(r"^CIPHERTEXT",line)
      x = match(r"^CIPHERTEXT = ([0-9A-Fa-f]+)", line)
      if encrypt
        testcase.output = x.captures[1]
      else
        testcase.input = x.captures[1]
      end

    end

    if testcase.key != "" &&
       testcase.input != "" &&
       testcase.output != ""

      testcase.encrypt = encrypt
      push!(testcases, testcase)
      testcase = Testcase("","","","",true)
    end
  end
  return testcases
end

for testfile in readdir("NIST-AES-Vectors")
  fn = identity

  if occursin(r"^(CBC|CFB|ECB|OFB)", testfile)
    testcases = loadtests("NIST-AES-Vectors/$testfile")

    for testcase in testcases
      if occursin(r"^CBC", testfile)
        @test AESCBC(testcase.input, testcase.key, testcase.iv, testcase.encrypt) == testcase.output
      elseif occursin(r"^CFB", testfile)
        @test AESCFB(testcase.input, testcase.key, testcase.iv, testcase.encrypt) == testcase.output
      elseif occursin(r"^ECB", testfile)
        @test AESECB(testcase.input, testcase.key, testcase.encrypt) == testcase.output
      elseif occursin(r"^OFB", testfile)
        @test AESOFB(testcase.input, testcase.key, testcase.iv) == testcase.output
      end
    end
  else
    println("skipping: $testfile")
  end
end

ext: jl


text:
"""

Helper infrastructure to compile and sample models using Stan's cmdstan binary.

$(SIGNATURES)

# Extended help

Exports:
```Julia
* `SampleModel`                        : Model structure to sample a Stan language model
* `StanModelError`                     : Exception error on cmdstan compilation failure
* `stan_sample`                        : Sample the model
* `read_samples`                       : Read the samples from .csv files
* `read_summary`                       : Read the cmdstan summary .csv file
* `stan_summary`                       : Create the stansummary .csv file
* `stan_generate_quantities`           : Simulate generated_quantities
* `read_generated_quantities`          : Read generated_quantities values
```
"""
module StanSample

using Reexport

@reexport using StanBase, DataFrames, MonteCarloMeasurements
@reexport using Tables, TableOperations, CSV

using Requires
using DocStringExtensions: FIELDS, SIGNATURES, TYPEDEF

import StanBase: stan_sample, get_cmdstan_home
import StanBase: cmdline, read_summary, stan_summary
import StanBase: RandomSeed, Init, Output, StanModelError

function __init__()
  @require MCMCChains="c7f686f2-ff18-58e9-bc7b-31028e88f75d" include("utils/mcmcchains.jl")
end

include("stanmodel/sample_types.jl")
include("stanmodel/SampleModel.jl")
include("stanrun/cmdline.jl")
include("stanrun/diagnose.jl")
include("stanrun/stan_generated_quantities.jl")
include("stansamples/read_samples.jl")
include("stansamples/read_csv_files.jl")
include("stansamples/convert_a3d.jl")
include("stansamples/read_generated_quantities.jl")
include("utils/dataframes.jl")
include("utils/namedtuples.jl")
include("utils/particles.jl")
include("utils/tables.jl")

export
  SampleModel,
  StanModelError,
  stan_sample,
  read_samples,
  read_summary,
  stan_summary,
  stan_generate_quantities,
  read_generated_quantities,
  convert_a3d,
  diagnose

end # module


text:
function Phase(arg0::GroundStation, arg1::AbsoluteDate, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::ObservableSatellite)
    return Phase((GroundStation, AbsoluteDate, jdouble, jdouble, jdouble, jdouble, ObservableSatellite), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end

function add_modifier(obj::AbstractMeasurement, arg0::EstimationModifier)
    return jcall(obj, "addModifier", void, (EstimationModifier,), arg0)
end

function compare_to(obj::ComparableMeasurement, arg0::ComparableMeasurement)
    return jcall(obj, "compareTo", jint, (ComparableMeasurement,), arg0)
end

function estimate(obj::AbstractMeasurement, arg0::jint, arg1::jint, arg2::Vector{SpacecraftState})
    return jcall(obj, "estimate", EstimatedMeasurement, (jint, jint, Vector{SpacecraftState}), arg0, arg1, arg2)
end

function get_base_weight(obj::AbstractMeasurement)
    return jcall(obj, "getBaseWeight", Vector{jdouble}, ())
end

function get_coordinates(::Type{AbstractMeasurement}, arg0::SpacecraftState, arg1::jint, arg2::jint)
    return jcall(AbstractMeasurement, "getCoordinates", TimeStampedFieldPVCoordinates, (SpacecraftState, jint, jint), arg0, arg1, arg2)
end

function get_date(obj::AbstractMeasurement)
    return jcall(obj, "getDate", AbsoluteDate, ())
end

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

function get_modifiers(obj::AbstractMeasurement)
    return jcall(obj, "getModifiers", List, ())
end

function get_observed_value(obj::AbstractMeasurement)
    return jcall(obj, "getObservedValue", Vector{jdouble}, ())
end

function get_parameters_drivers(obj::AbstractMeasurement)
    return jcall(obj, "getParametersDrivers", List, ())
end

function get_satellites(obj::AbstractMeasurement)
    return jcall(obj, "getSatellites", List, ())
end

function get_station(obj::Phase)
    return jcall(obj, "getStation", GroundStation, ())
end

function get_theoretical_standard_deviation(obj::AbstractMeasurement)
    return jcall(obj, "getTheoreticalStandardDeviation", Vector{jdouble}, ())
end

function get_wavelength(obj::Phase)
    return jcall(obj, "getWavelength", jdouble, ())
end

function is_enabled(obj::AbstractMeasurement)
    return jcall(obj, "isEnabled", jboolean, ())
end

function set_enabled(obj::AbstractMeasurement, arg0::jboolean)
    return jcall(obj, "setEnabled", void, (jboolean,), arg0)
end

function signal_time_of_flight(::Type{AbstractMeasurement}, arg0::TimeStampedFieldPVCoordinates, arg1::FieldVector3D, arg2::FieldAbsoluteDate)
    return jcall(AbstractMeasurement, "signalTimeOfFlight", RealFieldElement, (TimeStampedFieldPVCoordinates, FieldVector3D, FieldAbsoluteDate), arg0, arg1, arg2)
end

function signal_time_of_flight(::Type{AbstractMeasurement}, arg0::TimeStampedPVCoordinates, arg1::Vector3D, arg2::AbsoluteDate)
    return jcall(AbstractMeasurement, "signalTimeOfFlight", jdouble, (TimeStampedPVCoordinates, Vector3D, AbsoluteDate), arg0, arg1, arg2)
end


text:
function combine_nodes(nodes1::Union{Array{T,1},Array{T,2}},nodes2::Array{T,1}) where {T<:AbstractFloat}  # nodes1 can be a 1d or 2d array; nodes2 is a 1d array

  n1 = size(nodes1,1)
  n2 = size(nodes1,2)
  n3 = length(nodes2)

  combined_nodes = Array{T,2}(undef,n1*n3,n2+1)

  for i = 1:n3
    combined_nodes[(i-1)*n1+1:i*n1,1:n2] = nodes1
  end
  for i = 1:n1
    for j = 1:n3
      combined_nodes[(j-1)*n1+i,n2+1] = nodes2[j]
    end
  end

  return combined_nodes

end

#function combine_nodes(nodes1::Union{Array{T,1},Array{T,2}},nodes2::Array{T,1}) where {T<:AbstractFloat}  # nodes1 can be a 1d or 2d array; nodes2 is a 1d array

#  combined_nodes = [kron(ones(size(nodes2,1)),nodes1) kron(nodes2,ones(size(nodes1,1)))]

#  return combined_nodes

#end


text:
(x,DM,DeM,DoM) =cheb.chebdif(25,2);
De = DeM[:,:,1];
Do = DoM[:,:,1];
xp = x[1:13]
up = 1.0./(xp.-1.1).^4
A = [Do  Do*0.0;De*0.0 De];
B = [Do*0.0  -Do; -De -diagm(ones(13))];
AA = A*diagm([up;up])

max_stars_repo_name: diakosmos/razor

max_issues_repo_name: diakosmos/razor

max_forks_repo_name: diakosmos/razor


text:
# Lock state of a specific Component
# also has metadata regarding whether its
# locked only for specific operations
struct LockState
    locked::Bool
    _read::Bool # true by default (if false no read permission even if locked)
    _write::Bool # true by default (if false no write permission even if locked)
end

# A smart, distributed smart lock built out of bamboo
# and sticks I found in the streets of los angeles
# This basically locks used components
struct SmartLock
    r::Dict{ComponentPoolIndex, LockState}
end


text:
""" Algorithms operating on wiring diagrams.
"""
module WiringDiagramAlgorithms
export topological_sort, normalize_cartesian!, normalize_copy!,
  normalize_delete!, crossing_minimization_by_sort

using DataStructures
import LightGraphs
using Statistics: mean

using ..DirectedWiringDiagrams
import ..DirectedWiringDiagrams: set_box

# Traversal
###########

""" Topological sort of boxes in wiring diagram.

Returns a list of box IDs, excluding the outer box's input and output IDs.
"""
function topological_sort(d::WiringDiagram)::Vector{Int}
  vs = LightGraphs.topological_sort_by_dfs(graph(d))
  filter(v -> !(v in outer_ids(d)), vs)
end

# Normalization
###############

""" Put a wiring diagram for a cartesian category into normal form.

This function puts a wiring diagram representing a morphism in a free cartesian
category into normal form. Copies and deletions are simplified as much as
possible.
"""
function normalize_cartesian!(d::WiringDiagram)
  # TODO: By duality, the same algorithm will normalize a cocartesian diagram.
  # We should support cartesian, cocartesian, and bicartesian diagrams.
  normalize_copy!(normalize_delete!(d))
end

""" Normalize copies in a wiring diagram.

This function maximizes sharing of intermediate computations in a wiring
diagram where copies are natural.

This algorithm is basically the same as the congruence closure algorithm on
term graphs, in the special case of the empty relation R = ∅
(Baader & Nipkow, 1998, Term Rewriting and All That, Sec. 4.4).
The main difference is the possibility of zero or many function outputs.
"""
function normalize_copy!(d::WiringDiagram)
  # Compute equivalence classes of boxes (without modifying the diagram).
  sets = IntDisjointSets(nboxes(d)+2)
  initial = filter(box_ids(d)) do v
    all(u == input_id(d) for u in inneighbors(d,v))
  end
  for v1 in initial
    for v2 in initial
      merge_if_congruent!(d, sets, v1, v2)
    end
  end

  # Keep only the designated representative of each equivalence class.
  extra = Int[]
  for v in box_ids(d)
    rep = find_root(sets, v)
    if v != rep
      for wire in out_wires(d, v)
        add_wire!(d, Wire(set_box(wire.source, rep), wire.target))
      end
      push!(extra, v)
    end
  end
  rem_boxes!(d, extra)
  d
end

function merge_if_congruent!(d::WiringDiagram, sets::IntDisjointSets, v1::Int, v2::Int)
  if v1 == v2 || (!in_same_set(sets, v1, v2) && is_congruent(d, sets, v1, v2))
    union!(sets, v1, v2)
    for out1 in filter(v -> v != output_id(d), outneighbors(d, v1))
      for out2 in filter(v -> v != output_id(d), outneighbors(d, v2))
        merge_if_congruent!(d, sets, out1, out2)
      end
    end
  end
end

function is_congruent(d::WiringDiagram, sets::IntDisjointSets, v1::Int, v2::Int)::Bool
  box(d, v1) == box(d, v2) && all(eachindex(input_ports(box(d,v1)))) do port
    wires1, wires2 = in_wires(d,v1,port), in_wires(d,v2,port)
    n1, n2 = length(wires1), length(wires2)
    @assert n1 <= 1 && n2 <= 1 # TODO: Handle merges?
    n1 == n2 && all(zip(wires1, wires2)) do pair
      src1, src2 = pair[1].source, pair[2].source
      rep1, rep2 = find_root(sets, src1.box), find_root(sets, src2.box)
      set_box(src1, rep1) == set_box(src2, rep2)
    end
  end
end

""" Normalize deletions in a wiring diagram.

This function removes all unused intermediate computations in a wiring diagram
where deletion is natural.
"""
function normalize_delete!(d::WiringDiagram)
  unused = Set{Int}()
  for v in reverse(topological_sort(d))
    if all(wire.target.box in unused for wire in out_wires(d, v))
      push!(unused, v)
    end
  end
  rem_boxes!(d, unused)
  d
end

# Layout
########

""" Crossing minimization by sorting a univariate statistic.

The boxes in `sources` and/or `targets` are fixed and the boxes in `vs` are
permuted. A permutation `σ` of the latter is returned, such that `vs[σ]` are the
sorted box IDs. Both one-sided and two-sided crossing minimization are
supported, depending on whether just one, or both, of `sources` and `targets`
are given.

In this simple but popular heuristic algorithm, the boxes are permuted by
sorting a univariate statistic of the positions of incoming and/or outgoing
wires. Typical choices are:

- `mean`: the sample mean, yielding the "barycenter method"
- `median`: the sample median

In both cases, this algorithm has the property that if there is a permutation
with no crossings, it will find it.
"""
function crossing_minimization_by_sort(d::WiringDiagram, vs::Vector{Int};
    sources::Vector{Int}=Int[], targets::Vector{Int}=Int[],
    statistic::Function=mean)::Vector{Int}
  @assert allunique(vs) && allunique(sources) && allunique(targets)
  if isempty(sources) && isempty(targets)
    # Degenerate case: nothing to sort, so preserve original order.
    return collect(eachindex(vs))
  end

  source_coord = port_coords(d, sources, OutputPort)
  target_coord = port_coords(d, targets, InputPort)
  stats = map(vs) do v
    source_coords = mapreduce(vcat, sources; init=Int[]) do source
      Int[ source_coord(wire.source) for wire in wires(d, source, v) ]
    end
    target_coords = mapreduce(vcat, targets; init=Int[]) do target
      Int[ target_coord(wire.target) for wire in wires(d, v, target) ]
    end
    statistic(vcat(source_coords, target_coords))
  end
  sortperm(stats)
end

""" Make function mapping ports to logical coordinates.
"""
function port_coords(d::WiringDiagram, vs::Vector{Int}, kind::PortKind)
  get_ports = kind == InputPort ? input_ports : output_ports
  index = Dict(vs[i] => i for i in eachindex(vs))
  sizes = [ length(get_ports(d,v)) for v in vs ]
  offsets = cumsum(vcat([0], sizes))
  (port::Port -> offsets[index[port.box]] + port.port)
end

end


text:
#=0 -- C-SVC
1 -- nu-SVC
2 -- one-class SVM
3 -- epsilon-SVR
4 -- nu-SVR=#
baremodule SvmType
   const CSVC = 0 # Classification
   const NuSVC = 1 #Classification
   const OneClassSVM = 2 #Classification
   const EpsilonSVR = 3 # Regression
   const NuSVR =4  # Regression
end

#=0 -- linear: u'*v
1 -- polynomial: (gamma*u'*v + coef0)^degree
2 -- radial basis function: exp(-gamma*|u-v|^2)
3 -- sigmoid: tanh(gamma*u'*v + coef0)=#
baremodule SvmKernel
   const Linear = 0
   const Polynomial = 1
   const RBF = 2
   const Sigmoid = 3
end

type Svm
  nativeModel
  x
  y
  _svmType::Int32
  _kernel::Int32
  _degree::Integer
  _gamma::Float64
  _coef0::Float64
  _C::Float64
  _nu::Float64
  _p::Float64
  _cache_size::Float64
  _eps::Float64
  _shrinking::Bool
  _probabilityEstimates::Bool
  _weights::Union{Dict{Any, Float64}, Void}
  _verbose::Bool
  Svm() = new(nothing,# NaoviteModel
                nothing,# X
                nothing,# Y
                SvmType.CSVC, # SvmType
                SvmKernel.RBF, #Kernel
                Integer(3), # Degrre
                typemin(Float64),#gamma
                0.0, # coef0
                Float64(1.0), # C (Cost)
                Float64(0.5), # nu
                Float64(0.1),#p
                Float64(100.00), # cache_size
                Float64(0.001), # eps,
                true, #shrinking
                false,#probability_estimates
                nothing,#weights
                false) # verbosity
end
#=
options:
-s svm_type : set type of SVM (default 0)
	0 -- C-SVC
	1 -- nu-SVC
	2 -- one-class SVM
	3 -- epsilon-SVR
	4 -- nu-SVR
-t kernel_type : set type of kernel function (default 2)
	0 -- linear: u'*v
	1 -- polynomial: (gamma*u'*v + coef0)^degree
	2 -- radial basis function: exp(-gamma*|u-v|^2)
	3 -- sigmoid: tanh(gamma*u'*v + coef0)
-d degree : set degree in kernel function (default 3)
-g gamma : set gamma in kernel function (default 1/num_features)
-r coef0 : set coef0 in kernel function (default 0)
-c cost : set the parameter C of C-SVC, epsilon-SVR, and nu-SVR (default 1)
-n nu : set the parameter nu of nu-SVC, one-class SVM, and nu-SVR (default 0.5)
-p epsilon : set the epsilon in loss function of epsilon-SVR (default 0.1)
-m cachesize : set cache memory size in MB (default 100)
-e epsilon : set tolerance of termination criterion (default 0.001)
-h shrinking: whether to use the shrinking heuristics, 0 or 1 (default 1)
-b probability_estimates: whether to train a SVC or SVR model for probability estimates, 0 or 1 (default 0)
-wi weight: set the parameter C of class i to weight*C, for C-SVC (default 1)
=#
function train!{U<:Real,T}(model::Svm, x::Matrix{U}, y::Vector{T})
  model.y = y
  model.x = x
  if (model._svmType==nothing)
    if (eltype(y)<:Real)
      model._svm_type = SvmType.EpsilonSVR # eps-SVR
    else
      model._svm_type = SvmType.CSVC # C-SVR
    end
    if (model._gamma==typemin(Float64))
      model._gamma=1.0/size(x, 1)
    end
  end
  model.nativeModel= svmtrain(y, transpose(x);
          svm_type=model._svmType,
          kernel_type=model._kernel, degree=model._degree,
          gamma=model._gamma, coef0=model._coef0,
          C=model._C, nu=model._nu, p=model._p,
          cache_size=model._cache_size, eps=model._eps, shrinking=model._shrinking,
          probability_estimates=model._probabilityEstimates,
          weights=model._weights,
          verbose=model._verbose)
end

function predict{U<:Real}(model::Svm,x::Matrix{U})
  (predicted_labels, decision_values) = svmpredict(model.nativeModel,transpose(x))
  return predicted_labels
end


text:
"""
    function get_mean_m_h_vecs(data::Array{<:Real,2}, target::Array{<:Integer,1}, dist_func::Any,
                               sig::Real; f_type::String)::Tuple{Array{<:AbstractFloat}, Array{<:AbstractFloat}}

Get mean m and mean h values (see reference, auxiliary function). The f_type argument specifies whether the
features are continuous or discrete and can either have the value of "continuous" or "discrete".
"""
function get_mean_m_h_vecs(data::Array{<:Real,2}, target::Array{<:Integer,1}, dist_func::Any,
                           sig::Real; f_type::String)::Tuple{Array{<:AbstractFloat}, Array{<:AbstractFloat}}

    # Allocate matrices for storing results.
    mean_m = Array{Float64}(undef, size(data))
    mean_h = Array{Float64}(undef, size(data))

    # Go over samples and compute weighted mean m and h vectors.
    @inbounds for idx = 1:size(data, 1)

        if f_type == "continuous"
            # If features continuous.
            # Compute m and h matrices.
            m_vecs = abs.(data[idx:idx, :] .- data[target .!= target[idx], :])
            h_vecs = abs.(data[idx:idx, :] .- data[(1:size(data, 1) .!= idx) .& (target .== target[idx]), :])

        elseif f_type == "discrete"
            # If features discrete.
            # Compute m and h matrices.
            m_vecs = Int64.(data[idx:idx, :] .!= data[target .!= target[idx], :])
            h_vecs = Int64.(data[idx:idx, :] .!= data[(1:size(data, 1) .!= idx) .& (target .== target[idx]), :])
        else
            throw(DomainError(f_type, "f_type can only be equal to \"continuous\" or \"discrete\"."))
        end

        # Compute mean m vector for next sample.
        f_dists_other_m_i = exp.(-dist_func(data[idx:idx, :], data[target .!= target[idx], :])/sig)
        f_dists_other_sum_m_i = sum(f_dists_other_m_i)
        alpha_i = f_dists_other_m_i / (f_dists_other_sum_m_i + eps(Float64))    # Compute alpha vector.
        mean_m[idx, :] = m_vecs'*alpha_i

        # Compute mean h vector for next sample.
        f_dists_other_h_i = exp.(-dist_func(data[idx:idx, :], data[(1:size(data, 1) .!= idx) .& (target .== target[idx]), :])/sig)
        f_dists_other_sum_h_i = sum(f_dists_other_h_i)
        beta_i = f_dists_other_h_i / (f_dists_other_sum_h_i + eps(Float64))     # Compute beta vector.
        mean_h[idx, :] = h_vecs'*beta_i
    end

    # Return results.
    return mean_m, mean_h
end

"""
    get_gamma_vals(data::Array{<:Real,2}, target::Array{<:Integer,1}, dist_func::Any, sig::Real; f_type::String)

Get gamma values (see reference, auxiliary function). The f_type argument specifies whether the
features are continuous or discrete and can either have the value of "continuous" or "discrete".
"""
function get_gamma_vals(data::Array{<:Real,2}, target::Array{<:Integer,1}, dist_func::Any, sig::Real; f_type::String)

    # Allocate array for storing results.
    po_vals = Array{Float64}(undef, size(data, 1))

    # Go over samples and compute probabilities of sample being an outlier.
    @inbounds for idx = 1:size(data, 1)

        # Compute probability of n-th example being an outlier.
        if f_type == "continuous"
            # If features continuous.
            m_vecs = abs.(data[idx:idx, :] .- data[target .!= target[idx], :])
            d_vals = abs.(data[idx:idx, :] .- data)
        elseif f_type == "discrete"
            # If features discrete.
            m_vecs = Int64.(data[idx:idx, :] .!= data[target .!= target[idx], :])
            d_vals = Int64.(data[idx:idx, :] .!= data)
        else
            throw(DomainError(f_type, "f_type can only be equal to \"continuous\" or \"discrete\"."))
        end

        # Compute P_{o} value for next sample.
        f_dists_other_m_i = exp.(-dist_func(data[idx:idx, :], data[target .!= target[idx], :])/sig)
        f_dists_other_d_i = exp.(-dist_func(data[idx:idx, :], data[1:size(data, 1) .!= idx, :])/sig)
        po_vals[idx] = sum(f_dists_other_m_i)/(sum(f_dists_other_d_i) + eps(Float64))
    end

    # Return probabilities of examples being inliers.
    return permutedims(1 .- po_vals)'
end

"""
    irelief(data::Array{<:Real,2}, target::Array{<:Integer,1}, max_iter::Integer, k_width::Real, conv_condition::Real,
                 initial_w_div::Real; f_type::String="continuous")::Array{<:AbstractFloat}

Compute feature weights using I-Relief algorithm. The f_type argument specifies whether the features are continuous or discrete
and can either have the value of "continuous" or "discrete".

---
# Reference:
- Yijun Sun and Jian Li. Iterative RELIEF for feature weighting. In ICML
2006 - Proceedings of the 23rd International Conference on Machine
Learning, volume 2006, pages 913–920, 2006.
"""
function irelief(data::Array{<:Real,2}, target::Array{<:Integer,1}, max_iter::Integer=1000, k_width::Real=2.0, conv_condition::Real=1.0e-6,
                 initial_w_div::Real=-1.0, dist_func::Any=(e1, e2, w) -> sum(abs.(w.*(e1 .- e2)), dims=2); f_type::String="continuous")::Array{<:AbstractFloat}

    # If initial weight divisor argument has signal value of -1.0, set to I (number of features).
    if initial_w_div == -1.0
        initial_w_div = size(data, 2)
    end

    # Intialize convergence indicator and distance weights for features.
    convergence = false
    dist_weights = ones(Float64, 1, size(data, 2))/initial_w_div

    # Initialize iteration counter.
    iter_count = 0

    ### Main iteration loop. ###
    @inbounds while iter_count < max_iter && !convergence

        # Get mean m and mean h vals for all examples.
        mean_m_vecs, mean_h_vecs = get_mean_m_h_vecs(data, target, (e1, e2) -> dist_func(e1, e2, dist_weights), k_width, f_type=f_type)

        # Get gamma values and compute nu.
        gamma_vals = get_gamma_vals(data, target, (e1, e2) -> dist_func(e1, e2, dist_weights), k_width, f_type=f_type)

        # Get nu vector.
        nu = ((mean_m_vecs - mean_h_vecs)'*gamma_vals)'

        # Update distance weights.
        nu_clipped = clamp.(nu, 0, Inf)
        dist_weights_nxt = nu_clipped/(LinearAlgebra.norm(nu_clipped) + eps(Float64))

        # Check if convergence criterion satisfied. If not, continue with next iteration.
        if sum(abs.(dist_weights_nxt .- dist_weights)) < conv_condition
            dist_weights = dist_weights_nxt
            convergence = true
        else
            dist_weights = dist_weights_nxt
            iter_count += 1
        end
    end

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

    # Return feature ranks and last distance weights.
    return vec(dist_weights)
end


text:
function build_architecture(c::PMConstructor, json_dict)
    # Get the architecture from the options dictionary.
    options = json_dict[_options_path_]
    arch = options[:architecture]

    # If a custom function is passed - use that as an architecture
    # constructor. Otherwise, parse through the passed string to decode the
    # architecture.
    if isa(arch, Function)
        @debug "Dispatching Custom Architecture"
        return arch()
    end

    num_links = options[:num_links]

    # Perform manual dispatch based on the string.
    if arch == "Array_Asap3"
        toplevel =  asap3(Rectangular(num_links, 1))
    elseif arch == "Array_Asap4"
        toplevel =  asap4(Rectangular(num_links, 1))
    elseif arch == "Array_Asap2"
        toplevel = asap2(Rectangular(num_links, 1))
    else
        error("Unrecognized Architecture: $arch_string")
    end

    # Run some transformations on the generated architecture
    name_mappables(toplevel, json_dict)

    return toplevel
end

function ismatch(c::Component, pm_base_type)
    haskey(c.metadata, typekey()) || (return false)
    # Check the two implemented mapper attributes for processors
    if pm_base_type == "Processor_Core"
        return isproc(c)
    elseif pm_base_type == "Memory_Core"
        return ismemory(c)
    elseif pm_base_type == "Input_Handler_Core"
        return isinput(c)
    elseif pm_base_type == "Output_Handler_Core"
        return isoutput(c)
    end

    return false
end

function name_mappables(toplevel::TopLevel, json_dict)
    warnings_given = 0
    warning_limit = 10
    for core in json_dict["array_cores"]
        # Get the address for the core
        addr = CartesianIndex(core["address"]...)
        base_type = core["base_type"]

        # Spit out a warning is the address is not in the model.
        if !isaddress(toplevel, addr)
            # Suppress if to many warnings have been generated.
            if warnings_given < warning_limit
                # Check if core name is "nothing".
                if core["name"] === nothing
                    msg = "Address $addr not found"
                else
                    msg = "No address $addr found for core $(core["name"])."
                end
                @warn msg

                warnings_given += 1
            elseif warnings_given == warning_limit
                @warn "Suppressing further address warnings."
                warnings_given = warning_limit + 1
            end

            continue
        end

        parent = toplevel[addr]
        for path in walk_children(parent)
            component = parent[path]
            # Try to match the base_type of the Project_Manager core with
            # its type in the Mapper's model.
            #
            # This will break if multiple cores of the same type exist in the
            # same tile.
            if ismatch(component, base_type)
                # Set the metadata for this component and exit
                component.metadata["pm_name"] = get(core, "name", missing)
                component.metadata["pm_type"] = get(core, "type", missing)
                break
            end
        end
    end
end


text:
# see also https://github.com/tlienart/Franklin.jl/issues/330
@testset "locvar" begin
    s = raw"""
        @def va = 5
        @def vb = 7
        ```julia:ex
        #hideall
        println(locvar("va")+locvar("vb"))
        ```
        \output{ex}
        """ |> fd2html_td
    @test isapproxstr(s, """

12

        """)
end

@testset "shortcut" begin
    a = """
        A
        ```!
        x = 1
        ```
        B
        ```!
        print(x)
        ```
        """ |> fd2html_td
    @test isapproxstr(a, """

A

            x = 1

1

B

print(x)

1

            """)
end

@testset "#697" begin
    a = """
        A
        ```!
        x = :ab
        ```
        B
        ```!
        x = :(a+b)
        ```
        C
        """ |> fd2html_td
    @test isapproxstr(a, """

A

x = :ab

:ab

B

x = :(a+b)

:(a + b)

C

        """)
end


text:
# Helper functions

"""
    expw = normexp!(logw,expw)

Compute expw .= exp.(logw)/sum(exp.(logw)). Uses the log-sum-exp trick to control for overflow in exp function.

Smoothing adds constant value to each value prior to normalization (useful to avoid zero probabilities).
"""
function normexp!(logw,expw,smoothing=0.0)
    offset = maximum(logw)
    expw .= exp.(logw .- offset) .+ smoothing
    s = sum(expw)
    expw .*= 1/s
end
# function logsumexp!(w,we)
#     offset = maximum(w)
#     we .= exp.(w .- offset)
#     s = sum(we)
#     w .-= log(s) + offset # this computes logw
#     we .*= 1/s
# end
"""
    logsumexp(x)

Returns `log(sum(exp.(x)))`. Uses a numerically stable algorithm.

# References:
- https://arxiv.org/pdf/1412.8695.pdf eq 3.8
- https://discourse.julialang.org/t/fast-logsumexp/22827/7?u=baggepinnen for stable logsumexp
"""
function logsumexp(x::AbstractVector{<:Real})
    offset, maxind = findmax(x) # offset controls for overflow
    w = exp.(x .- offset) # Note: maximum(w) = 1
    Σ = sum_all_but(w, maxind) # Σ = sum(w)-1
    log1p(Σ) + offset #, Σ+1 # log1p controls for underflow
end
"""
    sum_all_but(x,i)

Computes `sum(x) - x[i]`.
"""
function sum_all_but(x,i)
    x[i] -= 1
    s = sum(x)
    x[i] += 1
    s
end


text:
"""
    cache_name(vars...)

Make a unique name based on the vars names.

# Examples

```julia
cache_name("test","var")
```
"""
function cache_name(vars...)
    "_cache_" * bytes2hex(sha1(join([vars...])))
end


text:
# Before running this, please make sure to activate and instantiate the environment
# corresponding to [this `Project.toml`](https://raw.githubusercontent.com/alan-turing-institute/MLJTutorials/master/Project.toml) and [this `Manifest.toml`](https://raw.githubusercontent.com/alan-turing-institute/MLJTutorials/master/Manifest.toml)
# so that you get an environment which matches the one used to generate the tutorials:
#
# ```julia
# cd("MLJTutorials") # cd to folder with the *.toml
# using Pkg; Pkg.activate("."); Pkg.instantiate()
# ```

# ## Machine type vs Scientific Type## ### Why make a distinction?## When analysing data, it is important to distinguish between## * *how the data is encoded* (e.g. `Int`), and# * *how the data should be interpreted* (e.g. a class label, a count, ...)## How the data is encoded will be referred to  as the **machine type** whereas how the data should be interpreted  will  be referred to   as the **scientific type** (or `scitype`).## In some cases, this may be un-ambiguous, for instance if you have a vector of floating point values, this should usually be interpreted as a continuous feature (e.g.: weights, speeds, temperatures, ...).## In many other cases however, there may be ambiguities, we list a few examples below:## * A vector  of `Int` e.g. `[1, 2, ...]` which should be interpreted as categorical labels,# * A vector of `Int` e.g. `[1, 2,  ...]` which should be interpreted as count data,# * A vector of `String` e.g. `["High", "Low", "High", ...]` which should  be  interpreted as ordered categorical labels,# * A vector of `String` e.g. `["John", "Maria", ...]` which should not interpreted as informative data,# * A vector of floating points `[1.5,  1.5, -2.3, -2.3]` which should be interpreted as categorical data (e.g. the few possible values of some setting), etc.## ### The Scientific Types## The package [ScientificTypes.jl](https://github.com/alan-turing-institute/ScientificTypes.jl) defines a barebone type hierarchy which can be used to indicate how a particular feature should be interpreted; in particular:## ```plaintext# Found# ├─ Known# │  ├─ Textual# │  ├─ Finite# │  │  ├─ Multiclass# │  │  └─ OrderedFactor# │  └─ Infinite# │     ├─ Continuous# │     └─ Count# └─ Unknown# ```## A *scientific type convention* is a specific implementation indicating how machine types can be related to scientific types. It may also provide helper functions to convert data to a given scitype.## The convention used in MLJ is implemented in [MLJScientificTypes.jl](https://github.com/alan-turing-institute/ScientificTypes.jl).# This is what we will use throughout; you never need to use ScientificTypes.jl# unless you intend to implement your own scientific type convention.## ### Inspecting the scitype## The `schema` function
using RDatasets, MLJScientificTypes
boston = dataset("MASS", "Boston")
sch = schema(boston)

# In this cases, most of the variables have a (machine) type `Float64` and# their default  interpretation is `Continuous`.# There is also `:Chas`, `:Rad` and `:Tax` that have a (machine) type  `Int64`# and their default interpretation is `Count`.## While the interpretation as `Continuous` is usually fine, the interpretation# as `Count` needs a bit more attention.# For instance note that:
unique(boston.Chas)

# so even  though it's got a machine type of `Int64` and consequently a# default  interpretation of `Count`, it would be more appropriate to interpret# it as an `OrderedFactor`.## ### Changing the scitype## In order to re-specify the scitype(s) of  feature(s) in a dataset, you can# use the `coerce` function and  specify pairs of variable name and  scientific# type:
boston2 = coerce(boston, :Chas => OrderedFactor);

# the effect of this is to convert the `:Chas` column to an ordered categorical# vector:
eltype(boston2.Chas)

# corresponding to the `OrderedFactor` scitype:
elscitype(boston2.Chas)

# You can also specify multiple pairs in one shot with `coerce`:
boston3 = coerce(boston, :Chas => OrderedFactor, :Rad => OrderedFactor);

# ### String and Unknown## If a feature in  your dataset has String elements, then the  default scitype# is `Textual`; you can either choose to  drop  such columns or to coerce them# to categorical:
feature = ["AA", "BB", "AA", "AA", "BB"]
elscitype(feature)

# which you can coerce:
feature2 = coerce(feature, Multiclass)
elscitype(feature2)

# ## Tips and tricks## ### Type to Type coercion## In  some cases you will want to reinterpret all features currently# interpreted as some scitype `S1` into some other scitype `S2`.# An example  is if some features are currently interpreted as `Count` because# their original type was `Int` but you  want  to  consider all such as# `Continuous`:
data = select(boston, [:Rad, :Tax])
schema(data)

# let's coerce from `Count` to `Continuous`:
data2 = coerce(data, Count => Continuous)
schema(data2)

# ### Autotype## A last useful tool is `autotype` which allows you to specify *rules* to# define the interpretation of features automatically.# You can code your own rules but there are three useful ones that are pre-# coded:## * the `:few_to_finite` rule which checks how many unique entries are present# in a vector and if there are "few" suggests a categorical type,# * the `:discrete_to_continuous` rule converts `Integer` or `Count` to# `Continuous`# * the `:string_to_multiclass` which returns `Multiclass` for any string-like# column.## For instance:
boston3 = coerce(boston, autotype(boston, :few_to_finite))
schema(boston3)

# You can also specify multiple rules, see [the docs](https://alan-turing-institute.github.io/MLJScientificTypes.jl/stable/#Automatic-type-conversion-for-tabular-data-1) for more information.
# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
module CLI # Bukdu

import ..Bukdu: Routing, Naming

"""
    CLI.routes()

Showing the routing table.
"""
function routes()
    A = Routing.store[:routing_tables]
    isempty(A) && return
    ncols = 5 # verb url C action pipe
    nrows = Int(length(A)/ncols)
    rt = reshape(A, ncols, nrows)
    paddings = maximum((length ∘ string).(rt), dims=2) .+ 2
    function f(idx, el, lastcolumn)
        if idx == lastcolumn
            el
        else
            rpad(el, paddings[idx])
        end
    end
    for rowidx in 1:nrows
        row = rt[:, rowidx]
        lastcolumn = isempty(row[ncols]) ? ncols-1 : ncols
        print.([f(idx, el, lastcolumn) for (idx, el) in enumerate(row[1:lastcolumn])])
        println()
    end
end

end # module Bukdu.CLI

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
struct Record{Tf<:Tuple, To}
    fields::Tf
end

function Record(t::T) where T<:Tuple
    To = Tuple{map(fieldtype, t)...}
    #Tov = Tuple{map(s->Vector{s},map(fieldtype, t))...}
    Record{T, To}(t)
end

# for dispatch on N
include_string(TextParse, "const RecN{N,U} = Record{T,U} where T<:NTuple{N, Any}")

@generated function tryparsenext(r::RecN{N, To}, str, i, len, opts=default_opts) where {N, To}
    quote
        R = Nullable{To}
        i > len && @goto error

        Base.@nexprs $N j->begin
            @chk2 (val_j, i) = tryparsenext(r.fields[j], str, i, len, opts)
        end

        @label done
        return R(Base.@ntuple $N val), i

        @label error
        R(), i
    end
end

const PARSE_SUCCESS = 0x00
const PARSE_ERROR   = 0x01

function gen_1parsesetindex(j, fieldexpr, colexpr)
    val_j = Symbol(:val, j)
    quote
        err_field = $j
        @chk2 $val_j, ii = tryparsenext($fieldexpr, str, i, len, opts)
        err = setcell!($colexpr, row, $val_j, str)
        if err != PARSE_SUCCESS
            err_code = err
            @goto error
        end
        i = ii
    end
end

@generated function tryparsesetindex(r::RecN{N,To}, str::AbstractString, i::Int, len::Int, columns::Tuple, row::Int, opts) where {N,To}
    fldtypes = Base.fieldtype(r, 1).parameters
    coltypes = columns.parameters
    fieldparsers = []
    j = 1
    while j <= N
        ft = fldtypes[j]
        ct = coltypes[j]
        rl = 1
        for k = j+1:N
            if fldtypes[k] == ft && coltypes[k] == ct
                rl += 1
            else
                break
            end
        end
        if rl > 2
            body = gen_1parsesetindex(:jj, :(r.fields[jj]::($ft)), :(columns[jj]::($ct)))
            push!(fieldparsers,
                  quote
                  for jj = $j:$(j+rl-1); $body; end
                  end)
            j += rl
        else
            push!(fieldparsers, gen_1parsesetindex(j, :(r.fields[$j]), :(columns[$j])))
            j += 1
        end
    end
    R = Result{Int, Tuple{Int,Int,Int,UInt8}}
    quote
        err_field = 1
        ii = i
        err_code = PARSE_ERROR
        i > len && @goto error

        $(fieldparsers...)

        @label done
        return $R(true, i)

        @label error
        $R(false, (ii, i, err_field, err_code)) # error char, start of error field, error field
    end
end

@inline function setcell!(col, i, val, str)
    col[i] = val
    PARSE_SUCCESS
end

@inline function setcell!(col::Array{String,1}, i, val::StrRange, str)
    col[i] = alloc_string(str, val)
    PARSE_SUCCESS
end

@inline Base.@propagate_inbounds function setcell!(col::StringVector, i, val::StrRange, str)
    col[i] = WeakRefString(pointer(str, val.offset + 1), val.length)
    PARSE_SUCCESS
end

# Weird hybrid of records and fields

struct UseOne{T,R<:Record,use} <: AbstractToken{T}
    record::R
end

fieldtype(::UseOne{T}) where {T} = T

function UseOne(fields::Tuple, use)
    r = Record(fields)
    UseOne{fieldtype(fields[use]), typeof(r), use}(r)
end
getthing(x, ::Type{Val{n}}) where {n} = x[n]
function tryparsenext(f::UseOne{T,S,use}, str, i, len, opts=default_opts) where {T,S,use}
    R = Nullable{T}
    @chk2 xs, i = tryparsenext(f.record, str, i, len, opts)

    @label done
    return R(getthing(xs, Val{use})), i

    @label error
    return R(), i
end

struct Repeated{F, T, N}
    field::F
end

Repeated(f::F, n) where {F} = Repeated{F, fieldtype(f), n}(f)

fieldtype(::Repeated{F,T,N}) where {F,T,N} = NTuple{N,T}

@generated function tryparsenext(f::Repeated{F,T,N}, str, i, len, opts=default_opts) where {F,T,N}
    quote
        R = Nullable{NTuple{N,T}}
        i > len && @goto error

        # pefect candidate for #11902
        Base.@nexprs $N j->begin
            @chk2 (val_j, i) = tryparsenext(f.field, str, i, len, opts)
        end

        @label done
        return R(Base.@ntuple $N val), i

        @label error
        R(), i
    end
end


text:
using Documenter, NeuralQuantumState

makedocs(modules = [NeuralQuantumState],
    authors= "Manu Compen",
    sitename="NeuralQuantumState.jl")

deploydocs(
    repo = "github.com/mcompen/NeuralQuantumState.jl.git",
    versions = ["stable" => "v^", "v#.#"]
)


text:
module Cluster
export cluster

import Random: randperm, MersenneTwister

srand = MersenneTwister

"""
    cluster(x, neigh, eps=1e-5, seed=42)

Find a partition of x such that
`abs(x[i] - x[j]) < eps`
for all edges `j in neighidx[i]` where `i` and `j` are in the same
partition.
"""
function cluster(
    x::Array,
    neighidx::N;
    eps::Float64=1e-5,
    seed::Int=42,
)::PartitionIndex where {N}
    n = length(x)
    s = Vector{Int}()
    sizehint!(s, n)
    partition = fill(-1, n)
    nextp = 0
    for i in randperm(srand(seed), n)
        if partition[i] < 0     # unexplored?
            push!(s, i)
            nextp += 1
            while !isempty(s)
                v = pop!(s)
                partition[v] = nextp
                for u in neighidx[v]
                    if abs(x[v] - x[u]) < eps && partition[u] < 0
                        push!(s, u)
                    end
                end
            end
        end
    end
    return PartitionIndex(partition)
end

struct PartitionIndex
    I::Vector{Int}
    V::Vector{Int}
end

function PartitionIndex(ps::Vector{Int})
    k = maximum(ps)
    n = length(ps)
    V = Vector{Int}(undef, n)

    I = zeros(Int, k+1)
    for p in ps
        I[p] += 1
    end

    # prefix sum
    acc = 1
    j = 0
    jj = I[1]
    for i = 1:k
        I[i] = acc
        acc += j
        j, jj = jj, I[i+1]
    end
    I[k+1] = acc

    for (i, p) in enumerate(ps)
        V[I[p+1]] = i
        I[p+1] += 1
    end

    return PartitionIndex(I, V)
end

@inline Base.getindex(c::PartitionIndex, j::Int) =
    view(c.V, c.I[j]:c.I[j+1]-1)

Base.length(c::PartitionIndex) = length(c.I) -1

end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
__precompile__()

module PAPA

export papa_reconstruction, papa_reconstruction_debug

using LinearAlgebra, LeastSquaresOptim

greet() = print("Welcome to PAPA world!")
include("Process_Solver.jl")
include("SupportFunctions.jl")

"""
    papa_reconstruction(N, Npairs, pair_order, initial_chilist, sigmaVec[, CP_penalty, f_tol, iter])

Function for PAPA bootstrapping from a list of χ matrices to a larger χ matrix.

# Arguments
- `N::Integer`: the number of qubits
- `Npairs::Integer`: the number of qubit pairs

# Examples
```jldoctest
julia>

```
"""
function papa_reconstruction(N::Int,Npairs::Int,pair_order::Array{Int,2},initial_chilist::Vector{Float64},sigmaVec::Vector{Float64};CP_penalty = 1.0,x_tol=1.0e-7,f_tol=1.0e-7,iter=1000000,alg=LevenbergMarquardt(),flag=true)

    setup_tuple = setup_problem(N,Npairs,pair_order)

    Nel_tot = convert(Int,2*Npairs*16*17/2 + 16*2*Npairs + Npairs)

    function min_fun!(F::Vector{Float64},x::Vector{Float64})
        process_solver!(F,x,sigmaVec,setup_tuple,CP_penalty)
    end

    # function min_fun(x)
    #     return process_solver(x,sigmaVec,setup_tuple,CP_penalty)
    # end

    # result  = nlsolve(min_fun!,initial_chilist;xtol=x_tol,ftol=f_tol,iterations=iter)
    # chi_list_final = result.zero
    if flag
        print("Optimization Beginning\n")
    end
    lsa_papa = LeastSquaresProblem(x = initial_chilist,f! = min_fun!,output_length = Nel_tot)
    result = optimize!(lsa_papa,alg;x_tol=x_tol,f_tol=f_tol,iterations=iter)
    chilist = result.minimizer
    # optimize(min_fun,initial_chilist)

    dim2q = 16
    chi_PAPA = zeros(ComplexF64,Npairs,dim2q,dim2q)
    chitemp = zeros(ComplexF64,16,16)
    Nel = 136
    for nn = 1:1:Npairs
        fill!(chitemp,0.)
        start = 1 + (nn-1)*Nel
        stop = nn*Nel;
        chitemp[triu(trues(dim2q,dim2q))] = chilist[start:1:stop] + 1im*chilist[(start+Npairs*Nel):1:(stop+Npairs*Nel)]
        chitemp[:,:] = chitemp+chitemp'
        chitemp[:,:] = chitemp - diagm(0 => diag(chitemp))/2
        chi_PAPA[nn,:,:] = chitemp[:,:]
    end

    return (chi_PAPA, result.ssr)
end

function papa_reconstruction_debug(N::Int,Npairs::Int,pair_order::Array{Int,2},initial_chilist::Vector{Float64},sigmaVec::Vector{Float64};CP_penalty = 1.0,x_tol=1.0e-7,f_tol=1.0e-7,iter=1000000,alg=LevenbergMarquardt(),flag=true)

    setup_tuple = setup_problem(N,Npairs,pair_order)

    Nel_tot = convert(Int,2*Npairs*16*17/2 + 16*2*Npairs + Npairs)

    function min_fun!(F::Vector{Float64},x::Vector{Float64})
        process_solver!(F,x,sigmaVec,setup_tuple,CP_penalty)
    end

    # function min_fun(x)
    #     return process_solver(x,sigmaVec,setup_tuple,CP_penalty)
    # end

    # result  = nlsolve(min_fun!,initial_chilist;xtol=x_tol,ftol=f_tol,iterations=iter)
    # chi_list_final = result.zero
    if flag
        print("Optimization Beginning\n")
    end
    lsa_papa = LeastSquaresProblem(x = initial_chilist,f! = min_fun!,output_length = Nel_tot)
    result = optimize!(lsa_papa,alg;x_tol=x_tol,f_tol=f_tol,iterations=iter)
    chilist = result.minimizer
    # optimize(min_fun,initial_chilist)

    dim2q = 16
    chi_PAPA = zeros(ComplexF64,Npairs,dim2q,dim2q)
    chitemp = zeros(ComplexF64,16,16)
    Nel = 136
    for nn = 1:1:Npairs
        fill!(chitemp,0.)
        start = 1 + (nn-1)*Nel
        stop = nn*Nel;
        chitemp[triu(trues(dim2q,dim2q))] = chilist[start:1:stop] + 1im*chilist[(start+Npairs*Nel):1:(stop+Npairs*Nel)]
        chitemp[:,:] = chitemp+chitemp'
        chitemp[:,:] = chitemp - diagm(0 => diag(chitemp))/2
        chi_PAPA[nn,:,:] = chitemp[:,:]
    end

    return (chi_PAPA, result)
end

pair_order_1 = [1 2]

# Experimental data
sigma12_I = zeros(16,16)
max_ent_state!(4;mrho=sigma12_I)
sigma12_I = convert(Array{ComplexF64,2},sigma12_I)

# Conver to vector of data

sigmaVec_I = [real(sigma12_I[triu(trues(16,16))]);imag(sigma12_I[triu(trues(16,16))])]

# Initial state
chi12_idI = zeros(ComplexF64,16,16)

for i1 = 1:1:4
    for i2 = 1:1:4
        for i3 = 1:1:4
            for i4 = 1:1:4
                chi12_idI[4*(i1-1)+i2,4*(i3-1)+i4] = 4*sigma12_I[4*(i2-1)+i1,4*(i4-1)+i3]
            end
        end
    end
end

initial_chilist_I = [real(chi12_idI[triu(trues(16,16))]);imag(chi12_idI[triu(trues(16,16))])]

function __init__()
    print("Precompiling for speed\n")
    papa_reconstruction(2,1,pair_order_1,initial_chilist_I,sigmaVec_I,CP_penalty = 1.0,x_tol=1.0e-7,f_tol=1.0e-7,iter=1000000,alg=LevenbergMarquardt(),flag=false);
    papa_reconstruction(2,1,pair_order_1,initial_chilist_I,sigmaVec_I,CP_penalty = 1.0,x_tol=1.0e-7,f_tol=1.0e-7,iter=1000000,alg=Dogleg(),flag=false);
    nothing
end

end # module


text:
## AbstractSparseTensor type and functions/methods

"""
    AbstractSparseTensor{Tv,Ti<:Integer,N} <: AbstractArray{Tv,N}

Abstract supertype for `N`-dimensional sparse tensors
with elements of type `Tv` and indices of type `Ti`.
"""
abstract type AbstractSparseTensor{Tv,Ti<:Integer,N} <: AbstractArray{Tv,N} end

## Overloads for specializing outputs

show(io::IO, A::AbstractSparseTensor) = invoke(show, Tuple{IO,Any}, io, A)
function show(io::IO, ::MIME"text/plain", A::AbstractSparseTensor)
    nstored, N = numstored(A), ndims(A)

    # Print summary
    summary(io, A)
    iszero(nstored) && return
    print(io, ":")

    # Print stored entries
    entrylines = get(io, :limit, false) ? displaysize(io)[1] - 4 : typemax(Int)
    pad = map(ndigits, size(A))
    if entrylines >= nstored                    # Enough space to print all the stored entries
        for (ind, val) in storedpairs(A)
            _print_ln_entry(io, pad, ind, val)
        end
    elseif entrylines <= 0                      # No space to print any of the stored entries
        print(io, " \u2026")
    elseif entrylines == 1                      # Only space to print vertical dots
        print(io, '\n', " \u22ee")
    elseif entrylines == 2                      # Only space to print first stored entry
        ind, val = first(storedpairs(A))
        _print_ln_entry(io, pad, ind, val)
        print(io, '\n', ' '^(3 + sum(pad) + 2 * (N - 1) + 3), '\u22ee')
    else                                        # Print the stored entries in two chunks
        # Fetch vectors of entries
        inds, vals = storedindices(A), storedvalues(A)

        # First chunk
        prechunk = div(entrylines - 1, 2, RoundUp)
        for ptr in 1:prechunk
            _print_ln_entry(io, pad, inds[ptr], vals[ptr])
        end

        # Dots
        print(io, '\n', ' '^(3 + sum(pad) + 2 * (N - 1) + 3), '\u22ee')

        # Second chunk
        postchunk = div(entrylines - 1, 2, RoundDown)
        for ptr in nstored-postchunk+1:nstored
            _print_ln_entry(io, pad, inds[ptr], vals[ptr])
        end
    end
end
function _print_ln_entry(io::IO, pad::NTuple{N,Int}, ind::NTuple{N,<:Integer}, val) where {N}
    print(io, '\n', "  [")
    for k in 1:N
        print(io, lpad(Int(ind[k]), pad[k]))
        k == N || print(io, ", ")
    end
    print(io, "]  =  ", val)
end

function summary(io::IO, A::AbstractSparseTensor)
    invoke(summary, Tuple{IO,AbstractArray}, io, A)
    nstored = numstored(A)
    print(io, " with ", nstored, " stored ", nstored == 1 ? "entry" : "entries")
end

## Overloads for improving efficiency

findall_pure(f::Function, A::AbstractSparseTensor) =
    f(zero(eltype(A))) ? invoke(findall, Tuple{Function,AbstractArray}, f, A) : findall_stored(f, A)
findall(A::AbstractSparseTensor{Bool}) = findall_pure(identity, A)
findall(f::typeof(!iszero), A::AbstractSparseTensor) = findall_pure(f, A)
findall(f::Base.Fix2{typeof(in)}, A::AbstractSparseTensor) = findall_pure(f, A)

## AbstractSparseTensor interface

"""
    dropstored!(f::Function, A::AbstractSparseTensor{Tv,Ti,N})

Update sparse tensor `A`, dropping stored entries for which `f` is true.
The function `f` is passed one argument: the value of the entry.

For example, use `dropstored!(iszero, A)` to drop all numerical zeros.
"""
function dropstored! end

"""
    numstored(A::AbstractSparseTensor{Tv,Ti,N})

Return the number of stored entries.
Includes stored numerical zeros; use `count(!iszero,A)`
to count the number of nonzeros.
"""
function numstored end

"""
    storedindices(A::AbstractSparseTensor{Tv,Ti,N})

Return a `Vector{NTuple{N,Ti}}` of all the stored indices.
May share underlying data with `A`.
"""
function storedindices end

"""
    storedvalues(A::AbstractSparseTensor{Tv,Ti,N})

Return a `Vector{Tv}` of all the stored values.
May share underlying data with `A`.
"""
function storedvalues end

"""
    storedpairs(A::AbstractSparseTensor{Tv,Ti,N})

Return an iterator over index => value pairs for all the stored entries.
May share underlying data with `A`.
"""
function storedpairs end

## AbstractSparseTensor optional interface (internal)

"""
    findall_stored(f::Function, A::AbstractSparseTensor)

Variant of `findall(f, A)` that searches over only the stored entries.
Useful when `f(0) == false` is known in advance.
"""
findall_stored(f::Function, A::AbstractSparseTensor) =
    sort!([convert(keytype(A), CartesianIndex(ind)) for (ind, val) in storedpairs(A) if f(val)])

## Generic methods

"""
    indtype(T::Type{<:AbstractSparseTensor})
    indtype(A::AbstractSparseTensor)

Return the index type of a sparse tensor.

# Examples
```julia-repl
julia> indtype(AbstractSparseTensor{Float64,Int8,2})
Int8
```
"""
indtype(::Type{<:AbstractSparseTensor}) = Integer
indtype(::Type{<:AbstractSparseTensor{<:Any,Ti}}) where {Ti<:Integer} = Ti
indtype(A::AbstractSparseTensor) = indtype(typeof(A))

## Utilities

"""
    check_Ti(dims, Ti)

Check that the `dims` tuple and `Ti` index type are valid:
+ `dims` are nonnegative and fit in `Ti` (`0 ≤ dims[k] ≤ typemax(Ti)`)
+ corresponding length fits in `Int` (`prod(dims) ≤ typemax(Int)`)
If not, throw an `ArgumentError`.
"""
function check_Ti(dims::Dims{N}, Ti::Type) where {N}
    # Check that dims are nonnegative and fit in Ti
    maxTi = typemax(Ti)
    for k in 1:N
        dim = dims[k]
        dim >= 0 || throw(ArgumentError("the size along dimension $k (dims[$k] = $dim) is negative"))
        dim <= maxTi ||
            throw(ArgumentError("the size along dimension $k (dims[$k] = $dim) does not fit in Ti = $(Ti) (typemax($Ti) = $(typemax(Ti)))"))
    end

    # Check that corresponding length fits in Int
    len = reduce(widemul, dims)
    len <= typemax(Int) ||
        throw(ArgumentError("number of elements (length = $len) does not fit in Int (prevents linear indexing)"))
    # do not need to check that dims[k] <= typemax(Int) for CartesianIndex since eltype(dims) == Int

    return nothing
end

"""
    checkbounds_dims(Bool, dims, I...)

Return `true` if the specified indices `I` are in bounds for an array
with the given dimensions `dims`. Useful for checking the inputs to constructors.
"""
function checkbounds_dims(::Type{Bool}, dims::Dims{N}, I::Vararg{Integer,N}) where {N}
    for k in 1:N
        (1 <= I[k] <= dims[k]) || return false
    end
    return true
end

"""
    checkbounds_dims(dims, I...)

Throw an error if the specified indices `I` are not in bounds for an array
with the given dimensions `dims`. Useful for checking the inputs to constructors.
"""
function checkbounds_dims(dims::Dims{N}, I::Vararg{Integer,N}) where {N}
    checkbounds_dims(Bool, dims, I...) ||
        throw(ArgumentError("index (= $I) out of bounds (dims = $dims)"))
    return nothing
end


text:
# Download BedMachine
using Pkg
Pkg.activate(".")
using NCDatasets, Plots, Parameters

smooth_surface = 2000.0 # smoothing radius, in meters

const rw = 1000.0
const ri = 910.0

!ispath("../data") && mkdir("../data")
!ispath("../output") && mkdir("../output")

url = "file://" * homedir() * "/itet-stor/glazioarch/GlacioData/BedMachine_Antarctica/168596330/BedMachineAntarctica_2019-11-05_v01.nc"

filename = if !isfile(joinpath("../data", splitdir(url)[2]))
    print("Downloading data ...")
    download(url, joinpath("../data", splitdir(url)[2]))
    println("done.")
else
    joinpath("../data", splitdir(url)[2])
end

if !(@isdefined ds)
    print("Loading NC file ... ")
    ds = NCDataset(filename)

    # downscale = 2 finest possible on 16GB
    # downscale = 1 ok on 32GB ram (needs around 16GB)
    #
    downscale = 200

    x = ds["x"][1:downscale:end];
    y = ds["y"][1:downscale:end];
    x = x[1]:x[2]-x[1]:x[end]
    y = y[1]:y[2]-y[1]:y[end]
    mask_ = ds["mask"][1:downscale:end,1:downscale:end];
    mask = (mask_.==2) .| (mask_.==4) .| (mask_.==1); # route water over(under): grounded ice or lake Vostok or ice-free-land
    surface = ds["surface"][1:downscale:end,1:downscale:end];
    bed = convert(Matrix{Float32}, ds["bed"][1:downscale:end,1:downscale:end]);
    println("done")
end

if smooth_surface>0
    dn = round(Int, smooth_surface/step(x))
    if dn>0
        #smooth
        @warn "todo"
    end
end

surf = copy(surface);
surf[.!mask] .= NaN;

## Now run flow routing with surface and bed

using WhereTheWaterFlows, PyPlot
const WWF = WhereTheWaterFlows
# inds = [45:55,45:55]
# mask = mask[inds...]
# hydpot = (rw*bed .+ (surf.-bed)*ri)[inds...]

# hydpot = (rw*bed .+ (surf.-bed)*ri)
# @time area, slen, dir, nout, nin, pits, c, bnds = waterflows( hydpot, drain_pits=true, bnd_as_pits=true);
# area[.!mask] .= NaN;
# WWF.plotarea(1:size(area,1), 1:size(area,2), area, pits)
# if !all(getindex.(Ref(mask), pits).==0)
#     @warn "Not all pits were removed"
# end
#

using WhereTheWaterAlsoFlows
#const WWAF = WhereTheWaterAlsoFlows

D = flow_routing2D(x, reverse(y), bed, surface.-bed, mask, mask*100, plotyes=false);

nothing


text:
module TranscriptionOpt

# Import the necessary packages.
import JuMP
import MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities
using ..InfiniteOpt

include("model.jl")
include("transcribe.jl")
include("optimize.jl")

# Export transcription datatypes
export TranscriptionData, TranscriptionModel

# Export transcription methods
export is_transcription_model, transcription_data, transcription_variable,
transcription_constraint, transcription_model

end # end module


text:
"""
    HeterogeneousMetaplex{T} <: AbstractEpidemicModel

A type representing a metaplex, i.e. a metapopulation model with embedded individual networks, with arbitrary dynamics.
"""
struct HeterogeneousMetaplex{T} <: AbstractEpidemicModel
    g::Vector{<:AbstractGraph} # individual relation network
    h::AbstractSimpleGraph # population network
    D::Vector{Float64}
    dynamics::T
end

HeterogeneousMetaplex(g, h, D::Real, dynamics) = HeterogeneousMetaplex(
    g,
    h,
    fill(float(D), num_states(dynamics)),
    dynamics
    )

Metapopulation(mpx::HeterogeneousMetaplex) = Metapopulation(mpx.h, mpx.D, mpx.dynamics)

#ContactProcess(mpx::HeterogeneousMetaplex) = ContactProcess(mpx.g, mpx.dynamics)

function init_rates(mpx::HeterogeneousMetaplex, state)
    N = nv(mpx.g[1])
    a = zeros(N)
    init_rates!(a, state, mpx)
    return a
end

output_type(::HeterogeneousMetaplex) = Int

progressbar(::HeterogeneousMetaplex) = true

struct HeterogeneousMetapopulation{T} <: AbstractEpidemicModel
    h::AbstractSimpleGraph
    N::Int
    ks::Vector{Int}
    D::Vector{Float64}
    dynamics::T
end

HeterogeneousMetapopulation(h, N, ks, D::Real, dynamics) = HeterogeneousMetapopulation(h, N, ks, fill(float(D), num_states(dynamics)), dynamics)

HeterogeneousMetapopulation(h, N, ks::Int, D, dynamics) = HeterogeneousMetapopulation(h, N, fill(ks, nv(h)), D, dynamics)

function init_rates(mp::HeterogeneousMetapopulation, state)
    a = zeros(nv(mp.h))
    init_rates!(a, state, mp)
    return a
end

output_type(::HeterogeneousMetapopulation) = Int


text:
const BTGdict = GenericIncidenceList{ExVertex,Edge{ExVertex},Array{ExVertex,1},Array{Array{Edge{ExVertex},1},1}}

# BayesTree declarations
"""
$(TYPEDEF)

Data structure for the Bayes (Junction) tree, which is used for inference and constructed from a given `::FactorGraph`.
"""
mutable struct BayesTree
  bt::BTGdict
  btid::Int
  cliques::Dict{Int,Graphs.ExVertex}
  frontals::Dict{Symbol,Int}
end

function emptyBayesTree()
    bt =   BayesTree(Graphs.inclist(Graphs.ExVertex,is_directed=true),
                     0,
                     Dict{Int,Graphs.ExVertex}(),
                     #[],
                     Dict{AbstractString, Int}())
    return bt
end

"""
    $TYPEDEF

Container for upward tree solve / initialization.

TODO
- remove proceed
- more direct clique access (cliq, parent, children), for multi-process solves
"""
mutable struct CliqStateMachineContainer{BTND, T <: AbstractDFG, InMemG <: InMemoryDFGTypes}
  dfg::T
  cliqSubFg::InMemG
  tree::BayesTree
  cliq::Graphs.ExVertex
  parentCliq::Vector{Graphs.ExVertex}
  childCliqs::Vector{Graphs.ExVertex}
  forceproceed::Bool # TODO: bad flag that must be removed by refactoring sm
  incremental::Bool
  drawtree::Bool
  dodownsolve::Bool
  delay::Bool
  opts::SolverParams
  refactoring::Dict{Symbol, String}
  oldcliqdata::BTND
  logger::SimpleLogger
  CliqStateMachineContainer{BTND}() where {BTND} = new{BTND, DFG.GraphsDFG, DFG.GraphsDFG}() # NOTE JT - GraphsDFG as default?
  CliqStateMachineContainer{BTND}(x1::G,
                                  x2::InMemoryDFGTypes,
                                  x3::BayesTree,
                                  x4::Graphs.ExVertex,
                                  x5::Vector{Graphs.ExVertex},
                                  x6::Vector{Graphs.ExVertex},
                                  x7::Bool,
                                  x8::Bool,
                                  x9::Bool,
                                  x10a::Bool,
                                  x10aa::Bool,
                                  x10aaa::SolverParams,
								  x10b::Dict{Symbol,String}=Dict{Symbol,String}(),
                                  x11::BTND=emptyBTNodeData(),
                                  x13::SimpleLogger=SimpleLogger(Base.stdout) ) where {BTND, G <: AbstractDFG} = new{BTND, G, typeof(x2)}(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10a,x10aa,x10aaa,x10b,x11, x13)
end

function CliqStateMachineContainer(x1::G,
                                   x2::InMemoryDFGTypes,
                                   x3::BayesTree,
                                   x4::Graphs.ExVertex,
                                   x5::Vector{Graphs.ExVertex},
                                   x6::Vector{Graphs.ExVertex},
                                   x7::Bool,
                                   x8::Bool,
                                   x9::Bool,
                                   x10::Bool,
                                   x10aa::Bool,
                                   x10aaa::SolverParams,
                                   x11::BTND=emptyBTNodeData(),
                                   x13::SimpleLogger=SimpleLogger(Base.stdout) ) where {BTND, G <: AbstractDFG}
  #
  CliqStateMachineContainer{BTND}(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x10aa,x10aaa,Dict{Symbol,String}(),x11,x13)
end

const CSMHistory = Vector{Tuple{DateTime, Int, Function, CliqStateMachineContainer}}

"""
$(TYPEDEF)

Data structure for each clique in the Bayes (Junction) tree.
"""
mutable struct BayesTreeNodeData
  frontalIDs::Vector{Symbol}
  conditIDs::Vector{Symbol}
  inmsgIDs::Vector{Symbol} # Int
  potIDs::Vector{Symbol} # Int # this is likely redundant TODO -- remove
  potentials::Vector{Symbol}
  partialpotential::Vector{Bool}

  dwnPotentials::Vector{Symbol}
  dwnPartialPotential::Vector{Bool}

  cliqAssocMat::Array{Bool,2}
  cliqMsgMat::Array{Bool,2}
  directvarIDs::Vector{Symbol} # Int
  directFrtlMsgIDs::Vector{Symbol} # Int
  msgskipIDs::Vector{Symbol} # Int
  itervarIDs::Vector{Symbol} # Int
  directPriorMsgIDs::Vector{Symbol} # Int
  debug
  debugDwn

  # future might concentrate these four fields down to two
  # these should become specialized BeliefMessage type
  upMsg::TempBeliefMsg # Dict{Symbol, BallTreeDensity}
  dwnMsg::TempBeliefMsg # Dict{Symbol, BallTreeDensity}
  upInitMsgs::Dict{Int, TempBeliefMsg}
  downInitMsg::TempBeliefMsg

  allmarginalized::Bool
  initialized::Symbol
  upsolved::Bool
  downsolved::Bool
  initUpChannel::Channel{Symbol}
  initDownChannel::Channel{Symbol}
  solveCondition::Condition
  lockUpStatus::Channel{Int}
  lockDwnStatus::Channel{Int}
  solvableDims::Channel{Dict{Symbol, Float64}}
  statehistory::Vector{Tuple{DateTime, Int, Function, CliqStateMachineContainer}}
  BayesTreeNodeData() = new()
  BayesTreeNodeData(x...) = new(x[1],x[2],x[3],x[4],x[5],x[6],x[7],x[8],x[9],x[10],
                                x[11],x[12],x[13],x[14],x[15],x[16],x[17],x[18],x[19],x[20],
                                x[21], x[22], x[23], x[24], x[25], x[26], x[27], x[28], x[29], x[30], x[31],
                                Vector{Tuple{DateTime, Int, Function, CliqStateMachineContainer}}() )
end

# TODO -- this should be a constructor
function emptyBTNodeData()
  BayesTreeNodeData(Symbol[],Symbol[],Symbol[],
                    Symbol[],Symbol[],Bool[], # 6
                    Symbol[],Bool[],
                    Array{Bool}(undef, 0,0),
                    Array{Bool}(undef, 0,0),
                    Int[],Int[],             # 10+2
                    Int[],Int[],Int[],       # 13+2
                    nothing, nothing,        # 15+2
                    Dict{Symbol, BallTreeDensity}(),  # :null => AMP.manikde!(zeros(1,1), [1.0;], (:Euclid,))),
                    Dict{Symbol, BallTreeDensity}(),  # :null => AMP.manikde!(zeros(1,1), [1.0;], (:Euclid,))),
                    Dict{Int, TempBeliefMsg}(),
                    TempBeliefMsg(),         # 19+2
                    false, :null,
                    false, false,            # 23+2
                    Channel{Symbol}(1), Channel{Symbol}(1), Condition(), # 26+2
                    Channel{Int}(1), Channel{Int}(1),
                    Channel{Dict{Symbol,Float64}}(1) )
end

"""
$(TYPEDEF)
"""
mutable struct PotProd
    Xi::Symbol # Int
    prev::Array{Float64,2}
    product::Array{Float64,2}
    potentials::Array{BallTreeDensity,1}
    potentialfac::Vector{Symbol}
end
"""
$(TYPEDEF)
"""
mutable struct CliqGibbsMC
    prods::Array{PotProd,1}
    lbls::Vector{Symbol}
    CliqGibbsMC() = new()
    CliqGibbsMC(a,b) = new(a,b)
end
"""
$(TYPEDEF)
"""
mutable struct DebugCliqMCMC
  mcmc::Union{Nothing, Array{CliqGibbsMC,1}}
  outmsg::NBPMessage
  outmsglbls::Dict{Symbol, Symbol} # Int
  priorprods::Vector{CliqGibbsMC}
  DebugCliqMCMC() = new()
  DebugCliqMCMC(a,b,c,d) = new(a,b,c,d)
end

"""
$(TYPEDEF)
"""
mutable struct UpReturnBPType
  upMsgs::NBPMessage
  dbgUp::DebugCliqMCMC
  IDvals::Dict{Symbol, EasyMessage}
  keepupmsgs::TempBeliefMsg # Dict{Symbol, BallTreeDensity} # TODO Why separate upMsgs?
  totalsolve::Bool
  UpReturnBPType() = new()
  UpReturnBPType(x1,x2,x3,x4,x5) = new(x1,x2,x3,x4,x5)
end

"""
$(TYPEDEF)

TODO refactor msgs into only a single variable
"""
mutable struct DownReturnBPType
  dwnMsg::NBPMessage
  dbgDwn::DebugCliqMCMC
  IDvals::Dict{Symbol,EasyMessage} # Int
  keepdwnmsgs::TempBeliefMsg # Dict{Symbol, BallTreeDensity}
end

"""
$(TYPEDEF)
"""
mutable struct FullExploreTreeType{T, T2, T3 <:InMemoryDFGTypes}
  fg::T3
  bt::T2
  cliq::Graphs.ExVertex
  prnt::T
  sendmsgs::Vector{NBPMessage}
end

const ExploreTreeType{T} = FullExploreTreeType{T, BayesTree}
const ExploreTreeTypeLight{T} = FullExploreTreeType{T, Nothing}

function ExploreTreeType(fgl::G,
                         btl::BayesTree,
                         vertl::Graphs.ExVertex,
                         prt::T,
                         msgs::Array{NBPMessage,1} ) where {G <: AbstractDFG, T}
  #
  ExploreTreeType{T}(fgl, btl, vertl, prt, msgs)
end

"""
$(TYPEDEF)
"""
mutable struct MsgPassType
  fg::GraphsDFG
  cliq::Graphs.ExVertex
  vid::Symbol # Int
  msgs::Array{NBPMessage,1}
  N::Int
end


text:
using Echidna
using Test

# TEST ZDT1
x = [.1 for i in 1:30]
@show obj = ZDT1(x)
@test obj[1] == 0.1
@test obj[2] == 1.464110105645933

nobjs = 2
nvars = 30
zdt1_problem = Problem(
    2,
    [false for i in 1:nobjs],
    nvars,
    [MOGA_Real(0.0, 1.0) for i in 1:30],
    ZDT1
)
sol = Solution(
    zdt1_problem,
    [.1 for i in 1:nvars],
    false,
    [0.0, 0.0],
    0.0,
    0
)
# @show sol.objectives = ZDT1(sol)
@show evaluate!(sol)
@test sol.objectives[1] == 0.1
@test sol.objectives[2] == 1.464110105645933
@test sol.evaluated

#### TEST WITH TOURNAMENT
pop_size = 100
pop = [random_candidate(zdt1_problem) for i in 1:pop_size]
for i in 1:pop_size
    evaluate!(pop[i])
end

max_stars_repo_name: phawthorne/echidna

max_issues_repo_name: phawthorne/echidna

max_forks_repo_name: phawthorne/echidna


text:
import IterTools

"""
    list_psp(element; functional, family, core, datadir_psp)

List the pseudopotential files known to DFTK. Allows various ways
to restrict the displayed files.

# Examples
```julia-repl
julia> list_psp(family="hgh")
```
will list all HGH-type pseudopotentials and
```julia-repl
julia> list_psp(family="hgh", functional="lda")
```
will only list those for LDA (also known as Pade in this context)
and
```julia-repl
julia> list_psp(:O, core=:semicore)
```
will list all oxygen semicore pseudopotentials known to DFTK.
"""
function list_psp(element=nothing; family=nothing, functional=nothing, core=nothing,
                  datadir_psp=datadir_psp())
    # Normalize input keys
    isnothing(element)    || (element = Symbol(PeriodicTable.elements[element].symbol))
    isnothing(functional) || (functional = lowercase(functional))
    isnothing(family)     || (family = lowercase(family))

    psplist = []
    for (root, _, files) in walkdir(datadir_psp)
        root = relpath(root, datadir_psp)
        for file in files
            base, ext = splitext(file)
            ext == ".sh" && continue                    # Ignore scripts
            count(isequal('-'), base) == 1 || continue  # Need exactly one '-' in filename
            count(isequal('/'), root) == 1 || continue  # family/functional

            f_family, f_functional = split(root, "/")
            f_element, f_nvalence = split(base, "-")
            f_nvalence[1] == 'q' || continue              # Need 'q' before valence number
            f_element = Symbol(uppercase(f_element[1]) * f_element[2:end])
            haskey(PeriodicTable.elements, Symbol(f_element)) || continue

            push!(psplist, (identifier=joinpath(root, file), family=f_family,
                  functional=f_functional, element=f_element,
                  n_elec_valence=parse(Int, f_nvalence[2:end]),
                  path=joinpath(datadir_psp, root, file)))
        end
    end

    # Annotate "core" property (:other, :semicore, :fullcore)
    per_elem = IterTools.groupby(psp -> (psp.family, psp.functional, psp.element), psplist)
    psp_per_element = map(per_elem) do elgroup
        @assert length(elgroup) > 0
        if length(elgroup) == 1
            cores = [(core=:fullcore, )]
        else
            cores = append!(fill((core=:other, ), length(elgroup) - 2),
                            [(core=:fullcore, ), (core=:semicore, )])
        end
        merge.(sort(elgroup, by=psp -> psp.n_elec_valence), cores)
    end

    filter(collect(Iterators.flatten(psp_per_element))) do psp
        !isnothing(element)    && psp.element    != element    && return false
        !isnothing(functional) && psp.functional != functional && return false
        !isnothing(family)     && psp.family     != family     && return false
        !isnothing(core)       && psp.core       != core       && return false
        true
    end
end


text:
@testset "api_retry" begin
    idx = 1
    fails = [true, false]

    result = CompatHelper.api_retry() do
        retval = fails[idx]
        idx += 1

        if retval
            throw(Exception())
        else
            return idx
        end
    end

    # We ran through the loop twice
    @test result == 3
end


text:
function a_res(model::SAFTVRMieFamily, z,Vol,Temp)
    return a_mono(model,z,Vol,Temp)+a_chain(model,z,Vol,Temp)+a_assoc(model,z,Vol,Temp)
end

function a_mono(model::SAFTVRMieFamily, z,Vol,Temp)
    return a_hs(model,z,Vol,Temp)+a_disp(model,z,Vol,Temp)
end
function a_disp(model::SAFTVRMieFamily, z,Vol,Temp)
    return a_1(model,z,Vol,Temp)+a_2(model,z,Vol,Temp)+a_3(model,z,Vol,Temp)
end

function a_hs(model::SAFTVRMieFamily, z,Vol,Temp)
    ζ0   = ζn(model, z,Vol,Temp, 0)
    ζ1   = ζn(model, z,Vol,Temp, 1)
    ζ2   = ζn(model, z,Vol,Temp, 2)
    ζ3   = ζn(model, z,Vol,Temp, 3)
    NParticles = N_A*sum(z[i] for i in model.components)
    return 6*Vol/π/NParticles*(3ζ1*ζ2/(1-ζ3) + ζ2^3/(ζ3*(1-ζ3)^2) + (ζ2^3/ζ3^2-ζ0)*log(1-ζ3))
end

function ζn(model::SAFTVRMieFamily, z,Vol,Temp,n)
    return π/6*ρs(model,z,Vol,Temp)*sum(xS(model,z,Vol,Temp,i)*d(model,z,Vol,Temp,i)^n for i in model.components)
end

function ρs(model::SAFTVRMieFamily, z,Vol,Temp)
    NParticles = N_A*sum(z[i] for i in model.components)
    x   = z/sum(z[i] for i in model.components)
    m   = model.params.segment
    m̄   = sum(x[i]*m[i] for i in model.components)
    return NParticles/Vol*m̄
end

function xS(model::SAFTVRMieFamily, z,Vol,Temp,i)
    x   = z/sum(z[i] for i in model.components)
    m   = model.params.segment
    m̄   = sum(x[j]*m[j] for j in model.components)
    return x[i]*m[i]/m̄
end

function d(model::SAFTVRMieFamily, z,Vol,Temp,i)
    Temp        = Temp
    ϵ           = model.params.epsilon[union(i,i)]
    σ           = model.params.sigma[union(i,i)]
    λR          = model.params.lambdaR[union(i,i)]
    λA          = model.params.lambdaA[union(i,i)]
    u           = [0.26356031971814109102031,1.41340305910651679221800,3.59642577104072208122300,7.08581000585883755692200,12.6408008442757826594300]
    w           = [0.5217556105828086524759,0.3986668110831759274500,7.5942449681707595390e-2,3.6117586799220484545e-3,2.3369972385776227891e-5]
    θ           = CMie(model,z,Vol,Temp,i,i)*ϵ/Temp
    return σ*(1-sum(w[j]*(θ./(θ+u[j]))^(1/λR)*(exp(θ*(1/(θ./(θ+u[j]))^(λA/λR)-1))/(u[j]+θ)/λR) for j in 1:5))
end

function CMie(model::SAFTVRMieFamily, z,Vol,Temp,i,j)
    λR          = model.params.lambdaR
    λA          = model.params.lambdaA
    return (λR[union(i,j)]/(λR[union(i,j)]-λA[union(i,j)]))*(λR[union(i,j)]/λA[union(i,j)])^(λA[union(i,j)]/(λR[union(i,j)]-λA[union(i,j)]))
end

function ζ_x(model::SAFTVRMieFamily, z,Vol,Temp)
    return pi/6*ρs(model,z,Vol,Temp)*sum(xS(model,z,Vol,Temp,i)*xS(model,z,Vol,Temp,j)*(d(model,z,Vol,Temp,i)+d(model,z,Vol,Temp,j))^3/8 for i in model.components for j in model.components)
end

function a_1(model::SAFTVRMieFamily, z,Vol,Temp)
    x    = z/sum(z[i] for i in model.components)

    m    = model.params.segment
    m̄    = sum(x[j]*m[j] for j in model.components)
    return m̄/Temp*sum(xS(model,z,Vol,Temp,i)*xS(model,z,Vol,Temp,j)*a_1ij(model,z,Vol,Temp,i,j) for i in model.components for j in model.components)
end

function a_1ij(model::SAFTVRMieFamily, z,Vol,Temp,i,j)
    ϵ           = model.params.epsilon[union(i,j)]
    λR          = model.params.lambdaR[union(i,j)]
    λA          = model.params.lambdaA[union(i,j)]

    HSd         = (d(model,z,Vol,Temp,i)+d(model,z,Vol,Temp,j))/2
    x_0         = x0(model,z,Vol,Temp,i,j)
    C           = CMie(model,z,Vol,Temp,i,j)
    return 2*π*ϵ*HSd^3*C*ρs(model,z,Vol,Temp)*(x_0^λA*(a1s(model,z,Vol,Temp,λA)+B(model,z,Vol,Temp,λA,x_0))-x_0^λR*(a1s(model,z,Vol,Temp,λR)+B(model,z,Vol,Temp,λR,x_0)))
end

function a1s(model::SAFTVRMieFamily, z,Vol,Temp,λ)
    ζeff = ζ_eff(model,z,Vol,Temp,λ)
    return -1/(λ-3)*(1-ζeff/2)/(1-ζeff)^3
end

function ζ_eff(model::SAFTVRMieFamily, z,Vol,Temp,λ)
    ζx = ζ_x(model,z,Vol,Temp)
    return sum(c(model,z,Vol,Temp,λ,n)*ζx^n for n in 1:4)
end

function c(model::SAFTVRMieFamily, z,Vol,Temp,λ,n)
    A           = [[0.81096   1.7888  -37.578   92.284],
                   [1.02050  -19.341   151.26  -463.50],
                  [-1.90570   22.845  -228.14   973.92],
                   [1.08850  -6.1962   106.98  -677.64]]
    return sum(A[n][m]/λ^(m-1) for m in 1:4)
end

function B(model::SAFTVRMieFamily, z,Vol,Temp,λ,x0)
    I  = (1-x0^(3-λ))/(λ-3)
    J  = (1-(λ-3)*x0^(4-λ)+(λ-4)*x0^(3-λ))/((λ-3)*(λ-4))
    ζx = ζ_x(model,z,Vol,Temp)
    return I*(1-ζx/2)/(1-ζx)^3-9*J*ζx*(ζx+1)/(2*(1-ζx)^3)
end

function x0(model::SAFTVRMieFamily, z,Vol,Temp,i,j)
    σ           = model.params.sigma
    HSd         = (d(model,z,Vol,Temp,i)+d(model,z,Vol,Temp,j))/2
    return σ[union(i,j)]/HSd
end

function a_2(model::SAFTVRMieFamily, z,Vol,Temp)
    x    = z/sum(z[i] for i in model.components)

    m    = model.params.segment
    m̄    = sum(x[j]*m[j] for j in model.components)
    return m̄/Temp^2*sum(xS(model,z,Vol,Temp,i)*xS(model,z,Vol,Temp,j)*a_2ij(model,z,Vol,Temp,i,j) for i in model.components for j in model.components)
end

function a_2ij(model::SAFTVRMieFamily, z,Vol,Temp,i,j)
    ϵ           = model.params.epsilon[union(i,j)]
    λR          = model.params.lambdaR[union(i,j)]
    λA          = model.params.lambdaA[union(i,j)]

    HSd         = (d(model,z,Vol,Temp,i)+d(model,z,Vol,Temp,j))/2
    x_0         = x0(model,z,Vol,Temp,i,j)
    C           = CMie(model,z,Vol,Temp,i,j)
    ζx          = ζ_x(model,z,Vol,Temp)
    KHS         = (1-ζx)^4/(1+4*ζx+4*ζx^2-4*ζx^3+ζx^4)
    return π*KHS*(1+χ(model,z,Vol,Temp,i,j))*ρs(model,z,Vol,Temp)*ϵ^2*HSd^3*C^2*(x_0^(2*λA)*(a1s(model,z,Vol,Temp,2*λA)+B(model,z,Vol,Temp,2*λA,x_0))-2*x_0^(λA+λR)*(a1s(model,z,Vol,Temp,λA+λR)+B(model,z,Vol,Temp,λA+λR,x_0))+x_0^(2*λR)*(a1s(model,z,Vol,Temp,2*λR)+B(model,z,Vol,Temp,2*λR,x_0)))
end

function χ(model::SAFTVRMieFamily, z,Vol,Temp,i,j)
    λR          = model.params.lambdaR[union(i,j)]
    λA          = model.params.lambdaA[union(i,j)]

    ζstar = ζ_star(model, z,Vol,Temp)
    α     = CMie(model,z,Vol,Temp,i,j)*(1/(λA-3)-1/(λR-3))
    return f(α,1)*ζstar+f(α,2)*ζstar^5+f(α,3)*ζstar^8
end

function f(α,m)
    ϕ           = [[7.5365557, -359.440,  1550.9, -1.199320, -1911.2800,  9236.9,  10.0],
                  [-37.604630,  1825.60, -5070.1,  9.063632,  21390.175, -129430,  10.0],
                   [71.745953, -3168.00,  6534.6, -17.94820, -51320.700,  357230,  0.57],
                  [-46.835520,  1884.20, -3288.7,  11.34027,  37064.540, -315530, -6.70],
                  [-2.4679820,- 0.82376, -2.7171,  20.52142,  1103.7420,  1390.2, -8.00],
                  [-0.5027200, -3.19350,  2.0883, -56.63770, -3264.6100, -4518.2,   NaN],
                  [8.0956883,  3.70900,  0.0000,  40.53683,  2556.1810,  4241.6,   NaN]];
    return sum(ϕ[i+1][m]*α^i for i in 0:3)/(1+sum(ϕ[i+1][m]*α^(i-3) for i in 4:6))
end

function ζ_star(model::SAFTVRMieFamily, z,Vol,Temp)
    σ  = model.params.sigma
    return ρs(model,z,Vol,Temp)*π/6*sum(xS(model,z,Vol,Temp,i)*xS(model,z,Vol,Temp,j)*(σ[union(i,j)])^3 for i in model.components for j in model.components)
end

function a_3(model::SAFTVRMieFamily, z,Vol,Temp)
    x    = z/sum(z[i] for i in model.components)

    m    = model.params.segment
    m̄    = sum(x[j]*m[j] for j in model.components)
    return m̄/Temp^3*sum(xS(model,z,Vol,Temp,i)*xS(model,z,Vol,Temp,j)*a_3ij(model,z,Vol,Temp,i,j) for i in model.components for j in model.components)
end

function a_3ij(model::SAFTVRMieFamily, z,Vol,Temp,i,j)
    ϵ     = model.params.epsilon[union(i,j)]
    λR    = model.params.lambdaR[union(i,j)]
    λA    = model.params.lambdaA[union(i,j)]
    ζstar = ζ_star(model, z,Vol,Temp)
    α     = CMie(model,z,Vol,Temp,i,j)*(1/(λA-3)-1/(λR-3))
    return -ϵ^3*f(α,4)*ζstar*exp(f(α,5)*ζstar+f(α,6)*ζstar^2)
end

function a_chain(model::SAFTVRMieFamily, z,Vol,Temp)
    x       = z/sum(z[i] for i in model.components)
    m       = model.params.segment
    return -sum(x[i]*(log(g_Mie(model,z,Vol,Temp,i))*(m[i]-1)) for i in model.components)
end

function g_Mie(model::SAFTVRMieFamily,z,Vol,Temp,i)
    ϵ    = model.params.epsilon[union(i,i)]
    gHS  = g_HS(model,z,Vol,Temp,i)
    g1   = g_1(model,z,Vol,Temp,i)
    g2   = g_2(model,z,Vol,Temp,i)
    gMie = gHS*exp(ϵ/Temp*g1/gHS+(ϵ/Temp)^2*g2/gHS);
    return gMie
end

function g_HS(model::SAFTVRMieFamily,z,Vol,Temp,i)
    ζx  = ζ_x(model,z,Vol,Temp)
    x_0 = x0(model,z,Vol,Temp,i,i)
    k0  = -log(1-ζx)+(42ζx-39ζx^2+9ζx^3-2ζx^4)/(6*(1-ζx)^3)
    k1  = (ζx^4+6*ζx^2-12*ζx)/(2*(1-ζx)^3)
    k2  = -3*ζx^2/(8*(1-ζx)^2)
    k3  = (-ζx^4+3*ζx^2+3*ζx)/(6*(1-ζx)^3)
    gHS = exp(k0+x_0*k1+x_0^2*k2+x_0^3*k3)
    return gHS
end

function g_1(model::SAFTVRMieFamily,z,Vol,Temp,i)
    λR          = model.params.lambdaR[union(i,i)]
    λA          = model.params.lambdaA[union(i,i)]
    x_0 = x0(model,z,Vol,Temp,i,i)
    g1 = 3*da1(model,z,Vol,Temp,i)-CMie(model,z,Vol,Temp,i,i)*(λA*x_0^λA*(a1s(model,z,Vol,Temp,λA)+B(model,z,Vol,Temp,λA,x_0))-λR*x_0^λR*(a1s(model,z,Vol,Temp,λR)+B(model,z,Vol,Temp,λR,x_0)))
    return g1
end

function da1(model::SAFTVRMieFamily,z,Vol,Temp,i)
    λR  = model.params.lambdaR[union(i,i)]
    λA  = model.params.lambdaA[union(i,i)]
    x_0 = x0(model,z,Vol,Temp,i,i)
    return CMie(model,z,Vol,Temp,i,i)*(x_0^λA*(da1s(model,z,Vol,Temp,λA)+dB(model,z,Vol,Temp,λA,x_0))-x_0^λR*(da1s(model,z,Vol,Temp,λR)+dB(model,z,Vol,Temp,λR,x_0)))
end

function da1s(model::SAFTVRMieFamily,z,Vol,Temp,λ)
    ζeff  = ζ_eff(model,z,Vol,Temp,λ)
    dζeff = dζ_eff(model,z,Vol,Temp,λ)
    return -1/(λ-3)*((1-ζeff/2)/(1-ζeff)^3+ρs(model,z,Vol,Temp)*((3*(1-ζeff/2)*(1-ζeff)^2-0.5*(1-ζeff)^3)/(1-ζeff)^6)*dζeff);
end

function dζ_eff(model::SAFTVRMieFamily,z,Vol,Temp,λ)
    ζx = ζ_x(model,z,Vol,Temp)
    return sum(n*c(model,z,Vol,Temp,λ,n)*ζx^(n-1) for n in 1:4)*ζx/ρs(model,z,Vol,Temp)
end

function dB(model::SAFTVRMieFamily,z,Vol,Temp,λ,x_0)
    I  = (1-x_0^(3-λ))/(λ-3)
    J  = (1-(λ-3)*x_0^(4-λ)+(λ-4)*x_0^(3-λ))/((λ-3)*(λ-4))
    ζx = ζ_x(model,z,Vol,Temp)
    return (((1-ζx/2)*I/(1-ζx)^3-9*ζx*(1+ζx)*J/(2*(1-ζx)^3))+ρs(model,z,Vol,Temp)*((3*(1-ζx/2)*(1-ζx)^2
            -0.5*(1-ζx)^3)*I/(1-ζx)^6-9*J*((1+2*ζx)*(1-ζx)^3+ζx*(1+ζx)*3*(1-ζx)^2)/(2*(1-ζx)^6))*ζx/ρs(model,z,Vol,Temp));
end

function g_2(model::SAFTVRMieFamily,z,Vol,Temp,i)
        γ    = γ_corr(model::SAFTVRMieFamily,z,Vol,Temp,i)
        gMCA = g_MCA(model::SAFTVRMieFamily,z,Vol,Temp,i)
    return (1+γ)*gMCA
end

function γ_corr(model::SAFTVRMieFamily,z,Vol,Temp,i)
    ϕ    = [[7.5365557, -359.440,  1550.9, -1.199320, -1911.2800,  9236.9,  10.0],
           [-37.604630,  1825.60, -5070.1,  9.063632,  21390.175, -129430,  10.0],
            [71.745953, -3168.00,  6534.6, -17.94820, -51320.700,  357230,  0.57],
           [-46.835520,  1884.20, -3288.7,  11.34027,  37064.540, -315530, -6.70],
           [-2.4679820,- 0.82376, -2.7171,  20.52142,  1103.7420,  1390.2, -8.00],
           [-0.5027200, -3.19350,  2.0883, -56.63770, -3264.6100, -4518.2,   NaN],
            [8.0956883,  3.70900,  0.0000,  40.53683,  2556.1810,  4241.6,   NaN]];
   ϵ     = model.params.epsilon[union(i,i)]
   λR    = model.params.lambdaR[union(i,i)]
   λA    = model.params.lambdaA[union(i,i)]
   ζstar = ζ_star(model, z,Vol,Temp)
   α     = CMie(model,z,Vol,Temp,i,i)*(1/(λA-3)-1/(λR-3))
   return ϕ[1][7]*(1-tanh(ϕ[2][7]*(ϕ[3][7]-α)))*ζstar*(exp(ϵ/Temp)-1)*exp(ϕ[4][7]*ζstar+ϕ[5][7]*ζstar^2)
end

function g_MCA(model::SAFTVRMieFamily,z,Vol,Temp,i)
    λR  = model.params.lambdaR[union(i,i)]
    λA  = model.params.lambdaA[union(i,i)]
    x_0 = x0(model,z,Vol,Temp,i,i)
    ζx  = ζ_x(model,z,Vol,Temp)
    KHS = (1-ζx)^4/(1+4*ζx+4*ζx^2-4*ζx^3+ζx^4)
    return 3*da2(model,z,Vol,Temp,i)-KHS*CMie(model,z,Vol,Temp,i,i)^2*(
          λR*x_0^(2*λR)*(a1s(model,z,Vol,Temp,2*λR)+B(model,z,Vol,Temp,2*λR,x_0))-
          (λA+λR)*x_0^(λA+λR)*(a1s(model,z,Vol,Temp,λA+λR)+B(model,z,Vol,Temp,λA+λR,x_0))+
          λA*x_0^(2*λA)*(a1s(model,z,Vol,Temp,2*λA)+B(model,z,Vol,Temp,2*λA,x_0)))
end
function da2(model::SAFTVRMieFamily,z,Vol,Temp,i)
    λR   = model.params.lambdaR[union(i,i)]
    λA   = model.params.lambdaA[union(i,i)]
    x_0  = x0(model,z,Vol,Temp,i,i)
    ζx   = ζ_x(model,z,Vol,Temp)
    KHS  = (1-ζx)^4/(1+4*ζx+4*ζx^2-4*ζx^3+ζx^4)
    dKHS = -((4*(1-ζx)^3*(1+4*ζx+4*ζx^2-4*ζx^3+ζx^4)+(1-ζx)^4*(4+8*ζx-12*ζx^2+4*ζx^3))/(1+4*ζx+4*ζx^2-4*ζx^3+ζx^4)^2)*ζx/ρs(model,z,Vol,Temp);
    return 0.5*CMie(model,z,Vol,Temp,i,i)^2*(
           ρs(model,z,Vol,Temp)*dKHS*(x_0^(2*λA)*(a1s(model,z,Vol,Temp,2*λA)+B(model,z,Vol,Temp,2*λA,x_0))-
           2*x_0^(λA+λR)*(a1s(model,z,Vol,Temp,λA+λR)+B(model,z,Vol,Temp,λA+λR,x_0))+
           x_0^(2*λR)*(a1s(model,z,Vol,Temp,2*λR)+B(model,z,Vol,Temp,2*λR,x_0)))+
           KHS*(x_0^(2*λA)*(da1s(model,z,Vol,Temp,2*λA)+dB(model,z,Vol,Temp,2*λA,x_0))-
           2*x_0^(λA+λR)*(da1s(model,z,Vol,Temp,λA+λR)+dB(model,z,Vol,Temp,λA+λR,x_0))+
           x_0^(2*λR)*(da1s(model,z,Vol,Temp,2*λR)+dB(model,z,Vol,Temp,2*λR,x_0))))
end

function a_assoc(model::SAFTVRMieFamily, z, v, T)
    x = z/sum(z[i] for i in model.components)
    n_sites = model.params.n_sites
    X_iA = X_assoc(model,z,v,T)
    return sum(x[i]*sum(n_sites[i][a]*(log(X_iA[i,a])+(1-X_iA[i,a])/2) for a in keys(model.params.n_sites[i])) for i in model.components)
end

function X_assoc(model::SAFTVRMieFamily, z, v, T)
    x = z/sum(z[i] for i in model.components)
    ρ = N_A*sum(z[i] for i in model.components)/v
    n_sites = model.params.n_sites
    X_iA = Dict()
    X_iA_old = Dict()
    tol = 1.
    iter = 1

    while tol > 1e-12 && iter < 100
        for i in model.components
            for a in keys(model.params.n_sites[i])
                A = 0.
                for j in model.components
                    B = 0
                    for b in keys(model.params.n_sites[j])
                        if haskey(model.params.epsilon_assoc,Set([(i,a),(j,b)]))
                            if iter!=1
                                B+=n_sites[j][b]*X_iA_old[j,b]*Δ(model,z,v,T,i,j,a,b)
                            else
                                B+=n_sites[j][b]*Δ(model,z,v,T,i,j,a,b)
                            end
                        end
                    end
                    A += ρ*x[j]*B
                end
                if iter == 1
                    X_iA[i,a] =0.5+0.5*(1+A)^-1
                else
                    X_iA[i,a] =0.5*X_iA_old[i,a]+0.5*(1+A)^-1
                end
            end
        end
        if iter == 1
            tol = sqrt(sum(sum((1. -X_iA[i,a])^2 for a in keys(model.params.n_sites[i])) for i in model.components))
        else
            tol = sqrt(sum(sum((X_iA_old[i,a] -X_iA[i,a])^2 for a in keys(model.params.n_sites[i])) for i in model.components))
        end
        X_iA_old = deepcopy(X_iA)
        iter += 1
    end

    return X_iA
end

function Δ(model::SAFTVRMieFamily, z, v, T, i, j, a, b)
    ρR = ρs(model,z,v,T)*σx3(model,z,v,T)
    TR = T/model.params.epsilon[union(i,j)]

    c  = [0.0756425183020431	-0.128667137050961	 0.128350632316055	-0.0725321780970292	   0.0257782547511452  -0.00601170055221687	  0.000933363147191978  -9.55607377143667e-05  6.19576039900837e-06 -2.30466608213628e-07 3.74605718435540e-09
          0.134228218276565	    -0.182682168504886 	 0.0771662412959262	-0.000717458641164565 -0.00872427344283170	0.00297971836051287	 -0.000484863997651451	 4.35262491516424e-05 -2.07789181640066e-06	4.13749349344802e-08 0
         -0.565116428942893	     1.00930692226792   -0.660166945915607	 0.214492212294301	  -0.0388462990166792	0.00406016982985030	 -0.000239515566373142	 7.25488368831468e-06 -8.58904640281928e-08	0	                 0
         -0.387336382687019	    -0.211614570109503	 0.450442894490509	-0.176931752538907	   0.0317171522104923  -0.00291368915845693	  0.000130193710011706  -2.14505500786531e-06  0	                0	                 0
          2.13713180911797	    -2.02798460133021 	 0.336709255682693	 0.00118106507393722  -0.00600058423301506	0.000626343952584415 -2.03636395699819e-05	 0	                   0	                0	                 0
         -0.300527494795524	     2.89920714512243   -0.567134839686498	 0.0518085125423494	  -0.00239326776760414	4.15107362643844e-05  0	                     0	                   0	                0                    0
         -6.21028065719194	    -1.92883360342573	 0.284109761066570	-0.0157606767372364	   0.000368599073256615	0 	                  0	                     0	                   0	                0	                 0
          11.6083532818029	     0.742215544511197  -0.0823976531246117	 0.00186167650098254   0	                0	                  0	                     0	                   0	                0	                 0
         -10.2632535542427	    -0.125035689035085	 0.0114299144831867	 0	                   0	                0	                  0	                     0	                   0	                0	                 0
          4.65297446837297	    -0.00192518067137033 0	                 0	                   0	                0	                  0	                     0	                   0	                0	                 0
         -0.867296219639940	     0	                 0	                 0	                   0	                0	                  0	                     0	                   0	                0	                 0]
    I = sum(sum(c[n+1,m+1]*ρR^n*TR^m for m in 0:(10-n)) for n in 0:10)
    ϵ_assoc = model.params.epsilon_assoc[Set([(i,a),(j,b)])]
    F = (exp(ϵ_assoc/T)-1)
    K = model.params.bond_vol[Set([(i,a),(j,b)])]
    return F*K*I
end

function σx3(model::SAFTVRMieFamily, z, v, T)
    σ = model.params.sigma
    return sum(sum(xS(model,z,v,T,i)*xS(model,z,v,T,j)*σ[union(i,j)]^3 for j in model.components) for i in model.components)
end


text:
################################################# FILE DESCRIPTION #########################################################

# This file contains tests for the Vertex DataFrame.

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

import Graft: eltypes

@testset "Vertex DataFrame" begin
   g = completegraph(10)
   float_arr = rand(10)
   int_arr = rand(Int, 10)
   bool_arr = rand(Bool, 10)
   str_arr = map(string, 1:10)
   void_arr = fill(nothing, 10)

   V = g.vdata

   # Mutation on emtpy DataFrame
   addvertex!(V)
   @test size(V) == (0,0)
   rmvertex!(V, 11)
   @test size(V) == (0,0)
   rmvertex!(V, 1:10)
   @test size(V) == (0,0)

   ###
   # SETVPROP
   ###
   # Try to create a property only on a subset of properties
   try setvprop!(g, rand(1:10), 1, :p1); @test false catch @test true end
   try setvprop!(g, rand(1:10, rand(1:10)), 1, :p1); @test false catch @test true end
   try setvprop!(g, :, rand(rand(1:9)), :p1); @test false catch @test true end

   # Create new properties
   setvprop!(g, :, float_arr, :p1)
   setvprop!(g, :, int_arr, :p2)
   setvprop!(g, :, bool_arr, :p3)
   setvprop!(g, :, str_arr, :p4)
   setvprop!(g, :, nothing, :p5)

   # Test property names and types
   @test listvprops(g) == [:p1, :p2, :p3, :p4, :p5]
   @test eltypes(g.vdata) == [Float64, Int, Bool, String, Void]

   # Test property values
   @test getvprop(g, :, :p1) == float_arr
   @test getvprop(g, :, :p2) == int_arr
   @test getvprop(g, :, :p3) == bool_arr
   @test getvprop(g, :, :p4) == str_arr
   @test getvprop(g, :, :p5) == void_arr

   # Mutation on DataFrame
   @test addvertex!(g) == 11
   for vprop in listvprops(g)
      isequal(getvprop(g, 11, vprop), NA)
   end

   rmvertex!(g, 11)
   @test size(g.vdata) == (10,5)

   # Modify single vertex
   v = rand(1:10)

   setvprop!(g, v, 0.5, :p1)
   @test getvprop(g, v, :p1) == 0.5

   setvprop!(g, v, "p4", :p4)
   @test getvprop(g, v, :p4) == "p4"

   # Modify a range of vertices
   setvprop!(g, 1:5, 0, :p1)
   @test getvprop(g, 1:5, :p1) == zeros(5)

   setvprop!(g, 2:7, 2:7, :p2)
   @test getvprop(g, 2:7, :p2) == collect(2:7)
end