text:
module ComputationalHomology

using LinearAlgebra
using Hungarian: hungarian
using SmithNormalForm: smith
using Distances: Metric, Euclidean, pairwise, colwise, pairwise!
using Random: shuffle!, randperm
using BoundingSphere: boundingsphere
using SparseArrays: SparseMatrixCSC, spzeros, findnz
using Combinatorics: combinations

import Base: ==, -, +, *, union, keys, values, hash, first, last, isless, show, length,
             eltype, valtype, getindex, setindex!, size, iterate, push!, append!, in, similar,
             read, write, vec, complex, iszero, convert, reduce, keytype, copy, map, map!,
             minimum, maximum, parse, isempty
import Base.Iterators: pairs
import SparseArrays: sparse
import Statistics: mean
import LinearAlgebra: diag

export AbstractCell,
       dim, faces, cofaces, volume,
       AbstractSimplex,
       indices, vertices,
       Simplex, Cube, Cell,

       AbstractChain,
       Chain, simplify,

       AbstractComplex,
       boundary, coboundary, cells,
       SimplicialComplex, addsimplex!, addsimplices!,
       CWComplex, BitmapComplex,
       vietorisrips, witness, cech, čech,

       AbstractHomology, grouptype, group,
       Homology, homology, withgenerators, generators,
       betti, euler,

       Filtration, filtration, order, simplices, similarity,

       AbstractPersistenceReduction,
       StandardReduction, TwistReduction,
       pairs, reduce!,
       PersistentHomology, persistenthomology,
       PersistentCocycleReduction, persistentcohomology, PersistentCocycleReduction,

       AbstractInterval, birth, death,
       Interval, AnnotatedInterval,
       PersistenceDiagram, diagram,
       Landscape, landscape, mean,
       PersistentImage,
       wasserstein,
       tropic

include("abstractchain.jl")
include("chains.jl")
include("abstractcell.jl")
include("simplex.jl")
include("cube.jl")
include("cwcell.jl")
include("complex.jl")
include("simplicialcomplex.jl")
include("cwcomplex.jl")
include("constructions.jl")
include("filtration.jl")
include("bitmapcomplex.jl")
include("homology.jl")
include("intervals.jl")
include("persistence.jl")
include("landscape.jl")
include("pimage.jl")
include("examples.jl")
include("distances.jl")
include("iterators.jl")
include("tropic.jl")

@deprecate intervals(d, ps...) diagram(d, ps...)
@deprecate celltype(cplx) eltype(cplx)

end # module

ext: jl


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

@testset "TestCompat" begin

#---

using JuLIP, ACE
import ZipFile, JSON

include("../artifacts.jl")
testdatadir = joinpath(datadir, "tests")

#---

@info("Test auto-generated testset `v07_compat`")
zipname = joinpath(testdatadir, "v07_compat.zip")
zipdir = ZipFile.Reader(zipname)

for fptr in zipdir.files
   print("    $(fptr.name): ")
   D = JSON.parse(fptr)
   calc = read_dict(D)
   ACE.Testing.runtests(calc, D["_tests"])
   if fptr.name[1:8] == "rpibasis"
      basis = ACE.Testing.test_basis(D["_params"])
      print(" / ")
      ACE.Testing.runtests(basis, D["_tests"])
   end
   println()
end

close(zipdir)

#---

@info("Test auto-generated testset `v08_compat`")
zipname = joinpath(testdatadir, "v08_compat.zip")
zipdir = ZipFile.Reader(zipname)

for fptr in zipdir.files
   print("    $(fptr.name): ")
   D = JSON.parse(fptr)
   calc = read_dict(D)
   ACE.Testing.runtests(calc, D["_tests"])
   if fptr.name[1:8] == "rpibasis"
      basis = ACE.Testing.test_basis(D["_params"])
      print(" / ")
      ACE.Testing.runtests(basis, D["_tests"])
   end
   println()
end

close(zipdir)

end


