text:
###################################################
# EXEMPLO - resolver PL com a interface CPLEX do Julia
# * Requer CPLEX instalado e CPLEX.jl configurado no Julia
#
# Autor: Leonardo D. Secchin
# Data : 28/06/2021
# leonardosecchin.github.io/topicospo/
###################################################

using CPLEX, JuMP

function cplex_barrier(A, b, c; lb=[], ub=[], eps=1e-8)

    # Lê dimensões
    m,n = size(A)

    # Verifica se as dimensões dos dados são compatíveis
    if (length(b) != m) | (length(c) != n) | (m > n)
        error("Dados com dimensões incompatíveis!")
    end

    # limitantes padrão das variáveis, caso não fornecidos
    if isempty(lb)
        lb = zeros(n,1)
    end
    if isempty(ub)
        ub = 1.0e308 * ones(n,1)
    end

    # estrutura JuMP associada com CPLEX
    P = Model(CPLEX.Optimizer)

    ### configura CPLEX (ver https://www.ibm.com/docs/en/icos/20.1.0?topic=cplex-list-parameters)
    ### os parâmetros são da interface C

    # desliga parte do presolve
    set_optimizer_attribute(P, "CPXPARAM_Preprocessing_Presolve", CPX_OFF)
    set_optimizer_attribute(P, "CPX_PARAM_DEPIND", 0)

    # saídas na tela
    set_optimizer_attribute(P, "CPXPARAM_ScreenOutput", 1)

    # escolhe o algoritmo de barreiras (pontos interiores) do CPLEX
    set_optimizer_attribute(P, "CPXPARAM_LPMethod", CPX_ALG_BARRIER)

    # desliga crossover para Simplex (não purifica a solução para uma solução básica)
    set_optimizer_attribute(P, "CPXPARAM_SolutionType", CPX_NONBASIC_SOLN)

    # eps para convergência (o padrão já é 1e-8)
    set_optimizer_attribute(P, "CPXPARAM_Barrier_ConvergeTol", eps)

    # adiciona variáveis
    @variable(P, x[1:n])

    # define limitantes das variáveis
    for j in 1:n
        set_lower_bound(x[j], lb[j])
        set_upper_bound(x[j], ub[j])
    end

    # adiciona função objetivo
    @objective(P, Min, sum(c[j]*x[j] for j in 1:n))

    # adiciona restrições
    for i in 1:m
        @constraint(P, sum(A[i,j]*x[j] for j in 1:n) == b[i])
    end

    # resolve P
    optimize!(P)

    # retorna solução
    # como capturar número de iterações?
    return getvalue.(x),getobjectivevalue(P),0,((termination_status(P)) == 1 ? 0 : 1)
end

max_stars_repo_name: leonardosecchin/temp

max_issues_repo_name: leonardosecchin/temp

max_forks_repo_name: leonardosecchin/temp


text:
using BenchmarkTools, ResumableFunctions
using ResumableFunctions

const n = 93

const N = BigInt

function direct(a::N, b::N)
  b, a+b
end

function test_direct(n::Int)
  a, b = zero(N), one(N)
  for _ in 1:n-1
    a, b = direct(a, b)
  end
  a
end

@resumable function fibonacci_resumable(n::Int)
  a, b = zero(N), one(N)
  for _ in 1:n
    @yield a
    a, b = b, a + b
  end
end

@noinline function test_resumable(n::Int)
  a = 0
  for v in fibonacci_resumable(n)
    a = v
  end
  a
end

function fibonacci_channel(n::Int, ch::Channel)
  a, b = zero(N), one(N)
  for _ in 1:n
    put!(ch, a)
    a, b = b, a + b
  end
end

@noinline function test_channel(n::Int, csize::Int)
  fib_channel = Channel(c -> fibonacci_channel(n, c); ctype=N, csize=csize)
  a = 0
  for v in fib_channel
    a = v
  end
  a
end

struct Generator
  callee :: Task
  function Generator(f::Function, args...; kwargs...)
      ct = current_task()
      new(@task begin task_local_storage(:caller, ct); f(args...; kwargs...); yield(ct) end)
  end
end

function Base.iterate(gen::Generator, state=nothing)
  ret = yieldto(gen.callee)
  if ret === nothing return nothing end
  ret, nothing
end

function consume(gen::Generator, val=nothing)
  yieldto(gen.callee, val)
end

function produce(val=nothing)
  t = task_local_storage(:caller)
  yieldto(t, val)
end

function fibonacci_task(n::Int)
  a,b = zero(N), one(N)
  for _ in 1:n
      produce(a)
      a, b = b, a+b
  end
end

@noinline function test_task(n::Int)
  a = 0
  for v in Generator(fibonacci_task, n)
    a = v
  end
  a
end

function fibonacci_closure()
  a, b = zero(N), one(N)
  function()
    tmp = a
    a, b = b, a + b
    tmp
  end
end

@noinline function test_closure(n::Int)
  fib_closure = fibonacci_closure()
  a = 0
  for _ in 1:n
    a = fib_closure()
  end
  a
end

function fibonacci_closure_opt()
  a = Ref(zero(N))
  b = Ref(one(N))
  function()
    tmp = a[]
    a[], b[] = b[], a[] + b[]
    tmp
  end
end

@noinline function test_closure_opt(n::Int)
  fib_closure = fibonacci_closure_opt()
  a = 0
  for _ in 1:n
    a = fib_closure()
  end
  a
end

function fibonacci_closure_stm(n::Int)
  _state = Ref(0x00)
  a = Ref{N}()
  b = Ref{N}()
  _iterstate_1 = Ref{Int}()
  _iterator_1 = Ref{UnitRange{Int}}()

  function(_arg::Any=nothing)
    _state[] === 0x00 && @goto _STATE_0
    _state[] === 0x01 && @goto _STATE_1
    error("@resumable function has stopped!")
    @label _STATE_0
    _state[] = 0xff
    _arg isa Exception && throw(_arg)
    a[], b[] = zero(Int), one(Int)
    _iterator_1[] = 1:n
    _iteratornext_1 = iterate(_iterator_1[])
    while _iteratornext_1 !== nothing
      (_, _iterstate_1[]) = _iteratornext_1
      _state[] = 0x01
      return a[]
      @label _STATE_1
      _state[] = 0xff
      _arg isa Exception && throw(_arg)
      (a[], b[]) = (b[], a[] + b[])
      _iteratornext_1 = iterate(_iterator_1[], _iterstate_1[])
    end
  end
end

const FibClosure = typeof(fibonacci_closure_stm(n))

function Base.iterate(f::FibClosure, state=nothing)
  a = f()
  f._state[] === 0xff && return nothing
  a, nothing
end

@noinline function test_closure_stm(n::Int)
  fib_closure = fibonacci_closure_stm(n)
  a = 0
  for v in fibonacci_closure_stm(n)
    a = v
  end
  a
end

struct FibN
  n::Int
end

function Base.iterate(f::FibN, state::Tuple{N, N, Int}=(zero(N), one(N), 1))
  a, b, iters = state
  iters > f.n && return nothing
  a, (b, a + b, iters + 1)
end

@noinline function test_iteration_protocol(n::Int)
  a = 0
  for v in FibN(n)
    a = v
  end
  a
end

isinteractive() || begin
  println("Direct: ")
  @btime test_direct($n)
  @assert test_direct(n) == 7540113804746346429

  println("ResumableFunctions: ")
  @btime test_resumable($n)
  @assert test_resumable(n) == 7540113804746346429

  println("Channels csize=0: ")
  @btime test_channel($n, $0)
  @assert test_channel(n, 0) == 7540113804746346429

  println("Channels csize=1: ")
  @btime test_channel($n, $1)
  @assert test_channel(n, 1) == 7540113804746346429

  println("Channels csize=20: ")
  @btime test_channel($n, $20)
  @assert test_channel(n, 20) == 7540113804746346429

  println("Channels csize=100: ")
  @btime test_channel($n, $100)
  @assert test_channel(n, 100) == 7540113804746346429

  println("Task scheduling")
  @btime test_task($n)
  @assert test_task(n) == 7540113804746346429

  println("Closure: ")
  @btime test_closure($n)
  @assert test_closure(n) == 7540113804746346429

  println("Closure optimised: ")
  @btime test_closure_opt($n)
  @assert test_closure_opt(n) == 7540113804746346429

  println("Closure statemachine: ")
  @btime test_closure_stm($n)
  @assert test_closure_stm(n) == 7540113804746346429

  println("Iteration protocol: ")
  @btime test_iteration_protocol($n)
  @assert test_iteration_protocol(n) == 7540113804746346429
end

ext: jl


text:
using DataFrames

const datafile = "test.csv"
const nloops = 10

function load_pkgs()
    println("loading packages...")
    @everywhere using Blocks
    @everywhere using DataFrames
end

if nprocs() < 4
    addwrkrs = 4 - nprocs()
    println("adding $addwrkrs more processors...")
    addprocs(addwrkrs)
end
println("\tnprocs: $(nprocs())")
load_pkgs()

df = dreadtable(datafile, header=false)
colnames(df)
colnames!(df, ["c1","c2","c3","c4","c5","c6","c7","c8","c9","c10"])

sum_result = df+df
mul_result = 2*df
eq_result = (sum_result .== mul_result)
@assert all(eq_result)

df1 = dreadtable(open(datafile), 1000)
@assert nrow(df1) == nrow(df)
@assert ncol(df1) == ncol(df)
@assert isapprox(sum(matrix(colsums(df1))), sum(matrix(colsums(df))))


text:
print("Testing mac: ")

v1 = [-2im, 1-im]
v2 = [1, 2]
m = mac(v1,v2)

if m > 0
    println("OK")
end


text:
include("Robustness.jl")

index_C3a = 19
index_C5a = 15
index_C5 = 14
index_C3 = 8

alpha_array_final_50 = zeros(65)
alpha_array_final_90 = zeros(65)
alpha_array_final_99 = zeros(65)

# 50% reduction C3a and C5a -
alpha_array_C3a_C3_50 = main(index_C3a,[index_C3],[0.5])
alpha_array_C3a_C5_50 = main(index_C3a,[index_C5],[0.5])

alpha_array_C5a_C3_50 = main(index_C5a,[index_C3],[0.5])
alpha_array_C5a_C5_50 = main(index_C5a,[index_C5],[0.5])

alpha_array_C3a_C3_C5_50 = main(index_C3a,[index_C3,index_C5],[0.5,0.5])
alpha_array_C5a_C3_C5_50 = main(index_C5a,[index_C3,index_C5],[0.5,0.5])

alpha_array_final_50 = [alpha_array_final_50 alpha_array_C3a_C3_50 alpha_array_C3a_C5_50 alpha_array_C5a_C3_50 alpha_array_C5a_C5_50 alpha_array_C3a_C3_C5_50 alpha_array_C5a_C3_C5_50]

# 90% -
alpha_array_C3a_C3_90 = main(index_C3a,[index_C3],[0.1])
alpha_array_C3a_C5_90 = main(index_C3a,[index_C5],[0.1])

alpha_array_C5a_C3_90 = main(index_C5a,[index_C3],[0.1])
alpha_array_C5a_C5_90 = main(index_C5a,[index_C5],[0.1])

alpha_array_C3a_C3_C5_90 = main(index_C3a,[index_C3,index_C5],[0.1,0.1])
alpha_array_C5a_C3_C5_90 = main(index_C5a,[index_C3,index_C5],[0.1,0.1])

alpha_array_final_90 = [alpha_array_final_90 alpha_array_C3a_C3_90 alpha_array_C3a_C5_90 alpha_array_C5a_C3_90 alpha_array_C5a_C5_90 alpha_array_C3a_C3_C5_90 alpha_array_C5a_C3_C5_90]

# 99% -
alpha_array_C3a_C3_99 = main(index_C3a,[index_C3],[0.01])
alpha_array_C3a_C5_99 = main(index_C3a,[index_C5],[0.01])

alpha_array_C5a_C3_99 = main(index_C5a,[index_C3],[0.01])
alpha_array_C5a_C5_99 = main(index_C5a,[index_C5],[0.01])

alpha_array_C3a_C3_C5_99 = main(index_C3a,[index_C3,index_C5],[0.01,0.01])
alpha_array_C5a_C3_C5_99 = main(index_C5a,[index_C3,index_C5],[0.01,0.01])

alpha_array_final_99 = [alpha_array_final_99 alpha_array_C3a_C3_99 alpha_array_C3a_C5_99 alpha_array_C5a_C3_99 alpha_array_C5a_C5_99 alpha_array_C3a_C3_C5_99 alpha_array_C5a_C3_C5_99]

# Write the results files -
writedlm("./robustness/alpha_Z1_50.dat",alpha_array_final_50[:,2:end])
writedlm("./robustness/alpha_Z1_90.dat",alpha_array_final_90[:,2:end])
writedlm("./robustness/alpha_Z1_99.dat",alpha_array_final_99[:,2:end])


text:
using LinearAlgebra
using NamedDims
using NamedDims: matrix_prod_names, names, symmetric_names
using Test

@testset "+" begin
    nda = NamedDimsArray{(:a,)}(ones(3))

    @testset "standard case" begin
        @test +(nda) == ones(3)
        @test names(+(nda)) == (:a,)

        @test +(nda, nda) == 2ones(3)
        @test names(+(nda, nda)) == (:a,)

        @test +(nda, nda, nda) == 3ones(3)
        @test names(+(nda, nda, nda)) == (:a,)
    end

    @testset "partially named dims" begin
        ndx = NamedDimsArray{(:x, :_)}(ones(3, 5))
        ndy = NamedDimsArray{(:_, :y)}(ones(3, 5))

        lhs = ndx + ndy
        rhs = ndy + ndx
        @test names(lhs) == (:x, :y) == names(rhs)
        @test lhs == 2ones(3, 5) == rhs
    end

    @testset "Dimension disagreement" begin
        @test_throws DimensionMismatch +(
            NamedDimsArray{(:a, :b, :c, :d)}(zeros(3, 3, 3, 3)),
            NamedDimsArray{(:w, :x, :y, :z)}(ones(3, 3, 3, 3))
        )

        @test_throws DimensionMismatch +(
            NamedDimsArray{(:time,)}(zeros(3,)),
            NamedDimsArray{(:time, :value)}(ones(3, 3))
        )
    end

    @testset "Mixed array types" begin
        lhs_sum = +(
            NamedDimsArray{(:a, :b, :c, :d)}(zeros(3, 3, 3, 3)),
            ones(3, 3, 3, 3)
        )
        @test lhs_sum == ones(3, 3, 3, 3)
        @test names(lhs_sum) == (:a, :b, :c, :d)

        rhs_sum = +(
            zeros(3, 3, 3, 3),
            NamedDimsArray{(:w, :x, :y, :z)}(ones(3, 3, 3, 3))
        )
        @test rhs_sum == ones(3, 3, 3, 3)
        @test names(rhs_sum) == (:w, :x, :y, :z)

        casts = (NamedDimsArray{(:foo, :bar)}, identity)
        for (T1, T2, T3, T4) in Iterators.product(casts, casts, casts, casts)
            all(isequal(identity), (T1, T2, T3, T4)) && continue
            total = T1(ones(3, 6)) + T2(2ones(3, 6)) + T3(3ones(3, 6)) + T4(4ones(3, 6))
            @test total == 10ones(3, 6)
            @test names(total) == (:foo, :bar)
        end
    end
end

@testset "-" begin
    # This is actually covered by the tests for + above, since that uses the same code
    # just one extra as a sensability check
    nda = NamedDimsArray{(:a, :b)}(ones(3, 100))
    @test nda - nda == zeros(3, 100)
    @test names(nda - nda) == (:a, :b)
end

@testset "scalar product" begin
    nda = NamedDimsArray{(:a, :b, :c, :d, :e)}(ones(10, 20, 30, 40, 50))
    @test 10nda == 10ones(10, 20, 30, 40, 50)
    @test names(10nda) == (:a, :b, :c, :d, :e)
end

@testset "matmul" begin
    @testset "matrix_prod_names" begin
        @test matrix_prod_names((:foo, :bar), (:bar, :buzz)) == (:foo, :buzz)
        @test matrix_prod_names((:foo, :bar), (:_, :buzz)) == (:foo, :buzz)
        @test matrix_prod_names((:foo, :_), (:bar, :buzz)) == (:foo, :buzz)
        @test matrix_prod_names((:foo, :_), (:_, :buzz)) == (:foo, :buzz)
        @test_throws DimensionMismatch matrix_prod_names((:foo, :bar), (:nope, :buzz))

        @test matrix_prod_names((:foo,), (:bar, :buzz)) == (:foo, :buzz)
        @test matrix_prod_names((:foo,), (:_, :buzz)) == (:foo, :buzz)
        # No error case with name mismatch here, as a Vector has "virtual" wildcard second dimension

        @test matrix_prod_names((:foo, :bar), (:bar,)) == (:foo,)
        @test matrix_prod_names((:foo, :bar), (:_, )) == (:foo,)
        @test matrix_prod_names((:foo, :_), (:bar,)) == (:foo,)
        @test matrix_prod_names((:foo, :_), (:_,)) == (:foo,)
        @test_throws DimensionMismatch matrix_prod_names((:foo, :bar), (:nope,))
    end

    @testset "Matrix-Matrix" begin
        nda = NamedDimsArray{(:a, :b)}(ones(2, 3))
        ndb = NamedDimsArray{(:b, :c)}(ones(3, 2))

        @testset "standard case" begin
            @test nda * ndb == 3ones(2, 2)
            @test names(nda * ndb) == (:a, :c)

            @test ones(4, 3) * ndb == 3ones(4, 2)
            @test names(ones(4, 3) * ndb) == (:_, :c)

            @test nda * ones(3, 7) == 3ones(2, 7)
            @test names(nda * ones(3,7)) == (:a, :_)
        end

        @testset "Dimension disagreement" begin
            @test_throws DimensionMismatch ndb * nda
        end
    end

    @testset "Matrix-Vector" begin
        ndm = NamedDimsArray{(:a, :b)}(ones(1, 1))
        ndv = NamedDimsArray{(:b, )}(ones(1))

        @test ndm * ndv == ones(1)
        @test names(ndm * ndv) == (:a,)
    end

    @testset "Vector-Matrix" begin
        ndm = NamedDimsArray{(:a, :b)}(ones(1, 1))
        ndv = NamedDimsArray{(:a, )}(ones(1))

        @test ndv * ndm == ones(1, 1)
        @test names(ndv * ndm) == (:a, :b)
    end

    @testset "Vector-Vector" begin
        v = [1, 2, 3]
        ndv = NamedDimsArray{(:vec,)}(v)
        @test_throws MethodError ndv * ndv
        @test ndv' * ndv == 14
        @test ndv' * ndv == adjoint(ndv) * v == transpose(ndv) * v
        @test ndv' * ndv == adjoint(v) * ndv == transpose(v) * ndv
        @test ndv * ndv' == [1 2 3; 2 4 6; 3 6 9]

        ndv2 = NamedDimsArray{(:b,)}([3, 2, 1])
        @test_throws DimensionMismatch ndv' * ndv2
    end
end
@testset "allocations: matmul names" begin
    @test 0 == @allocated (() -> matrix_prod_names((:foo, :bar), (:bar,)))()
    @test 0 == @allocated (() -> symmetric_names((:foo, :bar), 1))()
end

@testset "Mutmul with special types" begin
    nda = NamedDimsArray{(:a, :b)}(ones(5,5))
    @testset "$T" for T in (Diagonal, Symmetric, Tridiagonal, SymTridiagonal, BitArray,)
        x = T(ones(5,5))
        @test names(x * nda) == (:_, :b)
        @test names(nda * x) == (:a, :_)
    end
end

@testset "inv" begin
    nda = NamedDimsArray{(:a, :b)}([1.0 2; 3 4])
    @test names(inv(nda)) == (:b, :a)
    @test nda * inv(nda) ≈ NamedDimsArray{(:a, :a)}([1.0 0; 0 1])
    @test inv(nda) * nda ≈ NamedDimsArray{(:b, :b)}([1.0 0; 0 1])
end

@testset "cov/cor" begin
    @testset "symmetric_names" begin
        @test symmetric_names((:a, :b), 1) == (:b, :b)
        @test symmetric_names((:a, :b), 2) == (:a, :a)
        @test symmetric_names((:a, :b), 5) == (:_, :_)

        @test_throws MethodError symmetric_names((:a, :b, :c), 2)
    end
    @testset "$f" for f in (cov, cor)
        @testset "matrix input, matrix result" begin
            A = rand(3, 5)
            nda = NamedDimsArray{(:a, :b)}(A)
            @test f(nda; dims=:a) == f(A, dims=1)
            @test names(f(nda; dims=:a)) == (:b, :b)
            @test names(f(nda, dims=:b)) == (:a, :a)
            # `Statistic.cov/cor(A, dims=p)` for `p > 2` is allowed but returns NaNs.
            @test names(f(nda, dims=3)) == (:_, :_)
        end
        @testset "vector input, scalar result" begin
            a = rand(4)
            nda = NamedDimsArray{(:a,)}(a)
            @test f(nda) isa Number
            @test f(nda) == f(a)
        end
        @testset "high dimensional input" begin
            @test_throws MethodError f(NamedDimsArray(rand(3, 4, 5), (:a, :b, :c)))
        end
    end
    @testset "cov corrected=$bool" for bool in (true, false)
        # test that kwargs get passed on correctly
        A = rand(2, 4)
        nda = NamedDimsArray{(:a, :b)}(A)
        @test cov(nda; corrected=bool) == cov(A; corrected=bool)
        @test cov(nda; corrected=bool, dims=:b)  == cov(A; corrected=bool, dims=2)
    end
end


text:
module CUPTI

using ..APIUtils

using ..CUDA
using ..CUDA: libcupti
using ..CUDA: CUuuid, CUcontext, CUstream, CUdevice, CUdevice_attribute,
              CUgraph, CUgraphNode, CUgraphNodeType, CUgraphExec, CUaccessPolicyWindow

using CEnum

using Memoize

# core library
include("libcupti_common.jl")
include("error.jl")
include("libcupti.jl")

include("wrappers.jl")

end


text:
using PyCall
using StatsBase

abstract type PyODmodel end

function StatsBase.fit!(model::PyODmodel, X::Array{T, 2}) where T<:Real
	# transposition since Python models are row major
    model.model.fit(Array(transpose(X)))
end

function StatsBase.predict(model::PyODmodel, X::Array{T, 2}) where T<:Real
	model.model.decision_function(Array(transpose(X)))
end

"""
	LODA(contamination=0.1, n_bins=10, n_random_cuts=100)

The LODA model.
"""
mutable struct LODA <: PyODmodel
	model

	function LODA(;kwargs...)
		py"""
		from pyod.models.loda import LODA

		def construct_loda(kwargs):
			return LODA(**kwargs)
		"""
		new(py"construct_loda"(kwargs))
	end
end

"""
	ABOD(n_neighbors=5, method='fast')

The ABOD model.
"""
mutable struct ABOD <: PyODmodel
	model

	function ABOD(;kwargs...)
		py"""
		from pyod.models.abod import ABOD

		def construct_abod(kwargs):
			return ABOD(**kwargs)
		"""
		new(py"construct_abod"(kwargs))
	end
end

"""
	HBOS(n_bins=10, alpha=0.1, tol=0.5)

The HBOS model.
"""
mutable struct HBOS <: PyODmodel
	model

	function HBOS(;kwargs...)
		py"""
		from pyod.models.hbos import HBOS

		def construct_hbos(kwargs):
			return HBOS(**kwargs)
		"""
		new(py"construct_hbos"(kwargs))
	end
end

"""
	MO_GAAL(k=10, stop_epochs=20, lr_d=0.01, lr_g=0.0001, decay=1e-06, momentum=0.9, contamination=0.1)

The MO_GAAL model. (GAN based)
"""
mutable struct MO_GAAL <: PyODmodel
	model

	function MO_GAAL(;kwargs...)
		py"""
		from pyod.models.mo_gaal import MO_GAAL

		def construct_mogaal(kwargs):
			return MO_GAAL(**kwargs)
		"""
		new(py"construct_mogaal"(kwargs))
	end
end


text:
isdefined(Base, :__precompile__) && __precompile__()

#__precompile__()

module JLib

using JSON, Compat, DataArrays, GLM, DataFrames, Distributions, NLopt, Showoff, StatsBase, DataStructures, StatsFuns, JuMP, MixedModels, NamedArrays, NLsolve
#using JSON, Requests, HttpParser

import DataArrays
#import DBAPI, JavaCall, JDBC
import StatsBase: coef, coeftable, df, deviance, fit!, fitted, loglikelihood, model_response, nobs, vcov
import Base: cond, std
import Distributions: Bernoulli, Binomial, Poisson, Gamma
import GLM: LogitLink, LogLink, InverseLink
import DataFrames: @~

export
       @~,
       lowercase,
       grep,
       grep1,
       lowercase!,
       loadCFG,
       df2dict,
       FS_singleLevel,
       genF,
       retypeBools!,
       poolit!, unpool,
       getColswithType,
       corDFD,
       cmparr,
       vif,
       vifDF, vifDF2,
       coefDF,
       checksingularity,
       xResiduals,
       gethostModelDistribution,
       gethostworkers,
       genF,
       raneffect,
       relevel,
       profileθ,
       readModels, expandM,  #saveModels,
       dfcor1, dfcor,
       read_dfx, save_dfx, read_modelsDict, save_modelsDict, read_dfd, save_dfd, read_dfr, save_dfr, read_orig,
       isnumeric, isfloat, isint, isString, isBool,
       isDict, isArray, isNum, isFloat, isInt, isSymbol,
       filler,
       ZDict, CIs_O, calcPValue_Opt     #, np_val,
       dict2json, dict_Sym, mergeDict,  savejson, save_cfg, read_cfg,
       blank

import Base: ==, *

include("common.jl")
include("featureSelection.jl")
include("CIs.jl")
#include("diagnostics.jl")
##include("misc.jl")
#include("modsel.jl")

""" x
# vif(m::DataFrameRegressionModel)  --  The variance inflation factors for the terms in a model.

function vif(m)
    v = vcov(m)
    assign = m.mm.assign
    nms = coefnames(m.mf)
    if (ipos = findfirst(nms, "(Intercept)")) > 0
        inds = deleteat!(collect(eachindex(nms)), ipos)
        v = view(v, inds, inds)
        assign = view(assign, inds)
    else
        warning("No intercept: vifs may not be sensible.")
    end
    terms = convert(Vector{Symbol}, m.mf.terms.terms)
    n_terms = length(terms)
    if n_terms < 2
        error("model contains fewer than 2 terms")
    end
    R = cov2cor(v)
    detR = det(Symmetric(R))
    result = NamedArray(Array(eltype(R), (length(terms), 3)))
    setnames!(result, string.(terms), 1)
    setnames!(result, ["GVIF", "Df", "GVIF^(1/(2*Df))"], 2)
    inds = collect(eachindex(terms))
    ainds = collect(eachindex(assign))
    for i in inds
        subs = find(x -> x == i, assign)
        notsubs = setdiff(ainds, subs)
        result[i, 1] = det(Symmetric(R[notsubs, notsubs])) / detR
        if (result[i, 2] = length(subs)) > 1
            result[i, 1] *= det(Symmetric(R[subs, subs]))
        end
    end
    if all(result[:, 2] .== 1)
        return result[:, 1:1]
    end
    result
end
"""

function vifDF(m::RegressionModel)
    sdf=coefDF(m)
    vdf=vif(m)
    vdf[:parameter] = convert(Array{Symbol}, vdf[:variable])
    join(sdf,vdf[[:parameter,:vif]], on = :parameter, kind=:outer)
end

#dfrm = glm( exposed_flag ~ 1 + Prd_2_Qty_POS, df_in  , Poisson(), LogLink() )

function vif(dfrm::RegressionModel)
    X = dfrm.mf.df[2:end]
    rhs = extractrhs(dfrm)
    result = DataFrame(variable=rhs.rhsarray, vif=0.0)
    i = 1
    for var in rhs.rhsarray
        lhs = parse(var)
        rhsnew = replace(rhs.rhsstring, "+"*var, "")
        rhsnew = rhsnew[2:end]
        rhsnew = parse(rhsnew)
        fnew = Formula(lhs, rhsnew)
        newfit = fit(LinearModel, fnew, X)
        r2 = rsquared(newfit)
        result[:vif][i] = 1 / (1-r2)
        i = i + 1
    end
    result
end

function rsquared(dfrm::RegressionModel)
    SStot = sum((dfrm.model.rr.y - mean(dfrm.model.rr.y)).^2)
    SSres = sum((dfrm.model.rr.y - dfrm.model.rr.mu).^2)
    return (1-(SSres/SStot))
end

type extractrhs
    rhsarray::Array
    rhsstring::AbstractString
end

# Extract the right hand side of a RegressionModel formula
function extractrhs(dfrm::RegressionModel)
    regressors = dfrm.mf.terms.terms
    rhsarray = similar(regressors, AbstractString)
    for i in 1:length(regressors)
        rhsarray[i] = string(regressors[i])
        rhsarray[i] = replace(rhsarray[i], " ", "")
    end
    rhsstring = ""
    for var in rhsarray
        rhsstring = rhsstring*"+"*var
    end
    if !dfrm.mf.terms.intercept
        rhsarray = [0, rhsarray]
        rhsstring = "0+"*rhsstring
    end
    return extractrhs(rhsarray, rhsstring)
end

end # module

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
# load packages
using Plots
#pyplot()
gr()

# load modules
# select a model later
include("HexaLatticeModel03a.jl")
include("SquareLatticeModel03a.jl")
# MCMC
include("IsingMCMC03a.jl")

###################################################################################################
# view functions
function drawEnergy(simE, trial, T, d=d, N=N)
    xs = [t for t in 0:1000:trial]
    ys = simE
    p = Plots.plot(xs, ys, label="Energy history with T=$(T)")
end

function drawSim(Sf, d=d, N=N)
    snap = Sf

    p=Plots.plot([],[],legend=false, size=(500,500))

    npoly = length(genPoly(1,d)[1])
    nbuff = (npoly+1)*N
    p0xs = fill(NaN, nbuff);p0ys = fill(NaN, nbuff)
    p1xs = fill(NaN, nbuff);p1ys = fill(NaN, nbuff)
    i1 = 1; i0 = 1
    @inbounds @simd for i in 1:N
        # bits
        (xs,ys) = genPoly(i, d)
        if (snap[i] == 1) # up-spin
            p1xs[i1:(i1+npoly-1)] .= xs
            p1ys[i1:(i1+npoly-1)] .= ys
            i1 += npoly+1
        else # down-spin
            p0xs[i0:(i0+npoly-1)] .= xs
            p0ys[i0:(i0+npoly-1)] .= ys
            i0 += npoly+1
        end
    end
    poly0 = Shape(p0xs, p0ys) # creating polygon
    p = Plots.plot!(poly0, fillcolor = plot_color(:blue)  ,legend=false)
    poly1 = Shape(p1xs, p1ys)
    p = Plots.plot!(poly1, fillcolor = plot_color(:yellow),legend=false)
end

###################################################################################################
# prepare
# setting for energy function

function genInitialJB(d=d)
    Jd = [Dict{Int,Int8}() for i in 1:d.N]
    Bd = zeros(Int8, d.N)
    (Jd, Bd)
end

function setJB!( Jd=Jd, Bd=Bd, d=d)
    J0 = 1
    B0 = 0
    for i in 1:d.N
        for j in getNeighbors(i,d)
            @inbounds Jd[i][j] = J0
        end
        @inbounds Bd[i] = B0
    end
end

###################################################################################################
# do simulation

# setting for model
d = SquareLatticeModel.ParamDict(128,128,0)
#d = SquareLatticeModel.ParamDict(64,64,0)
#d = SquareLatticeModel.ParamDict(32,32,0)
#d = SquareLatticeModel.ParamDict(8,8,0)
N = d.nrows * d.ncols
d.N = N

# main function
# mdlno: 1(square), 2(hexa)
# T: temperature e.g. T=2.26
function main(mdlno, T)
    # setting for simulation
    trial = 100000 # MCMC trial

    # gen functors as global
    global getNeighbors, genPoly
    if mdlno == 1
        (getNeighbors, genPoly)= SquareLatticeModel.genFunctor(d)
    elseif mdlno == 2
        (getNeighbors, genPoly)= HexaLatticeModel.genFunctor(d)
    end

    # generate Jd, Bd
    (Jd, Bd) = genInitialJB(d)
    setJB!(Jd, Bd)

    # simulation
    global E, MCMC
    (E, MCMC) = IsingMCMC.genFunctor(Jd, Bd)
    print("start MCMC...")
    @time (simE, Sf) = MCMC(T, N, trial);

    # save images
    # energy history
    print("start drawEnergy...")
    @time p = drawEnergy(simE, trial, T)
    Plots.savefig(p, "Ising_energy_M$(mdlno)_T$(T).png")

    # cell state
    print("start drawSim...")
    @time p = drawSim(Sf, d)
    Plots.savefig(p, "Ising_sim_M$(mdlno)_T$(T).png")
    print("done.")
end
###################################################################################################


