text:
module SimplexFacesTests

using Test
using STLCutters
using Gridap
using Gridap.ReferenceFEs

using STLCutters: center
using STLCutters: distance
using STLCutters: projection
using STLCutters: measure
using STLCutters: normal
using STLCutters: signed_distance
using STLCutters: contains_projection
using STLCutters: voxel_intersection
using STLCutters: simplex_face
using STLCutters: min_height
using STLCutters: Plane

# Face{1,2}

p1 = Point(0.0,1.0)
p2 = Point(1.0,0.0)
f = simplex_face(p1,p2)

@test num_dims(f) == 1
@test num_point_dims(f) == 2
@test center(f) == Point(0.5,0.5)
@test measure(f) ≈ √2
@test normal(f) ≈ (√2)/2 * Point(1,1)
@test normal(f,1) ≈ (√2)/2 * Point(-1,1)
@test normal(f,2) ≈ (√2)/2 * Point(1,-1)

pmin,pmax = get_bounding_box(f)

@test measure(pmin,pmax) == 1

# Face{1,3}

p1 = Point(0.0,1.0,0.0)
p2 = Point(1.0,0.0,1.0)
f = simplex_face(p1,p2)

@test num_dims(f) == 1
@test num_point_dims(f) == 3
@test center(f) == Point(0.5,0.5,0.5)
@test measure(f) ≈ √3
@test normal(f,1) ≈ 1/(√3) * Point(-1,1,-1)
@test normal(f,2) ≈ 1/(√3) * Point(1,-1,1)

pmin,pmax = get_bounding_box(f)

@test measure(pmin,pmax) == 1

# Face{2,2}

p1 = Point(0.0,1.0)
p2 = Point(1.0,0.0)
p3 = Point(0.0,0.0)
f = simplex_face(p1,p2,p3)

@test num_dims(f) == 2
@test num_point_dims(f) == 2
@test center(f) ≈ Point(1/3,1/3)
@test measure(f) ≈ 1/2
@test min_height(f) ≈ 1/√2

pmin,pmax = get_bounding_box(f)

@test measure(pmin,pmax) == 1

# Face{2,3}

p1 = Point(0.0,1.0,0.0)
p2 = Point(1.0,0.0,0.0)
p3 = Point(0.0,0.0,1.0)
f = simplex_face(p1,p2,p3)

@test center(f) ≈ Point(1/3,1/3,1/3)
@test measure(f) ≈ (√3)/2

@test num_dims(f) == 2
@test num_point_dims(f) == 3
@test normal(f) ≈ 1/√3 * Point(-1,-1,-1)
@test min_height(f) ≈ √(3/2)

pmin,pmax = get_bounding_box(f)

@test measure(pmin,pmax) == 1

# Face{3,3}

p1 = Point(0.0,0.0,0.0)
p2 = Point(1.0,0.0,0.0)
p3 = Point(0.0,1.0,0.0)
p4 = Point(0.0,0.0,1.0)
f = simplex_face(p1,p2,p3,p4)

@test center(f) ≈ Point(1/4,1/4,1/4)
@test measure(f) ≈ 1/6
@test min_height(f) ≈ 1/√3

pmin,pmax = get_bounding_box(f)

@test measure(pmin,pmax) == 1

# Plane{2}

p1 = Point(0.0,1.0)
p2 = Point(1.0,0.0)
f = simplex_face(p1,p2)

Π = Plane(f)
p = Point(0.0,0.0)

@test signed_distance(p,Π) ≈ -1/√2

# Plane{3}

p1 = Point(0.0,1.0,0.0)
p2 = Point(1.0,0.0,0.0)
p3 = Point(0.0,0.0,1.0)
f = simplex_face(p1,p2,p3)

Π = Plane(f)
p = Point(0.0,0.0,0.0)

@test signed_distance(p,Π) ≈ 1/√3

# Distance Point - Point

p1 = Point(0.0,0.0)
p2 = Point(1.0,1.0)

distance(p1,p2) ≈ √2

p1 = Point(0.0,0.0,0.0)
p2 = Point(1.0,1.0,1.0)

distance(p1,p2) ≈ √3

# Distance Point - Face{1}

p1 = Point(0.0,1.0)
p2 = Point(1.0,0.0)
f = simplex_face(p1,p2)

a = Point(0.0,0.0)
b = Point(2.0,0.0)
c = Point(0.5,0.5)

@test distance(a,f) ≈ 1/√(2)
@test distance(b,f) ≈ 1.0
@test abs( distance(c,f) ) < 1e-9

# Distance Point - Face{2}

p1 = Point(0.0,1.0,0.0)
p2 = Point(1.0,0.0,0.0)
p3 = Point(0.0,0.0,1.0)
f = simplex_face(p1,p2,p3)

a = Point(0.0,0.0,0.0)
b = Point(2.0,0.0,0.0)
c = Point(1/3,1/3,1/3)

@test distance(a,f) ≈ 1/√(3)
@test distance(b,f) ≈ 1.0
@test abs( distance(c,f) ) < 1e-9

# Distance Point - Face{3}

p1 = Point(0.0,0.0,0.0)
p2 = Point(1.0,0.0,0.0)
p3 = Point(0.0,1.0,0.0)
p4 = Point(0.0,0.0,1.0)
f = simplex_face(p1,p2,p3,p4)

a = Point(0.0,0.0,0.0)
b = Point(2.0,0.0,0.0)
c = Point(1/4,1/4,1/4)
d = Point(1.0,1.0,1.0)

@test abs(distance(a,f)) < 1e-9
@test distance(b,f) ≈ 1.0
@test distance(c,f) == 0
@test distance(d,f) ≈ 2/√3

# Projection Point - Point

a = Point(0.0,1.0)
b = Point(1.0,0.0)

@test projection(a,b) == b

# Projection Point - Face{1}

p1 = Point(0.0,1.0)
p2 = Point(1.0,0.0)
f = simplex_face(p1,p2)

a = Point(0.0,0.0)
b = Point(0.5,0.0)

@test projection(a,f) ≈ Point(0.5,0.5)
@test projection(b,f) ≈ Point(0.75,0.25)

# Projection Point - Face{2}

p1 = Point(0.0,1.0,0.0)
p2 = Point(1.0,0.0,0.0)
p3 = Point(0.0,0.0,1.0)
f = simplex_face(p1,p2,p3)

a = Point(0.0,0.0,0.0)
@test projection(a,f) ≈ Point(1/3,1/3,1/3)

# Contains projection Point - Face{1}

p1 = Point(0.0,1.0)
p2 = Point(1.0,0.0)
f = simplex_face(p1,p2)

a = Point(0.0,0.0)
b = Point(2.0,0.0)

@test contains_projection(f,a)
@test !contains_projection(f,b)

# Contains projection Point - Face{2}

p1 = Point(0.0,1.0,0.0)
p2 = Point(1.0,0.0,0.0)
p3 = Point(0.0,0.0,1.0)
f = simplex_face(p1,p2,p3)

a = Point(0.0,0.0,0.0)
b = Point(2.0,0.0,0.0)

@test contains_projection(f,a)
@test !contains_projection(f,b)

# Contains projection Point - Face{3}

p1 = Point(0.0,0.0,0.0)
p2 = Point(1.0,0.0,0.0)
p3 = Point(0.0,1.0,0.0)
p4 = Point(0.0,0.0,1.0)
f = simplex_face(p1,p2,p3,p4)

a = Point(1/4,1/4,1/4)
b = Point(1.0,1.0,1.0)

@test contains_projection(f,a)
@test !contains_projection(f,b)

# Voxel intersection Face{2}

p = QUAD

pmin,pmax = Point(0.0,0.0),Point(1.0,1.0)
pmin -= 0.1
pmax += 0.1

p1 = Point(0.5,0.5)
p2 = Point(2.0,5.0)
a = simplex_face(p1,p2)

p1 = Point(0.0,1.0)
p2 = Point(1.0,0.0)
b = simplex_face(p1,p2)

p1 = Point(0.0,2.0)
p2 = Point(2.0,0.0)
c = simplex_face(p1,p2)

p1 = Point(1.0,2.0)
p2 = Point(2.0,1.0)
d = simplex_face(p1,p2)

p1 = Point(1.5,0.5)
p2 = Point(2.0,0.0)
e = simplex_face(p1,p2)

@test voxel_intersection(a,pmin,pmax,p)
@test voxel_intersection(b,pmin,pmax,p)
@test voxel_intersection(c,pmin,pmax,p)
@test !voxel_intersection(d,pmin,pmax,p)
@test !voxel_intersection(e,pmin,pmax,p)

# Voxel intersection Face{3}

p = HEX

pmin,pmax = Point(0.0,0.0,0.0),Point(1.0,1.0,1.0)
pmin -= 0.1
pmax += 0.1

p1 = Point(1.0,0.0,0.0)
p2 = Point(0.0,1.0,0.0)
p3 = Point(0.0,0.0,1.0)
a = simplex_face(p1,p2,p3)

p1 = Point(9.0,-11.0,-11.0)
p2 = Point(-11.0,9.0,9.0)
p3 = Point(10.0,10.0,10.0)
b = simplex_face(p1,p2,p3)

p1 = Point(0.5,0.5,0.5)
p2 = Point(0.0,2.0,0.0)
p3 = Point(0.0,0.0,2.0)
c = simplex_face(p1,p2,p3)

p1 = Point(4.0,0.0,0.0)
p2 = Point(0.0,4.0,0.0)
p3 = Point(0.0,0.0,4.0)
d = simplex_face(p1,p2,p3)

p1 = Point(2.0,-0.2,-0.2)
p2 = Point(1.5,0.5,0.0)
p3 = Point(1.5,0.0,0.5)
e = simplex_face(p1,p2,p3)

@test voxel_intersection(a,pmin,pmax,p)
@test voxel_intersection(b,pmin,pmax,p)
@test voxel_intersection(c,pmin,pmax,p)
@test !voxel_intersection(d,pmin,pmax,p)
@test !voxel_intersection(e,pmin,pmax,p)

end # module


text:
#= =====================================================
LINK patterns, see html link fixer
===================================================== =#
# here we're looking for [id] or [id][] or [stuff][id] or ![stuff][id] but not [id]:
# 1 > (!)? == either ! or nothing
# 2 > [(.*?)] == [...] inside of the brackets
# 3 > (?:[(.*?)])? == [...] inside of second brackets if there is such
const ESC_LINK_PAT = r"(!)?[(.*?)](?!:)(?:[(.*?)])?"

#= =====================================================
HBLOCK patterns, see html blocks
NOTE: the for block needs verification (matching parens)
===================================================== =#

const HBLOCK_IF_PAT        = r"{{\s*if\s+([a-zA-Z_]\S*)\s*}}"
const HBLOCK_ELSE_PAT      = r"{{\s*else\s*}}"
const HBLOCK_ELSEIF_PAT    = r"{{\s*else\s*if\s+([a-zA-Z_]\S*)\s*}}"
const HBLOCK_END_PAT       = r"{{\s*end\s*}}"

const HBLOCK_ISDEF_PAT     = r"{{\s*i(?:s|f)def\s+([a-zA-Z_]\S*)\s*}}"
const HBLOCK_ISNOTDEF_PAT  = r"{{\s*i(?:s|f)n(?:ot)?def\s+([a-zA-Z_]\S*)\s*}}"
const HBLOCK_ISPAGE_PAT    = r"{{\s*ispage\s+((.|\n)+?)}}"
const HBLOCK_ISNOTPAGE_PAT = r"{{\s*isnotpage\s+((.|\n)+?)}}"

const HBLOCK_FOR_PAT = r"{{\s*for\s+(\(?(?:\s*[a-zA-Z_][^\r\n\t\f\v,]*,\s*)*[a-zA-Z_]\S*\s*\)?)\s+in\s+([a-zA-Z_]\S*)\s*}}"

const HBLOCK_TOC_PAT = r"{{\s*toc\s*}}"

#= =====================================================
Pattern checkers
===================================================== =#

"""
    check_for_pat(v)

Check that we have something like `{{for v in iterate}}` or
`{for (v1,v2) in iterate}}` but not something with unmached parens.
"""
function check_for_pat(v)
    op = startswith(v, "(")
    cp = endswith(v, ")")
    xor(op, cp) &&
        throw(HTMLBlockError("Unbalanced expression in {{for ...}}"))
    !op && occursin(",", v) &&
        throw(HTMLBlockError("Missing parens in {{for ...}}"))
    return nothing
end

ext: jl


text:
"""
    g = gaussianize(z)

Gaussianize field `z` using a sorted random sequence.
"""
function gaussianize(z)
    (m,n) = size(z)
    nn = m*n
    zr = z[:]
    zs = sort(zr)
    i = sortperm(zr)
    rrk = zeros(Int16,nn,1)
    rrk[i] = collect(1:nn)
    gg = randn(nn)
    w = sort(gg)
    gr = w[rrk]
    gs = Statistics.std(gr)
    if gs==0
        gs=1
    end
    gr = gr/gs
    g = reshape(gr,m,n)
end


text:
## sample space/domain

abstract VariateForm
type Univariate    <: VariateForm end
type Multivariate  <: VariateForm end
type Matrixvariate <: VariateForm end

abstract ValueSupport
type Discrete   <: ValueSupport end
type Continuous <: ValueSupport end

Base.eltype(::Type{Discrete}) = Int
Base.eltype(::Type{Continuous}) = Float64

## Sampleable

abstract Sampleable{F<:VariateForm,S<:ValueSupport}

Base.length(::Sampleable{Univariate}) = 1
Base.length(s::Sampleable{Multivariate}) = throw(MethodError(length, (s,)))
Base.length(s::Sampleable) = prod(size(s))

Base.size(s::Sampleable{Univariate}) = ()
Base.size(s::Sampleable{Multivariate}) = (length(s),)

Base.eltype{F,S}(s::Sampleable{F,S}) = eltype(S)
Base.eltype{F}(s::Sampleable{F,Discrete}) = Int
Base.eltype{F}(s::Sampleable{F,Continuous}) = Float64

nsamples{D<:Sampleable{Univariate}}(::Type{D}, x::Number) = 1
nsamples{D<:Sampleable{Univariate}}(::Type{D}, x::AbstractArray) = length(x)
nsamples{D<:Sampleable{Multivariate}}(::Type{D}, x::AbstractVector) = 1
nsamples{D<:Sampleable{Multivariate}}(::Type{D}, x::AbstractMatrix) = size(x, 2)
nsamples{D<:Sampleable{Matrixvariate}}(::Type{D}, x::Number) = 1
nsamples{D<:Sampleable{Matrixvariate},T<:Number}(::Type{D}, x::Array{Matrix{T}}) = length(x)

## Distribution <: Sampleable

abstract Distribution{F<:VariateForm,S<:ValueSupport} <: Sampleable{F,S}

typealias UnivariateDistribution{S<:ValueSupport}   Distribution{Univariate,S}
typealias MultivariateDistribution{S<:ValueSupport} Distribution{Multivariate,S}
typealias MatrixDistribution{S<:ValueSupport}       Distribution{Matrixvariate,S}
typealias NonMatrixDistribution Union(UnivariateDistribution, MultivariateDistribution)

typealias DiscreteDistribution{F<:VariateForm}   Distribution{F,Discrete}
typealias ContinuousDistribution{F<:VariateForm} Distribution{F,Continuous}

typealias DiscreteUnivariateDistribution     Distribution{Univariate,    Discrete}
typealias ContinuousUnivariateDistribution   Distribution{Univariate,    Continuous}
typealias DiscreteMultivariateDistribution   Distribution{Multivariate,  Discrete}
typealias ContinuousMultivariateDistribution Distribution{Multivariate,  Continuous}
typealias DiscreteMatrixDistribution         Distribution{Matrixvariate, Discrete}
typealias ContinuousMatrixDistribution       Distribution{Matrixvariate, Continuous}

variate_form{VF<:VariateForm,VS<:ValueSupport}(::Type{Distribution{VF,VS}}) = VF
variate_form{T<:Distribution}(::Type{T}) = variate_form(super(T))

value_support{VF<:VariateForm,VS<:ValueSupport}(::Type{Distribution{VF,VS}}) = VS
value_support{T<:Distribution}(::Type{T}) = value_support(super(T))

## TODO: the following types need to be improved
abstract SufficientStats
abstract IncompleteDistribution

typealias DistributionType{D<:Distribution} Type{D}
typealias IncompleteFormulation Union(DistributionType,IncompleteDistribution)


text:
module MLJFlux

export CUDALibs, CPU1

import Flux
using MLJModelInterface
using MLJModelInterface.ScientificTypesBase
import Base.==
using ProgressMeter
using CategoricalArrays
using Tables
using Statistics
using ColorTypes
using ComputationalResources
using Random

include("penalizers.jl")
include("core.jl")
include("builders.jl")
include("types.jl")
include("regressor.jl")
include("classifier.jl")
include("image.jl")
include("mlj_model_interface.jl")

### Package specific model traits:
MLJModelInterface.metadata_pkg.((NeuralNetworkRegressor,
                                 MultitargetNeuralNetworkRegressor,
                                 NeuralNetworkClassifier,
                                 ImageClassifier),
              name="MLJFlux",
              uuid="094fc8d1-fd35-5302-93ea-dabda2abf845",
              url="https://github.com/alan-turing-institute/MLJFlux.jl",
              julia=true,
              license="MIT")

export NeuralNetworkRegressor, MultitargetNeuralNetworkRegressor
export NeuralNetworkClassifier, ImageClassifier

end #module


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

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

module Hyper1F1

using Nemo

export ModuleHyper1F1
export GammaHyp, V000255, V000262, V001339, V007060, V033815, V099022, V251568

"""
GammaHyp: ``(a, b, c, d)`` ↦ ``Γ(a) `` Hypergeometric``1F1(b, c, d).``

* GammaHyp, V000255, V000262, V001339, V007060, V033815, V099022, V251568
"""
const ModuleHyper1F1 = ""

# Numerical evaluation based on hypergeometric functions.
# Nemo.hyp1f1(a::acb, b::acb, x::acb)
# May fail if the required precicison is > 10000.
# A Python implementation is here:
# http://luschny.de/math/seq/GammaHypPython.html

"""
Return ``Γ(a) `` Hypergeometric``1F1(b, c, d).``
"""
function GammaHyp(a, b, c, d)
    prec = 64
    while prec <= 10000
        CC = AcbField(prec)
        g = gamma(CC(a)) * hyp1f1(CC(b), CC(c), CC(d))
        b, i = unique_integer(g)
        b && return i
        prec = div(8*prec, 5)
    end

    error("GammaHyp with $a $b $c $d gives an InexactError!")
end

"""
Return ``(n+1)!`` Hypergeometric1F1``[-n, -n-1, -1]``. Number of fixedpoint-free permutations beginning with 2. (L. Euler).
"""
V000255(n) = n == 0 ? ZZ(1) : GammaHyp(n + 2, -n, -n-1, -1)

"""
Return ``n!`` Hypergeometric1F1``[1-n, 2, -1]``. Number of partitions of ``{1,...,n}`` into any number of ordered subsets.
"""
V000262(n) = n == 0 ? ZZ(1) : GammaHyp(n + 1, 1 - n, 2, -1)

"""
Return ``(n+1)!`` Hypergeometric1F1``[-n, -n-1, 1]``. Number of arrangements of ``{1, 2, ..., n+1}`` containing the element 1.
"""
V001339(n) = n == 0 ? ZZ(1) : GammaHyp(n + 2, -n, -n-1, 1)

"""
Return ``(2n)!`` Hypergeometric1F1``[-n, -2n, -2]``. Number of ways ``n`` couples can sit in a row without any spouses next to each other.
"""
V007060(n) = n == 0 ? ZZ(1) : GammaHyp(2n + 1, -n, -2n, -2)

"""
Return ``(2n)!`` Hypergeometric1F1``[-n, -2n, -1]``. Number of acyclic orientations of the Turán graph ``T(2n,n)``.
"""
V033815(n) = n == 0 ? ZZ(1) : GammaHyp(2n + 1, -n, -2n, -1)

"""
Return ``(2n)!`` Hypergeometric1F1``[-n, -2n, 1]``.
"""
V099022(n) = n == 0 ? ZZ(1) : GammaHyp(2n + 1, -n, -2n, 1)

"""
Return ``((2n)!/(n+1)!)`` Hypergeometric1F1``[1-n, n+2, -1]``. Egf. exp ``(x C(x)^2)`` where ``C(x) = 1 + xC(x)^2`` is the generating function of the Catalan numbers.
"""
function V251568(n::Int)
    n == 0 && return fmpz(1)
    prec = 64
    while prec <= 10000
        CC = ComplexField(prec)
        c = gamma(CC(2 * n + 1)) * hyp1f1(CC(1 - n), CC(n + 2), CC(-1)) / gamma(CC(n + 2))
        b, i = unique_integer(c)
        b && return i
        prec *= 2
    end
    error("n = $n gives an InexactError!")
end

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

using Test, SeqTests

function test()

    @testset "Hyper1F1" begin
        @test V000255(10) == ZZ(16019531)
        @test V000262(10) == ZZ(58941091)
        @test V001339(10) == ZZ(98641011)
        @test V007060(10) == ZZ(871804170613555200)
        @test V033815(10) == ZZ(1465957162768492800)
        @test V099022(10) == ZZ(3984884716852972800)
        @test V251568(10) == ZZ(123320412181)

        if is_oeis_installed()
            V = [V000255, V000262, V001339, V007060, V033815,
                 V099022, V251568]
            for v in V SeqTest(v, 'V') end
        end
    end
end

function demo()
    a, b, c, d =  11,  -9,     2, -1; GammaHyp(a, b, c, d) |> println
    a, b, c, d =  21,  -10,  -20, -2; GammaHyp(a, b, c, d) |> println
    a, b, c, d =  21,  -10,  -20, -1; GammaHyp(a, b, c, d) |> println
    a, b, c, d =  21,  -10,  -20,  1; GammaHyp(a, b, c, d) |> println
    a, b, c, d =  74,  -72,    2, -1; GammaHyp(a, b, c, d) |> println
    a, b, c, d = 147,  -73, -146, -2; GammaHyp(a, b, c, d) |> println
    a, b, c, d = 147,  -73, -146, -1; GammaHyp(a, b, c, d) |> println
    a, b, c, d = 147,  -73, -146,  1; GammaHyp(a, b, c, d) |> println
    a, b, c, d = 151, -149,    2, -1; GammaHyp(a, b, c, d) |> println
    a, b, c, d = 301, -150, -300, -2; GammaHyp(a, b, c, d) |> println
    a, b, c, d = 301, -150, -300, -1; GammaHyp(a, b, c, d) |> println
    a, b, c, d = 301, -150, -300,  1; GammaHyp(a, b, c, d) |> println
end

"""
(for n in 0:150 V000262(n) end) ::
    0.010073 seconds (3.75 k allocations: 166.422 KB)
"""
function perf()
    GC.gc()
    @time (for n in 0:150 V000262(n) end)
end

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

main()

end # module


text:
# todo: skip not assigned fields

"""Get an overview of the fields from MERA composite types:
```julia
viewfields(object)
```
"""
function viewfields(object::InfoType)

    list_field = propertynames(object)

    for i=list_field
        if i == :scale
            println()
            printstyled(i, " ==> subfields: ", propertynames(object.scale), "\n", bold=true, color=:normal)
            println()
        elseif i == :grid_info
            printstyled(i, " ==> subfields: ", propertynames(object.grid_info), "\n", bold=true, color=:normal)
            println()
        elseif i == :part_info
            printstyled(i, " ==> subfields: ", propertynames(object.part_info), "\n", bold=true, color=:normal)
            println()
        elseif i == :compilation
            printstyled(i, " ==> subfields: ", propertynames(object.compilation), "\n", bold=true, color=:normal)
            println()
        elseif i == :constants
            printstyled(i, " ==> subfields: ", propertynames(object.constants), "\n", bold=true, color=:normal)
            println()
        elseif i == :fnames
            printstyled(i, " ==> subfields: ", propertynames(object.fnames), "\n", bold=true, color=:normal)
            println()
        elseif i == :descriptor
            printstyled(i, " ==> subfields: ", propertynames(object.descriptor), "\n", bold=true, color=:normal)
            println()
        elseif i == :namelist_content
            printstyled(i, " ==> dictionary: ", tuple(keys(object.namelist_content)...), "\n", bold=true, color=:normal)
            println()
        elseif i == :files_content
            printstyled(i, " ==> subfields: ", propertynames(object.files_content), "\n", bold=true, color=:normal)
            println()
        else
            println(i, "\t= ", getfield(object, i))
        end
    end

    println()
    return
end

function viewfields(object::ScalesType)

    list_field = propertynames(object)

    println()
    printstyled("[Mera]: Fields to scale from user/code units to selected units\n", bold=true, color=:normal)
    printstyled("=======================================================================\n", bold=true, color=:normal)
    for i=list_field
        #fname = fieldname(field, i)
        println(i, "\t= ", getfield(object, i)) #, "\t\t", "[",typeof(getfield(object, i)),"]")
    end
    println()
    return
end

function viewfields(object::PartInfoType)

    list_field = propertynames(object)

    println()
    printstyled("[Mera]: Particle overview\n", bold=true, color=:normal)
    printstyled("===============================\n", bold=true, color=:normal)
    for i=list_field
        #fname = fieldname(field, i)
        println(i, "\t= ", getfield(object, i)) #, "\t\t", "[",typeof(getfield(object, i)),"]")
    end
    println()
    return
end

function viewfields(object::GridInfoType)

    list_field = propertynames(object)

    println()
    printstyled("[Mera]: Grid overview \n", bold=true, color=:normal)
    printstyled("============================\n", bold=true, color=:normal)
    for i=list_field
        if i == :ngridmax || i == :nstep_coarse || i == :nx || i == :ny ||
            i == :nz || i == :nlevelmax || i == :nboundary ||
            i == :ngrid_current
            println(i, "\t= ", getfield(object, i) )

        elseif i == :bound_key
            array_length = length(object.bound_key)
            println(i, " ==> length($array_length)" )
        elseif i == :cpu_read
            array_length = length(object.bound_key)
            println(i, " ==> length($array_length)" )
        end
    end

    println()
    return
end

function viewfields(object::CompilationInfoType)

    list_field = propertynames(object)

    println()
    printstyled("[Mera]: Compilation file overview\n", bold=true, color=:normal)
    printstyled("========================================\n", bold=true, color=:normal)
    for i=list_field
        #fname = fieldname(field, i)
        println(i, "\t= ", getfield(object, i)) #, "\t\t", "[",typeof(getfield(object, i)),"]")
    end
    println()
    return
end

function viewfields(object::FileNamesType)

    list_field = propertynames(object)

    println()
    printstyled("[Mera]: Paths and file-names\n", bold=true, color=:normal)
    printstyled("=================================\n", bold=true, color=:normal)
            for i=list_field
                println(i, "\t= ", getfield(object, i) )
            end

    println()
    return
end

function viewfields(object::DescriptorType)

    list_field = propertynames(object)

    println()
    printstyled("[Mera]: Descriptor overview\n", bold=true, color=:normal)
    printstyled("=================================\n", bold=true, color=:normal)
            for i=list_field
                println(i, "\t= ", getfield(object, i) )
            end

    println()
    return
end

function viewfields(object::FilesContentType)

    list_field = propertynames(object)

    println()
    printstyled("[Mera]: List of files-content\n", bold=true, color=:normal)
    printstyled("=================================\n", bold=true, color=:normal)
    println(propertynames(object))

    println()
    return
end

function viewfields(object::PhysicalUnitsType)

    list_field = propertynames(object)

    println()
    printstyled("[Mera]: Constants given in cgs units\n", bold=true, color=:normal)
    printstyled("=========================================\n", bold=true, color=:normal)
            for i=list_field
                println(i, "\t= ", getfield(object, i) )
            end

    println()
    return
end

# todo: check
function viewfields(object::DataSetType)
    list_field = propertynames(object)
    println()
    for i=list_field
        if i== :data
            printstyled(i, " ==> JuliaDB table: ", colnames(object.data), "\n", bold=true, color=:normal)
            println()
        elseif i== :info
            printstyled(i, " ==> subfields: ", propertynames(object.info), "\n", bold=true, color=:normal)
            println()
        elseif i== :scale
            println()
            printstyled(i, " ==> subfields: ", propertynames(object.scale), "\n", bold=true, color=:normal)
            println()
        elseif i == :used_descriptors
            if object.info.descriptor.usehydro
                println(i, "\t= ", getfield(object, i) )
            end
        else
            println(i, "\t= ", getfield(object, i) )
        end

    end
    println()
    return
end

function namelist(object::InfoType)
    println()
    printstyled("[Mera]: Namelist file content\n", bold=true, color=:normal)
    printstyled("=================================\n", bold=true, color=:normal)
    keylist_header = keys(object.namelist_content)
    for i in keylist_header
        println(i)
        icontent = object.namelist_content[i]
        keylist_parameters = keys(icontent)
        for j in keylist_parameters
            println(j, "  \t=", icontent[j] )
        end
    println()
    end

    return