text:
# Client for FedProx with linear local model
mutable struct FedProxClient{T1<:Int64, T2<:Float64, T3<:SparseMatrixCSC{Float64, Int64}, T4<:Matrix{Float64}, T5<:Vector{Int64}, T6<:Function} <: AbstractClient
    id::T1                                  # client index
    Xtrain::T3                              # training data
    XtrainT::T3                             # Row copy
    Ytrain::T5                              # training label
    W::T4                                   # (model) primal variable
    lr::T2                                  # learning rate
    λ::T2                                   # regularization parameter
    μ::T2                                   # proximal parameter
    numLocalEpochs::T1                      # number of local steps
    function FedProxClient(id::Int64, Xtrain::SparseMatrixCSC{Float64, Int64}, Ytrain::Vector{Int64}, config::Dict{String,Real})
        num_classes = config["num_classes"]
        λ = config["lambda"]
        learning_rate = config["learning_rate"]
        numLocalEpochs = config["numLocalEpochs"]
        μ = config["mu"]
        d = size(Xtrain, 2)
        W = zeros(Float64, d, num_classes)
        XtrainT = copy(Xtrain')
        new{Int64, Float64, SparseMatrixCSC{Float64, Int64}, Matrix{Float64}, Vector{Int64}, Function}(id, Xtrain, XtrainT, Ytrain, W, learning_rate, λ, μ, numLocalEpochs)
    end
end

# Model update on local device
function update!(
    client::FedProxClient
)
    # Implement K epochs SGD, using lazy update to avoid dense update from the regularization.
    n, d = size(client.Xtrain)
    _, K = size(client.W)
    lr = client.lr
    Wold = copy(client.W)
    λ = client.λ
    μ = client.μ
    hitTime = zeros(Int64, d, K)
    perm = collect(1:n)
    for epoch = 1:client.numLocalEpochs
        fill!(hitTime, 0)
        shuffle!(perm)
        timeStep = 1
        for i in perm
            g = getStochasticGrad(client.XtrainT, client.Ytrain, client.W, i)
            I, J, V = findnz(g)
            for j = 1:length(I)
                idx1 = I[j]
                idx2 = J[j]
                # Lazy update
                delay = timeStep-hitTime[idx1, idx2]
                client.W[idx1, idx2] *= (1 - lr*(λ+μ))^delay
                client.W[idx1, idx2] += lr*μ*delay*Wold[idx1, idx2]
                client.W[idx1, idx2] -= lr*V[j]
                # Update hitTime to the current time
                hitTime[idx1, idx2] = timeStep
            end
            timeStep += 1
        end
        # Lazy update for staled coordinates
        timeStep -= 1
        for j = 1:d
            for k = 1:K
                if hitTime[j, k] < timeStep
                    delay = timeStep-hitTime[j, k]
                    client.W[j, k] *= (1 - lr*λ)^delay
                    client.W[j, k] += lr*μ*delay*Wold[j, k]
                end
            end
        end
    end
    return nothing
end

# Get objective value
function getObjValue(
    client::FedProxClient
)
    objValue = obj(client.Xtrain, client.Ytrain, client.W, client.λ)
    return objValue
end

# Client for FedProx with multilayer perceptron local model
mutable struct FedProxClientNN{T1<:Int64, T2<:Float64, T3<:SparseMatrixCSC{Float64, Int64}, T4<:Zygote.Params, T5<:Flux.OneHotArray, T6<:Function, T7<:Flux.Chain} <: AbstractClient
    id::T1                                  # client index
    Xtrain::T3                              # training data
    XtrainT::T3                             # Row copy
    Ytrain::T5                              # training label
    W::T7                                   # (model) primal variable
    lr::T2                                  # learning rate
    λ::T2                                   # regularization parameter
    μ::T2                                   # proximal parameter
    numLocalEpochs::T1                      # number of local steps
    function FedProxClientNN(id::Int64, Xtrain::SparseMatrixCSC{Float64, Int64}, Ytrain::Vector{Int64}, dim::Int64, config::Dict{String,Real})
        numClasses = config["num_classes"]
        Ytrain = Flux.onehotbatch(Ytrain, 1:numClasses)
        λ = config["lambda"]
        learning_rate = config["learning_rate"]
        numLocalEpochs = config["numLocalEpochs"]
        μ = config["mu"]
        W = Chain( Dense(780, dim, relu, bias=false), Dense(dim, 10, bias=false), NNlib.softmax)
        XtrainT = copy(Xtrain')
        new{Int64, Float64, SparseMatrixCSC{Float64, Int64}, Zygote.Params, Flux.OneHotArray, Function, Flux.Chain}(id, Xtrain, XtrainT, Ytrain, W, learning_rate, λ, μ, numLocalEpochs)
    end
end

# Model update on local device
function update!(
    client::FedProxClientNN
)
    sgd!(client.XtrainT, client.Ytrain, client.W, client.λ, client.μ, client.lr, client.numLocalEpochs)
end


text:
function center_test(x1, x2, x3, y1, y2, y3)
    arr1 = randn((x1, x2, x3))
    arr2 = zeros((y1, y2, y3))

    FFTResampling.center_set!(arr2, arr1)
    arr3 = FFTResampling.center_extract(arr2, (x1, x2, x3))
    @test arr1 == arr3
end

 # test center set and center extract methods
@testset "center methods" begin
    center_test(4, 4, 4, 6,7,4)
    center_test(5, 4, 4, 7, 8, 4)
    center_test(5, 4, 4, 8, 8, 8)
    center_test(6, 4, 4, 7, 8, 8)

    @test 1 == FFTResampling.center_pos(1)
    @test 2 == FFTResampling.center_pos(2)
    @test 2 == FFTResampling.center_pos(3)
    @test 3 == FFTResampling.center_pos(4)
    @test 3 == FFTResampling.center_pos(5)
    @test 513 == FFTResampling.center_pos(1024)

    @test FFTResampling.get_indices_around_center((5), (2)) == (2, 3)
    @test FFTResampling.get_indices_around_center((5), (3)) == (2, 4)
    @test FFTResampling.get_indices_around_center((4), (3)) == (2, 4)
    @test FFTResampling.get_indices_around_center((4), (2)) == (2, 3)
end

@testset "slice" begin

    x = randn((1,2,3,4))
    y = FFTResampling.slice(x, 2, 2)
    @test x[:, 2:2, :, :] == y

    x = randn((5,2,3,4))
    y = FFTResampling.slice(x, 1, 4)
    @test x[4:4, :, :, :] == y

    x = randn((5))
    y = FFTResampling.slice(x, 1, 5)
    @test x[5:5] == y

end

@testset "slice indices" begin
    x = randn((1,2,3))
    y = FFTResampling.slice_indices(x, 1, 1)
    @test y == [1:1, 1:2, 1:3]

    x = randn((20,4,20, 1, 2))
    y = FFTResampling.slice_indices(x, 2, 3)
    @test y == [1:20, 3:3, 1:20, 1:1, 1:2]
end

@testset "reverse function" begin
    reverse_all = FFTResampling.reverse_all
    reverse_all_indices = FFTResampling.reverse_all_indices
    function test_reverse(x, y)
        x_rev = reverse_all(x)
        @test x == reverse_all(x_rev)
        @test x_rev == y

        @test x_rev == x[reverse_all_indices(x)...]
    end

    test_reverse([1,2,3], [3,2,1])
    test_reverse([1 2; 3 4], [4 3; 2 1])

    x = randn((2,3,4,1, 5,6))

    x_rev = x[end:-1:1, end:-1:1,end:-1:1,end:-1:1,end:-1:1, end:-1:1]
    test_reverse(x, x_rev)

    y = [3:-1:1, 5:-1:1, 1:-1:1, 20:-1:1]
    @test y == reverse_all_indices(randn((3,5,1,20)))
end

 # Just a reference implementation I don't want to delete
 # is slower, but works
@gg function make_hermitian_ref(arr::AbstractArray{T, N}) where {T, N}
    N2 = N - 1
    quote
        # if the size is odd, we need to add a slice to that dimension
        # at this slice we then fill new values, so that `arr` is hermitian
        mf(x) = x % 2 == 0 ? x +1 : x
        size_new = map(mf, size(arr))
        # new array with new size
        arr_new = zeros(eltype(arr), size_new)

        # fill the new array with the old one
        same_ind = map(x -> 1:x, size(arr))
        arr_new[same_ind...] = arr

        # copy again to have a out array (might be necessary because of reassigning)
        arr_out = copy(arr_new)

        # now modify values at the edges to get hermitian property
        # and Parsveval's theorem correct
        for d = 1:$N
            # special case when singleton dimension
            if size(arr_new)[d] == 1 || size(arr)[d] % 2 == 1
                continue
            end

            # the N2 dimensional loops, loop over this
            iter_nloops = map(x -> 1:x, collect(size_new))
            # delete at dimension d because we have only N2 dimensions to loop over
            deleteat!(iter_nloops, d)

            # @show iter_nloops
            @nloops $N2 i (k -> iter_nloops[k]) begin
                # indices to access the left slice
                # left slice is the existing slice which will be copied
                inds_l = convert(Array{Int},collect(@ntuple $N2 i))
                insert!(inds_l, d, 1)

                # right slice which will be copied from the left slice
                inds_r = similar(inds_l)

                # because of the hermitian property, we reverse the slices
                for j in 1:length(inds_r)
                    # in this case, we want to get the right hand slice
                    if j == d
                        inds_r[j] = size_new[j]
                    else
                        inds_r[j] = size_new[j] - inds_l[j] + 1
                    end
                end

                # assign the complex conjugate
                arr_out[inds_l...] = 0.5 * arr_new[inds_l...]
                arr_out[inds_r...] = 0.5 * conj(arr_new[inds_l...])
            end

        end
        return arr_out
    end
end

@testset "make hermitian" begin
    a = FFTResampling.make_hermitian([1im 2.0; 3.0im 4.0im; 5.0 6.0im])

    ã = [0.0+0.5im  2.0+0.0im  2.5-0.0im; 0.0+1.5im  0.0+4.0im  0.0-1.5im; 2.5+0.0im  0.0+6.0im  0.0-0.5im]

    @test a ≈ ã

    a = FFTResampling.make_hermitian([1.0; 2.0; 3.0im])
    ã = [1.0; 2.0; 3.0im]
    @test a ≈ ã

    a = FFTResampling.make_hermitian([1.0+1.0im; 2.0; 3.0; 3.0im])
    ã = [0.5+0.5im; 2.0; 3.0; 3.0im; 0.5-0.5im]
    @test a ≈ ã

    function test_sum(s)
        a = randn(s)
        ã = FFTResampling.make_hermitian(a)
        @test sum(a) ≈ sum(ã)
    end
    test_sum((12, 3, 12, 11, 1))
    test_sum((123, 312))
    test_sum((1))
    test_sum((2))
    test_sum((100, 101))
    test_sum((101, 101))
    test_sum((101, 100))
    test_sum((101, 101))

    function test_symmetry(s)
        a = randn(s)
        b = FFTResampling.make_hermitian(a)

        for k1 = 1:size(b)[1]
            for k2 = 1:size(b)[2]
                for k3 = 1:size(b)[3]
                    if k1 != 1 && k2 != 1 && k3 != 1
                        continue
                    end
                    @test b[k1, k2, k3] ≈
                        conj(b[size(b)[1] + 1 - k1, size(b)[2] + 1 - k2, size(b)[3] + 1 - k3])
                end
            end
        end
    end

    test_symmetry((2,2,2))
    test_symmetry((1,1,1))
    test_symmetry((4,4,4))
    test_symmetry((12,8,12))

    function test_symmetry_4D(s)
        a = randn(s)
        b = FFTResampling.make_hermitian(a)

        for k1 = 1:size(b)[1]
            for k2 = 1:size(b)[2]
                for k3 = 1:size(b)[3]
                    for k4 = 1:size(b)[4]
                        if k1 != 1 && k2 != 1 && k3 != 1 && k4 != 1
                            continue
                        end
                        @test b[k1, k2, k3, k4] ≈ conj(b[size(b)[1] + 1 - k1, size(b)[2] + 1 - k2, size(b)[3] + 1 - k3, size(b)[4] + 1 - k4])
                    end
                end
            end
        end
    end

    test_symmetry_4D((2,4, 6, 8))

    function compare_ref(x)
        a = FFTResampling.make_hermitian(x)
        b = make_hermitian_ref(x)
        @test a ≈ b
    end

    compare_ref(randn((2, 1, 2,3)))
    compare_ref(randn((2, 1, 2,3, 2, 4, 8)))
    compare_ref(randn((1022, 4)))
    compare_ref(randn((1022, 1022)))
    compare_ref(randn((5, 4)))
end

@testset "test 1D dft" begin
    x = randn(1)
    @test fft(x) ≈ FFTResampling.dft_1D(x)

    x = randn(122)
    @test fft(x) ≈ FFTResampling.dft_1D(x)

    x = randn(123)
    @test fft(x) ≈ FFTResampling.dft_1D(x)
end


text:
#
# Correctness Tests
#
using Test

# Hack: force the loading of MbedTLS as otherwise it conflicts with the version in LibCURL
# Entropy() was chosen as a fairly quick exported function
using MbedTLS
MbedTLS.Entropy()

using XMLDict
using BioMedQuery.Processes
using BioMedQuery.PubMed
using BioMedQuery.DBUtils
using BioMedQuery
using MySQL
using SQLite
using DataStreams

#For now this corresponds to JULIACIBot... since we aren't testing anywhere else
global CI = get(ENV, "CI", "false")=="true"
global TRAVIS = get(ENV, "TRAVIS", "false")=="true"
println("CI = ", CI, ", TRAVIS = ", TRAVIS)
global CI_SKIP_MYSQL = false

if CI == true && TRAVIS == false
     CI_SKIP_MYSQL = true
     @warn("MySQL tests not running")
end

all_tests = [
            ("dbutils_sqlite.jl",   "       Testing: DBUtils SQLite"),
            ("pubmed.jl",           "       Testing: Eutils/PubMed"),
            ("export_citations.jl", "       Testing: Export Citations"),
            ("pubmed_parse.jl",     "       Testing: Entrez Parsing"),
            ("ct.jl",               "       Testing: CLINICAL TRIALS"),
            ("processes_sqlite.jl", "       Testing: Processes SQLite"),
            ("processes_df.jl",     "       Testing: Processes DataFrame")
            ]
if !CI_SKIP_MYSQL
    # on MacOS host name must be localhost
    if Sys.isapple()
        global const MYSQL_HOST = "localhost"
    else
        global const MYSQL_HOST = "127.0.0.1"
    end

    if Sys.isapple()
        global const MYSQL_USER = get(ENV, "MYSQL_USER", get(ENV, "LOGNAME", "root"))
    else
        global const MYSQL_USER = get(ENV, "MYSQL_USER", "root")
    end

    global const MYSQL_PASSWORD = get(ENV, "MYSQL_PASSWORD", "")

    push!(all_tests, ("dbutils_mysql.jl",    "       Testing: DBUtils MySQL"))
    push!(all_tests, ("processes_mysql.jl",  "       Testing: Processes MySQL"))
end
println("Running tests:")

for (test, test_string) in all_tests
    println("-----------------------------------------")
    println("-----------------------------------------")
    println(test_string)
    println("-----------------------------------------")
    println("-----------------------------------------")

    include(test)
end


text:
# window.jl: contains numerical test and windowed version of our simple 2D BEM
# About
#   Author       	- Peter Opsomer (Peter.Opsomer@cs.kuleuven.be)
#   History      	- Created December 2014
# References
#   Smooth cutoff function	- Huybrechs, Daan, "Phase extraction without phase extraction: asymptotic compression of dense BEM matrices", talk on high frequency scattering workshop, Reading, 16/12/2014.
include("bea.jl")

# Test the solution on BC, Helmholtz and Sommerfeld and save some density values
# Input
#	xyz	- Mesh
#	elmt	- Element connectivity
#	k	- Wave number
#   	v	- Known density
function checkSolution(xyz,elmt,k,c)
	rests = zeros(6,7)
	mrt = zeros(6,7)
	hh = 1e-5
	cdf = zeros(Complex,3,2)
	sommer = zeros(6,7)
	for ti = 1:size(sommer,1), ppi = 1:size(sommer,2)
		theta = ti*pi/(size(sommer,2)-1)
		phi = ppi*pi/(size(sommer,1)/2)
		x = sin(theta)*cos(phi)
		y = sin(theta)*sin(phi)
		z = cos(theta)
		rests[ti,ppi] = abs( (exp(1im*k*x)+potential(xyz,elmt,k,[x y z],c))/exp(1im*k*x))
		# Dont take pt on boundary when checking central differences
		pt = [x y z]*1.2
		potpt = exp(1im*k*pt[1])+potential(xyz,elmt,k,pt,c)
		for dim = 1:3, si = -1:2:1
			pt[dim] = pt[dim] + si*hh
			cdf[dim,int((si+3)/2)] = exp(1im*k*pt[1])+potential	(xyz,elmt,k,pt,c)
			pt[dim] = pt[dim] - si*hh
		end
		mrt[ti,ppi] = abs( (( (sum(sum(cdf,1),2)-6*potpt)[1])/hh^2+k^2*potpt)/(k^2*potpt))
		if mrt[ti,ppi] > 1
			print("Err ", ti, ppi, ", potpt and cdf are ", potpt, cdf, ", su = ", ( (sum(sum(cdf,1),2)-6*potpt)[1]), '\n')
		end
		pt = [sin(theta) sin(theta) cos(theta)]
		sommer[ti,ppi] = abs(potential(xyz,elmt,k,2^ppi*pt,c) )*2^ppi
	end
	print(norm(rests), "is Rel. error on boundary conditions and Rel. errors on central differences and going away are ", norm(mrt), sommer[3,:], '\n')

	nPointsTh = 4
	nPointsPh = 4
	theta = linspace(0.1,pi-0.1,nPointsTh)
	phi = linspace(pi/2,3*pi/2-0.1,nPointsPh)

	results= zeros(Complex,3+nPointsTh^2)
	results[1:3] = interpolateGrid( ([1,-1,-1/sqrt(3)],[0,0,-1/sqrt(3)], [0,0,1/sqrt(3)]), c, xyz)

	f = open(string("resultsBEJ", nrMesh, "k", k), "w")
	writedlm(f,k)
	write(f,"was k, ansatz is x, x&y&z are \n")
	writedlm(f,[1 0 0; -1 0 0; -1/sqrt(3) -1/sqrt(3) 1/sqrt(3)])
	for nt = 1:nPointsTh
		x = sin(theta[nt])*cos(phi)
		y = sin(theta[nt])*sin(phi)
		z = cos(theta[nt])*ones(size(phi))
		results[(4+(nt-1)*nPointsTh):(3+(nt)*nPointsTh)] = interpolateGrid((x, y, z ), c, xyz)
		writedlm(f,[x y z])
	end
	write(f,"were x&y&z, results are \n")
	writedlm(f,[real(results) imag(results)])
	close(f)
end

# Apply window function
# Input
#	r	- Distance
#	m	- Cutoff
# Output
#		- 0 if r>m, 1 if r<0.7*m and else C^\infty-smooth
function window(r,m)
	if (r <= 0.7*m) return 1;
	elseif (r < m) return exp(2*exp(-(0.3*m)/(r-0.7*m) )/( (r-0.7*m)/(0.3*m)-1) );
	else return 0;
	end
end

# Make the compressed system matrix for triangular elements.
# Input
#   xyz		- Mesh
#   elmts 	- Element connectivity
#   k 		- Wave number
#   m		- Cutoff
# Output
#   A		- The compressed system matrix
function getComprSystemMatrix(xyz, elmts, k, m)
	nnode = size(xyz,1)
	A = zeros(Complex{Float64}, nnode, nnode)
	for ni = 1:nnode
		if (ni%10 == 0)
			print(ni/nnode, " ")
		end
		if(xyz[ni,1] < -0.2) # Illuminated region
			updateComprA(xyz,elmts,k,ni,A,m);
		else
			updateA(xyz,elmts,k,ni,A);
		end
	end
	return A
end

# Update row of the compressed system matrix
# Input
#   xyz		- Mesh
#   elmts 	- Element connectivity
#   k 		- Wave number
#   ni		- Row to update
#   A		- Compressed system matrix
#   m		- Cutoff
function updateComprA(xyz, elmts,k, ni, A, m)
	nelem = size(elmts,1)
	xk = [xyz[ni,:]; xyz[ni,:]; xyz[ni,:]]
	for ei = 1:nelem
		xm = xyz[vec(elmts[ei,:]),:]
		# Sort so that third node of the element (xi1=0=xi2) is closest to or on singularity (xm = xyz[ni,:])
		(val,idx) = findmin(sqrt( sum((xm-xk).^2,2) ) )
		indir = [1, 2, 3]
		indir[idx] = 3
		indir[3] = idx
		J = jacob(xm[indir,1], xm[indir,2], xm[indir,3])/4
		for qi = 1:3, ri=1:3, si = -1:2:1, ti = -1:2:1
			dist = norm(triang([(1+si*qno[qi])/2*(1-(1+ti*qno[ri])/2), (1+si*qno[qi])*(1+ti*qno[ri])/4],xm[indir,:]) -xk[1,:])
			common = wei[qi]*wei[ri]*exp(1im*k*dist)/(4*pi)/sqrt(sum( (xm[indir[1],:]+(1+ti*qno[ri])/2*(xm[indir[2],:]-xm[indir[1],:]) -xm[indir[3],:]+(xm[indir[3],:]-xk[1,:])/(1+si*qno[qi])*2).^2,2)[1] )*J*window(dist,m)

			A[ni,elmts[ei,indir[1]] ] += common*((1+si*qno[qi])/2*(1-(1+ti*qno[ri])/2))
			A[ni,elmts[ei,indir[2]] ] += common*(1+si*qno[qi])*(1+ti*qno[ri])/4
			A[ni,elmts[ei,indir[3]] ] += common*(1-(1+si*qno[qi])/2*(1-(1+ti*qno[ri])/2)- (1+si*qno[qi])*(1+ti*qno[ri])/4)
		end
	end
end

max_stars_repo_name: popsomer/asyBEM

max_issues_repo_name: popsomer/asyBEM

max_forks_repo_name: popsomer/asyBEM


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using DirectTrajectoryOptimization
using LinearAlgebra
using Plots

# ## horizon
T = 51

# ## car
nx = 3
nu = 2
nw = 0

function car(x, u, w)
    [u[1] * cos(x[3]); u[1] * sin(x[3]); u[2]]
end

function midpoint_implicit(y, x, u, w)
    h = 0.1 # timestep
    y - (x + h * car(0.5 * (x + y), u, w))
end

# ## model
dt = Dynamics(midpoint_implicit, nx, nx, nu, nw=nw)
dyn = [dt for t = 1:T-1]
model = DynamicsModel(dyn)

# ## initialization
x1 = [0.0; 0.0; 0.0]
xT = [1.0; 1.0; 0.0]

# ## objective
ot = (x, u, w) -> 0.0 * dot(x - xT, x - xT) + 1.0 * dot(u, u)
oT = (x, u, w) -> 0.0 * dot(x - xT, x - xT)
ct = Cost(ot, nx, nu, nw, [t for t = 1:T-1])
cT = Cost(oT, nx, 0, nw, [T])
obj = [ct, cT]

# ## constraints
ul = -0.5 * ones(nu)
uu = 0.5 * ones(nu)
bnd1 = Bound(nx, nu, [1], xl=x1, xu=x1, ul=ul, uu=uu)
bndt = Bound(nx, nu, [t for t = 2:T-1], ul=ul, uu=uu)
bndT = Bound(nx, 0, [T], xl=xT, xu=xT)

p_obs = [0.5; 0.5]
r_obs = 0.1
function obs(x, u, w)
    e = x[1:2] - p_obs
    return [dot(e, e) - r_obs^2.0]
end
cont = StageConstraint(obs, nx, nu, nw, [t for t = 1:T-1], :inequality)
conT = StageConstraint(obs, nx, 0, nw, [T], :inequality)
cons = ConstraintSet([bnd1, bndt, bndT], [cont, conT])

# ## problem
trajopt = TrajectoryOptimizationProblem(obj, model, cons)
s = Solver(trajopt)

# ## initialize
x_interpolation = linear_interpolation(x1, xT, T)
u_guess = [0.001 * randn(nu) for t = 1:T-1]
z0 = zeros(s.p.num_var)
for (t, idx) in enumerate(s.p.trajopt.model.idx.x)
    z0[idx] = x_interpolation[t]
end
for (t, idx) in enumerate(s.p.trajopt.model.idx.u)
    z0[idx] = u_guess[t]
end
initialize!(s, z0)

# ## solve
@time solve!(s)

# ## solution
@show trajopt.x[1]
@show trajopt.x[T]

# ## state
plot(hcat(trajopt.x...)[1, :], hcat(trajopt.x...)[2, :], label = "", color = :orange, width=2.0)
pts = Plots.partialcircle(0.0, 2.0 * π, 100, r_obs)
cx, cy = Plots.unzip(pts)
plot!(Shape(cx .+ p_obs[1], cy .+ p_obs[2]), color = :black, label = "", linecolor = :black)

# ## control
plot(hcat(trajopt.u[1:end-1]..., trajopt.u[end-1])', linetype = :steppost)


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -2.3024137037198598 <= q <= 2.756640266423038)

                     add_NL_constraint(m, :(gelu(0.8745011162287222 + -0.4747792679508138*gelu(0.20979439694353186 + 0.2523963619940135*$(x[1]) + 0.6039440929888729*$(x[2])) + 0.06822380173369114*gelu(-0.192203947089018 + 0.35618645260012016*$(x[1]) + -0.9380253880003773*$(x[2]))) + gelu(-0.5284568968024108 + 0.8782860084347535*gelu(0.20979439694353186 + 0.2523963619940135*$(x[1]) + 0.6039440929888729*$(x[2])) + 0.990986450920027*gelu(-0.192203947089018 + 0.35618645260012016*$(x[1]) + -0.9380253880003773*$(x[2]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
@testset "runtestExtractInQuoteInside" begin
    println("Testing runtestExtractInQuoteInside")
    include("../src/bfgaExtractInQuoteInside.jl")
    using .bfgaExtractInQuoteInside

    include("../src/Bf.jl")
    using .BfInterpreter

    @time model = test_serial()
    #@time model = test_parallel()
    prems = model.refSet[length(model.refSet)]
    println("Generation : $(model.params.currentGeneration)")
    println("& fitness : $(model.specific_fitness.fitness(prems, model.instructionsSet))")
    #bfcode :: Array{Float64,1}  = prems.dna
    println("code :  $(prems.program)") #$(join(genesToBfInstructions(bfcode), "") ) ")

    #bft = bfType(prems.dna)
    input1 = "One Two \"One Two, repeat\" after me"
    input2 = "The \"RMA had\" been infected"
    n1 , n2 = length(input1), length(input2)

    output1, _ = execute(prems.program, input1, model.instructionsSet)
    output2, _ = execute(prems.program, input2, model.instructionsSet)

    @test length(output1) >= 15
    @test collect("One Two, repeat") == output1[1:15]
    @test length(output2) >= 7
    @test collect("RMA had") == output2[1:7]

    # do the test_parallel after ...
    #test_parallel()
end


text:
# ---
# title: 706. Design HashMap
# id: problem706
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Hash Table, Design
# link: <https://leetcode.com/problems/design-hashmap/description/>
# hidden: true
# ---
#
# Design a HashMap without using any built-in hash table libraries.
#
# To be specific, your design should include these functions:
#
#   * `put(key, value)` : Insert a (key, value) pair into the HashMap. If the value already exists in the HashMap, update the value.
#   * `get(key)`: Returns the value to which the specified key is mapped, or -1 if this map contains no mapping for the key.
#   * `remove(key)` : Remove the mapping for the value key if this map contains the mapping for the key.
#
#
# **Example:**
#
#
#
#     MyHashMap hashMap = new MyHashMap();
#     hashMap.put(1, 1);
#     hashMap.put(2, 2);
#     hashMap.get(1);            // returns 1
#     hashMap.get(3);            // returns -1 (not found)
#     hashMap.put(2, 1);          // update the existing value
#     hashMap.get(2);            // returns 1
#     hashMap.remove(2);          // remove the mapping for 2
#     hashMap.get(2);            // returns -1 (not found)
#
#
#
# **Note:**
#
#   * All keys and values will be in the range of `[0, 1000000]`.
#   * The number of operations will be in the range of `[1, 10000]`.
#   * Please do not use the built-in HashMap library.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using OptimalApplication
using DataFrames
using PrettyTables
using Random
using StatsBase
using Statistics
using Base.Threads
import Printf: @sprintf
using UnicodePlots

const fullscale = false

# A long benchmark; tweak parameters with caution.
# Set fullscale = false to run a smaller benchmark to check formatting etc.

# Number of times to repeat each computation, where min of these is reported as time
const n_reps = fullscale ? 3 : 2

# Number of markets to test at each intersection of the experimental variables
const n_markets = fullscale ? 50 : 5

# Cutoffs to exclude certain large computations
const bnbcutoff = fullscale ? 33 : 10
const twottbnbcutoff = 2^bnbcutoff
const fptascutoff_m = fullscale ? 300 : 200
const fptascutoff_eps = fullscale ? 0.0 : 0.1

# Market sizes
marketsizes_SCM = fullscale ? 4 .^ (2:7) : 2 .^ (3:6)
marketsizes_VCM = fullscale ? 2 .^ (3:8) : 2 .^ (3:6)

function makeunicodeplots(df::DataFrame)
    large_idx = df[!, :m] .== maximum(df[!, :m])
    return [
        histogram(
            Vector{Float64}(df[large_idx, c]),
            ylabel=c,
            title="Time when m = $(marketsizes_SCM[end])",
            canvas=BlockCanvas
        )
        for c in setdiff(names(df), ("m", "time_bnb")) if !(Inf in df[!, c])
    ]
end

function collectmeanstd(df::DataFrame)
    kv_pairs = Pair[]
    for c in names(df)
        if c != "m"
            push!(kv_pairs, c => mean)
            push!(kv_pairs, c => std)
        end
    end
    return combine(groupby(df, :m), kv_pairs...)
end

function benchmark1(marketsizes_SCM = marketsizes_SCM)
    @info "Benchmark 1: Homogeneous-cost algorithms"
    mkts_SCM = SameCostsMarket[SameCostsMarket(m) for m in marketsizes_SCM, i in 1:n_markets]

    sizes = Int[m for m in marketsizes_SCM, _ in 1:n_markets]
    times_list = fill(Inf, length(marketsizes_SCM), n_markets)
    times_heap = fill(Inf, length(marketsizes_SCM), n_markets)

    @threads for j in 1:n_markets
        println("  j = ", @sprintf("%2d", j), " of $n_markets")
        for (i, _) in enumerate(marketsizes_SCM), _ in 1:n_reps
            times_list[i, j] = min(times_list[i, j], @elapsed applicationorder_list(mkts_SCM[i, j]))
            times_heap[i, j] = min(times_heap[i, j], @elapsed applicationorder_heap(mkts_SCM[i, j]))
        end
    end

    df = DataFrame(
        "m" => sizes[:],
        "time_list" => times_list[:],
        "time_heap" => times_heap[:]
    )

    plts = makeunicodeplots(df)
    meanstd = collectmeanstd(df)

    return meanstd, plts, df
end

function benchmark2(marketsizes_VCM = marketsizes_VCM)
    @info "Benchmark 2: Heterogeneous-cost algorithms"
    epsilons = [0.5, 0.05]
    mkts_VCM = VariedCostsMarket[VariedCostsMarket(m) for m in marketsizes_VCM, i in 1:n_markets]

    sizes = Int[m for m in marketsizes_VCM, _ in 1:n_markets]
    times_bnb = fill(Inf, length(marketsizes_VCM), n_markets)
    times_dp = fill(Inf, length(marketsizes_VCM), n_markets)
    times_fptas = fill(Inf, length(epsilons), length(marketsizes_VCM), n_markets)

    @threads for j in 1:n_markets
        println("  j = ", @sprintf("%2d", j), " of $n_markets")
        for (i, m) in enumerate(marketsizes_VCM), _ in 1:n_reps
            if m ≤ bnbcutoff
                times_bnb[i, j] = min(times_bnb[i, j], @elapsed optimalportfolio_branchbound(mkts_VCM[i, j]; maxit=twottbnbcutoff))
                # times_bnb[i, j] =
                #     @belapsed optimalportfolio_branchbound($(mkts_VCM[i, j]); maxit=$twottbnbcutoff)
            end

            times_dp[i, j] = min(times_dp[i, j], @elapsed optimalportfolio_dynamicprogram(mkts_VCM[i, j]))
            # times_dp[i, j] = @belapsed optimalportfolio_dynamicprogram($(mkts_VCM[i, j]))

            for (k, epsilon) in enumerate(epsilons)
                if m ≤ fptascutoff_m || epsilon > fptascutoff_eps
                    times_fptas[k, i, j] = min(times_fptas[k, i, j], @elapsed optimalportfolio_fptas(mkts_VCM[i, j], epsilon))
                    # times_fptas[i, j, k] = @belapsed optimalportfolio_fptas($(mkts_VCM[i, j]), $epsilon)
                end
            end
        end
    end

    df = DataFrame(
        "m" => sizes[:],
        "time_bnb" => times_bnb[:],
        "time_dp" => times_dp[:],
        ["time_fptas_$epsilon" => times_fptas[k, :, :][:] for (k, epsilon) in enumerate(epsilons)]...
    )

    plts = makeunicodeplots(df)
    meanstd = collectmeanstd(df)

    return meanstd, plts, df
end

function fmter(v, i, j)
    j == 1 && return v

    v = ismissing(v) || isinf(v) || isnan(v) ? "—" : @sprintf "%1.2f" 1000 * v

    return isodd(j) ? "($v)" : v
end

function doeverything()
    @info "Benchmark parameters" fullscale n_markets n_reps nthreads()
    println()
    bm1, plts1, = benchmark1()
    display(pretty_table(bm1, formatters=fmter))
    for pl in plts1
        display(pl)
        println()
    end
    bm2, plts2, = benchmark2()
    display(pretty_table(bm2, formatters=fmter))
    for pl in plts2
        display(pl)
        println()
    end
end

@time doeverything()

max_stars_repo_name: maxkapur/OptimalApplication

max_issues_repo_name: maxkapur/OptimalApplication

max_forks_repo_name: maxkapur/OptimalApplication


text:
module JungleHelperSpinyPlant

using ..Ahorn, Maple

@mapdef Entity "JungleHelper/SpinyPlant" SpinyPlant(x::Integer, y::Integer, height::Integer=16, color::String="Blue", sprite::String="")

colors = ["Blue", "Pink", "Yellow", "Orange"]

const placements = Ahorn.PlacementDict(
    "Spiny Plant ($(color)) (Jungle Helper)" => Ahorn.EntityPlacement(
        SpinyPlant,
        "rectangle",
        Dict{String, Any}(
            "color" => color
        )
    ) for color in colors
)

Ahorn.editingOptions(entity::SpinyPlant) = Dict{String, Any}(
    "color" => colors
)

Ahorn.minimumSize(entity::SpinyPlant) = 0, 16
Ahorn.resizable(entity::SpinyPlant) = false, true

function Ahorn.selection(entity::SpinyPlant)
    x, y = Ahorn.position(entity)
    height = Int(get(entity.data, "height", 16))

    return Ahorn.Rectangle(x, y, 24, height)
end

sections = ["Top", "Mid", "Bottom"]

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::SpinyPlant, room::Maple.Room)
    color = get(entity.data, "color", "Blue")

    # Values need to be system specific integer
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    height = Int(get(entity.data, "height", 16))
    left = get(entity.data, "left", true)

    startX = div(x, 8) + 1
    startY = div(y, 8) + 1
    stopY = startY + div(height, 8) - 1

    len = stopY - startY

    if len == 1
        topOpen = get(room.fgTiles.data, (startY - 1, startX + 1), false) == '0'
        bottomOpen = get(room.fgTiles.data, (stopY + 1, startX + 1), false) == '0'

        if topOpen
            if bottomOpen
                section = "Solo"
            else
                section = "Top"
            end
        else
            if bottomOpen
                section = "Bottom"
            else
                section = "Mid"
            end
        end

        Ahorn.drawImage(ctx, "JungleHelper/SpinyPlant/Spiny$(color)$(section)", section == "Solo" ? 4 : 0, 0)
    else
        for i in 0:2:len
            if i == len
                i = len - 1
            end

            section = "Mid"
            if i == 0
                if get(room.fgTiles.data, (startY - 1, startX + 1), false) == '0'
                    section = "Top"
                end
            elseif i == len - 1
                if get(room.fgTiles.data, (stopY + 1, startX + 1), false) == '0'
                    section = "Bottom"
                end
            end

            Ahorn.drawImage(ctx, "JungleHelper/SpinyPlant/Spiny$(color)$(section)", 0, 8 * i)
        end
    end
end

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
using BinDeps

@BinDeps.setup

@linux_only begin
    hdf5 = library_dependency("libhdf5")
    provides(AptGet, "hdf5-tools", hdf5)
    provides(Yum, "hdf5", hdf5)
end

@windows_only begin
    hdf5 = library_dependency("hdf5")

    const OS_ARCH = WORD_SIZE == 64 ? "x86_64" : "x86"

    if WORD_SIZE == 32
            provides(Binaries,URI("https://ia601003.us.archive.org/29/items/julialang/windows/hdf5-1.8-win32.7z"),hdf5,os = :Windows)
    else
            provides(Binaries,URI("https://ia601003.us.archive.org/29/items/julialang/windows/hdf5-1.8-win64.7z"),hdf5,os = :Windows)
    end
    push!(DL_LOAD_PATH, joinpath(Pkg.dir("HDF5/deps/usr/lib/"), OS_ARCH))
end

@osx_only begin
    if Pkg.installed("Homebrew") === nothing
            error("Homebrew package not installed, please run Pkg.add(\"Homebrew\")")
    end
    using Homebrew
    hdf5 = library_dependency("libhdf5")
    provides( Homebrew.HB, "hdf5", hdf5, os = :Darwin )
end

@BinDeps.install


text:
# ------------------------------------------
#    Data to fit

export data_params

"""
data_params(; kwargs...)` : returns a dictionary containing the
complete set of parameters required to read data from an .xyz file. j
All parameters are passed as keyword arguments.

### Parameters
* `fname` : a "*.xyz" file with atomistic data (mandatory).
* `energy_key = "dft_energy` : ASE's `Atoms.info` key to read energy
for fitting.
* `force_key = "dft_force` : ASE's `Atoms.arrays` key to read forces
for fitting.
* `virial_key = "dft_virial` : ASE's `Atoms.info` key to read virial
for fitting.
"""
function data_params(;
    fname = nothing,
    energy_key = "dft_energy",
    force_key = "dft_force",
    virial_key = "dft_virial"
)
    @assert !isnothing(fname) "`fname` must be given. "

    return Dict(
        "fname" => fname,
        "energy_key" => energy_key,
        "force_key" => force_key,
        "virial_key" => virial_key
    )

end

function read_data(params::Dict)
    return IPFitting.Data.read_xyz(
        params["fname"];
        energy_key = params["energy_key"],
        force_key = params["force_key"],
        virial_key = params["virial_key"])
end


text:
module Shoco

let depsfile = joinpath(@__DIR__, "..", "deps", "deps.jl")
    if !isfile(depsfile)
        error("Shoco is not properly installed. Please run Pkg.build(\"Shoco\")",
              " and restart Julia.")
    end
    include(depsfile)
end

__init__() = check_deps()

export compress, decompress

function compress(s::AbstractString)
    isempty(s) && return ""
    # The output should be no longer than the input
    compressed = Vector{UInt8}(undef, sizeof(s))
    # The function modifies `compressed` and returns the number of bytes written
    nbytes = ccall((:shoco_compress, libshoco), Int,
                   (Ptr{Cchar}, Csize_t, Ptr{UInt8}, Csize_t),
                   s, 0, compressed, sizeof(s))
    nbytes > 0 || error("Compression failed for input $s")
    resize!(compressed, nbytes)
    return String(compressed)
end

function decompress(s::AbstractString)
    isempty(s) && return ""
    # The decompressed string will be at most twice as long as the input
    decompressed = Vector{UInt8}(undef, 2 * sizeof(s))
    nbytes = ccall((:shoco_decompress, libshoco), Int,
                   (Ptr{Cchar}, Csize_t, Ptr{UInt8}, Csize_t),
                   s, sizeof(s), decompressed, 2 * sizeof(s))
    nbytes > 0 || error("Decompression failed for input $s")
    resize!(decompressed, nbytes)
    return String(decompressed)
end

end # module


text:
using CpuId
using Logging
using Pkg
using Printf

using StructuralIdentifiability
using StructuralIdentifiability: _runtime_logger, ODE

logger = Logging.SimpleLogger(stdout, Logging.Info)
global_logger(logger)

include("benchmarks.jl")

runtimes = Dict()
TIME_CATEGORIES = [:loc_time, :glob_time, :ioeq_time, :wrnsk_time, :rank_time, :check_time, :total]
NUM_RUNS = 5

for bmark in benchmarks
    name = bmark[:name]
    if (:skip in keys(bmark)) && bmark[:skip]
        @info "Skipping $name"
        continue
    end
    @info "Processing $name"
    runtimes[bmark[:name]] = Dict(c => 0. for c in TIME_CATEGORIES)
    for _ in 1:NUM_RUNS
        runtimes[name][:total] += @elapsed assess_identifiability(bmark[:ode])
        for cat in TIME_CATEGORIES[1:end - 1]
            runtimes[name][cat] += _runtime_logger[cat]
        end
    end
    for k in keys(runtimes[name])
        runtimes[name][k] = runtimes[name][k] / NUM_RUNS
    end
end

resulting_md = ""

resulting_md *= "|Model|" * join(TIME_CATEGORIES, "|") * "|\n"
resulting_md *= "|-----|" * join(["---" for _ in TIME_CATEGORIES], "|") * "|\n"
for (name, times) in runtimes
    global resulting_md *= "|$name|"
    for c in TIME_CATEGORIES
        resulting_md *= @sprintf("%.2f", times[c]) * "|"
    end
    resulting_md *= "\n"
end

resulting_md *= "\n*Benchmarking environment:*\n\n"
resulting_md *= "* Total RAM (Mb): $(Sys.total_memory() / 2^20)\n"
resulting_md *= "* Processor: $(cpubrand())\n"
resulting_md *= "* Julia version: $(VERSION)\n\n"
resulting_md *= "Versions of the dependencies:\n\n"

deps = Pkg.dependencies()
stid_info = deps[findfirst(x -> x.name == "StructuralIdentifiability", deps)]
for (s, uid) in stid_info.dependencies
    if deps[uid].version != nothing
        global resulting_md *= "* $s : $(deps[uid].version)\n"
    end
end

open("benchmark_result.md", "w") do io
    write(io, resulting_md)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
The `AbstractChar` type is the supertype of all character implementations
in Julia. A character represents a Unicode code point, and can be converted
to an integer via the [`codepoint`](@ref) function in order to obtain the
numerical value of the code point, or constructed from the same integer.
These numerical values determine how characters are compared with `<` and `==`,
for example.  New `T <: AbstractChar` types should define a `codepoint(::T)`
method and a `T(::UInt32)` constructor, at minimum.

A given `AbstractChar` subtype may be capable of representing only a subset
of Unicode, in which case conversion from an unsupported `UInt32` value
may throw an error. Conversely, the built-in [`Char`](@ref) type represents
a *superset* of Unicode (in order to losslessly encode invalid byte streams),
in which case conversion of a non-Unicode value *to* `UInt32` throws an error.
The [`isvalid`](@ref) function can be used to check which codepoints are
representable in a given `AbstractChar` type.

Internally, an `AbstractChar` type may use a variety of encodings.  Conversion
via `codepoint(char)` will not reveal this encoding because it always returns the
Unicode value of the character. `print(io, c)` of any `c::AbstractChar`
produces an encoding determined by `io` (UTF-8 for all built-in `IO`
types), via conversion to `Char` if necessary.

`write(io, c)`, in contrast, may emit an encoding depending on
`typeof(c)`, and `read(io, typeof(c))` should read the same encoding as `write`.
New `AbstractChar` types must provide their own implementations of
`write` and `read`.
"""
AbstractChar

"""
    Char(c::Union{Number,AbstractChar})

`Char` is a 32-bit [`AbstractChar`](@ref) type that is the default representation
of characters in Julia. `Char` is the type used for character literals like `'x'`
and it is also the element type of [`String`](@ref).

In order to losslessly represent arbitrary byte streams stored in a `String`,
a `Char` value may store information that cannot be converted to a Unicode
codepoint — converting such a `Char` to `UInt32` will throw an error.
The [`isvalid(c::Char)`](@ref) function can be used to query whether `c`
represents a valid Unicode character.
"""
Char

(::Type{T})(x::Number) where {T<:AbstractChar} = T(UInt32(x))
(::Type{AbstractChar})(x::Number) = Char(x)
(::Type{T})(x::AbstractChar) where {T<:Union{Number,AbstractChar}} = T(codepoint(x))
(::Type{T})(x::T) where {T<:AbstractChar} = x

"""
    ncodeunits(c::Char) -> Int

Return the number of code units required to encode a character as UTF-8.
This is the number of bytes which will be printed if the character is written
to an output stream, or `ncodeunits(string(c))` but computed efficiently.

!!! compat "Julia 1.1"
    This method requires at least Julia 1.1. In Julia 1.0 consider
    using `ncodeunits(string(c))`.
"""
ncodeunits(c::Char) = write(devnull, c) # this is surprisingly efficient

"""
    codepoint(c::AbstractChar) -> Integer

Return the Unicode codepoint (an unsigned integer) corresponding
to the character `c` (or throw an exception if `c` does not represent
a valid character). For `Char`, this is a `UInt32` value, but
`AbstractChar` types that represent only a subset of Unicode may
return a different-sized integer (e.g. `UInt8`).
"""
function codepoint end

codepoint(c::Char) = UInt32(c)

struct InvalidCharError{T<:AbstractChar} <: Exception
    char::T
end
struct CodePointError{T<:Integer} <: Exception
    code::T
end
@noinline invalid_char(c::AbstractChar) = throw(InvalidCharError(c))
@noinline code_point_err(u::Integer) = throw(CodePointError(u))

function ismalformed(c::Char)
    u = reinterpret(UInt32, c)
    l1 = leading_ones(u) << 3
    t0 = trailing_zeros(u) & 56
    (l1 == 8) | (l1 + t0 > 32) |
    (((u & 0x00c0c0c0) ⊻ 0x00808080) >> t0 != 0)
end

@inline is_overlong_enc(u::UInt32) = (u >> 24 == 0xc0) | (u >> 24 == 0xc1) | (u >> 21 == 0x0704) | (u >> 20 == 0x0f08)

function isoverlong(c::Char)
    u = reinterpret(UInt32, c)
    is_overlong_enc(u)
end

# fallback: other AbstractChar types, by default, are assumed
#           not to support malformed or overlong encodings.

"""
    ismalformed(c::AbstractChar) -> Bool

Return `true` if `c` represents malformed (non-Unicode) data according to the
encoding used by `c`.  Defaults to `false` for non-`Char` types.  See also
[`show_invalid`](@ref).
"""
ismalformed(c::AbstractChar) = false

"""
    isoverlong(c::AbstractChar) -> Bool

Return `true` if `c` represents an overlong UTF-8 sequence. Defaults
to `false` for non-`Char` types.  See also [`decode_overlong`](@ref)
and [`show_invalid`](@ref).
"""
isoverlong(c::AbstractChar) = false

function UInt32(c::Char)
    # TODO: use optimized inline LLVM
    u = reinterpret(UInt32, c)
    u < 0x80000000 && return u >> 24
    l1 = leading_ones(u)
    t0 = trailing_zeros(u) & 56
    (l1 == 1) | (8l1 + t0 > 32) |
    ((((u & 0x00c0c0c0) ⊻ 0x00808080) >> t0 != 0) | is_overlong_enc(u)) &&
        invalid_char(c)::Union{}
    u &= 0xffffffff >> l1
    u >>= t0
    ((u & 0x0000007f) >> 0) | ((u & 0x00007f00) >> 2) |
    ((u & 0x007f0000) >> 4) | ((u & 0x7f000000) >> 6)
end

"""
    decode_overlong(c::AbstractChar) -> Integer

When [`isoverlong(c)`](@ref) is `true`, `decode_overlong(c)` returns
the Unicode codepoint value of `c`. `AbstractChar` implementations
that support overlong encodings should implement `Base.decode_overlong`.
"""
function decode_overlong end

function decode_overlong(c::Char)
    u = reinterpret(UInt32, c)
    l1 = leading_ones(u)
    t0 = trailing_zeros(u) & 56
    u &= 0xffffffff >> l1
    u >>= t0
    ((u & 0x0000007f) >> 0) | ((u & 0x00007f00) >> 2) |
    ((u & 0x007f0000) >> 4) | ((u & 0x7f000000) >> 6)
end

function Char(u::UInt32)
    u < 0x80 && return reinterpret(Char, u << 24)
    u < 0x00200000 || code_point_err(u)::Union{}
    c = ((u << 0) & 0x0000003f) | ((u << 2) & 0x00003f00) |
        ((u << 4) & 0x003f0000) | ((u << 6) & 0x3f000000)
    c = u < 0x00000800 ? (c << 16) | 0xc0800000 :
        u < 0x00010000 ? (c << 08) | 0xe0808000 :
                         (c << 00) | 0xf0808080
    reinterpret(Char, c)
end

function (T::Union{Type{Int8},Type{UInt8}})(c::Char)
    i = reinterpret(Int32, c)
    i ≥ 0 ? ((i >>> 24) % T) : T(UInt32(c))
end

function Char(b::Union{Int8,UInt8})
    0 ≤ b ≤ 0x7f ? reinterpret(Char, (b % UInt32) << 24) : Char(UInt32(b))
end

convert(::Type{AbstractChar}, x::Number) = Char(x) # default to Char
convert(::Type{T}, x::Number) where {T<:AbstractChar} = T(x)
convert(::Type{T}, x::AbstractChar) where {T<:Number} = T(x)
convert(::Type{T}, c::AbstractChar) where {T<:AbstractChar} = T(c)
convert(::Type{T}, c::T) where {T<:AbstractChar} = c

rem(x::AbstractChar, ::Type{T}) where {T<:Number} = rem(codepoint(x), T)

typemax(::Type{Char}) = reinterpret(Char, typemax(UInt32))
typemin(::Type{Char}) = reinterpret(Char, typemin(UInt32))

size(c::AbstractChar) = ()
size(c::AbstractChar, d::Integer) = d < 1 ? throw(BoundsError()) : 1
ndims(c::AbstractChar) = 0
ndims(::Type{<:AbstractChar}) = 0
length(c::AbstractChar) = 1
IteratorSize(::Type{Char}) = HasShape{0}()
firstindex(c::AbstractChar) = 1
lastindex(c::AbstractChar) = 1
getindex(c::AbstractChar) = c
getindex(c::AbstractChar, i::Integer) = i == 1 ? c : throw(BoundsError())
getindex(c::AbstractChar, I::Integer...) = all(x -> x == 1, I) ? c : throw(BoundsError())
first(c::AbstractChar) = c
last(c::AbstractChar) = c
eltype(::Type{T}) where {T<:AbstractChar} = T

iterate(c::AbstractChar, done=false) = done ? nothing : (c, true)
isempty(c::AbstractChar) = false
in(x::AbstractChar, y::AbstractChar) = x == y

==(x::Char, y::Char) = reinterpret(UInt32, x) == reinterpret(UInt32, y)
isless(x::Char, y::Char) = reinterpret(UInt32, x) < reinterpret(UInt32, y)
hash(x::Char, h::UInt) =
    hash_uint64(((reinterpret(UInt32, x) + UInt64(0xd4d64234)) << 32) ⊻ UInt64(h))

first_utf8_byte(c::Char) = (reinterpret(UInt32, c) >> 24) % UInt8

# fallbacks:
isless(x::AbstractChar, y::AbstractChar) = isless(Char(x), Char(y))
==(x::AbstractChar, y::AbstractChar) = Char(x) == Char(y)
hash(x::AbstractChar, h::UInt) = hash(Char(x), h)
widen(::Type{T}) where {T<:AbstractChar} = T

@inline -(x::AbstractChar, y::AbstractChar) = Int(x) - Int(y)
@inline -(x::T, y::Integer) where {T<:AbstractChar} = T(Int32(x) - Int32(y))
@inline +(x::T, y::Integer) where {T<:AbstractChar} = T(Int32(x) + Int32(y))
@inline +(x::Integer, y::AbstractChar) = y + x

# `print` should output UTF-8 by default for all AbstractChar types.
# (Packages may implement other IO subtypes to specify different encodings.)
# In contrast, `write(io, c)` outputs a `c` in an encoding determined by typeof(c).
print(io::IO, c::Char) = (write(io, c); nothing)
print(io::IO, c::AbstractChar) = print(io, Char(c)) # fallback: convert to output UTF-8

const hex_chars = UInt8['0', '1', '2', '3', '4', '5', '6', '7', '8', '9',
                        'a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i',
                        'j', 'k', 'l', 'm', 'n', 'o', 'p', 'q', 'r',
                        's', 't', 'u', 'v', 'w', 'x', 'y', 'z']

function show_invalid(io::IO, c::Char)
    write(io, 0x27)
    u = reinterpret(UInt32, c)
    while true
        a = hex_chars[((u >> 28) & 0xf) + 1]
        b = hex_chars[((u >> 24) & 0xf) + 1]
        write(io, 0x5c, UInt8('x'), a, b)
        (u <<= 8) == 0 && break
    end
    write(io, 0x27)
end

"""
    show_invalid(io::IO, c::AbstractChar)

Called by `show(io, c)` when [`isoverlong(c)`](@ref) or
[`ismalformed(c)`](@ref) return `true`.   Subclasses
of `AbstractChar` should define `Base.show_invalid` methods
if they support storing invalid character data.
"""
show_invalid

# show c to io, assuming UTF-8 encoded output
function show(io::IO, c::AbstractChar)
    if c <= '\\'
        b = c == '\0' ? 0x30 :
            c == '\a' ? 0x61 :
            c == '\b' ? 0x62 :
            c == '\t' ? 0x74 :
            c == '\n' ? 0x6e :
            c == '\v' ? 0x76 :
            c == '\f' ? 0x66 :
            c == '\r' ? 0x72 :
            c == '\e' ? 0x65 :
            c == '\'' ? 0x27 :
            c == '\\' ? 0x5c : 0xff
        if b != 0xff
            write(io, 0x27, 0x5c, b, 0x27)
            return
        end
    end
    if isoverlong(c) || ismalformed(c)
        show_invalid(io, c)
    elseif isprint(c)
        write(io, 0x27)
        print(io, c) # use print, not write, to use UTF-8 for any AbstractChar
        write(io, 0x27)
    else # unprintable, well-formed, non-overlong Unicode
        u = codepoint(c)
        write(io, 0x27, 0x5c, u <= 0x7f ? 0x78 : u <= 0xffff ? 0x75 : 0x55)
        d = max(2, 8 - (leading_zeros(u) >> 2))
        while 0 < d
            write(io, hex_chars[((u >> ((d -= 1) << 2)) & 0xf) + 1])
        end
        write(io, 0x27)
    end
    return
end

function show(io::IO, ::MIME"text/plain", c::T) where {T<:AbstractChar}
    show(io, c)
    get(io, :compact, false) && return
    if !ismalformed(c)
        print(io, ": ")
        if isoverlong(c)
            print(io, "[overlong] ")
            u = decode_overlong(c)
            c = T(u)
        else
            u = codepoint(c)
        end
        h = uppercase(string(u, base = 16, pad = 4))
        print(io, (isascii(c) ? "ASCII/" : ""), "Unicode U+", h)
    else
        print(io, ": Malformed UTF-8")
    end
    abr = Unicode.category_abbrev(c)
    str = Unicode.category_string(c)
    print(io, " (category ", abr, ": ", str, ")")
end

max_stars_repo_name: naveenarun/julia

max_issues_repo_name: naveenarun/julia

max_forks_repo_name: naveenarun/julia


text:
ComplexGradientConfig(f::F,v::Vector{<:Element},w...) where {F} =
      ComplexGradientConfig(f,position(v),w...)

"""
    seed_zeros(v::Vector{Element},cfg::ComplexGradientConfig)

Given a collection `v` of points or blobs, create a copy of the collection with
the element's positions and strengths replaced by complex duals with zeros.
`cfg` specifies the type and size of duals. The use of complex duals (in spite of the real-valued
strength) ensures correct treatment for gradient calculation, since some
calculations will be complex-valued.
"""
function seed_zeros(strength::F,v::Vector{<:Element},cfg::ComplexGradientConfig) where F
  posduals = convert.(eltype(cfg.duals),position(v))
  strduals = convert.(eltype(cfg.duals),complex(strength.(v)))
  return posduals, real(strduals)
end

function seed!(posduals::AbstractArray{<:ComplexDual{T,V,M}},strduals::AbstractArray{<:ComplexDual{T,V,M}},strfunc::F,v::Vector{<:Element}) where {F,T,V,M}
  seed!(posduals,position(v))
  seed!(strduals,complex(strfunc.(v)))
  return posduals, real(strduals)
end

function seed!(posduals::AbstractArray{<:ComplexDual{T,V,M}},strduals::AbstractArray{<:ComplexDual{T,V,M}},strfunc::F,v::Vector{<:Element},index) where {F,T,V,M}
  seed!(posduals,position(v),index)
  seed!(strduals,complex(strfunc.(v)),index)
  return posduals, real(strduals)
end

"""
    seed_position(v::Vector{Element},cfg::ComplexGradientConfig)

Given a collection `v` of points or blobs, create a copy of the collection with
the element's positions replaced by complex duals with unit seeds and with
the element's strengths replaced by complex duals with zeros. `cfg` specifies
the type and size of duals. The length of the partials in these duals is twice
the number of elements, to enable computation of gradients with respect to
all positions.
"""
function seed_position!(posduals::AbstractArray{<:ComplexDual{T,V,M}},strduals::AbstractArray{<:ComplexDual{T,V,M}},
                        strfunc::F,v::Vector{<:Element},rseeds::NTuple{N,Partials{M,V}},iseeds::NTuple{N,Partials{M,V}}) where {T,F,N,M,V}
  seed!(posduals,position(v),rseeds,iseeds)
  seed!(strduals,complex(strfunc.(v)))
  return posduals, real(strduals)
end

function seed_position!(posduals::AbstractArray{<:ComplexDual{T,V,M}},strduals::AbstractArray{<:ComplexDual{T,V,M}},
                        strfunc::F,v::Vector{<:Element},index,rseeds::NTuple{N,Partials{M,V}},iseeds::NTuple{N,Partials{M,V}},
                        chunksize=N) where {T,F,N,M,V}
  seed!(posduals,position(v),index,rseeds,iseeds,chunksize)
  seed!(strduals,complex(strfunc.(v)),index,zero(Partials{M,V}),zero(Partials{M,V}),chunksize)
  return posduals, real(strduals)
end

seed_position(strfunc::F,v::Vector{<:Element},cfg::ComplexGradientConfig) where F =
    seed_position!(copy(cfg.duals),copy(cfg.duals),strfunc,v,cfg.rseeds,cfg.iseeds)

"""
    seed_strength(v::Vector{Element},cfg::ComplexGradientConfig)

Given a collection `v` of points or blobs, create a copy of the collection with
the element's strengths replaced by complex duals with unit seeds and with
the element's positions replaced by complex duals with zeros. `cfg` specifies
the type and size of duals. The use of complex duals (in spite of the real-valued
strength) ensures correct treatment for gradient calculation, since some
calculations will be complex-valued.
"""
function seed_strength!(posduals::AbstractArray{<:ComplexDual{T,V,M}},strduals::AbstractArray{<:ComplexDual{T,V,M}},
                        strfunc::F,v::Vector{<:Element},rseeds::NTuple{N,Partials{M,V}},iseeds::NTuple{N,Partials{M,V}}) where {T,F,N,M,V}
  seed!(posduals,position(v))
  seed!(strduals,complex(strfunc.(v)),rseeds,iseeds)
  return posduals, real(strduals)
end

function seed_strength!(posduals::AbstractArray{<:ComplexDual{T,V,M}},strduals::AbstractArray{<:ComplexDual{T,V,M}},
                        strfunc::F,v::Vector{<:Element},index,
                        rseeds::NTuple{N,Partials{M,V}},iseeds::NTuple{N,Partials{M,V}},chunksize=N) where {T,F,N,M,V}
  seed!(posduals,position(v),index,zero(Partials{M,V}),zero(Partials{M,V}),chunksize)
  seed!(strduals,complex(strfunc.(v)),index,rseeds,iseeds,chunksize)
  return posduals, real(strduals)
end

seed_strength(strfunc::F,v::Vector{<:Element},cfg::ComplexGradientConfig) where F =
    seed_strength!(copy(cfg.duals),copy(cfg.duals),strfunc,v,cfg.rseeds,cfg.iseeds)

## Gradient

function gradient_position(f,v::Vector{<:Element},cfg::ComplexGradientConfig{T} = ComplexGradientConfig(f, v)) where {T}
    checktag(T, f, property_type(eltype(v)))
    if chunksize(cfg) == length(v)
      vector_mode_gradient(f,v,cfg,vector_mode_dual_eval_position)
    else
      chunk_mode_gradient_position(f,v,cfg)
    end
end

function gradient_strength(f,v::Vector{<:Element},cfg::ComplexGradientConfig{T} = ComplexGradientConfig(f, v)) where {T}
    checktag(T, f, property_type(eltype(v)))
    if chunksize(cfg) == length(v)
      dz, dzstar = vector_mode_gradient(f,v,cfg,vector_mode_dual_eval_strength)
    else
      dz, dzstar = chunk_mode_gradient_strength(f,v,cfg)
    end
    return dz .+ dzstar
end

@inline vector_mode_dual_eval_position(f::F, v::Vector{<:Element},
          cfg::ComplexGradientConfig) where {F} = f(seed_position(v,cfg))

@inline vector_mode_dual_eval_strength(f::F, v::Vector{<:Element},
          cfg::ComplexGradientConfig) where {F} = f(seed_strength(v,cfg))

@inline vector_mode_dual_eval_param(f::F, v::Tuple{Vector{<:Element},T},
          cfg::ComplexGradientConfig) where {F,T} = f(seed_zeros(v[1],cfg),seed(v[2],cfg))

function jacobian_position(f,v::Vector{<:Element},cfg::ComplexGradientConfig{T} = ComplexGradientConfig(f, v)) where {T}
    checktag(T, f, property_type(eltype(v)))
    if chunksize(cfg) == length(v)
      vector_mode_jacobian(f,v,cfg,vector_mode_dual_eval_position)
    else
      chunk_mode_jacobian_position(f,v,cfg)
    end
end

function jacobian_strength(f,v::Vector{<:Element},cfg::ComplexGradientConfig{T} = ComplexGradientConfig(f, v)) where {T}
    checktag(T, f, property_type(eltype(v)))
    if chunksize(cfg) == length(v)
      dz, dzstar = vector_mode_jacobian(f,v,cfg,vector_mode_dual_eval_strength)
    else
      dz, dzstar = chunk_mode_jacobian_strength(f,v,cfg)
    end
    return dz .+ dzstar
end

function jacobian_param(f,v::Tuple{Vector{<:Element},S},
        cfg::ComplexGradientConfig{T} = ComplexGradientConfig(f, ComplexF64[v[2]])) where {S,T}
    checktag(T, f, v[2])
    dz, dzstar = vector_mode_jacobian(f,v,cfg,vector_mode_dual_eval_param)
    return real(dz .+ dzstar)
end

for ftype = (:position, :strength)
  fname = Symbol("chunk_mode_gradient_",ftype)
  sname = Symbol("seed_",ftype,"!")
  @eval function $fname(f::F, z::Vector{<:Element}, cfg::ComplexGradientConfig{T,V,N}) where {F,T,V,N}
     $(chunk_mode_gradient_expr(quote
                                  posduals = copy(cfg.duals)
                                  strduals = copy(cfg.duals)
                                  zdual = seed_zeros(z,cfg)
                                  seed!(zdual,z,posduals,strduals)
                                end,
                                quote
                                  dz = similar(z, Complex{valtype(ydual)})
                                  dzstar = similar(z, Complex{valtype(ydual)})
                                end,
                                :(dz, dzstar),
                                :(ydual = f(zdual)),
                                :(),
                                :($sname(zdual, z, posduals, strduals, index, rseeds, iseeds, chunksize)),
                                :(seed!(zdual, z, posduals, strduals, index)),
                                :(extract_gradient_chunk!(T, dz, dzstar, ydual, index, chunksize)),
                                :()))
  end
end

for ftype = (:position, :strength)
  fname = Symbol("chunk_mode_jacobian_",ftype)
  sname = Symbol("seed_",ftype,"!")
  @eval function $fname(f::F, z::Vector{<:Element}, cfg::ComplexGradientConfig{T,V,N}) where {F,T,V,N}
     $(chunk_mode_gradient_expr(quote
                                  posduals = copy(cfg.duals)
                                  strduals = copy(cfg.duals)
                                  zdual = seed_zeros(z,cfg)
                                  seed!(zdual,z,posduals,strduals)
                                end,
                                quote
                                  dz = similar(z, Complex{valtype(ydual)},length(ydual), zlen)
                                  dzstar = similar(z, Complex{valtype(ydual)},length(ydual), zlen)
                                  dz_reshaped = reshape_jacobian(dz, ydual, zdual)
                                  dzstar_reshaped = reshape_jacobian(dzstar, ydual, zdual)
                                end,
                                :(dz, dzstar),
                                :(ydual = f(zdual)),
                                :(ydual isa AbstractArray || throw(JACOBIAN_ERROR)),
                                :($sname(zdual, z, posduals, strduals, index, rseeds, iseeds, chunksize)),
                                :(seed!(zdual, z, posduals, strduals, index)),
                                :(extract_jacobian_chunk!(T, dz_reshaped, dzstar_reshaped, ydual, index, chunksize)),
                                :()))
  end
end


text:
module Lenny

#--- Fundamentals

abstract type AbstractContinuationProblem{T} end
abstract type AbstractToolbox{T} end
abstract type AbstractCovering{T} end
abstract type AbstractSolver{T} end
abstract type AbstractZeroProblem{T} end

function initialise end
function finalise end
function toolboxes end
function covering end
function zeroproblem end
function solver end

#--- Modules

include("FiniteStateMachines.jl")
include("Coverings.jl")
include("ZeroProblems.jl")
include("ODEs.jl")

#--- Basic continuation problem structure (user facing)

mutable struct ContinuationProblem{T} <: AbstractContinuationProblem{T}
	options::Dict{String, Any}
	toolboxes::Vector{AbstractToolbox{T}}
	covering::AbstractCovering{T}
	efuncs
	solver::AbstractSolver{T}
end

ContinuationProblem{T}() where T = ContinuationProblem(Dict{String, Any}(), AbstractToolbox{T}[], ZeroProblems.emptyCoverings.DefaultCovering())

Base.push!(prob::ContinuationProblem{T}, tbx::AbstractToolbox{T}) where T = (push!(prob.toolboxes, tbx); prob)
Base.push!(prob::ContinuationProblem{T}, cov::AbstractCovering{T}) where T = (prob.covering = cov; prob)
Base.setindex!(prob::ContinuationProblem, idx, val) = setindex!(prob.options, idx, val)

function initialise(prob::ContinuationProblem{T}) where T
	# Initialise toolboxes
	for i in eachindex(prob.toolboxes)
		prob.toolboxes[i] = initialise(prob.toolboxes[i], prob)
	end
	# Initialise solver
	# prob.solver = initialise(prob.solver, prob)
	# Initialise covering
	prob.covering = initialise(prob.covering, prob)
	# Initialise efuncs
	prob.efuncs = initialise(prob.efuncs, prob)
	# Return the final problem structure (TODO: specialise at this point)
	return prob
end

toolboxes(prob::ContinuationProblem) = prob.toolboxes
covering(prob::ContinuationProblem) = prob.covering
efuncs(prob::ContinuationProblem) = prob.efuncs
solver(prob::ContinuationProblem) = prob.solver

end # module


text:
macro c(ret_type, func, arg_types, lib)
  local args_in = Any[ symbol(string('a',x)) for x in 1:length(arg_types.args) ]
  quote
    $(esc(func))($(args_in...)) = ccall( ($(string(func)), $(Expr(:quote, lib)) ), $ret_type, $arg_types, $(args_in...) )
  end
end

macro ctypedef(fake_t,real_t)
  quote
    typealias $fake_t $real_t
  end
end

@ctypedef MD5_CTX Void
@ctypedef HMAC_CTX Void
@ctypedef SHA_CTX Void
@ctypedef SHA256_CTX Void
@ctypedef SHA512_CTX Void
@ctypedef evp_sign_method Void
@ctypedef evp_verify_method Void
@ctypedef EVP_CIPHER_INFO Void
@ctypedef EVP_ENCODE_CTX Void
@ctypedef EVP_PBE_KEYGEN Void
@ctypedef EVP_PKEY_gen_cb Void


text:
"""
    default_obsdim(data)

Specify the default observation dimension for `data`.
Falls back to `nothing` when an observation dimension is undefined.

By default, the following implementations are provided:
```julia
default_obsdim(x::nothing) = nothing
default_obsdim(x::AbstractArray) = ndims(x)
````
"""
default_obsdim(data) = nothing
default_obsdim(A::AbstractArray{T,N}) where {T,N} = N

"""
   getobs(data, idx; obsdim = default_obsdim(data))

Return the observations corresponding to the observation-index `idx`.
Note that `idx` can be of type `Int` or `AbstractVector`.
*Both options must be supported by a custom type.*

The returned observation(s) should be in the form intended to
be passed as-is to some learning algorithm. There is no strict
interface requirement on how this "actual data" must look like.
Every author behind some custom data container can make this
decision himself/herself. We do, however, expect it to be consistent
for `idx` being an integer, as well as `idx` being a vector, respectively.

If it makes sense for the type of `data`, `obsdim` can be used
to indicate which dimension of `data` denotes the observations.
See [`default_obsdim`](@ref) for defining a default dimension.

Data types implementing `LearnBase.getobs` typically also
implement `StatsBase.nobs`.

# Examples

Let's see how to implement a dataset interface for a dataset
represented by an array:
```julia
using LearnBase

function LearnBase.getobs(x::AbstractArray{T,N}, idx; obsdim=default_obsdim(x)) where {T,N}
    _idx = ntuple(i->  i == obsdim ? idx : Colon(), N)
    return x[_idx...]
end

# LearnBase imports nobs from StatsBase
LearnBase.nobs(x::AbstractArray; obsdim=default_obsdim(x)) = size(x, obsdim)

X = rand(2,3)

nobs(X) # == 3

getobs(X, 2) # same as X[:,2]
```
In a supervised learning setting, it can be convenient to
interpret a tuple of arrays as the inputs and the tagets:
```julia
# Here we use Ref to protect idx against broadcasting
LearnBase.getobs(t::Tuple, idx) = getobs.(t, Ref(idx))

# Assume all elements have the same nummber of observations.
# It would be safer to check explicitely though.
LearnBase.nobs(t::Tuple) = nobs(t[1])

# A dataset with 3 observations, each with 2 input features
X, Y = rand(2, 3), rand(3)
dataset = (X, Y)

getobs(dataset, 2) # -> (X[:,2], Y[2])
getobs(dataset, 1:2) # -> (X[:,1:2], Y[1:2])
```
"""
function getobs end

function getobs(data::AbstractArray{T,N}, idx; obsdim::Union{Int,Nothing}=nothing) where {T, N}
   od = obsdim === nothing ? default_obsdim(data) : obsdim
   _idx = ntuple(i -> i == od ? idx : Colon(), N)
   data[_idx...]
end

function getobs(data::Union{Tuple, NamedTuple}, i; obsdim::Union{Int,Nothing}=default_obsdim(data))
   # We don't force users to handle the obsdim keyword if not necessary.
   fobs = obsdim === nothing ? Base.Fix2(getobs, i) : x -> getobs(x, i; obsdim=obsdim)
   map(fobs, data)
end

function getobs(data::D, i; obsdim::Union{Int,Nothing}=default_obsdim(data)) where {D<:AbstractDict}
   fobs = obsdim === nothing ? Base.Fix2(getobs, i) : x -> getobs(x, i; obsdim=obsdim)
   # Cannot return D because the value type can change
   Dict(k => fobs(v) for (k, v) in pairs(data))
end

"""
   getobs!(buffer, data, idx; obsdim = default_obsdim(obsdim))

Inplace version of `getobs(data, idx; obsdim)`. If this method
is defined for the type of `data`, then `buffer` should be used
to store the result, instead of allocating a dedicated object.

Implementing this function is optional. In the case no such
method is provided for the type of `data`, then `buffer` will be
*ignored* and the result of `getobs` returned. This could be
because the type of `data` may not lend itself to the concept
of `copy!`. Thus, supporting a custom `getobs!(::MyType, ...)`
is optional and not required.

If it makes sense for the type of `data`, `obsdim` can be used
to indicate which dimension of `data` denotes the observations.
See [`default_obsdim`](@ref) for defining a default dimension.
"""
function getobs! end

"""
   gettargets(data, idx; obsdim = default_obsdim(data))

Return the targets corresponding to the observation-index `idx`.
Note that `idx` can be of type `Int` or `AbstractVector`.
Implementing this function for a custom type of `data` is
optional. It is particularly useful if the targets in `data` can
be provided without invoking [`getobs`](@ref). For example if you have a
remote data-source where the labels are part of some metadata
that is locally available.
If it makes sense for the type of `data`, `obsdim` can be used
to indicate which dimension of `data` denotes the observations.
See [`default_obsdim`](@ref) for defining a default dimension.
"""
function gettargets end

"""
     datasubset(data, [idx]; obsdim = default_obsdim(data))

Return a lazy subset of the observations in `data` that correspond
to the given `idx`. No data should be copied except of the
indices. Note that `idx` can be of type `Int` or `AbstractVector`.
Both options must be supported by a custom type.
If it makes sense for the type of `data`, `obsdim` can be used
to disptach on which dimension of `data` denotes the observations.
"""
function datasubset end

# We don't own nobs but pirate it for basic types
"""
   nobs(data; [obsdim])

Return the number of observations in the dataset `data`.

If it makes sense for the type of `data`, `obsdim` can be used
to indicate which dimension of `data` denotes the observations.
See [`default_obsdim`](@ref) for defining a default dimension.
"""
function StatsBase.nobs(data::AbstractArray; obsdim::Union{Int,Nothing}=nothing)
    od = obsdim === nothing ? default_obsdim(data) : obsdim
    size(data, od)
end

function StatsBase.nobs(data::Union{Tuple, NamedTuple, AbstractDict}; obsdim::Union{Int,Nothing} = default_obsdim(data))
    length(data) > 0 || throw(ArgumentError("Need at least one data input"))

    # We don't force users to handle the obsdim
    # keyword if not necessary.
    fnobs = obsdim === nothing ? nobs : x -> nobs(x; obsdim=obsdim)

    n = fnobs(data[first(keys(data))])
    for i in keys(data)
        ni = fnobs(data[i])
        n == ni || throw(DimensionMismatch("All data inputs should have the same number of observations, i.e. size in the last dimension. "))
    end
    return n
end

# todeprecate
function target end
function gettarget end


text:
using DataFrames, PrettyTables
using StatsBase, StatsBase.Statistics

"""
    function print_decision_strategy(S::States, Z::DecisionStrategy)

Print decision strategy.

# Examples
```julia
print_decision_strategy(S, Z)
```
"""
function print_decision_strategy(S::States, Z::DecisionStrategy)
    for (d, Z_j) in zip(Z.D, Z.Z_j)
        a1 = vec(collect(paths(S[d.I_j])))
        a2 = [Z_j(s_I) for s_I in a1]
        labels = fill("States", length(a1))
        df = DataFrame(labels = labels, a1 = a1, a2 = a2)
        pretty_table(df, ["Nodes", "$((d.I_j...,))", "$(d.j)"])
    end
end

"""
    function print_utility_distribution(udist::UtilityDistribution; util_fmt="%f", prob_fmt="%f")

Print utility distribution.

# Examples
```julia
udist = UtilityDistribution(S, P, U, Z)
print_utility_distribution(udist)
```
"""
function print_utility_distribution(udist::UtilityDistribution; util_fmt="%f", prob_fmt="%f")
    df = DataFrame(Utility = udist.u, Probability = udist.p)
    formatters = (
        ft_printf(util_fmt, [1]),
        ft_printf(prob_fmt, [2]))
    pretty_table(df; formatters = formatters)
end

"""
    function print_state_probabilities(sprobs::StateProbabilities, nodes::Vector{Node}; prob_fmt="%f")

Print state probabilities with fixed states.

# Examples
```julia
sprobs = StateProbabilities(S, P, U, Z)
print_state_probabilities(sprobs, [c.j for c in C])
print_state_probabilities(sprobs, [d.j for d in D])
```
"""
function print_state_probabilities(sprobs::StateProbabilities, nodes::Vector{Node}; prob_fmt="%f")
    probs = sprobs.probs
    fixed = sprobs.fixed

    prob(p, state) = if 1≤state≤length(p) p[state] else NaN end
    fix_state(i) = if i∈keys(fixed) string(fixed[i]) else "" end

    # Maximum number of states
    limit = maximum(length(probs[i]) for i in nodes)
    states = 1:limit
    df = DataFrame()
    df[!, :Node] = nodes
    for state in states
        df[!, Symbol("State $state")] = [prob(probs[i], state) for i in nodes]
    end
    df[!, Symbol("Fixed state")] = [fix_state(i) for i in nodes]
    pretty_table(df; formatters = ft_printf(prob_fmt, (first(states)+1):(last(states)+1)))
end

"""
function print_statistics(udist::UtilityDistribution; fmt = "%f")

Print statistics about utility distribution.
"""
function print_statistics(udist::UtilityDistribution; fmt = "%f")
    u = udist.u
    w = ProbabilityWeights(udist.p)
    names = ["Mean", "Std", "Skewness", "Kurtosis"]
    statistics = [mean(u, w), std(u, w, corrected=false), skewness(u, w), kurtosis(u, w)]
    df = DataFrame(Name = names, Statistics = statistics)
    pretty_table(df, formatters = ft_printf(fmt, [2]))
end

"""
    function print_risk_measures(udist::UtilityDistribution, αs::Vector{Float64}; fmt = "%f")

Print risk measures.
"""
function print_risk_measures(udist::UtilityDistribution, αs::Vector{Float64}; fmt = "%f")
    u, p = udist.u, udist.p
    VaR = [value_at_risk(u, p, α) for α in αs]
    CVaR = [conditional_value_at_risk(u, p, α) for α in αs]
    df = DataFrame(α = αs, VaR = VaR, CVaR = CVaR)
    pretty_table(df, formatters = ft_printf(fmt))
end


text:
__precompile__(true)

"""
The module for [Random123 Family](@ref).

Provide 8 RNG types:

- [`Threefry2x`](@ref)
- [`Threefry4x`](@ref)
- [`Philox2x`](@ref)
- [`Philox4x`](@ref)
- [`AESNI1x`](@ref)
- [`AESNI4x`](@ref)
- [`ARS1x`](@ref)
- [`ARS4x`](@ref)
"""
module Random123

using RandomNumbers

export set_counter!
include("common.jl")

export Threefry2x, Threefry4x
include("threefry.jl")

export Philox2x, Philox4x
include("philox.jl")

export R123_USE_AESNI

"True when AES-NI has been enabled."
const R123_USE_AESNI = try
    cmd = Base.julia_cmd()
    push!(
        cmd.exec, "-e",
        "const __m128i = NTuple{2, VecElement{UInt64}};" *
        "@assert ccall(\"llvm.x86.aesni.aeskeygenassist\", " *
        "llvmcall, __m128i, (__m128i, UInt8), " *
        "__m128i((0x0123456789123450, 0x9876543210987654)), 0x1) ≡ " *
        "__m128i((0x857c266f7c266e85, 0x2346382146382023))"
    )
    success(cmd)
catch e
    false
end

if R123_USE_AESNI
    export AESNI1x, AESNI4x
    export ARS1x, ARS4x
    include("./aesni_common.jl")
    include("./aesni.jl")
    include("./ars.jl")
else
    @warn "AES-NI is not enabled, so AESNI and ARS are not available."
end

end


text:
### contains functions to compute probabilities in subsets
### which we used to call partitions
### the functions are in the process of being rewritten for
### multisets

"""
	matrix_phi(k, U, occupancy_vector, n)
"""
function matrix_phi(k, U, occupancy_vector, n)

	"""matrix P(phi_k) in the notes

	n = number of photons
	occupancy_vector is a Int64 vector of {0,1}^m with one if the mode belongs to the partition, zero otherwise

	remark : opposite sign convention from the notes called notes_singlemode_output
	by default for the r = (1,1, ..., 1) input """

	if k > n
		throw(ArgumentError("k > n"))
	end

	check_at_most_one_particle_per_mode(occupancy_vector)

	m = size(U,1)
	mat = Matrix{eltype(U)}(I, size(U))

	for i in 1:m
		if occupancy_vector[i] == 1
			mat[i,i] = exp(-2im * pi * k / (n+1))
		end
	end

	mat

end

"""
	proba_partition_partial(;U, S, occupancy_vector, input_state, checks=true)

Return a ``n+1`` sized array giving the probability to find `[0,1,...]`, photons
inside the bins given by `occupancy_vector` at the output of `U`.

!!! note
	- We take `U` of dimension ``m`` while `M` is the scattering matrix, as in
	  Tichy, ``M_ij = U_{d_i}, _{d_j}``.
	- Given ``n`` photons, generally in the first modes, the distinguishability
	  matrix is defined as in Tichy, ``S_{ij} = <phi_{d_i}|phi_{d_j}>``.
	  This is not a problem as it does not depend on the output partition but be aware of it.
"""
function proba_partition_partial(; U, S, occupancy_vector, input_state, checks = true)

	"""returns a n+1 sized array giving the probability of have [zero, one, ...]
	photons inside the bins given by occupancy_vector for the interferometer
	given by the matrix U (perfectly indistinguishable photons)

	like proba_partition but with partial distinguishability defined through the S
	matrix, with conventions as in Tichy (note that we take the S matrix to be n*n
	while the interferometer has m modes)

	NOTE :
	in the following, we take U to be m*m
	while M is the scattering matrix, as in Tichy, M_ij = U_{d_i}, _{d_j}
	and have n photons, generally in the first n modes
	the distinguishability matrix is defined as in tichy, to be n*n,
	S_{ij} = <phi_{d_i}|phi_{d_j}>
	this is not a problem as it does not depend on the output partition but be aware of it
	"""

	m = size(U,1)
	n = sum(input_state)

	if n == 0
		throw(ArgumentError("number of photons ill defined, zero cannot be computed but trivial"))
	end

	if size(S,1) != n
		throw(ArgumentError("S matrix does not have the size required for this number of photons"))
	end

	if checks
		if !is_unitary(U)
			throw(ArgumentError("U not unitary"))
		elseif !is_a_gram_matrix(S)
			throw(ArgumentError("S not gram"))
		else
			println("Note : checking unitarity of U, that S is gram, slows down significantly probabilities computations")
		end
	end

	function scattering_matrix_amplitudes(U, input_state, occupancy_vector, k)

	    """U tilde in the handwritten notes, corresponds to U^dagger Lambda(eta) U"""

	    scattering_matrix(U' * matrix_phi(k, U, occupancy_vector, sum(input_state)) * U, input_state, input_state)

	end

	proba_fourier = 1/vector_factorial(input_state) .* [ryser(copy(transpose(S)) .* scattering_matrix_amplitudes(U, input_state, occupancy_vector, k)) for k in 0:n]

	if length(proba_fourier) == 0
		throw(ArgumentError("cannot compute the idft of an empty array"))
	end

	p_recovered = [1/(n+1) * sum([exp(2im * pi * j * k / (n+1)) * proba_fourier[k+1] for k in 0:n]) for j in 0:n] # direct inverse fourier transform

	if any([!isapprox(imag(p_recovered[i]),  0., atol = 1e-7, rtol = 1e-5) for i in 1:length(p_recovered)])
		println("WARNING : proba[$j] has a significant imaginary part of $(imag(tmp))")
	elseif !(isapprox(sum(real.(p_recovered)), 1., rtol = 1e-5))
		println("WARNING : probabilites do not sum to one ($(sum(real.(tmp)))), maybe be erreneous")
	end

	real.(p_recovered)

end

"""
	proba_partition_bosonic(;U, occupancy_vector, input_state=ones(Int, size(U,1)), checks=true)

Indistinguishable version of [`proba_partition_partial`](@ref).
"""
function proba_partition_bosonic(;U, occupancy_vector, input_state = ones(Int, size(U,1)), checks = true)

	"""indistinguishable version of proba_partition_partial"""

	n = sum(input_state)

	proba_partition_partial(U = U , S = ones(n,n), occupancy_vector = occupancy_vector, input_state = input_state, checks = checks)

end

function proba_partition_distinguishable(;U, occupancy_vector, input_state = ones(Int, size(U,1)), checks = true)

	"""indistinguishable version of proba_partition_partial"""

	n = sum(input_state)

	proba_partition_partial(U = U , S = Matrix{eltype(U)}(I,n,n), occupancy_vector = occupancy_vector, input_state = input_state, checks = checks)

end

"""
	partition_probability_distribution_distinguishable_rand_walk(part, U)

Generate a vector giving the probability to have ``k`` photons in the partition
`part` at the output of the interferomter `U`.
"""
function partition_probability_distribution_distinguishable_rand_walk(part, U)

	"""generates a vector giving the probability to have k photons in
	the partition part for the interferometer U by the random walk method
	discussed in a 22/02/21 email
	part is written like (1,2,4) if it is the output modes 1, 2 and 4"""

	function proba_photon_in_partition(i, part, U)
		"""returns the probability that the photon i ends up in the set of outputs part, for distinguishable particles only"""

		sum([abs(U[i, j])^2 for j in part]) # note the inversion line column
	end

	function walk(probability_vector, walk_number, part, U)

		new_probability_vector = similar(probability_vector)
		n = size(U)[1]
		proba_this_walk = proba_photon_in_partition(walk_number, part, U)

		for i in 1 : n+1
			new_probability_vector[i] = (1-proba_this_walk) * probability_vector[i] + proba_this_walk * probability_vector[i != 1 ? i-1 : n+1]
		end

		new_probability_vector

	end

	n = size(U)[1]
	probability_vector = zeros(n+1)
	probability_vector[1] = 1

	for walk_number in 1 : n
		probability_vector = walk(probability_vector, walk_number, part, U)
	end

	if !(isapprox(sum(probability_vector), 1., rtol = 1e-8))
		println("WARNING : probabilites do not sum to one ($(sum(probability_vector)))")
	end

	probability_vector

end


text:
function count(source::Enumerable, filter::Function, filter_expr::Expr)
    return Base.count(filter, source)
end

function count(source::Enumerable)
    return Base.count(i->true, source)
end


text:
#! /usr/bin/julia

# Rosetta Code, Order two numerical lists

function isallreal{T<:AbstractArray}(a::T)
    all(map(x->isa(x, Real), a))
end

function islexfirst{T<:AbstractArray,U<:AbstractArray}(a::T, b::U)
    isallreal(a) && isallreal(b) || throw(DomainError())
    for i in 1:min(length(a), length(b))
        x = a[i]
        y = b[i]
        x != y || continue
        return x < y
    end
    return length(a) < length(b)
end

tests = {[1, 2, 3],
         primes(10),
         0:2:6,
         [-Inf, 0.0, Inf],
         [π, e, φ, catalan],
         [2015, 5],
         [-sqrt(50.0), 50.0^2],
         }

println("Testing islexfirst:")
for (a, b) in combinations(tests, 2)
    tres = islexfirst(a, b) ? " is " : " is not "
    tres *= "lexically prior to\n    "
    println("\n    ", a, tres, b)
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
module TestSwaps
  using Test
  using EcologicalNetworks
  using Random

  N = convert(BinaryNetwork, web_of_life("A_HP_001"))
  Ud = shuffle(N; constraint=:degree)
  Uf = shuffle(N; constraint=:fill)
  Ui = shuffle(N; constraint=:vulnerability)
  Uo = shuffle(N; constraint=:generality)

  @test links(N) == links(Ud)
  @test links(N) == links(Uf)
  @test links(N) == links(Ui)
  @test links(N) == links(Uo)

  @test degree(N) == degree(Ud)
  @test degree(N; dims=1) == degree(Ud; dims=1)
  @test degree(N; dims=2) == degree(Ud; dims=2)
  @test degree(N; dims=1) == degree(Uo; dims=1)
  @test degree(N; dims=2) == degree(Ui; dims=2)

  N = convert(BinaryNetwork, nz_stream_foodweb()[1])
  Ud = shuffle(N; constraint=:degree)
  Uf = shuffle(N; constraint=:fill)
  Ui = shuffle(N; constraint=:vulnerability)
  Uo = shuffle(N; constraint=:generality)

  @test links(N) == links(Ud)
  @test links(N) == links(Uf)
  @test links(N) == links(Ui)
  @test links(N) == links(Uo)

  @test degree(N) == degree(Ud)
  @test degree(N; dims=1) == degree(Ud; dims=1)
  @test degree(N; dims=2) == degree(Ud; dims=2)
  @test degree(N; dims=1) == degree(Uo; dims=1)
  @test degree(N; dims=2) == degree(Ui; dims=2)

end


text:
type DividendSchedule
  dividends::Vector{Dividend}
end

DividendSchedule() = DividendSchedule(Vector{Dividend}(0))

type BondMixin
  settlementDays::Int
  # schedule::Schedule
  # cashflows::L
  # dc::DC
  issueDate::Date
  maturityDate::Date
  # pricingEngine::P
end

get_settlement_days(bond::Bond) = bond.bondMixin.settlementDays
get_issue_date(bond::Bond) = bond.bondMixin.issueDate
get_maturity_date(bond::Bond) = bond.bondMixin.maturityDate

type FixedRateBond{DC <: DayCount, P <: PricingEngine} <: Bond
  lazyMixin::LazyMixin
  bondMixin::BondMixin
  faceAmount::Float64
  schedule::Schedule
  cashflows::FixedRateLeg
  dc::DC
  redemption::Float64
  startDate::Date
  pricingEngine::P
  settlementValue::Float64
end

function FixedRateBond{DC <: DayCount, B <: BusinessDayConvention, C <: BusinessCalendar, P <: PricingEngine}(settlementDays::Int, faceAmount::Float64, schedule::Schedule, coup_rate::Float64, dc::DC, paymentConvention::B,
  redemption::Float64, issueDate::Date, calendar::C, pricing_engine::P)
  maturityDate = schedule.dates[end]

  # num_payments = length(schedule.dates)
  # coups = zeros(num_payments + 1)
  #
  # # build coupons
  # for i = 1:num_payments
  #   fixed_leg = FixedRateCoupon(schedule, faceAmount, coup_rate, calendar, paymentConvention)
  #   coups[i] = fixed_leg
  # end

  # add redemption
  # coups[end] = SimpleCashFlow(redemption, maturityDate)

  coups = FixedRateLeg(schedule, faceAmount, coup_rate, calendar, paymentConvention, dc)

  return FixedRateBond{DC, P}(LazyMixin(), BondMixin(settlementDays, issueDate, maturityDate), faceAmount, schedule, coups, dc, redemption, schedule.dates[1], pricing_engine, 0.0)
end

type FloatingRateBond{X <: InterestRateIndex, DC <: DayCount, P <: PricingEngine} <: Bond
  lazyMixin::LazyMixin
  bondMixin::BondMixin
  faceAmount::Float64
  schedule::Schedule
  cashflows::IborLeg
  iborIndex::X
  dc::DC
  fixingDays::Int
  gearings::Vector{Float64}
  spreads::Vector{Float64}
  caps::Vector{Float64}
  floors::Vector{Float64}
  inArrears::Bool
  redemption::Float64
  pricingEngine::P
  settlementValue::Float64
end

function FloatingRateBond{X <: InterestRateIndex, DC <: DayCount, B <: BusinessDayConvention, P <: PricingEngine}(settlementDays::Int, faceAmount::Float64, schedule::Schedule, iborIndex::X, dc::DC,
                          convention::B, fixingDays::Int, issueDate::Date, pricingEngine::P, inArrears::Bool = false, redemption::Float64 = 100.0,
                          gearings::Vector{Float64} = ones(length(schedule.dates) - 1), spreads::Vector{Float64} = zeros(length(schedule.dates) - 1),
                          caps::Vector{Float64} = Vector{Float64}(), floors::Vector{Float64} = Vector{Float64}();
                          cap_vol::OptionletVolatilityStructure=NullOptionletVolatilityStructure())
  maturityDate = schedule.dates[end]
  fixingDaysVect = fill(fixingDays, length(schedule.dates) - 1)

  coups = IborLeg(schedule, faceAmount, iborIndex, dc, convention, fixingDaysVect, gearings, spreads, caps, floors, inArrears;
                  add_redemption=true, cap_vol = cap_vol)
  return FloatingRateBond{X, DC, P}(LazyMixin(), BondMixin(settlementDays, issueDate, maturityDate), faceAmount, schedule, coups, iborIndex, dc, fixingDays, gearings, spreads, caps, floors,
                          inArrears, redemption, pricingEngine, 0.0)
end

type ZeroCouponBond{BC <: BusinessCalendar, P <: PricingEngine} <: Bond
  lazyMixin::LazyMixin
  bondMixin::BondMixin
  faceAmount::Float64
  redemption::Float64
  cashflows::ZeroCouponLeg
  calendar::BC
  settlementValue::Float64
  pricingEngine::P
end

function ZeroCouponBond{B <: BusinessCalendar, C <: BusinessDayConvention, P <: PricingEngine}(settlementDays::Int, calendar::B, faceAmount::Float64, maturityDate::Date,
                        paymentConvention::C=Following(), redemption::Float64=100.0, issueDate::Date=Date(), pe::P = DiscountingBondEngine())
  # build redemption CashFlow
  redemption_cf = ZeroCouponLeg(SimpleCashFlow(redemption, maturityDate))
  return ZeroCouponBond{B, P}(LazyMixin(), BondMixin(settlementDays, issueDate, maturityDate), faceAmount, redemption, redemption_cf, calendar, 0.0, pe)
end

get_calendar(zcb::ZeroCouponBond) = zcb.calendar

get_calendar(b::Bond) = b.schedule.cal

get_settlement_date(b::Bond) = get_issue_date(b)

function notional(bond::Bond, d::Date)
  if d > get_maturity_date(bond)
    return 0.0
  else
    return bond.faceAmount
  end
end

# Calculation method
function perform_calculations!{B <: Bond}(bond::B)
  bond.settlementValue = 0.0 # reset - TODO this will be expanded
  _calculate!(bond.pricingEngine, bond)

  return bond
end

# bond methods
accrued_amount(bond::Bond, settlement::Date) = accrued_amount(bond.cashflows, settlement, false) * 100.0 / notional(bond, settlement)

maturity_date(bond::Bond) = maturity_date(bond.cashflows)

function yield(bond::Bond, clean_price::Float64, dc::DayCount, compounding::CompoundingType, freq::Frequency, settlement::Date, accuracy::Float64 = 1.0e-10,
              max_iter::Int = 100, guess::Float64 = 0.05)
  dirty_price = clean_price + accrued_amount(bond, settlement)
  dirty_price /= 100.0 / notional(bond, settlement)

  return yield(bond.cashflows, dirty_price, dc, compounding, freq, false, settlement, settlement, accuracy, max_iter, guess)
end

yield(bond::Bond, dc::DayCount, compounding::CompoundingType, freq::Frequency, settlement::Date, accuracy::Float64 = 1.0e-10, max_iter::Int = 100) =
      yield(bond, clean_price(bond), dc, compounding, freq, settlement, accuracy, max_iter)

yield(bond::Bond, dc::DayCount, compounding::CompoundingType, freq::Frequency, accuracy::Float64 = 1.0e-10, max_iter::Int = 100) =
      yield(bond, clean_price(bond), dc, compounding, freq, settlement_date(bond), accuracy, max_iter)

function duration(bond::Bond, yld::InterestRate, duration_::Duration, dc::DayCount, settlement_date::Date)
  return duration(duration_, bond.cashflows, yld, dc, false, settlement_date)
end

function duration(bond::Bond, yld::Float64, dc::DayCount, compounding::CompoundingType, freq::Frequency, duration_::Duration, settlement_date::Date)
  y = InterestRate(yld, dc, compounding, freq)
  return duration(bond, y, duration_, dc, settlement_date)
end

# function npv{B <: Bond, P <: PricingEngine}(bond::B, pe::P)
#   calculate!(pe, bond)
#   return bond.settlementValue
# end
function npv(bond::Bond)
  calculate!(bond)

  return bond.settlementValue
end

clean_price(bond::Bond, settlement_value::Float64, settlement_date::Date) = dirty_price(bond, settlement_value, settlement_date) - accrued_amount(bond, settlement_date)
dirty_price(bond::Bond, settlement_value::Float64, settlement_date::Date) = settlement_value * 100.0 / bond.faceAmount # replace with notionals

function clean_price(bond::Bond)
  calculate!(bond)
  clean_price(bond, bond.settlementValue, settlement_date(bond))
end

dirty_price(bond::Bond) = dirty_price(bond, bond.settlementValue, settlement_date(bond))

function settlement_date(bond::Bond, d::Date = Date())
  if d == Date()
    d = settings.evaluation_date
  end

  # return d + Base.Dates.Day(get_settlement_days(bond))
  return advance(Dates.Day(get_settlement_days(bond)), get_calendar(bond), d)
end

get_redemption(b::Bond) = b.cashflows.coupons[end]
get_frequency(b::Bond) = b.schedule.tenor.freq

## clone methods ##
function clone(bond::FixedRateBond, pe::PricingEngine = bond.pricingEngine)
  # if everything is the same, we keep the calculation status etc
  lazyMixin, settlementValue = pe == bond.pricingEngine ? (bond.lazyMixin, bond.settlementValue) : (LazyMixin(), 0.0)

  return FixedRateBond(lazyMixin, bond.bondMixin, bond.faceAmount, bond.schedule, bond.cashflows, bond.dc, bond.redemption, bond.startDate, pe, settlementValue)
end


text:
using Test
using LinearAlgebra
using ProximalAlgorithms: Broyden, initialize, BroydenOperator, update!, reset!

@testset "Broyden ($R)" for R in [Float32, Float64]
    H = R[
        0.63287    0.330934   -0.156908   -0.294776    0.10761;
        0.330934   0.673201    0.0459778   0.231011   -0.235265;
       -0.156908   0.0459778   0.635812   -0.232261   -0.388775;
       -0.294776   0.231011   -0.232261    0.726854   -0.0691783;
        0.10761   -0.235265   -0.388775   -0.0691783   0.336262;
    ]
    l = R[1.0, 2.0, 3.0, 4.0, 5.0]

    n = length(l)

    f(x) = dot(x, H * x) / 2 + dot(x, l)
    grad_f(x) = H * x + l

    x_star = -H \ l
    f_star = f(x_star)

    Lip = opnorm(H)
    gamma = 1 / Lip
    x = zeros(R, n)

    @inferred initialize(Broyden(), x)

    acc = BroydenOperator(x)

    grad_f_x = grad_f(x)

    for it = 1:10
        d = @inferred acc * grad_f_x
        x = x - d

        grad_f_x_prev = grad_f_x
        grad_f_x = grad_f(x)

        update!(acc, -d, grad_f_x - grad_f_x_prev)
    end

    @test f(x) <= f_star + (1 + abs(f_star)) * sqrt(eps(R))

    reset!(acc)

    @test acc * x == x
end


text:
using Flux.Tracker

function simple_differential()
    f(x)=3x^2 + 2x + 1
    f′(x) = Tracker.gradient(f,x)[1]
    @show f′(2)
    f′′(x) = Tracker.gradient(f′,x)[1]
    @show f′′(2)

    f(W,b,x) = W*x+b
    # returns ∂f/∂W ∂f/∂b ∂f/∂x with respect to (W,b,x)=(2,3,4)
    @show Tracker.gradient(f,2,3,4)
end

function differential_with_params()
    W = param(2) # 2.0 (tracked)
    b = param(3) # 3.0 (tracked)

    f(x) = W * x + b

    params = Params([W, b])
    grads = Tracker.gradient(() -> f(4), params)

    @show grads[W] # 4.0
    @show grads[b] # 1.0
end

function main()
    simple_differential()
    differential_with_params()
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using MDPTools, NestedFixedPoint
using DataFrames

# 1 Data
df = readtable(Pkg.dir("StructuralBase"*"/data/data1234.csv"))
n = 90
omax = 450000
Xval = linspace(0, omax/1000, n)
x = ceil(df[:odo]*n/omax)
data = DataFrame()
data[:a] = df[:replace]
data[:x] = x
data[:dx] = x[1:end]-[0;x[1:end-1]]
for i = eachindex(data[:dx])
    data[:dx][i] = data[:dx][i]*(1-data[:a][i])+data[:x][i]*data[:a][i]
end
remove_first_row_index=df[:id]-[0;df[:id][1:end-1]]
data=data[remove_first_row_index .== 0,:]
data[:a]=[data[:a][2:end];0]
maxi = Int64(maximum(data[:dx]))

# 2 Model
#function transition(p)

p=zeros(maxi+1);
for i = 1:1:maxi+1
    p[i]=sum(data[:dx] .== i-1)/length(data[:dx])
end

F = zeros(n, n)
for i = 1:n-length(p)+1
    F[i,i:i+length(p)-1] = p
end
for i=1:length(p)-1
    F[end-length(p)+i+1,end-length(p)+i+1:end]=[p[1:end-i-1];sum(p[end-i:end])]
end
F = [F, F[ones(Int64, n),:]]

S = State(Xval, F)

Z1 = [zeros(n) -0.001*Xval]
Z2 = [-ones(n) zeros(n)]

U = LinearUtility([Z1, Z2], 0.9999, [11., 3.])

a = vec(round(Int64,data[:a]))
x = vec(round(Int64,data[:x]))

D = Data(a, a+1, x)
U.Θ*=0
rust_nfxp = fit_nfxp(U, S, D; verbose = false);


text:
abstract type AbstractWignerCache end

struct WignerCache <: AbstractWignerCache
    Wigner3j::Dict{Tuple{UInt,UInt,UInt,Int,Int},Tuple{Rational{BigInt},Rational{BigInt}}}
    Wigner6j::Dict{NTuple{6,UInt},Tuple{Rational{BigInt},Rational{BigInt}}}

    primetable::Array{Int64,1}
    factortable::Array{Array{UInt8,1},1}
    factorialtable::Array{Array{UInt32,1},1}
    bigprimetable::Array{Array{BigInt,1},1}
    bigone::Base.RefValue{BigInt}
end

function WignerCache()
    cache = WignerCache(
        Dict{Tuple{UInt,UInt,UInt,Int,Int},Tuple{Rational{BigInt},Rational{BigInt}}}(),
        Dict{NTuple{6,UInt},Tuple{Rational{BigInt},Rational{BigInt}}}(),
        [2,3,5],
        [UInt8[], UInt8[1], UInt8[0,1], UInt8[2], UInt8[0,0,1]],
        [UInt32[], UInt32[], UInt32[1], UInt32[1,1], UInt32[3,1], UInt32[3,1,1]],
        [[big(2)], [big(3)], [big(5)]],
        Ref{BigInt}(big(1))
    )
    return cache
end


text:
#clzctz.jl
#leading_zeros and trailing_zeros operations, stored as global function variables.

doc"""
  `clz(::UInt64)` and 'clz(::ArrayNum)' count the leading zeros and return a
  UInt16 value (instead of the Int64 standard value for leading_zeros.)
"""
clz(n::UInt64) = to16(leading_zeros(n))  #NB:  This should be shimmed with a 'fast' version
                                           #that goes directly to UInt16
clz(n::UInt16) = to16(leading_zeros(n))
function clz{FSS}(n::ArrayNum{FSS})
  res::UInt16 = z16
  cellvalue::UInt64 = z64
  #iterate down the array starting from the most significant cell
  for idx = 1:__cell_length(FSS)
    @inbounds (cellvalue = n.a[idx])
    res += clz(cellvalue)
    (cellvalue != z64) && return res
  end
  res
end

doc"""
  `ctz(::UInt64)` and 'ctz(::ArrayNum)' count the trailing zeros and return a
  UInt16 value (instead of the Int64 standard value for trailing_zeros.)
"""
ctz(n::UInt64) = to16(trailing_zeros(n))
#for when it's a superint (that's not a straight Uint64)
function ctz{FSS}(n::ArrayNum{FSS})
  res::UInt16 = z16
  should_continue::Bool = true
  cellvalue::UInt64 = z64
  #iterate down the array starting from the least significant cell (highest index)
  for idx = __cell_length(FSS):-1:1
    @inbounds (cellvalue = n.a[idx])
    res += ctz(cellvalue)
    (cellvalue != z64) && return res
  end
  res
end

doc"""
  `cto(::UInt64)` and 'cto(::ArrayNum)' count the trailing ones and return a
  UInt16 value (instead of the Int64 standard value for trailing_zeros.)
"""
cto(n::UInt64) = to16(trailing_ones(n))
#for when it's a superint (that's not a straight Uint64)
function cto{FSS}(n::ArrayNum{FSS})
  res::UInt16 = z16
  should_continue::Bool = true
  cellvalue::UInt64 = z64
  #iterate down the array starting from the least significant cell (highest index)
  for idx = __cell_length(FSS):-1:1
    @inbounds (cellvalue = n.a[idx])
    res += cto(cellvalue)
    (cellvalue != f64) && return res
  end
  res
end

export clz, ctz, cto

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl