text:
"""
    abstract type Sector end

Abstract type representing unitary and pivotal (pre-)fusion categories.

The instances of a concrete subtype are simple objects of this concrete type.

Every concret subtype `I<:Sector` should implement the following methods:
*   `Base.one(::Type{I})`: unit object of `I`
*   `Base.conj(a::I)`: give ``a̅``, the conjugate label of `a`
*   `Base.isreal(::Type{I})`: whether the Fsymbol and Rsymbol of the sector is real
*   `Base.isless(a::I, b::I)`: give a partial order between simple objects
*   `FusionStyle(::Type{I})`: `UniqueFusion()`, `SimpleFusion()` or `GenericFusion()`
*   `⊗(a::I, b::I)`: iterable with unique fusion outputs of ``a ⊗ b``
    (i.e. don't repeat in case of multiplicities)
*   `Nsymbol(a::I, b::I, c::I)`: number of times `c` appears in `a ⊗ b`; return a Bool in
    the MultiplicityFreeFusion case.
*   `Fsymbol(a::I, b::I, c::I, d::I, e::I, f::I)`: scalar (`MultiplicityFreeFusion`)
    or rank-4 Array (`GenericFusion`)
*   (optionally)`vertex_ind2label(i::Int, a::I, b::I, c::I)`: a custom label for the `i`th
    copy of `c` appearing in `a ⊗ b`
*   (optioanlly)`vertex_labeltype(::Type{I})`: the type of labels for the fusion vertices
*   (optionally)`dim(a::I)`: quantum dimension of simple object `a`; sqrtdim(), isqrtdim()
*   (optionally)`frobeniusschur(a::I)`: Frobenius-Schur indicator of `a`
*   (optionally)`Bsymbol(a::I, b::I, c::I)`: B-symbol: scalar (`MultiplicityFreeFusion`)
    or matrix (`GenericFusion`)
*   `BraidingStyle(::Type{I})`: `Bosonic()`, `Fermionic()` or `Anyonic()`
*   `Rsymbol(a::I, b::I, c::I)`: R-symbol: scalar (`MultiplicityFreeFusion`) or matrix
    (`GenericFusion`)
*   (optionally)`twist(a::I)`: twist of simple object `a`

Furthermore, an iterator over all simple objects of the sector `I` is implemented by the
singleton type [`SectorValues{I}`](@ref).
"""
abstract type Sector end

"""
    struct SectorValues{I<:Sector}

Singleton type to represent an iterator over all simple objects of sector `I`, whose
instance is obtained by `values(I)`.

For a new `I::Sector`, the following should be defined:
*   `Base.iterate(::SectorValues{I}[, state])`: iterate over the simple objects
*   `Base.IteratorSize(::Type{SectorValues{I}})`: `HasLenght()` or `IsInfinite()`

If `IteratorSize(I) == HasLength()`, also the following must be implemented:
*   `Base.length(::SectorValues{I})`: the number of different simple objects
*   `Base.getindex(::SectorValues{I}, i::Int)`: map from an index `i` to a simple object
*   `findindex(::SectorValues{I}, a::I)`: map from a simple object `a::I` to an index
    `i::Integer ∈ 1:length(values(I))`
"""
struct SectorValues{I<:Sector} end
Base.IteratorEltype(::Type{<:SectorValues}) = HasEltype()
Base.eltype(::Type{SectorValues{I}}) where {I<:Sector} = I

"""
    Base.values(::Type{I}) where {I<:Sector}

Return the iterator that generate all simple objects of sector `I`.
"""
Base.values(::Type{I}) where {I<:Sector} = SectorValues{I}()

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

Return the unit object of the sector.
"""
Base.one(a::Sector) = one(typeof(a))

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

Return the conjugate label `conj(a)`.

In the language of category, `conj(a)` is `a̅`, which is in the set of simple objects and
is isomorphic to the dual of `a`, i.e. `a*`. Since we only represent the simple objects in
the codes, we don't have an instance that refer to `a*`, and we identify `dual(a)` and
`conj(a)`.

See also: [`conj`](@ref)
"""
dual(a::Sector) = conj(a)

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

Extend `Base.conj`. Return `a̅`, which is the conjugate or dual label of `a` which is in the
set of simple objects.
"""
Base.conj(a::Sector) = conj(a)

"""
    isless(a::I,b::I) where {I<:Sector} -> Bool

Extend `Base.isless()`. Give a partial order between simple objects in Sector `I`.
"""
Base.isless(a::I, b::I) where {I<:Sector} = isless(a,b)

# Fusion
abstract type FusionStyle end
struct UniqueFusion <: FusionStyle end # unique fusion output when fusion two sectors

abstract type MultipleFusion <: FusionStyle end
struct SimpleFusion <: MultipleFusion end # multiple fusion but multiplicity free
struct GenericFusion <: MultipleFusion end # multiple fusion with multiplicities

const MultiplicityFreeFusion = Union{UniqueFusion, SimpleFusion}

Base.:&(f::F, ::F) where {F<:FusionStyle} = f
Base.:&(f1::FusionStyle, f2::FusionStyle) = f2 & f1

Base.:&(::SimpleFusion, ::UniqueFusion) = SimpleFusion()
Base.:&(::GenericFusion, ::UniqueFusion) = GenericFusion()
Base.:&(::GenericFusion, ::SimpleFusion) = GenericFusion()

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

Return the type of fusion behavior of sector `I`, which can be either
*   `UniqueFusion()`: single fusion output when fusing any two simple objects;
*   `SimpleFusion()`: multiple outputs, but every output occurs at most one,
    also known as multiplicity-free (e.g. irreps of ``SU(2)``);
*   `GenericFusion()`: multiple outputs that can occur more than once (e.g. irreps
    of ``SU(3)``).
There is an abstract supertype `MultipleFusion` of which both `SimpleFusion` and
`GenericFusion` are subtypes.

Each secteor should implement this function to indicate which fusion style it has.
"""
FusionStyle(a::Sector) = FusionStyle(typeof(a))

"""
    fusiontensor(a::I, b::I, c::I) where {I<:AbstractIrrep{G<:Group}}

Return the fusiontensor ``X^{ab}_{c,μ}: c → a ⊗ b`` as a rank-4 tensor with size
`(dim(a),dim(b),dim(c),Int(Nsymbol(a,b,c)))`.
"""
function fusiontensor end

"""
    ⊗(a::I, b::I) where {I<:Sector}
    ⊗(a::I, b::I, c::I, rest::Vararg{I}) where {I<:Sector}

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

Note that every element `c` should appear at most once, fusion degeneracies (if
`FusionStyle(I) == GenericFusion()`) should be accessed via `Nsymbol(a, b, c)`.

The arguments can be more than two.
"""
@inline function ⊗(a::I, b::I, c::I, rest::Vararg{I}) where {I<:Sector}
    if FusionStyle(I) isa UniqueFusion
        return a ⊗ first(⊗(b, c, rest...))
    else
        s = Set{I}()
        for d in ⊗(b, c, rest...)
            for e in a ⊗ d
                push!(s, e)
            end
        end
        return s
    end
end

"""
    Nsymbol(a::I, b::I, c::I) where {I<:Sector} -> Bool or Int

Return a `Bool` if `FusionStyle(I) == UniqueFusion()` or `SimpleFusion()`.

Return an `Int` representing the number of times `c` appears in the fusion product `a ⊗ b`
if `FusionStyle(I) == GenericFusion()`.
"""
function Nsymbol end

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

The F-symbol ``F^{abc}_d`` that associates the two different fusion orders of simple objects
`a`, `b` and `c` into an ouput simple object `d`, using either an intermediate object
``a ⊗ b → e`` or ``b ⊗ c → f``:
```
a-<-μ-<-e-<-ν-<-d                                     a-<-λ-<-d
    ∨       ∨       -> Fsymbol(a,b,c,d,e,f)[μ,ν,κ,λ]      ∨
    b       c                                             f
                                                          v
                                                      b-<-κ
                                                          ∨
                                                          c
```
Return a number, if `FusionStyle(I) == UniqueFusion()` or `SimpleFusion()`.

Return a rank-4 array of size
`(Nsymbol(a, b, e), Nsymbol(e, c, d), Nsymbol(b, c, f), Nsymbol(a, f, d))` if
`FusionStyle(I) == GenericFusion()`.
"""
function Fsymbol end

"""
    vertex_ind2label(k::Int, a::I, b::I, c::I) where {I<:Sector}

Convert the index `k` of the fusion vertex ``a ⊗ b → c`` into a label.

For `FusionStyle(I) == UniqueFusion()` or `SimpleFusion()`, the default is to suppress
vertex labels by setting them equal to `nothing`.

For `FusionStyle(I) == GenericFusion()`, the default is to just use `k`, unless a
specialized method is provided.
"""
vertex_ind2label(k::Int, a::I, b::I, c::I) where {I<:Sector}=
    _ind2label(FusionStyle(I), k::Int, a::I, b::I, c::I)
_ind2label(::UniqueFusion, k, a, b, c) = nothing
_ind2label(::SimpleFusion, k, a, b, c) = nothing
_ind2label(::GenericFusion, k, a, b, c) = k

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

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

# properties that can be determined by Fsymbol and Nsymbol
"""
    dim(a::Sector)

Return the (quantum) dimension of the sector `a`.
"""
function dim(a::Sector)
    if FusionStyle(a) isa UniqueFusion
        1
    elseif FusionStyle(a) isa SimpleFusion
        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

"""
    sqrtdim(a::Sector)

Return the square root of the quantum dimension of the sector `a`.
"""
sqrtdim(a::Sector) = (FusionStyle(a) isa UniqueFusion) ? 1 : sqrt(dim(a))

"""
    isqrtdim(a::Sector)

Return the inverse of the square root of the quantum dimension of the sector `a`.
"""
isqrtdim(a::Sector) = (FusionStyle(a) isa UniqueFusion) ? 1 : inv(sqrt(dim(a)))

"""
    frobeniusschur(a::Sector)

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

"""
    Bsymbol(a::I, b::I, c::I) where {I<: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(I)` is `UniqueFusion()` or `SimpleFusion()`, 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::I, b::I, c::I) where {I<:Sector}
    if FusionStyle(I) isa UniqueFusion || FusionStyle(I) isa SimpleFusion
        (sqrtdim(a)*sqrtdim(b)*isqrtdim(c))*Fsymbol(a, b, dual(b), a, c, one(a))
    else
        reshape((sqrtdim(a)*sqrtdim(b)*isqrtdim(c))*Fsymbol(a, b, dual(b), a, c, one(a)),
            (Nsymbol(a, b, c), Nsymbol(c, dual(b), a)))
    end
end

# Braiding
abstract type BraidingStyle end # generic braiding
abstract type HasBraiding <: BraidingStyle end
struct NoBraiding <: BraidingStyle end

abstract type SymmetricBraiding <: HasBraiding end
struct Bosonic <: SymmetricBraiding end # all twists are one
struct Fermionic <: SymmetricBraiding end # twists one and minus one

struct Anyonic <: HasBraiding 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()
Base.:&(::Bosonic, ::NoBraiding) = NoBraiding()
Base.:&(::Fermionic, ::NoBraiding) = NoBraiding()
Base.:&(::Anyonic, ::NoBraiding) = NoBraiding()

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

Return the type of braiding and twist behavior of simple objects of type `I`, which can be
*   `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 `SimpleFusion` or
    `GenericFusion` 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:
`isone(Rsymbol(b,a,c)*Rsymbol(a,b,c)) == true`) and permutations are uniquely defined.

Each secteor should implement this function to indicate which braiding style it has.
"""
BraidingStyle(a::Sector) = BraidingStyle(typeof(a))

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

The R-symbol ``R^{ab}_c`` that maps from ``c → a ⊗ b`` to ``c → b ⊗ a`` as in
```
a -<-μ-<- c                                 b -<-ν-<- c
     ∨          -> Rsymbol(a,b,c)[μ,ν]           v
     b                                           a
```
Return a number if `FusionStyle(I) == UniqueFusion()` or `SimpleFusion()`.

Return a square matrix with row and column size `Nsymbol(a,b,c) == Nsymbol(b,a,c)` if
`FusionStyle(I) == GenericFusion()`.
"""
function Rsymbol end

# properties that can be determined by Rsymbol and Fsymbol
"""
    twist(a::Sector)

Return the twist of a simple object `a`.
"""
twist(a::Sector) = sum(dim(b)/dim(a)*tr(Rsymbol(a,a,b)) for b in a ⊗ a)

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

Return whether the topological data (Fsymbol, Rsymbol) of the sector is real or not (in
which case it is complex).
"""
function Base.isreal(I::Type{<:Sector})
    u = one(I)
    if BraidingStyle(I) isa HasBraiding
        return (eltype(Fsymbol(u, u, u, u, u, u))<:Real) && (eltype(Rsymbol(u, u, u))<:Real)
    else
        return (eltype(Fsymbol(u, u, u, u, u, u))<:Real)
    end
end

"""
    struct SectorSet{I<:Sector, F, S}
        f::F
        set::S
    end

An iterator that applies x->convert(I, f(x)) on the elements of set; if f is not provided it
is just taken as the function identity.
"""
struct SectorSet{I<:Sector, F, S}
    f::F
    set::S
end
SectorSet{I}(::Type{F}, set::S) where {I<:Sector, F, S} = SectorSet{I, Type{F}, S}(F, set)
SectorSet{I}(f::F, set::S) where {I<:Sector, F, S} = SectorSet{I, F, S}(f, set)
SectorSet{I}(set) where {I<:Sector} = SectorSet{I}(identity, set)

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

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

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

# possible sectors
include("trivial.jl")
include("groups.jl")
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


text:
# # 220: 2D Nonlinear Poisson with boundary reaction and boundary species
# ([source code](SOURCE_URL))

module Example220_NonlinearPoisson2D_BoundarySpecies

using Printf
using VoronoiFVM

function main(;n=10,Plotter=nothing,verbose=false,unknown_storage=:sparse)

    h=1.0/convert(Float64,n)
    X=collect(0.0:h:1.0)
    Y=collect(0.0:h:1.0)

    grid=VoronoiFVM.Grid(X,Y)

    k=1.0
    eps=1.0
    physics=VoronoiFVM.Physics(
    num_species=3,
    breaction=function(f,u,node)
        if  node.region==2
            f[1]=k*(u[1]-u[3])
            f[3]=k*(u[3]-u[1])+ k*(u[3]-u[2])
            f[2]=k*(u[2]-u[3])
        end
    end,

    bstorage=function(f,u,node)
        if  node.region==2
            f[3]=u[3]
        end
    end,

    flux=function(f,u,edge)
        f[1]=eps*(u[1,1]-u[1,2])
        f[2]=eps*(u[2,1]-u[2,2])
    end,

    source=function(f,node)
        x1=node[1]-0.5
        x2=node[2]-0.5
        f[1]=exp(-20.0*(x1^2+x2^2))
    end,

    storage=function(f,u,node)
        f[1]=u[1]
        f[2]=u[2]
    end
    )

    sys=VoronoiFVM.System(grid,physics,unknown_storage=unknown_storage)

    enable_species!(sys,1,[1])
    enable_species!(sys,2,[1])
    enable_boundary_species!(sys,3,[2])

    function tran32!(a,b)
        a[1]=b[2]
    end

    bgrid2=subgrid(grid,[2],boundary=true,transform=tran32!)

    inival=unknowns(sys)
    inival.=0.0
    U=unknowns(sys)

    eps=1.0e-2

    control=VoronoiFVM.NewtonControl()
    control.verbose=verbose
    control.tol_linear=1.0e-5
    control.tol_relative=1.0e-5
    control.max_lureuse=0
    tstep=0.01
    time=0.0
    istep=0
    u5=0
    while time<1
        time=time+tstep
        solve!(U,inival,sys,control=control,tstep=tstep)
        inival.=U
        if verbose
            @printf("time=%g\n",time)
        end
        tstep*=1.0
        istep=istep+1
        U_bound=view(U[3,:],bgrid2)
        u5=U_bound[5]
        if isplots(Plotter)
            p1=Plotter.contourf(X,Y,reshape(U[1,:],length(X),length(Y)),levels=collect(0:0.1:0.6),clim=(0,0.6),colorbar=:right,color=:viridis,title=@sprintf("max1=%g max2=%g maxb=%g\n",maximum(U[1,:]),maximum(U[2,:]),maximum(U_bound[3,:])))
            p2=Plotter.contourf(X,Y,reshape(U[2,:],length(X),length(Y)),levels=collect(0:0.0001:0.002),clim=(0,0.002), colorbar=:right,color=:viridis)
            p3=Plotter.plot(grid=true,ylims=(0,0.0025))
            plot(bgrid2,U_bound,p=p3,show=false,Plotter=Plotter)
            p=Plotter.plot(p1,p2,p3,layout=(3,1) )
            Plotter.gui(p)
        end
    end
    return u5
end

function test()
    main(unknown_storage=:sparse) ≈ 0.0020781361856598
    main(unknown_storage=:dense) ≈ 0.0020781361856598
end
end

ext: jl


text:
using YaoBlocks, YaoArrayRegister, BitBasis
using YaoBlocks: eigenbasis
using Random, Test
using YaoBase: QubitMismatchError

function check_eigenbasis(op)
    E, V = eigenbasis(op)
    mat(V) * mat(E) * mat(V') ≈ mat(op)
end

@testset "eigen basis" begin
    for op in [
        X,
        Y,
        H,
        Z,
        put(4, 2 => X),
        kron(X, Y, Z),
        repeat(3, H, 2:3),
        cache(Rx(0.5)),
        Daggered(Rx(0.5)),
        rot(SWAP, 0.5),
        time_evolve(kron(X, X), 0.5),
        3 * put(3, 2 => X),
        chain(10, put(2 => X), chain(put(4 => Y))),
        chain(10, put(2 => X), chain(put(2 => Y))),
    ]
        @show op
        @test check_eigenbasis(op)
    end
end

@testset "better operator measure" begin
    Random.seed!(21)
    nbit = 3
    for op in [
        put(nbit, 2 => X),
        control(nbit, 2, 1 => X),
        kron(X, X, X),
        put(nbit, (3, 1, 2) => kron(X, Y, Z)),
        repeat(nbit, X, 2:3),
        +(put(nbit, 2 => X), put(nbit, 1 => Rx(π))),
        2.8 * put(nbit, 1 => X),
        chain(nbit, put(nbit, 3 => X), put(nbit, 1 => Z)),
        cache(put(nbit, 2 => X)),
        Daggered(put(nbit, 2 => Rx(0.8))),
    ]
        reg = rand_state(nbit)
        reg2 = copy(reg)
        @show op
        @test isapprox(
            sum(measure(op, reg2; nshots = 100000)) / 100000,
            expect(op, reg),
            rtol = 0.1,
        )
        @test reg ≈ reg2

        reg = rand_state(nbit; nbatch = 10)
        reg2 = copy(reg)
        @test isapprox(
            dropdims(sum(measure(op, reg2; nshots = 100000), dims = 1), dims = 1) / 100000,
            expect(op, reg),
            rtol = 0.1,
        )
        @test reg ≈ reg2
    end
end

@testset "measure op, specify locs" begin
    Random.seed!(21)
    Nbit = 4
    nbit = 3
    locs = (1, 3, 2)
    @test_throws QubitMismatchError measure(X, rand_state(5), locs)
    @test_throws QubitMismatchError measure(X, rand_state(5), AllLocs())
    for op in [
        put(nbit, 2 => X),
        control(nbit, 2, 1 => X),
        kron(X, Y, Z),
        put(nbit, (3, 1, 2) => kron(X, Y, Z)),
        repeat(nbit, X, 2:3),
        +(put(nbit, 2 => X), put(nbit, 1 => Rx(π))),
        2.8 * put(nbit, 1 => X),
        chain(nbit, put(nbit, 3 => X), put(nbit, 1 => Z)),
        cache(put(nbit, 2 => X)),
        Daggered(put(nbit, 2 => Rx(0.8))),
    ]
        reg = rand_state(Nbit)
        reg2 = copy(reg)
        @show op
        @test isapprox(
            sum(measure(op, reg2, locs; nshots = 100000)) / 100000,
            expect(put(Nbit, locs => op), reg),
            rtol = 0.2,
        )
        @test reg ≈ reg2
    end
end

@testset "better operator measure!" begin
    Random.seed!(9)
    nbit = 3
    for op in [
        put(nbit, 2 => X),
        control(nbit, 2, 1 => X),
        kron(X, X, X),
        put(nbit, (3, 1, 2) => kron(X, X, X)),
        repeat(nbit, X, 2:3),
        +(put(nbit, 2 => X), put(nbit, 1 => im * Rx(π))),
        2.8 * put(nbit, 1 => X),
        chain(nbit, put(nbit, 3 => X), put(nbit, 1 => Z)),
        chain(nbit, put(nbit, 3 => X), put(nbit, 3 => Z)),
        cache(put(nbit, 2 => X)),
        Daggered(put(nbit, 2 => im * Rx(π))),
    ]
        @show op
        reg = rand_state(nbit)
        res = measure!(op, reg)
        @test isapprox(res, sum(measure(op, reg; nshots = 100))[1] / 100, rtol = 0.1)

        reg = rand_state(nbit; nbatch = 10)
        res = measure!(op, reg)
        @test isapprox(
            res,
            vec(sum(measure(op, reg; nshots = 100); dims = 1)) / 100,
            rtol = 0.1,
        )
    end
end

@testset "commute to eachother" begin
    @test YaoBlocks.simple_commute_eachother([put(5, 2 => X), put(5, 3 => Y)])
    @test !YaoBlocks.simple_commute_eachother([put(5, 2 => X), put(5, 2 => Y)])
end

@testset "fix measuring add" begin
    circ_wfn = zero_state(2) |> put(2, 1=>H)
    sp = put(2, 1=>(X + Val(1im)*Y)/2)
    X_diff = sp + sp'
    @test measure(X_diff, circ_wfn, nshots=100) |> sum ≈ 100.0
end

@testset "measure operator with post processing" begin
    # measure! and reset
    reg = rand_state(8; nbatch=32)
    op = repeat(5, X, 1:5)
    @test_throws ArgumentError measure!(ResetTo(0), op, reg, 2:6)
    @test_throws ArgumentError measure!(RemoveMeasured(), op, reg, 2:6)
end


text:
module FlexibleFunctors

import ConstructionBase: constructorof

export destructure, ffunctor, fieldmap

include("param_functors.jl")

end # module


text:
# Example taken from Golightly and Wilkinson 2005 `Bayesian Inference for Stochastic
# Kinetic Models Using a Diffusion Approximation`

using Bridge
using StaticArrays, LinearAlgebra
import Bridge: b, σ, B, β, a, constdiff

struct Prokaryote{T} <: ContinuousTimeProcess{ℝ{4,T}}
    c₁::T
    c₂::T
    c₃::T
    c₄::T
    c₅::T
    c₆::T
    c₇::T
    c₈::T
    K::Float64

    function Prokaryote(c₁::T, c₂::T, c₃::T, c₄::T, c₅::T, c₆::T, c₇::T,
                        c₈::T, K) where T
        new{T}(c₁, c₂, c₃, c₄, c₅, c₆, c₇, c₈, K)
    end
end

function b(t, x, P::Prokaryote{T}) where T # x <-> (RNA, P, P₂, DNA)
    ℝ{4}(P.c₃*x[4] - P.c₇*x[1],
         P.c₄*x[1] + 2.0*P.c₆*x[3] - P.c₅*x[2]*(x[2]-1)-P.c₈*x[2],
         P.c₂*(P.K - x[4]) + 0.5*P.c₅*x[2]*(x[2]-1) - P.c₁*x[3]*x[4] - P.c₆*x[3],
         P.c₂*(P.K - x[4]) - P.c₁*x[3]*x[4])
end

function σ(t, x, P::Prokaryote{T}) where T
    σσᵀ = a(t, x, P)
    cholesky(σσᵀ).U'
end

#=
 c₃DNA+c₇RNA |         0          |           0          |          0
--------------------------------------------------------------------------------
     0       |  c₄RNA+2c₅P(P-1)   |  -c₅P(P-1)-2c₆P₂     |          0
             |      +4c₆P₂+c₈P    |                      |
--------------------------------------------------------------------------------
     0       | -c₅P(P-1)-2c₆P₂    |  c₁P₂DNA+c₂(K-DNA)   |   c₁P₂DNA
             |                    |    +0.5c₅P(P-1)+c₆P₂ |    +c₂(K-DNA)
--------------------------------------------------------------------------------
     0       |         0          |   c₁P₂DNA+c₂(K-DNA)  |    c₁P₂DNA
                                                                +c₂(K-DNA)
=#
# x <-> (RNA, P, P₂, DNA)
function a(t, x, P::Prokaryote{T}) where T
    _a₃₂ = P.c₅*x[2]*(x[2]-1.0)+2.0*P.c₆*x[3]
    a₄₄ = P.c₁*x[3]*x[4]+P.c₂*(P.K-x[4])

    @SMatrix [P.c₃*x[4]+P.c₇*x[1]  0.0  0.0  0.0;
              0.0  P.c₄*x[1]+2.0*_a₃₂+P.c₈*x[2]  -_a₃₂  0.0;
              0.0  -_a₃₂  a₄₄+0.5*_a₃₂  a₄₄;
              0.0  0.0  a₄₄  a₄₄]
end

domain(P::Prokaryote) = BoundedDomain( LowerBoundedDomain((0.0, 1.0, 0.0, 0.0), (1,2,3,4)),
                                       UpperBoundedDomain((P.K,), (4,)) )
constdiff(::Prokaryote) = false
clone(P::Prokaryote, θ) = Prokaryote(θ..., P.K)
params(P::Prokaryote) = [P.c₁, P.c₂, P.c₃, P.c₄, P.c₅, P.c₆, P.c₇, P.c₈]

struct ProkaryoteAux{O,R,S1,S2} <: ContinuousTimeProcess{ℝ{4,R}}
    c₁::R
    c₂::R
    c₃::R
    c₄::R
    c₅::R
    c₆::R
    c₇::R
    c₈::R
    K::Float64
    t::Float64
    u::S1
    T::Float64
    v::S2

    function ProkaryoteAux(c₁::R, c₂::R, c₃::R, c₄::R, c₅::R, c₆::R, c₇::R,
                           c₈::R, K, t, u::S1, T, v::S2, ::O) where {O,R,S1,S2}
        new{O,R,S1,S2}(c₁, c₂, c₃, c₄, c₅, c₆, c₇, c₈, K, t, u, T, v)
    end
end

observables(::ProkaryoteAux{O}) where O = O()

function B(t, P::ProkaryoteAux{Val{(true,true,true,true)}})
    @SMatrix [-P.c₇  0.0  0.0  P.c₃;
              P.c₄  P.c₅-P.c₈-2*P.c₅*P.v[2]  2.0*P.c₆  0.0;
              0.0  P.c₅*P.v[2]-0.5*P.c₅  -P.c₁*P.v[4]-P.c₆  -P.c₂-P.c₁*P.v[3];
              0.0  0.0  -P.c₁*P.v[4]  -P.c₂-P.c₁*P.v[3]]
end

function β(t, P::ProkaryoteAux{Val{(true,true,true,true)}})
    @SVector [0.0,
              P.c₅*P.v[2]^2,
              P.c₂*P.K + P.c₁*P.v[3]*P.v[4] - 0.5*P.c₅*P.v[2]^2,
              P.c₂*P.K + P.c₁*P.v[3]*P.v[4]]
end

function σ(t, P::ProkaryoteAux{Val{(true,true,true,true)}})
    σσᵀ = a(t, P)
    cholesky(σσᵀ).U'
end

function a(t, P::ProkaryoteAux{Val{(true,true,true,true)}})
    _a₃₂ = P.c₅*P.v[2]*(P.v[2]-1.0)+2.0*P.c₆*P.v[3]
    a₄₄ = P.c₁*P.v[3]*P.v[4]+P.c₂*(P.K-P.v[4])

    @SMatrix [P.c₃*P.v[4]+P.c₇*P.v[1]  0.0  0.0  0.0;
              0.0  P.c₄*P.v[1]+2.0*_a₃₂+P.c₈*P.v[2]  -_a₃₂  0.0;
              0.0  -_a₃₂  a₄₄+0.5*_a₃₂  a₄₄;
              0.0  0.0  a₄₄  a₄₄]
end

a(t, x, P::ProkaryoteAux) = a(t, P)
σ(t, x, P::ProkaryoteAux) = σ(t, P)
b(t, x, P::ProkaryoteAux) = B(t, P)*x + β(t, P)

constdiff(::ProkaryoteAux) = true
clone(P::ProkaryoteAux, θ) = ProkaryoteAux(θ..., P.K, P.t, P.u, P.T, P.v, observables(P))
clone(P::ProkaryoteAux, θ, v) = ProkaryoteAux(θ..., P.K, P.t, zero(v), P.T, v, observables(P))
params(P::ProkaryoteAux) = [P.c₁, P.c₂, P.c₃, P.c₄, P.c₅, P.c₆, P.c₇, P.c₈]
dependsOnParams(::ProkaryoteAux) = (1,2,3,4,5,6,7,8)


text:
# ### m14.7t.jl

using Pkg, DrWatson

using MCMCChains
using StanSample
using StatisticalRethinking

Dmat = CSV.read(sr_datadir("islandDistMatrix.csv"), DataFrame);
df = CSV.read(sr_datadir("Kline2.csv"), DataFrame);
df.society = 1:10

stan14_7 = "
functions{
    matrix cov_GPL2(matrix x, real sq_alpha, real sq_rho, real delta) {
        int N = dims(x)[1];
        matrix[N, N] K;
        for (i in 1:(N-1)) {
          K[i, i] = sq_alpha + delta;
          for (j in (i + 1):N) {
            K[i, j] = sq_alpha * exp(-sq_rho * square(x[i,j]) );
            K[j, i] = K[i, j];
          }
        }
        K[N, N] = sq_alpha + delta;
        return K;
    }
}
data{
    int T[10];
    int society[10];
    int P[10];
    matrix[10,10] Dmat;
}
parameters{
    vector[10] k;
    real<lower=0> g;
    real<lower=0> b;
    real<lower=0> a;
    real<lower=0> etasq;
    real<lower=0> rhosq;
}
model{
    vector[10] lambda;
    matrix[10,10] SIGMA;
    rhosq ~ exponential( 0.5 );
    etasq ~ exponential( 2 );
    a ~ exponential( 1 );
    b ~ exponential( 1 );
    g ~ exponential( 1 );
    SIGMA = cov_GPL2(Dmat, etasq, rhosq, 0.01);
    k ~ multi_normal( rep_vector(0,10) , SIGMA );
    for ( i in 1:10 ) {
        lambda[i] = (a * P[i]^b/g) * exp(k[society[i]]);
    }
    T ~ poisson( lambda );
}
";

m14_7s_data = Dict(
    :T => df.total_tools,
    :P => df.population,
    :society => df.society,
    :Dmat => Array(Dmat[:, 2:end])
)

m14_7s = SampleModel("m14.7s", stan14_7)
rc14_7s = stan_sample(m14_7s; data = m14_7s_data)
if success(rc14_7s)
  chns14_7s = read_samples(m14_7s, :mcmcchains)
  chns14_7s |> display

  println()

  read_summary(m14_7s) |> display
end

m14_7s_results = "
       mean   sd  5.5% 94.5% n_eff Rhat4
k[1]  -0.16 0.34 -0.71  0.33   615  1.01
k[2]  -0.02 0.33 -0.55  0.47   550  1.01
k[3]  -0.07 0.32 -0.58  0.40   566  1.01
k[4]   0.35 0.29 -0.05  0.81   557  1.01
k[5]   0.08 0.29 -0.36  0.51   509  1.01
k[6]  -0.38 0.30 -0.88  0.03   585  1.00
k[7]   0.14 0.28 -0.29  0.56   567  1.00
k[8]  -0.21 0.29 -0.66  0.20   526  1.01
k[9]   0.26 0.27 -0.15  0.65   534  1.01
k[10] -0.17 0.36 -0.73  0.36   768  1.00
g      0.61 0.60  0.07  1.70  1321  1.00
b      0.28 0.09  0.14  0.42  1011  1.00
a      1.38 1.05  0.25  3.30  1800  1.00
etasq  0.21 0.22  0.03  0.61   982  1.00
rhosq  1.28 1.56  0.08  4.32  2020  1.00
";

# End of m14_7s.jl


text:
include("Include.jl")

# extra:
using Flux
using Flux: @epochs
using BSON: @load

# load training set -
full_training_data_frame = load_tga_training_data()

# filter -
has_TM_flag = 0
experimental_data_table = filter([:visitid, :TM] => (x, y) -> ((x == 2 || x == 3) && y == has_TM_flag), full_training_data_frame)
input_data = convert.(Float32, Matrix(experimental_data_table[!, 3:14]))
validation_data_array = convert.(Float32, Matrix(experimental_data_table[!, 15:19]))

# what is the size of the system?
(P, D) = size(experimental_data_table)

# How are we going to do this?
δ = 0.90
number_of_samples = 1000
indim = 12
outdim = 5
input_range_array = Array{Float32,2}(undef, indim, 2)
for i ∈ 1:indim
    L = (1 - δ) * minimum(input_data[:, i])
    U = (1 + δ) * maximum(input_data[:, i])
    input_range_array[i, 1] = L
    input_range_array[i, 2] = U
end

# archive -
simulation_dict = Dict()
ℳ = 24

# for now - load a particular model -
m = 2
model_name = "deep_coag_model-L$(m)O-TF-TM-$(has_TM_flag)-V2-V3.bson"
model_file_path = joinpath(_PATH_TO_MODELS, "network-1-epoch-12k", model_name)
@load model_file_path deep_coag_model

# run the model on the unperturbed data -
base_line = sample(deep_coag_model, input_data)

oi = 5
sensitivity_array = Array{Float64,2}(undef, (indim - 1), number_of_samples)
for i ∈ 1:indim-1

    # get the range for this input -
    L = input_range_array[i, 1]
    U = input_range_array[i, 2]

    # initialize space for the simulation data -
    simulation_archive = Array{Matrix{Float32},1}()
    input_archive = Array{Matrix{Float32},1}()

    # how many samples?
    for j ∈ 1:number_of_samples

        # compute a new input for channel i
        r = rand(P)
        v = L .+ r .* (U - L)

        # build the input data set -
        perturbed_input_data = copy(input_data)
        for k ∈ 1:P
            perturbed_input_data[k,i] = v[k]
        end

        # run the models on the pertubed input, collect the output -
        perturbed_output_data = sample(deep_coag_model, perturbed_input_data)
        Δ_out = perturbed_output_data .- base_line
        Δ_in = perturbed_input_data .- input_data
        sensitivity_array[i,j] = var(Δ_out[:,oi])/var(Δ_in[:,i])
    end
end

max_stars_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_issues_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_forks_repo_name: varnerlab/UVM-TopDown-LegacyModel


text:
#
# Test specific error methods
#
using ManifoldsBase
using Test

@testset "Test specific Errors and their format." begin
    e = DomainError(1.0, "Norm not zero.") # a dummy
    e2 = ComponentManifoldError(1, e)

    s1 = sprint(showerror, e)
    s2 = sprint(showerror, e2)
    @test s2 == "At #1: $(s1)"

    e3 = CompositeManifoldError()
    s3 = sprint(showerror, e3)
    @test s3 == "CompositeManifoldError()\n"
    @test length(e3) == 0
    @test isempty(e3)
    @test repr(e3) == "CompositeManifoldError()"

    e4 = CompositeManifoldError([e2])
    @test repr(e4) == "CompositeManifoldError([$(repr(e2)), ])"
    s4 = sprint(showerror, e4)
    @test s4 == "CompositeManifoldError: $(s2)"

    eV = [e2, e2]
    e5 = CompositeManifoldError(eV)
    @test repr(e5) == "CompositeManifoldError([$(repr(e2)), $(repr(e2)), ])"
    s5 = sprint(showerror, e5)
    @test s5 == "CompositeManifoldError: $(s2)\n\n...and $(length(eV)-1) more error(s).\n"
end


text:
using Flux
using MLDataPattern
using Random
using CSV
using DataFrames
using Images
using ImageMagick
using EfficientNet

using Flux: flatten
import MLDatasets

datasetdir(args...) = datadir("datasets", args...)

# -------------------------------------------------------------------------------
# Utility functions
# -------------------------------------------------------------------------------
abstract type Dataset; end

function load(D::Type{<:Dataset}, T = Float32; kwargs...)
    return reshape_dataset.(load_raw(D, T); kwargs...)
end

function reshape_dataset((x,y)::Tuple; labelmap = identity)
    return (reshape_samples(x), reshape_labels(labelmap.(y)))
end

reshape_labels(y) = y
reshape_labels(y::AbstractVector) = Array(reshape(y, 1, :))
reshape_samples(x) = x

function reshape_samples(x::AbstractArray{T, 3}) where T
    return Array(reshape(x, size(x, 1), size(x, 2), 1, size(x, 3)))
end

# -------------------------------------------------------------------------------
# MNIST-like datasets
# -------------------------------------------------------------------------------
abstract type AbstractMNIST <: Dataset end
abstract type MNIST <: AbstractMNIST end

function load_raw(::Type{MNIST}, T)
    train = MLDatasets.MNIST.traindata(T)
    test = MLDatasets.MNIST.testdata(T)

    return train, test
end

abstract type FashionMNIST <: AbstractMNIST end

function load_raw(::Type{FashionMNIST}, T)
    train = MLDatasets.FashionMNIST.traindata(T)
    test = MLDatasets.FashionMNIST.testdata(T)

    return train, test
end

function build_network(::Type{<:AbstractMNIST}; seed = 1234)
    Random.seed!(seed)

    return Chain(
        # First convolution
        Conv((5, 5), 1=>20, stride=(1,1), relu),
        MaxPool((2,2)),

        # Second convolution
        Conv((5, 5), 20=>50, stride=(1,1), relu),
        MaxPool((2,2)),

        flatten,
        Dense(800, 500),
        Dense(500, 1)
    )
end

# -------------------------------------------------------------------------------
# CIFAR dataset
# -------------------------------------------------------------------------------
abstract type AbstractCIFAR <: Dataset end
abstract type CIFAR10 <: AbstractCIFAR end

function load_raw(::Type{CIFAR10}, T)
    train = MLDatasets.CIFAR10.traindata(T)
    test = MLDatasets.CIFAR10.testdata(T)

    return train, test
end

abstract type CIFAR20 <: AbstractCIFAR end

function load_raw(::Type{CIFAR20}, T)
    train = MLDatasets.CIFAR100.traindata(T)[[1,2]]
    test = MLDatasets.CIFAR100.testdata(T)[[1,2]]

    return train, test
end

abstract type CIFAR100 <: AbstractCIFAR end

function load_raw(::Type{CIFAR100}, T)
    train = MLDatasets.CIFAR100.traindata(T)[[1,3]]
    test = MLDatasets.CIFAR100.testdata(T)[[1,3]]

    return train, test
end

function build_network(::Type{<:AbstractCIFAR}; seed = 1234)
    Random.seed!(seed)

    return Chain(
        # First convolution, operating upon a 32x32 image
        Conv((3, 3), 3=>64, pad=(1,1), relu),
        MaxPool((2,2)),

        # Second convolution, operating upon a 16x16 image
        Conv((3, 3), 64=>128, pad=(1,1), relu),
        MaxPool((2,2)),

        # Third convolution, operating upon a 4x4 image
        Conv((3, 3), 128=>128, pad=(1,1), relu),
        MaxPool((2,2)),

        flatten,
        Dense(2048, 1)
    )
end

# -------------------------------------------------------------------------------
# SVHN2 dataset
# -------------------------------------------------------------------------------
abstract type AbstractSVHN2 <: Dataset end
abstract type SVHN2 <: AbstractSVHN2 end

function load_raw(::Type{SVHN2}, T)
    train = MLDatasets.SVHN2.traindata(T)
    test = MLDatasets.SVHN2.testdata(T)

    return train, test
end

abstract type SVHN2Full <: AbstractSVHN2 end

function load_raw(::Type{SVHN2Full}, T)
    x1, y1 = MLDatasets.SVHN2.traindata(T)
    x2, y2 = MLDatasets.SVHN2.extradata(T)

    train = (cat(x1, x2; dims = ndims(x1)), cat(y1, y2; dims = ndims(y1)))
    test = MLDatasets.SVHN2.testdata(T)

    return train, test
end

function build_network(::Type{<:AbstractSVHN2}; seed = 1234)
    Random.seed!(seed)

    return Chain(
        # First convolution, operating upon a 32x32 image
        Conv((3, 3), 3=>64, pad=(1,1), relu),
        MaxPool((2,2)),

        # Second convolution, operating upon a 16x16 image
        Conv((3, 3), 64=>128, pad=(1,1), relu),
        MaxPool((2,2)),

        # Third convolution, operating upon a 4x4 image
        Conv((3, 3), 128=>128, pad=(1,1), relu),
        MaxPool((2,2)),

        flatten,
        Dense(2048, 1)
    )
end

# -------------------------------------------------------------------------------
# ImageNet dataset
# -------------------------------------------------------------------------------
struct ImageFiles{T<:Real, W, H, C}
    filenames::Vector{String}

    function ImageFiles(filenames; T = Float32)
        C, W, H = size(Images.channelview(Images.load(filenames[1])))
        return new{T, W, H, C}(filenames)
    end
end

function Base.show(io::IO, d::ImageFiles{T, W, H, C}) where {T, W, H, C}
    sz = (W, H, C, length(d.filenames))
    print(io, join(sz, "x"), " ImageArray{$T}")
end

Base.length(d::ImageFiles) = length(d.filenames)
Base.ndims(::ImageFiles) = 4
Base.getindex(d::ImageFiles, ind::Int) = getindex(d, [ind])

function Base.getindex(d::ImageFiles{T, W, H, C}, inds) where {T, W, H, C}
    X = Array{T}(undef, W, H, C, length(inds))

    for (j, i) in enumerate(inds)
        img = Images.load(d.filenames[i])
        X[:, :, :, j] .= PermutedDimsArray(channelview(img), (2,3,1))
    end
    return X
end

function getdim(imgs::ImageFiles, d::Integer, inds)
    d == ndims(imgs) || error("not supported slicing")
    return imgs[inds]
end

reshape_samples(imgs::ImageFiles) = imgs

abstract type ImageNet <: Dataset end

function load_raw(::Type{ImageNet}, T)
    d_train = CSV.read(datasetdir("ImageNet224", "data_train.csv"), DataFrame)
    d_test = CSV.read(datasetdir("ImageNet224", "data_test.csv"), DataFrame)

    train = (ImageFiles(d_train.files; T), d_train.labels)
    test = (ImageFiles(d_test.files; T), d_test.labels)
    return train, test
end

function build_network(::Type{<:ImageNet}; seed = 1234)
    Random.seed!(seed)
    model = EfficientNet.from_pretrained("efficientnet-b0")

    return EffNet(
        model.stem,
        model.blocks,
        model.head,
        model.pooling,
        Flux.Chain(model.top[1], Dense(1280, 1)),
        model.drop_connect,
        model.model_name,
    )
end

# -------------------------------------------------------------------------------
# Molecules dataset
# -------------------------------------------------------------------------------
abstract type Molecules <: Dataset end

function load_raw(::Type{Molecules}, T)
    d = BSON.load(datasetdir("molecules.bson"))
    train = (T.(Array(d[:train][:x]')), d[:train][:y])
    test = (T.(Array(d[:test][:x]')), d[:test][:y])

    return train, test
end

function load_activity(::Type{Molecules}, T)
    d = BSON.load(datasetdir("molecules.bson"))
    return T.(d[:train][:activity]), T.(d[:test][:activity])
end

function build_network(::Type{Molecules}; seed = 1234)
    Random.seed!(seed)

    return Chain(
        Dense(100, 50, sigmoid),
        BatchNorm(50, sigmoid),
        Dense(50, 25, sigmoid),
        BatchNorm(25, sigmoid),
        Dense(25, 1)
    )
end


text:
include("./pure.jl")
using .pure
using Test, Lazy, LazyJSON

const io = open("mode1.json", "r")
const mode1s = read(io)
const mode1 = LazyJSON.value(mode1s)
close(io)
const mode1_illusts = mode1["illusts"]

@test split_backslash_last("one/two/three") == "three"
@test generate_filepath("myfile") == "$(homedir())/Downloads/myfile"
@test prefix_filename("one/two/three.png", "four", 1) == "001_four.png"
@test prefix_artist_name("name", 1) == "01\n                   name"

const fnm = [find_number_map(x, y)
             for y in 1:6
             for x in 1:5]
@test fnm == [i for i in 0:29]

# Actually not broken
@test_broken print_multiple_imgs(mode1_illusts)  # #23 has 8 pages

@test url_given_size(mode1["illusts"][1], "medium") == "https://i.pximg.net/c/540x540_70/img-master/img/2020/05/14/06/45/24/81547984_p0_master1200.jpg"
@test post_title(mode1_illusts, 1) == "みこっちゃん"

@test medium_urls(mode1_illusts) == ["https://i.pximg.net/c/540x540_10_webp/img-master/img/2020/05/14/06/45/24/81547984_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2020/05/12/06/36/27/81501385_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2020/05/10/23/10/38/81468125_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2020/05/09/04/31/04/81416496_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2020/05/07/05/48/03/81368866_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2020/05/04/05/17/46/81276257_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2020/04/20/04/52/26/80923496_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2020/04/11/01/49/35/80701898_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2020/03/10/04/07/08/80017594_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2020/02/29/19/09/35/79799236_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2020/02/22/21/21/50/79658392_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2020/02/17/05/31/08/79549991_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2020/01/10/21/50/56/78823485_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2020/01/01/02/50/56/78628383_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2019/12/21/20/13/12/78403815_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2019/12/20/03/00/33/78378594_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2019/12/08/19/09/16/78201587_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2019/11/14/00/47/58/77804404_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2019/10/31/06/25/38/77565309_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2019/10/24/22/50/31/77460464_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2019/10/18/06/06/27/77347697_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2019/10/02/07/26/09/77068750_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2019/09/09/04/32/38/76695217_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2019/09/01/01/39/35/76561671_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2019/08/08/04/11/57/76138362_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2019/07/27/16/34/45/75933779_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2019/07/20/02/29/47/75810852_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2019/07/13/17/12/09/75698678_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2019/07/06/17/40/44/75579060_p0_square1200.jpg", "https://i.pximg.net/c/540x540_10_webp/img-master/img/2019/06/29/09/44/52/75457783_p0_square1200.jpg"]

@test post_titles_in_page(mode1_illusts) == ["みこっちゃん", "おりじなる", "0510", "5.3", "おりじなる", "ミコ誕オメ画！", "5.2", "5.1", "310", "Midnight Sun", "222", "バレンタイン", "祝！！！", "あけましておめでとうございます", "ミコサンタ", "C97告知", "ミコバニー", "たちかわ楽市2019仕様4人組", "ハロミコ", "夏服", "御坂美琴写真集１０用", "常盤台中学指定体操服改", "ツイッターまとめ", "スクミズミコクロ", "ミズミコ", "ミコニャン", "とある画帖", "御坂美琴写真集９", "ジャンプ！", "シャワミコ"]

@test page_urls_in_post(mode1_illusts[23]) == (8, ["https://i.pximg.net/c/540x540_70/img-master/img/2019/09/09/04/32/38/76695217_p0_master1200.jpg", "https://i.pximg.net/c/540x540_70/img-master/img/2019/09/09/04/32/38/76695217_p1_master1200.jpg", "https://i.pximg.net/c/540x540_70/img-master/img/2019/09/09/04/32/38/76695217_p2_master1200.jpg", "https://i.pximg.net/c/540x540_70/img-master/img/2019/09/09/04/32/38/76695217_p3_master1200.jpg", "https://i.pximg.net/c/540x540_70/img-master/img/2019/09/09/04/32/38/76695217_p4_master1200.jpg", "https://i.pximg.net/c/540x540_70/img-master/img/2019/09/09/04/32/38/76695217_p5_master1200.jpg", "https://i.pximg.net/c/540x540_70/img-master/img/2019/09/09/04/32/38/76695217_p6_master1200.jpg", "https://i.pximg.net/c/540x540_70/img-master/img/2019/09/09/04/32/38/76695217_p7_master1200.jpg"])

@test change_url_to_full("https://i.pximg.net/c/540x540_70/img-master/img/2019/09/09/04/32/38/76695217_p0_master1200.jpg") == "https://i.pximg.net/c/540x540_70/img-master/img/2019/09/09/04/32/38/76695217_p0.jpg"

@test process_user_url("https://www.pixiv.net/en/users/2232374") == ("2232374", "1")
@test process_user_url("2232374") == ("2232374", "1")
@test process_artwork_url("https://www.pixiv.net/en/artworks/76695217") == ("76695217", "2")
@test process_artwork_url("76695217") == ("76695217", "2")


text:
###################################################
# METODO DO GRADIENTE COM BUSCA LINEAR ARMIJO + INTERPOLAÇÃO QUADRÁTICA + SALVAGUARDAS
# Veja
#   leonardosecchin.github.io/files/otim1/4.1.Metodos_descida_gerais.pdf
#
# Aqui, a busca linear segue o proposto na implementação de referência do SPG:
#   www.ime.usp.br/~egbirgin/tango/codes.php
#
# Autor: Leonardo D. Secchin
# Data : 06/03/2021
#
# Exemplos de uso:
#   gradiente_interp(nlp);
#   gradiente_interp(nlp, x0=[1;1]);
#   gradiente_interp(nlp, eps=1e-8);
#   x, f, gsupn, iter, status = gradiente_interp(nlp, x0=[1;1], eps=1e-8, maxiter=2000, saidas=true);
# onde 'nlp' é a estrutura MathOptNLPModel do problema.
###################################################

using Printf, LinearAlgebra, NLPModels

# FUNÇÃO PRINCIPAL
function gradiente_interp(nlp; x0=nothing, eps=1.0e-6, maxiter=10000, saidas=true)

    # DADOS DE ENTRADA
    # nlp     : estrutura MathOptNLPModel do problema
    # x0      : ponto inicial (opcional)
    # eps     : precisão para convergencia (opcional, padrão=1.0e-6)
    # maxiter : número máximo de iterações (opcional, padrão=1000)
    # saidas  : mensagens na tela?

    # DADOS DE SAÍDA
    # x       : último iterando
    # f       : f(x)
    # gsupn   : norma do supremo do gradiente de f
    # iter    : número de iterações
    # status  : 0=sucesso, 1=falha

    iter   = 0
    status = 1

    # parâmetro busca linear inexata (Armijo)
    eta = 1e-4

    # define ponto inicial caso não fornecido
    if x0 == nothing
        # captura ponto inicial da estrutura nlp
        if nlp.meta.x0 != nothing
            x0 = nlp.meta.x0
        else
            # caso não exista, seta x0=origem
            x0 = zeros(Float64,nlp.meta.nvar)
        end
    end

    # f, gradiente e norma
    x     = float(x0)
    f     = obj(nlp, x)
    g     = grad(nlp, x)
    gsupn = norm(g, Inf)

    # imprimi cabeçalho saídas
    if saidas
        @printf("\nit     \tf         |grad|\n==========================")
        @printf("\n%d\t%8.2e  %8.2e", iter, f, gsupn)
    end

    while gsupn > eps && iter < maxiter
        # direção
        d = -g

        # retorna novo iterando após busca linear (Armijo + interpolação quadrática)
        x, f = buscalinear(nlp, x, f, f, g, d, eta)

        # atualiza dados do iterando
        g     = grad(nlp, x)
        gsupn = norm(g, Inf)

        iter += 1

        # imprimi iteração corrente
        if saidas
            if mod(iter,20) == 0
                @printf("\n\nit     \tf         |grad|\n==========================")
            end
            @printf("\n%d\t%8.2e  %8.2e", iter, f, gsupn)
        end
    end

    # status de sucesso
    if gsupn <= eps
        status = 0
    end

    # saída
    if saidas
        println("\n\n*******************************")
        if status == 0
            println("PROBLEMA RESOLVIDO COM SUCESSO!")
        else
            println("FALHA NA RESOLUÇÃO.")
        end
        println("*******************************\n")
    end

    return x, f, gsupn, iter, status
end

# BUSCA LINEAR COM INTERPOLAÇÃO QUADRÁTICA E SALVAGUARDAS
#
# Código implementado como no SPG disponível em
#   www.ime.usp.br/~egbirgin/tango/codes.php
#
# Aqui, a interpolação quadrática é tentada sucessivas vezes nos pontos x+t*d.
# Salvaguardas são empregadas, isto é, caso o passo da
# interpolação seja próximo de 0 ou 1, usa backtracking (divide por 2).
# Essa estratégia une a aceleração da interpolação quadrática
# com um melhor controle de passos extremos.
#
# Ao calcular o passo da interpolação no ponto x+t*d, aparece um t^2.
# Em
#   leonardosecchin.github.io/files/otim1/4.1.Metodos_descida_gerais.pdf
# só calculamos a interpolação referente ao ponto x+d. Você pode
# refazer as contas dessa aula e verificar que aparece o t^2 quando o ponto
# de referência é x+t*d.
#
# A variável "fmax" neste código será igual à "f". No SPG, "fmax" receberá o
# máximo dos últimos valores de f, enquanto "f" receberá o valor de f no ponto
# corrente.
function buscalinear(nlp, x, f, fmax, g, d, eta)

    # calcula g' * d
    gtd = g' * d

    # passo inicial
    t = 1.0

    xnew = x + d
    fnew = obj(nlp, xnew)

    # enquanto Armijo não é satisfeito, atualiza t...
    while fnew > fmax + t*eta*gtd
        # se t for próximo de 0, usa backtracking
        if t <= 0.1
            t /= 2.0
        else
            # passo da interpolação quadrática
            tquad = - 0.5*gtd*(t^2) / (fnew - f - gtd)

            # se tquad for próximo de 1, usa backtracking
            if (tquad < 0.1) || (tquad > 0.9)
                tquad = t / 2.0
            end

            t = tquad
        end

        # nova tentativa
        xnew = x + t*d
        fnew = obj(nlp, xnew)
    end

    # retorna novo iterando
    return xnew, fnew
end

max_stars_repo_name: leonardosecchin/temp

max_issues_repo_name: leonardosecchin/temp

max_forks_repo_name: leonardosecchin/temp


text:
module Strand1994

import sbp
using sbp.Operator: OperatorData, to_matrix

default_path = string(@__DIR__, "/")

function load_description(;path=""::String, order=2::Int32, verbose=false::Bool)
        if path == ""
                path = default_path
        end

        if order==2
                desc_D = sbp.Operator.read_operator(string(path, "D_21.txt"))
                desc_H = sbp.Operator.read_operator(string(path, "H_21.txt"))
        elseif order==4
                desc_D = sbp.Operator.read_operator(string(path, "D_42.txt"))
                desc_H = sbp.Operator.read_operator(string(path, "H_42.txt"))
        else
           error("No operators implemented for order = ", order)
        end

        return desc_D, desc_H
end

function build_operators(desc_D::OperatorData, desc_H::OperatorData, m::Int64;
                                               is_sparse=true::Bool)
        x, h = sbp.Grid.grid_xp(m)

        D = to_matrix(desc_D, m, m, is_sparse) / h
        H = to_matrix(desc_H, m, m, is_sparse) * h

        return x, h, D, H

end

end

max_stars_repo_name: ooreilly/sbpjl

max_issues_repo_name: ooreilly/sbpjl

max_forks_repo_name: ooreilly/sbpjl


text:
"Parameters for right-hand side function of Specialised Partitioned Additive Runge-Kutta methods for Variational systems."
mutable struct AbstractParametersVSPARK{IT,DT,TT,D,S,R,P,FT,PT,UT,GT,ϕT,tabType} <: Parameters{DT,TT}
    f_f::FT
    f_p::PT
    f_u::UT
    f_g::GT
    f_ϕ::ϕT

    Δt::TT

    tab::tabType

    t::TT
    q::Vector{DT}
    p::Vector{DT}
    λ::Vector{DT}

    function AbstractParametersVSPARK{IT,DT,D,S,R,P}(f_f::FT, f_p::PT, f_u::UT, f_g::GT, f_ϕ::ϕT, Δt::TT, tableau::tabType) where {IT,DT,TT,D,S,R,P,FT,PT,UT,GT,ϕT,tabType}
        # create solution vectors
        q = zeros(DT,D)
        p = zeros(DT,D)
        λ = zeros(DT,D)

        new{IT,DT,TT,D,S,R,P,FT,PT,UT,GT,ϕT,tabType}(f_f, f_p, f_u, f_g, f_ϕ, Δt, tableau, zero(TT), q, p, λ)
    end
end

function update_params!(params::AbstractParametersVSPARK, sol::AtomicSolutionPDAE)
    # set time for nonlinear solver and copy previous solution
    params.t  = sol.t
    params.q .= sol.q
    params.p .= sol.p
    params.λ .= sol.λ
end

equation(int::AbstractIntegratorVSPARK) = int.equation
timestep(int::AbstractIntegratorVSPARK) = int.params.Δt
tableau(int::AbstractIntegratorVSPARK) = int.tableau
pstages(int::AbstractIntegratorVSPARK) = int.tableau.r

function initialize!(int::AbstractIntegratorVSPARK, sol::AtomicSolutionPDAE)
    sol.t̅ = sol.t - timestep(int)

    equation(int).v(sol.t, sol.q, sol.p, sol.v)
    equation(int).f(sol.t, sol.q, sol.p, sol.f)

    initialize!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f,
                            sol.t̅, sol.q̅, sol.p̅, sol.v̅, sol.f̅)
end

function initial_guess!(int::AbstractIntegratorVSPARK, sol::AtomicSolutionPDAE)
    for i in eachstage(int)
        evaluate!(int.iguess, sol.q, sol.p, sol.v, sol.f,
                              sol.q̅, sol.p̅, sol.v̅, sol.f̅,
                              int.cache.q̃, int.cache.p̃, int.cache.ṽ, int.cache.f̃,
                              tableau(int).q.c[i], tableau(int).p.c[i])

        for k in eachdim(int)
            int.solver.x[3*(ndims(int)*(i-1)+k-1)+1] = (int.cache.q̃[k] - sol.q[k])/timestep(int)
            int.solver.x[3*(ndims(int)*(i-1)+k-1)+2] = (int.cache.p̃[k] - sol.p[k])/timestep(int)
            int.solver.x[3*(ndims(int)*(i-1)+k-1)+3] = int.cache.ṽ[k]
        end
    end

    for i in 1:pstages(int)
        evaluate!(int.iguess, sol.q, sol.p, sol.v, sol.f,
                              sol.q̅, sol.p̅, sol.v̅, sol.f̅,
                              int.cache.q̃, int.cache.p̃, int.cache.ṽ, int.cache.f̃,
                              tableau(int).q̃.c[i], tableau(int).p̃.c[i])

        for k in eachdim(int)
            int.solver.x[3*ndims(int)*nstages(int)+3*(ndims(int)*(i-1)+k-1)+1] = (int.cache.q̃[k] - sol.q[k])/timestep(int)
            int.solver.x[3*ndims(int)*nstages(int)+3*(ndims(int)*(i-1)+k-1)+2] = (int.cache.p̃[k] - sol.p[k])/timestep(int)
            int.solver.x[3*ndims(int)*nstages(int)+3*(ndims(int)*(i-1)+k-1)+3] = 0
        end
    end

    if isdefined(tableau(int), :λ) && tableau(int).λ.c[1] == 0
        for k in eachdim(int)
            int.solver.x[3*ndims(int)*nstages(int)+3*(k-1)+3] = int.cache.λ[k]
        end
    end

    if isdefined(tableau(int), :d) && length(tableau(int).d) > 0
        for k in eachdim(int)
            int.solver.x[3*ndims(int)*nstages(int)+3*ndims(int)*pstages(int)+k] = 0
        end
    end
end


text:
using ContinuumArrays, DifferentialEquations, Plots

##
# natural boundary conditions
##

L = LinearSpline(range(-1,1; length=20))
x = axes(L,1)
D = Derivative(x)
M = L'L
Δ = -((D*L)'D*L)
u0 = copy(L \ exp.(x))

heat(u,(M,Δ),t) = M\(Δ*u)
prob = ODEProblem(heat,u0,(0.0,1.0),(M,Δ))
sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)

g = range(-1,1;length=1000)
@gif for t in 0.0:0.01:1.0
    plot(g, L[g,:]*sol(t); ylims=(0,3))
end

##
# Dirichlet boundary conditions
##

S = L[:,2:end-1]
M = S'S
Δ = -((D*S)'D*S)
u0 = (L \ broadcast(x -> (1-x^2)*exp(x), x))[2:end-1]
prob = ODEProblem(heat,u0,(0.0,1.0),(M,Δ))
sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)

g = range(-1,1;length=1000)
@gif for t in 0.0:0.01:1.0
    plot(g, S[g,:]*sol(t); ylims=(0,3))
end


text:
module GmshDiscreteModelsTests

using Test
using GridapGmsh
using Gridap

d = mktempdir()
f = joinpath(d,"model")
#f = "model"

mshfile = joinpath(@__DIR__,"t1.msh")
model = GmshDiscreteModel(mshfile)
writevtk(model,f)

@test filesize("$(f)_0.vtu") > 0
@test filesize("$(f)_1.vtu") > 0
@test filesize("$(f)_2.vtu") > 0

mshfile = joinpath(@__DIR__,"../demo/","demo.msh")
model = GmshDiscreteModel(mshfile)
writevtk(model,f)

@test filesize("$(f)_0.vtu") > 0
@test filesize("$(f)_1.vtu") > 0
@test filesize("$(f)_2.vtu") > 0
@test filesize("$(f)_3.vtu") > 0

rm(d,recursive=true)

mshfile = joinpath(@__DIR__,"../demo/","demo.msh")
model = GmshDiscreteModel(mshfile)

order = 3

fespace = FESpace(
  reffe=:Lagrangian, order=order, valuetype=Float64,
  conformity=:H1, model=model)

U = TrialFESpace(fespace)

u(x) = x[1]*x[2]*x[3]

uh = interpolate(U,u)

trian = Triangulation(model)
quad = CellQuadrature(trian,degree=3*order)

e = u - uh

ce = integrate( inner(e,e), trian, quad)

#writevtk(trian,"trian",celldata=["ce"=>ce],cellfields=["e"=>e,"uh"=>uh,"u"=>CellField(trian,u)])

tol = 1.e-8
r = sum( ce  )
@test r < tol

end # module


text:
const version = "v0.7"

"""
    struct UnstructuredDiscreteModel{Dc,Dp,Tp,B} <: DiscreteModel{Dc,Dp}
      grid::UnstructuredGrid{Dc,Dp,Tp,B}
      grid_topology::UnstructuredGridTopology{Dc,Dp,Tp,B}
      face_labeling::FaceLabeling
    end
"""
struct UnstructuredDiscreteModel{Dc,Dp,Tp,B} <: DiscreteModel{Dc,Dp}
  grid::UnstructuredGrid{Dc,Dp,Tp,B}
  grid_topology::UnstructuredGridTopology{Dc,Dp,Tp,B}
  face_labeling::FaceLabeling
end

"""
    UnstructuredDiscreteModel(grid::Grid)
"""
function UnstructuredDiscreteModel(grid::Grid)
  _grid = UnstructuredGrid(grid)
  topo = UnstructuredGridTopology(_grid)
  nfaces = [num_faces(topo,d) for d in 0:num_cell_dims(topo)]
  labels = FaceLabeling(nfaces)
  UnstructuredDiscreteModel(_grid,topo,labels)
end

"""
"""
function UnstructuredDiscreteModel(model::DiscreteModel)
  grid = UnstructuredGrid(get_grid(model))
  topo = UnstructuredGridTopology(get_grid_topology(model))
  labels = get_face_labeling(model)
  UnstructuredDiscreteModel(grid,topo,labels)
end

function UnstructuredDiscreteModel(model::UnstructuredDiscreteModel)
  model
end

# Implementation of the interface

get_grid(model::UnstructuredDiscreteModel) = model.grid

get_grid_topology(model::UnstructuredDiscreteModel) = model.grid_topology

get_face_labeling(model::UnstructuredDiscreteModel) = model.face_labeling

# Io

function to_dict(model::UnstructuredDiscreteModel)
  dict = Dict{Symbol,Any}()
  grid = get_grid(model)
  labels = get_face_labeling(model)
  dict[:grid] = to_dict(grid)
  dict[:labeling] = to_dict(labels)
  dict[:vertex_node] = get_vertex_node(model)
  dict[:version] = version
  D = num_cell_dims(model)
  dict[:D] = D
  topo = get_grid_topology(model)
  for d in 1:D
    k = Symbol("face_vertices_$d")
    dict[k] = to_dict(get_faces(topo,d,0))
  end
  dict
end

function check_dict(::Type{UnstructuredDiscreteModel},dict::Dict{Symbol,Any})

  if ! haskey(dict,:version)
    @unreachable "Cannot convert Dict to UnstructuredDiscreteModel for legacy format"
  end

  v = VersionNumber(dict[:version])
  if v < VersionNumber(version)
    @unreachable "Cannot convert Dict to UnstructuredDiscreteModel for outdated format"
  end

end

function from_dict(T::Type{UnstructuredDiscreteModel},dict::Dict{Symbol,Any})

  check_dict(T,dict)

  grid = from_dict(UnstructuredGrid,dict[:grid])
  labeling = from_dict(FaceLabeling,dict[:labeling])
  D::Int = dict[:D]
  vertex_to_node::Vector{Int} = dict[:vertex_node]
  orientation = OrientationStyle(grid)
  polytopes = map(get_polytope, get_reffes(grid))
  cell_type = get_cell_type(grid)
  vertex_coordinates = get_node_coordinates(grid)[vertex_to_node]

  nvertices = length(vertex_to_node)
  d_dface_to_vertices = Vector{Table{Int,Vector{Int},Vector{Int32}}}(undef,D+1)
  d_dface_to_vertices[0+1] = identity_table(Int,Int32,nvertices)
  for d in 1:D
    k = Symbol("face_vertices_$d")
    dface_to_vertices = from_dict(Table{Int,Vector{Int},Vector{Int32}},dict[k])
    d_dface_to_vertices[d+1] = dface_to_vertices
  end

  topo = UnstructuredGridTopology(
    vertex_coordinates,
    d_dface_to_vertices,
    cell_type,
    polytopes,
    orientation)

  UnstructuredDiscreteModel(grid,topo,labeling)
end

function simplexify(model::UnstructuredDiscreteModel)

  grid = get_grid(model)
  grid_topology = get_grid_topology(model)
  labels = get_face_labeling(model)

  tgrid = simplexify(grid)
  tgrid_topology = GridTopology(tgrid)

  reffes = get_reffes(grid)
  tfacelabels = _generate_tfacelabels(grid_topology, tgrid_topology, reffes, labels)

  UnstructuredDiscreteModel(tgrid,tgrid_topology,tfacelabels)

end

const UNSET_ID = Int32(0)

function _generate_tfacelabels(grid_topology, tgrid_topology, reffes, facelabels)

  @notimplementedif length(reffes) != 1
  reffe = first(reffes)
  p = get_polytope(reffe)
  ltcell_to_lnodes, simplex = simplexify(p)

  D = num_cell_dims(grid_topology)

  dim_to_tface_to_label = [ fill(UNSET_ID,num_faces(tgrid_topology,d)) for d in 0:D ]
  dim_to_face_to_label = facelabels.d_to_dface_to_entity

  _fill_dim_to_tface_to_label!(
    dim_to_tface_to_label,
    dim_to_face_to_label,
    tgrid_topology,
    grid_topology,
    ltcell_to_lnodes)

  FaceLabeling(
    dim_to_tface_to_label, facelabels.tag_to_entities, facelabels.tag_to_name)

end

function  _fill_dim_to_tface_to_label!(
  dim_to_tface_to_label,
  dim_to_face_to_label,
  tgrid_topology,
  grid_topology,
  ltcell_to_lnodes)

  tpolytope = first(get_polytopes(tgrid_topology))
  polytope = first(get_polytopes(grid_topology))

  D = length(dim_to_face_to_label)-1
  d = 0
  dim_to_tface_to_label[d+1] = dim_to_face_to_label[d+1]

  for d in 1:(D-1)

    cell_to_faces = get_faces(grid_topology,D,d)
    tcell_to_tfaces = get_faces(tgrid_topology,D,d)

    ntfaces = maximum(tcell_to_tfaces.data)

    ltface_to_ltnodes = get_faces(tpolytope,d,0)
    lface_to_lnodes = get_faces(polytope,d,0)

    tface_to_face = _generate_tface_to_face(
      cell_to_faces.data,
      cell_to_faces.ptrs,
      tcell_to_tfaces.data,
      tcell_to_tfaces.ptrs,
      ltcell_to_lnodes,
      ltface_to_ltnodes,
      lface_to_lnodes,
      ntfaces)

    _update_labels!(
      dim_to_tface_to_label[d+1],dim_to_face_to_label[d+1],tface_to_face)

  end

  d = D
  ncells = length(dim_to_face_to_label[d+1])
  nltcells = length(ltcell_to_lnodes)
  tcell_to_cell = _generate_tcell_to_cell(ncells,nltcells)
  _update_labels!(
    dim_to_tface_to_label[d+1],dim_to_face_to_label[d+1],tcell_to_cell)

  for d = 1:(D-1)

    for j in (d+1):D

      dface_to_jfaces = get_faces(tgrid_topology,d,j)
      dface_to_label = dim_to_tface_to_label[d+1]
      jface_to_label = dim_to_tface_to_label[j+1]
      _fix_dface_to_label!(dface_to_label,jface_to_label,dface_to_jfaces)

    end

  end

end

function _fix_dface_to_label!(dface_to_label,jface_to_label,dface_to_jfaces)

  ndfaces = length(dface_to_label)
  @assert ndfaces == length(dface_to_jfaces)

  for dface in 1:ndfaces

    dlabel = dface_to_label[dface]
    if dlabel != UNSET_ID
      continue
    end

    jfaces = dface_to_jfaces[dface]
    for jface in jfaces
      jlabel = jface_to_label[jface]
      if jlabel != UNSET_ID
        dface_to_label[dface] = jlabel
        break
      end
    end

  end

end

function _update_labels!(tface_to_label,face_to_label,tface_to_face)
  for tface in 1:length(tface_to_label)
    face = tface_to_face[tface]
    if face != 0
      tface_to_label[tface] = face_to_label[face]
    end
  end
end

function _generate_tcell_to_cell(ncells,nltcells)
  ntcells = ncells * nltcells
  tcell_to_cell = Vector{Int}(undef,ntcells)
  tcell = 1
  for cell in 1:ncells
    for ltcell in 1:nltcells
      tcell_to_cell[tcell] = cell
      tcell += 1
    end
  end
  tcell_to_cell
end


text:
"""
---yaml
title: Setup
author: Jonah E Havel
...
"""

using FileIO, Flux, Metalhead, Pages
using Metalhead: classify
###

#
TEST_IMG = ""
vgg = VGG19()
img = load(TEST_IMG)
###

classify(vgg, img)

###
using Metalhead
download(CIFAR10)

###

max_stars_repo_name: jempyre/Dionysus

max_issues_repo_name: jempyre/Dionysus

max_forks_repo_name: jempyre/Dionysus


text:
using nanoJulia
using Documenter

makedocs(;
    modules=[nanoJulia],
    authors="godkin1211 <jiajyun.ciou@gmail.com> and contributors",
    repo="https://github.com/godkin1211/nanoJulia.jl/blob/{commit}{path}#L{line}",
    sitename="nanoJulia.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://godkin1211.github.io/nanoJulia.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/godkin1211/nanoJulia.jl",
)


text:
mutable struct EChart
   id::String
   options::Dict{String,Any}
   width::Int64
   height::Int64
end

function EChart(p::Plot)

    id=randstring(10)
    width=800
    height=600

    options_d=dict(p)
    if haskey(options_d,"width")
      width=options_d["width"]
      delete!(options_d,"width")
    end
    if haskey(options_d,"height")
      height=options_d["height"]
      delete!(options_d,"height")
    end

    return EChart(id,options_d,width,height)
end

const BASE_OPTIONS=["title","legend","grid","xAxis","yAxis","radiusAxis","angleAxis","dataZoom","visualMap","tooltip","axisPointer","toolbox","brush","parallel","parallelAxis","singleAxis","timeline",
"graphic","aria","color","backgroundColor","textStyle","animation","width","height"]

function EChart(kind::String,args...;kwargs...)
    data=[]
    for r in args
       push!(data,DataAttrs(r))
    end

    series_options=Dict()
    chart_options=Dict()

    for (k,v) in kwargs
       k=string(k)
       if k in BASE_OPTIONS
           chart_options[k]=v
       else
           series_options[k]=v
       end
    end

    series_name="Series1"
    series=Series(kind,series_name,data,series_options)
    p=Plot([series],chart_options)

    return EChart(p)
end

import Base.getindex
import Base.setindex!

function Base.getindex(ec::EChart, i::String)

    if i=="width"
       return ec.width
    elseif i=="height"
        return ec.height
    else
        return Base.getindex(ec.options, i)
    end
end

function Base.setindex!(ec::EChart, v,i::String)
    if i=="width"
       ec.width=v
    elseif i=="height"
        ec.height=v
    else

      Base.setindex!(ec.options, v,i)
    end
    return nothing
end


text:
#=
    Helper functions to store constraint mappings
=#
cmap(model::LinQuadOptimizer) = model.constraint_mapping
# dummy method for isvalid
constrdict(::LinQuadOptimizer, ::CI{F, S}) where {F, S} = Dict{CI{F, S}, Any}()

function Base.getindex(model::LinQuadOptimizer, index::CI)
    return constrdict(model, index)[index]
end

# Internal function: delete the constraint name.
function delete_constraint_name(model::LinQuadOptimizer, index::CI)
    model.name_to_constraint = nothing
    delete!(model.constraint_to_name, index)
    return
end

"""
    has_integer(model::LinQuadOptimizer)::Bool

A helper function to determine if `model` has any integer components (i.e.
binary, integer, special ordered sets, semicontinuous, or semi-integer
variables).
"""
function has_integer(model::LinQuadOptimizer)
    constraint_map = cmap(model)
    return length(constraint_map.integer) > 0 ||
        length(constraint_map.binary) > 0 ||
        length(constraint_map.sos1) > 0 ||
        length(constraint_map.sos2) > 0 ||
        length(constraint_map.semicontinuous) > 0 ||
        length(constraint_map.semiinteger) > 0
end

#=
    MOI.is_valid
=#

function MOI.is_valid(model::LinQuadOptimizer, index::CI{F,S}) where F where S
    dict = constrdict(model, index)
    return haskey(dict, index)
end

"""
    __assert_valid__(model::LinQuadOptimizer, index::MOI.Index)

Throw an MOI.InvalidIndex error if `MOI.is_valid(model, index) == false`.
"""
function __assert_valid__(model::LinQuadOptimizer, index::MOI.Index)
    if !MOI.is_valid(model, index)
        throw(MOI.InvalidIndex(index))
    end
end

"""
    __assert_supported_constraint__(model::LinQuadOptimizer, ::Type{F}, ::Type{S})

Throw an `UnsupportedConstraint{F, S}` error if the model cannot add constraints
of type `F`-in-`S`.
"""
function __assert_supported_constraint__(model::LinQuadOptimizer, ::Type{F}, ::Type{S}) where {F<:MOI.AbstractFunction, S<:MOI.AbstractSet}
    if !((F,S) in supported_constraints(model))
        throw(MOI.UnsupportedConstraint{F, S}())
    end
end

#=
    Get number of constraints
=#

function MOI.get(model::LinQuadOptimizer, attribute::MOI.NumberOfConstraints{F, S}) where F where S
    __assert_supported_constraint__(model, F, S)
    return length(constrdict(model, CI{F,S}(UInt(0))))
end

#=
    Get list of constraint references
=#

function MOI.get(model::LinQuadOptimizer, attribute::MOI.ListOfConstraintIndices{F, S}) where F where S
    __assert_supported_constraint__(model, F, S)
    dict = constrdict(model, CI{F,S}(UInt(0)))
    indices = collect(keys(dict))
    return sort(indices, by=x->x.value)
end

#=
    Get list of constraint types in model
=#

function MOI.get(model::LinQuadOptimizer, ::MOI.ListOfConstraints)
    return [(F, S) for (F, S) in supported_constraints(model) if
        MOI.get(model, MOI.NumberOfConstraints{F,S}()) > 0]
end

#=
    Get constraint names
=#

function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintName, index::CI)
    return get(model.constraint_to_name, index, "")
end

function MOI.supports(
        ::LinQuadOptimizer, ::MOI.ConstraintName, ::Type{<:MOI.ConstraintIndex})
    return true
end

function MOI.set(model::LinQuadOptimizer, ::MOI.ConstraintName,
                  index::MOI.ConstraintIndex, name::String)
    if name != ""
        model.constraint_to_name[index] = name
    else
        delete!(model.constraint_to_name, index)
    end
    model.name_to_constraint = nothing
    return
end

#=
    Get constraint by name. Rebuild the lazy reverse-lookup if required.
=#

function _rebuild_name_to_constraint(model::LinQuadOptimizer)
    model.name_to_constraint = Dict{String, MOI.ConstraintIndex}()
    for (constr, c_name) in model.constraint_to_name
        c_name == "" && continue
        if haskey(model.name_to_constraint, c_name)
            duuplicate_con = model.name_to_constraint[c_name]
            model.name_to_constraint = nothing
            error("Constraints $(constr) and $(duuplicate_con) have the ",
                  "same name: ", c_name)
        else
            model.name_to_constraint[c_name] = constr
        end
    end
    return
end

function MOI.get(
        model::LinQuadOptimizer, ::Type{<:MOI.ConstraintIndex}, name::String)
    if model.name_to_constraint === nothing
        _rebuild_name_to_constraint(model)
    end
    return get(model.name_to_constraint, name, nothing)
end

# this covers the type-stable get(m, ConstraintIndex{F,S}, name)::CI{F,S} case
function MOI.get(model::LinQuadOptimizer, ::Type{MOI.ConstraintIndex{F,S}},
                 name::String) where {F, S}
    index = MOI.get(model, MOI.ConstraintIndex, name)
    if isa(index, MOI.ConstraintIndex{F, S})
        return index::MOI.ConstraintIndex{F, S}
    end
    return nothing
end

"""
    CSRMatrix{T}

Type alias for `Transpose{T,SparseMatrixCSC{T,Int}} where T`, which allows for the
    efficient storage of compressed sparse row (CSR) sparse matrix

Please use the following interface to access this data in case model is changed in future
 - `row_pointers(mat)` returns that sparse row_pointer vector see description
    which has a length of rows+1
 - `colvals` this is the column version of `rowvals(SparseMatrixCSC)`
 - `row_nonzeros` this is the equivilent to `nonzeros(SparseMatrixCSC)`
"""

const CSRMatrix{T} = Transpose{T,SparseMatrixCSC{T,Int}} where T

function CSRMatrix{T}(row_pointers, columns, coefficients) where T
    @assert length(columns) == length(coefficients)
    @assert length(columns) + 1 == row_pointers[end]
    transpose(SparseMatrixCSC{T,Int}(maximum(columns), length(row_pointers)-1, row_pointers, columns, coefficients))
end

#   Move access to an interface so that if we want to change type definition we can.
colvals(mat::CSRMatrix) = rowvals(transpose(mat))
row_pointers(mat::CSRMatrix) = transpose(mat).colptr  # This is the only way to get at colptrs without recreating it
row_nonzeros(mat::CSRMatrix) = nonzeros(transpose(mat))

#=
    Below we add constraints.
=#

include("constraints/singlevariable.jl")
include("constraints/vectorofvariables.jl")
include("constraints/scalaraffine.jl")
include("constraints/vectoraffine.jl")
include("constraints/scalarquadratic.jl")


text:
# https://sudonull.com/post/7909-Julia-Scripts-and-parsing-command-line-arguments

using DelimitedFiles, SharedArrays, Random, LinearAlgebra, Clustering, Base
using DataFrames, CSV
using Statistics, ArgParse
using ProgressMeter

nmats = 100 # iteraciones para samplear los eigenvalues de las correlaciones
n_kmedoids = 100 # iteraciones para calcular el minimo de kmedoids
n_kmedoids_out = 10 # Output de de las iteraciones del minimo de los kmedoids

############# FUNCTIONS
function max_eigv(gene)
    #returns largest eigenvalue of a randomized matrix generated from data.
    gran = zeros(size(gene))
    for i = 1:size(gene,2)
        gran[:,i] = gene[shuffle(1:end),i]
    end
 #    @time mr = corr_mat(gran)
    # println(mr[[end-1],end])
    mr = cor(transpose(gran))
    # println(mr[[end-1],end])
    return maximum(eigvals(mr))
end
#------------------------

############### MAIN
## JULIA_NUM_THREADS=$(nproc)
s = ArgParseSettings()
@add_arg_table! s begin
    "-c"
        help = "Flag to calculate number of clusters"
        action = :store_true
    "expfile"
        help = "Matriz de expresion (read a clean format exp matrix)"
        arg_type = String
        required = true
end
parsed_args = parse_args(s)
# parsed_args = Dict([("c", true), ("expfile", "Data/clean/Basal-clean.tsv")])
println(parsed_args)
fname = basename(parsed_args["expfile"])
subtype = split(fname,'-')[1]
println("Working on " * subtype)

println("Kamedoids clustering method on Cancer Expression Matrix")
flush(stdout)

Random.seed!(34568)

### INPUT FILES
#! gname s1 s2 ... gstart chromname
expr = readdlm(parsed_args["expfile"],skipstart=1)
#! cut cols: gname, gstart, chromname
gene = convert(Array{Float64,2},expr[:, 2:end-2])

###################
println("First correlation (m)")
@time  m = cor(transpose(gene))
gnames = expr[:,1]
gstarts = expr[:,end-1]
chroms= expr[:,end]
df = DataFrame(src = [], dst = [], src_gstart = Int64[],
    dst_gstart = Int64[], src_chrom = [], dst_chrom = [],
    pearson = Float64[])
@showprogress for i in 1:length(gnames)
    for j in i+1:length(gnames)
        # println(i," ",j," ",gnames[i]," ",gnames[j]," ",
        #     gstarts[i]," ",gstarts[j]," ",chroms[i]," ",
        #     chroms[j]," ", m[i,j])
        row = [gnames[i] gnames[j] gstarts[i] gstarts[j] chroms[i] chroms[j] m[i,j]]
        push!(df, row)
    end
end
odir = "Data/Pearson/"
mkpath(odir)
oname = odir * subtype * "-allchroms.csv"
CSV.write(oname,df)

###################
println("Automatic number of clusters detection for all chromosomes")
if parsed_args["c"] == true
    println("Get Eigenvalues from shuffles")
    l_evmax = SharedArray{Float64, 1}((nmats))
    for i in 1:nmats #doing the surrogate matrices and saving max eigv.
        @time l_evmax[i] = max_eigv(gene)
        println("Iteration: ",i,", max eigenvalue: ",l_evmax[i])
        flush(stdout)
    end
    #! get eigenvalue of gene
    m_ev = filter(x->x>0.0000001, eigvals(m))
    #! Todos los eigenvalues de la correlacion que sean mayores al mayor de las
    #! correlaciones muestreadas; el numero de estos eigenvalues es el numero de cluesters
    #! nc = 13, for Basal-All.txt with nmats=100
    nc = length(filter(x->x>maximum(l_evmax),m_ev))
    println("Number of clusters (nc): ",nc)
end

###################
println("Starting Kmedoids sampling for all chromosomes")
# make m in [0,1] insted [-1,1]
# correalacion solo en valores positivos
dmat = map(x-> 1-abs(x),m)

ens_c = []
for i = 1:n_kmedoids
    push!(ens_c, kmedoids(dmat,nc))
    if (i % n_kmedoids_out) == 0
    	println(i)
    end
end

#! De todas los clustereos cual es el que tiene el costo total minimo
inx_c = findmin(map(x -> x.totalcost, ens_c))[2]
c = ens_c[inx_c]
println("Cluster counts: ",c.counts)

###################
#! Saving
# writedlm("out.tsv", [expr[:,1] c.assignments c.acosts expr[:,end-1] expr[:,end]])
# header = ["gname" "clusterid" "assignmentcost" "gstart" "chromosome"]
# data = [expr[:,1] c.assignments c.costs expr[:,end-1] expr[:,end]]
# writedlm("out.tsv", [header ; data])
df = DataFrame([expr[:,1],c.assignments,c.costs,expr[:,end-1],expr[:,end]])
header = ["gname","clusterid","assignmentcost","gstart","chromosome"]
rename!(df, Symbol.(header)) # Force set header names
odir = "Data/Clustered/" * subtype * "/"
mkpath(odir)
oname = odir * subtype * "-all-clusters.csv"
CSV.write(oname,df)

max_stars_repo_name: josemaz/gene-matrices

max_issues_repo_name: josemaz/gene-matrices

max_forks_repo_name: josemaz/gene-matrices


text:
include("CenterOfMass.jl")
include("setup_parameters.jl")

actuator_limit = 5.0
max_robots = 4

prior = initialize_prior(circle_ps, resolution, interior_q, masses)

all_robots = collect(1:10)
attached_robots = [1,4,7]

@time @show maximize_csqmi_combinations(attached_robots, 2, prior,
  circle_ws, attachment_ws, sigma, interior_q, actuator_limit)

@time @show maximize_csqmi_available_robots(attached_robots, prior,
  circle_ws, attachment_ws, sigma, interior_q, actuator_limit)

@time @show maximize_csqmi_additional_robot(attached_robots, prior,
  circle_ws, attachment_ws, sigma, interior_q, actuator_limit)

println("null tests")
no_robots = Array{Int64,1}()

@time @show maximize_csqmi_combinations(no_robots, 2, prior,
  circle_ws, attachment_ws, sigma, interior_q, actuator_limit)

@time @show maximize_csqmi_available_robots(no_robots, prior,
  circle_ws, attachment_ws, sigma, interior_q, actuator_limit)

@time @show maximize_csqmi_additional_robot(no_robots, prior,
  circle_ws, attachment_ws, sigma, interior_q, actuator_limit)

max_stars_repo_name: mcorah/CenterOfMass

max_issues_repo_name: mcorah/CenterOfMass

max_forks_repo_name: mcorah/CenterOfMass


text:
# Chemfiles.jl, a modern library for chemistry file reading and writing
# Copyright (C) Guillaume Fraux and contributors -- BSD license

export lengths, set_lengths!, angles, set_angles!, cell_matrix, shape,
set_shape!, volume, CellShape, wrap!

type CellShape
    value::lib.chfl_cellshape

"""
The possible shape for an unit cell are:

- ``Chemfiles.ORTHORHOMBIC`` for unit cell with the three angles are 90°
- ``Chemfiles.TRICLINIC`` for unit cell where the three angles may not be 90°
- ``Chemfiles.INFINITE`` for unit cells without boundaries
"""
    function CellShape(value)
        value = lib.chfl_cellshape(value)
        if value in [lib.CHFL_CELL_INFINITE, lib.CHFL_CELL_ORTHORHOMBIC, lib.CHFL_CELL_TRICLINIC]
            return new(value)
        else
            throw(ChemfilesError("Invalid value for conversion to CellShape: $value"))
        end
    end
end

Base.:(==)(x::CellShape, y::CellShape) = x.value == y.value

const ORTHORHOMBIC = CellShape(lib.CHFL_CELL_ORTHORHOMBIC)
const TRICLINIC = CellShape(lib.CHFL_CELL_TRICLINIC)
const INFINITE = CellShape(lib.CHFL_CELL_INFINITE)

"""
Create an ``UnitCell`` from the three lenghts, with all the angles equal to 90°.
"""
function UnitCell(a::Number, b::Number, c::Number)
    handle = lib.chfl_cell(Float64[a, b, c])
    return UnitCell(handle)
end

"""
Create an ``UnitCell`` from the three lenghts and three angles.
"""
function UnitCell(a::Number, b::Number, c::Number, α::Number, β::Number, γ::Number)
    handle = lib.chfl_cell_triclinic(Float64[a, b, c], Float64[α, β, γ])
    return UnitCell(handle)
end

"""
Get a copy of the ``UnitCell`` of a ``frame``.
"""
function UnitCell(frame::Frame)
    handle = lib.chfl_cell_from_frame(frame.handle)
    return UnitCell(handle)
end

"""
Free the allocated memory for the ``Cell`` object.
"""
function _free(cell::UnitCell)
    _check(
        lib.chfl_cell_free(cell.handle)
    )
    return nothing
end

"""
Get the unit ``cell`` volume
"""
function volume(cell::UnitCell)
    result = Ref{Float64}(0)
    _check(
        lib.chfl_cell_volume(cell.handle, result)
    )
    return result[]
end

"""
Get the three ``cell`` lenghts (a, b and c) in angstroms.
"""
function lengths(cell::UnitCell)
    result = Float64[0, 0, 0]
    _check(
        lib.chfl_cell_lengths(cell.handle, result)
    )
    return result
end

"""
Set the ``cell`` lenghts to ``a``, ``b`` and ``c``.

``a``, ``b`` and ``c`` should be in angstroms.
"""
function set_lengths!(cell::UnitCell, a::Real, b::Real, c::Real)
    _check(
        lib.chfl_cell_set_lengths(cell.handle, Float64[a, b, c])
    )
    return nothing
end

"""
Get the three ``cell`` angles (alpha, beta and gamma) in degrees.
"""
function angles(cell::UnitCell)
    result = Float64[0, 0, 0]
    _check(
        lib.chfl_cell_angles(cell.handle, result)
    )
    return result
end

"""
Set the `cell` angles to ``α``, ``β`` and ``γ``.

``α``, ``β`` and ``γ`` should be in degrees.
"""
function set_angles!(cell::UnitCell, α::Real, β::Real, γ::Real)
    _check(
        lib.chfl_cell_set_angles(cell.handle, Float64[α, β, γ])
    )
    return nothing
end

"""
Get the ``cell`` matricial representation, *i.e.* the representation of the
three base vectors as::

        | a_x   b_x   c_x |
        |  0    b_y   c_y |
        |  0     0    c_z |
"""
function cell_matrix(cell::UnitCell)
    matrix = Array{Float64}(3, 3)
    _check(
        lib.chfl_cell_matrix(cell.handle, pointer(matrix))
    )
    return Array{Float64,2}(matrix)
end

"""
Get the ``cell`` shape, as a ``CellShape`` value
"""
function shape(cell::UnitCell)
    result = Ref{lib.chfl_cellshape}(0)
    _check(
        lib.chfl_cell_shape(cell.handle, result)
    )
    return CellShape(result[])
end

"""
Set the ``cell`` shape to the given ``shape``.
"""
function set_shape!(cell::UnitCell, shape::CellShape)
    _check(
        lib.chfl_cell_set_shape(cell.handle, shape.value)
    )
    return nothing
end

"""
Wrap a `vector` in the unit `cell`.
"""
function wrap!(cell::UnitCell, vector::Vector{Float64})
    _check(
        lib.chfl_cell_wrap(cell.handle, vector)
    )
    return vector
end

"""
Make a deep copy of a ``cell``.
"""
function Base.deepcopy(cell::UnitCell)
    handle = lib.chfl_cell_copy(cell.handle)
    return UnitCell(handle)
end


text:
using MathProgBase

using NLPModels, SparseArrays
import NLPModels.increment!
import MathProgBase.SolverInterface

export MathProgNLPModel

mutable struct ModelReader <: MathProgBase.AbstractMathProgSolver
end

mutable struct MathProgModel <: MathProgBase.AbstractMathProgModel
  eval :: Union{MathProgBase.AbstractNLPEvaluator, Nothing}
  numVar :: Int
  numConstr :: Int
  x :: Vector{Float64}
  y :: Vector{Float64}
  lvar :: Vector{Float64}
  uvar :: Vector{Float64}
  lcon :: Vector{Float64}
  ucon :: Vector{Float64}
  sense :: Symbol
  status :: Symbol
end

MathProgBase.NonlinearModel(solver :: ModelReader) = MathProgModel(nothing,
                                                                   0,
                                                                   0,
                                                                   Float64[],
                                                                   Float64[],
                                                                   Float64[],
                                                                   Float64[],
                                                                   Float64[],
                                                                   Float64[],
                                                                   :Min,
                                                                   :Uninitialized)

function MathProgBase.loadproblem!(m :: MathProgModel,
                                   numVar, numConstr,
                                   l, u, lb, ub,
                                   sense,
                                   eval :: MathProgBase.AbstractNLPEvaluator)

  # TODO: :JacVec is not yet available.
  # [:Grad, :Jac, :JacVec, :Hess, :HessVec, :ExprGraph]
  MathProgBase.initialize(eval, [:Grad, :Jac, :Hess, :HessVec, :ExprGraph])
  m.numVar = numVar
  m.numConstr = numConstr
  m.x = zeros(numVar)
  m.y = zeros(numConstr)
  m.eval = eval
  m.lvar = l
  m.uvar = u
  m.lcon = lb
  m.ucon = ub
  m.sense = sense
end

MathProgBase.setwarmstart!(m :: MathProgModel, x) = (m.x = x)
MathProgBase.status(m :: MathProgModel) = m.status
MathProgBase.getsolution(m :: MathProgModel) = m.x
MathProgBase.getobjval(m :: MathProgModel) = MathProgBase.eval_f(m.eval, m.x)

mutable struct MathProgNLPModel <: AbstractNLPModel
  meta :: NLPModelMeta
  mpmodel :: MathProgModel
  counters :: Counters      # Evaluation counters.

  jrows :: Vector{Int}      # Jacobian sparsity pattern.
  jcols :: Vector{Int}
  jvals :: Vector{Float64}  # Room for the constraints Jacobian.

  hrows :: Vector{Int}      # Hessian sparsity pattern.
  hcols :: Vector{Int}
  hvals :: Vector{Float64}  # Room for the Lagrangian Hessian.
end

"""
    MathProgNLPModel(model, name="Generic")

Construct a `MathProgNLPModel` from a `MathProgModel`.
"""
function MathProgNLPModel(mpmodel :: MathProgModel; name :: String="Generic")

  nvar = mpmodel.numVar
  lvar = mpmodel.lvar
  uvar = mpmodel.uvar

  nlin = length(mpmodel.eval.m.linconstr)         # Number of linear constraints.
  nquad = length(mpmodel.eval.m.quadconstr)       # Number of quadratic constraints.
  nnln = length(mpmodel.eval.m.nlpdata.nlconstr)  # Number of nonlinear constraints.
  ncon = mpmodel.numConstr                        # Total number of constraints.
  lcon = mpmodel.lcon
  ucon = mpmodel.ucon

  jrows, jcols = MathProgBase.jac_structure(mpmodel.eval)
  hrows, hcols = MathProgBase.hesslag_structure(mpmodel.eval)
  nnzj = length(jrows)
  nnzh = length(hrows)

  meta = NLPModelMeta(nvar,
                      x0=mpmodel.x,
                      lvar=lvar,
                      uvar=uvar,
                      ncon=ncon,
                      y0=zeros(ncon),
                      lcon=lcon,
                      ucon=ucon,
                      nnzj=nnzj,
                      nnzh=nnzh,
                      lin=collect(1:nlin),  # linear constraints appear first in MPB
                      nln=collect(nlin+1:ncon),
                      minimize=(mpmodel.sense == :Min),
                      islp=MathProgBase.isobjlinear(mpmodel.eval) & (nlin == ncon),
                      name=name,
                      )

  return MathProgNLPModel(meta,
                      mpmodel,
                      Counters(),
                      jrows,
                      jcols,
                      zeros(nnzj),  # jvals
                      hrows,
                      hcols,
                      zeros(nnzh),  # hvals
                      )
end

import Base.show
show(nlp :: MathProgNLPModel) = show(nlp.mpmodel)

function NLPModels.obj(nlp :: MathProgNLPModel, x :: AbstractVector)
  increment!(nlp, :neval_obj)
  return MathProgBase.eval_f(nlp.mpmodel.eval, x)
end

function NLPModels.grad(nlp :: MathProgNLPModel, x :: AbstractVector)
  g = zeros(nlp.meta.nvar)
  return grad!(nlp, x, g)
end

function NLPModels.grad!(nlp :: MathProgNLPModel, x :: AbstractVector, g :: AbstractVector)
  increment!(nlp, :neval_grad)
  MathProgBase.eval_grad_f(nlp.mpmodel.eval, g, x)
  return g
end

function NLPModels.cons(nlp :: MathProgNLPModel, x :: AbstractVector)
  c = zeros(nlp.meta.ncon)
  return cons!(nlp, x, c)
end

function NLPModels.cons!(nlp :: MathProgNLPModel, x :: AbstractVector, c :: AbstractVector)
  increment!(nlp, :neval_cons)
  MathProgBase.eval_g(nlp.mpmodel.eval, c, x)
  return c
end

function NLPModels.jac_coord(nlp :: MathProgNLPModel, x :: AbstractVector)
  increment!(nlp, :neval_jac)
  MathProgBase.eval_jac_g(nlp.mpmodel.eval, nlp.jvals, x)
  return (nlp.jrows, nlp.jcols, nlp.jvals)
end

function NLPModels.jac(nlp :: MathProgNLPModel, x :: AbstractVector)
  return sparse(jac_coord(nlp, x)..., nlp.meta.ncon, nlp.meta.nvar)
end

function NLPModels.jprod(nlp :: MathProgNLPModel, x :: AbstractVector, v :: AbstractVector)
  Jv = zeros(nlp.meta.ncon)
  return jprod!(nlp, x, v, Jv)
end

function NLPModels.jprod!(nlp :: MathProgNLPModel,
                x :: AbstractVector,
                v :: AbstractVector,
                Jv :: AbstractVector)
  nlp.counters.neval_jac -= 1
  increment!(nlp, :neval_jprod)
  Jv[:] = jac(nlp, x) * v
  return Jv
end

function NLPModels.jtprod(nlp :: MathProgNLPModel, x :: AbstractVector, v :: AbstractVector)
  Jtv = zeros(nlp.meta.nvar)
  return jtprod!(nlp, x, v, Jtv)
end

function NLPModels.jtprod!(nlp :: MathProgNLPModel,
                x :: AbstractVector,
                v :: AbstractVector,
                Jtv :: AbstractVector)
  nlp.counters.neval_jac -= 1
  increment!(nlp, :neval_jtprod)
  Jtv[1:nlp.meta.nvar] = jac(nlp, x)' * v
  return Jtv
end

# Uncomment if/when :JacVec becomes available in MPB.
# "Evaluate the Jacobian-vector product at `x`."
# function NLPModels.jprod(nlp :: MathProgNLPModel, x :: AbstractVector, v :: AbstractVector)
#   jv = zeros(nlp.meta.ncon)
#   return jprod!(nlp, x, v, jv)
# end
#
# "Evaluate the Jacobian-vector product at `x` in place."
# function NLPModels.jprod!(nlp :: MathProgNLPModel, x :: AbstractVector, v :: AbstractVector, jv
# :: AbstractVector)
#   increment!(nlp, :neval_jprod)
#   MathProgBase.eval_jac_prod(nlp.mpmodel.eval, jv, x, v)
#   return jv
# end
#
# "Evaluate the transposed-Jacobian-vector product at `x`."
# function NLPModels.jtprod(nlp :: MathProgNLPModel, x :: AbstractVector, v :: AbstractVector)
#   jtv = zeros(nlp.meta.nvar)
#   return jtprod!(nlp, x, v, jtv)
# end
#
# "Evaluate the transposed-Jacobian-vector product at `x` in place."
# function NLPModels.jtprod!(nlp :: MathProgNLPModel, x :: AbstractVector, v :: AbstractVector,
# jtv :: AbstractVector)
#   increment!(nlp, :neval_jtprod)
#   MathProgBase.eval_jac_prod_t(nlp.mpmodel.eval, jtv, x, v)
#   return jtv
# end

function NLPModels.hess_coord(nlp :: MathProgNLPModel, x :: AbstractVector;
    obj_weight :: Float64=1.0, y :: AbstractVector=zeros(nlp.meta.ncon))
  increment!(nlp, :neval_hess)
  MathProgBase.eval_hesslag(nlp.mpmodel.eval, nlp.hvals, x, obj_weight, y)
  return (nlp.hrows, nlp.hcols, nlp.hvals)
end

function NLPModels.hess(nlp :: MathProgNLPModel, x :: AbstractVector;
    obj_weight :: Float64=1.0, y :: AbstractVector=zeros(nlp.meta.ncon))
  return sparse(hess_coord(nlp, x, y=y, obj_weight=obj_weight)..., nlp.meta.nvar, nlp.meta.nvar)
end

function NLPModels.hprod(nlp :: MathProgNLPModel, x :: AbstractVector, v :: AbstractVector;
    obj_weight :: Float64=1.0, y :: AbstractVector=zeros(nlp.meta.ncon))
  hv = zeros(nlp.meta.nvar)
  return hprod!(nlp, x, v, hv, obj_weight=obj_weight, y=y)
end

function NLPModels.hprod!(nlp :: MathProgNLPModel, x :: AbstractVector, v :: AbstractVector,
    hv :: AbstractVector;
    obj_weight :: Float64=1.0, y :: AbstractVector=zeros(nlp.meta.ncon))
  increment!(nlp, :neval_hprod)
  MathProgBase.eval_hesslag_prod(nlp.mpmodel.eval, hv, x, v, obj_weight, y)
  return hv
end


text:
using LinearAlgebra

@inline function forward(::typeof(*), A::Matrix, B::Matrix)
    A * B, function (Δ::Matrix)
        Base.@_inline_meta
        (nothing, Δ * B', A' * Δ)
    end
end

@inline function forward(::typeof(tr), A::Matrix)
    tr(A), function (Δ::Real)
        Base.@_inline_meta
        (nothing, Δ * Matrix(I, size(A)))
    end
end


text:
#=
   Football bet Simple.

   Netica model of a football bet.
   From Neapolitan?
   The Netica model use an utility node (u) and a decision node (accept_bet)
   which I try to model here...

   Distributions of variable weather
   wet        =>   10000  (1.000000)

   Distributions of variable accept_bet
   no         =>    5079  (0.507900)
   yes        =>    4921  (0.492100)

   Distributions of variable result
   melbwins   =>    6027  (0.602700)
   melbloses  =>    3973  (0.397300)

   Distributions of variable u (num:0)
   20.00000 =>    3154  (0.315400)
   40.00000 =>    2873  (0.287300)
   -20.00000 =>    2048  (0.204800)
   -5.00000 =>    1925  (0.192500)

   Distributions of variable u_accept_bet_yes (num:0)
   0.00000 =>    5079  (0.507900)
   40.00000 =>    2873  (0.287300)
   -20.00000 =>    2048  (0.204800)

   [accept_bet,u]
   Distributions of variable (num:0)
   [1, 20]	=>	3154 (0.3154)
   [2, 40]	=>	2873 (0.2873)
   [2, -20]	=>	2048 (0.2048)
   [1, -5]	=>	1925 (0.1925)

   Cf ~/blog/football_bet_simple.blog
      ~/webppl/football_bet_simple.wppl
=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function football_bet_simple()
    wet = 1
    dry = 2

    no = 1
    yes = 2

    melbwins = 1
    melbloses = 2

    weather ~ Categorical([0.30,0.70]) # [wet,dry]
    accept_bet ~ Categorical([0.50,0.50]) # [no,yes]

    result ~ weather == wet ? Categorical([0.60,0.40]) :  # [melbwins,melbloses]
        Categorical([0.25,0.75])

    u ~ Dirac(
        (result == melbwins  && accept_bet==yes) ?  40 :
        (result == melbwins  && accept_bet==no)  ?  20 :
        (result == melbloses && accept_bet==yes) ? -20 :
        (result == melbloses && accept_bet==no)  ? -5 : 0;
    )

    u_accept_bet_yes ~ Dirac(accept_bet == yes ? u : 0)
    u_accept_bet_no  ~ Dirac(accept_bet == no  ? u : 0)

    true ~ Dirac(weather == wet)
    #  true ~ Dirac(result == melbwins)

    return [accept_bet,u]

end

model = football_bet_simple()

num_chains = 4

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(5), 10_000)
chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns, :weather,["wet","dry"])
show_var_dist_pct(chns, :accept_bet,["no","yes"])
show_var_dist_pct(chns, :result,["melbwins","melbloses"])
show_var_dist_pct(chns, :u)
show_var_dist_pct(chns, :u_accept_bet_yes)

chains_params = Turing.MCMCChains.get_sections(chns, :parameters)
genq = generated_quantities(model, chains_params)
println("[accept_bet,u]")
show_var_dist_pct(genq)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using SparseIR
using Test

import PyCall: pyimport, PyNULL, PyVector

sparse_ir = pyimport("sparse_ir")

@testset "basis.FiniteTempBasis" begin
    lambda_ = 10.0
    beta = 1.0
    eps = 1e-5
    wmax = lambda_ / beta
    kernels = [SparseIR.LogisticKernel, SparseIR.RegularizedBoseKernel]

    # test points
    taus = collect(range(0, beta; length=100))
    omegas = collect(range(-wmax, -wmax; length=10))

    for K in kernels, stat in [fermion, boson]
        shift = (stat == fermion ? 1 : 0)
        v = 2 .* [-1, 1, 100, 1000] .+ shift

        k = K(lambda_)
        basis = FiniteTempBasis(k, stat, beta, eps)
        basis_py = sparse_ir.FiniteTempBasis(stat, beta, wmax; eps=eps, kernel=k.o)
        @test basis.beta == beta
        @test all(basis.u(taus) .== basis_py.u(taus))
        @test all(basis.v(omegas) .== basis_py.v(omegas))
        @test all(basis.uhat(v) .== basis_py.uhat(v))
        @test basis.statistics == stat
        @test basis.size == basis_py.size
        @test size(basis) == basis_py.size

        for l in 1:size(basis)
            @test all(basis.u[l](taus) .== basis_py.u.__getitem__(l - 1)(taus))
            @test all(basis.uhat[l](v) .== basis_py.uhat.__getitem__(l - 1)(v))
        end
    end
end

@testset "basis.FiniteTempBases" begin
    lambda_ = 10.0
    beta = 1.0
    eps = 1e-5
    wmax = lambda_ / beta
    basis_f, basis_b = finite_temp_bases(beta, wmax, eps)

    basis_f_ref = FiniteTempBasis(fermion, beta, wmax, eps)
    basis_b_ref = FiniteTempBasis(boson, beta, wmax, eps)

    @test all(basis_f.s == basis_f_ref.s)
    @test all(basis_b.s == basis_b_ref.s)
end

@testset "basis.FiniteTempBasisInt" begin
    lambda = 10
    beta = 1
    basis = FiniteTempBasis(LogisticKernel(lambda), fermion, beta, 1e-5)
    @test basis isa FiniteTempBasis
end

@testset "basis.FiniteTempBasisDefaultEPS" begin
    lambda = 10
    beta = 1
    basis = FiniteTempBasis(LogisticKernel(lambda), fermion, beta)
    @test basis isa FiniteTempBasis
end

@testset "basis.FiniteTempBasisNewConstructor" begin
    lambda = 10
    beta = 1
    basis = FiniteTempBasis(LogisticKernel(lambda), fermion, beta)
    basis2 = FiniteTempBasis(fermion, beta, lambda / beta; kernel=LogisticKernel(lambda))
    basis3 = FiniteTempBasis(fermion, beta, lambda / beta)
    @test all(basis.s .== basis2.s)
    @test all(basis.s .== basis3.s)
end