text:
const T = Float32

@testset "functions" for i = 1:5

@testset "activation" begin
    x = parameter(randn(T,10,5))
    for i = 1:length(x.data)
        abs(x.data[i]) < 0.1 && (x.data[i] += 1)
    end
    for f in (relu,sigmoid,tanh)
        checkgrad(()->f(x), x)
    end
end

@testset "cnn" begin
    x = parameter(randn(T,10,10))
    f = Conv1d(T, 5, 10, 7, padding=2)
    checkgrad(()->f(x,[3,7]), x, f.W, f.b)
end

@testset "loss" begin
    @testset "softmax_crossentropy" begin
        p = Var(rand(1:10,5))
        q = parameter(rand(T,10,5))
        checkgrad(()->softmax_crossentropy(p,q), p, q)
    end

    @testset "l2" begin
        x = parameter(rand(T,10,5))
        # @test_grad l2(x,0.01) x
    end

    @testset "mse" begin
        x1 = parameter(rand(T,10,5))
        x2 = parameter(rand(T,10,5))
        #@test_grad mse(x1,x2) x1 x2
    end
end

@testset "math" begin
    x = parameter(rand(T,10,5) .+ T(1))
    #@test_function exp x
    #@test_function log x

    x1 = parameter(randn(T,10,5))
    x2 = parameter(randn(T,10,5))
    checkgrad(()->x1+x2, x1, x2)
    checkgrad(()->x1-x2, x1, x2)
    checkgrad(()->x1.*x2, x1, x2)

    A = parameter(randn(T,10,5))
    B = parameter(randn(T,5,7))
    checkgrad(()->A*B, A, B)

    x1 = parameter(randn(T,10,5))
    x2 = parameter(randn(T,10,1))
    checkgrad(()->x1.+x2, x1, x2)
    checkgrad(()->x1.-x2, x1, x2)
    checkgrad(()->x1.*x2, x1, x2)
end

@testset "reduction" begin
    x = parameter(randn(T,10,10)*T(10))
    @testset "max" begin
        for dim = 1:2
            checkgrad(()->max(x,dim), x)
        end
        checkgrad(()->max(x,[2,5,3]), x)
    end
end

@testset "rnn" begin
    x = parameter(randn(T,10,10))
    for nlayers = 1:1
        f = LSTM(T, 10, 5, nlayers, 0.0, true)
        checkgrad(()->f(x,[5,3,2]), x, f.weights...)
    end
end

@testset "blas" begin
    A = parameter(randn(T,10,5))
    x = parameter(randn(T,10))
    checkgrad(()->gemv('T',1,A,x), A, x)

    A = parameter(randn(T,10,5))
    B = parameter(randn(T,10,5))
    checkgrad(()->gemm('N','T',1,A,B), A, B)

    A = parameter(randn(T,10,5,7))
    B = parameter(randn(T,10,5,7))
    #test_gradient(gemm_batch, 'N', 'T', 1, A, B)
    #test_cuda(gemm_batch, 'N', 'T', 1, A, B)
end

@testset "concat" begin
    x1 = parameter(randn(T,10,5,2))
    x2 = parameter(randn(T,10,5,2))
    x3 = parameter(randn(T,10,5,2))
    for dim = 1:3
        checkgrad(()->concat(dim,x1,x2,x3), x1, x2, x3)
    end
end

@testset "dropout" begin
    x = parameter(randn(T,10,5))
    dropout(x, 0.5)
end

@testset "getindex" begin
    x = parameter(randn(T,10,5,4))
    checkgrad(()->x[2:7,:,1:3], x)
    x = parameter(randn(T,10,5))
    checkgrad(()->x[:,3:3], x)
end

@testset "linear" begin
    x = parameter(randn(T,10,5))
    f = Linear(T, 10, 7, init_b=Uniform(-1,1))
    checkgrad(()->f(x), x, f.W, f.b)
end

@testset "lookup" begin
    w = parameter(randn(T,10,15))
    x = Var(rand(1:15,10))
    checkgrad(()->lookup(w,x), w, x)
end

@testset "pack" begin
    x = parameter(randn(T,10,10))
    dims = [2, 5, 3]
    checkgrad(()->pack(x,[2,5,3],0), x)
end

@testset "reshape" begin
    x = parameter(randn(T,10,1,5))
    checkgrad(()->reshape(x,5,10), x)
    checkgrad(()->vec(x), x)
    checkgrad(()->dropdims(x), x)
end

@testset "softmax" begin
    x = parameter(rand(T,10,5))
    checkgrad(()->softmax(x), x)
    logsoftmax(x)
end

@testset "window1d" begin
    x = parameter(randn(T,10,10))
    ksize, padding, stride, dilation = 5, 2, 1, 1
    checkgrad(()->window1d(x,[3,7],ksize,padding,stride,dilation), x)
end

end


text:
# xvecs_read and xvecs_write
@testset "xvecs" begin
  fn = tempname()
  d, n, m, h = 32, Int(1e3), 4, 256
  X, C, B = generate_random_dataset(Float32, Int16, d, n, m, h)

  # fvecs
  fvecs_write(X, fn)
  X2 = fvecs_read(n, fn)
  @test all(X .== X2)
  rm(fn)

  Xint = convert(Matrix{Int32}, floor.(X.-0.5f0)*1000)
  ivecs_write(Xint, fn)
  Xint2 = ivecs_read(n, fn)
  @test all(Xint .== Xint2)
  rm(fn)
end

ext: jl


text:
struct ALObjective{T,O<:AbstractObjective} <: TO.AbstractObjective
    obj::O
    conset::ALConstraintSet{T}
    alcost::Vector{T}
    function ALObjective{T}(obj::AbstractObjective) where T
        new{T,typeof(obj)}(obj, ALConstraintSet{T}(), zeros(T, length(obj)))
    end
end
RD.state_dim(alobj::ALObjective, k::Integer) = RD.state_dim(alobj.obj, k)
RD.control_dim(alobj::ALObjective, k::Integer) = RD.state_dim(alobj.obj, k)
RD.dims(alobj::ALObjective) = RD.dims(alobj.obj)

function TO.cost(alobj::ALObjective, Z::SampledTrajectory)
    # Calculate unconstrained cost
    J = TO.cost(alobj.obj, Z)

    # Calculate constraints
    # TODO: update trajectory if Z is not the cached one
    evaluate_constraints!(alobj.conset, Z)

    # Calculate AL penalty
    alobj.alcost .= 0
    alcost(alobj.conset)
    J += sum(alobj.alcost)

    return J
end

function cost_expansion!(alobj::ALObjective, E::CostExpansion, Z)
    # Calculate expansion of original cost
    cost_expansion!(alobj.obj, E, Z)

    # Calculate constraint Jacobians
    constraint_jacobians!(alobj.conset, Z)

    # Calculate expansion of AL penalty
    algrad!(alobj.conset)
    alhess!(alobj.conset)

    # Add to existing expansion
    # each ALConstraint has local alias to ilqr.Efull
    # this is a hack to avoid an allocation for each constraint
    # due to type instability
    @assert E === alobj.conset.Eref[]
    add_alcost_expansion!(alobj.conset, E)

    return nothing
end


text:
module CategoricalAlgebra

using Reexport

include("FreeDiagrams.jl")
include("Limits.jl")
include("FinSets.jl")
include("Matrices.jl")
include("FinRelations.jl")
include("PredicatedSets.jl")
include("CSets.jl")
include("StructuredCospans.jl")
include("Squares.jl")

@reexport using .FreeDiagrams
@reexport using .Limits
@reexport using .CSets
@reexport using .StructuredCospans
@reexport using .Squares

end


text:
function create_exodus_database(file_name::ExoFileName)
    """
    TODO: clean up this method, I don't think it's used or tested
    For some reason #define statements are not picked up by ccall
    """
    @suppress begin
        exo_id = ccall((:ex_create_int, libexodus), ExoID,
                    (Base.Cstring, Int64, Ref{Int64}, Ref{Int64}, Int64),
                    file_name, EX_CLOBBER, cpu_word_size, IO_word_size,
                    version_number_2)
        exodus_error_check(exo_id, "create_exodus_database")
        return exo_id
    end
end

function copy_exodus_database(exo_id::ExoID, new_exo_id::ExoID)
    @suppress begin
        new_exo_id = ccall((:ex_copy, libexodus), ExoID,
                        (ExoID, ExoID),
                        exo_id, new_exo_id)
        exodus_error_check(new_exo_id, "copy_exodus_database")
    end
end

function close_exodus_database(exo_id::ExoID)
    @suppress begin
        error = ccall((:ex_close, libexodus), ExodusError, (ExoID,), exo_id)
        exodus_error_check(error, "close_exodus_database")
    end
end

function open_exodus_database(file_name::ExoFileName)
    # TODO: maybe add multiple methods for different EX_* options
    @suppress begin
        exo_id = ccall((:ex_open_int, libexodus), ExoID,
                    (Base.Cstring, Int64, Ref{Int64},
                        Ref{Int64}, Ref{Float64}, Int64),
                    file_name, EX_CLOBBER, cpu_word_size, IO_word_size,
                    version_number, version_number_2)
        exodus_error_check(exo_id, "open_exodus_database")
        return exo_id
    end
end


text:
using BinaryBuilder

name = "NodeJS_16"
version = v"16.15.0"