end

function namelist(object::Dict{Any,Any})
    println()
    printstyled("[Mera]: Namelist file content\n", bold=true, color=:normal)
    printstyled("=================================\n", bold=true, color=:normal)
    keylist_header = keys(object)
    for i in keylist_header
        println(i)
        icontent = object[i]
        keylist_parameters = keys(icontent)
        for j in keylist_parameters
            println(j, "  \t=", icontent[j] )
        end
    println()
    end
    println()
    return
end

"""Get a printout of the makefile:
```julia
makefile(object::InfoType)
```
"""
function makefile(object::InfoType)
    println()
    printstyled("[Mera]: Makefile content\n", bold=true, color=:normal)
    printstyled("=================================\n", bold=true, color=:normal)

    if object.makefile
        for i in object.files_content.makefile
            println(i)
        end
    else
        println("[Mera]: No Makefile found!")
    end

    println()

    return
end

"""Get a printout of the timerfile:
```julia
timerfile(object::InfoType)
```
"""
function timerfile(object::InfoType)
    println()
    printstyled("[Mera]: Timer-file content\n", bold=true, color=:normal)
    printstyled("=================================\n", bold=true, color=:normal)

    if object.timerfile
        for i in object.files_content.timerfile
            println(i)
        end
    else
        println("[Mera]: No timer-file found!")
    end

    println()

    return
end

"""Get a printout of the patchfile:
```julia
patchfile(object::InfoType)
```
"""
function patchfile(object::InfoType)
    println()
    printstyled("[Mera]: Patch-file content\n", bold=true, color=:normal)
    printstyled("=================================\n", bold=true, color=:normal)

    if object.patchfile
        for i in object.files_content.patchfile
            println(i)
        end
    else
        println("[Mera]: No timer-file found!")
    end

    println()

    return
end

"""Get a detailed overview of many fields from the MERA InfoType:
```julia
viewallfields(dataobject::InfoType)
```
"""
function viewallfields(dataobject::InfoType)
    viewfields(dataobject)

    viewfields(dataobject.scale)
    viewfields(dataobject.constants)
    viewfields(dataobject.fnames)
    viewfields(dataobject.descriptor)
    namelist(dataobject)
    viewfields(dataobject.grid_info)
    viewfields(dataobject.part_info)
    viewfields(dataobject.compilation)
    makefile(dataobject)
    timerfile(dataobject)
    return
end


text:
function runsim(numworkers, runs)
    for count in 1:runs
        @sync begin
            for worker in 1:numworkers
                @async begin
                    tasktime = rand()
                    sleep(tasktime)
                    println("Worker $worker finished after $tasktime seconds")
                end
            end
        end
        println("Checkpoint reached for run $count.")
    end
    println("Finished all runs.\n")
end

const trials = [[3, 2], [4, 1], [2, 5], [7, 6]]
for trial in trials
    runsim(trial[1], trial[2])
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using UnsafeArrays
import Base: showarg, eltype
const DSYEVR_ = (BLAS.@blasfunc(dsyevr_),Base.liblapack_name)
const SSYEVR_ = (BLAS.@blasfunc(ssyevr_),Base.liblapack_name)

# ----------------------------------------------------
# Zero cone
# ----------------------------------------------------
"""
    ZeroSet(dim)

Creates the zero set ``\\{ 0 \\}^{dim}`` of dimension `dim`. If `x` ∈ `ZeroSet` then all entries of x are zero.
"""
struct ZeroSet{T} <: AbstractConvexCone{T}
	dim::Int
	function ZeroSet{T}(dim::Int) where {T}
		dim >= 0 ? new(dim) : throw(DomainError(dim, "dimension must be nonnegative"))
	end
end
ZeroSet(dim) = ZeroSet{DefaultFloat}(dim)

function project!(x::AbstractVector{T}, ::ZeroSet{T}) where{T}
	x .= zero(T)
	return nothing
end

function in_dual(x::AbstractVector{T}, ::ZeroSet{T}, tol::T) where{T}
	return true
end

function in_pol_recc(x::AbstractVector{T}, ::ZeroSet{T}, tol::T) where{T}
	!any( x-> (abs(x) > tol), x)
end

function allocate_memory!(cone::AbstractConvexSet{T}) where {T}
  return nothing
end

# ----------------------------------------------------
# Nonnegative orthant
# ----------------------------------------------------
"""
    Nonnegatives(dim)

Creates the nonnegative orthant ``\\{ x \\in \\mathbb{R}^{dim} : x \\ge 0 \\}``  of dimension `dim`.
"""
struct Nonnegatives{T} <: AbstractConvexCone{T}
	dim::Int
	constr_type::BitArray{1} #store rows of constraints that are loose (+1)
	function Nonnegatives{T}(dim::Int) where {T}
		dim >= 0 ? new(dim, falses(dim)) : throw(DomainError(dim, "dimension must be nonnegative"))
	end
end
Nonnegatives(dim) = Nonnegatives{DefaultFloat}(dim)

"Classify the inequality constraints as loose, if b very large, as s = -Ax + b >= 0."
function classify_constraints!(constr_type::BitArray, b::AbstractVector{T}, COSMO_INFTY::Real, MIN_SCALING::Real) where {T}
	for i = 1:length(b)
		if b[i] > COSMO_INFTY * MIN_SCALING
			constr_type[i] = true
		end
	end
	return nothing
end

function project!(x::AbstractVector{T}, C::Nonnegatives{T}) where{T}
	x .= max.(x, zero(T))
	return nothing
end

function in_dual(x::AbstractVector{T}, ::Nonnegatives{T}, tol::T) where{T}
	return !any( x-> (x < -tol), x)
end

function in_pol_recc(x::AbstractVector{T}, ::Nonnegatives{T}, tol::T) where{T}
	return !any( x-> (x > tol), x)
end

# ----------------------------------------------------
# Second Order Cone
# ----------------------------------------------------
"""
    SecondOrderCone(dim)

Creates the second-order cone (or Lorenz cone) ``\\{ (t,x) \\in \\mathrm{R}^{dim} : || x ||_2  \\leq t \\}``.
"""
struct SecondOrderCone{T} <: AbstractConvexCone{T}
	dim::Int
	function SecondOrderCone{T}(dim::Int) where {T}
		dim >= 0 ? new(dim) : throw(DomainError(dim, "dimension must be nonnegative"))
	end
end
SecondOrderCone(dim) = SecondOrderCone{DefaultFloat}(dim)

function project!(x::AbstractVector{T}, ::SecondOrderCone{T}) where{T}
	t = x[1]
	xt = view(x, 2:length(x))
	norm_x = norm(xt, 2)
	if norm_x <= t
		nothing
	elseif norm_x <= -t
		x[:] .= zero(T)
	else
		x[1] = (norm_x + t) / T(2)
		#x(2:end) assigned via view
		@. xt = (norm_x + t) / (T(2) * norm_x) * xt
	end
	return nothing
end

function in_dual(x::AbstractVector{T}, ::SecondOrderCone{T}, tol::T) where{T}
	@views norm(x[2:end]) <= (tol + x[1]) #self dual
end

function in_pol_recc(x::AbstractVector{T}, ::SecondOrderCone, tol::T) where{T}
	@views norm(x[2:end]) <= (tol - x[1]) #self dual
end

# ----------------------------------------------------
# Positive Semidefinite Cone
# ----------------------------------------------------

#a type to maintain internal workspace data for the BLAS syevr function
mutable struct PsdBlasWorkspace{T}
    m::Base.RefValue{BLAS.BlasInt}
    w::Vector{T}
    Z::Matrix{T}
    isuppz::Vector{BLAS.BlasInt}
    work::Vector{T}
    lwork::BLAS.BlasInt
    iwork::Vector{BLAS.BlasInt}
    liwork::BLAS.BlasInt
    info::Base.RefValue{BLAS.BlasInt}

    function PsdBlasWorkspace{T}(n::Integer) where{T}

        BlasInt = BLAS.BlasInt

        #workspace data for BLAS
        m      = Ref{BlasInt}()
        w      = Vector{T}(undef,n)
        Z      = Matrix{T}(undef,n,n)
        isuppz = Vector{BlasInt}(undef, 2*n)
        work   = Vector{T}(undef, 1)
        lwork  = BlasInt(-1)
        iwork  = Vector{BlasInt}(undef, 1)
        liwork = BlasInt(-1)
        info   = Ref{BlasInt}()

        new(m,w,Z,isuppz,work,lwork,iwork,liwork,info)
    end
end

for (syevr, elty) in
    ((DSYEVR_,:Float64),
     (SSYEVR_,:Float32))
   @eval begin
        function _syevr!(A::AbstractMatrix{$elty}, ws::PsdBlasWorkspace{$elty})

            #Float64 only support for now since we call dsyevr_ directly
            n       = size(A,1)
            ldz     = n
            lda     = stride(A,2)

                ccall($syevr, Cvoid,
                (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BLAS.BlasInt},
                Ptr{$elty}, Ref{BLAS.BlasInt}, Ref{$elty}, Ref{$elty},
                Ref{BLAS.BlasInt}, Ref{BLAS.BlasInt}, Ref{$elty}, Ptr{BLAS.BlasInt},
                Ptr{$elty}, Ptr{$elty}, Ref{BLAS.BlasInt}, Ptr{BLAS.BlasInt},
                Ptr{$elty}, Ref{BLAS.BlasInt}, Ptr{BLAS.BlasInt}, Ref{BLAS.BlasInt},
                Ptr{BLAS.BlasInt}),
                'V', 'A', 'U', n,
                A, max(1,lda), 0.0, 0.0,
                0, 0, -1.0,
                ws.m, ws.w, ws.Z, ldz, ws.isuppz,
                ws.work, ws.lwork, ws.iwork, ws.liwork,
                ws.info)
                LAPACK.chklapackerror(ws.info[])
        end
    end #@eval
end #for

function _project!(X::AbstractMatrix, ws::PsdBlasWorkspace{T}) where{T}

    #computes the upper triangular part of the projection of X onto the PSD cone

     #allocate additional workspace arrays if the ws
     #work and iwork have not yet been sized
     if ws.lwork == -1
         _syevr!(X,ws)
         ws.lwork = BLAS.BlasInt(real(ws.work[1]))
         resize!(ws.work, ws.lwork)
         ws.liwork = ws.iwork[1]
         resize!(ws.iwork, ws.liwork)
     end

	 # below LAPACK function does the following: w,Z  = eigen!(Symmetric(X))
	 	_syevr!(X, ws)
		# compute upper triangle of: X .= Z*Diagonal(max.(w, 0.0))*Z'
		rank_k_update!(X, ws)
end

function rank_k_update!(X::AbstractMatrix, ws::COSMO.PsdBlasWorkspace{T}) where {T}
  n = size(X, 1)
  @. X = zero(T)
  nnz_λ = 0
  for j = 1:length(ws.w)
    λ = ws.w[j]
    if λ > 0
      nnz_λ += 1
      @inbounds for i = 1:n
        ws.Z[i, j] = ws.Z[i, j] * sqrt(λ)
      end
    end
  end

  if nnz_λ > 0
    V = uview(ws.Z, :, (n - nnz_λ + 1):n)
    BLAS.syrk!('U', 'N', one(T), V, one(T), X)
  end
  return nothing
end

"""
    PsdCone(dim)

Creates the cone of symmetric positive semidefinite matrices ``\\mathcal{S}_+^{dim}``. The entries of the matrix `X` are stored column-by-column in the vector `x` of dimension `dim`.
Accordingly  ``X \\in \\mathbb{S}_+ \\Rightarrow x \\in \\mathcal{S}_+^{dim}``, where ``X = \\text{mat}(x)``.
"""
struct PsdCone{T} <: AbstractConvexCone{T}
	dim::Int
	sqrt_dim::Int
  work::PsdBlasWorkspace{T}
  tree_ind::Int  # tree number that this cone belongs to
  clique_ind::Int
	function PsdCone{T}(dim::Int, tree_ind::Int, clique_ind::Int) where{T}
		dim >= 0       || throw(DomainError(dim, "dimension must be nonnegative"))
		iroot = isqrt(dim)
		iroot^2 == dim || throw(DomainError(dim, "dimension must be a square"))
		new(dim, iroot, PsdBlasWorkspace{T}(iroot), tree_ind, clique_ind)
	end
end
PsdCone(dim) = PsdCone{DefaultFloat}(dim)
PsdCone{T}(dim::Int) where{T} = PsdCone{T}(dim, 0, 0)

struct DensePsdCone{T} <: AbstractConvexCone{T}
  dim::Int
  sqrt_dim::Int
  work::PsdBlasWorkspace{T}
  function DensePsdCone{T}(dim::Int) where{T}
    dim >= 0       || throw(DomainError(dim, "dimension must be nonnegative"))
    iroot = isqrt(dim)
    iroot^2 == dim || throw(DomainError(dim, "dimension must be a square"))
    new(dim, iroot, PsdBlasWorkspace{T}(iroot))
  end
end
DensePsdCone(dim) = DensePsdCone{DefaultFloat}(dim)

function project!(x::AbstractVector{T}, cone::Union{PsdCone{T}, DensePsdCone{T}}) where{T}
	n = cone.sqrt_dim

    # handle 1D case
    if length(x) == 1
        x .= max(x[1], zero(T))
    else
        # symmetrized square view of x
        X = reshape(x, n, n)
        symmetrize_upper!(X)
        _project!(X, cone.work)

        #fill in the lower triangular part
        for j=1:n, i=1:(j-1)
            X[j,i] = X[i,j]
        end
    end
    return nothing
end

# Notice that this is an in-place version that uses x as workspace
function in_dual!(x::AbstractVector{T}, cone::Union{PsdCone{T}, DensePsdCone{T}}, tol::T) where{T}
	n = cone.sqrt_dim
	X = reshape(x, n, n)
  return COSMO.is_pos_def!(X, tol)
end
in_dual(x::AbstractVector{T}, cone::Union{PsdCone{T}, DensePsdCone{T}}, tol::T) where{T} = in_dual!(copy(x), cone, tol)

function in_pol_recc!(x::AbstractVector{T}, cone::Union{PsdCone{T}, DensePsdCone{T}}, tol::T) where{T}
	n = cone.sqrt_dim
	X = reshape(x, n, n)
	return COSMO.is_neg_def!(X, tol)
end
in_pol_recc(x::AbstractVector{T}, cone::Union{PsdCone{T}, DensePsdCone{T}}, tol::T) where{T} = in_pol_recc!(copy(x), cone, tol)

# ----------------------------------------------------
# Positive Semidefinite Cone (Triangle)
# ----------------------------------------------------
# Psd cone given by upper-triangular entries of matrix
"""
    PsdConeTriangle(dim)

Creates the cone of symmetric positive semidefinite matrices. The entries of the upper-triangular part of matrix `X` are stored in the vector `x` of dimension `dim`.
A ``r \\times r`` matrix has ``r(r+1)/2`` upper triangular elements and results in a vector of ``\\mathrm{dim} = r(r+1)/2``.

### Examples
The matrix
```math
\\begin{bmatrix} x_1 & x_2 & x_4\\\\ x_2 & x_3 & x_5\\\\ x_4 & x_5 & x_6 \\end{bmatrix}
```
is transformed to the vector ``[x_1, x_2, x_3, x_4, x_5, x_6]^\\top `` with corresponding constraint  `PsdConeTriangle(6)`.

"""
mutable struct PsdConeTriangle{T} <: AbstractConvexCone{T}
    dim::Int #dimension of vector
    sqrt_dim::Int # side length of matrix
    X::Array{T,2}
    work::PsdBlasWorkspace{T}
    tree_ind::Int # tree number that this cone belongs to
    clique_ind::Int

    function PsdConeTriangle{T}(dim::Int, tree_ind::Int, clique_ind::Int) where{T}
        dim >= 0       || throw(DomainError(dim, "dimension must be nonnegative"))
        side_dimension = Int(sqrt(0.25 + 2 * dim) - 0.5);
        new(dim, side_dimension, zeros(side_dimension, side_dimension), PsdBlasWorkspace{T}(side_dimension), tree_ind, clique_ind)

    end
end
PsdConeTriangle(dim) = PsdConeTriangle{DefaultFloat}(dim)
PsdConeTriangle{T}(dim::Int) where{T} = PsdConeTriangle{T}(dim, 0, 0)

DecomposableCones{T} = Union{PsdCone{T}, PsdConeTriangle{T}}

mutable struct DensePsdConeTriangle{T} <: AbstractConvexCone{T}
    dim::Int #dimension of vector
    sqrt_dim::Int # side length of matrix
    X::Array{T,2}
    work::PsdBlasWorkspace{T}

    function DensePsdConeTriangle{T}(dim::Int) where{T}
        dim >= 0       || throw(DomainError(dim, "dimension must be nonnegative"))
        side_dimension = Int(sqrt(0.25 + 2 * dim) - 0.5);
        new(dim, side_dimension, zeros(side_dimension, side_dimension), PsdBlasWorkspace{T}(side_dimension))
    end
end
DensePsdConeTriangle(dim) = DensePsdConeTriangle{DefaultFloat}(dim)

function project!(x::AbstractArray, cone::Union{PsdConeTriangle{T}, DensePsdConeTriangle{T}}) where{T}
    # handle 1D case
    if length(x) == 1
        x .= max(x[1],zero(T))
    else
        populate_upper_triangle!(cone.X, x, 1 / sqrt(2))
        _project!(cone.X,cone.work)
        extract_upper_triangle!(cone.X, x, sqrt(2) )
    end
    return nothing
end

# Notice that we are using a (faster) in-place version that modifies the input
function in_dual!(x::AbstractVector{T}, cone::Union{PsdConeTriangle{T}, DensePsdConeTriangle{T}}, tol::T) where{T}
    n = cone.sqrt_dim
    populate_upper_triangle!(cone.X, x, 1 / sqrt(2))
    return COSMO.is_pos_def!(cone.X, tol)
end
in_dual(x::AbstractVector{T}, cone::Union{PsdConeTriangle{T}, DensePsdConeTriangle{T}}, tol::T) where {T} = in_dual!(x, cone, tol)

function in_pol_recc!(x::AbstractVector{T}, cone::Union{PsdConeTriangle{T}, DensePsdConeTriangle{T}}, tol::T) where{T}
    n = cone.sqrt_dim
    populate_upper_triangle!(cone.X, x, 1 / sqrt(2))
    Xs = Symmetric(cone.X)
    return COSMO.is_neg_def!(cone.X, tol)
end
in_pol_recc(x::AbstractVector{T}, cone::Union{PsdConeTriangle{T}, DensePsdConeTriangle{T}}, tol::T) where {T} = in_pol_recc!(x, cone, tol)

function allocate_memory!(cone::Union{PsdConeTriangle{T}, DensePsdConeTriangle{T}}) where {T}
  cone.X = zeros(cone.sqrt_dim, cone.sqrt_dim)
end

function populate_upper_triangle!(A::AbstractMatrix, x::AbstractVector, scaling_factor::Float64)
 	k = 0
  	for j in 1:size(A, 2)
     	for i in 1:j
        	k += 1
        	if i != j
        		A[i, j] = scaling_factor * x[k]
        	else
        		A[i, j] = x[k]
        	end
      	end
  	end
  	nothing
end

function extract_upper_triangle!(A::AbstractMatrix, x::AbstractVector, scaling_factor::Float64)
	k = 0
  	for j in 1:size(A, 2)
     	for i in 1:j
        	k += 1
        	if i != j
        		x[k] = scaling_factor * A[i, j]
        	else
        		x[k] = A[i, j]
        	end
      	end
  	end
	nothing
end

"""
    ExponentialCone(MAX_ITERS = 100, EXP_TOL = 1e-8)

Creates the exponential cone ``\\mathcal{K}_{exp} = \\{(x, y, z) \\mid y \\geq 0 ye^{x/y} ≤ z\\} \\cup \\{ (x,y,z) \\mid   x \\leq 0, y = 0, z \\geq 0 \\}``
"""
struct ExponentialCone{T} <: AbstractConvexCone{T}
  dim::Int
  v0::Vector{T}
  MAX_ITER::Int
  EXP_TOL::T

  function ExponentialCone{T}(dim = 3, MAX_ITERS = 100, EXP_TOL = T(1e-8)) where{T}
    new(3, zeros(T, 3), MAX_ITERS, EXP_TOL)
  end
end
ExponentialCone(args...) = ExponentialCone{DefaultFloat}(args...)

function project!(v::AbstractVector{T}, cone::ExponentialCone{T}) where{T}

  # Check the four different cases
  # 1. v in K_exp => v = v
  in_cone(v, cone, zero(T)) && return nothing

  # 2. -v in K_exp^* => v = 0
  if in_dual(-v, cone, zero(T))
    v .= zero(T)
    return nothing
  end

  # 3. x < 0 and y < 0 => v = (x, 0, max(z, 0))
  if v[1] < 0 && v[2] < 0
    v[2] = T(0)
    v[3] = max(v[3], zero(T))
    return nothing
  end

  # 4. Otherwise solve the following minimisation problem
  # min_w (1/2) ||v - v0||_2^2
  # s.t.  y * exp(x/y) == z
  #       y > 0
  project_exp!(v, cone)
end

# This is a modified version of the projection code used in SCS
# https://github.com/cvxgrp/scs/blob/master/src/cones.c
# We are solving the dual problem g(λ) via a bisection method
function project_exp!(v::AbstractVector{T}, cone::ExponentialCone{T}) where {T <: AbstractFloat}
  # save input vector and use v as working variable
  @. cone.v0 = v
  l, u = get_bisection_bounds(v, cone.v0, cone.EXP_TOL)

  for k = 1:cone.MAX_ITER
    λ = (u + l) / 2
    g = grad_dual!(λ, v, cone.v0, cone.EXP_TOL)
    g > 0 ? (l = λ) : (u = λ)
    u - l < cone.EXP_TOL && break
  end
end

function get_bisection_bounds(v::AbstractVector{T}, v0::Vector{T}, tol::T) where {T <: AbstractFloat}
  l = zero(T)
  λ = T(0.125)
  g = grad_dual!(λ, v, v0, tol)
  while g > 0
    l = λ
    λ *= 2
    g = grad_dual!(λ, v, v0, tol)
  end
  u = λ
  return l, u
end

function grad_dual!(λ::T, v::AbstractVector{T}, v0::Vector{T}, tol::T) where {T <: AbstractFloat}
  find_minimizers!(λ, v, v0, tol)
  v[2] == 0 ? (g = v[1]) : (g = v[1] + v[2] * log(v[2] / v[3]))
  return g
end

function find_minimizers!(λ::T, v::AbstractVector{T}, v0::Vector{T}, tol::T) where {T <: AbstractFloat}
  v[3] = find_min_t(λ, v0[2], v0[3], tol)
  # s* = (t - t0) * t / λ
  v[2] = (1 / λ) * (v[3] - v0[3]) * v[3]
  # r* = r0 - λ
  v[1] = v0[1] - λ
end

# use Newton method to find minimizer t* for given λ, i.e. find the zero of
# f(t) = t * (t - t0) / lambda - s0 + λ * log( t - t0 / λ) + λ
# Define Δt = t - t0
function find_min_t(λ::T, s0::T, t0::T, tol::T) where {T <: AbstractFloat}
  Δt = max(-t0, tol)
  for k = 1:150
    f = Δt * (Δt + t0) / λ^2 - s0 / λ + log(Δt / λ) + one(T)
    grad_f = (2 * Δt + t0) / λ^2 + one(T) / Δt
    Δt = Δt - f / grad_f

    if (Δt <= -t0)
      Δt = -t0
      break
    elseif (Δt <= 0)
      Δt = zero(T)
      break
    elseif abs(f) < tol
      break
    end
  end
  return Δt + t0
end

function in_cone(v::AbstractVector{T}, cone::ExponentialCone{T}, tol::T) where {T <: AbstractFloat}
  x = v[1]
  y = v[2]
  z = v[3]
  return (y > 0 && y * exp(x/y) <= z + tol) || (x <= tol &&  y == 0. && z >= -tol )
end
# Kexp^* = { (x,y,z) | x < 0, -xe^(y/x) <= e^1 z } cup { (0,y,z) | y >= 0,z >= 0 }
function in_dual(v::AbstractVector{T}, cone::ExponentialCone{T}, tol::T) where {T <: AbstractFloat}
  x = v[1]
  y = v[2]
  z = v[3]
  return (x < 0 && -x * exp(y / x) - exp(1) *  z <= tol) || (abs(x) <= tol && y >= -tol && z >= -tol)
end

function in_pol_recc(v::AbstractVector{T},cone::ExponentialCone{T}, tol::T) where {T <: AbstractFloat}
  return in_dual(-v, cone, tol)
end

"""
    PowerCone(alpha::Float64, MAX_ITERS::Int = 20, POW_TOL = 1e-8)

Creates the 3-d power cone ``\\mathcal{K}_{pow} = \\{(x, y, z) \\mid x^\\alpha y^{(1-\\alpha)} \\geq  \\|z\\|, x \\geq 0, y \\geq 0 \\}`` with ``0 < \\alpha < 1``
"""
struct PowerCone{T} <: AbstractConvexCone{T}
  dim::Int
  α::T
  MAX_ITER::Int
  POW_TOL::T

  function PowerCone{T}(alpha::Real, MAX_ITERS::Int = 20, POW_TOL::Real = 1e-8) where{T <: AbstractFloat}
    (alpha <= 0 || alpha >= one(T)) && throw(DomainError("The exponent α of the power cone has to be in (0, 1)."))
    new(3, alpha, MAX_ITERS, POW_TOL)
  end
end
PowerCone(alpha::T, args...) where {T <: AbstractFloat}= PowerCone{T}(alpha, args...)

# Allow precision conversion of PowerCones
function convert(dest_type::Type{COSMO.PowerCone{Ta}}, src::COSMO.PowerCone{Tb}) where {Ta <: AbstractFloat, Tb <: AbstractFloat}
	dest_type(Ta(src.α), src.MAX_ITER, Ta(src.POW_TOL))
end

# The projection onto the power cone is described in
# Hien - Differential properties of Euclidean projections onto power cone (2015)
function project!(v::AbstractVector{T}, cone::PowerCone{T}) where{T}
  # Check the special cases first
  # 1. v in K_pow => v = v
  in_cone(v, cone, zero(T)) && return nothing

  # 2. -v in K_pow^* => v .= 0
  if in_dual(-v, cone, zero(T))
    v .= zero(T)
    return nothing
  end

  # 3. v not in K_pow and -v not in K_pow^* and z == 0 => x = max(x, 0), y = max(y, 0)
  if abs(v[3]) <= cone.POW_TOL
    v[1] = max(v[1], zero(T))
    v[2] = max(v[2], zero(T))
    return nothing
  end

  # 4. Otherwise solve the following problem
  # find  r
  # s.t.  sigma(x0, y0, z0, r) == 0
  # and   0 < r < ||z0||
  #
  # x = 0.5 * (x0 + sqrt(x0^2 + 4α*r*(||z0||-r)))
  # y = 0.5 * (y0 + sqrt(y0^2 + 4(1-α)*r*(||z0||-r)))
  # z = z0 * r / ||z0||
  project_pow!(v, cone)
end

 # find the zero of above condition for r by applying Newton's method
  function project_pow!(v::AbstractVector{T}, cone::PowerCone{T}) where {T <: AbstractFloat}
  x0 = v[1]
  y0 = v[2]
  z0 = v[3]
  r = abs(z0) / T(2)
  ϕx = zero(T)
  ϕy = zero(T)
  # compute a zero of phi(v, r, α)
  for k = 1:cone.MAX_ITER
    ϕx = ϕc(x0, z0, r, cone.α)
    ϕy = ϕc(y0, z0, r, 1 - cone.α)
    phi = ϕ(ϕx, ϕy, r, cone.α)
    abs(phi) < cone.POW_TOL && break

    dϕx_dr = dϕc_dr(ϕx, x0, z0, r, cone.α)
    dϕy_dr = dϕc_dr(ϕy, y0, z0, r, 1 - cone.α)
    dphi_dr = dϕ_dr(ϕx, ϕy, dϕx_dr, dϕy_dr, r, cone.α)
    r = r - phi / dphi_dr

    # ensure 0 < r < abs(z0)
    r = min(max(r, 0), abs(z0))
  end

  # given a solution r, update the vector components (x, y, z)
  v[1] = ϕx
  v[2] = ϕy
  v[3] = z0 * r / abs(z0)
  return nothing
end

function ϕc(x0::T, z0::T, r::T, α::T) where{T <: Real}
  return max(T(0.5) * (x0 + sqrt(x0^2 + T(4) * α * r * (abs(z0) - r))), T(1e-10))
end

function dϕc_dr(ϕx::T, x0::T, z0::T, r::T, α::T) where{T <: Real}
  return α / (2 * ϕx - x0) * (abs(z0) - 2 * r)
end

function ϕ(ϕx::T, ϕy::T, r::T, α::T) where{T <: Real}
  return ϕx^α * ϕy^(1 - α) - r
end

# dϕ / dr = Π fi^αi * (Σ αi fi' / fi) - 1
function dϕ_dr(ϕx::T, ϕy::T, ϕx_dr::T, ϕy_dr::T, r::T, α::T) where{T <: Real}
  return ϕx^α * ϕy^(1-α) * (α * ϕx_dr / ϕx + (1 - α) * ϕy_dr / ϕy) - 1
end

function in_cone(v::AbstractVector{T}, cone::PowerCone{T}, tol::T) where{T}
  x = v[1]
  y = v[2]
  z = v[3]
  α = cone.α
  return x >= 0 && y >= 0 && x^α * y^(1 - α) >= abs(z) - tol
end

# Kpow^* = { (s, t, w) | (s / α)^α * (t/(1-α))^(1-α) >= abs(w), s >= 0, t >= 0 }
function in_dual(v::AbstractVector{T}, cone::PowerCone{T}, tol::T) where{T}
  s = v[1]
  t = v[2]
  w = v[3]
  α = cone.α
  return s >= -tol && t >= -tol && s^α * t^(1-α) >= abs(w) * α^α * (1-α)^(1-α) - tol
end

function in_pol_recc(v::AbstractVector{T},cone::PowerCone{T}, tol::T) where{T}
  return in_dual(-v, cone, tol)
end

"""
    DualExponentialCone(MAX_ITERS::Int = 100, EXP_TOL = 1e-8)

Creates the dual exponential cone ``\\mathcal{K}^*_{exp} = \\{(x, y, z) \\mid x < 0,  -xe^{y/x} \\leq e^1 z \\} \\cup \\{ (0,y,z) \\mid   y \\geq 0, z \\geq 0 \\}``
"""
struct DualExponentialCone{T} <: AbstractConvexCone{T}
  dim::Int
  v0::Vector{T}
  primal_cone::ExponentialCone{T}

  function DualExponentialCone{T}(dim::Int = 3, MAX_ITERS::Int = 100, EXP_TOL = T(1e-8)) where{T}
    new(3, zeros(T, 3), ExponentialCone{T}(dim, MAX_ITERS, EXP_TOL))
  end
end
DualExponentialCone(args...) = DualExponentialCone{DefaultFloat}(args...)

"""
    DualPowerCone(alpha::Float64, MAX_ITERS::Int = 20, POW_TOL = 1e-8)

Creates the 3-d dual power cone ``\\mathcal{K}^*_{pow} = \\{(u, v, w) \\mid \\left( \\frac{u}{\\alpha}\\right)^\\alpha \\left( \\frac{v}{1-\\alpha}\\right)^{(1-\\alpha)} \\geq  \\|w\\|, u \\geq 0, v \\geq 0 \\}`` with ``0 < \\alpha < 1``
"""
struct DualPowerCone{T} <: AbstractConvexCone{T}
  dim::Int
  v0::Vector{T}
  primal_cone::PowerCone{T}

  function DualPowerCone{T}(alpha::Real, MAX_ITERS::Int = 20, POW_TOL = T(1e-8)) where{T}
    (alpha <= 0 || alpha >= 1) && throw(DomainError, "The exponent α of the dual power cone has to be in (0, 1).")
    new(3, zeros(T,3), PowerCone{T}(alpha, MAX_ITERS, POW_TOL))
  end
end
DualPowerCone(alpha::T, args...) where {T <: AbstractFloat}= DualPowerCone{T}(alpha, args...)

DualCones{T} = Union{DualExponentialCone{T}, DualPowerCone{T}}
in_cone(v::AbstractVector{T}, cone::DualCones{T}, tol::T) where {T <: AbstractFloat} = in_dual(v, cone.primal_cone, tol)
in_dual(v::AbstractVector{T}, cone::DualCones{T}, tol::T) where {T <: AbstractFloat} = in_cone(v, cone.primal_cone, tol)
in_pol_recc(v::AbstractVector{T}, cone::DualCones{T}, tol::T) where {T <: AbstractFloat} = in_dual(-v, cone, tol)

# Project dual cones by using Moreau decomposition: Proj^*(v) = v + Proj(-v)
function project!(v::AbstractVector{T}, cone::DualCones{T}) where{T <: AbstractFloat}
  @. cone.v0 = v
  @. v *= -one(T)
  project!(v, cone.primal_cone)
  @. v += cone.v0
end

# Union for all types where the user has to provide extra information to create the cone
ArgumentCones = Union{PowerCone, DualPowerCone}

# ----------------------------------------------------
# Box
# ----------------------------------------------------
"""
    Box(l, u)

Creates a box or intervall with lower boundary vector ``l \\in  \\mathbb{R}^m \\cup \\{-\\infty\\}^m`` and upper boundary vector``u \\in \\mathbb{R}^m\\cup \\{+\\infty\\}^m``.
"""
struct Box{T} <: AbstractConvexSet{T}
	dim::Int
	constr_type::Vector{Int} #store type of constraint {-1: loose, 0: inequality, 1: equality}
	l::Vector{T}
	u::Vector{T}
	function Box{T}(dim::Int) where {T <: AbstractFloat}
		dim >= 0 || throw(DomainError(dim, "dimension must be nonnegative"))
		l = fill!(Vector{T}(undef, dim), -Inf)
		u = fill!(Vector{T}(undef, dim), +Inf)
		return new(dim, zeros(Int, dim), l, u)
	end
	function Box{T}(l::Vector{T}, u::Vector{T}) where {T <: AbstractFloat}
		length(l) == length(u) || throw(DimensionMismatch("bounds must be same length"))
		#enforce consistent bounds
        _box_check_bounds(l,u)
		return new(length(l), zeros(Int, length(l)), l, u)
	end
end
Box(dim) = Box{DefaultFloat}(dim)
Box(l::AbstractArray{T}, u::AbstractArray{T}) where {T <: AbstractFloat} = Box{T}(l, u)

function _box_check_bounds(l,u)
    for i in eachindex(l)
        l[i] > u[i] && error("Box set: inconsistent lower/upper bounds specified at index i = ", i, ": l[i] = ",l[i],", u[i] = ",u[i])
    end
end

"Classify the type of constraint of the box into loose inequality (-1), inequality (0) or equality (+1)."
function classify_box_constraints!(box::COSMO.Box{T}, COSMO_INFTY::Real, MIN_SCALING::Real, RHO_TOL::Real) where{T}
	@inbounds for i = 1:length(box.l)
		if box.l[i] < (-COSMO_INFTY * MIN_SCALING) && box.u[i] > (COSMO_INFTY * MIN_SCALING)
			box.constr_type[i] = -1
		elseif (box.u[i] - box.l[i]) < RHO_TOL
			box.constr_type[i] = 1
		else
			box.constr_type[i] = 0
		end
	end
	return nothing
end

function project!(x::AbstractVector{T}, box::Box{T}) where{T}
	@. x = clip(x, box.l, box.u)
	return nothing
end

function support_function(x::AbstractVector{T}, B::Box{T}, tol::T) where{T}
    s = 0.
    for i in eachindex(x)
        s+= ( abs(x[i]) > tol && x[i] > 0) ? x[i]*B.u[i] : x[i]*B.l[i]
    end
    return s
end
support_function!(x::AbstractVector{T}, B::Box{T}, tol::T) where{T} = support_function(x, B, tol)

function in_pol_recc(x::AbstractVector{T}, B::Box{T}, tol::T) where{T}
    !any(XU -> (XU[2] == Inf && XU[1] > tol), zip(x,B.u)) && !any(XL -> (XL[2] == -Inf && XL[1] < -tol), zip(x,B.l))
end

function scale!(box::Box{T}, e::AbstractVector{T}) where{T}
	@. box.l = box.l * e
	@. box.u = box.u * e
	return nothing
end

function Base.deepcopy(box::Box{T}) where {T}
  Box{T}(deepcopy(box.l), deepcopy(box.u))
end

# ----------------------------------------------------
# Composite Set
# ----------------------------------------------------

#struct definition is provided in projections.jl, since it
#must be available to SplitVector, which in turn must be
#available for most of the methods here.

CompositeConvexSet(args...) = CompositeConvexSet{DefaultFloat}(args...)

function project!(x::SplitVector{T}, C::CompositeConvexSet{T}) where{T}
    for i = 1:length(C.sets)
        project!(x.views[i],C.sets[i])
    end
	#foreach(xC -> project!(xC[1], xC[2]), zip(x.views, C.sets))
	return nothing
end

function support_function!(x::SplitVector{T}, C::CompositeConvexSet{T}, tol::T) where{T}
	sum(xC -> support_function!(xC[1], xC[2], tol), zip(x.views, C.sets))
end

function in_pol_recc(x::SplitVector{T}, C::CompositeConvexSet{T}, tol::T) where{T}
	all(xC -> in_pol_recc(xC[1], xC[2], tol), zip(x.views, C.sets))
end

function in_pol_recc!(x::SplitVector{T}, C::CompositeConvexSet{T}, tol::T) where{T}
  all(xC -> in_pol_recc!(xC[1], xC[2], tol), zip(x.views, C.sets))
end

function scale!(C::CompositeConvexSet{T}, e::SplitVector{T}) where{T}
	for i = eachindex(C.sets)
		scale!(C.sets[i], e.views[i])
	end
end

function rectify_scaling!(E::SplitVector{T},
	work::SplitVector{T},
	C::CompositeConvexSet{T}) where {T}
	any_changed = false
	for i = eachindex(C.sets)
		any_changed |= rectify_scaling!(E.views[i], work.views[i], C.sets[i])
	end
	return any_changed
end

#-------------------------
# general AbstractConvexCone operations
#-------------------------

# sup_{z in K_tilde_b = {-K} x {b} } <z,δy> = { <y,b> ,if y in Ktilde_polar
#                                                 +∞   ,else}

function support_function(y::SplitView{T}, cone::AbstractConvexCone{T}, tol::T) where{T}
  in_dual(-y, cone, tol) ? 0. : Inf;
end

# An in-place method that is faster, but uses the input variable y as workspace
function support_function!(y::SplitView{T}, cone::AbstractConvexCone{T}, tol::T) where{T}
  @. y *= - one(T)
  in_dual!(y, cone, tol) ? 0. : Inf;
end

# Notice: for every convex set apart from PsdCone and PsdConeTriangle use the normal non-modifying function
# for PsdCone and PsdConeTriangle we are using (faster) in-place functions.
function in_dual!(x::AbstractVector{T}, cone::AbstractConvexSet{T}, tol::T) where{T}
  return in_dual(x, cone, tol)
end

function in_pol_recc!(x::AbstractVector{T}, cone::AbstractConvexSet{T}, tol::T) where{T}
  return in_pol_recc(x, cone, tol)
end

function scale!(cone::AbstractConvexSet{T}, ::AbstractVector{T}) where{T}
  return nothing
end

# The fall-back method. Notice that this is conservative as not all custom cones that users might define have to be scalar-scaled.
rectify_scaling!(E, work, set::AbstractConvexSet{T}) where {T} = rectify_scalar_scaling!(E, work)

function rectify_scaling!(E, work, set::Union{SecondOrderCone{T}, PsdCone{T}, DensePsdCone{T}, PsdConeTriangle{T}, DensePsdConeTriangle{T}, PowerCone{T}, DualPowerCone{T}, ExponentialCone{T}, DualExponentialCone{T}}) where{T}
  return rectify_scalar_scaling!(E, work)
end
rectify_scaling!(E, work, set::Union{ZeroSet{<:Real}, Nonnegatives{<:Real}, Box{<:Real}}) = false

#-------------------------
# generic set operations
#-------------------------
# function Base.showarg(io::IO, C::AbstractConvexSet{T}, toplevel) where{T}
#    print(io, typeof(C), " in dimension '", A.dim, "'")
# end

eltype(::AbstractConvexSet{T}) where{T} = T
num_subsets(C::AbstractConvexSet{T}) where{T}  = 1
num_subsets(C::CompositeConvexSet{T}) where{T} = length(C.sets)

function get_subset(C::AbstractConvexSet, idx::Int)
	idx == 1 || throw(DimensionMismatch("Input only has 1 subset (itself)"))
	return C
end
get_subset(C::CompositeConvexSet, idx::Int) = C.sets[idx]

function rectify_scalar_scaling!(E, work)
	tmp = mean(E)
	work .= tmp ./ E
	return true
end

# computes the row indices of A,b for each convex set
function get_set_indices(sets::Array{COSMO.AbstractConvexSet, 1})
	sidx = 0
	indices = Array{UnitRange{Int}, 1}(undef, length(sets))
	for i = eachindex(sets)
		indices[i] = (sidx + 1) : (sidx + sets[i].dim)
		sidx += sets[i].dim
	end
	return indices
end


text:
using Plots

function bracket_minimum(f, x=0; s=1e-2, k=2.0)
    a, ya = x, f(x)
    b, yb = a + s, f(a + s)
    if yb > ya
        a, b = b, a
        ya, yb = yb, ya
        s = -s
    end
    while true
        c, yc = b + s, f(b + s)
        if yc > yb
            return a < c ? (a,b, c) : (c, b, a)
        end
        a, ya, b, yb = b, yb, c, yc
        s *= k
    end
end

function quadratic_fit_search(f, a, b, c, n)
    ya, yb, yc = f(a), f(b), f(c)
    for i in 1:n-3
        x = 0.5*(ya*(b^2-c^2)+yb*(c^2-a^2)+yc*(a^2-b^2)) /
            (ya*(b-c) +yb*(c-a) +yc*(a-b))
        yx = f(x)
        if x > b
            if yx > yb
                c, yc = x, yx
            else
                a, ya, b, yb = b, yb, x, yx
            end
        elseif x < b
            if yx > yb
                a, ya = x, yx
            else
                c, yc, b, yb = b, yb, x, yx
            end
        end
    end
    return (a, b, c)
end

#definicao da funcao
f = x-> -(0.4/((x^2+1))^(1/2))+((x^2+1)^(1/2))*(1-0.4/(x^2+1))-x    #multiplicada por -1 para encontrar o minimo
g = x-> -(-(0.4/((x^2+1))^(1/2))+((x^2+1)^(1/2))*(1-0.4/(x^2+1))-x) #funcao original

#definiciao do ponto inicial
x0 = 0
#busca do intervalo
a,b,c = bracket_minimum(f,x0)

#numero de iteracoes
φ = MathConstants.φ
n=(c-a)/(10^-3*log(φ))

am,bm,cm = quadratic_fit_search(f,a,b,c,n)
#(am,f(am))
#am,bm,cm

println("x* = ", cm)

plot(g,a-1.,c+1.,lw=2,draw_arrow="true",label="Function")
scatter!([a,b,c],[g(a),g(b),g(c)],label="interval")
scatter!([cm],[g(cm)],label="x*")


