text:
"""
    equal_weights(returns::AbstractMatrix)::AbstractVector
```math
    w_i = 1/N \\quad  i \\in [1, \\hdots, N]
```
External links
* DeMiguel, V., Garlappi, L., & Uppal, R. (2009).
  Optimal versus naive diversification: How inefficient is the 1/N portfolio strategy?.
  The review of Financial studies, 22(5), 1915-1953.
  doi: [10.1093/hhm075](https://doi.org/10.1093/rfs/hhm075)
"""
function equal_weights(returns::AbstractMatrix)::AbstractVector
    """
    Input: Returns over a specific period.
    Output: The weights of equal weights portfolio allocation.
    """
    num_assets = size(returns)[2]
    weights = fill(1/num_assets, (num_assets))
    return weights
end

"""
long_only_minimum_variance
    ```math
    \\textrm{minimize}_{\\frac{\\vec{w}^T \\vec{\\sigma}}{\\sqrt{\\vec{w}^T \\Sigma \\vec{w}}}, \\\\
    \\textrm{subject to} \\quad \\vec{1}^T \\vec{w} = 1 \\quad \\vec{w} \\geq 0.
```
External links
* DeMiguel, V., Garlappi, L., & Uppal, R. (2009).
  Optimal versus naive diversification: How inefficient is the 1/N portfolio strategy?.
  The review of Financial studies, 22(5), 1915-1953.
  doi: [10.1093/hhm075](https://doi.org/10.1093/rfs/hhm075)
"""
function long_only_minimum_variance(returns::AbstractMatrix)::AbstractVector
    cov_matrix = cov(returns)
    portfolio = Model(Ipopt.Optimizer)
    set_silent(portfolio) # Suppress the output
    num_assets = size(cov_matrix)[1]
    @variable(portfolio, x[1:num_assets] >= 0)
    @objective(portfolio, Min, x' * cov_matrix * x)
    @constraint(portfolio, sum(x[i] for i = 1:num_assets) == 1)
    optimize!(portfolio)
    weights = value.(x)
    return weights
end
"""
    random_weights
```math
f(x) = \\frac{1}{\\pi \\sqrt{(x - a) (b - x)}}, \\quad x \\in [a, b]
```
"""
function random_weights(returns::AbstractMatrix)::AbstractVector
    """
    Input: Returns over a specific period.
    Output: Randomized weights of a portfolio.
    """
    num_assets = size(returns)[2]
    weights = rand(num_assets)
    norm = sum(weights)
    weights = weights ./ norm
    return weights
end
"""
    mdp
```math
    \\textrm{maximize}_{\\vec{w}} \\quad  \\frac{\\vec{w}^T \\vec{\\sigma}}{\\sqrt{\\vec{w}^T \\Sigma \\vec{w}}}, \\\\
    \\textrm{subject to} \\quad \\vec{1}^T \\vec{w} = 1.
```
External links
* Choueifaty, Y., Froidure, T., & Reynier, J. (2013).
  Properties of the most diversified portfolio.
  Journal of investment strategies, 2(2), 49-70.
  doi: [10.2139/ssrn.1895459](https://dx.doi.org/10.2139/ssrn.1895459)
"""
function mdp(returns::AbstractMatrix)::AbstractVector
    cov_matrix = cov(returns)
    var_vector = var(returns)
    portfolio = Model(Ipopt.Optimizer)
    set_silent(portfolio) # Suppress the output
    num_assets = size(cov_matrix)[1]
    @variable(portfolio, x[1:num_assets] >= 0)
    @objective(portfolio, Min, 0.5 * x' * cov_matrix * x)
    @constraint(portfolio, sum(x[i] * var_vector[i] for i = 1:num_assets) == 1)
    optimize!(portfolio)
    weights = value.(x)
    weights = weights ./ sum(weights)
    return weights
end

function inverse_variance(returns::AbstractMatrix)::AbstractVector
    return 0
end
function inverse_volatility(returns::AbstractMatrix)::AbstractVector
    return 0
end


text:
# Superselection sectors (quantum numbers):
# for defining graded vector spaces and invariant subspaces of tensor products
#==============================================================================#
"""
    abstract type Sector end

Abstract type for representing the (isomorphism classes of) simple objects in (unitary
and pivotal) (pre-)fusion categories, e.g. the irreducible representations of a finite or
compact group.

Every new `G<:Sector` should implement the following methods:
*   `one(::Type{G})`: unit element of `G`
*   `conj(a::G)`: ``a̅``, conjugate or dual label of ``a``
*   `⊗(a::G, b::G)`: iterable with unique fusion outputs of ``a ⊗ b``
    (i.e. don't repeat in case of multiplicities)
*   `Nsymbol(a::G, b::G, c::G)`: number of times `c` appears in `a ⊗ b`, i.e. the
    multiplicity
*   `FusionStyle(::Type{G})`: `Abelian()`, `SimpleNonAbelian()` or
    `DegenerateNonAbelian()`
*   `BraidingStyle(::Type{G})`: `Bosonic()`, `Fermionic()`, `Anyonic()`, ...
*   `Fsymbol(a::G, b::G, c::G, d::G, e::G, f::G)`: F-symbol: scalar (in case of
    `Abelian`/`SimpleNonAbelian`) or matrix (in case of `DegenerateNonAbelian`)
*   `Rsymbol(a::G, b::G, c::G)`: R-symbol: scalar (in case of
    `Abelian`/`SimpleNonAbelian`) or matrix (in case of `DegenerateNonAbelian`)
and optionally
*   `dim(a::G)`: quantum dimension of sector `a`
*   `frobeniusschur(a::G)`: Frobenius-Schur indicator of `a`
*   `Bsymbol(a::G, b::G, c::G)`: B-symbol: scalar (in case of
    `Abelian`/`SimpleNonAbelian`) or matrix (in case of `DegenerateNonAbelian`)
*   `twist(a::G)` -> twist of sector `a`
and optionally, if `FusionStyle(G) isa DegenerateNonAbelian`
*   `vertex_ind2label(i::Int, a::G, b::G, c::G)` -> a custom label for the `i`th copy of
    `c` appearing in `a ⊗ b`

Furthermore, `iterate` and `Base.IteratorSize` should be made to work for the singleton type
[`SectorValues{G}`](@ref).
"""
abstract type Sector end

# iterator over the values in the sector
"""
    struct SectorValues{G<:Sector}

Singleton type to represent an iterator over the possible values of type `G`, whose instance is obtained as `values(G)`. For a new `G::Sector`, the following should be defined
*   `Base.iterate(::SectorValues{G}[, state])`: iterate over the values
*   `Base.IteratorSize(::Type{SectorValues{G}})`: `HasLenght()`, `SizeUnkown()`
    or `IsInfinite()` depending on whether the number of values of type `G` is finite (and sufficiently small) or infinite; for a large number of values, `SizeUnknown()` is recommend because this will trigger the use of `GenericRepresentationSpace`.
If `IteratorSize(G) == HasLength()`, also the following must be implemented:
*   `Base.length(::SectorValues{G})`: the number of different values
*   `Base.getindex(::SectorValues{G}, i::Int)`: a mapping between an index `i` and an
    instance of `G`
*   `findindex(::SectorValues{G}, c::G)`: reverse mapping between a value `c::G` and an
    index `i::Integer ∈ 1:length(values(G))`
"""
struct SectorValues{G<:Sector} end
Base.IteratorEltype(::Type{<:SectorValues}) = HasEltype()
Base.eltype(::Type{SectorValues{G}}) where {G<:Sector} = G
Base.values(::Type{G}) where {G<:Sector} = SectorValues{G}()

# Define a sector for ungraded vector spaces
struct Trivial <: Sector
end
Base.show(io::IO, ::Trivial) = print(io, "Trivial()")
Base.show(io::IO, ::Type{Trivial}) = print(io, "Trivial")

"""
    one(::Sector) -> Sector
    one(::Type{<:Sector}) -> Sector

Return the unit element within this type of sector.
"""
Base.one(a::Sector) = one(typeof(a))
Base.one(::Type{Trivial}) = Trivial()

"""
    dual(a::Sector) -> Sector

Return the conjugate label `conj(a)`.
"""
dual(a::Sector) = conj(a)
Base.conj(::Trivial) = Trivial()

"""
    isreal(::Type{<:Sector}) -> Bool

Return whether the topological data (Fsymbol, Rsymbol) of the sector is real or not (in
which case it is complex).
"""
Base.@pure function Base.isreal(G::Type{<:Sector})
    u = one(G)
    return (eltype(Fsymbol(u,u,u,u,u,u))<:Real) && (eltype(Rsymbol(u,u,u))<:Real)
end
Base.@pure Base.isreal(::Type{Trivial}) = true

Base.isless(::Trivial, ::Trivial) = false

# FusionStyle: the most important aspect of Sector
#---------------------------------------------
"""
    ⊗(a::G, b::G) where {G<:Sector}

Return an iterable of elements of `c::G` that appear in the fusion product `a ⊗ b`.

Note that every element `c` should appear at most once, fusion degeneracies (if
`FusionStyle(G) == DegenerateNonAbelian()`) should be accessed via `Nsymbol(a,b,c)`.
"""
⊗(::Trivial, ::Trivial) = (Trivial(),)

"""
    Nsymbol(a::G, b::G, c::G) where {G<:Sector} -> Integer

Return an `Integer` representing the number of times `c` appears in the fusion product
`a ⊗ b`. Could be a `Bool` if `FusionStyle(G) == Abelian()` or `SimpleNonAbelian()`.
"""
function Nsymbol end
Nsymbol(::Trivial, ::Trivial, ::Trivial) = true

# trait to describe the fusion of superselection sectors
abstract type FusionStyle end
struct Abelian <: FusionStyle
end
abstract type NonAbelian <: FusionStyle end
struct SimpleNonAbelian <: NonAbelian # non-abelian fusion but multiplicity free
end
struct DegenerateNonAbelian <: NonAbelian # non-abelian fusion with multiplicities
end

"""
    FusionStyle(a::Sector) -> ::FusionStyle
    FusionStyle(G::Type{<:Sector}) -> ::FusionStyle

Return the type of fusion behavior of sectors of type G, which can be either
*   `Abelian()`: single fusion output when fusing two sectors;
*   `SimpleNonAbelian()`: multiple outputs, but every output occurs at most one,
    also known as multiplicity free (e.g. irreps of ``SU(2)``);
*   `DegenerateNonAbelian()`: multiple outputs that can occur more than once (e.g. irreps
    of ``SU(3)``).
There is an abstract supertype `NonAbelian` of which both `SimpleNonAbelian` and
`DegenerateNonAbelian` are subtypes.
"""
FusionStyle(a::Sector) = FusionStyle(typeof(a))
FusionStyle(::Type{Trivial}) = Abelian()

# NOTE: the following inline is extremely important for performance, especially
# in the case of Abelian, because ⊗(...) is computed very often
@inline function ⊗(a::G, b::G, c::G, rest::Vararg{G}) where {G<:Sector}
    if FusionStyle(G) isa Abelian
        return a ⊗ first(⊗(b, c, rest...))
    else
        s = Set{G}()
        for d in ⊗(b, c, rest...)
            for e in a ⊗ d
                push!(s, e)
            end
        end
        return s
    end
end

"""
    Fsymbol(a::G, b::G, c::G, d::G, e::G, f::G) where {G<:Sector}

Return the F-symbol ``F^{abc}_d`` that associates the two different fusion orders of sectors
`a`, `b` and `c` into an ouput sector `d`, using either an intermediate sector ``a ⊗ b → e``
or ``b ⊗ c → f``:
```
a-<-μ-<-e-<-ν-<-d                                     a-<-λ-<-d
    ∨       ∨       -> Fsymbol(a,b,c,d,e,f)[μ,ν,κ,λ]      ∨
    b       c                                         b-<-κ
                                                          ∨
                                                          c
```
If `FusionStyle(G)` is `Abelian` or `SimpleNonAbelian`, the F-symbol is a number. Otherwise
it is a rank 4 array of size
`(Nsymbol(a,b,e), Nsymbol(e,c,d), Nsymbol(b,c,f), Nsymbol(a,f,d))`.
"""
function Fsymbol end
Fsymbol(::Trivial, ::Trivial, ::Trivial, ::Trivial, ::Trivial, ::Trivial) = 1

"""
    Rsymbol(a::G, b::G, c::G) where {G<:Sector}

Returns the R-symbol ``R^{ab}_c`` that maps between ``a ⊗ b → c`` and ``b ⊗ a → c`` as in
```
a -<-μ-<- c                                 b -<-ν-<- c
     ∨          -> Rsymbol(a,b,c)[μ,ν]           ∧
     b                                           a
```
If `FusionStyle(G)` is `Abelian()` or `SimpleNonAbelian()`, the R-symbol is a number.
Otherwise it is a square matrix with row and column size `Nsymbol(a,b,c) == Nsymbol(b,a,c)`.
"""
function Rsymbol end
Rsymbol(::Trivial, ::Trivial, ::Trivial) = 1

# If a G::Sector with `fusion(G) == DegenerateNonAbelian` fusion wants to have custom vertex
# labels, a specialized method for `vertindex2label` should be added
"""
    vertex_ind2label(i::Int, a::G, b::G, c::G) where {G<:Sector}

Convert the index i of the fusion vertex (a,b)->c into a label. For
`FusionStyle(G) == Abelian()` or `FusionStyle(G) == NonAbelian()`, where every fusion
output occurs only once and `i == 1`, the default is to suppress vertex labels by setting
them equal to `nothing`. For `FusionStyle(G) == DegenerateNonAbelian()`, the default is to
just use `i`, unless a specialized method is provided.
"""
vertex_ind2label(i::Int, s1::G, s2::G, sout::G) where {G<:Sector}=
    _ind2label(FusionStyle(G), i::Int, s1::G, s2::G, sout::G)
_ind2label(::Abelian, i, s1, s2, sout) = nothing
_ind2label(::SimpleNonAbelian, i, s1, s2, sout) = nothing
_ind2label(::DegenerateNonAbelian, i, s1, s2, sout) = i

"""
    vertex_labeltype(G::Type{<:Sector}) -> Type

Return the type of labels for the fusion vertices of sectors of type `G`.
"""
Base.@pure vertex_labeltype(G::Type{<:Sector}) =
    typeof(vertex_ind2label(1, one(G), one(G), one(G)))

# combine fusion properties of tensor products of sectors
Base.:&(f::F, ::F) where {F<:FusionStyle} = f
Base.:&(f1::FusionStyle, f2::FusionStyle) = f2 & f1

Base.:&(::SimpleNonAbelian, ::Abelian) = SimpleNonAbelian()
Base.:&(::DegenerateNonAbelian, ::Abelian) = DegenerateNonAbelian()
Base.:&(::DegenerateNonAbelian, ::SimpleNonAbelian) = DegenerateNonAbelian()

# properties that can be determined in terms of the F symbol
# TODO: find mechanism for returning these numbers with custom type T<:AbstractFloat
"""
    dim(a::Sector)

Return the (quantum) dimension of the sector `a`.
"""
function dim(a::Sector)
    if FusionStyle(a) isa Abelian
        1
    elseif FusionStyle(a) isa SimpleNonAbelian
        abs(1/Fsymbol(a,conj(a),a,a,one(a),one(a)))
    else
        abs(1/Fsymbol(a,conj(a),a,a,one(a),one(a))[1])
    end
end

"""
    frobeniusschur(a::Sector)

Return the Frobenius-Schur indicator of a sector `a`.
"""
function frobeniusschur(a::Sector)
    if FusionStyle(a) isa Abelian || FusionStyle(a) isa SimpleNonAbelian
        sign(Fsymbol(a,conj(a),a,a,one(a),one(a)))
    else
        sign(Fsymbol(a,conj(a),a,a,one(a),one(a))[1])
    end
end

"""
    twist(a::Sector)

Return the twist of a sector `a`
"""
function twist(a::Sector)
    if FusionStyle(a) isa Abelian || FusionStyle(a) isa SimpleNonAbelian
        θ = sum(dim(b)/dim(a)*Rsymbol(a,a,b) for b in a ⊗ a)
    else
        # TODO: is this correct?
        # θ = sum(dim(b)/dim(a)*tr(Rsymbol(a,a,b)) for b in a ⊗ a)
        throw(MethodError(twist, (a,)))
    end
    return θ
end

"""
    Bsymbol(a::G, b::G, c::G) where {G<:Sector}

Return the value of ``B^{ab}_c`` which appears in transforming a splitting vertex
into a fusion vertex using the transformation
```
a -<-μ-<- c                                                    a -<-ν-<- c
     ∨          -> √(dim(c)/dim(a)) * Bsymbol(a,b,c)[μ,ν]           ∧
     b                                                            dual(b)
```
If `FusionStyle(G)` is `Abelian()` or `SimpleNonAbelian()`, the B-symbol is a number.
Otherwise it is a square matrix with row and column size
`Nsymbol(a, b, c) == Nsymbol(c, dual(b), a)`.
"""
function Bsymbol(a::G, b::G, c::G) where {G<:Sector}
    if FusionStyle(G) isa Abelian || FusionStyle(G) isa SimpleNonAbelian
        sqrt(dim(a)*dim(b)/dim(c))*Fsymbol(a, b, dual(b), a, c, one(a))
    else
        reshape(sqrt(dim(a)*dim(b)/dim(c))*Fsymbol(a,b,dual(b),a,c,one(a)),
            (Nsymbol(a,b,c), Nsymbol(c,dual(b),a)))
    end
end

# Not necessary
function Asymbol(a::G, b::G, c::G) where {G<:Sector}
    if FusionStyle(G) isa Abelian || FusionStyle(G) isa SimpleNonAbelian
        sqrt(dim(a)*dim(b)/dim(c))*conj(frobeniusschur(a)*Fsymbol(dual(a),a,b,b,one(a),c))
    else
        reshape(sqrt(dim(a)*dim(b)/dim(c))*
                conj(frobeniusschur(a)*Fsymbol(dual(a),a,b,b,one(a),c)),
                (Nsymbol(a,b,c), Nsymbol(dual(a),c,b)))
    end
end

# Braiding:
#-------------------------------------------------
# trait to describe type to denote how the elementary spaces in a tensor product space
# interact under permutations or actions of the braid group
abstract type BraidingStyle end # generic braiding
abstract type SymmetricBraiding <: BraidingStyle end # symmetric braiding => actions of permutation group are well defined
struct Bosonic <: SymmetricBraiding end # trivial under permutations
struct Fermionic <: SymmetricBraiding end
struct Anyonic <: BraidingStyle end

Base.:&(b::B,::B) where {B<:BraidingStyle} = b
Base.:&(B1::BraidingStyle, B2::BraidingStyle) = B2 & B1

Base.:&(::Bosonic,::Fermionic) = Fermionic()
Base.:&(::Bosonic,::Anyonic) = Anyonic()
Base.:&(::Fermionic,::Anyonic) = Anyonic()

"""
    BraidingStyle(::Sector) -> ::BraidingStyle
    BraidingStyle(G::Type{<:Sector}) -> ::BraidingStyle

Return the type of braiding and twist behavior of sectors of type `G`, which can be either
*   `Bosonic()`: symmetric braiding with trivial twist (i.e. identity)
*   `Fermionic()`: symmetric braiding with non-trivial twist (squares to identity)
*   `Anyonic()`: general ``R_(a,b)^c`` phase or matrix (depending on `SimpleNonAbelian` or
    `DegenerateNonAbelian` fusion) and arbitrary twists

Note that `Bosonic` and `Fermionic` are subtypes of `SymmetricBraiding`, which means that
braids are in fact equivalent to crossings (i.e. braiding twice is an identity:
`Rsymbol(b,a,c)*Rsymbol(a,b,c) = I`) and permutations are uniquely defined.
"""
BraidingStyle(a::Sector) = BraidingStyle(typeof(a))
BraidingStyle(::Type{Trivial}) = Bosonic()

# SectorSet:
#-------------------------------------------------------------------------------
# Custum generator to represent sets of sectors with type inference
struct SectorSet{G<:Sector,F,S}
    f::F
    set::S
end
SectorSet{G}(::Type{F}, set::S) where {G<:Sector,F,S} = SectorSet{G,Type{F},S}(F, set)
SectorSet{G}(f::F, set::S) where {G<:Sector,F,S} = SectorSet{G,F,S}(f, set)
SectorSet{G}(set) where {G<:Sector} = SectorSet{G}(identity, set)

Base.IteratorEltype(::Type{<:SectorSet}) = HasEltype()
Base.IteratorSize(::Type{SectorSet{G,F,S}}) where {G<:Sector,F,S} = Base.IteratorSize(S)

Base.eltype(::SectorSet{G}) where {G<:Sector} = G
Base.length(s::SectorSet) = length(s.set)
Base.size(s::SectorSet) = size(s.set)

function Base.iterate(s::SectorSet{G}, args...) where {G<:Sector}
    next = iterate(s.set, args...)
    next === nothing && return nothing
    val, state = next
    return convert(G, s.f(val)), state
end

# possible sectors
include("irreps.jl") # irreps of symmetry groups, with bosonic braiding
# include("fermions.jl") # irreps with defined fermionparity and fermionic braiding
include("anyons.jl") # non-group sectors
include("product.jl") # direct product of different sectors

ext: jl


text:
using MPSKit,TensorKit,DelimitedFiles

let
    #defining the hamiltonian
    th = nonsym_ising_ham(lambda = 0.3)
    (sxt,syt,szt) = nonsym_spintensors(1//2);

    #Center gauging a random mps
    ts=InfiniteMPS([ℂ^2],[ℂ^12]);

    #Finding the groundstate
    (ts,pars,_)=find_groundstate(ts,th,Vumps(maxiter=400));

    len=20;deltat=0.05;totaltime=3.0

    #apply a single spinflip at the middle site
    mpco = MPSComoving(ts,copy.([ts.AC[1];ts.AR[2:len]]),ts)
    @tensor mpco.AC[Int(round(len/2))][-1 -2;-3]:=ts.AC[Int(round(len/2))][-1,1,-3]*sxt[-2,1]

    pars = params(mpco,th)

    szdat=[expectation_value(mpco,szt)]

    for i in 1:(totaltime/deltat)
        (mpco,pars) = changebonds(mpco,th,OptimalExpand()&SvdCut(trschemes=[truncdim(20)])) # grow the bond dimension by 1, and truncate at bond dimension 20
        (mpco,pars) = timestep(mpco,th,deltat,Tdvp(),pars)
        push!(szdat,expectation_value(mpco,szt))
    end

    szdat=real.(reduce((a,b)->[a b],szdat))
    #writedlm("spinflip_evolution.csv",szdat)
end


text:
function myf4(;bin_data=nothing, arg1=10, arg2=31.4, int32=nothing)
   if !isnothing(bin_data)
      bin_lng = length(bin_data)
   else
      bin_lng = 0
   end
   a1 = arg1 * 10
   a2 = arg2 / 10.0
   return (; bin_lng, a1, a2)
end


text:
using ITensors
using ITensorMPOCompression

N=5 #5 sites
NNN=2 #Include 2nd nearest neighbour interactions
sites = siteinds("S=1/2",N);
H=make_transIsing_MPO(sites,NNN);
pprint(H[2]) #H[1] is a row vector, so let's see what H[2] looks like
Q,L,iq=block_qx(H[2]); #Block respecting QL
pprint(Q)
pprint(L,iq) #we need to tell pprint which index is the row index.


text:
using Heuristics
using BenchmarkTools

# Schwefel Function
function schaffer(x)
    obj = 0.5 + (sin(x[1]^2 + x[2]^2)^2 - 0.5)/(1 + 0.001*(x[1]^2+x[2]^2))^2
    return obj
end

function waveDrop(x)
    obj = -(1 + cos(12*sqrt(x[1]^2 + x[2]^2)))/(0.5*(x[1]^2 + x[2]^2) + 2.0)
    return obj
end

# Setup Problem
d = 2
LB = -5.12*ones(d)
UB = 5.12*ones(d)

prob = Problem(waveDrop, LB, UB)
pso  = PSO(prob; numParticles = 100)

# optimize
opts = Options(;display = false, maxStallIters = 100)
@benchmark optimize!($pso, $opts)
#optimize!(pso, opts)


text:
const jobs = 12

mutable struct Environment
    seq::Int
    cnt::Int
    Environment() = new(0, 0)
end

const env = [Environment() for i in 1:jobs]
const currentjob = [1]

seq() = env[currentjob[1]].seq
cnt() = env[currentjob[1]].cnt
seq(n) = (env[currentjob[1]].seq = n)
cnt(n) = (env[currentjob[1]].cnt = n)

function hail()
    print(lpad(seq(), 4))
    if seq() == 1
        return
    end
    cnt(cnt() + 1)
    seq(isodd(seq()) ? 3 * seq() + 1 : div(seq(), 2))
end

function runtest()
    for i in 1:jobs
        currentjob[1] = i
        env[i].seq = i
    end
    computing = true
    while computing
        for i in 1:jobs
            currentjob[1] = i
            hail()
        end
        println()
        for j in 1:jobs
            currentjob[1] = j
            if seq() != 1
                break
            elseif j == jobs
                computing = false
            end
        end
    end
    println("\nCOUNTS:")
    for i in 1:jobs
        currentjob[1] = i
        print(lpad(cnt(), 4))
    end
    println()
end

runtest()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
export CuModule, CuFunction

type CuModule
    ptr::Ptr{Void}

    function CuModule(ptr)
        m = new(ptr)
        finalizer(m, cuModuleUnload)
        m
    end
end

Base.unsafe_convert(::Type{Ptr{Void}}, m::CuModule) = m.ptr

function CuModule(image::Vector{UInt8})
    p = Ptr{Void}[0]
    cuModuleLoadData(p, image)
    #cuModuleLoadDataEx(p, image, 0, CUjit_option[], Ptr{Void}[])
    CuModule(p[1])
end

type CuFunction
    m::CuModule # avoid CuModule gc-ed
    ptr::Ptr{Void}
end

Base.unsafe_convert(::Type{Ptr{Void}}, f::CuFunction) = f.ptr

function CuFunction(m::CuModule, name::String)
    p = CUfunction[0]
    cuModuleGetFunction(p, m, name)
    CuFunction(m, p[1])
end

function CuFunction(code::String)
    code = replace(code, "Float32", "float")
    code = replace(code, "Float64", "double")
    code = replace(code, "Int32", "int")
    contains(code, "Array<") && (code = "$(Interop.array_h)\n$code")
    #contains(code, "Ranges<") && (code = "$range_h\n$code")

    ptx = NVRTC.compile(code)
    p = Ptr{Void}[0]
    cuModuleLoadData(p, pointer(ptx))
    mod = CuModule(p[1])

    fnames = []
    for line in split(ptx,'\n')
        m = match(r".visible .entry (.+)\(", line) # find function name
        m == nothing && continue
        push!(fnames, String(m[1]))
    end
    length(fnames) > 1 && throw("Multiple functions are found.")
    CuFunction(mod, fnames[1])
end

box(x) = x
box(x::Int) = Cint(x)
box{N}(t::NTuple{N,Int}) = map(Cint, t)
box(x::Vector{Int}) = ntuple(i -> Cint(x[i]), length(x))

function (f::CuFunction)(args...;
    dx=1, dy=1, dz=1, bx=128, by=1, bz=1, sharedmem=0, stream=C_NULL)
    argptrs = Ptr{Void}[pointer_from_objref(box(a)) for a in args]
    gx = ceil(dx / bx)
    gy = ceil(dy / by)
    gz = ceil(dz / bz)
    cuLaunchKernel(f, gx, gy, gz, bx, by, bz, sharedmem, stream, argptrs, stream)
end

#=
macro compile(expr)
    expr.head == :string || throw("expr is not string")
    idx = findfirst(expr.args) do a
        isa(a,String) && match(r"__global__ void", a) != nothing
    end
    idx == 0 && throw("Cannot find \"__global__ void\".")
    dict = Dict()
    for i = idx+1:length(expr.args)
        arg = expr.args[i]
        isa(arg, Symbol) && (dict[arg] = arg)
        isa(arg, String) && contains(arg, ")") && break # end of function declaration
    end
    syms = Expr(:tuple, keys(dict)...)
    dict = Dict()

    quote
        local dict = $dict
        local key = $(esc(syms))
        get!(dict, key) do
            local code = $(esc(expr))
            local ptx = NVRTC.compile(code)
            p = Ptr{Void}[0]
            cuModuleLoadData(p, pointer(ptx))
            mod = CuModule(p[1])

            # TODO: multi-device
            for line in split(ptx,'\n')
                m = match(r".visible .entry (.+)\(", line) # find function name
                m == nothing && continue
                fname = Symbol(m[1])
                return CuFunction(mod, string(fname))
            end
        end
    end
end
=#


text:
tropomi_tle_string = """
       S5P
       1 42969U 17064A   21193.45435625 -.00000014  00000-0  14036-4 0  9992
       2 42969  98.7381 133.0236 0001130  90.2950 269.8356 14.19543373194087
       """

function CO2M_uncertainty_function(TOA)

    return 0.068055560 + 0.0017075595 * TOA - 1.0989804e-5 * TOA^2 + 2.5627241e-8 * TOA^3

end

function intersect_with_wgs84!(
    origin,
    look_vector,
    output,
    a=6378137.0, # m
    b=6356752.314245 # m
    )

    x, y, z = origin
    u, v, w = look_vector

    t = -(1.0/(b^2 * (u^2 + v^2) + (a^2 * w^2))) * (b^2 * (u*x + v*y) +
        (a^2 * w*z) + 0.5*sqrt(4 * (b^2 * (u*x + v*y) + (a^2 * w*z))^2 -
        4 * (b^2 * (u^2 + v^2) + (a^2 * w^2)) * (b^2 * (-a^2 + x^2 + y^2) + (a^2 * z^2))))

    output[1] = x + t * u
    output[2] = y + t * v
    output[3] = z + t * w

end

function NadirSwathOrbiterSampling(
    target_lon::Real,
    target_lat::Real,
    radius::Real,
    start_time::DateTime,
    end_time::DateTime,
    tle_string_OR_kepler,
    field_of_view::Real,
    N_pixel::Integer,
    exposure_time::Real,
    uncertainty_function::Function;
    instrument::String="instrument",
    return_only_locations::Bool=false,
    )

    locarray = Geolocation[]
    loctimearray = GeolocationTime[]
    scenearray = Scene[]

    #wgs84_a = 6378137.0 # m
    #wgs84_b = 6356752.314245 # m

    #wgs84_a2 = wgs84_a^2
    #wgs84_b2 = wgs84_b^2

    # This is needed for whatever reason?
    eop = get_iers_eop()

    # Overall strategy:
    # First produce orbit, propagate orbit and generate
    # sounding locations, after which we can kick out
    # any scenes that are nighttime etc.

    sat_lons = Float32[]
    sat_lats = Float32[]

    sounding_lons = Float32[]
    sounding_lats = Float32[]

    # Check the type of
    if tle_string_OR_kepler isa String

        tle_string = tle_string_OR_kepler
        tle = read_tle_from_string(tle_string)
        orbp = init_orbit_propagator(Val(:sgp4), tle[1])

    elseif tle_string_OR_kepler isa KeplerianElements

        k = tle_string_OR_kepler
        orbp = init_orbit_propagator(Val(:J4),
                                     k.t, # epoch
                                     k.a, # semi-major axis [m]
                                     k.e, # eccentricity
                                     k.i, # inclination [rad]
                                     k.Ω, # right ascension of ascending node [rad]
                                     k.ω, # argument of perigee / periapsis [rad],
                                     k.f  # true anomaly [rad]
                                     )

    end

    # Array of angles (swath angles away from full nadir)
    theta_array_deg = collect(-N_pixel//2:1:N_pixel//2) .* (field_of_view / N_pixel) # TROPOMI: 216 pixels with 0.5 deg each
    theta_array = deg2rad.(theta_array_deg)
    cos_theta_array = cos.(theta_array)
    sin_theta_array = sin.(theta_array)

    # Rotation around x, y and z
    Rx = Matrix(1.0I, 3, 3)
    Ry = Matrix(1.0I, 3, 3)
    Rz = Matrix(1.0I, 3, 3)

    # temp matrices
    tmp1 = Matrix(1.0I, 3, 3)
    tmp2 = Matrix(1.0I, 3, 3)
    tmp3 = ones(3)
    tmp4 = ones(3)
    look_vector = ones(3)

    r_intersection = ones(3)

    # First - construct a sensible bounding box to load the surface data
    wgs84 = Projection(Proj4.epsg[4326])
    xy = lonlat2xy([Float64(target_lon), Float64(target_lat)], wgs84)
    r = radius * 1000.0 * 1.10 # Add x percent as buffer, radius is given in [km] but we need [m]

    # If the user only wants location arrays, we can skip this portion
    # needed for the construction of SIF scenes.
    if ~return_only_locations

        north = xy2lonlat(geod_destination(xy, 0.0, r, wgs84), wgs84)
        east = xy2lonlat(geod_destination(xy, 90.0, r, wgs84), wgs84)
        south = xy2lonlat(geod_destination(xy, 180.0, r, wgs84), wgs84)
        west = xy2lonlat(geod_destination(xy, 270.0, r, wgs84), wgs84)

        lon_min = minimum([west[1], south[1]])
        lat_min = minimum([west[2], south[2]])
        lon_max = maximum([east[1], north[1]])
        lat_max = maximum([east[2], north[2]])

        vnp_sd = create_VNP_SD_from_locbounds(
            lon_min,
            lat_min,
            lon_max,
            lat_max
        )

    end

    # this does not create a copy, however if you add a time period
    # it *will* create a new object, such that "start_time" will not
    # be modified.
    current_time = start_time

    pasted_swath_width = false

    while current_time <= end_time

        # Increment by exposure time
        current_time += Dates.Millisecond(exposure_time * 1000.0)

        # Turn the current time into a julian date
        current_epoch = date_to_jd(
            Dates.year(current_time),
            Dates.month(current_time),
            Dates.day(current_time),
            Dates.hour(current_time),
            Dates.minute(current_time),
            Dates.second(current_time) + Dates.millisecond(current_time) / 1000.0
        )

        # Propagate the orbit to the requested date
        r_teme, v_teme = propagate_to_epoch!(orbp, current_epoch)
        D_ITRF_TEME = r_eci_to_ecef(TEME(), ITRF(), current_epoch, eop)

        r_itrf = D_ITRF_TEME * r_teme
        x, y, z = r_itrf
        x2 = x*x
        y2 = y*y
        z2 = z*z

        v_itrf = D_ITRF_TEME * v_teme

        # Calculate the off-nadir footprints from the current
        # satellite position!
        this_lat, this_lon, this_alt = ecef_to_geodetic(r_itrf)
        sat_lon = rad2deg(this_lon)
        sat_lat = rad2deg(this_lat)

        # Normalized vector pointing from the satellite down
        # to center of the Earth
        r_sat_to_origin = normalize(-r_itrf)

        # Normalized velocity vector gives us flight direction
        v_norm = normalize(v_itrf)

        # Compute the rotation matrices for this point in time
        u, v, w = r_sat_to_origin # The vector to be rotated around rotation axis
        a, b, c = v_norm # The rotation axis
        d = sqrt(b^2 + c^2)

        # Rotations around x and y are independent of the swath angle
        Rx[2,2] = c / d
        Rx[3,3] = c / d
        Rx[2,3] = -b / d
        Rx[3,2] = b / d

        Ry[1,1] = d
        Ry[3,3] = d
        Ry[3,1] = a
        Ry[1,3] = -a

        @views tmp1[:,:] = Rx' * Ry'
        #=
        tmp1[1,1] = d
        tmp1[1,3] = a
        tmp1[2,1] = -a*b/d
        tmp1[2,2] = c/d
        tmp1[2,3] = b
        tmp1[3,1] = -a*c/d
        tmp1[3,2] = -b/d
        tmp1[3,3] = c
        =#

        @views tmp2[:,:] = Ry * Rx
        #=
        @views tmp2[:,:] = tmp1[:,:]
        tmp2[3,1] *= -1
        tmp2[3,2] *= -1
        tmp2[1,3] *= -1
        tmp2[2,3] *= -1
        =#

        # tmp3 = Ry * Rx * r_sat_to_origin
        @views tmp3[:] = tmp2 * r_sat_to_origin
        #=
        tmp3[1] = u*d - a*b*v/d - a*c*w/d #u*d - a*w
        tmp3[2] = c*v/d - b*w/d #-a*b*u/d + c*v/d - b*w
        tmp3[3] = a*u + b*v + c*w #a*c*u/d + b*v/d + c*w
        =#

        # Some smart guessing method here, to determine whether we actually
        # need to calculate the sounding locations of the entire swath - or
        # whether we can just move on to the next frame.

        # If the distance between ROI center and subsattelite point is larger than twice the
        # swath extremal point distance, skip this!

        swath_end_lons = zeros(2)
        swath_end_lats = zeros(2)

        for i in 1:2

            if i == 1
                idx = 1
            else
                idx = length(theta_array)
            end

            Rz[1,1] = cos_theta_array[idx]
            Rz[2,2] = cos_theta_array[idx]
            Rz[1,2] = -sin_theta_array[idx]
            Rz[2,1] = sin_theta_array[idx]

            #tmp4 = Rz * tmp3
            tmp4[1] = Rz[1,1] * tmp3[1] + Rz[1,2] * tmp3[2]
            tmp4[2] = Rz[2,1] * tmp3[1] + Rz[2,2] * tmp3[2]
            tmp4[3] = tmp3[3]

            look_vector[1] = tmp1[1,1] * tmp4[1] + tmp1[1,2] * tmp4[2] + tmp1[1,3] * tmp4[3]
            look_vector[2] = tmp1[2,1] * tmp4[1] + tmp1[2,2] * tmp4[2] + tmp1[2,3] * tmp4[3]
            look_vector[3] = tmp1[3,1] * tmp4[1] + tmp1[3,2] * tmp4[2] + tmp1[3,3] * tmp4[3]

            # Intersect look vector with WGS84 geoid
            intersect_with_wgs84!(r_itrf, look_vector, r_intersection)
            this_lat, this_lon, this_altitude = ecef_to_geodetic(r_intersection)

            swath_end_lons[i] = rad2deg(this_lon)
            swath_end_lats[i] = rad2deg(this_lat)

        end

        # Distance between first and last footprint in a swath
        distance_swath = calculate_distance(swath_end_lons[1], swath_end_lats[1],
                                            swath_end_lons[2], swath_end_lats[2])

        if ~pasted_swath_width
            println("Swath width: $(distance_swath)")
            pasted_swath_width = true
        end

        # Distance between subsatellite point and target ROI center
        distance_ss_to_ROI_center = calculate_distance(
            sat_lon, sat_lat,
            target_lon, target_lat
        )

        # If our target location is too far away from this particular swath, just
        # move on to the next one. If the user wants locations, then skip this check
        # as well.
        if (distance_ss_to_ROI_center - radius > distance_swath) & ~return_only_locations
            continue
        end

        # Process the full swath
        for i in 1:length(theta_array)

            Rz[1,1] = cos_theta_array[i]
            Rz[2,2] = cos_theta_array[i]
            Rz[1,2] = -sin_theta_array[i]
            Rz[2,1] = sin_theta_array[i]

            #tmp4 = Rz * tmp3
            tmp4[1] = Rz[1,1] * tmp3[1] + Rz[1,2] * tmp3[2]
            tmp4[2] = Rz[2,1] * tmp3[1] + Rz[2,2] * tmp3[2]
            tmp4[3] = tmp3[3]

            look_vector[1] = tmp1[1,1] * tmp4[1] + tmp1[1,2] * tmp4[2] + tmp1[1,3] * tmp4[3]
            look_vector[2] = tmp1[2,1] * tmp4[1] + tmp1[2,2] * tmp4[2] + tmp1[2,3] * tmp4[3]
            look_vector[3] = tmp1[3,1] * tmp4[1] + tmp1[3,2] * tmp4[2] + tmp1[3,3] * tmp4[3]

            # @views look_vector[:] = Rx' * Ry' * Rz * Ry * Rx * r_sat_to_origin
            # @views look_vector[:] = tmp1 * tmp4

            # Intersect look vector with WGS84 geoid
            intersect_with_wgs84!(r_itrf, look_vector, r_intersection)
            this_lat, this_lon, this_altitude = ecef_to_geodetic(r_intersection)

            this_loc = Geolocation(rad2deg(this_lon), rad2deg(this_lat))

            if return_only_locations
                # If the user only wants locations, push the current location into the array
                # and move on.
                this_loctime = GeolocationTime(this_loc, current_time)
                push!(loctimearray, this_loctime)

            else

                # If the user wants the full SIF scene, do necessary calculations here.

                # This footprint is indeed within our requested target
                if check_location_within_radius(target_lon, target_lat, radius, this_loc)

                    this_loctime = GeolocationTime(this_loc, current_time)

                    # solar azimuth (saa) at this point unused!
                    this_sza, this_saa = calculate_solar_angles(this_loctime)

                    # If we are looking at nighttime - skip
                    if isnan(this_sza) | (this_sza > 90.0)
                        continue
                    end

                    # Calculate viewing zenith angles
                    # This should be the angle between the vector pointing
                    # from measurement location (r_intersections) and the
                    # satellite position (r_itrf) in ECEF, and the measurement location
                    # and the normal (going from earth center to measurement location)

                    r_norm = normalize(r_intersection)
                    r_loc_to_sat_norm = normalize(r_itrf - r_intersection)

                    _tmp = dot(r_loc_to_sat_norm, r_norm)
                    if (_tmp > 1) & (_tmp < 1 + 1e-6)
                        _tmp = 1.0
                    end

                    this_vza = rad2deg(acos(_tmp))

                    this_scene = create_SIF_scene(
                        instrument,
                        "nadir",
                        this_loctime,
                        this_vza,
                        uncertainty_function,
                        vnp_sd
                    )

                    push!(locarray, this_loc)
                    push!(scenearray, this_scene)

                end

            end

        end

        # Compute the new look vectors going from center of spacecraft
        # in ECEF frame.
        #new_look_vectors = rotate_vec_arbitrary.(
        #    Ref(r_sat_to_origin), # point to be rotated
        #    Ref(r_itrf), # position of rotation origin
        #    Ref(v_norm), # rotation axis relative to origin
        #    theta_array) # rotation angle

    end

    if return_only_locations

        return loctimearray

    else

        return NadirSwathOrbiterSampling(
            "",
            ["instrument"],
            locarray,
            scenearray
        )

    end

end


text:
module Normalize

using Optim: optimize, minimizer
using Statistics: mean, var
using StatsBase: geomean, Histogram
using Distributions: fit, Normal, pdf

include("types.jl")
include("chisq.jl")
include("boxcox.jl")
include("yeojohnson.jl")

end


text:
module ChainRulesTest

using LinearAlgebra
using ChainRulesCore
using DiffResults
using ReverseDiff
using Test

f(x) = sum(4x .+ 1)

function ChainRulesCore.rrule(::typeof(f), x)
    r = f(x)
    function back(d)
        #=
        The proper derivative of `f` is 4, but in order to
        check if `ChainRulesCore.rrule` had taken over the compuation,
        we define a rrule that returns 3 as `f`'s derivative.

        After importing this rrule into ReverseDiff, if we get 3
        rather than 4 when we compute the derivative of `f`, it means
        the importing mechanism works.
        =#
        return ChainRulesCore.NoTangent(), fill(3 * d, size(x))
    end
    return r, back
end

ReverseDiff.@grad_from_chainrules f(x::ReverseDiff.TrackedArray)

g(x, y) = sum(4x .+ 4y)

function ChainRulesCore.rrule(::typeof(g), x, y)
    r = g(x, y)
    function back(d)
        # same as above, use 3 and 5 as the derivatives
        return ChainRulesCore.NoTangent(), fill(3 * d, size(x)), fill(5 * d, size(x))
    end
    return r, back
end

ReverseDiff.@grad_from_chainrules g(x::ReverseDiff.TrackedArray, y)
ReverseDiff.@grad_from_chainrules g(x, y::ReverseDiff.TrackedArray)
ReverseDiff.@grad_from_chainrules g(x::ReverseDiff.TrackedArray, y::ReverseDiff.TrackedArray)

@testset "rrule in ChainRules and ReverseDiff" begin
    ## ChainRules
    # function f
    input = rand(3, 3)
    output, back = ChainRulesCore.rrule(f, input);
    _, d = back(1)
    @test output == f(input)
    @test d == fill(3, size(input))
    # function g
    inputs = rand(3, 3), rand(3, 3)
    output, back = ChainRulesCore.rrule(g, inputs...);
    _, d1, d2 = back(1)
    @test output == g(inputs...)
    @test d1 == fill(3, size(inputs[1]))
    @test d2 == fill(5, size(inputs[2]))

    ## ReverseDiff
    #function f
    inputs = (rand(3, 3), )

    results = (similar(inputs[1]),)
    f_tape = ReverseDiff.GradientTape(x -> f(x) + 2, (rand(3, 3),))
    ReverseDiff.gradient!(results, f_tape, inputs)

    @test results[1] == fill(3, size(inputs[1]))

    results = (similar(inputs[1]),)
    compiled_tape = ReverseDiff.CompiledTape(f_tape)
    ReverseDiff.gradient!(results, compiled_tape, inputs)
    @test results[1] == fill(3, size(inputs[1]))

    # function g
    inputs = rand(3, 3), rand(3, 3)

    results = (similar(inputs[1]), similar(inputs[2]))
    f_tape = ReverseDiff.GradientTape((x, y) -> g(x, y) + 2, (rand(3, 3), rand(3, 3)))
    ReverseDiff.gradient!(results, f_tape, inputs)

    @test results[1] == fill(3, size(inputs[1]))
    @test results[2] == fill(5, size(inputs[2]))

    results = (similar(inputs[1]), similar(inputs[2]),)
    compiled_tape = ReverseDiff.CompiledTape(f_tape)
    ReverseDiff.gradient!(results, compiled_tape, inputs)
    @test results[1] == fill(3, size(inputs[1]))
    @test results[2] == fill(5, size(inputs[2]))

end

### Tape test
@testset "Tape test: Ensure ordinary call is not tracked" begin
    tp = ReverseDiff.InstructionTape()

    f(x) = sum(2x .+ g([1, 2], [3, 4]))
    x = rand(3, 3)
    xt = ReverseDiff.track(copy(x), tp)
    # record
    yt = f(xt)
    @test length(tp) == 3 # sum, broadcast+, broadcast*, but not `g`
end

### Functions with varargs and kwargs
# Varargs
f_vararg(x, args...) = sum(4x .+ sum(args))

function ChainRulesCore.rrule(::typeof(f_vararg), x, args...)
    r = f_vararg(x, args...)
    function back(d)
        return ChainRulesCore.NoTangent(), fill(3 * d, size(x))
    end
    return r, back
end

ReverseDiff.@grad_from_chainrules f_vararg(x::ReverseDiff.TrackedArray, args...)

@testset "Function with Varargs" begin
    inputs = (rand(3, 3), )

    results = (similar(inputs[1]),)
    f_tape = ReverseDiff.GradientTape(x -> f_vararg(x, 1, 2, 3) + 2, (rand(3, 3),))
    ReverseDiff.gradient!(results, f_tape, inputs)

    @test results[1] == fill(3, size(inputs[1]))
end

# Vargs and kwargs
f_kw(x, args...; k=1, kwargs...) = sum(4x .+ sum(args) .+ (k + kwargs[:j]))

function ChainRulesCore.rrule(::typeof(f_kw), x, args...; k=1, kwargs...)
    r = f_kw(x, args...; k=k, kwargs...)
    function back(d)
        return ChainRulesCore.NoTangent(), fill(3 * d, size(x))
    end
    return r, back
end

ReverseDiff.@grad_from_chainrules f_kw(x::ReverseDiff.TrackedArray, args...; k=1, kwargs...)

@testset "Function with Varargs and kwargs" begin
    inputs = (rand(3, 3), )

    results = (similar(inputs[1]),)
    f_tape = ReverseDiff.GradientTape(x -> f_kw(x, 1, 2, 3; k=2, j=3) + 2, (rand(3, 3),))
    ReverseDiff.gradient!(results, f_tape, inputs)

    @test results[1] == fill(3, size(inputs[1]))
end

### Mix @grad and @grad_from_chainrules

h(x) = 10x
h(x::ReverseDiff.TrackedArray) = ReverseDiff.track(h, x)
ReverseDiff.@grad function h(x)
    xv = ReverseDiff.value(x)
    return h(xv), Δ -> (Δ * 7,) # use 7 asits derivatives
end

@testset "ReverseDiff and ChainRules Mixed" begin
    t(x) = g(x, h(x))
    inputs = (rand(3, 3), )
    results = (similar(inputs[1]),)

    g_tape = ReverseDiff.GradientTape(t, (rand(3, 3),))
    ReverseDiff.gradient!(results, g_tape, inputs)
    @test results[1] == fill(38, size(inputs[1])) # 38 = 3 + 5 * 7
end

### Isolated Scope
module IsolatedModuleForTestingScoping
using ChainRulesCore
using ReverseDiff: @grad_from_chainrules

f(x) = sum(4x .+ 1)

function ChainRulesCore.rrule(::typeof(f), x)
    r = f(x)
    function back(d)
        # return a distinguishable but improper grad
        return ChainRulesCore.NoTangent(), fill(3 * d, size(x))
    end
    return r, back
end

@grad_from_chainrules f(x::TrackedArray)

module SubModule
using Test
using ReverseDiff: TrackedArray, GradientTape, gradient!
using ..IsolatedModuleForTestingScoping: f
@testset "rrule in Isolated Scope" begin
    inputs = (rand(3, 3), )

    results = (similar(inputs[1]),)
    f_tape = GradientTape(x -> f(x) + 2, (rand(3, 3),))
    gradient!(results, f_tape, inputs)

    @test results[1] == fill(3, size(inputs[1]))
end

end # end of SubModule
end # end of IsolatedModuleForTestingScoping

end

max_stars_repo_name: jrevels/ReverseDiffPrototype

max_issues_repo_name: jrevels/ReverseDiffPrototype

max_forks_repo_name: jrevels/ReverseDiffPrototype


text:
# дз (!)
# 1. графики
# 2. базовые операции с данными
# 3. операции с несколькими табличками
# 4. группировка табличных данных
# 5. регрессия
# 6. PCA - SVD
# 7. гетероскедастичность
# предложить своё :)

versioninfo()

using RDatasets # импортировать обучающие данные из R
using Econometrics
using GLM
using CSV # чтение-запись .csv
using Gadfly # графики
using LinearAlgebra

# если чего-то из этого нет
# using Pkg; Pkg.add("GLM")
# using Pkg ~ import Pkg
# Pkg.add("GLM")

diamonds = dataset("ggplot2", "diamonds")

x = 5
x = 5;

pwd() # посмотреть
cd("/home/boris/Downloads/") # сменить

CSV.write("diam.csv", diamonds)
d2 = CSV.read("diam.csv")

# кусочек
diamonds[2:5, ["Cut", "Table"]]
diamonds[2:5, [:Cut, :Table]]

diamonds[!, ["Cut"]]

# столбец
diamonds.Cut
diamonds["Cut"]
diamonds[:Cut]
diamonds["Cut_copy"] = diamonds["Cut"]
diamonds[:Cut_copy2] = diamonds["Cut"]

diamonds.Cut_copy3 = diamonds["Cut"]

# описательное :)
describe(diamonds)
head(diamonds)
tail(diamonds)

# разница
[5, 6] == 5 # правда ли что вектор и число это одно и то же?
[5, 6] .== 5 # поточечно есть ли равенство?

# отбор строк

diamonds[diamonds.X > 6, ["Cut"]] # не сработает
diamonds[diamonds.X .> 6, ["Cut"]] # сработает (!)

# идея как сделать: robust_model = model(vce=HC1)
?model

model = fit(EconometricModel, @formula(Price ~ Carat + X + Y + Z), diamonds) # единичный столбец добавляется автоматически!
model_nointercept = fit(EconometricModel, @formula(Price ~ 0 + Carat + X + Y + Z), diamonds) # единичный столбец не добавляем

vcov(model) # ковариационная матрица обычная, неустойчивая к гетероскедастичности
# RSS/(n-k) * (X'X)^{-1}

vcov(model, HC1) # устойчивая к гетероскедастичности
vcov(model, HC0) # устойчивая к гетероскедастичности, оценки Уайта
vcov(model, HC3) # устойчивая к гетероскедастичности! выдаст ошибку

# \hat\beta = (X'X)^{-1}X'y
# \hat\Var\hat\beta = (X'X)^{-1}X' \hat \Sigma X(X'X)^{-1}
# [k x k] = [5 x 5]
# \hat \Sigma = [50000 x 50000] => out of memory

# извлечем что-то из модели :)
coef(model) # оценки коэффициентов
residuals(model) # остатки
fitted(model) # прогнозы
coeftable(model) # табличка
sqrt.(diag(vcov(model))) # стандартные ошибки коэффициентов
sqrt.([5, 6]) # поэлементный корень
sqrt([5, 6]) # не сработает!
sqrt(vcov(model)) # матричный корень!

# доверительные интервалы для коэффициентов неустойчивые к гетероскедастичности
confint(model, level=0.9)

# доверительные интервалы для коэффициентов устойчивые к гетероскедастичности
confint(model, level=0.9, se=sqrt.(diag(vcov(model, HC1))))

# графики
iris = dataset("datasets", "iris")
iris
plot(iris, x=:SepalLength, y=:PetalLength,
    Guide.xlabel("Длина чашелистика в мм"))

plot(iris, x=:SepalLength, Geom.bar,
    Guide.xlabel("Длина чашелистика в мм"))

# Gadfly

# создать новую переменную
diamonds["price2"] = diamonds[:Price] .^ 2
diamonds[:price3] = diamonds[:Price] .^ 2

# тест Уайта
# (!) не подвешивайте выбор робастные-обычные стандартные ошибки на результат теста Уайта
# (!) сразу используйте робастные стандартные ошибки
# (!) тест Уайта только для проверки гипотезы, есть гетероскедастиность или нет (!)
diamonds["resid2"] = residuals(model) .^ 2
aux_model = fit(EconometricModel, @formula(resid2 ~ Carat + X + Carat^2 + X^2 + X*Carat), diamonds)
# PATSY model syntax

# LR-test likelihood ratio test = тест отношения правдоподобия
# (H0, restricted) model_simple ∈ model_complex (Ha, unrestricted)
# LR = 2 * (ln L_ur - ln L_r) → χ^2 (d), где d — разница в числе параметров между simple и complex model
# H0: @formula(resid2 ~ 1) - нет гетероскедастичности
# Ha: @formula(resid2 ~ Carat + X + Carat^2 + X^2 + X*Carat) - есть гетероскедастичность
# H0 отвергается

β = 5
χ = 6
β = 9
β

# штрафные критерии
aic(model)
bic(model)

# McReel Econometrics in julia
# https://github.com/mcreel/Econometrics/
# quantecon julia
# https://julia.quantecon.org/
# MIT, Computational Thinking, 2020
# https://computationalthinking.mit.edu/Fall20/
# 3blue1brown


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "libunistring"
version = v"0.9.10"

# Collection of sources required to build libunistring
sources = [
    "http://ftp.gnu.org/gnu/libunistring/libunistring-0.9.10.tar.xz" =>
    "eb8fb2c3e4b6e2d336608377050892b54c3c983b646c561836550863003c05d7",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd libunistring-0.9.10/
./configure --prefix=$prefix --host=$target LDFLAGS=-L$prefix/lib CFLAGS=-I$prefix/include
make
make install
exit

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libunistring", :libunistring)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/Libiconv-v1.15-0/build_Libiconv.v1.15.0.jl"
]

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

max_stars_repo_name: yukota/libunistringBuilder

max_issues_repo_name: yukota/libunistringBuilder

max_forks_repo_name: yukota/libunistringBuilder


text:
module TestStateSpace
using CustomTest
using Base.Test
using ControlSystems
# Naming convention:
# ------------------
# {type}[S]_{dims}[_d][_n]
# type: C: Continuous, D: Discrete
# scalar: add "S" if any input is scalar
# dims: "nxnuny"
# feedthrough: append "_d" if `D` is present
# names: append "_n" if some inputs/outputs/states are named

# SCALARS
a_2 = [-5 -3; 2 -9]
CS_111 = ss(-5, 2, 3, [0])
CS_111_d = ss([3], 2, 1, 1)
CS_211 = ss(a_2, [1; 2], [1 0], 0)
CS_221 = ss(a_2, [1 0; 0 2], [1 0], 0)
CS_222 = ss(a_2, [1 0; 0 2], eye(2), 0)

# CONTINUOUS
a_1 = [-5]
C_111 = ss(a_1, [2], [3], [0])
C_211 = ss(a_2, [1; 2], [1 0], [0])
C_212 = ss(a_2, [1; 2], eye(2), [0; 0])
C_221 = ss(a_2, [1 0; 0 2], [1 0], [0 0])
C_222 = ss(a_2, [1 0; 0 2], eye(2), zeros(2,2))
C_222_d = ss(a_2, [1 0; 0 2], eye(2), eye(2))
C_022 = ss(4*eye(2))

# DISCRETE
da_1 = [-0.5]
da_2 = [0.2 -0.8; -0.8 0.07]
D_111 = ss(da_1, [2], [3], [0], 0.005)
D_211 = ss(da_2, [1; 2], [1 0], [0], 0.005)
D_221 = ss(da_2, [1 0; 0 2], [1 0], [0 0], 0.005)
D_222 = ss(da_2, [1 0; 0 2], eye(2), zeros(2,2), 0.005)
D_222_d = ss(da_2, [1 0; 0 2], eye(2), eye(2), 0.005)
D_022 = ss(4*eye(2), 0.005)

# Definition of input, output and state names
C_222_d_n = ss(a_2, [1 0; 0 2], eye(2), eye(2),
  statenames=["i","u"],inputnames=String("e"),outputnames="theta")

# TESTS
# Contstuct with scalars
@test CS_111 == C_111
@test CS_111_d == ss([3],[2],[1],[1])
@test CS_211 == C_211
@test CS_221 == C_221
@test CS_222 == C_222

# Addition
@test C_111 + C_111 == ss([-5 0; 0 -5],[2; 2],[3 3],[0])
@test C_222 + C_222 == ss([-5 -3 0 0; 2 -9 0 0; 0 0 -5 -3;
        0 0 2 -9],[1 0; 0 2; 1 0; 0 2], [1 0 1 0; 0 1 0 1],[0 0; 0 0])
@test C_222 + 1 == ss([-5 -3; 2 -9],[1 0; 0 2],[1 0; 0 1],[1 1; 1 1])
@test D_111 + D_111 == ss([-0.5 0; 0 -0.5],[2; 2],[3 3],[0], 0.005)

# Subtraction
@test C_111 - C_211 == ss([-5 0 0; 0 -5 -3; 0 2 -9],[2; 1; 2],[3 -1 -0],[0])
@test 1 - C_222 == ss([-5 -3; 2 -9],[1 0; 0 2],[-1 -0; -0 -1],[1 1; 1 1])
@test D_111 - D_211 == ss([-0.5 0 0; 0 0.2 -0.8; 0 -0.8 0.07],[2; 1; 2],
        [3 -1 -0],[0], 0.005)

# Multiplication
@test C_111 * C_221 == ss([-5 2 0; 0 -5 -3; 0 2 -9],
        [0 0; 1 0; 0 2],[3 0 0],[0 0])
@test C_212 * C_111 == ss([-5 -3 3; 2 -9 6; 0 0 -5],
        [0; 0; 2],[1 0 0; 0 1 0],[0; 0])
@test 4*C_222 == ss([-5 -3; 2 -9],[1 0; 0 2],[4 0; 0 4],[0 0; 0 0])
@test D_111 * D_221 == ss([-0.5 2 0; 0 0.2 -0.8; 0 -0.8 0.07],
        [0 0; 1 0; 0 2],[3 0 0],[0 0],0.005)

# Division
@test 1/C_222_d == ss([-6 -3; 2 -11],[1 0; 0 2],[-1 0; -0 -1],[1 -0; 0 1])
@test C_221/C_222_d == ss([-5 -3 -1 0; 2 -9 -0 -2; 0 0 -6 -3;
        0 0 2 -11],[1 0; 0 2; 1 0; 0 2],[1 0 0 0],[0 0])
@test 1/D_222_d == ss([-0.8 -0.8; -0.8 -1.93],[1 0; 0 2],[-1 0; -0 -1],
        [1 -0; 0 1],0.005)

# Indexing
@test size(C_222) == (2, 2)
@test size(C_212) == (2, 1)
@test C_222[1,1] == ss([-5 -3; 2 -9],[1; 0],[1 0],[0])
@test C_222[1:1,1] == ss([-5 -3; 2 -9],[1; 0],[1 0],[0])
@test C_222[1,1:2] == C_221
@test size(C_222[1,[]]) == (1,0)

# Naming signals
@test C_222_d_n.statenames == String["i","u"]
@test C_222_d_n.inputnames == String["e1","e2"]
@test C_222_d_n.outputnames == String["theta1","theta2"]

# Printing
res = ("StateSpace:\nA = \n          x1      x2 \n  x1   -5.0    -3.0  \n  x2"*
       "    2.0    -9.0  \nB = \n         u1     u2 \n  x1   1.0    0.0  \n"*
       "  x2   0.0    2.0  \nC = \n         x1     x2 \n  y1   1.0    0.0  \n"*
       "  y2   0.0    1.0  \nD = \n         u1     u2 \n  y1   0.0    0.0  \n"*
       "  y2   0.0    0.0  \n\nContinuous-time state-space model")
@test sprint(show, C_222) == res
res = ("StateSpace:\nA = \n          x1      x2 \n  x1    0.2    -0.8  \n  x2"*
       "   -0.8     0.07 \nB = \n         u1     u2 \n  x1   1.0    0.0  \n"*
       "  x2   0.0    2.0  \nC = \n         x1     x2 \n  y1   1.0    0.0  \n"*
       "  y2   0.0    1.0  \nD = \n         u1     u2 \n  y1   0.0    0.0  \n"*
       "  y2   0.0    0.0  \n\nSample Time: 0.005 (seconds)\n"*
       "Discrete-time state-space model")
@test sprint(show, D_222) == res
res = ("StateSpace:\nD = \n         u1     u2 \n  y1   4.0    0.0  \n  y2  "*
       " 0.0    4.0  \n\nStatic gain")
@test sprint(show, C_022) == res
res = ("StateSpace:\nD = \n         u1     u2 \n  y1   4.0    0.0  \n  y2  "*
       " 0.0    4.0  \n\nSample Time: 0.005 (seconds)\nStatic gain")
@test sprint(show, D_022) == res

# Errors
@test_throws ErrorException C_111 + C_222             # Dimension mismatch
@test_throws ErrorException C_111 - C_222             # Dimension mismatch
@test_throws ErrorException C_111 * C_222             # Dimension mismatch
@test_throws ErrorException D_111 + C_111             # Sampling time mismatch
@test_throws ErrorException D_111 - C_111             # Sampling time mismatch
@test_throws ErrorException D_111 * C_111             # Sampling time mismatch
D_diffTs = ss([1], [2], [3], [4], 0.1)
@test_throws ErrorException D_111 + D_diffTs            # Sampling time mismatch
@test_throws ErrorException D_111 - D_diffTs            # Sampling time mismatch
@test_throws ErrorException D_111 * D_diffTs            # Sampling time mismatch
@test_throws ErrorException 1/C_222                     # Not invertible
@test_throws ErrorException 1/C_212                     # Not invertible
@test_throws ErrorException ss([1 2], [1], [2], [3])      # Not square A
@test_throws ErrorException ss([1], [2 0], [1], [2])      # I/0 dim mismatch
@test_throws ErrorException ss([1], [2], [3 4], [1])      # I/0 dim mismatch
@test_throws ErrorException ss([1], [2], [3], [4], -0.1)  # Negative samping time
@test_throws ErrorException ss(eye(2), eye(2), eye(2), [0]) # Dimension mismatch
end


text:
# LSL.jl: Julia interface for Lab Streaming Layer
# Copyright (C) 2019 Samuel Powell

# LSL.jl: main module definition
module LSL

using liblsl_jll

using Libdl
using CEnum
using LightXML

export protocol_version, library_version, library_info, local_clock

# Include wrapper
module lib
  using CEnum
  import ..liblsl
  include(joinpath(@__DIR__, "..", "gen", "liblsl_common.jl"))
  include(joinpath(@__DIR__, "..", "gen", "liblsl_api.jl"))
end

using .lib

# Constants
const FOREVER = lib.LSL_FOREVER
const IRREGULAR_RATE = lib.LSL_IRREGULAR_RATE
const DEDUCED_TIMESTAMP = lib.LSL_DEDUCED_TIMESTAMP
const NO_PREFERENCE = lib.LSL_NO_PREFERENCE

# Type maps
_lsl_channel_format(::Type{Float32})    = lib.lsl_channel_format_t(1)
_lsl_channel_format(::Type{Float64})    = lib.lsl_channel_format_t(2)
_lsl_channel_format(::Type{String})     = lib.lsl_channel_format_t(3)
_lsl_channel_format(::Type{Int32})      = lib.lsl_channel_format_t(4)
_lsl_channel_format(::Type{Int16})      = lib.lsl_channel_format_t(5)
_lsl_channel_format(::Type{Cchar})      = lib.lsl_channel_format_t(6)
_lsl_channel_format(::Type{Int64})      = lib.lsl_channel_format_t(7)
_lsl_channel_format(::Type{T}) where T  = lib.lsl_channel_format_t(0)

function _jl_channel_format(lsl_format::lib.lsl_channel_format_t)

  if lsl_format == 1
    return Float32
  elseif lsl_format == 2
    return Float64
  elseif lsl_format == 3
    return String
  elseif lsl_format == 4
    return Int32
  elseif lsl_format == 5
    return Int16
  elseif lsl_format == 6
    return Cchar
  elseif lsl_format == 7
    return Int64
  else
    return Cvoid
  end

end

_jl_channel_format(::Type{Float32})    = lib.lsl_channel_format_t(1)
_jl_channel_format(::Type{Float64})    = lib.lsl_channel_format_t(2)
_jl_channel_format(::Type{String})     = lib.lsl_channel_format_t(3)
_jl_channel_format(::Type{Int32})      = lib.lsl_channel_format_t(4)
_jl_channel_format(::Type{Int16})      = lib.lsl_channel_format_t(5)
_jl_channel_format(::Type{Int64})      = lib.lsl_channel_format_t(7)
_jl_channel_format(::Type{T}) where T  = lib.lsl_channel_format_t(0)

# Error handling

function handle_error(errcode)
  if errcode >= 0
    return errcode
  elseif errcode == -1
    error("operation failed due to timeout")
  elseif errcode == -2
    error("the stream has been lost")
  elseif errcode == -3
    error("an argument was incorrectly specified")
  elseif errcode == -4
    error("an internal error has occurred")
  else
    error("an unknown error has occurred")
  end
  return 0
end

# Basic functions
"""
    protocol_version()

Return the version number of the protocol.

Clients with different minor versions are protocol-compatible with each
other while clients with different major versions will refuse to work
together.
"""
function protocol_version()
    major, minor = divrem(lib.lsl_protocol_version(), 100)
    return VersionNumber(major, minor)
end

"""
    library_version()

Return the version number of the underlying liblsl library.
"""
function library_version()
    major, minor = divrem(lib.lsl_library_version(), 100)
    return VersionNumber(major, minor)
end

"""
    library_info()

Return string containing library information.

Contents are for debugging purposes and should not be replied upon by a user application.
"""
library_info() = unsafe_string(lib.lsl_library_info())

"""
    local_clock()

Return local system time stamp in seconds.

The resolution is better than a milisecond. This reading can be used to assign time stamps
to samples as they are being acquired. If the "age" of a sample is known at a particular
time (e.g., from USB transmission delays), it can be used as an offset to local_clock() to
obtain a better estimate of when a sample was actually captured. See push_sample() for a use
case.
"""
local_clock() = lib.lsl_local_clock()

# High level object API
include("XMLElement.jl")
include("StreamInfo.jl")
include("StreamOutlet.jl")
include("StreamInlet.jl")
include("Resolver.jl")

end # module


text:
using KissABC
using AbstractMCMC
using Statistics
using Test
using Random
Random.seed!(1)

@testset "Factored" begin
    d = Factored(Uniform(0, 1), Uniform(100, 101))
    @test all((0, 100) .<= rand(d) .<= (1, 101))
    @test pdf(d, (0.0, 0.0)) == 0.0
    @test pdf(d, (0.5, 100.5)) == 1.0
    @test logpdf(d, (0.5, 100.5)) == 0.0
    @test logpdf(d, (0.0, 0.0)) == -Inf
    @test length(d) == 2
    m = Factored(Uniform(0.00, 1.0), DiscreteUniform(1, 2))
    sample = rand(m)
    @test 0 < sample[1] < 1
    @test sample[2] == 1 || sample[2] == 2
    @test pdf(m, sample) == 0.5
    @test logpdf(m, sample) ≈ log(0.5)
end

@testset "Push" begin
    push_p = KissABC.push_p
    a /′ b = (typeof(a) == typeof(b)) && all(a .== b)
    @test push_p(Normal(), 1) /′ 1.0
    @test push_p(DiscreteUniform(), 1.0) /′ 1
    @test push_p(Factored(Normal(), DiscreteUniform()), (2, 1.0)) /′ (2.0, 1)
    @test push_p(Product([Normal(), Normal()]), [2, 1]) /′ [2.0, 1.0]
end

@testset "Tiny Data, Approximate Bayesian Computation and the Socks of Karl Broman" begin
    function model((n_socks, prop_pairs), consts)
        n_picked = 11
        n_pairs = round(Int, prop_pairs * floor(n_socks / 2))
        n_odd = n_socks - 2 * n_pairs
        socks = sort([repeat(1:n_pairs, 2); (n_pairs+1):(n_pairs+n_odd)])
        picked_socks = socks[randperm(n_socks)][1:min(n_socks, n_picked)]
        lu = length(unique(picked_socks))
        sample_pairs = min(n_socks, n_picked) - lu
        sample_odds = lu - sample_pairs
        sample_pairs, sample_odds
    end

    prior_mu = 30
    prior_sd = 15
    prior_size = -prior_mu^2 / (prior_mu - prior_sd^2)

    pr_socks = NegativeBinomial(prior_size, prior_size / (prior_mu + prior_size))
    pr_prop = Beta(15, 2)

    pri = Factored(pr_socks, pr_prop)

    tinydata = (0, 11)
    nparticles = 5000
    modelabc = ApproxPosterior(pri, x -> sum(abs, model(x, 0) .- tinydata), 0.1)
    results_st = sample(modelabc, AIS(500), 5000, ntransitions = 100, progress = false)
    @test results_st[1] ≈ 46.2
    @test results_st[2] ≈ 0.866

    P =
        smc(
            pri,
            x -> sum(abs, model(x, 0) .- tinydata),
            nparticles = 5000,
            verbose = false,
            alpha = 0.99,
            r_epstol = 0,
            epstol = 0.01,
        ).P

    @test P[1] ≈ 46.2
    @test P[2] ≈ 0.866
end

@testset "Normal dist -> Dirac Delta inference" begin
    pri = Normal(1, 0.2)
    sim(μ) = μ * μ + 1
    cost(x) = abs(sim(x) - 1.5)
    abc = ApproxKernelizedPosterior(pri, cost, 0.001)
    res = sample(abc, AIS(12), 500, discard_initial = 1000, progress = false)

    @test sim(res) ≈ 1.5
    @test smc(pri, cost, epstol = 0.1).P ≈ 0.707
end

@testset "Normal dist -> Dirac Delta inference, MCMCThreads" begin
    pri = Normal(1, 0.2)
    sim(μ) = μ * μ + 1
    cost(x) = abs(sim(x) - 1.5)
    abc = ApproxKernelizedPosterior(pri, cost, 0.001)
    res = sample(
        abc,
        AIS(12),
        MCMCThreads(),
        100,
        50,
        discard_initial = 50 * 12,
        progress = false,
    )
    @show sim(res)
    @test sim(res) ≈ 1.5
end

@testset "Normal dist + Uniform Distr inference" begin
    pri = Factored(Normal(1, 0.5), DiscreteUniform(1, 10))
    sim((n, du)) = (n * n + du) * (n + randn() * 0.01)
    cost(x) = abs(sim(x) - 5.5)
    model_abc = ApproxPosterior(pri, cost, 0.01)
    res = sample(model_abc, AIS(100), 1000, discard_initial = 10000, progress = false)
    @test sim(Tuple(res)) ≈ 5.5
    @test smc(pri, cost).P[2] ≈ 5
end

function brownianrms((μ, σ), N, samples = 200)
    t = 0:N
    #    rand()<1/20 && sleep(0.001)
    @.(sqrt(μ * μ * t * t + σ * σ * t)) .* (0.95 + 0.1 * rand())
end

@testset "Inference on drifted Wiener Process" begin
    params = (0.5, 2.0)
    tdata = brownianrms(params, 30, 10000)
    prior = Factored(Uniform(0, 1), Uniform(0, 4))
    cost(x) = sum(abs, brownianrms(x, 30) .- tdata) / length(tdata)
    modelabc = ApproxPosterior(prior, cost, 0.1)
    sim = sample(modelabc, AIS(50), 100, discard_initial = 50000, progress = false)
    @test all(sim .≈ params)
    @test all(smc(prior, cost, verbose=false,min_r_ess=0.55).P .≈ params)
end

@testset "Classical Mixture Model 0.1N+N" begin
    st(res) =
        ((quantile(res, 0.1:0.1:0.9)-reverse(quantile(res, 0.1:0.1:0.9)))/2)[1+(end-1)÷2:end]
    st_n = [
        0.0,
        0.04680825481526908,
        0.1057221226763449,
        0.2682111969397526,
        0.8309228020477986,
    ]

    prior = Uniform(-10, 10)
    sim(μ) = μ + rand((randn() * 0.1, randn()))
    cost(x) = abs(sim(x) - 0.0)
    plan = ApproxPosterior(prior, cost, 0.01)
    res = sample(
        plan,
        AIS(50),
        2000,
        ntransitions = 100,
        discard_initial = 5000,
        progress = false,
    )
    plan = ApproxKernelizedPosterior(prior, cost, 0.01 / sqrt(2))
    resk = sample(
        plan,
        AIS(50),
        2000,
        ntransitions = 100,
        discard_initial = 5000,
        progress = false,
    )
    ressmc = smc(prior, cost, nparticles = 2000, alpha = 0.9, epstol = 0.01,mcmc_retrys=500,mcmc_tol=0.9,verbose=false).P
    testst(alg, r) = begin
        m = mean(abs, st(r) - st_n)
        println(":", alg, ": testing m = ", m)
        @show r
        m < 0.1
    end
    @test testst("Hard threshold", res)
    @test testst("Kernelized threshold", resk)
    @test testst("SMC", ressmc)
end

@testset "Usecase of issue #10" begin
    plan = ApproxPosterior(Normal(0, 1), x -> abs(x - 1.5), 0.01)
    res = sample(plan, AIS(20), 100, discard_initial = 2000, progress = false)
    @show res
    @test res ≈ 1.5
end

@testset "MVNormal vector test + 4 MCMCThreads" begin
    plan =
        ApproxPosterior(MultivariateNormal(4, 1.0), x -> abs(sum(abs2, x)^0.5 - 1.5), 0.01)
    res = sample(
        plan,
        AIS(20),
        MCMCThreads(),
        100,
        4,
        discard_initial = 10000,
        ntransitions = 40,
        progress = false,
    )
    @test mean(plan.cost(res)) < 0.01
end

@testset "CommonLogDensity: rosenbrock banana density" begin
    D = CommonLogDensity(
        2,
        rng -> randn(rng, 2),
        x -> -100 * (x[1] - x[2]^2)^2 - (x[2] - 1)^2,
    )
    @test length(D) == 2
    @test typeof(KissABC.unconditional_sample(Random.GLOBAL_RNG, D)) <: KissABC.Particle
    res = sample(
        D,
        AIS(50),
        1000,
        ntransitions = 100,
        discard_initial = 2000,
        progress = false,
    )
    @show res
    @test quantile(D.lπ(res), 0.97) > -0.69
end

@testset "Handling of ∞ costs" begin
    D = CommonLogDensity(
        2,
        rng -> rand(2) .* (2, 1) .- (1, 0),
        x -> ifelse(sum(abs2, x) <= 1, 0.0, -Inf),
    )
    D2 = CommonLogDensity(2, rng -> rand(2) .* (2, 1) .- (1, 0), x -> -Inf)
    res = sample(
        D,
        AIS(50),
        1000,
        ntransitions = 100,
        discard_initial = 5000,
        progress = false,
    )
    @test D.lπ(res) == 0
    @test_throws ErrorException sample(D2, AIS(50), 10, progress = false)
end

@testset "SMC" begin
    pp = Factored(Normal(0, 5), Normal(0, 5))
    cc((x, y)) = 50 * (x + randn() * 0.01 - y^2)^2 + (y - 1 + randn() * 0.01)^2

    R = smc(pp, cc, verbose = false, alpha = 0.9, nparticles = 500, epstol = 0.01, parallel=true).P
    @test R[1] ≈ 1
    @test R[2] ≈ 1

    cc2((x, y)) = rand((50 * (x + randn() * 0.01 - y^2)^2 + (y - 1 + randn() * 0.01)^2,Inf))

    R = smc(pp, cc2, verbose = false, alpha = 0.9, nparticles = 1000, epstol = 0.01, parallel=true).P

    @test R[1] ≈ 1
    @test R[2] ≈ 1
end

#benchmark
#=
using KissABC, Distributions, Random
function cost((u1, p1); n=10^6, raw=false)
    u2 = (1.0 - u1*p1)/(1.0 - p1)
    x = randexp(n) .* ifelse.(rand(n) .< p1, u1, u2)
    raw && return x
    sqrt(sum(abs2,[std(x)-2.2, median(x)-0.4]./[2.2,0.4]))
end

plan=ApproxPosterior(Factored(Uniform(0,1), Uniform(0.5,1)), cost, 0.01)

@show res=sample(plan, AIS(100),100,discard_initial=10000)

early stop:
 [0.49006664933267297, 0.49313860531909304, 0.49497013116625105]
 [0.8804136291097875, 0.8819843728641816, 0.8834306754737902]

=#
#plotting stuff

#=
using KissABC
using Distributions

function cost((μ,σ))
    x=randn(1000) .* σ .+ μ
    d1=mean(x)-2.0
    d2=std(x)-0.04
    hypot(d1,d2*50)
end

prior=Factored(Uniform(1,3),Truncated(Normal(0,0.05),0,100))
plan=ApproxKernelizedPosterior(prior,cost,0.005)
res=sample(plan,AIS(10),10000,ntransitions=50)

prsample=[rand(prior) for i in 1:10000]
μ_pr=getindex.(prsample,1)
σ_pr=getindex.(prsample,2)

μ_p=vec(res[:,1,:])
σ_p=vec(res[:,2,:])

mean(μ_p),std(μ_p)
mean(σ_p),std(σ_p)

cd(@__DIR__); pwd()

using PyPlot
pygui(true)
figure(figsize=1.5 .*(7.5,7.5).*(1,(sqrt(5)-1)/2),dpi=200)
subplot(2,2,1)
title("PRIOR")
hist(μ_pr,150,histtype="step",label=L" π(μ)",density=true,range=(0.95,3.05))

legend()
xlabel(L"\mu")
subplot(2,2,2)
title("POSTERIOR")
hist(μ_p,150,histtype="step",label=L" P(μ|{\rm data})",density=true,range=(1.9,2.1))

legend()
xlabel(L"\mu")
subplot(2,2,3)
hist(σ_pr,150,histtype="step",label=L" π(σ)",density=true,range=(-0.005,0.3))

xlabel(L"\sigma")
legend()
subplot(2,2,4)
hist(σ_p,150,histtype="step",label=L" P(σ|{\rm data})",density=true,range=(0.02,0.06))

xlabel(L"\sigma")
legend()
tight_layout()
PyPlot.savefig("../images/inf_normaldist.png")

=#
#=

using ApproxInferenceProblems,Distributions,KissABC
normalizep(X)=X./sum(X)
problem = ApproxInferenceProblem(BlowFly,T=1000,statistics = y -> normalizep(diff([count(x->x<=α,y) for α in 14:16:16014])) )
problem.model(rand(problem.prior))
function cost(X)
    s=problem.model(X)
    costs=0.0
    N=0.0
    for i in eachindex(s)
        if s[i]>0 && problem.data[i]>0
            costs+=abs(s[i]-problem.data[i])/(problem.data[i]+s[i])
            N+=1
        end
    end
    costs/N
end
cost(rand(problem.prior))
approx_density = ApproxPosterior(problem.prior,cost,0.01)

mcmc(approx_density, nparticles=100, generations=5000, parallel=true)

problem.target
=#


text:
function pred_accuracy(preds, obsv)
    return sum(preds .== obsv) / size(preds, 1)
end

function accuracy_measures(obsv::Array{T,1}, preds::Array{T,1}) where {T}
    cm = confusion_matrix(obsv, preds)
    n = sum(cm.matrix)
    if isone(length(cm.classes))
        @warn "Only one class in confusion matrix"
        TN, FN, FP, TP = 0.0, 0.0, 0.0, 1.0
    else
        TN, FN, FP, TP = cm.matrix
    end
    sensitivity = TP / (TP + FN)
    FN_rate = 1 - sensitivity
    specificity = TN / (TN + FP)
    FP_rate = 1 - specificity
    accuracy = (TP + TN) / n
    positive_predictive_power = TP / (TP + FP)
    odds_ratio = (TP * TN) / (FP * FN)
    kappa =
        ((TP + TN) - (((TP + FN) * (TP + FP) + (FP + TN) * (FN + TN)) / n)) /
        (n - (((TP + FN) * (TP + FP) + (FP + TN) * (FN + TN)) / n))
    TSS = 1 - (sensitivity + specificity)
    acc_mes = (
        cm=cm.matrix,
        TN=TN,
        FN=FN,
        FP=FP,
        TP=TP,
        sensitivity=sensitivity,
        FN_rate=FN_rate,
        specificity=specificity,
        FP_rate=FP_rate,
        accuracy=accuracy,
        positive_predictive_power=positive_predictive_power,
        odds_ratio=odds_ratio,
        kappa=kappa,
        TSS=TSS,
    )
    return acc_mes
end

function auc(
    model::Ensemble{S,T}, vld_features::Array{S}, vld_labels::Array{T,1}; plot=true, kw...
) where {S,T}
    presence_proba = apply_forest_proba(model, vld_features, [0, 1])[:, 2]
    thresholds = collect(0.0:0.01:1.0)
    thrsh_preds = [
        Int64.(replace(prob -> prob .>= thrsh ? 1 : 0, presence_proba)) for
        thrsh in thresholds
    ]
    acc_mes = [accuracy_measures(vld_labels, pred) for pred in thrsh_preds]
    FP_rates = map(x -> x.FP_rate, acc_mes)
    sensitivities = map(x -> x.sensitivity, acc_mes)
    score = auc_score(FP_rates, sensitivities)
    if plot
        p = auc_plot(FP_rates, sensitivities, score; kw...)
        return (score=score, FP_rates=FP_rates, sensitivities=sensitivities, plot=p)
    else
        return (score=score, FP_rates=FP_rates, sensitivities=sensitivities)
    end
end
function auc_score(FP_rates, sensitivities)
    FPrev, sensrev = FP_rates[end:-1:1], sensitivities[end:-1:1]
    deltas = FPrev[2:end] .- FPrev[1:(end - 1)]
    areas = deltas .* sensrev[2:end]
    score = sum(areas)
    return score
end
function auc_plot(FP_rates, sensitivities, score; kw...)
    p = plot(
        [0, 1],
        [0, 1];
        xlim=(0, 1),
        ylim=(0, 1),
        title="Receiver operating curve",
        xlabel="False-positive rate (1 - Specificity)",
        ylabel="True-positive rate (Sensitivity)",
        legend=:none,
        aspect_ratio=1,
    )
    plot!(
        p,
        FP_rates,
        sensitivities,
        ann=(0.8, 0.1, "AUC = $(round(score, digits = 3))"),
        kw...,
    )
    return p
end

# From https://github.com/marubontan/MLUtils.jl/blob/master/src/utils.jl
#=
function auc2(yTruth::Array{Int}, yScore::Array{Float64})
    dIndex = findall(1 .== yTruth)
    dnIndex = findall(0 .== yTruth)
    score = 0.0
    for ydn in dnIndex
        for yd in dIndex
            if yScore[yd] > yScore[ydn]
                score += 1
            elseif yScore[ydn] == yScore[yd]
                score += 0.5
            end
        end
    end
    return score / (length(dIndex) * length(dnIndex))
end
yTruth = vld_labels[:, nsp]
yScore = apply_forest_proba(model, vld_features, [0,1])[:,2]
auc2(yTruth, yScore)
=#

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
"""
    q4

Compute the solution of the Poisson equation of heat conduction with a nonzero
heat source. Quadrilateral four-node elements are used.
"""
module q4

using LinearAlgebra
using SparseArrays
using MeshCore.Exports
using MeshCore: @_check
using MeshSteward.Exports
using Elfel.Exports
using UnicodePlots

A = 1.0 # length of the side of the square
kappa =  1.0; # conductivity matrix
Q = -6.0; # internal heat generation rate
tempf(x, y) =(1.0 + x^2 + 2.0 * y^2);#the exact distribution of temperature
N = 1000;# number of subdivisions along the sides of the square domain

function genmesh()
    conn = Q4block(A, A, N, N)
    mesh = Mesh()
    attach!(mesh, conn)
    return mesh
end

function assembleKF(fesp, kappa, Q)
    function integrate!(am, av, geom, elit, qpit, kappa, Q)
        nedof = ndofsperel(elit)
        ke = LocalMatrixAssembler(nedof, nedof, 0.0)
        fe = LocalVectorAssembler(nedof, 0.0)
        for el in elit
            init!(ke, eldofs(el), eldofs(el))
            init!(fe, eldofs(el))
            for qp in qpit
                Jac, J = jacjac(el, qp)
                gradN = bfungrad(qp, Jac)
                JxW = J * weight(qp)
                N = bfun(qp)
                for j in 1:nedof
                    for i in 1:nedof
                        ke[i, j] += dot(gradN[i], gradN[j]) * (kappa * JxW)
                    end
                    fe[j] += N[j] * Q * JxW
                end
            end
            assemble!(am, ke)
            assemble!(av, fe)
        end
        return am, av
    end

    elit = FEIterator(fesp)
    qpit = QPIterator(fesp, (kind = :Gauss, order = 2))
    geom = geometry(fesp.mesh)
    am = start!(SysmatAssemblerSparse(0.0), ndofs(fesp), ndofs(fesp))
    av = start!(SysvecAssembler(0.0), ndofs(fesp))

    @time integrate!(am, av, geom, elit, qpit, kappa, Q)

    return finish!(am), finish!(av)
end

function solve!(T, K, F, nu)
    @time KT = K * T
    @time T[1:nu] = K[1:nu, 1:nu] \ (F[1:nu] - KT[1:nu])
end

function checkcorrectness(fesp)
    geom = geometry(fesp.mesh)
    ir = baseincrel(fesp.mesh)
    T = attribute(ir.right, "T")
    std = 0.0
    for i in 1:length(T)
        std += abs(T[i][1] - tempf(geom[i]...))
    end
    @_check (std / length(T)) <= 1.0e-9
end

function run()
    mesh = genmesh()
    fesp = FESpace(Float64, mesh, FEH1_Q4())
    bir = boundary(mesh);
    vl = connectedv(bir);
    locs = geometry(mesh)
    for i in vl
        setebc!(fesp, 0, i, 1, tempf(locs[i]...))
    end
    numberfreedofs!(fesp)
    numberdatadofs!(fesp)
    @show nunknowns(fesp)
    K, F = assembleKF(fesp, kappa, Q)
    T = fill(0.0, ndofs(fesp))
    gathersysvec!(T, fesp)
    solve!(T, K, F, nunknowns(fesp))
    scattersysvec!(fesp, T)
    makeattribute(fesp, "T", 1)
    checkcorrectness(fesp)
    vtkwrite("q4-T", baseincrel(mesh), [(name = "T",)])
end

end

@time q4.run()
# q4.run()


text:
function solve()
    n = parse(Int, readline())
    ret = 0
    for a=1:n-1
        ret += floor(Int, (n-1) / a)
    end
    ret
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
#
#   Allows the PETSc dynamic library to be used from Julia. (http::/julialang.org)
#     PETSc must be configured with --with-shared-libraries
#     You can build with or without MPI, but cannot run on more than one process
#
#  Load the PETSc dynamic library
#

module PETSc

using MPI
using ArrayViews
include("petsc_constants.jl")  # typedefs and constants
include("petsc_error.jl")  # error handling

import Base.show

# export all names
export PetscInitialize, PetscInitialized, getPETSC_COMM_SELF, PetscView, PetscIS, PetscDestroy, PetscISSetType, PetscISGetSizePetscISGetIndices, PetscDataTypeFromString, PetscDataTypeGetSize, PetscFinalize

# -------------------------------------
function echodemo(filename)
  f = open(filename)
  h = readall(f)
  close(f)
  pos = 0
  while (pos <= length(h))
    (ex,pos)=parse(h,pos)
    str = string(ex.args)
    if _jl_have_color
      print("\033[1m\033[30m")
    end
    #  the next line doesn't work right for multiple line commands like for loops
    println(str[2:strlen(str)-1])
    if _jl_have_color
      print(_jl_answer_color())
      println(" ") # force a color change, otherwise answer color is not used
    end
    e = eval(ex)
    println(" ")
  end
end

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

#=
PETSC_INSERT_VALUES = 1;
PETSC_ADD_VALUES    = 2;
PETSC_COPY_VALUES   = 0;

PETSC_NORM_1         = 0;
PETSC_NORM_2         = 1;
PETSC_NORM_FROBENIUS = 2;
PETSC_NORM_INFINITY  = 3;
PETSC_NORM_MAX       = PETSC_NORM_INFINITY;
=#

# -------------------------------------
#    These are the Julia interface methods for all the visible PETSc functions. Julia datatypes for PETSc objects simply contain the C pointer to the
#    underlying PETSc object.
#
# -------------------------------------
function PetscInitialized()
  init = Array(PetscBool, 1);
  err = ccall( (:PetscInitialized,  libpetsclocation),Int32,(Ptr{PetscBool},), init);

  return init[1]
end

function PetscFinalize()
  println("finalizing PETSc")
  gc() # call garbage collection to force all PETSc objects be destroy that are queued up for destruction
  return ccall( (:PetscFinalize,  libpetsclocation),Int32,());
end

function PetscInitialize()
  PetscInitialize([])
end

function PetscInitialize(args)
  PetscInitialize(args,"","")
end

function PetscInitialize(args,filename,help)
  # argument list starts with program name
  args = ["julia";args];
  println("typeof(args) = ", typeof(args))
  #
  #   If the user forgot to PetscFinalize() we do it for them, before restarting PETSc
  #

 init = PetscInitialized()
 if (init != 0)
#    gc() # call garbage collection to force all PETSc objects be destroy that are queued up for destruction
#    err = ccall( (:PetscFinalize,  libpetsclocation),Int32,());if (err != 0) return err; end
    err = PetscFinalize()

    if err != 0
      return err
    end
  end

  # convert arguments to Cstring
  arr = Array(ByteString, length(args))
  for i = 1:length(args)
    arr[i] = Base.cconvert(Cstring, args[i])
    # = cstring(args[i])
  end
#  ptrs = _jl_pre_exec(arr)
#  err = ccall(Libdl.dlsym(libpetsc, :PetscInitializeNoPointers),Int32,(Int32,Ptr{Ptr{UInt8}},Ptr{UInt8},Ptr{UInt8}), length(ptrs), ptrs,cstring(filename),cstring(help));

  err = ccall( (:PetscInitializeNoPointers,  libpetsclocation),Int32,(Int32,Ptr{Ptr{UInt8}},Cstring,Cstring), length(arr), arr,filename,help);
  return err
end

function getPETSC_COMM_SELF()
  comm = Array(MPI.CComm, 1)
  err = ccall( (:PetscGetPETSC_COMM_SELF,  libpetsclocation),Int32,(Ptr{comm_type},),comm);
  return convert(MPI.Comm, comm[1])
#   return MPI.COMM_WORLD.val
end

function PetscDataTypeFromString(name::AbstractString)
    ptype = Array(Cint, 1)
    found = Array(PetscBool, 1)
    ccall((:PetscDataTypeFromString,petsc),PetscErrorCode,(Cstring,Ptr{PetscDataType},Ptr{PetscBool}), name, ptype, found)

    return ptype[1], convert(Bool, found[1])
end

function PetscDataTypeGetSize(dtype::PetscDataType)
    datasize = Array(Csize_t, 1)
    ccall((:PetscDataTypeGetSize,petsc),PetscErrorCode,(PetscDataType,Ptr{Csize_t}), dtype, datasize)

    return datasize[1]
end

# -------------------------------------
#
#abstract PetscObject

function PetscView(obj)
  PetscView(obj,0)
end

type PetscIS
  pobj::Ptr{Void}
  function PetscIS(comm::MPI_Comm)
#    comm = PETSC_COMM_SELF();
    is = Array(Int64,1)
    err = ccall( (:ISCreate,  libpetsclocation),Int32,(comm_type,Ptr{Void}),comm,is)
#    if (err != 0)  # return type stability
#      return err
#    end
    is_ = new(is[1])
    finalizer(is_,PetscDestroy)
    # does not seem to be called immediately when is is no longer visible, is it called later during garbage collection?
    return is_
  end
end

  function PetscDestroy(is::PetscIS)
    if (is.pobj != 0) then
      err = ccall( ( :ISDestroy,  libpetsclocation),PetscErrorCode,(Ptr{Void},), &is.pobj);
    end
    is.pobj = 0
    println("ISDestroy called")
    return 0
  end

  function PetscIS(indices::Array{PetscInt})
    is = PetscIS()
    err = ccall( ( :ISSetType,  libpetsclocation),PetscErrorCode,(Ptr{Void},Cstring), is.pobj,"general");
    err = ccall( (:ISGeneralSetIndices,  libpetsclocation), PetscErrorCode, (Ptr{Void}, PetscInt, Ptr{PetscInt},Int32),is.pobj,length(indices), indices, PETSC_COPY_VALUES)
    return is
  end

  function PetscISSetType(vec::PetscIS, name)
    err = ccall( (:ISSetType,  libpetsclocation),PetscErrorCode,(Ptr{Void},Cstring), vec.pobj,name);
  end

  function PetscView(obj::PetscIS,viewer)
  # what is the viewer argument used for?
   err = ccall( ( :ISView,  libpetsclocation),PetscErrorCode,(Ptr{Void},Int64),obj.pobj,0);
  end

  function PetscISGetSize(obj::PetscIS)
    n = Array(PetscInt, 1)
    err = ccall( ( :ISGetSize,  libpetsclocation), PetscErrorCode,(Ptr{Void}, Ptr{PetscInt}), obj.pobj, n);
    return n[1]
  end

  # this should really take the indices array as an arguments
  function PetscISGetIndices(obj::PetscIS)
    len = PetscISGetSize(obj)
    indices = Array(PetscInt,len)
    err = ccall( (:ISGetIndicesCopy,  libpetsclocation), PetscErrorCode,(Ptr{Void},Ptr{PetscInt}),obj.pobj,indices);

    return indices
  end

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

include("petsc_vec.jl")

# -------------------------------------
include("petsc_mat.jl")

# -------------------------------------
include("petsc_ksp.jl")

# -------------------------------------
include("petsc_pc.jl")

#--------------------------------------
include("petsc_options.jl")

end  # end module


text:
using OMEinsum
using LinearAlgebra

function boundary_vector(A, O)
    @assert size(A, 2)^4 == size(O, 1)
    X = kron(A, A, A, A)
    return vec(X * O * X')
end

A = rand(3, 2)
O = rand(2^4, 2^4)
V = ein"ab,cd,ef,gh,ij,kl,mn,op,jlnpbdfh->acegikmo"(A,A,A,A,A,A,A,A,reshape(O, ntuple(_->2, 8)));
dot(V, V)

st = vec(transpose(A) * A)
psi = kron(st, st, st, st)
psi = vec(permutedims(reshape(psi, ntuple(_->2, 8)), (1,3,5,7,2,4,6,8)))
psi' * kron(O, O) * psi

X = kron(A, A, A, A)
V = X * O * X'
dot(V, V)
(X * O * transpose(X))[]
st
dot(A[:, 1], A[:, 1])
dot(A[:, 1], A[:, 2])
dot(A[:, 2], A[:, 1])
dot(A[:, 2], A[:, 2])


text:
import Pkg
Pkg.activate(".")
Pkg.instantiate()

using DelimitedFiles
using LinearAlgebra
using Printf
using SparseArrays

using BandedMatrices
using DataStructures: top
using Interpolations
using JSON
using OrdinaryDiffEq, DiffEqCallbacks
using Parameters
using RecursiveArrayTools
using Roots
using StaticArrays
using SummationByPartsOperators

import FFTW; FFTW.set_num_threads(1)

using PyCall, LaTeXStrings; import PyPlot; plt=PyPlot
inset_locator = pyimport("mpl_toolkits.axes_grid.inset_locator")

cycler = pyimport("cycler").cycler
line_cycler   = (cycler(color=["#E69F00", "#56B4E9", "#009E73", "#0072B2", "#D55E00", "#CC79A7", "#F0E442"]) +
                 cycler(linestyle=["-", "--", "-.", ":", "-", "--", "-."]))
marker_cycler = (cycler(color=["#E69F00", "#56B4E9", "#009E73", "#0072B2", "#D55E00", "#CC79A7", "#F0E442"]) +
                 cycler(linestyle=["none", "none", "none", "none", "none", "none", "none"]) +
                 cycler(marker=["4", "2", "3", "1", "+", "x", "."]))

plt.rc("axes", prop_cycle=line_cycler)
plt.rc("text", usetex=true)
plt.rc("text.latex", preamble="\\usepackage{newpxtext}\\usepackage{newpxmath}\\usepackage{commath}\\usepackage{mathtools}")
plt.rc("font", family="serif", size=18.)
plt.rc("savefig", dpi=100)
plt.rc("legend", loc="best", fontsize="medium", fancybox=true, framealpha=0.5)
plt.rc("lines", linewidth=2.5, markersize=10, markeredgewidth=2.5)

function convergence_orders(Ns, errors)
    @assert length(Ns) == length(errors)

    orders = similar(errors, length(errors)-1)
    orders[1] = NaN
    for i in 2:length(errors)
        orders[i-1] = - log(errors[i-1] / errors[i]) / log(Ns[i-1] / Ns[i])
    end
    orders
end

function linear_regression(x, y)
    A = hcat(one.(x), x)
    param = A \ y
    param
end

function even_values(val_N)
    res = Vector{eltype(val_N)}()
    for N in val_N
        N1 = iseven(N) ? N : N-1
        push!(res, N1)
    end
    res
end

function odd_values(val_N)
    res = Vector{eltype(val_N)}()
    for N in val_N
        N1 = iseven(N) ? N-1 : N
        push!(res, N1)
    end
    res
end

function evenodd_values(val_N)
    res = Vector{eltype(val_N)}()
    for N in val_N
        N1, N2 = iseven(N) ? (N, N+1) : (N-1, N)
        push!(res, N1)
        push!(res, N2)
    end
    res
end


text:
@testset "typeassert" begin
    test_rrule(typeassert, 1.1, Float64)
    test_frule(typeassert, 1.1, Float64)
end

@testset "ifelse" begin
    test_rrule(ifelse, true, 1.1, 2.0)
    test_frule(ifelse, false, 1.1, 2.0)

    test_rrule(ifelse, true, [1.1], [2.0]; check_inferred=false)
    test_frule(ifelse, false, [1.1], [2.0]; check_inferred=false)
end


text:
# Falta adaptar a função para diferentes resoluções radiométricas
# Para tanto, deve-se mudar o zeros(256)

function otsu_threshold_(in_img,t="image")
    tipo = typeof(in_img)

    img_orig = convert(Array{Float64},in_img)
    img = convert(Array{Int32},round(Int,in_img*255))

    d = size(img)
    hist = zeros(256)

    # Making the histogram
    for i=1:d[1]
        for j=1:d[2]
            hist[img[i,j]+1] += 1
        end
    end
    # end of the making

    total = d[1]*d[2]

    sumB = 0
    wB = 0
    max = 0.0
    sum1 = sum(img)

    for i=1:156
        wB += hist[i]
        wF = total - wB
        if wB == 0 || wF == 0
            continue
        end
        sumB = sumB + (i-1) * hist[i]
        mF = (sum1 - sumB) / wF
        between = wB * wF * ((sumB / wB) - mF) * ((sumB / wB) - mF)
        if between >= max
            level = i
            maximum = between
        end
    end

    level = level/255.0

    if t=="lvl"
        return(level)
    elseif t=="image"
        img_result = zeros(d)
        for i=1:d[1]
            for j=1:d[2]
                if img_orig[i,j] >= level
                    img_result[i,j] = 1.0
                end
            end
        end
        return(img_result)
    end
end


text:
const _TS_DATA_TO_METADATA_MAP = Dict(
    Deterministic => DeterministicMetadata,
    DeterministicSingleTimeSeries => DeterministicMetadata,
    AbstractDeterministic => DeterministicMetadata,
    Probabilistic => ProbabilisticMetadata,
    Scenarios => ScenariosMetadata,
    SingleTimeSeries => SingleTimeSeriesMetadata,
)

const _TS_METADATA_TO_DATA_MAP = Dict(
    # DeterministicMetadata is used for two types, and so this cannot be used for it.
    ProbabilisticMetadata => Probabilistic,
    ScenariosMetadata => Scenarios,
    SingleTimeSeriesMetadata => SingleTimeSeries,
)

function time_series_data_to_metadata(::Type{T}) where {T <: TimeSeriesData}
    return _TS_DATA_TO_METADATA_MAP[T]
end

function time_series_metadata_to_data(ts_metadata::TimeSeriesMetadata)
    return _TS_METADATA_TO_DATA_MAP[typeof(ts_metadata)]
end

function time_series_metadata_to_data(ts_metadata::DeterministicMetadata)
    return ts_metadata.time_series_type
end

is_time_series_sub_type(::Type{<:TimeSeriesMetadata}, ::Type{<:TimeSeriesData}) = false
is_time_series_sub_type(::Type{SingleTimeSeriesMetadata}, ::Type{StaticTimeSeries}) = true
is_time_series_sub_type(::Type{DeterministicMetadata}, ::Type{AbstractDeterministic}) = true
is_time_series_sub_type(::Type{DeterministicMetadata}, ::Type{Forecast}) = true
is_time_series_sub_type(::Type{ProbabilisticMetadata}, ::Type{Forecast}) = true
is_time_series_sub_type(::Type{ScenariosMetadata}, ::Type{Forecast}) = true


text:
"
    use an idmrg2 step to truncate/expand the bond dimension
"
@with_kw struct VumpsSvdCut <: Algorithm
    tol_gauge = Defaults.tolgauge
    tol_galerkin = Defaults.tol
    tol_eigenval = Defaults.tol
    trscheme = notrunc()
end

function changebonds_1(state::InfiniteMPS, H::Hamiltonian,alg::VumpsSvdCut,pars=params(state,H)) #would be more efficient if we also repeated pars
    #the unitcell==1 case is unique, because there you have a sef-consistency condition

    #expand the one site to two sites
    nstate = InfiniteMPS(repeat(state.AL,2))
    nH = repeat(H,2)

    nstate,npars = changebonds(nstate,nH,alg)

    A1 = nstate.AL[1]; A2 = nstate.AL[2]
    D1 = space(A1,1); D2 = space(A2,1)

    #collapse back to 1 site
    if dim(D2) != dim(D1)
        (nstate,npars) = changebonds(nstate,nH,SvdCut(trschemes = [truncdim(min(dim(D1),dim(D2)))]),npars)
    end

    collapsed = InfiniteMPS([nstate.AL[1]];tol=alg.tol_gauge)
    return collapsed, pars
end

function changebonds_n(state::InfiniteMPS, H::Hamiltonian,alg::VumpsSvdCut,pars=params(state,H))
    meps=0.0
    for loc in 1:length(state)
        @tensor AC2[-1 -2;-3 -4] := state.AC[loc][-1,-2,1]*state.AR[loc+1][1,-3,-4]

        (vals,vecs,_) = let state=state,pars=pars
            eigsolve(x->ac2_prime(x,loc,state,pars),AC2, 1, :SR, tol = alg.tol_eigenval; ishermitian=false )
        end
        nAC2=vecs[1]

        (vals,vecs,_)  = let state=state,pars=pars
            eigsolve(x->c_prime(x,loc+1,state,pars),state.CR[loc+1], 1, :SR, tol = alg.tol_eigenval; ishermitian=false )
        end
        nC2=vecs[1]

        #svd ac2, get new AL1 and S,V ---> AC
        (AL1,S,V,eps) = tsvd(nAC2, (1,2), (3,4), trunc=alg.trscheme)
        @tensor AC[-1,-2,-3]:=S[-1,1]*V[1,-2,-3]
        meps=max(eps,meps)

        #find AL2 from AC and C as in vumps paper
        QAC,_ = leftorth( AC,(1,2,),(3,), alg=TensorKit.Polar())
        QC ,_ = leftorth(nC2 ,(1,),(2,) , alg=TensorKit.Polar())
        @tensor AL2[-1,-2,-3] := QAC[-1,-2,1]*conj(QC[-3,1])

        #make a new state using the updated A's
        allAls = copy(state.AL)
        allAls[loc]   = permute(AL1,(1,2),(3,))
        allAls[loc+1] = permute(AL2,(1,2),(3,))

        state = InfiniteMPS(allAls; tol = alg.tol_gauge)
    end

    return state, pars
end

changebonds(state::InfiniteMPS,H,alg::VumpsSvdCut,pars=params(state,H)) = (length(state) == 1) ? changebonds_1(state,H,alg,pars) : changebonds_n(state,H,alg,pars);


text:
example_filepath(example_name, examples_dir="../examples/") =
    joinpath(examples_dir, example_name * ".jl")

@testset "Examples" begin
    @info "Testing examples..."

    @testset "One-dimensional diffusion" begin
        @info "  Testing one-dimensional diffusion example"

        replace_strings = [
            ("size = (1, 1, 128)", "size = (1, 1, 16)"),
            ("stop_iteration = 1000", "stop_iteration = 1"),
            ("simulation.stop_iteration += 10000", "simulation.stop_iteration += 100"),
            ("mp4(", "# mp4(")
        ]

        @test run_script(replace_strings, "one_dimensional_diffusion", example_filepath("one_dimensional_diffusion"))
    end

    @testset "Two-dimensional turbulence example" begin
        @info "  Testing two-dimensional turbulence example"

        replace_strings = [
            ("size=(128, 128, 1)", "size=(16, 16, 1)"),
            ("for i=1:100", "for i=1:1"),
            ("stop_iteration += 10", "stop_iteration += 1"),
            ("mp4(", "# mp4(")
        ]

        @test run_script(replace_strings, "two_dimensional_turbulence", example_filepath("two_dimensional_turbulence"))
    end

    for arch in archs
        @testset "Wind and convection-driven mixing example [$(typeof(arch))]" begin
            @info "  Testing wind and convection-driven mixing example [$(typeof(arch))]"

            replace_strings = [
                ("Nz = 32", "Nz = 16"),
                ("iteration_interval=10", "iteration_interval=1"),
                ("for i in 1:100", "for i in 1:1"),
                ("stop_iteration += 10", "stop_iteration += 1"),
                ("mp4(", "# mp4(")
            ]

            if arch isa GPU
                push!(replace_strings, ("architecture = CPU()", "architecture = GPU()"))
            end

            @test run_script(replace_strings, "ocean_wind_mixing_and_convection",
                             example_filepath("ocean_wind_mixing_and_convection"),
                             string(typeof(arch)))

            rm("ocean_wind_mixing_and_convection.jld2", force=true)
        end
    end

    @testset "Ocean convection with plankton example" begin
        @info "  Testing ocean convection with plankton example"

        replace_strings = [
            ("Nz = 128", "Nz = 16"),
            ("iteration_interval=100", "iteration_interval=1"),
            ("for i = 1:100", "for i = 1:1"),
            ("stop_iteration += 100", "stop_iteration += 1"),
            ("mp4(", "# mp4(")
        ]

        @test run_script(replace_strings, "ocean_convection_with_plankton",
                         example_filepath("ocean_convection_with_plankton"))
    end

    @testset "Internal wave example" begin
        @info "  Testing internal wave example"

        replace_strings = [
            ("Nx = 128", "Nx = 16"),
            ("iteration_interval = 20", "iteration_interval = 1"),
            ("for i=0:100", "for i=1:1"),
            ("stop_iteration += 20", "stop_iteration += 1"),
            ("mp4(", "# mp4(")
        ]

        @test run_script(replace_strings, "internal_wave", example_filepath("internal_wave"))
    end

    @testset "Eady turbulence" begin
        @info "  Testing Eady turbulence example"

        replace_strings = [
            ("Nh = 64", "Nh = 16"),
            ("Nz = 32", "Nz = 16"),
            ("end_time = 3day", "end_time = 1"),
            # Get rid of anything PyPlot/PyCall related
            ("using PyPlot, PyCall", ""),
            ("GridSpec =", "#"),
            ("fig =", "#"),
            ("gs =", "#"),
            ("fig.add_subplot", "#"),
            ("gcf()", "#"),
            ("makeplot!", "#makeplot!"),
            ("function #makeplot!(axs, model)", "function makeplot!(axs, model)")
        ]

        @test run_script(replace_strings, "eady_turbulence", example_filepath("eady_turbulence"))
    end
end