url_prefix = "https://nodejs.org/dist/v$version/node-v$version"
sources = [
    ArchiveSource("$(url_prefix)-linux-x64.tar.gz", "d1c1de461be10bfd9c70ebae47330fb1b4ab0a98ad730823fb1340e34993edee"; unpack_target = "x86_64-linux-gnu"),
    ArchiveSource("$(url_prefix)-linux-arm64.tar.gz", "2aa387e6a57ade663849efdc4fabf7431a38d975db98dcc79293840e6894d28b"; unpack_target = "aarch64-linux-gnu"),
    ArchiveSource("$(url_prefix)-linux-ppc64le.tar.gz", "625bf1f6cc2d608c51fc5b412ca162251871d14eb795cb006107d743c1da200c"; unpack_target = "powerpc64le-linux-gnu"),
    ArchiveSource("$(url_prefix)-linux-armv7l.tar.gz", "3b54c8f57a8ab211b5e969cdf6d20b3bcd7f30f7e0444e00c409f78b90486d30"; unpack_target = "arm-linux-gnueabihf"),

    ArchiveSource("$(url_prefix)-linux-x64.tar.gz", "d1c1de461be10bfd9c70ebae47330fb1b4ab0a98ad730823fb1340e34993edee"; unpack_target = "x86_64-linux-musl"),
    ArchiveSource("$(url_prefix)-linux-arm64.tar.gz", "2aa387e6a57ade663849efdc4fabf7431a38d975db98dcc79293840e6894d28b"; unpack_target = "aarch64-linux-musl"),
    ArchiveSource("$(url_prefix)-linux-armv7l.tar.gz", "3b54c8f57a8ab211b5e969cdf6d20b3bcd7f30f7e0444e00c409f78b90486d30"; unpack_target = "arm-linux-musleabihf"),

    ArchiveSource("$(url_prefix)-darwin-x64.tar.gz", "a6bb12bbf979d32137598e49d56d61bcddf8a8596c3442b44a9b3ace58dd4de8"; unpack_target = "x86_64-apple-darwin14"),
    ArchiveSource("$(url_prefix)-darwin-arm64.tar.gz", "ad8d8fc5330ef47788f509c2af398c8060bb59acbe914070d0df684cd2d8d39b"; unpack_target = "aarch64-apple-darwin20"),

    ArchiveSource("$(url_prefix)-win-x64.zip", "dbe04e92b264468f2e4911bc901ed5bfbec35e0b27b24f0d29eff4c25e428604"; unpack_target = "x86_64-w64-mingw32"),
    ArchiveSource("$(url_prefix)-win-x86.zip", "0d11a3844dad4ab679502495a4aa41041168a2caa81b8da9c7b5a14902c46986"; unpack_target = "i686-w64-mingw32"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/
cp -r ${target}/*/* ${prefix}/.
cd ${prefix}
if [[ "${target}" == *-mingw* ]]; then
    mkdir bin
    chmod +x {node.exe,npm,npm.cmd,npx,npx.cmd}
    mv {node.exe,npm,npm.cmd,npx,npx.cmd} bin
    mv node_modules bin
fi
install_license LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc="glibc", cxxstring_abi="cxx11"),
    Platform("aarch64", "linux"; libc="glibc", cxxstring_abi="cxx11"),
    Platform("powerpc64le", "linux"; libc="glibc", cxxstring_abi="cxx11"),
    Platform("armv7l", "linux"; libc="glibc", cxxstring_abi="cxx11"),

    Platform("x86_64", "linux"; libc="musl", cxxstring_abi="cxx11"),
    Platform("aarch64", "linux"; libc="musl", cxxstring_abi="cxx11"),
    Platform("armv7l", "linux"; libc="musl", cxxstring_abi="cxx11"),

    Platform("x86_64", "macos"),
    Platform("aarch64", "macos"),

    Platform("x86_64", "windows"),
    Platform("i686", "windows"),
]

# The products that we will ensure are always built
products = [
    ExecutableProduct("node", :node),
    FileProduct("bin/npm", :npm),
    FileProduct("bin/npx", :npx),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies, julia_compat = "1.6")


text:
## csgObject

abstract type csgObject <: Object end

export csgUnion
struct csgUnion <: csgObject
    objects :: Array{Object,1}

    function csgUnion(objects)
        if length(DeleteDuplicates(deleteat!(objects, isempty.(objects)))) == 0
            Empty()
        elseif length(DeleteDuplicates(deleteat!(objects, isempty.(objects)))) == 1
            objects[1]
        else
            new(DeleteDuplicates(deleteat!(objects, isempty.(objects))))
        end
    end
    function csgUnion(objects::Object...)
        return csgUnion([objects...])
    end
    function csgUnion(itr::Base.Generator)
        return csgUnion(collect(itr))
    end
end

export csgIntersection
struct csgIntersection <: csgObject
    objects :: Array{Object,1}

    function csgIntersection(objects)
        if Empty() ∈ objects
            Empty()
        elseif length(DeleteDuplicates(objects)) == 0
            Empty()
        elseif length(DeleteDuplicates(objects)) == 1
            objects[1]
        else
            new(DeleteDuplicates(objects))
        end
    end
    function csgIntersection(objects::Object...)
        return csgIntersection([objects...])
    end
    function csgIntersection(itr::Base.Generator)
        return csgIntersection(collect(itr))
    end
end

export csgMerge
struct csgMerge <: csgObject
    objects :: Array{Object,1}
    function csgMerge(objects)
        if length(DeleteDuplicates(deleteat!(objects, isempty.(objects)))) == 0
            Empty()
        elseif length(DeleteDuplicates(deleteat!(objects, isempty.(objects)))) == 1
            objects[1]
        else
            new(DeleteDuplicates(deleteat!(objects, isempty.(objects))))
        end
    end
    function csgMerge(objects::Object...)
        return csgMerge([objects...])
    end
    function csgMerge(itr::Base.Generator)
        return csgMerge(collect(itr))
    end
end

export csgDifference
struct csgDifference <: csgObject
    objects :: Array{Object,1}
    function csgDifference(objects)
        if length(objects) ≠ 2
            error("Too many objects.")
        elseif objects[1] == Empty()
            Empty()
        elseif objects[2] == Empty()
            objects[1]
        else
            new(objects)
        end
    end
    function csgDifference(object1::Object,object2::Object)
        return csgDifference([object1,object2])
    end
end

export csgClip
struct csgClip <: csgObject
    objects :: Array{Object,1}
    function csgClip(objects)
        if length(objects) ≠ 2
            error("Too many objects.")
        elseif objects[1] == Empty()
            Empty()
        elseif objects[2] == Empty()
            Empty()
        else
            new(objects)
        end
    end
    function csgClip(object1::Object,object2::Object)
        return csgClip([object1,object2])
    end
end

export csgBound
struct csgBound <: csgObject
    objects :: Array{Object,1}
    function csgBound(objects)
        if length(objects) ≠ 2
            error("Too many objects.")
        elseif objects[1] == Empty()
            Empty()
        elseif objects[2] == Empty()
            Empty()
        else
            new(objects)
        end
    end
    function csgBound(object1::Object,object2::Object)
        return csgBound([object1,object2])
    end
end

function povray_script(csg :: csgUnion; preindent = 0)
    script = "union{\n" * "  "^(preindent+1)
    for object in csg.objects
        script *= povray_script(object; preindent=preindent+1) * "\n" * "  "^(preindent+1)
    end
    script = script[1:end-2] * "}"
    return script
end

function povray_script(csg :: csgIntersection; preindent = 0)
    script = "intersection{\n" * "  "^(preindent+1)
    for object in csg.objects
        script *= povray_script(object; preindent=preindent+1) * "\n" * "  "^(preindent+1)
    end
    script = script[1:end-2] * "}"
    return script
end

function povray_script(csg :: csgMerge; preindent = 0)
    script = "merge{\n" * "  "^(preindent+1)
    for object in csg.objects
        script *= povray_script(object; preindent=preindent+1) * "\n" * "  "^(preindent+1)
    end
    script = script[1:end-2] * "}"
    return script
end

function povray_script(csg :: csgDifference; preindent = 0)
    script = "difference{\n" * "  "^(preindent+1)
    for object in csg.objects
        script *= povray_script(object; preindent=preindent+1) * "\n" * "  "^(preindent+1)
    end
    script = script[1:end-2] * "}"
    return script
end

function povray_script(csg :: csgClip; preindent = 0)
    return "object{"*povray_script(csg.objects[1])*"clipped_by{"*povray_script(csg.objects[2])*"}}"
end
function povray_script(csg :: csgBound; preindent = 0)
    return "object{"*povray_script(csg.objects[1])*"bounded_by{"*povray_script(csg.objects[2])*"}}"
end

max_stars_repo_name: hyrodium/JuliRay

max_issues_repo_name: hyrodium/JuliRay

max_forks_repo_name: hyrodium/JuliRay


text:
# ---
# title: 15. 3Sum
# id: problem15
# author: zhwang
# date: 2022-03-15
# difficulty: Medium
# categories: Array, Two Pointers
# link: <https://leetcode.com/problems/3sum/description/>
# hidden: true
# ---
#
# Given an array `nums` of _n_ integers, are there elements _a_ , _b_ , _c_ in
# `nums` such that _a_ \+ _b_ \+ _c_ = 0? Find all unique triplets in the array
# which gives the sum of zero.
#
# Notice that the solution set must not contain duplicate triplets.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [-1,0,1,2,-1,-4]
#     Output: [[-1,-1,2],[-1,0,1]]
#
#
# **Example 2:**
#
#
#
#     Input: nums = []
#     Output: []
#
#
# **Example 3:**
#
#
#
#     Input: nums = [0]
#     Output: []
#
#
#
#
# **Constraints:**
#
#   * `0 <= nums.length <= 3000`
#   * `-105 <= nums[i] <= 105`
#
#
## @lc code=start
using LeetCode

three_sum(nums::Vector{Int}) = three_sum!(copy(nums))
function three_sum!(nums::Vector{Int})::Vector{NTuple{3, Int}}
    sort!(nums)
    n, res, pre = length(nums), NTuple{3, Int}[], NaN
    for (i, first) in enumerate(nums)
        first == pre && continue
        pre, left, right = first, i + 1, n
        while left < right
            total = nums[left] + nums[right]
            if total < -first
                left += 1
            elseif total > -first
                right -= 1
            else
                push!(res, (first, nums[left], nums[right]))
                pos = findfirst(!=(nums[left]), @view(nums[left+1:right]))
                isnothing(pos) ? break : left += pos
                pos = findfirst(!=(nums[right]), @view(nums[right-1:-1:left]))
                isnothing(pos) ? break : right -= pos
            end
        end
    end
    res
end

## @lc code=end

## @lc test=start
@testset "15.3sum.jl" begin
    @test three_sum([-1,0,1,2,-1,-4]) == [(-1,-1,2),(-1,0,1)]
    @test three_sum([0]) == three_sum!(Int[]) == NTuple{3, Int}[]
    @test three_sum([1, 1, -2]) == [(-2, 1, 1)]
    @test three_sum([0, 0, 0, 0]) == [(0, 0, 0)]
end
## @lc test=end

## @lc info=start
# link: [solution to 15]()
## @lc info=end


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module ClausenNumbers
using Nemo, PrimesIterator, Products, NumberTheory

export ModuleClausenNumbers
export ClausenNumber, ClausenNumberList
export V002445, L002445, V027642

"""

* ClausenNumber, ClausenNumberList, V002445, L002445, V027642
"""
const ModuleClausenNumbers = ""

"""

Return the Clausen number ``C_n`` which is the denominator of the Bernoulli number ``B_{2n}``.
"""
function ClausenNumber(n::Int)
    n == 0 && return fmpz(1)
    m = [d + 1 for d ∈ Divisors(2n)]
    ∏([q for q ∈ m if isPrime(q)])
end

"""

Return the list of length len of Clausen numbers which are the denominators of the Bernoulli numbers ``B_{2n}``.
"""
function ClausenNumberList(len::Int)
    len ≤ 0 && return fmpz[]

    A = fill(fmpz(2), len)
    A[1] = 1
    m = len - 1
    m == 0 && return A

    for p ∈ Primes(3, 2m + 1)
        r = div(p - 1, 2)
        for k ∈ range(r, step = r, length = div(m, r))
            A[k+1] *= p
        end
    end
    A
end

"""

Return the Clausen number ``C(n)`` which is the denominator of the Bernoulli number ``B_{2n}``.
"""
V002445(n::Int) = ClausenNumber(n)

"""

Return the list of length len of Clausen numbers which are the denominators of the Bernoulli numbers ``B_{2n}``.
"""
L002445(len::Int) = ClausenNumberList(len)

"""

Return the denominator of Bernoulli number ``B_n``.
"""
function V027642(n::Int)
    isEven(n) && return ClausenNumber(div(n, 2))
    n == 1 && return fmpz(2)
    return fmpz(1)
end

#START-TEST-########################################################

using Test, SeqUtils, SeqTests

function test()
∈
    @testset "Clausen" begin
        C = ClausenNumberList(800)
        @test C[125] == 30
        @test C[781] == 32695402455500348373810
        @test C[794] == 6

        @test isa(C[781], fmpz)
        @test isa(ClausenNumber(10), fmpz)
    end

    if data_installed()
        SeqTest([L002445], 'L')
        SeqTest(V002445, 'V', 0)
        SeqTest(V027642, 'V', 0)
    end
end

function demo()
    println("\nReturn the Clausen number C(n) which is the denominator of the Bernoulli number B(2n).")
    for n ∈ 0:10
        println(n, " ↦ ", V002445(n))
    end

    println("\nReturn the denominator of Bernoulli number B_n.")
    for n ∈ 0:10
        println(n, " ↦ ", V027642(n))
    end

    println()
    for n ∈ 0:6
        print(n, " ↦ ")
        Println(ClausenNumberList(n))
    end

    println()
    Println(ClausenNumberList(9))
    println()
end

"""

# ClausenNumberList is ~ 16 times faster than a list of ClausenNumbers!
ClausenNumberList(10000) ::
    0.016498 seconds (141.08 k allocations: 2.252 MiB)
for n ∈ 0:10000 ClausenNumber(n) end ::
    0.737140 seconds (2.40 M allocations: 66.330 MiB, 53.08% gc time)
"""
function perf()
    GC.gc()
    @time ClausenNumberList(10000)
    @time (for n ∈ 0:10000
        ClausenNumber(n)
    end)
end

function main()
    test()
    demo()
    perf()
end

main()

end # module

#=

0 ↦ fmpz[]
1 ↦ fmpz[3]
2 ↦ fmpz[3, 30]
3 ↦ fmpz[3, 30, 42]
4 ↦ fmpz[3, 30, 42, 30]
5 ↦ fmpz[3, 30, 42, 30, 66]
6 ↦ fmpz[3, 30, 42, 30, 66, 2730]

0 ↦ 1
1 ↦ 2
2 ↦ 6
3 ↦ 1
4 ↦ 30
5 ↦ 1
6 ↦ 42
7 ↦ 1
8 ↦ 30
9 ↦ 1
10 ↦ 66

  0.009131 seconds (141.08 k allocations: 2.252 MiB)
  0.633522 seconds (2.40 M allocations: 66.330 MiB, 57.26% gc time)
=#

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
const modules = Array{Symbol, 1}()

################################################################################
## INTEGER VERSION.
## a thin shim over the binary version, but shows how effective method rewrite
## can be.

function inject!(f::Expr, expr)
  unshift!(f.args[2].args, expr)
end

function set_output_type!(f::Expr, output_type)
  #if the function already has an expression type
  if f.args[1].head == :(::)
    f.args[1].args[2] = output_type
  else
    call_list = f.args[1]
    f.args[1] = Expr(:(::))
    push!(f.args[1].args, call_list)
    push!(f.args[1].args, output_type)
  end
end

function substitute_wire_inputs!(f::Expr)
  #issue a list of places where we need extra wire macros.
  wire_macro_list = Dict{Symbol,Expr}()

  for argument in f.args[1].args[1].args[2:end]
    if isa(argument, Symbol)
      continue
    end

    #a bare wire statement.
    if argument.args[2] == :Wire
      #reset it to Unsigned.
      argument.args[2] = :Unsigned
    elseif argument.args[2] == :SingleWire
      wire_macro_list[argument.args[1]] = :(0:0v)
      argument.args[2] = :Unsigned
    elseif isa(argument.args[2], Expr) &&
      (argument.args[2].head == :curly) &&
      (argument.args[2].args[1] == :Wire) &&
      (argument.args[2].args[2].head == :(:))

      #save the needed symbol and wire range parameter.
      wire_macro_list[argument.args[1]] = argument.args[2].args[2]

      #then reset it to Unsigned.
      argument.args[2] = :Unsigned
    end
  end

  wire_macro_list
end

function inject_input_macro!(f::Expr, identifier, structure)
  inject!(f, :(@input $identifier $structure))
end

function integer_translate(f::Expr)
  f_integer = copy(f)

  set_output_type!(f_integer, :Unsigned)

  input_wire_list = substitute_wire_inputs!(f_integer)

  for key in keys(input_wire_list)
    inject_input_macro!(f_integer, key, input_wire_list[key])
  end

  inject!(f_integer, :(Verilog.@verimode :integermode))

  f_integer
end

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

include("module_rewrite.jl")

################################################################################
# strip module names if it's not in module mode.
function strip_nonmodule!(f::Expr)
  #generate a temporary block
  stripped_block = []
  for argument in f.args[2].args
    if isa(argument, Expr) &&
      (argument.head == :macrocall) &&
      (argument.args[1] == Symbol("@suffix"))
      continue
    end
    push!(stripped_block, argument)
  end
  #reset the block
  f.args[2].args = stripped_block
end

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

macro verilog(f)
  #first check to see if it's a function.
  if f.head == :function
    #make three copies of the function.

    #make the module version.
    f_module = module_translate(f)
    f_wiretext = wiretext_translate(f)

    #strip module instructions
    strip_nonmodule!(f)

    #make version that substitutes all wires with integers.
    f_integer = integer_translate(f)
    #make a version that substitutes wires with a fused integer.
    make_fused = false

    if (make_fused)
      f_fusedinteger = fused_integer_translate(f)
    end

    inject!(f, :(Verilog.@verimode :wiremode))

    esc(quote
      #release all three forms of the function.
      Base.@__doc__ $f
      $f_integer
      $f_module
      $f_wiretext
    end)
  else
    throw(ArgumentError("@verilog macro must be run on a function."))
  end
end

export @verilog


text:
# to_expr
to_expr(vi::MOI.VariableIndex) = :(x[$(vi.value)])

function to_expr(f::SAF)
    f = MOIU.canonical(f)
    if isempty(f.terms)
        return f.constant
    else
        linear_term_exprs = map(f.terms) do term
            :($(term.coefficient) * x[$(term.variable.value)])
        end
        expr = :(+($(linear_term_exprs...)))
        if !iszero(f.constant)
            push!(expr.args, f.constant)
        end
        return expr
    end
end

function to_expr(f::SQF)
    f = MOIU.canonical(f)
    linear_term_exprs = map(f.affine_terms) do term
        i = term.variable.value
        :($(term.coefficient) * x[$i])
    end
    quadratic_term_exprs = map(f.quadratic_terms) do term
        i = term.variable_1.value
        j = term.variable_2.value
        if i == j
            :($(term.coefficient / 2) * x[$i] * x[$j])
        else
            :($(term.coefficient) * x[$i] * x[$j])
        end
    end
    expr = :(+($(linear_term_exprs...), $(quadratic_term_exprs...)))
    if !iszero(f.constant)
        push!(expr.args, f.constant)
    end
    return expr
end

# check_variable_indices
function check_variable_indices(model::Optimizer, index::VI)
    @assert 1 <= index.value <= length(model.inner.variable_info)
end

function check_variable_indices(model::Optimizer, f::SAF)
    for term in f.terms
        check_variable_indices(model, term.variable)
    end
end

function check_variable_indices(model::Optimizer, f::SQF)
    for term in f.affine_terms
        check_variable_indices(model, term.variable)
    end
    for term in f.quadratic_terms
        check_variable_indices(model, term.variable_1)
        check_variable_indices(model, term.variable_2)
    end
end

function find_variable_info(model::Optimizer, vi::VI)
    check_variable_indices(model, vi)
    model.inner.variable_info[vi.value]
end

function check_constraint_indices(model::Optimizer, index::CI{MOI.VariableIndex})
    @assert 1 <= index.value <= length(model.inner.variable_info)
end

function check_constraint_indices(model::Optimizer, index::CI)
    @assert 1 <= index.value <= length(model.inner.constraint_info)
end


text:
abstract Language

# Languages currently supported by Snowball
abstract DanishLanguage     <: Language
abstract DutchLanguage      <: Language
abstract EnglishLanguage    <: Language
abstract FinnishLanguage    <: Language
abstract FrenchLanguage     <: Language
abstract GermanLanguage     <: Language
abstract HungarianLanguage  <: Language
abstract ItalianLanguage    <: Language
abstract NorwegianLanguage  <: Language
abstract PortugueseLanguage <: Language
abstract RomanianLanguage   <: Language
abstract RussianLanguage    <: Language
abstract SpanishLanguage    <: Language
abstract SwedishLanguage    <: Language
abstract TurkishLanguage    <: Language

# These are ISO 639-2T alpha-3 and ISO 639-3 codes
isocode(lang::Type{DanishLanguage})     = "dan"
isocode(lang::Type{DutchLanguage})      = "nld"
isocode(lang::Type{EnglishLanguage})    = "eng"
isocode(lang::Type{FinnishLanguage})    = "fin"
isocode(lang::Type{FrenchLanguage})     = "fra"
isocode(lang::Type{GermanLanguage})     = "deu"
isocode(lang::Type{HungarianLanguage})  = "hun"
isocode(lang::Type{ItalianLanguage})    = "ita"
isocode(lang::Type{NorwegianLanguage})  = "nor"
isocode(lang::Type{PortugueseLanguage}) = "por"
isocode(lang::Type{RomanianLanguage})   = "ron"
isocode(lang::Type{RussianLanguage})    = "rus"
isocode(lang::Type{SpanishLanguage})    = "spa"
isocode(lang::Type{SwedishLanguage})    = "swe"
isocode(lang::Type{TurkishLanguage})    = "tur"

name(lang::Type{DanishLanguage})     = "danish"
name(lang::Type{DutchLanguage})      = "dutch"
name(lang::Type{EnglishLanguage})    = "english"
name(lang::Type{FinnishLanguage})    = "finnish"
name(lang::Type{FrenchLanguage})     = "french"
name(lang::Type{GermanLanguage})     = "german"
name(lang::Type{HungarianLanguage})  = "hungarian"
name(lang::Type{ItalianLanguage})    = "italian"
name(lang::Type{NorwegianLanguage})  = "norwegian"
name(lang::Type{PortugueseLanguage}) = "portuguese"
name(lang::Type{RomanianLanguage})   = "romanian"
name(lang::Type{RussianLanguage})    = "russian"
name(lang::Type{SpanishLanguage})    = "spanish"
name(lang::Type{SwedishLanguage})    = "swedish"
name(lang::Type{TurkishLanguage})    = "turkish"


text:
isdefined(LibGit2, :discover) ||
@eval LibGit2 function discover(
    start_path::AbstractString = pwd();
    ceiling::Union{AbstractString,Vector} = "",
    across_fs::Bool = false,
)
    sep = @static is_windows() ? ";" : ":"
    ceil = ceiling isa AbstractString ? ceiling :
        join(convert(Vector{String}, ceiling), sep)
    buf_ref = Ref(Buffer())
    @check ccall((:git_repository_discover, :libgit2), Cint,
                 (Ptr{Buffer}, Cstring, Cint, Cstring),
                 buf_ref, start_path, across_fs, ceil)
    buf = buf_ref[]
    str = unsafe_string(buf.ptr, buf.size)
    free(buf_ref)
    return str
end


text:
"""
$(SIGNATURES)
    Get true number of constraints, accounting for minimum time and infeasible start constraints
    -p: total number of stage constraints (state and control)
    -pI: number of inequality stage constraints (state and control)
    -pE: number of equality stage constraints (stage and control)
"""
function get_num_constraints(solver::Solver,set::Symbol=:all)
    if solver.state.constrained
        if set == :custom
            pI = solver.obj.pI_custom
            pE = solver.obj.pE_custom
            p = pI + pE
        else
            if solver.obj isa ConstrainedObjective
                p = solver.obj.p
                pI = solver.obj.pI
                pE = p - pI
            else
                p,pI,pE = 0,0,0
            end
            if is_min_time(solver)
                pI += 2
                pE += 1
            end
            solver.state.infeasible ? pE += solver.model.n : nothing
            p = pI + pE
        end
        return p, pI, pE
    else
        return 0,0,0
    end
end

function get_num_terminal_constraints(solver::Solver,set::Symbol=:all)
    if solver.state.constrained
        if set == :custom
            pI_N = solver.obj.pI_N_custom
            pE_N = solver.obj.pE_N_custom
            p_N = pI_N + pE_N
        else
            if solver.obj isa ConstrainedObjective
                p_N = solver.obj.p_N
                pI_N = solver.obj.pI_N

                pE_N = p_N - pI_N
            else
                p_N,pI_N,pE_N = 0,0,0
            end
        end
        return p_N,pI_N,pE_N
    else
        return 0,0,0
    end
end

"""
$(SIGNATURES)
Evalutes all inequality and equality constraints (in place) for the current
state and control trajectories
    see: A Novel Augmented Lagrangian Approach for Inequalities and Convergent
    Any-Time Non-Central Updates (Toussaint)
"""
function update_constraints!(res::ConstrainedIterResults, solver::Solver, X=res.X, U=res.U)::Nothing
    n,m,N = get_sizes(solver)
    p,pI,pE = get_num_constraints(solver)
    p_N,pI_N,pE_N = get_num_terminal_constraints(solver)
    m̄,mm = get_num_controls(solver)
    n̄,nn = get_num_states(solver)

    # c_fun = solver.c_fun
    c_fun = constraint_function(solver)

    for k = 1:N-1
        # Update constraints
        c_fun(res.C[k], X[k], U[k])

        # Minimum time special case
        if solver.state.minimum_time
            if k == 1
                res.C[k][p] = 0.0
            end
        end

        # Get active constraint set
        get_active_set!(res,solver,p,pI,k)

        # Update penality-indicator matrices based on active set
        res.Iμ[k] = Diagonal(res.active_set[k].*res.μ[k])
    end

    # Terminal constraint
    c_fun(res.C[N],X[N][1:n])
    get_active_set!(res,solver,p_N,pI_N,N)

    res.Iμ[N] = Diagonal(res.active_set[N].*res.μ[N])

    return nothing
end

function update_constraints!(res::UnconstrainedIterResults, solver::Solver, X=res.X, U=res.U)::Nothing
    return nothing
end

"""
$(SIGNATURES)
    Determine active set for inequality constraints
"""
function get_active_set!(results::ConstrainedIterResults,solver::Solver,p::Int,pI::Int,k::Int)
    # Inequality constraints
    for j = 1:pI
        if active_set_criteria(solver,results.C[k][j], results.λ[k][j], results.μ[k][j])
            results.active_set[k][j] = true
        else
            results.active_set[k][j] = false
        end
    end
    # Equality constraints
    for j = pI+1:p
        results.active_set[k][j] = true
    end
    return nothing
end

function active_set_criteria(solver::Solver,c,λ,μ)::Bool
    if solver.opts.al_type == :default
        c > -solver.opts.active_constraint_tolerance || λ > 0.0
    elseif solver.opts.al_type == :algencan
        λ > 0
        c > -solver.opts.active_constraint_tolerance || λ > 0.0
        λ + μ*c > 0
    else
        error("al_type $(solver.opts.al_type) not recognized")
    end
end

"""
$(SIGNATURES)
    Count the number of constraints of each type from an objective
"""
function count_constraints(obj::ConstrainedObjective, constraints::Symbol=:all)
    n,m = get_sizes(obj)
    p = obj.p # number of constraints
    pI = obj.pI # number of inequality and equality constraints
    pE = p-pI # number of equality constraints

    pI_c = obj.pI_custom
    pE_c = obj.pE_custom

    # Terminal constraints
    p_N = obj.p_N
    pI_N = obj.pI_N
    pE_N = p_N - pI_N

    pI_N_c = obj.pI_N_custom
    pE_N_c = obj.pE_N_custom

    if constraints == :all
        return (pI, pI_c, pI_N, pI_N_c), (pE, pE_c, pE_N, pE_N_c)
    elseif constraints == :custom
        return (pI_c, pI_N_c), (pE_c, pE_N_c)
    elseif constraints == :total
        return (pI, pI_N), (pE, pE_N)
    end

end

"""
$(SIGNATURES)
    Generate the Jacobian of a general nonlinear constraint function
        -constraint function must be inplace
        -automatic differentition via ForwardDiff.jl
"""
function generate_general_constraint_jacobian(c::Function,p::Int,n::Int64,m::Int64)::Function
    c_aug! = f_augmented!(c,n,m)
    J = zeros(p,n+m)
    S = zeros(n+m)
    cdot = zeros(p)
    F(J,cdot,S) = ForwardDiff.jacobian!(J,c_aug!,cdot,S)

    function c_jacobian(cx,cu,x,u)
        S[1:n] = x[1:n]
        S[n+1:n+m] = u[1:m]
        F(J,cdot,S)
        cx[1:p,1:n] = J[1:p,1:n]
        cu[1:p,1:m] = J[1:p,n+1:n+m]
    end
    return c_jacobian
end

function generate_general_constraint_jacobian(c::Function,p::Int,n::Int64)::Function
    J_N = zeros(p,n)
    xdot = zeros(p)
    F_N(J_N,xdot,x) = ForwardDiff.jacobian!(J_N,c,xdot,x) # NOTE: terminal constraints can only be dependent on state x_N
    function c_jacobian(cx,x)
        F_N(J_N,xdot,x)
        cx .= J_N
    end
    return c_jacobian
end

"""
$(SIGNATURES)
Generate the constraints function C(x,u) and a function to compute the jacobians
Cx, Cu = Jc(x,u) from a `ConstrainedObjective` type. Automatically stacks inequality
and equality constraints and takes jacobians of custom functions with `ForwardDiff`.
Stacks the constraints as follows:

[upper control inequalities
 (√dt upper bound)
 lower control inequalities
 (√dt lower bound)
 upper state inequalities
 lower state inequalities
 general inequalities
 general equalities
 (control equalities for infeasible start)
 (dt equality)]
"""
function generate_constraint_functions(obj::ConstrainedObjective; max_dt::Float64=1.0, min_dt::Float64=1e-2)
    n,m = get_sizes(obj)

    # Key: I=> inequality,   E=> equality
    #     _c=> custom   (lack)=> box constraint
    #     _N=> terminal (lack)=> stage

    min_time = obj.tf == 0

    p = obj.p # number of constraints
    pI, pI_c, pI_N, pI_N_c = obj.pI, obj.pI_custom, obj.pI_N, obj.pI_N_custom
    pE, pE_c, pE_N, pE_N_c = p-obj.pI, obj.pE_custom, obj.p_N - obj.pI_N, obj.pE_N_custom
    m̄ = m
    n̄ = n
    if min_time
         m̄ += 1
         n̄ += 1
    end
    labels = String[]

    # Append on min time bounds
    u_max = obj.u_max
    u_min = obj.u_min
    if min_time
        u_max = [u_max; sqrt(max_dt)]
        u_min = [u_min; sqrt(min_dt)]
    end

    # Mask for active (state|control) constraints
    u_min_active = isfinite.(u_min)
    u_max_active = isfinite.(u_max)
    x_min_active = isfinite.(obj.x_min)
    x_max_active = isfinite.(obj.x_max)

    # Inequality on control
    pI_u_max = count(u_max_active)
    pI_u_min = count(u_min_active)
    pI_u = pI_u_max + pI_u_min
    function c_control_limits!(c,x,u)
        c[1:pI_u_max] = (u[1:m̄] - u_max)[u_max_active]
        c[pI_u_max+1:pI_u_max+pI_u_min] = (u_min - u[1:m̄])[u_min_active]
    end

    lbl_u_min = ["control (lower bound)" for i = 1:pI_u_min]
    lbl_u_max = ["control (upper bound)" for i = 1:pI_u_max]
    if min_time
        lbl_u_min[end] = "* √dt (lower bound)"
        lbl_u_max[end] = "* √dt (upper bound)"
    end

    # Inequality on state
    pI_x_max = count(x_max_active)
    pI_x_min = count(x_min_active)
    pI_x = pI_x_max + pI_x_min
    function c_state_limits!(c,x,u)
        c[1:pI_x_max] = (x[1:n] - obj.x_max )[x_max_active]
        c[pI_x_max+1:pI_x_max+pI_x_min] = (obj.x_min - x[1:n])[x_min_active]
    end
    lbl_x_max = ["state (upper bound)" for i = 1:pI_x_max]
    lbl_x_min = ["state (lower bound)" for i = 1:pI_x_min]

    # Update pI
    pI = pI_x + pI_u + pI_c

    # Form inequality constraint
    function cI!(c,x,u)
        c_control_limits!(view(c,1:pI_u),x,u)
        c_state_limits!(view(c,(1:pI_x).+pI_u),x,u)
        if pI_c > 0
            obj.cI(view(c,(1:pI_c).+pI_u.+pI_x),x,u)
        end
    end
    lbl_cI = ["custom inequality" for i = 1:pI_c]
    lbl_cE = ["custom equality" for i = 1:pE_c]

    function cI!(c,x)
        if obj.use_xf_equality_constraint
            nothing
        else
            if pI_x_max > 0
                c[1:pI_x_max] = (x[1:n] - obj.x_max )[x_max_active]
            end
            if pI_x_min > 0
                c[pI_x_max .+ (1:pI_x_min)] = (obj.x_min - x[1:n])[x_min_active]
            end
            if obj.pI_N_custom > 0
                c[pI_x .+ (1:obj.pI_N_custom)] = obj.cI_N(c,x[1:n])
            end
        end
    end

    # Construct labels
    c_labels = [lbl_u_max; lbl_u_min; lbl_x_max; lbl_x_min; lbl_cI; lbl_cE]

    # Augment functions together
    function c_function!(c,x,u)::Nothing
        infeasible = length(u) != m̄
        if pI > 0
            cI!(view(c,1:pI),x,u[1:m̄])
        end
        if pE_c > 0
            obj.cE(view(c,(1:pE_c).+pI),x,u[1:m])
        end
        if infeasible
            c[pI.+pE_c.+(1:n)] = u[m̄.+(1:n)]
        end
        if min_time
            c[pI+pE_c+(n*infeasible)+1] = u[m̄] - x[n̄]
        end
        return nothing
    end

    # Terminal Constraint
    """
    [xN-xf]
    or
    [x - x_max
     x_min - x
     cI_N_custom
     cI_E_custom]
    """
    function c_function!(c,x)
        if obj.use_xf_equality_constraint
            c[1:n] = x - obj.xf
        else
            if pI_x_max > 0
                c[1:pI_x_max] = (x[1:n] - obj.x_max)[x_max_active]
            end
            if pI_x_min > 0
                c[pI_x_max .+ (1:pI_x_min)] = (obj.x_min - x[1:n])[x_min_active]
            end
            if obj.pI_N_custom > 0
                obj.cI_N(view(c,pI_x .+ (1:obj.pI_N_custom)),x[1:n])
            end
            if obj.pE_N_custom > 0
                obj.cE_N(view(c,(obj.pI_N_custom + pI_x) .+ (1:pE_N)),x[1:n])
            end
        end
    end

    ### Jacobians ###
    # Declare known Jacobians
    In = Matrix(I,n,n)
    cx_control_limits = zeros(pI_u,n)
    cx_state_limits = zeros(pI_x,n)
    cx_state_limits[1:pI_x_max, :] = In[x_max_active,:]
    cx_state_limits[pI_x_max+1:end,:] = -In[x_min_active,:]

    Im = Matrix(I,m̄,m̄)
    cu_control_limits = zeros(pI_u,m̄)
    cu_control_limits[1:pI_u_max,:] = Im[u_max_active,:]
    cu_control_limits[pI_u_max+1:end,:] = -Im[u_min_active,:]
    cu_state_limits = zeros(pI_x,m̄)

    if pI_c > 0
        cI_custom_jacobian! = generate_general_constraint_jacobian(obj.cI, pI_c, n, m)
    end
    if pE_c > 0
        cE_custom_jacobian! = generate_general_constraint_jacobian(obj.cE, pE_c, n, m)
    end

    cx_infeasible = zeros(n,n)
    cu_infeasible = In

    function c_jacobian!(cx::AbstractMatrix, cu::AbstractMatrix, x::AbstractArray,u::AbstractArray)
        infeasible = length(u) != m̄
        let m = m̄
            cx[1:pI_u, 1:n] = cx_control_limits
            cx[(1:pI_x).+pI_u, 1:n] = cx_state_limits

            cu[1:pI_u, 1:m] = cu_control_limits
            cu[(1:pI_x).+pI_u, 1:m] = cu_state_limits
        end

        if pI_c > 0
            cI_custom_jacobian!(view(cx,pI_x+pI_u+1:pI_x+pI_u+pI_c,1:n), view(cu,pI_x+pI_u+1:pI_x+pI_u+pI_c,1:m), x, u[1:m])
        end
        if pE_c > 0
            cE_custom_jacobian!(view(cx,pI_x+pI_u+pI_c+1:pI_x+pI_u+pI_c+pE_c,1:n), view(cu,pI_x+pI_u+pI_c+1:pI_x+pI_u+pI_c+pE_c,1:m), x, u[1:m])
        end

        if infeasible
            cx[pI+pE_c+1:pI+pE_c+n,1:n] = cx_infeasible
            cu[pI+pE_c+1:pI+pE_c+n,m̄+1:m̄+n] = cu_infeasible
        end
        if min_time
            cx[pI+pE_c+(n*infeasible)+1,n̄] = -1.0
            cu[pI+pE_c+(n*infeasible)+1,m̄] = 1.0
        end
    end

    # Terminal Constraint
    """
    [In
     -In
    JI_N_custom;
     In;
     JI_E_custom]
    """

    if obj.pI_N_custom > 0
        cI_N_custom_jacobian! = generate_general_constraint_jacobian(obj.cI_N, obj.pI_N_custom, n)
    end
    if obj.pE_N_custom > 0
        cE_N_custom_jacobian! = generate_general_constraint_jacobian(obj.cE_N, obj.pE_N_custom, n)
    end

    function c_jacobian!(j::AbstractArray,x::AbstractArray)
        if obj.use_xf_equality_constraint
            j[1:n,1:n] = In
        else
            if pI_x_max > 0
                j[1:pI_x_max,1:n] = In[x_max_active,:]
            end
            if pI_x_min > 0
                j[pI_x_max .+ (1:pI_x_min),1:n] = -In[x_min_active,:]
            end
            if obj.pI_N_custom > 0
                cI_N_custom_jacobian!(view(j,pI_x .+(1:obj.pI_N_custom),1:n),x[1:n])
            end
            if obj.pE_N_custom > 0
                cE_N_custom_jacobian!(view(j,(obj.pI_N_custom + pI_x) .+ (1:pE_N),1:n),x[1:n])
            end
        end
    end

    cE!(c,x,u) = obj.cE(c,x,u)

    function cE!(c,x)
        if obj.use_xf_equality_constraint
            c[1:n] = x - obj.xf
        elseif obj.pE_N_custom > 0
            c[1:pE_N] = obj.cE_N(view(c,1:pE_N),x[1:n])
        end
    end

    return c_function!, c_jacobian!, c_labels, cI!, cE!
end

generate_constraint_functions(obj::UnconstrainedObjective; max_dt::Float64=1.0,min_dt=1.0e-2) = null_constraint, null_constraint_jacobian, String[], null_constraint, null_constraint

"""
$(SIGNATURES)
Compute the maximum constraint violation. Inactive inequality constraints are
not counted (masked by the Iμ matrix).
"""
function constraint_ℓ2_norm(results::ConstrainedIterResults)
    return norm(map((x)->x.>0, results.Iμ) .* results.C)
end

function contraint_ℓ2_norm(results::UnconstrainedIterResults)
    return 0.0
end

"""
$(SIGNATURES)
Compute the ℓ2-norm over entire constraint trajectory. Inactive inequality constraints are
not counted (masked by the Iμ matrix).
"""
function max_violation(results::ConstrainedIterResults)
    return maximum(norm.(map((x)->x.>0, results.Iμ) .* results.C, Inf))
end

function max_violation(results::UnconstrainedIterResults)
    return 0.0
end

## Simple constraint primitives
"""
$(SIGNATURES)
Circle constraint function (c ⩽ 0, negative is satisfying constraint)
"""
function circle_constraint(x,x0,y0,r)
	return -((x[1]-x0)^2 + (x[2]-y0)^2  - r^2)
end

circle_constraint(x,c,r) = circle_constraint(x,c[1],c[2],r)

"""
$(SIGNATURES)
Sphere constraint function (c ⩽ 0, negative is satisfying constraint)
"""
function sphere_constraint(x,x0,y0,z0,r)
	return -((x[1]-x0)^2 + (x[2]-y0)^2 + (x[3]-z0)^2  - r^2)
end

function sphere_constraint(x,x0,r)
	return -((x[1]-x0[1])^2 + (x[2]-x0[2])^2 + (x[3]-x0[3])^2-r^2)
end


text:
# Fonctions permettant de calculer des modèles types SIR/SEIR

export SIR ;
"""
Function SIR

"""
function SIR(β::float64,γ::float64,S0::float64)

    S = zeros(Float64,ParamSIR.NB_DT) ;
    I = zeros(Float64,ParamSIR.NB_DT) ;
    R = zeros(Float64,ParamSIR.NB_DT) ;
    R0 = β/γ ;
    println("Compute SEIR model with R_0 = ", R0) ;

    # Initial conditions
    S[1] = S0 ;
    I[1] = 1 - S0 ;

    for i=2:ParamSIR.NB_DT
        S[i+1] = (1 - β*I[i])*S[i] ;
        I[i+1] = (1 + β*S[i] - γ)*I[i] ;
        R[i+1] = γ*I[i]+R[i] ;
    end

    return(S,I,R)

end

export SEIR ;

function SEIR(α::float64,β::float64,γ::float64,S0::float64)

    S = zeros(Float64,ParamSEIR.NB_DT) ;
    E = zeros(Float64,ParamSEIR.NB_DT) ;
    I = zeros(Float64,ParamSEIR.NB_DT) ;
    R = zeros(Float64,ParamSEIR.NB_DT) ;

    R0 = β/γ ;
    println("Compute SEIR model with R_0 = ", R0) ;

    # Initial conditions
    S[1] = S0 ;
    I[1] = 1-S0 ;

    for i=2:ParamSEIR.NB_DT
        S[i+1] = (1 - β*I[i])*S[i] ;
        E[i+1] = β*S[i]*I[i] + (1-α)*E[i] ;
        I[i+1] = α*E[i] + (1-γ)*I[i] ;
        R[i+1] = γ*I[i]+R[i] ;
    end

    return()

end


text:
# Script to run a test case for understanding how M scales with β #

# Add paths
loc_gm_package = "../src/"
if !(loc_gm_package in LOAD_PATH)
    push!(LOAD_PATH, "../src/")
end

# Required packages
using LinearAlgebra
using Convex, SCS, Mosek, MosekTools

using GraphicalModelLearning
using LightGraphs

# Function definitions
using Printf, SparseArrays
include("beta_scaling.jl")
include("active_learning.jl")

using DelimitedFiles

# parameters
FLAG_create_struct_gm = true
FLAG_regular_random_gm = false
FLAG_lattice_gm = true
FLAG_weak_impurity = true
sampling_regime = M_regime()
learning_algo = NLP()

if FLAG_lattice_gm
    # text file to save adjacency matrix in
    file_adj_matrix_gm = "adj_matrix_ferro_lattice_gm_E_MEE2.txt"
    # name of picture to save graphical model in
    file_plot_gm = "ferro_lattice_gm_E_MEE2.eps"
    # File to save final results
    file_M_opt_gm = "M_opt_Ferro_Lattice_E_MEE2.txt"
elseif FLAG_regular_random_gm
    # text file to save adjacency matrix in
    file_adj_matrix_gm = "adj_matrix_spin_glass_gm.txt"
    # name of picture to save graphical model in
    file_plot_gm = "spin_glass_gm3.eps"
    # File to save final results
    file_M_opt_gm = "M_opt_M_Regime_Regular_Random_B.txt"
end

# Simple test
N = 16
d = 4
α = 0.4
β_array = [1.8,3.2]

# Create and plot the initial graphical graphical model
β = copy(β_array[1])
adj_matrix, struct_adj_matrix = ferro_lattice(N,α,β,FLAG_weak_impurity)
m = n = Int(sqrt(N))
open(file_adj_matrix_gm, "w") do io
    writedlm(io, adj_matrix)
end;

## Start the complexity studies
τ = α/2
L = 45
M_g = 25000
M_factor = 0.05

M_opt = Array{Int64,1}(undef,length(β_array))

for i = 1:length(β_array)
    # define β
    β = copy(β_array[i])
    println(β)

    if i==1
        M_g = 25000
    end
    println(M_g)

    # Create the adjacency matrix
    adj_matrix, struct_adj_matrix = ferro_lattice(N,α,β,FLAG_weak_impurity)

    # Get the optimum number of samples
    M_opt[i] = get_M_opt_glauber_dynamics_AL(adj_matrix, τ, L, M_g, M_factor)

    # Update the guess of M_opt
    @printf("beta=%f, M_opt=%d \n", β, M_opt[i])

    M_g = Int(floor(1.8*M_opt[i]))
    println(M_g)
end

writedlm(file_M_opt_gm,M_opt)


text:
"""
    module Atomdata

Holds the data for isolated atoms.
"""
module Atomdata
using ..AtomicMod:makeatom
#using ..Atomic:AtomicMod

"""
    atoms::Dict()

Periodic table information.
"""
atoms = Dict()

#      name                     name  Z  row col mass nval nsemi  orbitals total_energy    e_orb1, e_orb2
#                                                (amu)                      (Ryd)             (eV)
#
atoms["T"] = makeatom("T",    0, 0, 0, 0.0,  1.0,   0,   [:s], 0.0, [0.0], 0.0);
atoms["X"] = makeatom("X",    0, 0, 0, 0.0,  1.0,   0,   [:s], 0.0, [0.0], 0.0, 0.5);
atoms["Xa"] = makeatom("Xa",    0, 0, 0, 0.0,  1.0,   0,   [:s], 3.0, [3.0* 13.605693122], 0.0, 0.5);

atoms["Is"] = makeatom("Si", 14, 3, 4, 28.09,  4.0,   0,   [:s, :p], -9.15243175, [-10.6594, -3.9477],0.007075865418, 0.5300262399999837);

atoms["H" ] = makeatom("H" , 1, 1, 1, 1.0079,   1.0,   0,   [:s ],     -0.90531367, [-6.3233], 0.0, 0.9194820210526328);
#atoms["H" ] = makeatom("H" , 1, 1, 1, 1.0079,   1.0,   0,   [:s, :p ],     -0.90531367, [-6.3233, 0.3], 0.0, 0.9194820210526328);

atoms["Li"] = makeatom("Li", 3, 2, 1, 6.941,    1.0,   2,   [:s, :p], -14.23107803, [-2.7786, -1.0063], 0.0, 0.3666106694736778);
#atoms["Li"] = makeatom("Li", 3, 2, 1, 6.941,    1.0,   2,   [:s], -14.23107803, [-2.7786], 0.0, 0.3666106694736778);

atoms["Be"] = makeatom("Be", 4, 2, 2, 9.012,    2.0,   2,   [:s, :p], -27.80694399, [-5.5055, -1.9313], 0.0, 0.5967315368420951);
atoms["B" ] = makeatom("B" , 5, 2, 3, 10.81,    3.0,   0,   [:s, :p], -5.84632720, [-9.3213, -3.5217], 0.0, 0.6907072084210543);

atoms["C" ] = makeatom("C" , 6, 2, 4, 12.011,   4.0,   0,   [:s, :p], -10.78129974, [-13.5882, -5.1806],0.0, 0.8299248757894666 );
atoms["N" ] = makeatom("N" , 7, 2, 5, 14.007,   5.0,   0,   [:s, :p], -19.47573456, [-18.3876, -6.9917],0.0, 0.9670884968421092 );
atoms["O" ] = makeatom("O" , 8, 2, 6, 15.999,   6.0,   0,   [:s, :p], -31.71324769, [-23.7250, -8.9381],0.0, 1.100101094736835 );
atoms["F" ] = makeatom("F" , 9, 2, 7, 18.998,   7.0,   0,   [:s, :p], -48.22248233, [-29.6811, -11.0125],0.0, 1.2312884800000015 );

#atoms["C" ] = makeatom("C" , 6, 2, 4, 12.011,   4.0,   0,   [:s, :p], -10.77360203, [-13.6099, -5.2024],0.003505246254, 0.8299248757894666 );
#atoms["N" ] = makeatom("N" , 7, 2, 5, 14.007,   5.0,   0,   [:s, :p], -19.47575537, [-18.4071, -7.0112],0.0, 0.9670884968421092 );
#atoms["O" ] = makeatom("O" , 8, 2, 6, 15.999,   6.0,   0,   [:s, :p], -31.71327683, [-23.7413, -8.9545],0.0, 1.100101094736835 );
#atoms["F" ] = makeatom("F" , 9, 2, 7, 18.998,   7.0,   0,   [:s, :p], -48.22196210, [-29.6972, -11.0286],0.0, 1.2312884800000015 );

#atoms["F" ] = makeatom("F" , 9, 2, 7, 18.998,   5.0,   0,   [:p], -48.22196210, [-11.0286],0.0, 1.2312884800000015 );

atoms["Na"] = makeatom("Na", 11, 3, 1, 22.98,  1.0,   8,   [:s, :p], -95.06749946, [-2.6455, -0.6556], 0.0, 0.35322667789488926 );
atoms["Mg"] = makeatom("Mg", 12, 3, 2, 24.31,  2.0,   8,   [:s, :p], -124.99127898, [-4.6187, -1.2416], 0.0, 0.49684460631582567);
atoms["Al"] = makeatom("Al", 13, 3, 3, 26.98,  3.0,   0,   [:s, :p],  -6.37406439, [-7.6515, -2.6256], 0.0, 0.45958674526316);

#atoms["Si"] = makeatom("Si", 14, 3, 4, 28.09,  4.0,   0,   [:s, :p, :d], -9.15243181, [-10.7049, -3.9932, -1.0 ],0.0, 0.5459773389473677);
#atoms["Si"] = makeatom("Si", 14, 3, 4, 28.09,  4.0,   0,   [:s, :d, :p], -9.15243181, [-10.7049, -1.0, -3.9932],0.0, 0.5459773389473677);

atoms["Si"] = makeatom("Si", 14, 3, 4, 28.09,  4.0,   0,   [:s, :p], -9.15243181, [-10.7049, -3.9932],0.0, 0.5459773389473677);

#atoms["Is"] = makeatom("Si", 14, 3, 4, 28.09,  4.0,   0,   [:s, :p], -9.15243175, [-10.6594, -3.9477],0.007075865418, 0.5300262399999837);
atoms["P" ] = makeatom("P" , 15, 3, 5, 30.97,  5.0,   0,   [:s, :p], -15.07327848, [-13.8557, -5.4181], 0.0, 0.6280724631578809);
atoms["S" ] = makeatom("S" , 16, 3, 6, 32.06,  6.0,   0,   [:s, :p], -23.70562681, [-17.1360, -6.9201], 0.0, 0.7069649852631505);
atoms["Cl"] = makeatom("Cl", 17, 3, 7, 35.45,  7.0,   0,   [:s, :p], -33.06905525, [-20.5627, -8.5062], 0.0, 0.7840431242104959);

atoms["K" ] = makeatom("K" , 19, 4, 1, 39.098,  1.0,   8,   [:s, :d, :p], -56.84482202 , [-2.1989, -0.1312, -0.6747], 0.0, 0.29296732631577904);
#atoms["K" ] = makeatom("K" , 19, 4, 1, 39.098,  1.0,   8,   [:s, :d], -56.84482202 , [-2.1989, -0.1312], 0.0, 0.29296732631577904);
#atoms["K" ] = makeatom("K" , 19, 4, 1, 39.098,  1.0,   8,   [:s, :d, :p], -56.84482202 , [-2.1989, -0.1312, -0.6747], 0.0, 0.29296732631577904);
#atoms["K" ] = makeatom("K" , 19, 4, 1, 39.098,  1.0,   8,   [:s, :p], -56.84482202 , [-2.1989, -0.6747], 0.0, 0.29296732631577904);
atoms["Ca"] = makeatom("Ca", 20, 4, 2, 40.078,  2.0,   8,   [:s, :d, :p], -74.68919704, [-3.6547, -1.8454, -1.2718], 0.0, 0.39270917052637117 * 1.1);

atoms["Sc"] = makeatom("Sc", 21, 4, 2.05, 44.956,  3.0,   8,   [:s, :d, :p], -93.89389249,  [-4.0850, -3.1428, -1.3456],  0.0, 0.494323979080634  );
atoms["Ti"] = makeatom("Ti", 22, 4, 2.15, 47.867,  4.0,   8,   [:s, :d, :p], -118.63461912, [-4.3470, -4.0656, -1.3310],  0.0, 0.450666224264095   *1.2 );
atoms["V" ] = makeatom("V" , 23, 4, 2.25, 50.941,  5.0,    8,   [:s, :d, :p], -144.34826807, [-4.3727, -4.2223, -1.2153], 0.0, 0.44622367939588864 *1.2);
atoms["Cr"] = makeatom("Cr", 24, 4, 2.35, 51.996,  6.0,    8,   [:s, :d, :p], -174.69181247, [-4.3746, -4.2928, -1.0764], 0.0, 0.4399949632549368  *1.2);
atoms["Mn"] = makeatom("Mn", 25, 4, 2.45, 54.938,  7.0,    8,   [:s, :d, :p], -210.40629067, [-4.3786, -4.3505, -0.9746], 0.0, 0.44538861740151653 *1.2);
atoms["Fe"] = makeatom("Fe", 26, 4, 2.55, 55.845,  8.0,    8,   [:s, :d, :p], -249.56808439, [-4.3911, -4.4108, -0.8715], 0.0, 0.4851813000402511  *1.2);
atoms["Co"] = makeatom("Co", 27, 4, 2.65, 58.993,  9.0,    8,   [:s, :d, :p], -297.72249071, [-4.2335, -4.2995, -0.5167], 0.0, 0.4955992980049762  *1.2);
atoms["Ni"] = makeatom("Ni", 28, 4, 2.75, 58.963,  10.0,   8,   [:s, :d, :p], -342.65159683, [-4.2480, -4.3731, -0.3946], 0.0, 0.6142338677950465  *1.2);
atoms["Cu"] = makeatom("Cu", 29, 4, 2.85, 63.546,  11.0,   8,   [:s, :d, :p], -403.19763559, [-4.5506, -4.9911, -0.5855], 0.0, 0.5142569938899915  *1.2);
atoms["Zn"] = makeatom("Zn", 30, 4, 2.95, 65.38,   12.0,   8,   [:s, :d, :p], -461.12493484, [-5.9305, -10.1008, -1.0049], 0.0, 0.7190857528540386);

#atoms["Ga"] = makeatom("Ga", 30, 4, 3.0,  69.72,   3.0,   16,   [:s, :p], -414.75774465, [-8.8812, -2.4866], 0.0, 0.4697093810);

atoms["Ga"] = makeatom("Ga", 31, 4, 3.0,  69.72,   13.0,   6,   [:s,:d, :p], -414.75774465, [-8.8812,-19.0647, -2.4866], 0.0, 0.4697093810);

#atoms["Ge"] = makeatom("Ge", 31, 4, 4.0,  72.63,   4.0,    10,   [:s,:p, :d], -212.84384318, [-11.6643,  -3.8180, -1.0], 0.0, 0.5314231326);
#atoms["Ge"] = makeatom("Ge", 31, 4, 4.0,  72.63,   4.0,    10,   [:s,:d, :p], -212.84384318, [-11.6643, -1.0, -3.8180], 0.0, 0.5314231326);

atoms["Ge"] = makeatom("Ge", 32, 4, 4.0,  72.63,   4.0,    10,   [:s, :p], -212.84384318, [-11.6643, -3.8180], 0.0, 0.5314231326);

atoms["As"] = makeatom("As", 33, 4, 5.0,  74.92,   5.0,    0,   [:s, :p], -39.62501000, [-14.4349, -5.1279], 0.0,  0.5906811873);
atoms["Se"] = makeatom("Se", 34, 4, 6.0,  78.97,   6.0,    0,   [:s, :p],  -43.03651432, [-17.2230, -6.4409], 0.0, 0.6459585094);
atoms["Br"] = makeatom("Br", 35, 4, 7.0,  79.90,   7.0,    0,   [:s, :p], -40.38769338, [-20.0640, -7.7862], 0.0,  0.6992542568);

###

atoms["Rb"] = makeatom("Rb", 37, 5, 1, 85.468,  1.0,   8,   [:s, :d, :p], -53.05281238, [-2.0994, -0.3, -0.5921], 0.0, 0.282637);
#atoms["Rb"] = makeatom("Rb", 37, 5, 1, 85.468,  1.0,   8,   [:s, :p], -53.05281238, [-2.0994, -0.5921], 0.0, 0.282637);
atoms["Sr"] = makeatom("Sr", 38, 5, 2, 87.62,  2.0,   8,   [:s, :d, :p], -69.92604259, [-3.4017, -1.1543, -1.1318], 0.0,  0.365728 * 1.1);

atoms["Y" ] = makeatom("Y" , 39, 5, 2.05, 88.906,  3.0,   8,    [:s, :d, :p], -92.30423365,  [-3.9864, -2.3907, -1.2867], 0.0,  0.5538339);
atoms["Zr"] = makeatom("Zr", 40, 5, 2.15, 91.224,  4.0,   8,    [:s, :d, :p], -98.57421435,  [-4.3347, -3.4817, -1.3201], 0.0,  0.4710552712651725 *1.2);
atoms["Nb"] = makeatom("Nb", 41, 5, 2.25, 92.906,  5.0,    8,   [:s, :d, :p], -117.19066274 ,  [-4.5053, -4.3124, -1.3159],0.0, 0.43108825397374606*1.2);
atoms["Mo"] = makeatom("Mo", 42, 5, 2.35, 95.95,   6.0,    8,   [:s, :d, :p], -138.51328928,  [-4.2469, -4.2022, -1.0704], 0.0, 0.41259037004299504*1.2);
atoms["Tc"] = makeatom("Tc", 43, 5, 2.45, 98.0,    7.0,    8,   [:s, :d, :p], -173.39188188,  [-3.9325, -3.9874, -0.7985], 0.0, 0.4367177272916297*1.2);
atoms["Ru"] = makeatom("Ru", 44, 5, 2.55, 101.07,  8.0,    8,   [:s, :d, :p], -193.54935366,  [-3.5929, -3.7493, -0.4736], 0.0, 0.5163241436955993*1.2);
atoms["Rh"] = makeatom("Rh", 45, 5, 2.65, 102.91,  9.0,    6,   [:s, :d, :p], -170.64161619,  [-3.3079, -3.6154, -0.2568], 0.0, 0.5892015173558101*1.2);
atoms["Pd"] = makeatom("Pd", 46, 5, 2.75, 106.42,  10.0,   6,   [:s, :d, :p], -198.31288985,  [-3.2109, -3.9641, -0.1564], 0.0, 0.6171662403405886*1.2);
atoms["Ag"] = makeatom("Ag", 47, 5, 2.85, 107.87,  11.0,   8,   [:s, :d, :p], -295.31434496,  [-4.3268, -7.3569, -0.5745], 0.0, 0.5141000517496574*1.2);
atoms["Cd"] = makeatom("Cd", 48, 5, 2.95, 112.41,  12.0,   0,   [:s, :d, :p], -114.39800378,  [-5.6202, -11.5954, -1.0602], 0.0, 0.6749099495485315);

#atoms["In"] = makeatom("In", 49, 5, 3.0,  114.82,   3.0,   10,   [:s, :p], -132.97861648 , [-8.1551, -2.4193], 0.0, 0.4292669);
atoms["In"] = makeatom("In", 49, 5, 3.0,  114.82,   13.0,   0,   [:s,:d, :p], -132.97861648 , [-8.1551, -18.4828, -2.4193], 0.0, 0.4292669);
atoms["Sn"] = makeatom("Sn", 50, 5, 4.0,  118.71,   4.0,   10,   [:s, :p], -158.61463088, [-10.4877, -3.5995], 0.0,  0.4766457);
atoms["Sb"] = makeatom("Sb", 51, 5, 5.0,  121.76,   5.0,   10,   [:s, :p], -184.28336729, [-12.7656, -4.7332], 0.0,  0.5212358);
atoms["Te"] = makeatom("Te", 52, 5, 6.0,  127.60,   6.0,    0,   [:s, :p], -26.19141773, [-15.0608, -5.8767], 0.0,   0.5638125);
atoms["I" ] = makeatom("I" , 53, 5, 7.0,  126.90,   7.0,    0,   [:s, :p], -65.01087258 , [-17.3546, -7.0038], 0.0,  0.6039061);

###

atoms["Cs"] = makeatom("Cs", 55, 6, 1, 132.9,  1.0,   8,   [:s, :d, :p ], -62.83174426, [-1.9276, -0.3, -0.5431 ], 0.0, 0.262230);
#atoms["Cs"] = makeatom("Cs", 55, 6, 1, 132.9,  1.0,   8,   [:s, :p ], -62.83174426, [-1.9276, -0.5431 ], 0.0, 0.262230);
atoms["Ba"] = makeatom("Ba", 56, 6, 2, 139.3,  2.0,   8,   [:s, :d, :p], -70.01638694, [-3.0703, -1.8523, -1.0457], 0.0, 0.332378 * 1.1);

#atoms["La"] = makeatom("La", 57, 6, 2.05, 138.9,  3.0,   8,    [:s, :d, :p], -101.82373826,  [-3.0750,-2.2784,-0.6575], 0.0, 0.5); #f orbitals?
atoms["La"] = makeatom("La", 57, 6, 2.05, 138.9,  3.0,   8,    [:s, :d], -101.82373826,  [-3.0750,-2.2784], 0.0, 0.5); #f orbitals?
atoms["Hf"] = makeatom("Hf", 72, 6, 2.15, 178.5,  4.0,   8,    [:s, :d, :p], -157.97552412,  [-5.0187, -2.6052, -1.2994], 0.0, 0.5723412); #lan
atoms["Ta"] = makeatom("Ta", 73, 6, 2.25, 181.0,  5.0,    8,   [:s, :d, :p], -141.00656328,  [-5.3481, -3.4993, -1.2418], 0.0, 0.6127437);
atoms["W" ] = makeatom("W" , 74, 6, 2.35, 183.8,  6.0,    8,   [:s, :d, :p], -158.16868802,  [-5.6245, -4.3914, -1.2092], 0.0, 0.6426146);
atoms["Re"] = makeatom("Re", 75, 6, 2.45, 186.2,  7.0,    8,   [:s, :d, :p], -189.48871022,  [-5.8618, -5.2763, -1.1642], 0.0, 0.5348510*1.2);
atoms["Os"] = makeatom("Os", 76, 6, 2.55, 190.2,  8.0,    8,   [:s, :d, :p], -197.71975092,  [-5.9152, -5.8238, -1.0547], 0.0, 0.5172697301479442*1.2);
atoms["Ir"] = makeatom("Ir", 77, 6, 2.65, 192.2,  9.0,    6,   [:s, :d, :p], -180.96527345,  [-5.5876, -5.6416, -0.8041], 0.0, 0.5004484210685588*1.2);
atoms["Pt"] = makeatom("Pt", 78, 6, 2.75, 195.1,  10.0,   6,   [:s, :d, :p], -210.06169163,  [-5.3377, -5.5187, -0.5975], 0.0, 0.5793362);
atoms["Au"] = makeatom("Au", 79, 6, 2.85, 197.0,  11.0,   0,   [:s, :d, :p], -114.46631173,  [-5.7155, -6.8176, -0.6483], 0.0, 0.5508305);
atoms["Hg"] = makeatom("Hg", 80, 6, 2.95, 200.6,  12.0,   0,   [:s, :d, :p], -107.48183748,  [-9.7901, -6.7757, -0.8974], 0.0, 0.5915364);

atoms["Tl"] = makeatom("Tl", 81, 6, 3.0,  204.4,   13.0,   0,   [:s, :d, :p], -144.53037033, [ -9.4990,-15.3839,  -2.2651], 0.0, 0.4202546); #include d
#atoms["Tl"] = makeatom("Tl", 81, 6, 3.0,  204.4,   13.0,   0,   [:d, :s, :p], -144.53037033, [-15.3839, -9.4990,  -2.2651], 0.0, 0.4202546); #include d
#atoms["Tl"] = makeatom("Tl", 81, 6, 3.0,  204.4,   3.0,   10,   [:s, :p], -144.53037033, [-9.4990, -2.2651], 0.0, 0.4202546);

atoms["Pb"] = makeatom("Pb", 82, 6, 4.0,  207.2,   4.0,   10,   [:s, :p], -163.61141039, [-11.9754, -3.3964], 0.0, 0.4606776);
atoms["Bi"] = makeatom("Bi", 83, 6, 5.0,  209.0,   5.0,   10,   [:s, :p], -184.68674216, [-14.4080, -4.4700], 0.0, 0.4996092);

for key in keys(atoms)
    atoms[Symbol(key)] = atoms[key]
end

#U = IE - EA  (ionization energy - electron affinity

#stop printing
nothing

#approximate metallic radius / single covalent radius
#for use in prerelaxing structures in pm
"""
    atom_radius::Dict()

Atomic radius (metallic), in pm
"""
atom_radius = Dict()

atom_radius[ "X" ] =  200
atom_radius[ "Xa" ] =  200
atom_radius[ "T" ] =  200

atom_radius[ "H" ] =  60
atom_radius[ "He" ] =  60
atom_radius[ "Li" ] =  152
atom_radius[ "Be" ] =  112
atom_radius[ "B" ] =   82
atom_radius[ "C" ] =   77
atom_radius[ "N" ] =   75
atom_radius[ "O" ] =   73
atom_radius[ "F" ] =   71
atom_radius[ "Ne" ] =  80
atom_radius[ "Na" ] =  186
atom_radius[ "Mg" ] =  160
atom_radius[ "Al" ] =  143
atom_radius[ "Si" ] =  111
atom_radius[ "P" ] =   106
atom_radius[ "S" ] =   102
atom_radius[ "Cl" ] =  99
atom_radius[ "Ar" ] =  100
atom_radius[ "K" ] =   227
atom_radius[ "Ca" ] =  197
atom_radius[ "Sc" ] =  162
atom_radius[ "Ti" ] =  147
atom_radius[ "V" ] =   134
atom_radius[ "Cr" ] =  128
atom_radius[ "Mn" ] =  127
atom_radius[ "Fe" ] =  126
atom_radius[ "Co" ] =  125
atom_radius[ "Ni" ] =  124
atom_radius[ "Cu" ] =  128
atom_radius[ "Zn" ] =  134
atom_radius[ "Ga" ] =  135
atom_radius[ "Ge" ] =  122
atom_radius[ "As" ] =  119
atom_radius[ "Se" ] =  116
atom_radius[ "Br" ] =  114
atom_radius[ "Kr" ] =  100
atom_radius[ "Rb" ] =  248
atom_radius[ "Sr" ] =  215
atom_radius[ "Y" ] =  180
atom_radius[ "Zr" ] =  160
atom_radius[ "Nb" ] =  146
atom_radius[ "Mo" ] =  139
atom_radius[ "Tc" ] =  136
atom_radius[ "Ru" ] =  134
atom_radius[ "Rh" ] =  134
atom_radius[ "Pd" ] =  137
atom_radius[ "Ag" ] =  144
atom_radius[ "Cd" ] =  151
atom_radius[ "In" ] =  167
atom_radius[ "Sn" ] =  141
atom_radius[ "Sb" ] =  138
atom_radius[ "Te" ] =  135
atom_radius[ "I" ] =   133
atom_radius[ "Xe" ] =  130
atom_radius[ "Cs" ] =  265
atom_radius[ "Ba" ] =  222
atom_radius[ "La" ] =  187
atom_radius[ "Ce" ] =  181
atom_radius[ "Pr" ] =  181
atom_radius[ "Nd" ] =  181
atom_radius[ "Pm" ] =  181
atom_radius[ "Sm" ] =  181
atom_radius[ "Eu" ] =  181
atom_radius[ "Gd" ] =  181
atom_radius[ "Tb" ] =  181
atom_radius[ "Dy" ] =  181
atom_radius[ "Ho" ] =  181
atom_radius[ "Er" ] =  181
atom_radius[ "Tm" ] =  181
atom_radius[ "Yb" ] =  181
atom_radius[ "Lu" ] =  173
atom_radius[ "Hf" ] =  159
atom_radius[ "Ta" ] =  146
atom_radius[ "W" ] =  139
atom_radius[ "Re" ] =  137
atom_radius[ "Os" ] =  135
atom_radius[ "Ir" ] =  135
atom_radius[ "Pt" ] =  135
atom_radius[ "Au" ] =  144
atom_radius[ "Hg" ] =  151
atom_radius[ "Tl" ] =  170
atom_radius[ "Pb" ] =  147
atom_radius[ "Bi" ] =  146
atom_radius[ "Po" ] =  140
atom_radius[ "At" ] =  140

for key in keys(atom_radius)
    atom_radius[Symbol(key)] = atom_radius[key]
end

"""
    atom_prefered_oxidation::Dict()

Prefered oxidation states of atoms, in descending order of preference (approximate).
"""
atom_prefered_oxidation = Dict()

atom_prefered_oxidation[ "H" ] =  [1, -1, 0]
atom_prefered_oxidation[ "He" ] = [0]
atom_prefered_oxidation[ "Li" ] = [1, 0]
atom_prefered_oxidation[ "Be" ] = [2, 0]
atom_prefered_oxidation[ "B" ] =  [3, 0]
atom_prefered_oxidation[ "C" ] =  [4, -4, 0]
atom_prefered_oxidation[ "N" ] =  [-3, 0]
atom_prefered_oxidation[ "O" ] =  [-2, 0]
atom_prefered_oxidation[ "F" ] =  [-1, 0]
atom_prefered_oxidation[ "Ne" ] = [0]
atom_prefered_oxidation[ "Na" ] = [1, 0]
atom_prefered_oxidation[ "Mg" ] = [2, 0]
atom_prefered_oxidation[ "Al" ] = [3, 0]
atom_prefered_oxidation[ "Si" ] = [4, -4, 2, 0]
atom_prefered_oxidation[ "P" ] =  [-3, 3, 5, 0]
atom_prefered_oxidation[ "S" ] =  [-2, 6, 0]
atom_prefered_oxidation[ "Cl" ] = [-1, 0]
atom_prefered_oxidation[ "Ar" ] = [0]
atom_prefered_oxidation[ "K" ] =  [1, 0]
atom_prefered_oxidation[ "Ca" ] = [2, 0]
atom_prefered_oxidation[ "Sc" ] = [3, 2, 0]
atom_prefered_oxidation[ "Ti" ] = [4, 3, 2, 0]
atom_prefered_oxidation[ "V" ] =  [5, 3, 0]
atom_prefered_oxidation[ "Cr" ] = [4, 3, 6, 2, 0]
atom_prefered_oxidation[ "Mn" ] = [4, 3, 2, 0]
atom_prefered_oxidation[ "Fe" ] = [4, 3, 2, 0]
atom_prefered_oxidation[ "Co" ] = [3, 4, 2, 0]
atom_prefered_oxidation[ "Ni" ] = [2, 4, 1, 0]
atom_prefered_oxidation[ "Cu" ] = [2, 4, 1, 0]
atom_prefered_oxidation[ "Zn" ] = [2, 0]
atom_prefered_oxidation[ "Ga" ] = [3, 0]
atom_prefered_oxidation[ "Ge" ] = [4, -4, 2, 0]
atom_prefered_oxidation[ "As" ] = [-3, 3, 5, 0]
atom_prefered_oxidation[ "Se" ] = [-2, 2, 4, 6, 0]
atom_prefered_oxidation[ "Br" ] = [-1, 0]
atom_prefered_oxidation[ "Kr" ] = [0]
atom_prefered_oxidation[ "Rb" ] = [1, 0]
atom_prefered_oxidation[ "Sr" ] = [2, 0]
atom_prefered_oxidation[ "Y" ] =  [3, 0]
atom_prefered_oxidation[ "Zr" ] = [4, 2, 0]
atom_prefered_oxidation[ "Nb" ] = [5, 0]
atom_prefered_oxidation[ "Mo" ] = [6, 4, 0]
atom_prefered_oxidation[ "Tc" ] = [4, 3, 2, 0]
atom_prefered_oxidation[ "Ru" ] = [4, 3, 2, 0]
atom_prefered_oxidation[ "Rh" ] = [4, 3, 2, 0]
atom_prefered_oxidation[ "Pd" ] = [4, 2, 0]
atom_prefered_oxidation[ "Ag" ] = [1, 2, 0]
atom_prefered_oxidation[ "Cd" ] = [2,0]
atom_prefered_oxidation[ "In" ] = [3,0]
atom_prefered_oxidation[ "Sn" ] = [4, 2, -4, 0]
atom_prefered_oxidation[ "Sb" ] = [-3, 3, 5, 0]
atom_prefered_oxidation[ "Te" ] = [-2, 6, 0]
atom_prefered_oxidation[ "I" ] =  [-1, 0]
atom_prefered_oxidation[ "Xe" ] = [0]
atom_prefered_oxidation[ "Cs" ] = [1, 0]
atom_prefered_oxidation[ "Ba" ] = [2, 0]
atom_prefered_oxidation[ "La" ] = [3, 0]
atom_prefered_oxidation[ "Ce" ] = [3, 0]
atom_prefered_oxidation[ "Pr" ] = [3, 0]
atom_prefered_oxidation[ "Nd" ] = [3, 0]
atom_prefered_oxidation[ "Pm" ] = [3, 0]
atom_prefered_oxidation[ "Sm" ] = [3, 0]
atom_prefered_oxidation[ "Eu" ] = [2, 3, 0]
atom_prefered_oxidation[ "Gd" ] = [3, 0]
atom_prefered_oxidation[ "Tb" ] = [3, 0]
atom_prefered_oxidation[ "Dy" ] = [3, 0]
atom_prefered_oxidation[ "Ho" ] = [3, 0]
atom_prefered_oxidation[ "Er" ] = [3, 0]
atom_prefered_oxidation[ "Tm" ] = [3, 0]
atom_prefered_oxidation[ "Yb" ] = [3, 0]
atom_prefered_oxidation[ "Lu" ] = [3, 0]
atom_prefered_oxidation[ "Hf" ] = [4, 2, 0]
atom_prefered_oxidation[ "Ta" ] = [5, 3, 0]
atom_prefered_oxidation[ "W" ] =  [6, 4, 0]
atom_prefered_oxidation[ "Re" ] = [4, 3, 6, 0]
atom_prefered_oxidation[ "Os" ] = [4, 3, 0]
atom_prefered_oxidation[ "Ir" ] = [4, 3, 0]
atom_prefered_oxidation[ "Pt" ] = [2, 4, 0]
atom_prefered_oxidation[ "Au" ] = [3, 1, 0]
atom_prefered_oxidation[ "Hg" ] = [2, 1, 0]
atom_prefered_oxidation[ "Tl" ] = [3, 1, 0]
atom_prefered_oxidation[ "Pb" ] = [4, 2, -4, 0]
atom_prefered_oxidation[ "Bi" ] = [3, -3, 0]
atom_prefered_oxidation[ "Po" ] = [-2, 0]
atom_prefered_oxidation[ "At" ] = [-1, 0]

for key in keys(atom_prefered_oxidation)
    atom_prefered_oxidation[Symbol(key)] = atom_prefered_oxidation[key]
end

"""
    cutoff_dist::Dict()

Cutoff distance for Hamiltonian calculation
"""
cutoff_dist = Dict()

function get_cutoff(at1, at2)

    if (at1,at2) in keys(cutoff_dist)
        return cutoff_dist[(at1,at2)]
    else
        rad1 = atom_radius[at1] / 0.529 / 100.0
        rad2 = atom_radius[at2] / 0.529 / 100.0

        cutoff2X    = (rad1 + rad2) / 2.0 * 7.0
        cutoff_onX  = (rad1 + rad2) / 2.0 * 6.0

        cutoff2X   = max(min(cutoff2X,   19.01), 15.01) #2body
        cutoff_onX = max(min(cutoff_onX, 18.01), 14.51) #onsite

        cutoff_dist[(at1,at2)] = [cutoff2X, cutoff_onX]
        cutoff_dist[(at2,at1)] = [cutoff2X, cutoff_onX]

        return [cutoff2X, cutoff_onX]
    end
end

function get_cutoff(at1, at2, at3)
    if (at1,at2,at3) in keys(cutoff_dist)
        return cutoff_dist[(at1,at2,at3)]
    else
        rad1 = atom_radius[at1] / 0.529 / 100.0
        rad2 = atom_radius[at2] / 0.529 / 100.0
        rad3 = atom_radius[at3] / 0.529 / 100.0

        cutoff3bX = minimum([rad1, rad2, rad3])*5.0
#        cutoff3bX    = (rad1 + rad2 + rad3) / 3.0 * 5.0
        cutoff3bX  = max(min(cutoff3bX,  13.51),  10.01) #3body

        cutoff_dist[(at1,at2,at3)] = cutoff3bX
        cutoff_dist[(at1,at3,at2)] = cutoff3bX
        cutoff_dist[(at2,at1,at3)] = cutoff3bX
        cutoff_dist[(at2,at3,at1)] = cutoff3bX
        cutoff_dist[(at3,at1,at2)] = cutoff3bX
        cutoff_dist[(at3,at2,at1)] = cutoff3bX

        return cutoff3bX
    end
end

"""
    min_dimer_dist_dict::Dict()

Approximate minimum distance in current fitting data.
"""
min_dimer_dist_dict = Dict()

min_dimer_dist_dict["Pd"] =  3.73246
min_dimer_dist_dict["Si"] =  3.33537
min_dimer_dist_dict["C"] =  1.92837
min_dimer_dist_dict["P"] =  2.86897
min_dimer_dist_dict["Nb"] =  3.10416
min_dimer_dist_dict["Ag"] =  3.8093
min_dimer_dist_dict["Ru"] =  2.95446
min_dimer_dist_dict["Sb"] =  3.75256
min_dimer_dist_dict["Cs"] =  4.49968
min_dimer_dist_dict["Be"] =  3.2197
min_dimer_dist_dict["Sr"] =  3.95284
min_dimer_dist_dict["Ta"] =  3.2543
min_dimer_dist_dict["Ga"] =  3.94627
min_dimer_dist_dict["Te"] =  3.89208
min_dimer_dist_dict["Y"] =  4.13392
min_dimer_dist_dict["Ca"] =  5.03584
min_dimer_dist_dict["K"] =  4.29261
min_dimer_dist_dict["Au"] =  3.74318
min_dimer_dist_dict["V"] =  2.59038
min_dimer_dist_dict["Mo"] =  2.88349
min_dimer_dist_dict["Br"] =  3.46016
min_dimer_dist_dict["Zn"] =  3.88444
min_dimer_dist_dict["Sc"] =  3.49735
min_dimer_dist_dict["H"] =  1.16274
min_dimer_dist_dict["Al"] =  3.9947
min_dimer_dist_dict["S"] =  2.87718
min_dimer_dist_dict["Ge"] =  3.53994
min_dimer_dist_dict["Fe"] =  2.68548
min_dimer_dist_dict["Mg"] =  4.609
min_dimer_dist_dict["F"] =  2.12195
min_dimer_dist_dict["La"] =  3.92813
min_dimer_dist_dict["W"] =  3.03842
min_dimer_dist_dict["Li"] =  2.88861
min_dimer_dist_dict["O"] =  1.84197
min_dimer_dist_dict["B"] =  2.47166
min_dimer_dist_dict["Mn"] =  2.47598
min_dimer_dist_dict["Bi"] =  3.96946
min_dimer_dist_dict["Re"] =  3.06043
min_dimer_dist_dict["Hf"] =  3.71755
min_dimer_dist_dict["Cl"] =  3.00634
min_dimer_dist_dict["In"] =  4.58286
min_dimer_dist_dict["Rb"] =  4.50733
min_dimer_dist_dict["Se"] =  3.29481
min_dimer_dist_dict["Ba"] =  4.10738
min_dimer_dist_dict["Rh"] =  3.20676
min_dimer_dist_dict["Ti"] =  2.86562
min_dimer_dist_dict["Zr"] =  3.39044
min_dimer_dist_dict["Cr"] =  2.39985
min_dimer_dist_dict["Cu"] =  3.29852
min_dimer_dist_dict["Tc"] =  2.90788
min_dimer_dist_dict["Tl"] =  4.8252
min_dimer_dist_dict["Hg"] =  4.27212
min_dimer_dist_dict["N"] =  1.67319
min_dimer_dist_dict["Na"] =  3.50215
min_dimer_dist_dict["Os"] =  3.09582
min_dimer_dist_dict["Ni"] =  3.11813
min_dimer_dist_dict["Co"] =  2.91763
min_dimer_dist_dict["Pt"] =  3.47029
min_dimer_dist_dict["Ir"] =  3.25225
min_dimer_dist_dict["As"] =  3.18345
min_dimer_dist_dict["Sn"] =  4.14125
min_dimer_dist_dict["I"] =  4.02336
min_dimer_dist_dict["Cd"] =  4.29162
min_dimer_dist_dict["Pb"] =  4.35765

for key in keys(min_dimer_dist_dict)
    min_dimer_dist_dict[Symbol(key)] = min_dimer_dist_dict[key]
end

########
"""
    sub_list::Dict()

Most likely substitution list, approximate.
"""
sub_list = Dict()

sub_list[ "H" ] =  ["F","Li"]
sub_list[ "He" ] = []
sub_list[ "Li" ] = ["Na"]
sub_list[ "Be" ] = ["Mg"]
sub_list[ "B" ] =  ["Al"]
sub_list[ "C" ] =  ["Si"]
sub_list[ "N" ] =  ["P"]
sub_list[ "O" ] =  ["S"]
sub_list[ "F" ] =  ["Cl"]
sub_list[ "Ne" ] = []
sub_list[ "Na" ] = ["K", "Li"]
sub_list[ "Mg" ] = ["Be", "Ca"]
sub_list[ "Al" ] = ["Ga", "Sc"]
sub_list[ "Si" ] = ["Ge", "Ti"]
sub_list[ "P" ] = ["As", "N"]
sub_list[ "S" ] = ["Se", "O"]
sub_list[ "Cl" ] = ["F", "Br"]
sub_list[ "Ar" ] = []
sub_list[ "K" ] = ["Rb", "Na"]
sub_list[ "Ca" ] = ["Mg", "Sr"]
sub_list[ "Sc" ] = ["Al", "Ga", "Y"]
sub_list[ "Ti" ] = ["Zr", "Hf"]
sub_list[ "V" ] = ["Nb", "Ti"]
sub_list[ "Cr" ] = ["V", "Mn", "Mo"]
sub_list[ "Mn" ] = ["Cr", "Fe", "Ru"]
sub_list[ "Fe" ] = ["Mn", "Co", "Ru"]
sub_list[ "Co" ] = ["Fe", "Ni", "Rh"]
sub_list[ "Ni" ] = ["Co", "Pd", "Cu"]
sub_list[ "Cu" ] = ["Ni", "Zn", "Ag"]
sub_list[ "Zn" ] = ["Mg", "Cd", "Ca"]
sub_list[ "Ga" ] = ["Al", "In", "Sc"]
sub_list[ "Ge" ] = ["Si", "Sn", "Ti"]
sub_list[ "As" ] = ["Sb", "P", "Ga"]
sub_list[ "Se" ] = ["S", "Te"]
sub_list[ "Br" ] = ["Cl", "I"]
sub_list[ "Kr" ] = []
sub_list[ "Rb" ] = ["Cs", "K"]
sub_list[ "Sr" ] = ["Ca", "Ba", "Pb"]
sub_list[ "Y" ] = ["Cs", "Bi"]
sub_list[ "Zr" ] = ["Ti", "Hf", "Ge"]
sub_list[ "Nb" ] = ["Ta", "V", "Mo"]
sub_list[ "Mo" ] = ["Cr", "W", "Nb"]
sub_list[ "Tc" ] = ["Mn", "Re", "Ru"]
sub_list[ "Ru" ] = ["Fe", "Os"]
sub_list[ "Rh" ] = ["Co", "Ir"]
sub_list[ "Pd" ] = ["Ni", "Pt"]
sub_list[ "Ag" ] = ["Cu", "Au"]
sub_list[ "Cd" ] = ["Zn", "Hg"]
sub_list[ "In" ] = ["Ga", "Tl", "Sb"]
sub_list[ "Sn" ] = ["Ge", "Pb"]
sub_list[ "Sb" ] = ["As", "Bi"]
sub_list[ "Te" ] = ["Se"]
sub_list[ "I" ] = ["Br"]
sub_list[ "Xe" ] = []
sub_list[ "Cs" ] = ["Rb", "K"]
sub_list[ "Ba" ] = ["Sr"]
sub_list[ "La" ] = ["Y", "Sc"]
sub_list[ "Ce" ] = ["Ti", "Zr", "Hf"]
sub_list[ "Pr" ] = ["La",  "Y"]
sub_list[ "Nd" ] = ["La",  "Y"]
sub_list[ "Pm" ] = ["La",  "Y"]
sub_list[ "Sm" ] = ["La",  "Y"]
sub_list[ "Eu" ] = ["Sr", "Ba", "Y"]
sub_list[ "Gd" ] = ["La",  "Y", "Sc"]
sub_list[ "Tb" ] = ["La",  "Y", "Sc"]
sub_list[ "Dy" ] = ["La",  "Y", "Sc"]
sub_list[ "Ho" ] = ["La",  "Y", "Sc"]
sub_list[ "Er" ] = ["La",  "Y", "Sc"]
sub_list[ "Tm" ] = ["La",  "Y", "Sc"]
sub_list[ "Yb" ] = ["La",  "Y", "Sc"]
sub_list[ "Lu" ] = ["La", "Y", "Sc"]
sub_list[ "Hf" ] = ["Zr", "Ti"]
sub_list[ "Ta" ] = ["Nb", "V"]
sub_list[ "W" ] = ["Mo", "Re", "Ta", "Cr"]
sub_list[ "Re" ] = ["Tc", "Os", "W", "Mn"]
sub_list[ "Os" ] = ["Ru", "Ir", "Re", "Fe"]
sub_list[ "Ir" ] = ["Rh", "Os", "Pt", "Co"]
sub_list[ "Pt" ] = ["Pd", "Ir", "Au", "Ni"]
sub_list[ "Au" ] = ["Ag", "Pt", "Hg", "Cu"]
sub_list[ "Hg" ] = ["Cd", "Pb", "Zn"]
sub_list[ "Tl" ] = ["In", "Bi"]
sub_list[ "Pb" ] = ["Sn", "Sr"]
sub_list[ "Bi" ] = ["Sb", "In"]
sub_list[ "Po" ] = []
sub_list[ "At" ] = []

"""
    electronegativity::Dict()

Pauling scale.
"""
electronegativity = Dict()

electronegativity["H"] =  2.20
electronegativity["He"] =  0.98
electronegativity["Be"] =  1.57
electronegativity["B"] =  2.04
electronegativity["C"] =  2.55
electronegativity["N"] =  3.04
electronegativity["O"] =  3.44
electronegativity["F"] =  3.98
electronegativity["Ne"] = 0.0
electronegativity["Na"] = 0.93
electronegativity["Mg"] = 1.31
electronegativity["Al"] = 1.61
electronegativity["Si"] = 1.90
electronegativity["P"] = 2.19
electronegativity["S"] = 2.58
electronegativity["Cl"] = 3.16
electronegativity["K"] = 0.82
electronegativity["Ca"] = 1.00
electronegativity["Sc"] = 1.36
electronegativity["Ti"] = 1.54
electronegativity["V"] = 1.63
electronegativity["Cr"] = 1.66
electronegativity["Mn"] = 1.55
electronegativity["Fe"] = 1.83
electronegativity["Co"] = 1.88
electronegativity["Ni"] = 1.91
electronegativity["Cu"] = 1.90
electronegativity["Zn"] = 1.65
electronegativity["Ga"] = 1.81
electronegativity["Ge"] = 2.01
electronegativity["As"] = 2.18
electronegativity["Se"] = 2.55
electronegativity["Br"] = 2.96
electronegativity["Kr"] = 3.00
electronegativity["Rb"] = 0.82
electronegativity["Sr"] = 0.95
electronegativity["Y"] = 1.22
electronegativity["Zr"] = 1.33
electronegativity["Nb"] = 1.6
electronegativity["Mo"] = 2.16
electronegativity["Tc"] = 1.9
electronegativity["Ru"] = 2.2
electronegativity["Rh"] = 2.28
electronegativity["Pd"] = 2.20
electronegativity["Ag"] = 1.93
electronegativity["Cd"] = 1.69
electronegativity["In"] = 1.78
electronegativity["Sn"] = 1.96
electronegativity["Sb"] = 2.05
electronegativity["Te"] = 2.1
electronegativity["I"] = 2.66
electronegativity["Xe"] = 2.6
electronegativity["Cs"] = 0.79
electronegativity["Ba"] = 0.89
electronegativity["La"] = 1.10
electronegativity["Ce"] = 1.12
electronegativity["Pr"] = 1.13
electronegativity["Nd"] = 1.14
electronegativity["Pm"] = 1.13
electronegativity["Sm"] = 1.17
electronegativity["Eu"] = 1.2
electronegativity["Gd"] = 1.2
electronegativity["Tb"] = 1.22
electronegativity["Dy"] = 1.23
electronegativity["Ho"] = 1.24
electronegativity["Er"] = 1.24
electronegativity["Tm"] = 1.25
electronegativity["Yb"] = 1.1
electronegativity["Lu"] = 1.27
electronegativity["Hf"] = 1.3
electronegativity["Ta"] = 1.5
electronegativity["WT"] = 2.36
electronegativity["Re"] = 1.9
electronegativity["Os"] = 2.2
electronegativity["Ir"] = 2.2
electronegativity["Pt"] = 2.28
electronegativity["Au"] = 2.54
electronegativity["Hg"] = 2.00
electronegativity["Tl"] = 1.62
electronegativity["Pb"] = 2.33
electronegativity["Bi"] = 2.02
electronegativity["Po"] = 2.0
electronegativity["At"] = 2.2
electronegativity["Fr"] = 0.7
electronegativity["Ra"] = 0.89
electronegativity["Ac"] = 1.1
electronegativity["Th"] = 1.3
electronegativity["Pa"] = 1.5
electronegativity["U"] = 1.38
electronegativity["Np"] = 1.36
electronegativity["Pu"] = 1.28
electronegativity["Am"] = 1.3
electronegativity["Cm"] = 1.3
electronegativity["Bk"] = 1.3
electronegativity["Cf"] = 1.3
electronegativity["Es"] = 1.3

end #ends module


text:
using Documenter, Genie
using App, AppServer
using Cache, Commands, Configuration, Cookies, WebChannels
using FileTemplates, Generator
using Helpers, Inflector, Input, Loggers
using Renderer, REPL, Router, Sessions, Tester, Toolbox, Util
using FileCacheAdapter, FileSessionAdapter, Encryption

push!(LOAD_PATH,  "../../src",
                  "../../src/cache_adapters",
                  "../../src/session_adapters")

makedocs()


text:
module Make_AEM_Synth
using PyPlot, GP, Random, UseGA_AEM, TransD_GP

function geomprogdepth(n, dy, c)
    dy*(1.0-c^n)/(1-c)
end

function getn(z, dy, c)
    log(1 - z/dy*(1-c))/log(c)
end

function main()
    n, dz, extendfrac  = 52, 1.1, 1.05
    znrange            = 1.0:n
    zboundaries        = geomprogdepth.(znrange, dz, extendfrac)
    thickness          = [zboundaries[1]; diff(zboundaries)[1:end-1]]
    zall               = [zboundaries[1]/2; 0.5*(zboundaries[1:end-1] + zboundaries[2:end])]
    znall              = getn.(zall, dz, extendfrac)

    figure()
    plot(znall, zall)
    xlabel("depth index")
    ylabel("depth associated")
    grid()
    nicenup(gcf())

    f, ax = plt.subplots(1, 2, figsize=(10,5))
    ax[1].stem(zboundaries[1:end-1], zboundaries[1:end-1], markerfmt="")
    ax[1].stem(zall, zall, "k--", markerfmt=" ")
    ax[1].set_xlabel("depth m")
    ax[1].set_ylabel("depth m")

    ax[2].stem(znrange[1:end-1], znrange[1:end-1], markerfmt="")
    ax[2].stem(znall, znall, "k--", markerfmt=" ")
    ax[2].set_ylabel("depth index")
    ax[2].set_xlabel("depth index")
    nicenup(gcf())

    f, ax = plt.subplots(1, 2, figsize=(10,5), sharey=true)
    ax[1].stem(zall[1:end-1],thickness, "k--", markerfmt=" ")
    ax[1].set_ylabel("thickness m")
    ax[1].yaxis.grid(which="major")
    ax[2].stem(znall[1:end-1],thickness, "k--", markerfmt=" ")
    ax[2].set_xlabel("depth index")
    ax[2].yaxis.grid(which="major")
    nicenup(f)

    λ, δ = [2.0], 0.01
end

function nicenup(g::PyPlot.Figure;fsize=14)
    for ax in gcf().axes
        ax.tick_params("both",labelsize=fsize)
        ax.xaxis.label.set_fontsize(fsize)
        ax.yaxis.label.set_fontsize(fsize)
        ax.title.set_fontsize(fsize)
        if typeof(ax.get_legend_handles_labels()[1]) != Array{Any,1}
            ax.legend(loc="best", fontsize=fsize)
        end
    end
    g.tight_layout()
end

end


text:
"""
A highly optimized Business Days calculator written in Julia language. Also known as Working Days calculator.

Website: https://github.com/JuliaFinance/BusinessDays.jl
"""
module BusinessDays

import Dates

# exported types
export
    HolidayCalendar,
    CompositeHolidayCalendar,
    GenericHolidayCalendar

# exported functions
export
    isholiday,
    isweekday,
    isweekend,
    isbday,
    tobday,
    advancebdays,
    bdayscount,
    bdays,
    firstbdayofmonth,
    lastbdayofmonth,
    listholidays,
    listbdays

include("dateutils.jl")
include("holidaycalendar.jl")
include("bdayscache.jl")
include("bdays.jl")
include("bdaysvecfun.jl")
include("composite.jl")
include("query.jl")
include("calendars/calendars.jl")
include("generic.jl")

end # module BusinessDays


text:
function discretise(dist::ContinuousUnivariateDistribution, n::Integer=10)
    F(t) = cdf(dist, t)
    lb = minimum(dist)
    ub = maximum(dist)
    nodes = range(lb, ub; length = n + 1)
    pmf = [F(nodes[i + 1]) - F(nodes[i]) for i in 1:n]
    return DiscreteNonParametric(nodes[2:end], pmf)
end

discretise(dist::ContinuousUnivariateDistribution, n::Number) = discretise(dist, Integer(n))

max_stars_repo_name: carlocab/BPPrivInfo

max_issues_repo_name: carlocab/BPPrivInfo

max_forks_repo_name: carlocab/BPPrivInfo


text:
# We have n particles here.
# For each particle i (i ∈ N⁺)  :
# We have Mₘᵢₙ ≤ mᵢ ≤ Mₘₐₓ
# We have an three dimentional vector v⃗ᵢ for i's velocity.
# And an three dimentional point Pᵢ for i's locationn.
# We have Gravitational Constant G ≈ 6.67408*10⁻¹¹ (m³kg⁻¹s⁻²) .
# And a small time step Δt→0 (actually, it's a very small value).
# At the beginning, every particle will spread in an area randomly.
# At the beginning, every particle will have a velocity which direction and size are all given randomly but in a range.
#

# Modules
using LinearAlgebra

# Constants
n = 1000
G = 6.67408e-11
Δt = 60*12
Δtₕ = Δt/2
Mₘₐₓ = 1e22
Mₘᵢₙ = 1e20
Areaₘᵢₙ = [-1e11 -1e11 -1e11]
Areaₘₐₓ = [1e11 1e11 1e11]
Areaₜₑₘₚ = Areaₘₐₓ - Areaₘᵢₙ
v⃗ₘᵢₙ  = [-5e2 -5e2 -5e2]
v⃗ₘₐₓ  = [5e2 5e2 5e2]
v⃗ₜₑₘₚ  = v⃗ₘₐₓ  - v⃗ₘᵢₙ
SavePer = 5
Total = 1024*16

# Generate
MGΔt = (rand(n)*(Mₘₐₓ-Mₘᵢₙ) + ones(n)*Mₘᵢₙ)*G*Δt
P₀ = zeros(3, n)
P₁ = zeros(3, n)
v⃗  = zeros(3, n)
for i=1:3
  P₀[i, :] = rand(n)*Areaₜₑₘₚ[i] + ones(n)*Areaₘᵢₙ[i]
  v⃗[i,  :] = rand(n)*v⃗ₜₑₘₚ[i]  + ones(n)*v⃗ₘᵢₙ[i]
end

#Set an blackhole
P₀[:, 1] = [0.0 0.0 0.0]
v⃗[:, 1] = [0.0 0.0 0.0]
MGΔt[1] = 1e30*G*Δt

# Functions
#rsqrt(x) = 1/sqrt(x);
function update(P₀, P₁, v⃗)
  for i=1:n
    Δta⃗  = [0.0, 0.0, 0.0]
    for j=1:(i-1)
      Δta⃗  += MGΔt[j]/(norm(P₀[:, j] - P₀[:, i])^3)*(P₀[:, j] - P₀[:, i])
    end
    for j=(i+1):n
      Δta⃗  += MGΔt[j]/(norm(P₀[:, j] - P₀[:, i])^3)*(P₀[:, j] - P₀[:, i])
    end
    P₁[:, i] = P₀[:, i] + v⃗[:, i]*Δt + Δta⃗ *Δtₕ
    v⃗[:, i] = v⃗[:, i] + Δta⃗
  end
end

# Run
for j=0:Total-1
  @time if mod(j, 2)==0
    update(P₀, P₁, v⃗)
  else
    update(P₁, P₀, v⃗)
  end
  if mod(j, SavePer)==0
    println(P₀[:, 1])
  end
end

max_stars_repo_name: RESBI/N-Body-in-Julia

max_issues_repo_name: RESBI/N-Body-in-Julia

max_forks_repo_name: RESBI/N-Body-in-Julia


text:
function Base.show(io::IO, dfr::DataFrameRow;
                   allcols::Bool = !get(io, :limit, false),
                   rowlabel::Symbol = :Row,
                   eltypes::Bool = true,
                   truncate::Int = 32,
                   kwargs...)
    r, c = parentindices(dfr)
    _show(io, view(parent(dfr), [r], c); allcols=allcols, rowlabel=rowlabel,
          summary=false, rowid=r, eltypes=eltypes, truncate=truncate,
          title="DataFrameRow", kwargs...)
end

Base.show(io::IO, mime::MIME"text/plain", dfr::DataFrameRow;
          allcols::Bool = !get(io, :limit, false),
          rowlabel::Symbol = :Row,
          eltypes::Bool = true,
          truncate::Int = 32,
          kwargs...) =
    show(io, dfr; allcols=allcols, rowlabel=rowlabel, eltypes=eltypes,
         truncate=truncate, kwargs...)

Base.show(dfr::DataFrameRow;
          allcols::Bool = !get(stdout, :limit, true),
          rowlabel::Symbol = :Row,
          eltypes::Bool = true,
          truncate::Int = 32,
          kwargs...) =
    show(stdout, dfr; allcols=allcols, rowlabel=rowlabel, eltypes=eltypes,
         truncate=truncate, kwargs...)


text:
cd(Pkg.dir("uJSON", "deps"))
run(`make`)


text:
# typedef int Menu_Options;
# typedef int Item_Options;
#
# /* Menu options: */
# #define O_ONEVALUE      (0x01)
# #define O_SHOWDESC      (0x02)
# #define O_ROWMAJOR      (0x04)
# #define O_IGNORECASE    (0x08)
# #define O_SHOWMATCH     (0x10)
# #define O_NONCYCLIC     (0x20)
#
# /* Item options: */
# #define O_SELECTABLE    (0x01)
#
# typedef struct
# {
#   const char* str;
#   unsigned short length;
# } TEXT;
#
# typedef struct tagITEM
# {
#   TEXT           name;        /* name of menu item                         */
#   TEXT           description; /* description of item, optional in display  */
#   struct tagMENU *imenu;      /* Pointer to parent menu                    */
#   void           *userptr;    /* Pointer to user defined per item data     */
#   Item_Options   opt;         /* Item options                              */
#   short          index;       /* Item number if connected to a menu        */
#   short          y;           /* y and x location of item in menu          */
#   short          x;
#   bool           value;       /* Selection value                           */
#
#   struct tagITEM *left;       /* neighbor items                            */
#   struct tagITEM *right;
#   struct tagITEM *up;
#   struct tagITEM *down;
#
# } ITEM;
#
# typedef void (*Menu_Hook)(struct tagMENU *);
#
# typedef struct tagMENU
# {
#   short          height;                /* Nr. of chars high               */
#   short          width;                 /* Nr. of chars wide               */
#   short          rows;                  /* Nr. of items high               */
#   short          cols;                  /* Nr. of items wide               */
#   short          frows;                 /* Nr. of formatted items high     */
#   short          fcols;                 /* Nr. of formatted items wide     */
#   short          arows;                 /* Nr. of items high (actual)      */
#   short          namelen;               /* Max. name length                */
#   short          desclen;               /* Max. description length         */
#   short          marklen;               /* Length of mark, if any          */
#   short          itemlen;               /* Length of one item              */
#   short          spc_desc;              /* Spacing for descriptor          */
#   short          spc_cols;              /* Spacing for columns             */
#   short          spc_rows;              /* Spacing for rows                */
#   char          *pattern;               /* Buffer to store match chars     */
#   short          pindex;                /* Index into pattern buffer       */
#   WINDOW        *win;                   /* Window containing menu          */
#   WINDOW        *sub;                   /* Subwindow for menu display      */
#   WINDOW        *userwin;               /* User's window                   */
#   WINDOW        *usersub;               /* User's subwindow                */
#   ITEM          **items;                /* array of items                  */
#   short          nitems;                /* Nr. of items in menu            */
#   ITEM          *curitem;               /* Current item                    */
#   short          toprow;                /* Top row of menu                 */
#   chtype         fore;                  /* Selection attribute             */
#   chtype         back;                  /* Nonselection attribute          */
#   chtype         grey;                  /* Inactive attribute              */
#   unsigned char  pad;                   /* Pad character                   */
#
#   Menu_Hook      menuinit;              /* User hooks                      */
#   Menu_Hook      menuterm;
#   Menu_Hook      iteminit;
#   Menu_Hook      itemterm;
#
#   void          *userptr;               /* Pointer to menus user data      */
#   char          *mark;                  /* Pointer to marker string        */
#
#   Menu_Options   opt;                   /* Menu options                    */
#   unsigned short status;                /* Internal state of menu          */
#
# } MENU;
#
#
# /* Define keys */
#
# #define REQ_LEFT_ITEM           (KEY_MAX + 1)
# #define REQ_RIGHT_ITEM          (KEY_MAX + 2)
# #define REQ_UP_ITEM             (KEY_MAX + 3)
# #define REQ_DOWN_ITEM           (KEY_MAX + 4)
# #define REQ_SCR_ULINE           (KEY_MAX + 5)
# #define REQ_SCR_DLINE           (KEY_MAX + 6)
# #define REQ_SCR_DPAGE           (KEY_MAX + 7)
# #define REQ_SCR_UPAGE           (KEY_MAX + 8)
# #define REQ_FIRST_ITEM          (KEY_MAX + 9)
# #define REQ_LAST_ITEM           (KEY_MAX + 10)
# #define REQ_NEXT_ITEM           (KEY_MAX + 11)
# #define REQ_PREV_ITEM           (KEY_MAX + 12)
# #define REQ_TOGGLE_ITEM         (KEY_MAX + 13)
# #define REQ_CLEAR_PATTERN       (KEY_MAX + 14)
# #define REQ_BACK_PATTERN        (KEY_MAX + 15)
# #define REQ_NEXT_MATCH          (KEY_MAX + 16)
# #define REQ_PREV_MATCH          (KEY_MAX + 17)
#
# #define MIN_MENU_COMMAND        (KEY_MAX + 1)
# #define MAX_MENU_COMMAND        (KEY_MAX + 17)
#
# /*
#  * Some AT&T code expects MAX_COMMAND to be out-of-band not
#  * just for menu commands but for forms ones as well.
#  */
# #if defined(MAX_COMMAND)
# #  if (MAX_MENU_COMMAND > MAX_COMMAND)
# #    error Something is wrong -- MAX_MENU_COMMAND is greater than MAX_COMMAND
# #  elif (MAX_COMMAND != (KEY_MAX + 128))
# #    error Something is wrong -- MAX_COMMAND is already inconsistently defined.
# #  endif
# #else
# #  define MAX_COMMAND (KEY_MAX + 128)
# #endif
#
#
# /* --------- prototypes for libmenu functions ----------------------------- */
#
# extern NCURSES_EXPORT(ITEM **)	menu_items (const MENU *);
# extern NCURSES_EXPORT(ITEM *)	current_item (const MENU *);
# extern NCURSES_EXPORT(ITEM *)	new_item (const char *,const char *);
#
# extern NCURSES_EXPORT(MENU *)	new_menu (ITEM **);
#
# extern NCURSES_EXPORT(Item_Options)	item_opts (const ITEM *);
# extern NCURSES_EXPORT(Menu_Options)	menu_opts (const MENU *);
#
# extern NCURSES_EXPORT(Menu_Hook)	item_init (const MENU *);
# extern NCURSES_EXPORT(Menu_Hook)	item_term (const MENU *);
# extern NCURSES_EXPORT(Menu_Hook)	menu_init (const MENU *);
# extern NCURSES_EXPORT(Menu_Hook)	menu_term (const MENU *);
#
# extern NCURSES_EXPORT(WINDOW *)	menu_sub (const MENU *);
# extern NCURSES_EXPORT(WINDOW *)	menu_win (const MENU *);
#
# extern NCURSES_EXPORT(const char *)	item_description (const ITEM *);
# extern NCURSES_EXPORT(const char *)	item_name (const ITEM *);
# extern NCURSES_EXPORT(const char *)	menu_mark (const MENU *);
# extern NCURSES_EXPORT(const char *)	menu_request_name (int);
#
# extern NCURSES_EXPORT(char *)	menu_pattern (const MENU *);
#
# extern NCURSES_EXPORT(void *)	menu_userptr (const MENU *);
# extern NCURSES_EXPORT(void *)	item_userptr (const ITEM *);
#
# extern NCURSES_EXPORT(chtype)	menu_back (const MENU *);
# extern NCURSES_EXPORT(chtype)	menu_fore (const MENU *);
# extern NCURSES_EXPORT(chtype)	menu_grey (const MENU *);
#
# extern NCURSES_EXPORT(int)	free_item (ITEM *);
# extern NCURSES_EXPORT(int)	free_menu (MENU *);
# extern NCURSES_EXPORT(int)	item_count (const MENU *);
# extern NCURSES_EXPORT(int)	item_index (const ITEM *);
# extern NCURSES_EXPORT(int)	item_opts_off (ITEM *,Item_Options);
# extern NCURSES_EXPORT(int)	item_opts_on (ITEM *,Item_Options);
# extern NCURSES_EXPORT(int)	menu_driver (MENU *,int);
# extern NCURSES_EXPORT(int)	menu_opts_off (MENU *,Menu_Options);
# extern NCURSES_EXPORT(int)	menu_opts_on (MENU *,Menu_Options);
# extern NCURSES_EXPORT(int)	menu_pad (const MENU *);
# extern NCURSES_EXPORT(int)	pos_menu_cursor (const MENU *);
# extern NCURSES_EXPORT(int)	post_menu (MENU *);
# extern NCURSES_EXPORT(int)	scale_menu (const MENU *,int *,int *);
# extern NCURSES_EXPORT(int)	set_current_item (MENU *menu,ITEM *item);
# extern NCURSES_EXPORT(int)	set_item_init (MENU *, Menu_Hook);
# extern NCURSES_EXPORT(int)	set_item_opts (ITEM *,Item_Options);
# extern NCURSES_EXPORT(int)	set_item_term (MENU *, Menu_Hook);
# extern NCURSES_EXPORT(int)	set_item_userptr (ITEM *, void *);
# extern NCURSES_EXPORT(int)	set_item_value (ITEM *,bool);
# extern NCURSES_EXPORT(int)	set_menu_back (MENU *,chtype);
# extern NCURSES_EXPORT(int)	set_menu_fore (MENU *,chtype);
# extern NCURSES_EXPORT(int)	set_menu_format (MENU *,int,int);
# extern NCURSES_EXPORT(int)	set_menu_grey (MENU *,chtype);
# extern NCURSES_EXPORT(int)	set_menu_init (MENU *, Menu_Hook);
# extern NCURSES_EXPORT(int)	set_menu_items (MENU *,ITEM **);
# extern NCURSES_EXPORT(int)	set_menu_mark (MENU *, const char *);
# extern NCURSES_EXPORT(int)	set_menu_opts (MENU *,Menu_Options);
# extern NCURSES_EXPORT(int)	set_menu_pad (MENU *,int);
# extern NCURSES_EXPORT(int)	set_menu_pattern (MENU *,const char *);
# extern NCURSES_EXPORT(int)	set_menu_sub (MENU *,WINDOW *);
# extern NCURSES_EXPORT(int)	set_menu_term (MENU *, Menu_Hook);
# extern NCURSES_EXPORT(int)	set_menu_userptr (MENU *,void *);
# extern NCURSES_EXPORT(int)	set_menu_win (MENU *,WINDOW *);
# extern NCURSES_EXPORT(int)	set_top_row (MENU *,int);
# extern NCURSES_EXPORT(int)	top_row (const MENU *);
# extern NCURSES_EXPORT(int)	unpost_menu (MENU *);
# extern NCURSES_EXPORT(int)	menu_request_by_name (const char *);
# extern NCURSES_EXPORT(int)	set_menu_spacing (MENU *,int,int,int);
# extern NCURSES_EXPORT(int)	menu_spacing (const MENU *,int *,int *,int *);
#
#
# extern NCURSES_EXPORT(bool)	item_value (const ITEM *);
# extern NCURSES_EXPORT(bool)	item_visible (const ITEM *);
#
# extern NCURSES_EXPORT(void)	menu_format (const MENU *,int *,int *);
#
# #if NCURSES_SP_FUNCS
# extern NCURSES_EXPORT(MENU *)	NCURSES_SP_NAME(new_menu) (SCREEN*, ITEM **);
# #endif
#
# #ifdef __cplusplus
#   }
# #endif


text:
## Gegenbauer Polynomials
@registerN Gegenbauer AbstractCCOP1 α
export  Gegenbauer

"""
   Gegenbauer{α, T <: Number}

The Gegenbauer polynomials have weight function
`(1-x^2)^(α-1/2)` over the domain `[-1,1]`. The parameter `α` is
specified in the constructor. These are also called the ultra-spherical polynomials.
The Legendre polynomials are the specialization  `Gegenbauer{1/2}`.

```jldoctest
julia> using Polynomials, SpecialPolynomials

julia> p =  Gegenbauer{1/2}([1,2,3])
typename(Gegenbauer){0.5}(1⋅Cᵅ₀(x) + 2⋅Cᵅ₁(x) + 3⋅Cᵅ₂(x))

julia> convert(Polynomial, p)
Polynomials.Polynomial(-0.5 + 2.0*x + 4.5*x^2)
```

"""
Gegenbauer

basis_symbol(::Type{<:Gegenbauer{α}}) where {α} = "Cᵅ"
Polynomials.domain(::Type{<:Gegenbauer{α}}) where {α} = Polynomials.Interval(-1,1)

abcde(::Type{<:Gegenbauer{α}})  where {α} = NamedTuple{(:a,:b,:c,:d,:e)}((-1,0,1,-(2α+1),0))

k0(P::Type{<:Gegenbauer}) = one(eltype(P))
k1k0(P::Type{<:Gegenbauer{α}}, k::Int) where {α} =(one(eltype(P))*2*(α+k))/(k+1)

function norm2(::Type{<:Gegenbauer{α}}, n) where{α}
    pi * 2^(1-2α) * gamma(n + 2α) / (gamma(n+1) * (n+α) * gamma(α)^2)
end
function ω₁₀(::Type{<:Gegenbauer{α}}, n) where{α}
    val = gamma(n + 1 + 2α) / gamma(n + 2α)
    val /= (n+1)
    val *= (n + α)/(n+1+α)
    sqrt(val)

end

weight_function(::Type{<:Gegenbauer{α}}) where {α} = x -> (1-x^2)^(α-1/2)
generating_function(::Type{<:Gegenbauer{α}}) where {α} = (t,x) -> begin
    1/(1-2x*t +t^2)^α
end
function classical_hypergeometric(::Type{<:Gegenbauer{α}}, n, x) where {α}
    (α > -1/2 && !iszero(α)) || throw(ArgumentError("α > -1/2 and α ≠ 0 is necessary"))

    as = (-n, 2α+n)
    bs = (α + 1/2,)
    return Pochhammer_factorial(2α,n) * pFq(as, bs, (1-x)/2)

end

# overrides; big speedup if we avoid generating  from  a,b,c,d,e
B̃n(P::Type{<:Gegenbauer{α}}, n::Int) where {α} = zero(eltype(P))
C̃n(P::Type{<:Gegenbauer{α}}, n::Int) where {α} = (one(eltype(P))*n*(n + 2*α - 1))/(4*(n^2 + 2*n*α - n + α^2 - α))
b̂̃n(P::Type{<:Gegenbauer{α}}, n::Int) where {α} = zero(eltype(P))# one(S) *  4/α/(α+2)
b̂̃n(P::Type{<:Gegenbauer{α}}, ::Val{N}) where {α,N} = zero(eltype(P))# one(S) *  4/α/(α+2)
ĉ̃n(P::Type{<:Gegenbauer{α}}, ::Val{0}) where {α} = zero(eltype(P)) #one(S) *  4/α/(α+2)

@registerN OrthonormalGegenbauer AbstractCCOP1 α
export OrthonormalGegenbauer
ϟ(::Type{<:OrthonormalGegenbauer{α}}) where {α} = Gegenbauer{α}
ϟ(::Type{<:OrthonormalGegenbauer{α,T}}) where {α,T} = Gegenbauer{α,T}
@register_orthonormal(OrthonormalGegenbauer)


text:
# we write a function that applies Euler's method to solve y' = f(x,y)
# in the interval [x0, x1] with y(x0) = y0 with n steps
function eulerSolve(f, x0, y0, x1, n)

    # we compute the step size h
    h = (x1 - x0)/n

    # we initialize the outpur vectors
    x = zeros(n+1)
    y = zeros(n+1)
    x[1] = x0
    y[1] = y0

    # we interate over the interval in n steps of length h
    for i = 1:n
    	# and calculate x and y accordingly
        y[i+1] = y[i] + h * f(x[i], y[i])
        x[i+1] = x[i] + h
    end

    # we return the data points of the calculated function y(x)
    return (x, y)
end

max_stars_repo_path: 1a/Euler Method/eulerSolve.jl

max_issues_repo_path: 1a/Euler Method/eulerSolve.jl

max_forks_repo_path: 1a/Euler Method/eulerSolve.jl


text:
"""
    _osascript() -> Matrix{<:Colorant}

Paste an image from clipboard using `osascript`
"""
function _osascript()
    mktempdir() do dir
        # Define path
        path_png = joinpath(dir, "clipboard.png")

        # Compose command & run
        cmd = Cmd(["osascript", "-e", "try\nwrite (the clipboard as «class PNGf») to (open for access \"$(path_png)\" with write permission)\nend try"])
        run(cmd)

        # Paste from clipboard
        if isfile(path_png)
            img = load(path_png)
            return img
        else
            error("No image in clipboard")
        end
    end
end

"""
    _osascript(img::Matrix{<:Colorant})

Copy an image to clipboard using `osascript`
"""
function _osascript(img::Matrix{<:Colorant})
    mktempdir() do dir
        # Define path
        path_png = joinpath(dir, "clipboard.png")

        # Save image
        save(path_png, img)

        # Compose command & run
        cmd = Cmd(["osascript", "-e", "set the clipboard to (read \"$(path_png)\" as «class PNGf»)"])
        run(cmd)
    end
end