text:
using LinearAlgebra
using Random
using Distributions

# plot features require Gadfly, etc
using ApproxManifoldProducts

# more packages required
using TransformUtils

const LinAlg = LinearAlgebra
const KDE = KernelDensityEstimate
const TU = TransformUtils
const AMP = ApproxManifoldProducts

# temporarily required for experimental versions of on-manifold products
KDE.setForceEvalDirect!(true)

## create two densities

# two densities on a cylinder
p = manikde!(randn(2,100), (:Euclid, :Circular) )

pts2a = 3.0*randn(1,100).+5.0
pts2b = TransformUtils.wrapRad.(0.5*randn(1,100).+pi)
q = manikde!([pts2a;pts2b], (:Euclid, :Circular) )

# approximate the product between hybrid manifold densities
pq = manifoldProduct([p;q], (:Euclid, :Circular), Niter=2)
# getBW(pq)[:,1]

## plotting with Makie

using AbstractPlotting, FileIO
using Makie

N = 200
x = LinRange(-7.5, 7.5, N)
y = LinRange(-pi, pi, N)
zp = zeros(N,N)
zq = zeros(N,N)
zpq = zeros(N,N)

pt = zeros(2,1)
for i in 1:length(x), j in 1:length(y)
  pt[1,1] = x[i]
  pt[2,1] = y[j]
  zp[i,j] = p(pt,false, 1e-3, (+, addtheta), (-,difftheta))[1]
  zq[i,j] = q(pt,false, 1e-3, (+, addtheta), (-,difftheta))[1]
  zpq[i,j] = pq(pt,false, 1e-3, (+, addtheta), (-,difftheta))[1]
end

## Compound plot

scene = hbox(
        # contour(x, y, z, levels = 20, linewidth =3),
        contour(x, y, zpq, levels = 0, linewidth = 0, fillrange = true),
        contour(x, y, zq, levels = 0, linewidth = 0, fillrange = true),
        contour(x, y, zp, levels = 0, linewidth = 0, fillrange = true),
)

## simpler plot

scene = contour(x, y, zpq, levels = 0, linewidth = 0, fillrange = true)

# save("/tmp/test.png", scene)

## better wrap around

pts1a = 6.0*rand(1,100).-3.0
pts1b = TransformUtils.wrapRad.(0.3*randn(1,100).-0.6*pi)
p = manikde!([pts1a;pts1b], (:Euclid, :Circular) )

pts2a = 6.0*rand(1,100).-3.0
pts2b = TransformUtils.wrapRad.(0.3*randn(1,100).+0.6*pi)
q = manikde!([pts2a;pts2b], (:Euclid, :Circular) )

# approximate the product between hybrid manifold densities
pq = manifoldProduct([p;q], (:Euclid, :Circular), Niter=3)

N = 200
x = LinRange(-7.5, 7.5, N)
y = LinRange(-pi, pi, N)
zp = zeros(N,N)
zq = zeros(N,N)
zpq = zeros(N,N)

pt = zeros(2,1)
for i in 1:length(x), j in 1:length(y)
  pt[1,1] = x[i]
  pt[2,1] = y[j]
  zp[i,j] = p(pt,false, 1e-3, (+, addtheta), (-,difftheta))[1]
  zq[i,j] = q(pt,false, 1e-3, (+, addtheta), (-,difftheta))[1]
  zpq[i,j] = pq(pt,false, 1e-3, (+, addtheta), (-,difftheta))[1]
end

scene = hbox(
  # contour(x, y, z, levels = 20, linewidth =3),
  contour(x, y, zpq, levels = 0, linewidth = 0, fillrange = true),
  contour(x, y, zq, levels = 0, linewidth = 0, fillrange = true),
  contour(x, y, zp, levels = 0, linewidth = 0, fillrange = true),
)

## Save scene to file

# using AbstractPlotting, CairoMakie, FileIO
# AbstractPlotting.current_backend[] = CairoMakie.CairoBackend("/tmp/test.svg")
# open("/tmp/test.svg","w") do io
#   show(io, MIME"image/svg+xml"(), scene)
# end

#


text:
# models.jl
#
# author: Wade Shen
# swade@ll.mit.edu
# Copyright © 2014 Massachusetts Institute of Technology, Lincoln Laboratory
# version 0.1
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import Main.apply
export make_background, stats, vocab_size, apply, BKG

type BKG
  dict  :: Associative{AbstractString, Int32}
  index :: Array{AbstractString}
  stats :: Vector{Float64}
end
vocab_size(bkg::BKG) = length(bkg.index)
getindex(bkg::BKG, token :: AbstractString) = bkg.dict[token]
stats(bkg::BKG, s::AbstractString) = bkg.stats[bkg[s]]

function tfnorm(stats; cutoff = 1e10, squash :: Function = log)
  for i = 1:length(stats)
    stats[i] = min(cutoff, squash(1.0 / stats[i]))
  end

  return stats
end

function apply(bkg::BKG, counts)
  for i in indices(counts)
    counts[i] *= bkg.stats[i]
  end
  return counts
end

function make_background(features; mincount = 1, prune = 0.0, unk = true, norm = stats -> min(1.0 ./ stats, 1e10), logger = Log(STDERR))
  dict = DefaultDict(AbstractString, Int32, 0)

  @timer logger "building background dictionary" begin
  # Count
  for fv in features
    for f in fv
      dict[f] += 1
    end
  end

  # clean
  unkcount = 0
  total    = 0
  for (k, v) in dict
    if v < mincount
      unkcount += v
      delete!(dict, k)
    end
    total += v
  end
  dict[unk_token] = unkcount
  end # timer

  # index
  index          = (AbstractString)[unk_token]
  rev            = DefaultDict(AbstractString, Int32, 1)
  rev[unk_token] = 1
  i              = 2
  @timer logger "building index" begin
  for (k, v) in dict
    if k != unk_token
      push!(index, k)
      rev[k] = i
      i += 1
    end
  end
  end # timer

  # make model vector
  stats = Array(Float64, i-1)
  @timer logger "making bkg vector and pruning" begin
  for (k, v) in dict
    stats[rev[k]] = v
  end
  stats /= total

  # prune
  for i = 1:length(stats)
    if stats[i] < prune
      stats[i] = 0.0
    end
  end
  end # timer

  return BKG(rev, index, norm(stats))
end


text:
import Base.Random: rand

export bwmutate, bwmutate!

@doc """bwmutate(p::Population, ls::Landscape, mutprob::Float64)

Return a new population such that each locus of each individual in the original
population has been mutated with probability `mutprob`.
"""
function bwmutate(p::Population, mutprob::Float64)
  np = Population(p)
  bwmutate!(np, mutprob)
  return np
end

@doc """bwmutate(p::MetaPopulation, ls::Landscape, mutprob::Float64)
"""
function bwmutate(mp::MetaPopulation, mutprob::Float64)
  np = MetaPopulation(mp)
  bwmutate!(np, mutprob)
  return np
end

@doc """bwmutate(g::Genotype, mutprob::Float64)
"""
function bwmutate(g::Genotype, mutprob::Float64)
  ng = Genotype(g)
  bwmutate!(ng, mutprob)
  return ng
end

@doc """bwmutate!(g::Genotype, mutprob::Float64)

Mutate the genotype in-place by mutating each locus with probability `mutprob`.
"""
function bwmutate!(g::Genotype, mutprob::Float64)
  mask = AlleleMask(0)
  current = AlleleMask(1)
  for _ = 1:g.landscape.n
    if rand() < mutprob
      mask = mask | current
    end
    current = current << 1
  end
  g.alleles = g.alleles $ mask
end

@doc """bwmutate!(p::Population, ls::Landscape, mutprob::Float64)

Mutate the population in-place by mutating each locus of each individual in the
population with probability `mutprob`.
"""
function bwmutate!(p::Population, mutprob::Float64)
  for g = p.genotypes
    bwmutate!(g, mutprob)
  end
end

@doc """bwmutate!(mp::MetaPopulation, ls::Landscape, mutprob::Float64)
"""
function bwmutate!(mp::MetaPopulation, mutprob::Float64)
  for p = mp.populations
    bwmutate!(p, mutprob)
  end
end


text:
using ITensors
using ITensorPySCF

molecule_N2 = "
N 0 0 0;
N 0 0 1;
"

molecule_H20 = "
O 1.2091536548 1.7664118189 -0.0171613972;
H 2.1984800075 1.7977100627 0.0121161719;
H 0.9197881882 2.458018557 0.629793883;
"

molecule = molecule_H20

basis = "sto3g"

mol = gto.M(atom=molecule, basis=basis, verbose=3)

# Run HF
mf = scf.RHF(mol).run()
println("Nuclear Repulsion (Ha): ", mf.energy_nuc())
println("RHF Energy (Ha): ", mf.e_tot)

# Create shorthands for 1- and 2-body integrals in MO basis
mo = mf.mo_coeff
one_body = mo' * mf.get_hcore() * mo
two_body = mol.ao2mo(mf.mo_coeff, aosym=1)

# FCI (i.e. exact diagonalization)
cisolver = fci.FCI(mf)
cisolver.kernel()
println("FCI Energy (Ha): ", cisolver.e_tot)

t = one_body
V = two_body
n_occ = mf.mo_occ
e_nuclear = mf.energy_nuc()

@show e_nuclear

n = size(t, 1)
V = reshape(V, n, n, n, n)

# convert from chemistry convention to physics convention
V = 0.5 * permutedims(V, (3, 2, 1, 4))

os = OpSum()
os += e_nuclear, "Id", 1
for i in 1:n, j in 1:n
  os .+= t[i, j], "Cdagup", i, "Cup", j
  os .+= t[i, j], "Cdagdn", i, "Cdn", j
end
for i in 1:n, j in 1:n, k in 1:n, l in 1:n
  os .+= V[i, j, k, l], "Cdagup", i, "Cdagup", j, "Cup", k, "Cup", l
  os .+= V[i, j, k, l], "Cdagdn", i, "Cdagdn", j, "Cdn", k, "Cdn", l
  os .+= V[i, j, k, l], "Cdagup", i, "Cdagdn", j, "Cdn", k, "Cup", l
  os .+= V[i, j, k, l], "Cdagdn", i, "Cdagup", j, "Cup", k, "Cdn", l
end

s = siteinds("Electron", n; conserve_nf=true)
H = MPO(os, s)