text:
function compute_factor(Z::Array{Float64,2},
                        g::Array{Float64,1},
                        p::Array{Float64,1})::Array{Float64,1}
    return (Z ./(p'*Z))*g
end

function do_cover_iterations(Z::Array{Float64,2},
                             g::Array{Float64,1},
                             p::Array{Float64,1},
                             num_of_iter::Int64 = 10)::Array{Float64,1}
    for i in 1:num_of_iter   ## do num_of_iter iterations
        p = p.* compute_factor(Z,g,p)
    end
    return p
end

function cancel_zero_alternatives(Z::Array{Float64,2},
                                  p::Array{Float64,1},
                                  aTolSetChoiceProbToZero::Real=0)
    notToZero = findall(a->(!isapprox(a,0; atol = aTolSetChoiceProbToZero)), p)
    if length(notToZero)<length(p)   ## we'll set props close to zero to equal zero
        p=p[notToZero]
        p/=sum(p)
        Z=Z[notToZero,:]
    end
    return Z, p, notToZero
end

function run_cover_algorithm(Z::Array{Float64,2},
                         g::Array{Float64,1};
                         initialIterations::Int  = 10,
                         aTolSetChoiceProbToZero::Float64 = sqrt(eps()),
                         aTolerance4Stopping::Float64 = eps(),
                         maxIterations::Int = 9999)

    #initialize
    message = "Run Cover Algorithm"
    notToZeroLists = []
    localZ = copy(Z)
    m = size(Z)[1]
    p = ones(m)./m    ## initialize the uncondition probs column vector, i.e. arrays of dimension n x 1
    iterationsDone = 0

    for coverCall in 1:fld(maxIterations,initialIterations)
        p = do_cover_iterations(localZ, g, p, initialIterations)
        iterationsDone += initialIterations

        localZ, p, notToZero = cancel_zero_alternatives(localZ, p, aTolSetChoiceProbToZero)
        push!(notToZeroLists, notToZero)

        ## check necessary conditions (cover alg doesn't improve anymore beyond precision)
        all(  isapprox.( compute_factor(localZ,g,p), ones(length(p)); atol = aTolerance4Stopping )  ) && break
    end

    ### next, we put back the zeros. to this end, we infer the indices that survived
    idxs = collect(1:m)
    while length(notToZeroLists)>0
        idxs=idxs[notToZeroLists[1]]
        popfirst!(notToZeroLists)
    end
    p_withzeros = zeros(m)
    p_withzeros[idxs]=p
    return p_withzeros, ((log.(p'*localZ))*g)[1] , iterationsDone, message
    # optimal probs, value, number of iterations, message
end


text:
include("solver_preamble.jl")
import CDCS
# Iterations:
# dsos_concave_then_convex_cubic : > 2000, < 3000
factory = with_optimizer(CDCS.Optimizer, verbose=0, maxIter=3000)
config = MOI.Test.TestConfig(atol=1e-3, rtol=1e-3, query=false)
@testset "Linear" begin
    Tests.linear_test(factory, config)
end
@testset "SOC" begin
    Tests.soc_test(factory, config)
end
@testset "SDP" begin
    Tests.sd_test(factory, config)
end


text:
# SKCE

# predicted normal distributions with squared exponential kernel for the targets
function CalibrationErrors.unsafe_skce_eval_targets(
    kernel::SqExponentialKernel{Euclidean}, p::Normal, y::Real, p̃::Normal, ỹ::Real
)
    # extract parameters
    μ = p.μ
    σ = p.σ
    μ̃ = p̃.μ
    σ̃ = p̃.σ

    # compute scaling factors
    # TODO: use `hypot`?
    sqσ = σ^2
    sqσ̃ = σ̃^2
    α = inv(sqrt(1 + sqσ))
    β = inv(sqrt(1 + sqσ̃))
    γ = inv(sqrt(1 + sqσ + sqσ̃))

    return kernel(y, ỹ) - α * kernel(α * μ, α * ỹ) - β * kernel(β * y, β * μ̃) +
           γ * kernel(γ * μ, γ * μ̃)
end

# UCME

function CalibrationErrors.unsafe_ucme_eval_targets(
    kernel::SqExponentialKernel{Euclidean}, p::Normal, y::Real, ::Normal, testy::Real
)
    # compute scaling factor
    # TODO: use `hypot`?
    α = inv(sqrt(1 + p.σ^2))

    return kernel(y, testy) - α * kernel(α * p.μ, α * testy)
end

# kernels with input transformations
# TODO: scale upfront?

function CalibrationErrors.unsafe_skce_eval_targets(
    kernel::TransformedKernel{
        SqExponentialKernel{Euclidean},<:Union{ScaleTransform,ARDTransform}
    },
    p::Normal,
    y::Real,
    p̃::Normal,
    ỹ::Real,
)
    # obtain the transform
    t = kernel.transform

    return CalibrationErrors.unsafe_skce_eval_targets(
        SqExponentialKernel(), apply(t, p), t(y), apply(t, p̃), t(ỹ)
    )
end

function CalibrationErrors.unsafe_ucme_eval_targets(
    kernel::TransformedKernel{
        SqExponentialKernel{Euclidean},<:Union{ScaleTransform,ARDTransform}
    },
    p::Normal,
    y::Real,
    testp::Normal,
    testy::Real,
)
    # obtain the transform
    t = kernel.transform

    # `testp` is irrelevant for the evaluation and therefore not transformed
    return CalibrationErrors.unsafe_ucme_eval_targets(
        SqExponentialKernel(), apply(t, p), t(y), testp, t(testy)
    )
end

# utilities

# internal `apply` avoids type piracy of transforms
apply(t::Union{ScaleTransform,ARDTransform}, d::Normal) = Normal(t(d.μ), t(d.σ))


text:
using BinaryBuilder

src_version = v"3.0.0"  # also change in raw script string

# Collection of sources required to build GDAL
sources = [
    "https://download.osgeo.org/gdal/$src_version/gdal-$src_version.tar.xz" =>
    "ad316fa052d94d9606e90b20a514b92b2dd64e3142dfdbd8f10981a5fcd5c43e",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd gdal-3.0.0/

if [[ ${target} == *w64-mingw32* ]]; then
    # Symlink libproj for Windows, else configure couldn't find it
    # TODO fix in PROJBuilder or in GDAL configure?
    ln -s $prefix/lib/libproj_6_1.dll.a $prefix/lib/libproj.dll.a
    # Also symlink the library folder of mingw, so libstdc++ and
    # others can be found. The path is ignored by BinaryBuilder.
    mkdir $prefix/$target
    ln -s /opt/$target/$target/lib $prefix/$target/lib
elif [[ ${target} == *freebsd* ]]; then
    # FreeBSD's default Clang ran into issues with configure
    # which seemed to think it was the GNU compiler. Therefore,
    # let's just use the GNU compiler instead.
    CC=gcc
    CXX=g++
fi

# Show options in the log
./configure --help
./configure --prefix=$prefix --host=$target \
    --with-geos=$prefix/bin/geos-config \
    --with-proj=$prefix \
    --with-libz=$prefix \
    --with-sqlite3=$prefix \
    --with-curl=$prefix/bin/curl-config \
    --with-python=no \
    --enable-shared \
    --disable-static \
    "CC=$CC" \
    "CXX=$CXX"

make -j${nproc}
make install

# strip shared libraries to reduce filesize
if [[ ${target} == *w64-mingw32* ]]; then
    strip $prefix/bin/*.dll
elif [[ ${target} == *apple-darwin* ]]; then
    :
else
    strip $prefix/lib/*.so
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
platforms = expand_gcc_versions(platforms)

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libgdal", :libgdal),
    ExecutableProduct(prefix, "gdal_contour", :gdal_contour_path),
    ExecutableProduct(prefix, "gdal_grid", :gdal_grid_path),
    ExecutableProduct(prefix, "gdal_rasterize", :gdal_rasterize_path),
    ExecutableProduct(prefix, "gdal_translate", :gdal_translate_path),
    ExecutableProduct(prefix, "gdaladdo", :gdaladdo_path),
    ExecutableProduct(prefix, "gdalbuildvrt", :gdalbuildvrt_path),
    ExecutableProduct(prefix, "gdaldem", :gdaldem_path),
    ExecutableProduct(prefix, "gdalinfo", :gdalinfo_path),
    ExecutableProduct(prefix, "gdallocationinfo", :gdallocationinfo_path),
    ExecutableProduct(prefix, "gdalmanage", :gdalmanage_path),
    ExecutableProduct(prefix, "gdalsrsinfo", :gdalsrsinfo_path),
    ExecutableProduct(prefix, "gdaltindex", :gdaltindex_path),
    ExecutableProduct(prefix, "gdaltransform", :gdaltransform_path),
    ExecutableProduct(prefix, "gdalwarp", :gdalwarp_path),
    ExecutableProduct(prefix, "nearblack", :nearblack_path),
    ExecutableProduct(prefix, "ogr2ogr", :ogr2ogr_path),
    ExecutableProduct(prefix, "ogrinfo", :ogrinfo_path),
    ExecutableProduct(prefix, "ogrlineref", :ogrlineref_path),
    ExecutableProduct(prefix, "ogrtindex", :ogrtindex_path)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaGeo/GEOSBuilder/releases/download/v3.7.2-0/build_GEOS.v3.7.2.jl",
    "https://github.com/JuliaGeo/PROJBuilder/releases/download/v6.1.0-1/build_PROJ.v6.1.0.jl",
    "https://github.com/bicycle1885/ZlibBuilder/releases/download/v1.0.4/build_Zlib.v1.2.11.jl",
    "https://github.com/JuliaDatabases/SQLiteBuilder/releases/download/v0.10.0/build_SQLite.v3.28.0.jl",
    "https://github.com/JuliaWeb/LibCURLBuilder/releases/download/v0.5.1/build_LibCURL.v7.64.1.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "GDAL", src_version, sources, script, platforms, products, dependencies)


text:
#!/usr/bin/env julia
#
# @license Apache-2.0
#
# Copyright (c) 2017 The Stdlib Authors.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#    http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

# Lints a list of Julia files.
#
# The script is called with one or more arguments, where each argument is a
# filepath. Note that each provided filepath is resolved relative to the current
# working directory of the calling process.

import Base.Test
import Lint

# Lint each file separately...
@Base.Test.testset "Lint tests" begin
	for i in 1:length( ARGS )
		@Base.Test.test isempty( Lint.lintfile( ARGS[ i ] ) );
	end
end

quit();

max_stars_repo_name: andbmme/stdlib

max_issues_repo_name: andbmme/stdlib

max_forks_repo_name: andbmme/stdlib


text:
################################################################################
#
#  Integral norm equations
#
################################################################################

# Returns the solution of the norm equation and the number of the order, in which
# the solution was found.
function norm_equation(orders::Vector{<: AlgAssRelOrd}, a::NfAbsOrdElem)
  A = algebra(orders[1])
  if iszero(a)
    return A(), 1
  end
  d, i = norm_equation_fac_elem(orders, a)
  return evaluate(d), i
end

function norm_equation_fac_elem(orders::Vector{<: AlgAssRelOrd}, a::NfAbsOrdElem)
  A = algebra(orders[1])
  @assert !iszero(a) # We cannot represent 0 as a FacElem
  NC = NormCache(A, orders, a)
  while true
    # We test the orders alternately until we find a solution in one of them
    for i = 1:length(orders)
      b, d = _norm_equation_relative(NC, i)
      if b
        return d, i
      end
    end
  end
end

# This tests max_num_fields fields. If max_num_fields == 0, we try until we
# find a solution (although at some point NC.field_oracle does not return any
# new fields).
function _norm_equation_relative(NC::NormCache, order_num::Int; max_num_fields::Int = 10)
  A = NC.algebra
  n = NC.n
  primes = NC.primes
  vals = NC.valuations
  full_set = Set(1:length(primes))
  Ok = parent(NC.a)
  mUk = NC.UktoOk
  Uk = domain(mUk)
  partial_solutions = NC.partial_solutions[order_num]
  solutions_mod_units = NC.solutions_mod_units[order_num]
  GtoUk = NC.GtoUk[order_num]
  G = domain(GtoUk)
  fields_in_product = NC.fields_in_product[order_num]

  m = 0
  while max_num_fields == 0 || m < max_num_fields
    m += 1
    LtoA = get_new_field(NC.field_oracle, order_num)
    L = domain(LtoA)
    OL = maximal_order(L)

    K, KtoL = simplified_absolute_field(L)
    ktoK = restrict(inv(KtoL), base_field(L))
    OK = maximal_order_via_relative(K, KtoL)

    good_primes = __neq_find_good_primes(NC, OL)

    local vals2::Vector{Int}

    if !isempty(good_primes)
      if degree(L) != n
        n2 = degree(L)//n
        vals2 = [ Int(n2*vals[i]) for i = 1:length(vals) ]
      else
        vals2 = vals
      end

      remaining_primes = Set{Int}()
      # First search for solutions for single primes
      for p in good_primes
        sols = __neq_sunit(ktoK, eltype(primes)[ primes[p] ], Int[ vals2[p] ])
        if !isempty(sols)
          # We better evaluate s here: Some factors may not be integral which can
          # become a problem since the multiplication in A is non-commutative.
          # Also, multiplications and computing inverses in K are cheap.
          add_partial_solutions!(NC, order_num, Set(p), [ NC.fac_elem_mon(LtoA(KtoL(evaluate(s::FacElem{nf_elem, AnticNumberField})))) for s in sols ])
        else
          push!(remaining_primes, p)
        end
      end

      if !isempty(remaining_primes)
        # Now the primes together, for which have not found a solution yet
        sols = __neq_sunit(ktoK, eltype(primes)[ primes[i] for i in remaining_primes ], Int[ vals2[i] for i in remaining_primes ])
        if !isempty(sols)
          add_partial_solutions!(NC, order_num, remaining_primes, [ NC.fac_elem_mon(LtoA(KtoL(evaluate(ss::FacElem{nf_elem, AnticNumberField})))) for ss in sols ])
        elseif length(remaining_primes) < length(good_primes)
          # If this also failed, we test all primes together
          sols = __neq_sunit(ktoK, [ primes[i] for i in good_primes ], [ vals2[i] for i in good_primes ])
          if !isempty(sols)
            add_partial_solutions!(NC, order_num, good_primes, [ NC.fac_elem_mon(LtoA(KtoL(evaluate(ss::FacElem{nf_elem, AnticNumberField})))) for ss in sols ])
          end
        end
      end
    end

    if !NC.GtoUk_surjective[order_num] && degree(L) == n
      UK, mUK = unit_group(OK)
      N = hom(UK, Uk, [ mUk\(norm(ktoK, K(mUK(g)))) for g in gens(UK) ])
      if ngens(G) == 0
        G = UK
        GtoUk = N
        push!(fields_in_product, (LtoA, KtoL))
        if is_surjective(GtoUk)
          NC.GtoUk_surjective[order_num] = true
        end
      else
        # Check whether this group brings anything new
        skip = true
        for g in gens(UK)
          b = haspreimage(GtoUk, N(g))[1]
          if !b
            skip = false
            break
          end
        end

        if !skip
          push!(fields_in_product, (LtoA, KtoL))

          G, pi = direct_product(G, UK, task = :prod)::Tuple{GrpAbFinGen, Tuple{GrpAbFinGenMap, GrpAbFinGenMap}}
          GtoUk = hom(gens(G), [ GtoUk(pi[1](g)) + N(pi[2](g)) for g in gens(G) ])
          if is_surjective(GtoUk)
            NC.GtoUk_surjective[order_num] = true
          end
        end
      end
      NC.GtoUk[order_num] = GtoUk
    end

    if !haskey(partial_solutions, full_set)
      continue
    end

    sols = partial_solutions[full_set]
    for s in sols
      u = get!(solutions_mod_units, s) do
        return mUk\divexact(NC.a, Ok(normred(s)))
      end
      b, v = haspreimage(GtoUk, u)
      if !b
        continue
      end
      x = __neq_lift_unit(NC, order_num, v)
      return true, s*x
    end
  end
  return false, NC.fac_elem_mon()
end

################################################################################
#
#  Integral norm equations - valuations only
#
################################################################################

# A special case needed for principal_generator_eichler.
# Finds a in O such that v_{p_i}(normred(a)) == v_i where p_i = primes[i] and
# v_i = valuations[i] and such that v_q(normred(a)) == 0 for all other primes q,
# assuming that such an element exists.
function _norm_equation_valuations_only(O::AlgAssRelOrd, primes::Vector{<: NfAbsOrdIdl}, valuations::Vector{Int})
  @assert !isempty(primes)
  @assert length(primes) == length(valuations)
  A = algebra(O)
  _K = base_ring(A)
  Ok = order(primes[1])
  NC = NormCache(A, typeof(O)[ O ], Ok, primes, valuations)
  n = NC.n
  full_set = Set(1:length(primes))
  partial_solutions = NC.partial_solutions[1]

  while true
    LtoA = get_new_field(NC.field_oracle, 1, no_restriction = true)
    L = domain(LtoA)::_ext_type(elem_type(_K))
    OL = maximal_order(L)

    K, KtoL = simplified_absolute_field(L)
    ktoK = restrict(inv(KtoL), base_field(L))
    OK = maximal_order_via_relative(K, KtoL)

    good_primes = __neq_find_good_primes(NC, OL)

    if !isempty(good_primes)
      cache = Vector{Any}(undef, 3) # Used in __neq_find_sol_in_order
      if degree(L) != n
        n2 = degree(L)//n
        vals2 = Int[ Int(n2*valuations[i]) for i = 1:length(valuations) ]
      else
        vals2 = valuations
      end

      remaining_primes = Set{Int}()
      # First search for solutions for single primes
      for p in good_primes
        if haskey(partial_solutions, Set(p))
          continue
        end
        b, sss::elem_type(A) = __neq_find_sol_in_order(O, LtoA, KtoL, ktoK, eltype(primes)[ primes[p] ], eltype(vals2)[ vals2[p] ], cache)
        if b
          add_partial_solutions!(NC, 1, Set(p), [ NC.fac_elem_mon(sss) ])
        else
          push!(remaining_primes, p)
        end
      end

      if !isempty(remaining_primes) && !haskey(partial_solutions, remaining_primes)
        # Now the primes together, for which have not found a solution yet
        b, s = __neq_find_sol_in_order(O, LtoA, KtoL, ktoK, eltype(primes)[ primes[i] for i in remaining_primes ], [ vals2[i] for i in remaining_primes ], cache)
        if b
          add_partial_solutions!(NC, 1, remaining_primes, [ NC.fac_elem_mon(s) ])
        elseif length(remaining_primes) < length(good_primes)
          # If this also failed, we test all primes together
          b, s = __neq_find_sol_in_order(O, LtoA, KtoL, ktoK, eltype(primes)[ primes[i] for i in good_primes ], [ vals2[i] for i in good_primes ], cache)
          if b
            add_partial_solutions!(NC, 1, good_primes, [ NC.fac_elem_mon(s) ])
          end
        end
      end
    end

    if haskey(partial_solutions, full_set)
      s = first(partial_solutions[full_set])
      return s
    end
  end
end

# Finds any a \in O \cap L such that v_p(nr(a)) = vals[i] for p = primes[i].
function __neq_find_sol_in_order(O::AlgAssRelOrd, LtoA::NfRelToAbsAlgAssMor, KtoL::NfToNfRel, ktoK::NfToNfMor, primes_in_k::Vector{<: NfAbsOrdIdl}, vals::Vector{Int}, cache::Vector{Any})
  A = algebra(O)
  sols = __neq_sunit(ktoK, primes_in_k, vals)
  if isempty(sols)
    return false, A()
  end

  K = codomain(ktoK)
  L = domain(LtoA)
  sols_eval = Vector{elem_type(K)}()
  for s in sols
    s_eval = evaluate(s)
    t = LtoA(KtoL(s_eval))
    if t in O
      return true, t
    end
    push!(sols_eval, s_eval)
  end

  # None of the found solutions lies in O.
  # We assume that equation_order(L) \subseteq O and try to adjust one of the
  # elements of sols by a unit.
  OK = maximal_order(K)
  UK, mUK = unit_group(OK)
  if !isassigned(cache, 1)
    cache[1] = Order(K, [ KtoL\b for b in absolute_basis(equation_order(L)) ], check = false, isbasis = true)
  end
  OE = cache[1]
  if !isassigned(cache, 2)
    cache[2] = OO_mod_F_mod_O_mod_F(OE)
  end
  G, GtoQ, OKtoQ = cache[2]
  sols2 = Vector{elem_type(K)}()
  for s in s_eval
    sinQ = OKtoQ(OK(s))
    if !is_invertible(sinQ)[1]
      push!(s, sols2)
    end
    # s is coprime to the conductor

    g = GtoQ\sinQ
    h = hom(UK, G, [ GtoQ\(OKtoQ(mUK(UK[i]))) for i = 1:ngens(UK) ])
    b, u = haspreimage(h, g)
    if b
      Q, toQ = quo(UK, kernel(h)[1])
      u = toQ\(toQ(u)) # Reduce the coefficient size (hopefully)
      s = s*inv(elem_in_nf(mUK(u), copy = false))
      return true, LtoA(KtoL(s))
    end
  end

  if !isassigned(cache, 3)
    UE, mUE = unit_group(OE)
    UEinUK = [ mUK\(OK(elem_in_nf(mUE(UE[i]), copy = false))) for i = 1:ngens(UE) ]
    cache[3] = quo(UK, UEinUK)
  end
  Q, toQ = cache[3]
  for (i, g) in enumerate(Q)
    u = mUK(toQ\g)
    for s in sols2
      su = LtoA(KtoL(s*u))
      if su in O
        return true, su
      end
    end
  end
  return false, A()
end

################################################################################
#
#  Some helpers
#
################################################################################

# Let S be the set of primes lying over primes_in_k.
# This finds a set of representatives of the S-units of K modulo units of O_K,
# whose norm has the valuations vals at the primes primes_in_k.
function __neq_sunit(ktoK::NfToNfMor, primes_in_k::Vector{<: NfAbsOrdIdl}, vals::Vector{Int})
  K = codomain(ktoK)
  OK = maximal_order(K)
  primes_in_K = Vector{ideal_type(OK)}()
  for p in primes_in_k
    p1 = minimum(p)
    p2 = uniformizer(p)
    @assert p == p1*order(p) + p2*order(p)
    P = p1*OK + OK(ktoK(elem_in_nf(p2, copy = false)))*OK
    append!(primes_in_K, collect(keys(factor(P))))
  end
  class_group(K)
  SK, mSK = sunit_mod_units_group_fac_elem(primes_in_K)
  M = zero_matrix(FlintZZ, length(primes_in_k), ngens(SK))
  N = zero_matrix(FlintZZ, length(primes_in_K), ngens(SK))
  b = zero_matrix(FlintZZ, length(primes_in_k), 1)
  for c = 1:ngens(SK)
    g = mSK(SK[c])
    ng = norm(ktoK, g)
    for r = 1:length(primes_in_k)
      M[r, c] = valuation(ng, primes_in_k[r])
      b[r, 1] = vals[r]
    end
    for r = 1:length(primes_in_K)
      N[r, c] = valuation(g, primes_in_K[r])
    end
  end

  sols = try
    solve_mixed(M, b, N)
  catch e
    if e isa UndefVarError
      error("This function requires Polymake.jl (original error message follows)")
    else
      rethrow(e)
    end
  end
  return elem_type(codomain(mSK))[ mSK(SK(sols[i, :])) for i in 1:nrows(sols) ]
end

function __neq_lift_unit(NC::NormCache, order_num::Int, g::GrpAbFinGenElem)
  if isone(g)
    return one(NC.algebra)
  end

  fields_in_product = NC.fields_in_product[order_num]
  G = domain(NC.GtoUk[order_num])
  if length(fields_in_product) == 1
    LtoA, KtoL = fields_in_product[1]
    K = domain(KtoL)
    mUK = unit_group(maximal_order(K))[2]
    return LtoA(KtoL(elem_in_nf(mUK(g), copy = false)))
  end

  GtoH = flat(G)
  H = codomain(GtoH)
  h = GtoH(g)
  x = NC.fac_elem_mon()
  for j = 1:length(fields_in_product)
    LtoA, KtoL = fields_in_product[j]
    pi = canonical_projection(H, j)
    h2 = pi(h)
    K = domain(KtoL)
    mUK = unit_group(maximal_order(K))[2]
    x *= LtoA(KtoL(elem_in_nf(mUK(h2), copy = false)))
  end
  return x
end

function __neq_find_good_primes(NC::NormCache, OL::NfRelOrd)
  n = NC.n
  m = degree(nf(OL))
  mn = m//n
  good_primes = Set{Int}()
  for j = 1:length(NC.primes)
    v = NC.valuations[j]*mn
    if !isone(denominator(v))
      continue
    end
    p = NC.primes[j]
    pdec = prime_decomposition(OL, p, compute_uniformizer = false)
    for (q, e) in pdec
      if q.splitting_type[2] <= v
        push!(good_primes, j)
        break
      end
    end
  end
  return good_primes
end

################################################################################
#
#  Functions for FieldOracle and NormCache
#
################################################################################

# Let K = base_ring(A). Then this returns the field K(x) and a map to A.
function _as_subfield(A::AbsAlgAss{T}, x::AbsAlgAssElem{T}) where { T <: Union{ fmpq, nf_elem, NfRelElem } }
  return _as_subfield(A, x, minpoly(x))
end

function _as_subfield(A::AbsAlgAss{fmpq}, x::AbsAlgAssElem{fmpq}, f::PolyElem{fmpq})
  s = one(A)
  M = zero_matrix(FlintQQ, degree(f), dim(A))
  elem_to_mat_row!(M, 1, s)
  for i = 2:degree(f)
    s = mul!(s, s, x)
    elem_to_mat_row!(M, i, s)
  end
  K = number_field(f)[1]
  return K, NfAbsToAbsAlgAssMor(K, A, M)
end

function _as_subfield(A::AbsAlgAss{T}, x::AbsAlgAssElem{T}, f::PolyElem{T}) where { T <: Union{ nf_elem, NfRelElem } }
  s = one(A)
  M = zero_matrix(base_ring(A), degree(f), dim(A))
  elem_to_mat_row!(M, 1, s)
  for i = 2:degree(f)
    s = mul!(s, s, x)
    elem_to_mat_row!(M, i, s)
  end
  K = number_field(f)[1]
  return K, NfRelToAbsAlgAssMor(K, A, M)
end

# The matrices in modules should generate full lattices in A.
# This returns the numbers of the modules of which LtoA(O) is a submodule.
function _issubmodule(modules::Vector{<: PMat}, O::NfRelOrd, LtoA::NfRelToAbsAlgAssMor)
  L = domain(LtoA)
  A = codomain(LtoA)
  B = absolute_basis(O)
  M = zero_matrix(base_ring(A), length(B), dim(A))
  for i = 1:length(B)
    elem_to_mat_row!(M, i, LtoA(B[i]))
  end
  PM = pseudo_hnf_kb(PseudoMatrix(M), :lowerleft)
  for i = 1:nrows(PM)
    if !is_zero_row(PM.matrix, i)
      PM = sub(PM, i:nrows(PM), 1:ncols(PM))
      break
    end
  end
  result = Vector{Int}()
  for i = 1:length(modules)
    if _spans_subset_of_pseudohnf(PM, modules[i], :lowerleft)
      push!(result, i)
    end
  end
  return result
end

function _issubmodule(modules::Vector{<: AlgAssAbsOrd}, O::NfAbsOrd, LtoA::NfAbsToAbsAlgAssMor)
  L = domain(LtoA)
  A = codomain(LtoA)
  B = basis(O)
  result = Vector{Int}()
  for i = 1:length(modules)
    issub = true
    for b in B
      if !(LtoA(elem_in_nf(b, copy = false)) in modules[i])
        issub = false
        break
      end
    end
    if issub
      push!(result, i)
    end
  end
  return result
end

_issubmodule(FO::FieldOracle, OL::NfRelOrd, LtoA::NfRelToAbsAlgAssMor) = _issubmodule(FO.hnf_basis_pmats, OL, LtoA)
_issubmodule(FO::FieldOracle, OL::NfAbsOrd, LtoA::NfAbsToAbsAlgAssMor) = _issubmodule(FO.maximal_orders, OL, LtoA)

# Returns a LLL-reduced basis of O
function small_elements(O::AlgAssRelOrd)
  A = algebra(O)
  K = base_ring(A)
  B, BtoA = restrict_scalars(A, FlintQQ)
  pbO = pseudo_basis(O, copy = false)
  M = zero_matrix(FlintQQ, dim(B), dim(B))
  for i = 1:degree(O)
    for j = 1:degree(K)
      t = basis(pbO[i][2])[j]*pbO[i][1]
      elem_to_mat_row!(M, (i - 1)*degree(K) + j, BtoA\(t))
    end
  end
  MM = FakeFmpqMat(M)
  L = lll(MM.num)
  res = Vector{elem_type(A)}()
  for i = 1:dim(B)
    t = elem_from_mat_row(B, L, i, MM.den)
    push!(res, BtoA(t))
  end
  return res
end

function small_elements(O::AlgAssAbsOrd)
  A = algebra(O)
  M = basis_matrix(O, copy = false)
  L = lll(M.num)
  res = Vector{elem_type(A)}()
  for i = 1:dim(A)
    t = elem_from_mat_row(A, L, i, M.den)
    push!(res, t)
  end
  return res
end

# Adds a field for order number i (and possibly for other orders too)
function add_field(FO::FieldOracle, i::Int; no_restriction::Bool = false)
  A = FO.algebra
  function _add_field(x::AbsAlgAssElem)
    f = minpoly(x)
    L, LtoA = _as_subfield(A, x, f)

    if no_restriction
      # The users wants no restriction on the fields. Let's hope he/she knows
      # what he/she is doing.
      for j = 1:length(FO.fields)
        push!(FO.fields[j], LtoA)
      end
      return true
    end

    # Find the orders in which maximal_order(L) lies
    good_orders = _issubmodule(FO, maximal_order(L), LtoA)
    for o in good_orders
      push!(FO.fields[o], LtoA)
    end
    if i in good_orders
      # We found one in the order we wanted
      return true
    end
    return false
  end

  function _increase_counter(k::Int)
    j = FO.last_generated_field[k]
    if j != length(FO.small_elements) - FO.rounds + k
      FO.last_generated_field[k] += 1
      return true
    end

    if k == 1
      return false
    end

    b = _increase_counter(k - 1)
    if !b
      return false
    end
    FO.last_generated_field[k] = FO.last_generated_field[k - 1] + 1
    return true
  end

  while true
    if FO.phase == 1
      if !_increase_counter(FO.rounds)
        FO.phase = 2
        continue
      end
      x = sum([ j != 0 ? FO.small_elements[j] : zero(A) for j in FO.last_generated_field ])
    else
      coeffs = rand(-FO.rand_coeff_bound:FO.rand_coeff_bound, length(FO.small_elements))
      x = sum( coeffs[i]*FO.small_elements[i] for i = 1:length(FO.small_elements) )
    end
    if !is_integral(x)
      continue
    end
    if !is_irreducible(minpoly(x))
      continue
    end
    if _add_field(x)
      return nothing
    end
  end
end

# Returns a "new" field for order i
function get_new_field(FO::FieldOracle, i::Int; no_restriction::Bool = false)
  if isempty(FO.fields[i])
    add_field(FO, i, no_restriction = no_restriction)
  end

  return pop!(FO.fields[i])
end

function add_partial_solutions!(NC::NormCache{T1, T2, T3}, order_num::Int, S::Set{Int}, sols::Vector{FacElem{T3, T1}}) where { T1, T2, T3 }
  partial_solutions = NC.partial_solutions[order_num]
  @assert haskey(partial_solutions, Set{Int}())
  for T in keys(partial_solutions)
    if !isempty(intersect(S, T))
      continue
    end

    ST = union(S, T)
    sols2 = Vector{FacElem{T3, T1}}()
    for b in sols
      for c in partial_solutions[T]
        bc = b*c
        push!(sols2, bc)
      end
    end
    if haskey(partial_solutions, ST)
      append!(partial_solutions[ST], sols2)
    else
      partial_solutions[ST] = sols2
    end
  end
  return nothing
end


text:
"""
A getter function for objects of type ``Atom``.
It also provides information not physically stored in the ``Atom`` object
by standardized by the design of ``Atom``.

Arguments
----------
obj:Atom
    a Atom object

field:Symbol
    the requested field name
"""
function obtain{T<:AbstractMolecularContainer}(obj::T, field::Symbol)
    if field in fieldnames(obj)
        return getfield(obj, field)
    else
        return nothing
    end
end


text:
export logjac_forwarddiff, value_and_logjac_forwarddiff, CustomTransform

"""
$(SIGNATURES)

Calculate the log Jacobian determinant of `f` at `x` using `ForwardDiff.

# Note

`f` should be a bijection, mapping from vectors of real numbers to vectors of
equal length.

When `handleNaN = true` (the default), NaN log Jacobians are converted to -Inf.
"""
function logjac_forwarddiff(f, x; handleNaN = true,
                            chunk = ForwardDiff.Chunk(x),
                            cfg = ForwardDiff.JacobianConfig(f, x, chunk))
    lj = first(logabsdet(ForwardDiff.jacobian(f, x, cfg)))
    isnan(lj) && handleNaN && return oftype(lj, -Inf)
    lj
end

"""
$(SIGNATURES)

Calculate the value and the log Jacobian determinant of `f` at `x`. `flatten` is
used to get a vector out of the result that makes `f` a bijection.
"""
function value_and_logjac_forwarddiff(f, x; flatten = identity, handleNaN = true,
                                      chunk = ForwardDiff.Chunk(x),
                                      cfg = ForwardDiff.JacobianConfig(flatten ∘ f, x, chunk))
    f(x), logjac_forwarddiff(flatten ∘ f, x; handleNaN = handleNaN, cfg = cfg)
end

@calltrans struct CustomTransform{G <: AbstractTransform, F, H, C} <: VectorTransform
    g::G
    f::F
    flatten::H
    cfg::C
end

_custom_f(g, f) = x -> f(transform(g, x))

_custom_chunk(g) = ForwardDiff.Chunk(zeros(dimension(g)))

_custom_cfg(g, f, flatten, chunk = _custom_chunk(g)) =
    ForwardDiff.JacobianConfig(flatten ∘ _custom_f(g, f), zeros(dimension(g)), chunk)

"""
$(SIGNATURES)

Wrap a custom transform `y = f(transform(g, x))`` in a type that calculates the log Jacobian
of ``∂y/∂x`` using `ForwardDiff` when necessary.

Usually, `g::TransformReals`, but when an integer is used, it amounts to the identity
transformation with that dimension.

`flatten` should take the result from `f`, and return a flat vector with no redundant
elements, so that ``x ↦ y`` is a bijection. For example, for a covariance matrix the
elements below the diagonal should be removed.

`chunk` and `cfg` can be used to configure `ForwardDiff.JacobianConfig`. `cfg` is used
directly, while `chunk = ForwardDiff.Chunk{N}()` can be used to obtain a type-stable
configuration.
"""
function CustomTransform(g::AbstractTransform, f, flatten;
                         chunk = _custom_chunk(g),
                         cfg = _custom_cfg(g, f, flatten, chunk))
    CustomTransform(g, f, flatten, cfg)
end

CustomTransform(n::Integer, f, flatten; kwargs...) =
    CustomTransform(as(Array, n), f, flatten; kwargs...)

dimension(t::CustomTransform) = dimension(t.g)

function transform_with(flag::NoLogJac, t::CustomTransform, x::AbstractVector, index)
    @unpack g, f = t
    f(first(transform_with(flag, g, x, index))), flag, index + dimension(t)
end

function transform_with(flag::LogJac, t::CustomTransform, x::AbstractVector, index)
    @unpack g, f, flatten, cfg = t
    index = firstindex(x)
    index′ = index + dimension(g)
    y, ℓ = value_and_logjac_forwarddiff(_custom_f(g, f), x[index:(index′ - 1)];
                                        flatten = flatten, cfg = cfg)
    y, ℓ, index′
end


text:
"""
    ElementalFieldModule

Module for elemental fields.
"""
module ElementalFieldModule

using ..FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import ..FieldModule: AbstractField, nents
import ..FieldModule.@add_Field_fields

"""
    ElementalField{T<:Number} <: AbstractField

Elemental field, meaning the entities are finite elements.

The values in the field are indexed by the element number.  This means  that
there needs to be one field per finite element set.
"""
mutable struct ElementalField{T<:Number} <: AbstractField
	@add_Field_fields()
end

"""
   ElementalField(data::FMat{T}=[]) where {T<:Number}

Constructor of elemental field. The values of the field are given by the array
on input, `data`. This array needs to have as many rows as there are elements,
and as many columns as there are degrees of freedom per element.
"""
function ElementalField(data::FMat{T}=[]) where {T<:Number}
	values = deepcopy(data)
	dofnums = 0*similar(values,FInt)
	is_fixed = similar(values,Bool)
	fill!(is_fixed, 0)
	fixed_values = zeros(T,size(values))
	nfreedofs = 0
	return ElementalField(values, dofnums, is_fixed, fixed_values, nfreedofs)
end

"""
    ElementalField(data::FVec{T}) where {T<:Number}

Constructor of elemental field. The values of the field are given by the vector
on input, `data`. This vector needs to have as many entries as there are elements;
there is just one degree of freedom per element.
"""
function ElementalField(data::FVec{T}) where {T<:Number}
    return ElementalField(reshape(data, length(data), 1))
end

"""
    nelems(self::ElementalField)::FInt = nents(self)

Provide the number of elements  in the elemental field.
"""
nelems(self::ElementalField)::FInt = nents(self)

end


text:
"Hierarchical Dirichlet Process Mixture Model Hyperparameters"
immutable HDPHyperparam <: AbstractHyperparam

  γ_a::Float64
  γ_b::Float64

  α_a::Float64
  α_b::Float64

  "default values"
  HDPHyperparam() = new(1.0, 1.0, 1.0, 1.0)

end

"Hierarchical Dirichlet Process Mixture Model Data Object"
type HDPData <: AbstractModelData

  # Energy
  energy::Float64

  # Distributions
  distributions::Vector{ConjugatePostDistribution}

  # Assignments
  assignments::Vector{Vector{Int}}

	# Weights
  weights::Vector{Vector{Float64}}

end

type HDPBuffer{T <: Real} <: AbstractModelBuffer

  # samples
  X::Vector{Vector{T}}

  # number of groups
  N0::Int

  # indecies of groups
  N0idx::Vector{Int}

  # number of samples per group
  Nj::Vector{Int}

  # indecies of samples per group
  Njidx::Array{Vector}

  # assignments
  Z::Vector{Vector{Int}}

  # number of active cluster
  K::Int

  # number of samples per group per cluster
  C::Array{Int, 2}

  # total number of tables
  totalnt::Array{Int}

  # number of clusters per table
  classnt::Array{Int}

  # distributions
  G::Vector{ConjugatePostDistribution}

  # base distribution
  G0::ConjugatePostDistribution

  # concentration parameter
  α::Float64

  # concentration parameters of clusters
  β::Array{Float64}

  # gamma
  γ::Float64
end

"HDP initialization using random assignments."
function init_random_hdp{T <: Real}(samples::Vector{Vector{T}}, G0::ConjugatePostDistribution; k::Int = 10)

    # TODO make sure the code is type-stable and Julia 0.4 compatible

    Z = Vector{Vector{Int}}(length(samples))
    G = ConjugatePostDistribution[]

    for i in 1:length(samples)

        cc = 1 + (collect(1:length(samples[i])) % k)
        Z[i] = vec( cc[randperm(length(samples[i]))] )

        for c in 1:k
            idx = find(Z[i] .== c)

            if length(idx) > 0
              if size(G, 1) < c
                  Gc = add_data(G0, samples[i][idx])
                  push!(G, Gc)
              else
                  G[c] = add_data(G[c], samples[i][idx])
              end
            end
        end

    end

    (Z, G)

end

"Single Gibbs sweep of HDP training using beta variables."
function gibbs_beta_hdp!(B::HDPBuffer)

    shuffle!(B.N0idx)

    prob = zeros(B.K + 1) * -Inf

    for j in B.N0idx

        # get samples of group
        data = B.X[j]
        z = B.Z[j]

        shuffle!(B.Njidx[j])

        for i in B.Njidx[j]

            cluster = z[i]

            # remove data item from model
            remove_data!(B.G[cluster], data[i])
            B.C[cluster, j] -= 1

            if isdistempty(B.G[cluster])

                # remove cluster
                B.G = B.G[[1:cluster-1; cluster+1:end]]

                for jj in 1:B.N0
                    B.Z[jj][B.Z[jj] .> cluster] -= 1
                end

                B.β = B.β[[1:cluster-1; cluster+1:end]]

                B.classnt = B.classnt[:,[1:cluster-1; cluster+1:end]]
                B.totalnt = B.totalnt[[1:cluster-1; cluster+1:end]]'

                B.K -= 1

                B.C = B.C[[1:cluster-1; cluster+1:end],:]
                prob = zeros(B.K + 1) * -Inf
            end

            # compute log likelihood
            for k in 1:B.K
              llh = logpred(B.G[k], data[i])[1]
              prior = B.C[k, j] + B.β[k] * B.α
              prob[k] = llh + log( prior )
            end

            #prob[1:B.K] = pmap(k -> logpred(B.G[k], data[i])[1] + log(B.C[k, j] + B.β[k] * B.α), 1:B.K)

            prob[B.K + 1] = logpred(B.G0, data[i])[1] + log( B.β[B.K+1] * B.α )

            prob = exp(prob - maximum(prob))

            c = rand_indices(prob)

            # add data to model
            B.Z[j][i] = c

            if c > B.K
                # create new cluster
                B.K += 1
                B.G = cat(1, B.G, deepcopy(B.G0))
                b = rand(Dirichlet([1, B.γ]))
                b = b * B.β[end]
                B.β = cat(1, B.β, 1)
                B.β[end-1:end] = b
                B.C = cat(1, B.C, zeros(Int, 1, B.N0))
                prob = zeros(B.K + 1) * -Inf
            end

            B.C[c, j] += 1
            add_data!(B.G[c], [data[i]])

        end

    end

    # sample number of tables
    kk = maximum([0, B.K - length(B.totalnt)])
    B.totalnt = cat(2, B.totalnt - sum(B.classnt, 1), zeros(Int, 1, kk))
    B.classnt = randnumtable(B.α .* B.β[:,ones(Int, B.N0)]', B.C')
    B.totalnt = B.totalnt + sum(B.classnt, 1)

    # update beta weights
    a = zeros(B.K + 1)
    a[1:end-1] = B.totalnt
    a[end] = B.γ

    B.β = rand(Dirichlet(a))

    B
end

function compute_energy!{T}(B::HDPData, X::Vector{Vector{T}})

  E = 1e-10
  numl = 0

  for j in length(X)

      # get samples of group
      data = X[j]

      for i in length(data)
        pp = 0.0
        c = 0

        for k = 1:length(B.distributions)
          p = exp( logpred( B.distributions[k], data[i] )[1] ) * B.weights[j][k]

          # only sum over actual values (excluding nans)
          if p == p
            pp += p
            c += 1
          end
        end

        E += pp
        numl += 1

      end

  end

  B.energy = log( E / numl )

end

function train_gibbs_hdp{T}(X::Vector{Vector{T}}, G0::ConjugatePostDistribution, Z::Vector{Vector{Int}}, G::Vector{ConjugatePostDistribution}, hyper::HDPHyperparam, K::Int;
                           α = 1.0, γ = 1.0, burnin = 0, thinout = 1, maxiter = 100)

  N0 = length(X)
  N0idx = collect(1:N0)

  Nj = [length(x) for x in X]
  Njidx = [collect(1:N) for N in Nj]

  C = reduce(hcat, [StatsBase.counts(z, 1:K) for z in Z])

  # init step
  β = ones(K) ./ K
  classnt = randnumtable(α * β[:,ones(Int, N0)]', C')
  totalnt = sum(classnt, 1)

  # set alpha vector of Dirichlet Distribution to sample β
  a = zeros(K + 1)
  a[1:end-1] = totalnt
  a[end] = γ

  # update beta
  β = rand(Dirichlet(a));

  B = HDPBuffer{T}(
    X,
    N0,
    N0idx,
    Nj,
    Njidx,
    Z,
    K,
    C,
    totalnt,
    classnt,
    G,
    G0,
    α,
    β,
    γ)

  # burn-in phase
  for iter in 1:burnin
    gibbs_beta_hdp!(B)

    totalnt = sum(B.classnt, 1)
    B.γ = random_concentration_parameter(B.γ, hyper.γ_a, hyper.γ_b, sum(B.totalnt), B.K, maxiter = 10)
    B.α = random_concentration_parameter(B.α, hyper.α_a, hyper.α_b, B.Nj, B.totalnt)
  end

  results = HDPData[]

  #println("starting")
  #addprocs(CPU_CORES)

  for iter = 1:maxiter

    for t in 1:thinout
      gibbs_beta_hdp!(B)

      totalnt = sum(B.classnt, 1)
      B.γ = random_concentration_parameter(B.γ, hyper.γ_a, hyper.γ_b, sum(B.totalnt), B.K, maxiter = 10)
      B.α = random_concentration_parameter(B.α, hyper.α_a, hyper.α_b, B.Nj, B.totalnt)
    end

		# compute weights
		W = map(j -> map(k -> B.C[k, j] + B.β[k] * B.α, 1:B.K), B.N0idx)

    model = HDPData(0.0, deepcopy(B.G), deepcopy(B.Z), deepcopy(W) )

    # update energy
    compute_energy!(model, X)

    # record results
    push!( results, model )

  end

  return results

end


text:
#==============================================================================================#
#                              interface.jl -- TypeTree interface                              #
#==============================================================================================#

#----------------------------------------------------------------------------------------------#
#                                   Empty generic functions                                    #
#----------------------------------------------------------------------------------------------#

function tt end

#----------------------------------------------------------------------------------------------#
#                                            Export                                            #
#----------------------------------------------------------------------------------------------#

export tt


text:
function partialtrace(rho,dimOfSubsystemToKeep)
  # First check if the density matrix is Hermitian
  if maximum(abs(rho'-rho))>8E-15
    error("This density matrix is not Hermitian!")
  end

  dimOfSubsystemToTraceOut = div(size(rho)[1] , dimOfSubsystemToKeep);

  matrix = zeros(Complex128,dimOfSubsystemToKeep,dimOfSubsystemToKeep);

  for i=1:dimOfSubsystemToTraceOut
    vector = zeros(Complex128,dimOfSubsystemToTraceOut,1);
    vector[i] = 1;
    fullVector = kron(eye(dimOfSubsystemToKeep),vector);
    matrix = matrix + fullVector' * rho * fullVector;
  end

  return matrix
end


text:
import Mosek

export solveTakouda
function solveTakouda(prob::Problem)

    n   = prob.N
	wlb = prob.wlb
	hlb = prob.hlb
	wub = prob.wub
	hub = prob.hub
	c   = prob.c
	area = prob.area
	W  = prob.W / 2
	H  = prob.H / 2

    β = prob.aspect

    Qˣ = Array(Float64, n, n)
    Qʸ = Array(Float64, n, n)
    Uˣ = Array(Float64, n, n)
    Uʸ = Array(Float64, n, n)
    for i in 1:n, j in 1:n
        Qˣ[i,j] = min(W, wub[i]+wub[j])
        Qʸ[i,j] = min(H, hub[i]+hub[j])
        Uˣ[i,j] = W - 0.5*(wlb[i] + wlb[j])
        Uʸ[i,j] = W - 0.5*(hlb[i] + hlb[j])
    end

    model = Model(solver=Mosek.MosekSolver())

    @variable(model, x[1:n])
    @variable(model, y[1:n])
    @variable(model, wlb[i] <= w[i=1:n] <= wub[i])
    @variable(model, hlb[i] <= h[i=1:n] <= hub[i])
    @variable(model, dˣ[1:n,1:n] >= 0)
    @variable(model, dʸ[1:n,1:n] >= 0)
    @variable(model, Sˣ[1:n,1:n] >= 0)
    @variable(model, Sʸ[1:n,1:n] >= 0)
    tn  = convert(Int,  n*( n-1)/2)
    ttn = convert(Int, tn*(tn-1)/2)
    nvar = 1+2tn+2ttn
    @variable(model, bin[1:nvar,1:nvar], SDP)
    for i in 1:nvar
        @constraint(model, bin[i,i] == 1)
    end

    # functions for referencing elements of σα; #assume access of upper tri
    function full_into_tn(i::Int,j::Int)
        ( 1 <= i <= n && 1 <= j <= n) || throw(BoundsError())
        i == j && error("Whoops")
        if i > j
            i,j = j,i
        end
        return sum(Int[x for x in (n-1):-1:(n-i+1)]) + (j - i)
    end

    function full_into_ttn(i::Int,j::Int)
        (1 <= i <= tn && 1 <= j <= tn) || throw(BoundsError())
        i == j && error("Whoops")
        if i > j
            i,j = j,i
        end
        return sum(Int[x for x in (tn-1):-1:(tn-i+1)]) + (j - i)
    end
    σ(i::Int,j::Int)  = bin[1,1+(full_into_tn(i,j))]
    α(i::Int,j::Int)  = bin[1,1+tn+(full_into_tn(i,j))]
    σα(i,j,k,l) = bin[1+full_into_tn(i,j),1+tn+full_into_tn(k,l)]
    σσ(i,j,k,l) = bin[1+   full_into_tn(i,j), 1+   full_into_tn(k,l)]
    αα(i,j,k,l) = bin[1+tn+full_into_tn(i,j), 1+tn+full_into_tn(k,l)]
    σσαα(i,j,k,l,p,q,r,s) = bin[1+2tn+   full_into_ttn(full_into_tn(i,j),full_into_tn(k,l)),
                                 1+2tn+ttn+full_into_ttn(full_into_tn(p,q),full_into_tn(r,s))]

    for i in 1:tn
        for j in (i+1):tn
            @constraint(model, bin[1,1+2tn+    full_into_ttn(i,j)] == bin[1+i,1+j])
            @constraint(model, bin[1,1+2tn+ttn+full_into_ttn(i,j)] == bin[1+tn+i,1+tn+j])
        end
    end
    for i in 1:n
        areaconstr = [     h[i]  √(area[i])
                      √(area[i])      w[i] ]
        @SDconstraint(model, areaconstr >= 0)

        aspconstr = [β[i]    w[i]
                     w[i] area[i]]
        @SDconstraint(model, aspconstr >= 0)
        aspconstr = [β[i]    h[i]
                     h[i] area[i]]
        @SDconstraint(model, aspconstr >= 0)

        sitbconstr = [0.5*(W-w[i])        x[i]
                             x[i]  0.5*(W-w[i])]
        @SDconstraint(model, sitbconstr >= 0)
        sitbconstr = [0.5*(H-h[i])        y[i]
                             y[i]  0.5*(H-h[i])]
        @SDconstraint(model, sitbconstr >= 0)

        for j in (i+1):n
            @constraint(model, dˣ[i,j] >= 0.5*(w[i]+w[j]) - 0.5*(1-σ(i,j))*Qˣ[i,j]) # 12
            @constraint(model, dʸ[i,j] >= 0.5*(h[i]+h[j]) - 0.5*(1+σ(i,j))*Qʸ[i,j]) # 19
            @constraint(model, dˣ[i,j] - 2Sˣ[i,j] == x[j] - x[i]) # 13
            @constraint(model, Sˣ[i,j] <= 0.25*(3 - σ(i,j) - α(i,j) - σα(i,j,i,j))*Uˣ[i,j]) # 15
            @constraint(model, Sˣ[i,j] + x[j] - x[i] >= 0) # 16
            @constraint(model, Sˣ[i,j] + x[j] - x[i] <= 0.25*(3 - σ(i,j) + α(i,j) + σα(i,j,i,j))*Uˣ[i,j]) # 17
            @constraint(model, dʸ[i,j] - 2Sʸ[i,j] == y[j] - y[i]) # 20
            @constraint(model, Sʸ[i,j] <= 0.25*(3 + σ(i,j) - α(i,j) + σα(i,j,i,j))*Uʸ[i,j]) # 22
            @constraint(model, Sʸ[i,j] + y[j] - y[i] >= 0) # 23
            @constraint(model, Sʸ[i,j] + y[j] - y[i] <= 0.25*(3 + σ(i,j) + α(i,j) - σα(i,j,i,j))*Uʸ[i,j]) # 24

            for k in (j+1):n
                @constraint(model, αα(i,j,i,k) - αα(i,j,j,k) + αα(i,k,j,k) +
                                    -σσ(i,j,i,k) - σσ(i,j,j,k) - σσ(i,k,j,k) +
                                     σσαα(i,j,i,k,i,j,i,k) - σσαα(i,j,i,k,i,j,j,k) + σσαα(i,j,i,k,i,k,j,k) +
                                     σσαα(i,j,j,k,i,j,i,k) - σσαα(i,j,j,k,i,j,j,k) + σσαα(i,j,j,k,i,k,j,k) +
                                     σσαα(i,k,j,k,i,j,i,k) - σσαα(i,k,j,k,i,j,j,k) + σσαα(i,k,j,k,i,k,j,k) == 1)
            end
            # symmetry breaking
            @constraint(model, dˣ[i,j] >= 0.25*(wlb[i]+wlb[j]) * (1 + σ(i,j)))
            @constraint(model, dʸ[i,j] >= 0.25*(hlb[i]+hlb[j]) * (1 - σ(i,j)))
            # # S valid inequalities
            @constraint(model, x[i] + 0.5w[i] <= x[j] - 0.5w[j] + 0.5*(3 - σ(i,j) - α(i,j) - σα(i,j,i,j))*W)
            @constraint(model, x[j] + 0.5w[j] <= x[i] - 0.5w[i] + 0.5*(3 - σ(i,j) + α(i,j) + σα(i,j,i,j))*W)
            @constraint(model, y[i] + 0.5h[i] <= y[j] - 0.5h[j] + 0.5*(3 + σ(i,j) - α(i,j) - σα(i,j,i,j))*H)
            @constraint(model, y[j] + 0.5h[j] <= y[i] - 0.5h[i] + 0.5*(3 + σ(i,j) + α(i,j) + σα(i,j,i,j))*H)
        end
    end

    @objective(model, Min, dot(c,dˣ + dʸ))
    buildInternalModel(model)
    elapse = @elapsed (stat = solve(model))
    println("Objective value = $(getObjectiveValue(model)) in $elapse sec")
    return getObjectiveValue(model)
end


text:
using DataFrames
using CSV
using Plots
#%%
data_rk = open("rk_ode.dat")
data_rk4 = open("rk4_ode.dat")
data_eu = open("eu_ode.dat")

df_rk=CSV.read(data_rk, header=false, delim=" ")
df_rk4=CSV.read(data_rk4, header=false, delim=" ")
df_eu=CSV.read(data_eu, header=false, delim=" ")

#%%
plot(df_rk[1],df_rk[2], label="Runge-Kutta", dpi=600)
plot!(df_rk4[1],df_rk4[2], label="RK 4th order", dpi=600)
plot(df_eu[1],df_eu[2], label="Euler", dpi=600)
# xaxis!(:log10)
# yaxis!(:log10)

max_stars_repo_name: aurelio-amerio/corso-cpp


text:
struct FixedNesterovSequence{R} end

"""
    FixedNesterovSequence(R::Type)

Construct the standard sequence of extrapolation coefficients, of type `R`, used in Nesterov-accelerated gradient methods.

Argument `R` should be a floating point type, and specifies the type of the elements in the sequence.

See also: [`SimpleNesterovSequence`](@ref), [`ConstantNesterovSequence`](@ref).
"""
FixedNesterovSequence(R) = FixedNesterovSequence{R}()

function Base.iterate(::FixedNesterovSequence{R}, t=R(1)) where R
    t_next = (1 + sqrt(1 + 4 * t^2)) / 2
    return (t - 1) / t_next, t_next
end

Base.IteratorSize(::Type{<:FixedNesterovSequence}) = Base.IsInfinite()
Base.eltype(::Type{FixedNesterovSequence{R}}) where R = R

struct SimpleNesterovSequence{R} end

"""
    SimpleNesterovSequence(R::Type)

Construct an iterator yielding the sequence `(k - 1) / (k + 2)`, of type `R`, for `k >= 1`.

This is the "simplified" sequence of extrapolation coefficients used in Nesterov-accelerated gradient methods.
Argument `R` should be a floating point type, and specifies the type of the elements in the sequence.

See also: [`FixedNesterovSequence`](@ref), [`ConstantNesterovSequence`](@ref).
"""
SimpleNesterovSequence(R) = SimpleNesterovSequence{R}()

Base.iterate(::SimpleNesterovSequence{R}, k=1) where {R} = R(k - 1) / (k + 2), k + 1

Base.IteratorSize(::Type{<:SimpleNesterovSequence}) = Base.IsInfinite()
Base.eltype(::Type{SimpleNesterovSequence{R}}) where R = R

"""
    ConstantNesterovSequence(m::R, stepsize::R)

Construct a constant iterator yielding `(1 - sqrt(m * stepsize)) / (1 + sqrt(m * stepsize))` with type `R`.

This sequence of extrapolation coefficients is commonly used in Nesterov-accelerated gradient methods
when applied to strongly convex functions.

See also: [`FixedNesterovSequence`](@ref), [`SimpleNesterovSequence`](@ref).
"""
function ConstantNesterovSequence(m::R, stepsize::R) where R
    k_inverse = m * stepsize
    return repeated((1 - sqrt(k_inverse)) / (1 + sqrt(k_inverse)))
end

mutable struct AdaptiveNesterovSequence{R}
    m::R
    stepsize::R
    theta::R
end

"""
    AdaptiveNesterovSequence(m::R)

Construct an object that yields a sequence of extrapolation coefficients, commonly used in Nesterov-accelerated gradient methods.

Argument `m` is usually the convexity modulus of the interested function, so it should be zero for barely convex functions,
and greater than zero for strongly convex ones. The extrapolation coefficients can be extracted by calling [`next!`](@ref) on the
object.

Note that this type, as opposed to e.g. [`FixedNesterovSequence`](@ref), is not a standard Julia iterator, in that it requires
the sequence of stepsizes (input to [`next!`](@ref)) in order to produce the output sequence.

Special cases: if `next!` is invoked with a constant stepsize `gamma`, then
- If `m == 0`, the output sequence will be the same as for `FixedNesterovSequence(R)`.
- If `m > 0`, the output sequence will be the same as for `ConstantNesterovSequence(m, gamma)`.

See also: [`next!`](@ref).
"""
AdaptiveNesterovSequence(m::R) where R = AdaptiveNesterovSequence{R}(m, -R(1), -R(1))

"""
    next!(seq::AdaptiveNesterovSequence{R}, stepsize::R)

Push a new stepsize value into the sequence, and return the next extrapolation coefficient.

See also: [`AdaptiveNesterovSequence`](@ref).
"""
function next!(seq::AdaptiveNesterovSequence{R}, stepsize::R) where R
    if seq.stepsize < 0
        seq.stepsize = stepsize
        seq.theta = seq.m > 0 ? sqrt(seq.m * stepsize) : R(1)
    end
    b = seq.theta^2 / seq.stepsize - seq.m
    delta = b^2 + 4 * (seq.theta^2) / (seq.stepsize * stepsize)
    theta = stepsize * (-b + sqrt(delta)) / 2
    beta = stepsize * (seq.theta) * (1 - seq.theta) / (seq.stepsize * theta + stepsize * seq.theta ^ 2)
    seq.stepsize = stepsize
    seq.theta = theta
    return beta
end

struct NesterovExtrapolation{S} end

NesterovExtrapolation(S) = NesterovExtrapolation{S}()
NesterovExtrapolation() = NesterovExtrapolation{SimpleNesterovSequence}()

acceleration_style(::Type{<:NesterovExtrapolation}) = NesterovStyle()

initialize(::NesterovExtrapolation{S}, x) where S = Iterators.Stateful(S{real(eltype(x))}())


text:
function cos_test(x, y)
    cos(x) * cos(y)
end

function u_test(x, y)
    exp(cos(x) + cos(2y))
end

function dudx_test(x,y)
    - sin(x) * u_test(x,y)
end

function dudy_test(x,y)
    - 2sin(2y) * u_test(x,y)
end

function f_test(x, y)
    (sin(x)^2 - cos(x) + 4sin(2y)^2 - 4cos(2y)) * u_test(x,y)
end

function evaluate_function_on_grid(grid, ufunc, u)
    @assert size(u,1) == length(grid.x)
    @assert size(u,2) == length(grid.y)
    @inbounds for j in 1:size(u,2)
        for i in 1:size(u,1)
            u[i,j]  = ufunc(grid.x[i], grid.y[j])
        end
    end
end


text:
# functions in readme file
# Claudia October 2015

##using PhyloNetworks
using DataFrames

if !isdefined(:individualtest) individualtest = false; end

if(individualtest)
    include("../src/types.jl")
    include("../src/functions.jl")
end

if isdefined(:PhyloNetworks)
    PhyloNetworks.CHECKNET || error("need CHECKNET==true in PhyloNetworks to test snaq in test_correctLik.jl")
else
    CHECKNET || error("need CHECKNET==true to test snaq in test_correctLik.jl")
end

try
    d=readTrees2CF("../examples/treefile.txt", CFfile = "none", writeTab=false, writeSummary=false)
    T=readTopologyLevel1("../examples/startTree.txt");
    writeTopologyLevel1(T)
    # snaq
    net1=snaq!(T,d,filename="", runs=2)
    otherNet = undirectedOtherNetworks(net1)
    net2=snaq!(T,d,hmax=2, filename="", runs=2)
    d=readTableCF("../examples/tableCF.txt", summaryfile="");
    # snaq
    net1=snaq!(T,d,filename="", runs=2)
    net2=snaq!(T,d,hmax=2, filename="", runs=2)
    rootatnode!(net1, "4")
    writeTopologyLevel1(net1)
    writeTopologyLevel1(net1,di=true)
    # topologyMaxQ
    net1topo = readTopology("(2,(4,(3,(5,(6,#H1)))),(1)#H1);")
    net1par = topologyMaxQPseudolik!(net1topo,d)
    # topologyQ
    net1withBL = readTopology("(2,(4,(3,(5,(6,#H6:1.0::0.3):5.006):0.518):0.491):1.533,(1)#H6:1.0::0.7);");
    topologyQPseudolik!(net1withBL,d)
catch
    error("error in README with input trees")
end

try
    T=readTopologyLevel1("../examples/startTree.txt");
    d=readTableCF("../examples/tableCF.txt", summaryfile="")
    # snaq
    net1=snaq!(T,d,filename="", runs=2)
    net2=snaq!(T,d,hmax=2, filename="", runs=2)
    rootatnode!(net1, "4")
    writeTopologyLevel1(net1)
    writeTopologyLevel1(net1,di=true)
    # topologyMaxQ
    net1topo = readTopology("(2,(4,(3,(5,(6,#H1)))),(1)#H1);")
    net1par = topologyMaxQPseudolik!(net1topo,d)
    # topologyQ
    net1withBL = readTopology("(2,(4,(3,(5,(6,#H6:1.0::0.3):5.006):0.518):0.491):1.533,(1)#H6:1.0::0.7);");
    topologyQPseudolik!(net1withBL,d)
catch
    error("error in README with CF table")
end

println("NO ERRORS!")

if(false)
try
    df = readtable("../examples/tableCFCI.csv")
    bootnet = bootsnaq(T, df, hmax=1, nrep=3, runs=2) # need to make sure does not produce files
    BStable, tree1 = treeEdgesBootstrap(bootnet,net1)
    BSn, BSe, BSc, BSgam, BSedgenum = hybridBootstrapSupport(bootnet, net1)
catch
    error("error in README for bootstrap")
end
end

#plot obsCF vs expCF
## df = fittedQuartetCF(d)
## using Gadfly
## p = plot(df,layer(x="obsCF1",y="expCF1",Geom.point,Theme(default_color=colorant"orange")),layer(x="obsCF2",y="expCF2",Geom.point,Theme(default_color=colorant"purple")),layer(x="obsCF3",y="expCF3",Geom.point,Theme(default_color=color("blue"))),layer(x=0:1,y=0:1),Geom.line,Theme(default_color=color("black")))


text:
import Test: @inferred

@testset "Substitution" begin
    Mod.@polyvar x[1:3]

    @test subs(2, x[1]=>3) == 2
    @test iszero((x[1]-x[1])(x[1]=>x[2]))
    @test subs(CustomPoly(x[1]+x[2]), x[2]=>x[1]) == 2x[1]

    a = (x[1])(x[1]=>x[2])
    b = x[2]
    @test (x[1])(x[1]=>x[2]) == x[2]

    p = x[1]*x[2]*x[3]
    @test convert(Int, p(x => (1, 2, 3))) == 6

    p = x[1]^2 + x[1]*x[3] - 3
    @test convert(Int, p(x => (5, x[1]*x[2], 4))) == 42

    p = x[1]^2 + x[2]^2
    q = p(x[1:2] => [1 -1; 1 1] * vec(x[1:2]))
    @test q == 2p

    t = 2.0 * x[1] * x[2]^2
    @test t(x => [1, 2, 3]) == 8.0
    @test t((x[2], x[1]) => [1, 3]) == 6.0

    p = x[1] + x[2] + 2*x[1]^2 + 3*x[1]*x[2]^2
    #@inferred p((x[1], x[2]) => (1.0, 2.0))

    @inferred subs(p, x[2] => 2.0)
    @test subs(p, x[2] => 2.0) == 13x[1] + 2 + 2x[1]^2
    @inferred subs(p, x[2] => x[1])
    @inferred subs(p, x[2] => x[1]) == 2x[1] + 2x[1]^2 + 3x[1]^3

    p = x[1] - x[1]
    @test iszero(p)
    @inferred p(x[1] => 1)
    @test iszero(p(x[1] => 1))
    @inferred subs(p, x[1] => 1)
    @test iszero(subs(p, x[1] => 1))

    q = (x[1] + 1) / (x[1] + 2)
    @test isapproxzero(q(x[1] => -1))
    @test !isapproxzero(q(x[1] => 1))
    @test q(x[1] => 1) ≈ 2/3

    q = (x[1] + x[3]) / (x[2] - 1)
    @test subs(q, x[1:2] => (x[2], x[1])) == (x[2] + x[3]) / (x[1] - 1)

    # Taken from perf/runbenchmark.jl
    p = x[1] + x[2] + 2x[1]^2 + x[2]^3

    varst = (x[1], x[2])
    varsv = [x[1], x[2]]
    valst = (1.0, 2.0)
    valsv = [1.0, 2.0]

    @test variable.(Tuple(varsv)) == varst
    @test variable.(Tuple(varsv)) isa Tuple{<:AbstractVariable, <:AbstractVariable}

    @test p(variables(p) => valst) == 13
    @test p(variables(p) => valsv) == 13
    @test p(varst => valst) == 13
    @test p(varsv => valsv) == 13

    Mod.@polyvar x y
    @test subs([x^2 + y, x + y], x => y) == [y^2+y, 2y]
end


text:
using Flux
using Flux.Tracker, Test, NNlib
using Flux.Tracker: TrackedReal, gradient, gradcheck, grad, checkpoint, forwarddiff
using NNlib: conv, depthwiseconv
using Printf: @sprintf
using LinearAlgebra: diagm, dot, LowerTriangular, norm
using Statistics: mean, std
using Random
# using StatsBase

gradtest(f, xs::AbstractArray...) = gradcheck((xs...) -> sum(sin.(f(xs...))), xs...)
gradtest(f, dims...) = gradtest(f, rand.(Float64, dims)...)
@testset "Tracker" begin
@test gradtest((x, W, b) -> σ.(W*x .+ b), 5, (2,5), 2)
@test gradtest((x, W, b) -> σ.(W*x .+ b), (5,3), (2,5), 2)
@test gradtest((x, W, b) -> logσ.(W*x .+ b), 5, (2,5), 2)
@test gradtest((x, W, b) -> logσ.(W*x .+ b), (5,3), (2,5), 2)
@test gradtest((w, x) -> w'*x, randn(Float64,10, 2), randn(Float64,10))
@test gradtest((w, x) -> w*x', randn(Float64,5,5), randn(Float64,5,5))
@test gradtest(x -> sum(x, dims = (2, 3)), (3,4,5))
@test gradtest(x -> sum(x, dims = 1), randn(Float64,2,3))
@test gradtest(x -> sum(x, dims = [1,2]), randn(Float64,2,3))
@test gradtest(x -> sum(x), randn(Float64,2,3))
@test gradtest(x -> prod(x, dims=(2, 3)), (3,4,5))
@test gradtest(x -> prod(x), (3,4,5))

@test gradtest(x -> softmax(x).*(1:3), 3)
@test gradtest(x -> softmax(x).*(1:3), (3,5))
@test gradtest(x -> logsoftmax(x).*(1:3), 3)
@test gradtest(x -> logsoftmax(x).*(1:3), (3,5))

@test gradtest(Flux.mse, rand(5,5), rand(5, 5))
@test gradtest(Flux.crossentropy, rand(5,5), rand(5, 5))

@test gradtest(x -> x', rand(5))

@testset "indexing & slicing" begin
  gradtest(x->view(x, 1:2, 1:2), rand(4, 4))
end

function promotiontest(f, A, B, C)
  r0 = f(A, B, C)
  r1 = f(param(A), B, C)
  r2 = f(A, param(B), C)
  r3 = f(A, B, param(C))
  r4 = f(param(A), param(B), param(C))

  @test !isa(r0, TrackedArray)
  @test all(isa.([r1,r2,r3,r4], TrackedArray))
  @test r1 == r2 == r3 == r4
  @test r0 == Flux.data(r4)
end

@testset "concat" begin
  cat1(x...) = cat(x..., dims = 1)
  cat2(x...) = cat(x..., dims = 2)

  @testset for vcatf in [vcat, cat1]
    @test gradtest(vcatf, rand(5), rand(3))
    @test gradtest(vcatf, rand(5), rand(3), rand(8))
    @test gradtest(vcatf, rand(5)', rand(5)')
    @test gradtest(vcatf, rand(5,2), rand(3,2), rand(8,2))
    @test gradtest(vcatf, rand(5,2,3), rand(3,2,3), rand(8,2,3))
    @test gradtest(vcatf, rand(5), rand(3,1))
    @test gradtest(vcatf, rand(5)', rand(2,5))
  end

  @testset for hcatf in [hcat, cat2]
    @test gradtest(hcatf, rand(5), rand(5))
    @test gradtest(hcatf, rand(5)', rand(5)')
    @test gradtest(hcatf, rand(2,5), rand(2,3), rand(2,8))
    @test gradtest(hcatf, rand(2,5,3), rand(2,3,3), rand(2,8,3))
    @test gradtest(hcatf, rand(5), rand(5), rand(5,2))
    @test gradtest(hcatf, rand(5)', rand(1,3))
    @test gradtest(hcatf, rand(5), rand(5,2))
end

  @testset for catf in [vcat, cat1, hcat, cat2, (x...) -> cat(x..., dims = 3), (x...) -> cat(x..., dims = (1,2))]
    @test gradtest(catf, rand(5))
    @test gradtest(catf, rand(5)')
    @test gradtest(catf, rand(2,5))
    @test gradtest(catf, rand(2,5,3))
  end

  @test gradtest((x...) -> cat(x..., dims = 3), rand(2,5,2), rand(2,5,3), rand(2,5,4))

  @testset "cat($dim, ...)" for dim in 3:5
    catdim = (x...) -> cat(x..., dims = dim)
    @test gradtest(catdim, rand(5), rand(5), rand(5))
    @test gradtest(catdim, rand(2,5), rand(2,5), rand(2,5))
    @test gradtest(catdim, rand(2,5,3), rand(2,5,3), rand(2,5,3))
  end

  @test !isa(vcat(rand(2)), TrackedArray)
  @test !isa(hcat(rand(2)), TrackedArray)
  @test !isa(cat(rand(2), dims=1), TrackedArray)

  @test gradtest((a,b)->cat(a, b, dims = (2,3,5)), rand(2,3), rand(2,4,2,1))

  @testset "promotiontest" begin
    @testset for fcat in [hcat, vcat, (x...) -> cat(x..., dims = 3), (x...) -> cat(x..., dims = (1,2))]
      promotiontest(fcat, rand(2), rand(2), rand(2))
      promotiontest(fcat, rand(2)', rand(2)', rand(2)')
      promotiontest(fcat, rand(2,2), rand(2,2), rand(2,2))
      promotiontest(fcat, rand(2,2,2), rand(2,2,2), rand(2,2,2))
    end

    promotiontest(vcat, rand(1,2), rand(2)', rand(2,2))
    promotiontest(hcat, rand(2,1), rand(2), rand(2,2))
    promotiontest(vcat, rand(3,4,5), rand(1,4,5), rand(2,4,5))
    promotiontest(hcat, rand(4,3,5), rand(4,1,5), rand(4,2,5))
    promotiontest((x...) -> cat(x..., dims = 3), rand(4,5,3), rand(4,5,1), rand(4,5,2))
  end

end

@test gradtest(x -> permutedims(x, [3,1,2]), rand(4,5,6))

@test gradtest(x -> repeat(x; inner=2), rand(5))
@test gradtest(x -> repeat(x; inner=2, outer=3), rand(5))
@test gradtest(x -> repeat(x; inner=(2,2,1), outer=(1,1,3)), rand(5,4,3))

@test gradtest(kron, rand(5), rand(3))
@test gradtest(kron, rand(5), rand(3), rand(8))
@test gradtest(kron, rand(5,1), rand(3,1))
@test gradtest(kron, rand(5,1), rand(3,1), rand(8,1))
@test gradtest(kron, rand(5,2), rand(3,2), rand(8,2))

@test gradtest(x -> diagm(0 => x), rand(3))

@test gradtest(W -> inv(log.(W * W)), (5,5))
@test gradtest((A, B) -> A / B , (1,5), (5,5))
@test gradtest((A, B) -> log.(A * A) / exp.(B * B), (5,5), (5,5))
@test gradtest((A, B) -> log.(A * A) \ exp.(B * B), (5,5), (5,5))

@testset "mean" begin
  @test gradtest(mean, rand(2, 3))

  @test gradtest(x -> mean(x, dims=1), rand(2, 3))
  @test gradtest(x -> mean(x, dims=2), rand(2, 3))
  @test gradtest(x -> mean(x, dims=3), rand(2, 3, 4))

  @test gradtest(x -> mean(x, dims=[1, 2]), rand(2, 3, 4))
end

@testset "maximum" begin
  @test gradtest(maximum, rand(2, 3))

  @test gradtest(x -> maximum(x, dims=1), rand(2, 3))
  @test gradtest(x -> maximum(x, dims=2), rand(2, 3))
  @test gradtest(x -> maximum(x, dims=3), rand(2, 3, 4))

  @test gradtest(x -> maximum(x, dims=[1, 2]), rand(2, 3, 4))
end

@testset "minimum" begin
  @test gradtest(minimum, rand(2, 3))

  @test gradtest(x -> minimum(x, dims=1), rand(2, 3))
  @test gradtest(x -> minimum(x, dims=2), rand(2, 3))
  @test gradtest(x -> minimum(x, dims=3), rand(2, 3, 4))

  @test gradtest(x -> minimum(x, dims=[1, 2]), rand(2, 3, 4))
end

@test gradtest(x -> std(x), rand(5,5))
@test gradtest(x -> std(x, dims = 1), rand(5,5))

@test gradtest((x, y) -> x .* y, rand(5), rand(5))
@test gradtest(dot, rand(5), rand(5))

@test gradtest(norm, rand(5))

@test gradtest(rand(5)) do x
  y = x.^2
  2y + x
end

@test gradtest(conv, rand(10, 3, 2), randn(Float64,2, 3, 2))
@test gradtest(conv, rand(10, 10, 3, 2), randn(Float64,2, 2, 3, 2))
@test gradtest(conv, rand(10, 10, 10, 3, 2), randn(Float64,2, 2, 2, 3, 2))

@test gradtest(depthwiseconv, rand(10,10,3,2), randn(2, 2, 2, 3))

@test gradtest(x -> maxpool(x, (2,2)), rand(10, 10, 3, 2))
@test gradtest(x -> maxpool(x, (2,2,2)), rand(10, 10, 10, 3, 2))

@test gradtest(x -> meanpool(x, (2,2)), rand(10, 10, 3, 2))
@test gradtest(x -> meanpool(x, (2,2,2)), rand(5, 5, 5, 3, 2))

@test gradtest(x -> Float64.(x), 5)

@testset "equality & order" begin
    # TrackedReal
    @test param(2)^2 == param(4)
    @test param(2)^2 == 4
    @test 4 == param(2)^2

    @test param(2)^2 ≈ param(4)
    @test param(2)^2 ≈ 4
    @test 4 ≈ param(2)^2

    @test (param([1,2,3]) .< 2) == [true, false, false]
    @test (param([1,2,3]) .<= 2) == [true, true, false]
    @test (2 .> param([1,2,3])) == [true, false, false]
    @test (2 .>= param([1,2,3])) == [true, true, false]

    # TrackedArray
    @test param([1,2,3]).^2 == param([1,4,9])
    @test [1,2,3].^2 == param([1,4,9])
    @test param([1,2,3]).^2 == [1,4,9]

    @test param([1,2,3]).^2 ≈ param([1,4,9])
    @test [1,2,3].^2 ≈ param([1,4,9])
    @test param([1,2,3]).^2 ≈ [1,4,9]
end

@testset "reshape" begin
  x = reshape(param(rand(2,2,2)), 4, 2)
  @test x isa TrackedArray
  @test size(x) == (4,2)
  x = reshape(param([1]), (1,:))
  @test x isa TrackedArray
  @test size(x) == (1,1)
  x = reshape(param(rand(2)), (2,:))
  @test x isa TrackedArray
  @test size(x) == (2,1)
  x = reshape(param(rand(2,2)), (1,:,2))
  @test x isa TrackedArray
  @test size(x) == (1,2,2)
end

@testset "Intermediates" begin
  x = param([1])
  l = sum((x .+ x).^2)
  Flux.back!(l, once = false)
  @test x.grad == [8]
  x.grad .= 0
  Flux.back!(l, once = false)
  @test x.grad == [8]
end

@testset "Fallbacks" begin
  xs = param([1 2; 3 4])
  @test similar(xs) isa Matrix{Float64}
end

@test @sprintf("%.2f", sum(param([1,2,3]))) == "6.00"

@inferred NNlib.conv(param(rand(10,10,3,2)),randn(Float64,2,2,3,4))

b = param(rand())
Tracker.back!(b)
@test Tracker.grad(b) == 1

@testset "collect" begin
  x, y = param(2), param(3)
  xy = Tracker.collect([x, y])
  @test xy isa TrackedArray{Float64}
  z = xy[1]*xy[2]
  back!(z)
  @test grad.((x,y)) == (3, 2)

  @test gradient(2, 3) do x, y
    xy = Tracker.collect([x, y])
    xy[1]*xy[2]
  end == (3, 2)
end

# Gradient Hooks
@testset "Hooks" begin
  x = param(2)
  y = Tracker.hook(-, x)
  back!(y)
  @test grad(x) == -1
end

@testset "Checkpointing" begin
  count = 0
  function mul(a, b)
    count += 1
    a * b
  end
  @test gradient(x -> mul(5, x), 3)[1] == 5
  @test count == 1
  @test gradient(x -> checkpoint(mul, 5, x), 3)[1] == 5
  @test count == 3
end

@testset "Updates" begin
  xs = param([1, 2, 3])
  Tracker.update!(xs, param([4, 5, 6]))
  @test xs == [5, 7, 9]
  x = param(3)
  Tracker.update!(x, param(4))
  @test x == 7
end

@testset "Params" begin
  W = param(randn(5, 10))
  x = rand(10)
  dW = gradient(W -> sum(W*x), W)[1]
  gs = gradient(() -> sum(W*x), Tracker.Params([W]))
  @test gs[W] == dW
end

@testset "Forward" begin
  @test @inferred(Tracker.forward_jacobian(x -> [sum(x)], rand(5,5), Val(12)))[2] ==
    reshape(ones(25), :, 1)
  @test gradient([2, 3]) do x
    forwarddiff(x) do x
      x[1]*x[2]
    end
  end == ([3, 2],)
end

end #testset


text:
using Optim
using Calculus

struct CaseModelFitSEI3R <: CaseModelFitResult
    ρ::Float64
    E0::Float64
    info::Matrix{Float64}
    inputs::NamedTuple
    fit::Optim.MultivariateOptimizationResults
end

logit(x::Float64) = log(x/(1-x))
invlogit(x::Float64) = exp(x)/(1+exp(x))

function poissonLogLik( x::Vector{Int64}, λ::Vector{Float64})
    return sum(poissonLogLik.(x,λ))
end

function poissonLogLik(x::Int64,λ::Float64)
    return logpdf(Poisson(λ),x)
end

"""
function fitCaseModel(cases::Vector{Int64},N::Int64,
    ρ::Float64,ϕ0::Float64,E::Float64,d0::SEI3R)
"""
function fitCaseModel(cases::Vector{Int64},N::Int64,
    ρ::Float64,ϕ0::Float64,E::Float64,d0::SEI3R)

    d = deepcopy(d0)

    ntime = length(cases)
    # dvec = fill(d0,ntime)

    function f(param::Vector{Float64})
        d.β[1] = invlogit(param[2] + logit(d0.β[1]) )
        s = estimatedStates(ntime,N,cases[1],exp(param[3]),
                            invlogit(param[1]),d)
        caseIntensity = s.I1*invlogit(param[1])
        # println(caseIntensity)
        return -poissonLogLik(diff(cases),caseIntensity[2:end])
        # return caseIntensity
    end

    fit = optimize(f,[logit(ρ),logit(ϕ0),log(E)],NelderMead())
    d.β[1] = invlogit(fit.minimizer[2] + logit(d0.β[1]) )
    info = Calculus.hessian(f,fit.minimizer)
    return (cases=cases,N=N,ρ=invlogit(fit.minimizer[1]),
            E0=exp(fit.minimizer[3]),d=d,info=info,fit=fit)
end

###########
"""
function fitCaseModel(cases::Vector{Int64},x::Matrix{Float64},
    N::Int64,ρ::Float64,ϕ0::Float64,E::Float64,d0::SEI3R)

Fit a model to confirmed cases where the β1 transmission parameter is allowed to depend on covariates x.  The dependence is through a logit link function.
"""
function fitCaseModel(cases::Vector{Int64},x::Matrix{Float64},
    N::Int64,ρ::Float64,E::Float64,d0::SEI3R)

    d = deepcopy(d0)
    nc = size(x,2)

    ntime = length(cases)
    # dvec = fill(d0,ntime)

    function f(param::Vector{Float64})
        ϕ = hcat(param[3:end])
        betachange = (x*ϕ)[:,1]
        s = estimatedStates(N,cases[1],exp(param[2]),
                            invlogit(param[1]),betachange,d)
        caseIntensity = s.I1*invlogit(param[1])
        # println(caseIntensity)
        return -poissonLogLik(diff(cases),caseIntensity[2:end])
        # return caseIntensity
    end

    fit = optimize(f,[logit(ρ);log(E);zeros(nc)],NelderMead())
    # d.β[1] = invlogit(fit.minimizer[2] + logit(d0.β[1]))
    info = Calculus.hessian(f,fit.minimizer)
    return (cases=cases,N=N,ρ=invlogit(fit.minimizer[1]),
            E0=exp(fit.minimizer[3]),d=d0,info=info,fit=fit)
end

"""
predictCases(fit,ntime)

Predict number of cases using a fit object.
"""
function predictCases(fit,ntime::Int64)
    s = estimatedStates(ntime,fit.N,fit.cases[1],fit.E0,fit.ρ,fit.d)
    return cumsum(s.I1*fit.ρ)
end

function predictCases(fit,x::Matrix{Float64})
    ϕ = hcat(fit.fit.minimizer[3:end])
    betachange = (x*ϕ)[:,1]
    s = estimatedStates(fit.N,fit.cases[1],fit.E0,fit.ρ,
        betachange,fit.d)
    return cumsum(s.I1*fit.ρ)
end


text:
# deprecations should also include a comment stating the version
# for which the deprecation should be removed.

# Deprecated to fix spelling. Can be removed for version 2.0.
@deprecate simplecycles_hadwick_james simplecycles_hawick_james