occ_to_state = Dict([0 => 1, 2 => 4])
occupation_to_state(n) = occ_to_state[n]
ψmf = MPS(s, occupation_to_state.(n_occ))
@show inner(ψmf', H, ψmf)

ψ0 = randomMPS(s, occupation_to_state.(n_occ); linkdims=1)
@show inner(ψ0', H, ψ0)

sweeps = Sweeps(10)
setmaxdim!(sweeps, 100, 200)
setcutoff!(sweeps, 1e-6)
setnoise!(sweeps, 1e-6, 1e-7, 1e-8, 0.0)
e, ψ = dmrg(H, ψ0, sweeps)

@show e - cisolver.e_tot


text:
"Utility functions"

using Zygote: gradient

function Jacobian(F, x)
    y = F(x)
    n = length(y)
    m = length(x)
    T = eltype(y)
    j = Array{T, 2}(undef, n, m)
    for i in 1:n
        j[i, :] .= gradient(x -> F(x)[i], x)[1]
    end
    return j
end

function wMatrix(γ, order)
    mW = huge*Matrix{Float64}(I, order, order)
    mW[1, 1] = γ
    return mW
end

function transition(γ, order)
    V = zeros(order, order)
    V[1] = 1/γ
    return V
end

function shift(dim)
    S = Matrix{Float64}(I, dim, dim)
    for i in dim:-1:2
           S[i,:] = S[i-1, :]
    end
    S[1, :] = zeros(dim)
    return S
end

function uvector(dim, pos=1)
    u = zeros(dim)
    u[pos] = 1
    return dim == 1 ? u[pos] : u
end

function defineCS(order)
    uvector(order), shift(order)
end

max_stars_repo_name: wmkouw/NLARX


text:
module FermiContour

module Transforms

include("JacobiFunc.jl")

import Base.(|>)

immutable Sn{M}
    m::M
end
@inline function |>(p, s::Sn)
    w,z = p
    sn,cn,dn = JacobiFunc.jacobi(z, s.m)
    return w*cn*dn, sn
end

immutable Möbius{A,B,C,D}
    a::A; b::B; c::C; d::D
end
@inline function |>(p, m::Möbius)
    w,z = p
    return w*(m.a*(m.c*z+m.d ) - m.c*(m.a*z+m.b))/(m.c*z+m.d)^2,
          (m.a*z+m.b)/(m.c*z+m.d)
end

immutable Affine{A,B}
    a::A; b::B
end
@inline function |>(p, a::Affine)
    w,z = p
    return w*a.a, a.a*z+a.b
end

immutable Sqrt end
@inline function |>(p, s::Sqrt)
    w,z = p
    return w/(2*sqrt(z)), sqrt(z)
end

end # module Transforms

include("JacobiFunc.jl")

export fermicontour
"""
    fermicontour(E1,E2,β,n) -> w,z

Compute a quadrature rule for evaluating the Fermi-Dirac function through contour integration.

For any `x ∈ [-E2,-E1] ∪ [E1,E2]` and any `β ∈ [0,∞]` it holds
```julia
fermidirac(x,β) ≈ sum(
    real(w*fermidirac(z,β)/(z-x))
    for (w,z) in zip(fermicontour(E1,E2,β,n)...)
)
```
and the error converges exponentially for `n → ∞`. Note that both `E1 = 0` (single interval)
and `β = Inf` (zero temperature) are valid inputs.

Reference: Lin Lin et al., Pole-Based Approximation of the Fermi-Dirac Function, Chinese Annals of Mathematics, Series B
"""
function fermicontour(E1,E2,β,n)
    m = E1^2 + π^2/β^2
    M = E2^2 + π^2/β^2

    k = (sqrt(M/m)-1)/(sqrt(M/m)+1)
    K = JacobiFunc.K(k^2)
    iK = im*JacobiFunc.iK(k^2)

    t = -K + iK/2 + 4*K*(0.5:1:2*n-0.5)/(2*n)
    w = Vector{Complex{typeof(k)}}(length(t))
    z = Vector{Complex{typeof(k)}}(length(t))
    for i = 1:length(t)
        # Quadrature points and weights
        w[i],z[i] = (-1/(2π*im) * 4*K/n, t[i]) |>
            Transforms.Sn(k^2) |>
            Transforms.Möbius(1,1/k, -1,1/k) |>
            Transforms.Affine(-sqrt(m*M), π^2/β^2) |>
            Transforms.Sqrt() |>
            Transforms.Affine(im,0)
    end
    return w,z
end

export fermidirac
"""
    fermidirac(z,β)

Compute `1/(1+exp(β*z))` in a numerically stable way.
"""
function fermidirac(z,β)
    return abs(β*z) < log(realmax(real(typeof(z)))) ?
            1/(1 + exp(β*z)) :
            0.5*(1-sign(real(z)))
end

end # module


text:
using Clang

include("cbasic.jl")
include("functions.jl")
include("wrap_c.jl")


text:
data = readlines()

function inner_string_size(s)
    n = 0
    pass_next = 0
    for (c1, c2) in zip(s[2:end-1], s[3:end])
        if pass_next > 0
            pass_next -= 1
            continue
        end
        if c1 == '\\'
            if c2 == 'x'
                pass_next = 3
            else
                pass_next = 1
            end
        end
        n += 1
    end
    n
end

function solve_1()
    tot_len = data .|> length |> sum
    tot_inner = data .|> inner_string_size |> sum
    tot_len - tot_inner
end

function outer_string_size(s)
    count(occursin("\\\""), s) + 2
end

function solve_2()
    tot_out = data .|> outer_string_size |> sum
    tot_out
end

solve_1() |> println
solve_2() |> println

max_stars_repo_name: pLOPeGG/JuliaOfCode

max_issues_repo_name: pLOPeGG/JuliaOfCode

max_forks_repo_name: pLOPeGG/JuliaOfCode


text:
using Base.Threads, LinearAlgebra

export npairs, nsites, maxneigs, max_neighbours, neigs, neighbours, neigs!

PairList(X::Vector{SVec{T}}, cutoff::AbstractFloat, cell::AbstractMatrix, pbc;
            int_type::Type = Int32, fixcell = true) where {T} =
   _pairlist_(X, SMat{T}(cell), SVec{Bool}(pbc), T(cutoff), int_type,
              fixcell)

PairList(X::Matrix{T}, args...; kwargs...) where {T} =
   PairList(reinterpret(SVec{T}, X, (size(X,2),)), args...; varargs...)

npairs(nlist::PairList) = length(nlist.i)
nsites(nlist::PairList) = length(nlist.first) - 1
cutoff(nlist::PairList) = nlist.cutoff

# ====================== Cell Index Algebra =====================

"Map i back to the interval [0,n) by shifting by integer multiples of n"
@inline function bin_wrap(i::Integer, n::Integer)
   while i <= 0;  i += n; end
   while i > n;  i -= n; end
   return i;
end

"apply bin_wrap only if periodic bdry"
@inline bin_wrap(i::Integer, pbc::Bool, n::Integer) = pbc ? bin_wrap(i, n) : i

"Map i back to the interval [0,n) by assigning edge value if outside interval"
@inline function bin_trunc(i::TI, n::TI) where {TI <: Integer}
   if i <= 0;     i = TI(1)
   elseif i > n;  i = n
   end
   return i
end

"apply bin_trunc only if open bdry"
@inline bin_trunc(i::TI, pbc::Bool, n::TI) where {TI <: Integer} =
      pbc ? TI(i) : bin_trunc(i, n)

"applies bin_trunc to open bdry and bin_wrap to periodic bdry"
@inline bin_wrap_or_trunc(i::Integer, pbc::Integer, n::Integer) =
      pbc ? bin_wrap(i, n) : bin_trunc(i, n)

"Map particle position to a (cartesian) cell index"
@inline position_to_cell_index(inv_cell::SMat{T}, x::SVec{T}, ns::SVec{TI}
         ) where {T, TI <: Integer} =
   floor.(TI, ((inv_cell' * x) .* ns .+ 1))

# ------------ The next two functions are the only dimension-dependent
#              parts of the code!

# an extension of sub2ind for the case when i is a vector (cartesian index)
# @inline Base.sub2ind(dims::NTuple{3,TI}, i::SVec{TI}) where {TI <: Integer} =
#    sub2ind(dims, i[1], i[2], i[3])
# WARNING: this smells like a performance regression!
# WARNING: `LinearIndices` always returnsstores Int, not TI!!!
@inline _sub2ind(dims::NTuple{3,TI}, i::SVec{TI}) where {TI <: Integer} =
   TI( (LinearIndices(dims))[i[1], i[2], i[3]] )

lengths(C::SMat{T}) where {T} =
   det(C) ./ SVec{T}(norm(C[2,:]×C[3,:]), norm(C[3,:]×C[1,:]), norm(C[1,:]×C[2,:]))

# --------------------------------------------------------------------------

function analyze_cell(cell, cutoff, TI)
   @assert TI <: Integer
   # check the cell volume (allow only 3D volumes!)
   volume = abs(det(cell))
   if volume < 1e-12
      @warn("zero volume detected - proceed at your own risk")
      @show cell
      @show volume
      @show cutoff
   end
   # precompute inverse of cell matrix for coordiate transformation
   inv_cell = inv(cell)
   # Compute distance of cell faces
   lens = abs.(lengths(cell))
   # Number of cells for cell subdivision
   _t = floor.(TI, lens / cutoff)
   ns_vec = max.(_t, one(TI))
   return inv_cell, ns_vec, lens
end

# multi-threading setup

function setup_mt(niter::TI, maxnt = MAX_THREADS[1]) where TI <: Integer
   nt = minimum([6, nthreads(), ceil(TI, niter / 20), maxnt])
   # nn = ceil.(TI, linspace(1, niter+1, nt+1))
   # range(start, stop=stop, length=length)
   nn = ceil.(TI, range(1, stop=niter+1, length=nt+1))
   return nt, nn
end

# ==================== CellList Core  ================

function _celllist_(X::Vector{<:SVec}, cell::SMat, pbc::SVec{Bool},
            cutoff::T, TI) where {T <: Real}
   @assert TI <: Integer
   # ----- analyze cell -----
   nat = length(X)
   inv_cell, ns_vec, lens = analyze_cell(cell, cutoff, TI)
   ns = ns_vec.data

   if prod(BigInt.(ns_vec)) > typemax(TI)
      error("""Ratio of simulation cell size to cutoff is very large.
               Are you using a cell with lots of vacuum? To fix this
               use a larger integer type (e.g. Int128), a
               larger cut-off, or a smaller simulation cell.""")
   end

   # data structure to store a linked list for each bin
   ncells = prod(ns_vec)
   seed = fill(TI(-1), ncells)
   last = Vector{TI}(undef, ncells)
   next = Vector{TI}(undef, nat)
   nats = zeros(TI, ncells)

   for i = 1:nat
      # Get cell index
      c = position_to_cell_index(inv_cell, X[i], ns_vec)
      # Periodic/non-periodic boundary conditions
      c = bin_wrap_or_trunc.(c, pbc, ns_vec)
      # linear cell index  # (+1 due to 1-based indexing)
      ci = _sub2ind(ns, c)

      # Put atom into appropriate bin (list of linked lists)
      if seed[ci] < 0   #  ci contains no atom yet
         next[i] = -1
         seed[ci] = i
         last[ci] = i
         nats[ci] += 1
      else
         next[i] = -1
         next[last[ci]] = i
         last[ci] = i
         nats[ci] += 1
      end
   end

   return CellList(X, cell, inv_cell, pbc, cutoff, seed, last, next, nats)
end

function _pairlist_(clist::CellList{T, TI}) where {T, TI}

   X, cell, pbc, cutoff, seed, last, next =
         clist.X, clist.cell, clist.pbc, clist.cutoff,
         clist.seed, clist.last, clist.next
   nat = length(X)
   inv_cell, ns_vec, lens = analyze_cell(cell, cutoff, TI)

   # guess how many neighbours per atom
   #    atoms in cell x (8 cells) * (ball / cube)
   max_nat_cell = maximum(clist.nats)
   nneigs_guess = ceil(TI, 1.5 * max_nat_cell * (π^2/4))

   # allocate arrays for the many threads
   # set number of threads
   nt, nn = setup_mt(nat)
   # allocate arrays
   first_t = Vector{TI}[ Vector{TI}()  for n = 1:nt ]    # i
   secnd_t = Vector{TI}[ Vector{TI}()  for n = 1:nt ]    # j
   shift_t = Vector{SVec{TI}}[ Vector{SVec{TI}}()  for n = 1:nt ]  # ~ X[i] - X[j]
   # give size hints
   sz = (nat ÷ nt + nt) * nneigs_guess
   for n = 1:nt
      sizehint!(first_t[n], sz)
      sizehint!(secnd_t[n], sz)
      sizehint!(shift_t[n], sz)
   end

   # We need the shape of the bin ( bins[:, i] = cell[i,:] / ns[i] )
   # bins = cell' ./ ns_vec
   bins = hcat( cell[1,:]/ns_vec[1], cell[2,:]/ns_vec[2], cell[3,:] / ns_vec[3] )

   # Find out over how many neighbor cells we need to loop (if the box is small)
   nxyz = ceil.(TI, cutoff * (ns_vec ./ lens))
   # cxyz = CartesianIndex(nxyz.data)
   # WARNING : 3D-specific hack; also potential performance regression
   # WARNING: `CartesianIndices` always stores Int, not TI!!!
   xyz_range = CartesianIndices((-nxyz[1]:nxyz[1],
                                 -nxyz[2]:nxyz[2],
                                 -nxyz[3]:nxyz[3]))

   # Loop over threads
   # @threads for it = 1:nt
   for it = 1:nt
      for i = nn[it]:(nn[it+1]-1)
         # current atom position
         _find_neighbours_!(i, clist, ns_vec, bins, xyz_range,
                      first_t[it], secnd_t[it], shift_t[it])
      end # for i = 1:nat
   end # @threads

   # piece them back together
   sz = sum( length(first_t[i]) for i = 1:nt )
   first = first_t[1];     sizehint!(first, sz)
   secnd = secnd_t[1];     sizehint!(secnd, sz)
   shift = shift_t[1]; sizehint!(shift, sz)
   for it = 2:nt
      append!(first, first_t[it])
      append!(secnd, secnd_t[it])
      append!(shift, shift_t[it])
   end

   # Build return tuple
   return first, secnd, shift
end

function _find_neighbours_!(i, clist, ns_vec::SVec{TI}, bins, xyz_range,
                            first, secnd, shift) where TI
   inv_cell, X, pbc = clist.inv_cell, clist.X, clist.pbc
   seed, last, next = clist.seed, clist.last, clist.next
   xi = X[i]

   # funnily testing with cutoff^2 actually makes a measurable difference
   # which suggests we are pretty close to the performance limit
   cutoff_sq = clist.cutoff^2

   # cell index (cartesian) of xi
   ci0 = position_to_cell_index(inv_cell, xi, ns_vec)

   # Truncate if non-periodic and outside of simulation domain
   # (here, we don't yet want to wrap the pbc as well)
   ci = bin_trunc.(ci0, pbc, ns_vec)
   # dxi is the position relative to the lower left corner of the bin
   dxi = xi - bins * (ci .- 1)

   # Apply periodic boundary conditions as well now
   ci = bin_wrap_or_trunc.(ci0, pbc, ns_vec)

   for ixyz in xyz_range  # Integer
      # convert cartesian index to SVector
      xyz = SVec{TI}(ixyz.I)
      # get the bin index
      cj = bin_wrap.(ci + xyz, pbc, ns_vec)
      # skip this bin if not inside the domain
      all(TI(1) .<= cj .<= ns_vec) || continue
      # linear cell index
      ncj = _sub2ind(ns_vec.data, cj)
      # Offset of the neighboring bins
      off = bins * xyz

      # Loop over all atoms in neighbouring bin (all potential
      # neighbours in the bin with linear index cj1)
      j = seed[ncj] # the first atom in the ncj cell
      while j > 0
         if i != j || any(xyz .!= 0)
            xj = X[j] # position of current neighbour

            # we need to find the cell index again, because this is
            # not really the cell index, but it could be outside
            # the domain -> i.e. this only makes a difference for pbc
            cj = position_to_cell_index(inv_cell, xj, ns_vec)
            cj = bin_trunc.(cj, pbc, ns_vec)

            # drj is position relative to lower left corner of the bin
            dxj = xj - bins * (cj .- 1)
            # Compute distance between atoms
            dx = dxj - dxi + off
            norm_dx_sq = dot(dx, dx)

            # append to the list
            if norm_dx_sq < cutoff_sq
               push!(first, i)
               push!(secnd, j)
               push!(shift, (ci0 - cj + xyz) .÷ ns_vec)
            end
         end  # if i != j || any(xyz .!= 0)

         # go to the next atom in the current cell
         j = next[j];
      end # while j > 0 (loop over atoms in current cell)
   end # loop over neighbouring bins
end

function _pairlist_(X::Vector{SVec{T}}, cell::SMat{T}, pbc::SVec{Bool},
            cutoff::T, TI, fixcell::Bool) where {T}
   @assert TI <: Integer
   # temporary (?) fix to make sure all atoms are within the cell
   if fixcell
      X, cell = _fix_cell_(X, cell, pbc)
   end

   clist = _celllist_(X, cell, pbc, cutoff, TI)
   i, j, S = _pairlist_(clist)

   first = get_first(i, length(X))
   sort_neigs!(j, (S,), first)

   return PairList(X, cell, cutoff, i, j, S, first)
end

# ==================== Post-Processing =========================

"""
`get_first(i::Vector{TI}, nat::Integer) -> first::Vector{TI}
where TI <: Integer`

Assumes that `i` is sorted in ascending order.
For `n = 1, . . ., nat`, `first[n]` will be the index to the first
element of `i` with value `n`. Further, `first[nat+1]` will be
`length(i) + 1`.

If `first[n] == first[n+1]` then this means that `i` contains no element `n`.
"""
function get_first(i::Vector{TI}, nat::Integer = i[end]) where {TI}
   # compute the first index for each site
   first = Vector{TI}(undef, nat + 1)
   idx = 1
   n = 1
   while n <= nat && idx <= length(i)
      first[n] = idx
      while i[idx] == n
         idx += 1
         if idx > length(i)
            break
         end
      end
      n += 1
   end
   first[n:end] .= length(i)+1
   return first
end

"""
`sort_neigs!(j, arrays, first)`

sorts each sub-range of `j` corresponding to one site  in ascending order
and applies the same permutation to `r, R, S`.
"""
function sort_neigs!(j, arrays::Tuple, first)
   nat = length(first) - 1
   nt, nn = setup_mt(nat)
   @threads for it = 1:nt
      for n = nn[it]:(nn[it+1]-1)
         if first[n+1] > first[n] + 1
            rg = first[n]:first[n+1]-1
            jrg = j[rg]
            if !issorted(jrg)
               I = sortperm(j[rg])
               rg_perm = rg[I]
               j[rg] = j[rg_perm]
               for a in arrays
                  a[rg] .= a[rg_perm]
               end
            end
         end
      end
   end
end

"""
`_fix_cell_(X::Vector{SVec{T}}, C::SMat{T}, pbc)`

produces new `X`, `C` such that PBC are respected, but all positions are
inside the cell. This Potentially involves
 * wrapping atom positions in the pbc directions
 * shifting atom positions in the non-pbc directions
 * enlarging the cell
If either `X` or `C` are modified, then they will be copied.
"""
function _fix_cell_(X::Vector{SVec{T}}, C::SMat{T}, pbc) where {T}
   invCt = inv(C)'
   copy_X = false
   min_lam = @MVector zeros(3)
   max_lam = @MVector ones(3)
   for (ix, x) in enumerate(X)
      λ = Vector(invCt * x)
      for i = 1:3
         update_x = false
         if !(0.0 <= λ[i] < 1.0)
            if pbc[i]
               λ[i] = mod(λ[i], 1.0)
               update_x = true
            else
               min_lam[i] = min(min_lam[i], λ[i])
               max_lam[i] = max(max_lam[i], λ[i])
            end
         end
         if update_x
            if !copy_X; X = copy(X); copy_X = true end
            X[ix] = C' * SVec{T}(λ)
         end
      end
   end
   # check whether we need to adjust the non-PBC directions
   if (minimum(min_lam) < 0.0) || (maximum(max_lam) > 1.0)
      # shift vector:
      if !copy_X; X = copy(X); copy_X = true end
      t = - C' * min_lam
      for n = 1:length(X)
         X[n] += t
      end
      # the new min_lam is now zero and the new max_lam is max_lam - min_lam
      max_lam -= min_lam
      min_lam .= 0.0
      # we need to multiply the cell by the correct lambda
      C = hcat( max_lam[1] * C[1,:], max_lam[2] * C[2,:], max_lam[3] * C[3,:])'
   end
   return X, C
end

"""
`maxneigs(nlist::PairList) -> Integer`

returns the maximum number of neighbours that any atom in the
neighbourlist has.
"""
maxneigs(nlist::PairList) = maximum( nneigs(nlist, n) for n = 1:nsites(nlist) )

# retire max_neigs
const max_neigs = maxneigs

"""
`nneigs(nlist::PairList, i0::Integer) -> Integer` :
return number of neighbours of particle with index `i0`.
(within cutoff radius + buffer)
"""
nneigs(nlist::PairList, i0::Integer) = nlist.first[i0+1]-nlist.first[i0]

function _getR(nlist, idx)
   i = nlist.i[idx]
   j = nlist.j[idx]
   return _getR(nlist.X[j] - nlist.X[i], nlist.S[idx], nlist.C)
end

_getR(dX::SVec, S::SVec, C::SMat) = dX + C' * S

"""
`neigs!(Rtemp, nlist, i) -> j, R`

For `nlist::PairList` this returns the interaction neighbourhood of
the atom indexed by `i`. E.g., in the standard loop approach one
would have
```
Rtemp = zeros(JVecF, maxneigs(nlist))
for (i, j, R) in sites(nlist)
   (j, R) == neigs(nlist, i) == neigs!(Rtemp, nlist, i)
end
```

`R` is a view into `Rtemp` with the correct length, while `j` is a view
into `nlist.j`.
"""
function neigs!(Rs::AbstractVector{<: SVec}, nlist::PairList, i0::Integer)
   n1, n2 = nlist.first[i0], nlist.first[i0+1]-1
   _grow_array!(Rs, n2-n1+1)
   J = (@view nlist.j[n1:n2])
   for n = 1:length(J)
      Rs[n] = _getR(nlist, n1+n-1)
   end
   return J, (@view Rs[1:length(J)])
end

function neigs!(Js::AbstractVector{<: SVec},
                Rs::AbstractVector{<: SVec},
                nlist::PairList, i0::Integer)
   _grow_array!(Rs, n2-n1+1)
   _grow_array!(Js, n2-n1+1)
   j, Rs = neigs!(Rs, nlist, i0)
   N = length(j)
   copyto!(Js, j)
   return (@view Js[1:length(j)]), Rs
end

function _grow_array!(A::Vector{T}, N) where {T}
   if length(A) < N
      append!(A, zeros(T, N-length(A)))
   end
   return A
end

"""
`neigss!(Rs, nlist, i0) -> j, R, S` : return neighbourhood as in
`neigs!` as well as the corresponding cell shifts.

(`R` is a view into `Rs` with the correct length)
"""
function neigss!(Rs::AbstractVector{<: SVec}, nlist::PairList, i0::Integer)
   n1, n2 = nlist.first[i0], nlist.first[i0+1]-1
   _grow_array!(Rs, n2-n1+1)
   J = (@view nlist.j[n1:n2])
   for n = 1:length(J)
      Rs[n] = _getR(nlist, n1+n-1)
   end
   return J, (@view Rs[1:length(J)]), (@view nlist.S[n1:n2])
end

neigs(nlist::PairList{T}, i0::Integer) where {T} =
      neigs!( zeros(SVec{T}, nneigs(nlist, i0)), nlist, i0 )

neigss(nlist::PairList{T}, i0::Integer) where {T} =
      neigss!( zeros(SVec{T}, nneigs(nlist, i0)), nlist, i0 )

"""
alias for `neigs`
"""
neighbours = neigs

"""
alias for `max_neigs`
"""
max_neighbours = maxneigs


text:
include("../src/Dataclass.jl")
using .DataClass

@dataclass mutable struct Register
    value::Int64
    valuestr::String
    _initvars = [_init=true, _repr=true, _eq=true, _order=true, _unsafe_hash=false, _frozen=true]
end


text:
@testset "672.bulb-switcher-ii.jl" begin
    @test flip_lights(0, 0) == 1
    @test flip_lights(2, 1) == 3
    @test flip_lights(2, 2) == 4
    @test flip_lights(2, 3) == 4
end


text:
# misc constraint tests
@testset "node bc" begin
    grid = generate_grid(Triangle, (1, 1))
    addnodeset!(grid, "nodeset", x-> x[2] == -1 || x[1] == -1)
    dh = DofHandler(grid)
    push!(dh, :u, 2)
    push!(dh, :p, 1)
    close!(dh)
    ch = ConstraintHandler(dh)
    dbc1 = Dirichlet(:u, getnodeset(grid, "nodeset"), (x,t) -> x, [1, 2])
    dbc2 = Dirichlet(:p, getnodeset(grid, "nodeset"), (x,t) -> 0, 1)
    add!(ch, dbc1)
    add!(ch, dbc2)
    close!(ch)
    update!(ch)

    @test ch.prescribed_dofs == collect(1:9)
    @test ch.values == [-1, -1, 1, -1, -1, 1, 0, 0, 0]
end


text:
using TBLIS
using Test
#TBLIS.init()

O = 2
V = 2
TT = Float64

A = rand(TT,O,O,O,O)
B = rand(TT,O,O,V,V)
C = zeros(TT,O,O,V,V)
tijab = TBLIS.TTensor{TT}(A,1.0)
WmBeJ = TBLIS.TTensor{TT}(B,1.0)
T2 = TBLIS.TTensor{TT}(C,1.0)

#@tensor C[i,j,a,b] = A[i,j,k,l]*B[k,l,a,b]
for i=1:O,j=1:O,a=1:V,b=1:V
    for k=1:O,l=1:O
        C[i,j,a,b] += A[i,j,k,l]*B[k,l,a,b]
    end
end
_C = deepcopy(C)
C .= 0
TBLIS.mul!(T2,tijab,WmBeJ,"ijkl","klab","ijab")
@test isapprox(_C,C,rtol=1E-10)

_B = deepcopy(B) + deepcopy(C)
TBLIS.add!(WmBeJ,T2,"ijab","ijab")

@test T2.data ≈ C
@test isapprox(_B,C,rtol=1E-10)

@testset "Test changing the number of TBLIS threads" begin
  tblis_num_threads = TBLIS.get_num_threads()

  TBLIS.set_num_threads(4)
  @test TBLIS.get_num_threads() == 4

  TBLIS.set_num_threads(1)
  @test TBLIS.get_num_threads() == 1

  # Set the number of TBLIS threads back to the original value
  TBLIS.set_num_threads(tblis_num_threads)
  @test TBLIS.get_num_threads() == tblis_num_threads
end


text:
a1=solveQ([:x,:y,:z],[[:x は :y を 作る],[:z は :y です]],魔法会議)
f1=statements([:x,:y,:z],[[:x は :y を 作る],[:z は :y です]],a1)
u1=unique(f1)

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
"""
    Polylogarithms

Module containing functions to calculate the polylogarithm and associated functions

"""
module Polylogarithms
import SpecialFunctions

# not using MPFR for the moment
# using Base.MPFR: ROUNDING_MODE, big_ln2

export polylog, bernoulli, harmonic, stieltjes, dirichlet_beta
export Diagnostics
export parse

# @compat ComplexOrReal{T} = Union{T,Complex{T}}
# s::ComplexOrReal{Float64}
ComplexOrReal{T} = Union{T,Complex{T}}

# Some extra parsing routines for reading Mathematics output, but also, complex numbers
include("utilities.jl")

# Constants
include("constants.jl")

# Series
include("stieltjes.jl")
include("bernoulli_n.jl")
include("harmonic.jl")
# include("gamma_derivatives.jl") # this just generates a table that we include into the code

# Functions
include("beta.jl")
include("bernoulli_poly.jl")
include("polylog.jl")

end


text:
using Modia
using Tables
using DataFrames
using StatsPlots
using BenchmarkTools

SIR = Model(
    equations = :[
        N = S + I + R
        der(S) = -β*c*I/N*S
        der(I) = β*c*I/N*S - γ*I
        der(R) = γ*I
    ]
);

δt = 0.1
tmax = 40.0;

u0 = Map(S = Var(init=990.0),
         I = Var(init=10.0),
         R = Var(init=0.0));

p = Map(β = 0.05,
        c = 10.0,
        γ = 0.25);

model = SIR | u0 | p;

sir = @instantiateModel(model);

simulate!(sir,Tsit5(),stopTime=tmax,interval=δt);

result = sir.result_x.u # extract result
result = hcat(result...) # concatenate vectors
result = result' # transpose
result = Tables.table(result) # convert to table
df_modia = DataFrame(result) # convert to DataFrame
rename!(df_modia,["S","I","R"]) # rename
df_modia[!,:t] = sir.result_x.t; # add in time

@df df_modia plot(:t,
    [:S :I :R],
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number")

@benchmark simulate!(sir,Tsit5(),stopTime=tmax,interval=δt)


text:
mutable struct ObjectiveSurface
    param1_name::String
    param1_values::Array{Float64,1}
    param2_name::String
    param2_values::Array{Float64,1}
    surface_points::Int64
    replications::Int64
    params::Dict{String, Any}
    weight_matrix::Any
end

function objective_function(parameter_dict::Dict, os::ObjectiveSurface, seed, ind)
    try
        parameter_dict[os.param1_name] = os.param1_values[ind]
        parameter_dict[os.param2_name] = os.param2_values[ind]
        log_price_paths = log.(SLOB(parameter_dict)(seed))
        log_price_paths = vcat(log_price_paths[1:1,:], log_price_paths)
        log_returns = diff(log_price_paths, dims=1)
        objective_value = weighted_moment_distance(os.weight_matrix,
            log_returns)
        return objective_value
    catch e
        return NaN
    end
end

function (os::ObjectiveSurface)(seed, parallel=false)
    iterations = os.surface_points^2

    seeds = Int.(rand(MersenneTwister(seed), UInt32, iterations))
    parameter_vector = fill(os.params, iterations)
    if parallel==true
        os_values = SharedArray{Float64,1}(iterations)
        @sync @distributed for i in 1:iterations
            os_values[i] = objective_function(parameter_vector[i], os, seeds[i], i)
        end
        return Array(os_values)
    else
        os_values = Array{Float64,1}(undef, iterations)
        @showprogress 1 "Iteration: " for i in 1:iterations
            os_values[i] = objective_function(parameter_vector[i], os, seeds[i], i)
        end
        return os_values
    end
end


text:
struct LinearInterpolationData{uType,tType} <: AbstractDiffEqInterpolation
  timeseries::uType
  ts::tType
end

DiffEqBase.interp_summary(::LinearInterpolationData) = "1st order linear"
(interp::LinearInterpolationData)(tvals,idxs,deriv,p,continuity::Symbol=:left) = sde_interpolation(tvals,interp,idxs,deriv,p,continuity)
(interp::LinearInterpolationData)(val,tvals,idxs,deriv,p,continuity::Symbol=:left) = sde_interpolation!(val,tvals,interp,idxs,deriv,p,continuity)


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
mutable struct Mean_covar
    mu::Array{Float64,2}
    sigma::Array{Float64,2}
    w::Float64
end

const B = Array{Mean_covar,1}(undef, 5)
display(B)

# %%
B[1] = Mean_covar(zeros(2,2), zeros(2,2), 0.0)
display(B)

# %%

max_stars_repo_path: 0014/const array.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0014/const array.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0014/const array.jl

max_forks_repo_name: genkuroki/public


text:
const JLL_ENV_CLANG_TARGETS_MAPPING = Dict(
    "aarch64-apple-darwin20"=>"aarch64-apple-darwin20",
    "aarch64-linux-gnu"=>"aarch64-unknown-linux-gnu",
    "aarch64-linux-musl"=>"aarch64-unknown-linux-musl",
    "armv7l-linux-gnueabihf"=>"armv7l-unknown-linux-gnueabihf",
    "armv7l-linux-musleabihf"=>"armv7l-unknown-linux-musleabihf",
    "i686-linux-gnu"=>"i686-unknown-linux-gnu",
    "i686-linux-musl"=>"i686-unknown-linux-musl",
    "i686-w64-mingw32"=>"i686-w64-windows-gnu",
    "powerpc64le-linux-gnu"=>"powerpc64le-unknown-linux-gnu",
    "x86_64-apple-darwin14"=>"x86_64-apple-darwin14",
    "x86_64-linux-gnu"=>"x86_64-unknown-linux-gnu",
    "x86_64-linux-musl"=>"x86_64-unknown-linux-musl",
    "x86_64-unknown-freebsd11.1"=>"x86_64-unknown-freebsd11.1",
    "x86_64-w64-mingw32"=>"x86_64-w64-windows-gnu",
)

target(triple::String) = get(JLL_ENV_CLANG_TARGETS_MAPPING, triple, "unknown")
target(p::Platform) = target(__triplet(p))