text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

"""
    BoundingBox{T<:Real}

Axis-aligned three-dimensional bounding box.

```julia
BoundingBox(xmin::T, xmax::T, ymin::T, ymax::T, zmin::T, zmax::T)
BoundingBox(s::Surface{T})
BoundingBox(s::ParametricSurface{T,3}, transform::Transform{T} = identitytransform(T))
BoundingBox(c::CSGTree{T})
BoundingBox(tri::Triangle{T})
BoundingBox(triangles::AbstractVector{Triangle{T}})
BoundingBox(points::AbstractArray{SVector{3,T}})
BoundingBox(la::LensAssembly{T})
```
"""
struct BoundingBox{T<:Real}
    xmin::T
    ymin::T
    zmin::T
    xmax::T
    ymax::T
    zmax::T

    function BoundingBox(xmin::T, xmax::T, ymin::T, ymax::T, zmin::T, zmax::T) where {T<:Real}
        # if anything is NaN then we want to fall back to the infinite bounding box
        # if the input bounds are infinite then we also want to fall back to the infinite bounding box - this is necessary because the signs of the Infinity can sometimes get messed up
        # isinf(-Inf) == true so this works for any infinity
        xmin = isnan(xmin) || isinf(xmin) ? typemin(T) : xmin
        ymin = isnan(ymin) || isinf(ymin) ? typemin(T) : ymin
        zmin = isnan(zmin) || isinf(zmin) ? typemin(T) : zmin
        xmax = isnan(xmax) || isinf(xmax) ? typemax(T) : xmax
        ymax = isnan(ymax) || isinf(ymax) ? typemax(T) : ymax
        zmax = isnan(zmax) || isinf(zmax) ? typemax(T) : zmax
        if xmin <= xmax && ymin <= ymax && zmin <= zmax
            return new{T}(xmin, ymin, zmin, xmax, ymax, zmax)
        else
            throw(ErrorException("Invalid bounding box"))
        end
    end
end
export BoundingBox

function BoundingBox(s::ParametricSurface{T,3}, transform::Transform{T} = identitytransform(T)) where {T<:Real}
    # get the bounding box of a transformed bounding box
    bbox = BoundingBox(s)
    if transform == identitytransform(T)
        return bbox
    else
        p1 = transform * SVector(bbox.xmin, bbox.ymin, bbox.zmin)
        p2 = transform * SVector(bbox.xmin, bbox.ymax, bbox.zmin)
        p3 = transform * SVector(bbox.xmin, bbox.ymax, bbox.zmax)
        p4 = transform * SVector(bbox.xmin, bbox.ymin, bbox.zmax)
        p5 = transform * SVector(bbox.xmax, bbox.ymin, bbox.zmin)
        p6 = transform * SVector(bbox.xmax, bbox.ymax, bbox.zmin)
        p7 = transform * SVector(bbox.xmax, bbox.ymax, bbox.zmax)
        p8 = transform * SVector(bbox.xmax, bbox.ymin, bbox.zmax)
        return BoundingBox(SVector(p1, p2, p3, p4, p5, p6, p7, p8))
    end
end

function BoundingBox(tri::Triangle{T}) where {T<:Real}
    big = fill(typemin(T), SVector{3,T})
    small = fill(typemax(T), SVector{3,T})

    for vert in vertices(tri)
        small = min.(small, vert)
        big = max.(big, vert)
    end

    return BoundingBox(small[1], big[1], small[2], big[2], small[3], big[3])
end

function BoundingBox(triangles::AbstractVector{Triangle{T}}) where {T<:Real}
    big = fill(typemin(T), SVector{3,T})
    small = fill(typemax(T), SVector{3,T})

    for tri in triangles
        for vert in vertices(tri)
            small = min.(small, vert)
            big = max.(big, vert)
        end
    end

    return BoundingBox(small[1], big[1], small[2], big[2], small[3], big[3])
end

function BoundingBox(points::AbstractArray{SVector{3,T}}) where {T<:Real}
    xmax = maximum((x) -> x[1], points)
    ymax = maximum((x) -> x[2], points)
    zmax = maximum((x) -> x[3], points)
    xmin = minimum((x) -> x[1], points)
    ymin = minimum((x) -> x[2], points)
    zmin = minimum((x) -> x[3], points)
    return BoundingBox(xmin, xmax, ymin, ymax, zmin, zmax)
end

function area(a::BoundingBox{T}) where {T<:Real}
    dx = a.xmax - a.xmin
    dy = a.ymax - a.ymin
    dz = a.zmax - a.zmin
    return 2 * (dx * dy + dy * dz + dz * dx)
end

union(::Nothing, b::BoundingBox{T}) where {T<:Real} = b
union(a::BoundingBox{T}, ::Nothing) where {T<:Real} = a

function union(a::BoundingBox{T}, b::BoundingBox{T}) where {T<:Real}
    return BoundingBox(min(a.xmin, b.xmin), max(a.xmax, b.xmax), min(a.ymin, b.ymin), max(a.ymax, b.ymax), min(a.zmin, b.zmin), max(a.zmax, b.zmax))
end

function intersection(a::BoundingBox{T}, b::BoundingBox{T}) where {T<:Real}
    if a.xmax < b.xmin || a.ymax < b.ymin || a.zmax < b.zmin || b.xmax < a.xmin || b.ymax < a.ymin || b.zmax < b.zmin
        @info a, b
        return nothing
    else
        return BoundingBox(max(a.xmin, b.xmin), min(a.xmax, b.xmax), max(a.ymin, b.ymin), min(a.ymax, b.ymax), max(a.zmin, b.zmin), min(a.zmax, b.zmax))
    end
end

"""
    doesintersect(bbox::BoundingBox{T}, r::AbstractRay{T,3}) -> Bool

Tests whether `r` intersects an axis-aligned [`BoundingBox`](@ref), `bbox`.
"""
function doesintersect(a::BoundingBox{T}, r::AbstractRay{T,3}) where {T<:Real}
    if inside(a, origin(r)) || onsurface(a, origin(r))
        return true
    end

    d = direction(r)
    o = origin(r)

    # Infs and zeros get us into all kinds of problems with AutoDiff here..
    # work arounds are possible in all cases, the code just gets messy :(
    # really we are doing this:
    #     dfx = one(T) / d[1]
    #     dfy = one(T) / d[2]
    #     dfz = one(T) / d[3]
    #     t1 = (a.xmin - o[1]) * dfx
    #     t2 = (a.xmax - o[1]) * dfx
    #     t3 = (a.ymin - o[2]) * dfy
    #     t4 = (a.ymax - o[2]) * dfy
    #     t5 = (a.zmin - o[3]) * dfz
    #     t6 = (a.zmax - o[3]) * dfz
    # tmin = max(max(min(t1, t2), min(t3, t4)), min(t5, t6))
    # tmax = min(min(max(t1, t2), max(t3, t4)), max(t5, t6))

    tmin = typemin(T)
    tmax = typemax(T)

    if d[1] != zero(T)
        # don't really understand why the gradient fails here is xmin or xmax is Inf, but this fixes it
        tx1 = isinf(a.xmin) ? sign(d[1]) * a.xmin : (a.xmin - o[1]) / d[1]
        tx2 = isinf(a.xmax) ? sign(d[1]) * a.xmax : (a.xmax - o[1]) / d[1]
    else
        # avoiding divide by zero (when d[1] == 0) to preserve gradients
        tx1 = sign(a.xmin - o[1]) * typemax(T)
        tx2 = sign(a.xmax - o[1]) * typemax(T)
    end
    # avoid min/max on Infs to preserve gradients
    if tx1 > tx2
        if tmin < tx2
            tmin = tx2
        end
        if tmax > tx1
            tmax = tx1
        end
    else
        if tmin < tx1
            tmin = tx1
        end
        if tmax > tx2
            tmax = tx2
        end
    end

    # below uses the same NaN avoidance techniques...
    if d[2] != zero(T)
        ty1 = isinf(a.ymin) ? sign(d[2]) * a.ymin : (a.ymin - o[2]) / d[2]
        ty2 = isinf(a.ymax) ? sign(d[2]) * a.ymax : (a.ymax - o[2]) / d[2]
    else
        ty1 = sign(a.ymin - o[2]) * typemax(T)
        ty2 = sign(a.ymax - o[2]) * typemax(T)
    end
    if ty1 > ty2
        if tmin < ty2
            tmin = ty2
        end
        if tmax > ty1
            tmax = ty1
        end
    else
        if tmin < ty1
            tmin = ty1
        end
        if tmax > ty2
            tmax = ty2
        end
    end

    if d[3] != zero(T)
        tz1 = isinf(a.zmin) ? sign(d[3]) * a.zmin : (a.zmin - o[3]) / d[3]
        tz2 = isinf(a.zmax) ? sign(d[3]) * a.zmax : (a.zmax - o[3]) / d[3]
    else
        tz1 = sign(a.zmin - o[3]) * typemax(T)
        tz2 = sign(a.zmax - o[3]) * typemax(T)
    end
    if tz1 > tz2
        if tmin < tz2
            tmin = tz2
        end
        if tmax > tz1
            tmax = tz1
        end
    else
        if tmin < tz1
            tmin = tz1
        end
        if tmax > tz2
            tmax = tz2
        end
    end

    return tmax >= tmin && tmax > zero(T)
end
export doesintersect

function inside(a::BoundingBox{T}, p::SVector{3,T}) where {T<:Real}
    return a.xmin < p[1] < a.xmax && a.ymin < p[2] < a.ymax && a.zmin < p[3] < a.zmax
end

function onsurface(a::BoundingBox{T}, p::SVector{3,T}) where {T<:Real}
    return ((p[1] === a.xmin || p[1] === a.xmax) && a.ymin < p[2] < a.ymax && a.zmin < p[3] < a.zmax) || (a.xmin < p[1] < a.xmax && (p[2] === a.ymin || p[2] === a.ymax) && a.zmin < p[3] < a.zmax) || (a.xmin < p[1] < a.xmax && a.ymin < p[2] < a.ymax && (p[3] === a.zmin || p[3] === a.zmax))
end

"""
    surfaceintersection(bbox::BoundingBox{T}, r::AbstractRay{T,3}) -> Union{EmptyInterval{T},Interval{T}}

Calculates the intersection of `r` with an axis-aligned [`BoundingBox`](@ref), `bbox`.

Returns an [`EmptyInterval`](@ref) if there is no intersection or an [`Interval`](@ref) if there is one or two intersections.
Note that the uv of the returned intersection is always **0**.
"""
function surfaceintersection(a::BoundingBox{T}, r::AbstractRay{T,3}) where {T<:Real}
    d = direction(r)
    o = origin(r)

    # Infs and zeros get us into all kinds of problems with AutoDiff here..
    # work arounds are possible in all cases, the code just gets messy :(
    # really we are doing this:
    #     dfx = one(T) / d[1]
    #     dfy = one(T) / d[2]
    #     dfz = one(T) / d[3]
    #     t1 = (a.xmin - o[1]) * dfx
    #     t2 = (a.xmax - o[1]) * dfx
    #     t3 = (a.ymin - o[2]) * dfy
    #     t4 = (a.ymax - o[2]) * dfy
    #     t5 = (a.zmin - o[3]) * dfz
    #     t6 = (a.zmax - o[3]) * dfz
    # tmin = max(max(min(t1, t2), min(t3, t4)), min(t5, t6))
    # tmax = min(min(max(t1, t2), max(t3, t4)), max(t5, t6))

    tmin = typemin(T)
    tmax = typemax(T)
    upper_normal = nothing
    lower_normal = nothing

    if d[1] != zero(T)
        # don't really understand why the gradient fails here is xmin or xmax is Inf, but this fixes it
        tx1 = isinf(a.xmin) ? sign(d[1]) * a.xmin : (a.xmin - o[1]) / d[1]
        tx2 = isinf(a.xmax) ? sign(d[1]) * a.xmax : (a.xmax - o[1]) / d[1]
    else
        # avoiding divide by zero (when d[1] == 0) to preserve gradients
        tx1 = sign(a.xmin - o[1]) * typemax(T)
        tx2 = sign(a.xmax - o[1]) * typemax(T)
    end
    # avoid min/max on Infs to preserve gradients
    if tx1 > tx2
        if tmin < tx2
            lower_normal = SVector{3,T}(1, 0, 0)
            tmin = tx2
        end
        if tmax > tx1
            upper_normal = SVector{3,T}(-1, 0, 0)
            tmax = tx1
        end
    else
        if tmin < tx1
            lower_normal = SVector{3,T}(-1, 0, 0)
            tmin = tx1
        end
        if tmax > tx2
            upper_normal = SVector{3,T}(1, 0, 0)
            tmax = tx2
        end
    end

    # below uses the same NaN avoidance techniques...
    if d[2] != zero(T)
        ty1 = isinf(a.ymin) ? sign(d[2]) * a.ymin : (a.ymin - o[2]) / d[2]
        ty2 = isinf(a.ymax) ? sign(d[2]) * a.ymax : (a.ymax - o[2]) / d[2]
    else
        ty1 = sign(a.ymin - o[2]) * typemax(T)
        ty2 = sign(a.ymax - o[2]) * typemax(T)
    end
    if ty1 > ty2
        if tmin < ty2
            lower_normal = SVector{3,T}(0, 1, 0)
            tmin = ty2
        end
        if tmax > ty1
            upper_normal = SVector{3,T}(0, -1, 0)
            tmax = ty1
        end
    else
        if tmin < ty1
            lower_normal = SVector{3,T}(0, -1, 0)
            tmin = ty1
        end
        if tmax > ty2
            upper_normal = SVector{3,T}(0, 1, 0)
            tmax = ty2
        end
    end

    if d[3] != zero(T)
        tz1 = isinf(a.zmin) ? sign(d[3]) * a.zmin : (a.zmin - o[3]) / d[3]
        tz2 = isinf(a.zmax) ? sign(d[3]) * a.zmax : (a.zmax - o[3]) / d[3]
    else
        tz1 = sign(a.zmin - o[3]) * typemax(T)
        tz2 = sign(a.zmax - o[3]) * typemax(T)
    end
    if tz1 > tz2
        if tmin < tz2
            lower_normal = SVector{3,T}(0, 0, 1)
            tmin = tz2
        end
        if tmax > tz1
            upper_normal = SVector{3,T}(0, 0, -1)
            tmax = tz1
        end
    else
        if tmin < tz1
            lower_normal = SVector{3,T}(0, 0, -1)
            tmin = tz1
        end
        if tmax > tz2
            upper_normal = SVector{3,T}(0, 0, 1)
            tmax = tz2
        end
    end

    if !(tmax >= tmin && tmax > zero(T))
        return EmptyInterval(T)
    else
        if tmin <= zero(T)
            if tmax == typemax(T)
                return rayorigininterval(Infinity(T))
            else
                return rayorigininterval(Intersection(tmax, point(r, tmax), upper_normal, zero(T), zero(T), NullInterface(T)))
            end
        else
            lower = Intersection(tmin, point(r, tmin), lower_normal, zero(T), zero(T), NullInterface(T))
            if tmax == typemax(T)
                return positivehalfspace(lower)
            else
                return Interval(lower, Intersection(tmax, point(r, tmax), upper_normal, zero(T), zero(T), NullInterface(T)))
            end
        end
    end
end


text:
using ArgParse,JLD
include("preprocess.jl")

function main(args=ARGS)
    s = ArgParseSettings()
    s.exc_handler = ArgParse.debug_handler
    @add_arg_table s begin
        ("--trainfile"; required=true; help="To save the metadata about trainfile")
        ("--vocabfile"; required=true; help="To be able to use with different vocabulary sizes")
    end
    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true)
    for (k, v) in o
        println("$k => $v")
    end

    @time d = Data(o[:trainfile]; vocabfile=o[:vocabfile])
    ofs = create_offset(d)
    vocabsize = length(d.word_to_index)
    fname = split(o[:trainfile], '/')[end]
    savefile = string(fname, "_", vocabsize, "k", ".jld")
    info("saving to $savefile")
    save(savefile, "filename", o[:trainfile], "sequences", d.sequences, "word_to_index", d.word_to_index ,"offsets", ofs, "index_to_word", d.index_to_word)
end
!isinteractive() && main(ARGS)

ext: jl

max_stars_repo_name: kirnap/wordsubs

max_issues_repo_name: kirnap/wordsubs

max_forks_repo_name: kirnap/wordsubs


text:
using Catalyst, DiffEqBase, ModelingToolkit, Test

using ModelingToolkit: value

@parameters t k1 k2
@variables S(t) I(t) R(t)
rxs = [Reaction(k1, [S,I], [I], [1,1], [2]),
       Reaction(k2, [I], [R]) ]
rs = ReactionSystem(rxs, t, [S,I,R], [k1,k2])

specset = Set([value(S) =>1, value(I) => 2, value(R) => 3])
@test issetequal(specset, speciesmap(rs))

pset = Set([value(k1) => 1, value(k2) => 2])
@test issetequal(pset, paramsmap(rs))

rxs2 = [Reaction(k2, [I], [R], [1], [1]),
        Reaction(k1, [S,I], [I], [1,1], [2])]
rs2 = ReactionSystem(rxs2, t, [R,I,S], [k2,k1])
@test rs == rs2

rs3 = make_empty_network()
@parameters k3 k4
@variables D(t)
addspecies!(rs3, S)
addspecies!(rs3, D)
addparam!(rs3, k3)
addparam!(rs3, k4)
@test issetequal(species(rs3), [S, D])
@test issetequal(params(rs3), [k3, k4])
addreaction!(rs3, Reaction(k3, [S], [D]))
addreaction!(rs3, Reaction(k4, [S,I], [D]))
merge!(rs, rs3)
addspecies!(rs2, S)
addspecies!(rs2, D)
addparam!(rs2, k3)
addparam!(rs2, k4)
addreaction!(rs2, Reaction(k3, [S], [D]))
addreaction!(rs2, Reaction(k4, [S,I], [D]))
@test rs2 == rs

rxs = [Reaction(k1, [S,I], [I], [1,1], [2]),
       Reaction(k2, [I], [R]) ]
rs = ReactionSystem(rxs, t, [S,I,R], [k1,k2])
rs3 = make_empty_network()
addspecies!(rs3, S)
addspecies!(rs3, D)
addparam!(rs3, k3)
addparam!(rs3, k4)
addreaction!(rs3, Reaction(k3, [S], [D]))
addreaction!(rs3, Reaction(k4, [S,I], [D]))
rs4 = merge(rs, rs3)
@test rs2 == rs4

rxs = [Reaction(k1*S, [S,I], [I], [2,3], [2]),
       Reaction(k2*R, [I], [R]) ]
rs = ReactionSystem(rxs, t, [S,I,R], [k1,k2])
deps = dependents(rxs[2], rs)
@test isequal(deps, [R,I])
@test isequal(dependents(rxs[1], rs), dependants(rxs[1], rs))
addspecies!(rs, S)
@test numspecies(rs) == 3
addspecies!(rs, S, disablechecks=true)
@test numspecies(rs) == 4
addparam!(rs, k1)
@test numparams(rs) == 2
addparam!(rs, k1, disablechecks=true)
@test numparams(rs) == 3

rnmat = @reaction_network begin
       α, S + 2I --> 2I
       β, 3I --> 2R + S
   end α β

smat = [1 2 0;
        0 3 0]
pmat = [0 2 0;
        1 0 2]
@test all(smat .== substoichmat(rnmat))
@test all(pmat .== prodstoichmat(rnmat))


text:
using Test

import CryptoGroups: ElGamal, PGroup, Enc, Dec, specialize, ECGroup, generator, <|
import CryptoGroups

import ShuffleProofs: prove, verify, Simulator, gen_shuffle, Verifier, PoSChallenge, Shuffle, shuffle, VShuffleProof, PoSProof

import ShuffleProofs: step, challenge, PoSChallenge, gen_roprg

@enum VState Config Init PermCommit PoSCommit

###
struct HonestVerifier{T} <: Verifier
    challenge::PoSChallenge
end

HonestVerifier(challenge::PoSChallenge) = HonestVerifier{Config}(challenge)
HonestVerifier{T}(verifier::HonestVerifier) where T = HonestVerifier{T}(verifier.challenge)

PoSChallenge(verifier::HonestVerifier{PoSCommit}) = verifier.challenge

step(verifier::HonestVerifier{Config}, proposition::Shuffle) = HonestVerifier{Init}(verifier)
step(verifier::HonestVerifier{Init}, 𝐜) = HonestVerifier{PermCommit}(verifier)
step(verifier::HonestVerifier{PermCommit}, 𝐜̂, t) = HonestVerifier{PoSCommit}(verifier)
#step(verifier::HonestVerifier{PoSCommit}, s) = HonestVerifier{End}(verifier)

challenge(verifier::HonestVerifier{Init}) = (verifier.challenge.𝐡, verifier.challenge.𝐡[1])
challenge(verifier::HonestVerifier{PermCommit}) = verifier.challenge.𝐮
challenge(verifier::HonestVerifier{PoSCommit}) = verifier.challenge.c

function test_prover(g)

    sk = 5
    pk = g^sk

    enc = Enc(pk, g)
    dec = Dec(sk)

    𝐦 = [g^4, g^2, g^3]
    𝐞 = enc(𝐦, [2, 3, 7])

    N = length(𝐞)

    𝐡 = [g^i for i in 2:N+1]

    𝐫′ = [4, 2, 3]
    proposition, secret = gen_shuffle(enc, 𝐞, 𝐫′) # In practice total of random factors can't match as it reveals
    @test verify(proposition, secret)
    @test verify(proposition, sk)

    (; 𝛙) = secret
    (; 𝐞, 𝐞′) = proposition
    @test dec(𝐞)[𝛙] == dec(𝐞′)

    𝐡 = [g^i for i in 2:N+1]
    𝐮 = [3, 4, 5]
    c = 9

    chg = PoSChallenge(𝐡, 𝐮, c)

    verifier = HonestVerifier(chg)

    # Since the group is small
    # chances that at least one group element will point to 1 are large
    #roprg = gen_roprg(UInt8[7]) # 14, 27, 152
    roprg = gen_roprg(reinterpret(UInt8, Int[304])) # 14, 27, 152, 204, 689, 961
    simulator = prove(proposition, secret, verifier; roprg)
    @test verify(simulator)

    roprg = gen_roprg(reinterpret(UInt8, Int[409])) # 14, 27, 152, 204, 689, 961
    simulator2 = shuffle(𝐞, g, pk, verifier; roprg)
    @test verify(simulator2)

    ### Testing proof translation and verification with Verificatum notation written verifier

    (; proof) = simulator
    vproof = VShuffleProof(proof)
    @test verify(proposition, vproof, chg)

    ### To make it easier I need to type vproof
    @test proof == PoSProof(vproof)

end

p = 23
q = 11
G = specialize(PGroup, p, q)
g = G(3)

test_prover(g)

spec = CryptoGroups.MODP_1024
G = specialize(PGroup, spec)
g = G <| generator(spec)

test_prover(g)

spec = CryptoGroups.Curve_P_256
G = specialize(ECGroup, spec; name = :P_256)
g = G <| generator(spec)

test_prover(g)


text:
using Yota
using Base.Iterators
using Statistics
using MLDataUtils
using Distributions
using CUDA
import NNlib
import NNlibCUDA
import ChainRulesCore: rrule, rrule_via_ad, NoTangent, ZeroTangent, @thunk, unthunk
using Tullio, KernelAbstractions, LoopVectorization, CUDAKernels

include("yota_ext.jl")
include("utils.jl")
include("bcast.jl")
include("init.jl")
include("conv.jl")
include("rnn.jl")
include("activations.jl")
include("losses.jl")
include("layers.jl")
include("batchnorm.jl")
include("optim.jl")
include("fit.jl")
include("metrics.jl")
include("devices.jl")
include("cuda.jl")


text:
Base.parent(A::StarAlgebra{<:Group}) = A.object

function StarAlgebras.AlgebraElement(
    χ::AbstractClassFunction,
    RG::StarAlgebra{<:Group},
)
    G = parent(RG)
    @assert G === parent(χ)
    b = basis(RG)

    dim = degree(χ)
    ord = order(Int, G)
    c = dim // ord

    T = Base._return_type(*, Tuple{typeof(c), eltype(χ)})

    x = AlgebraElement(zeros(T, length(b)), RG)

    for (v, cc) in zip(values(χ), conjugacy_classes(χ))
        for g in cc
            x[inv(g)] = c * v
        end
    end
    return x
end

function algebra_elt_from_support(support, RG::StarAlgebra)
    b = basis(RG)
    I = [b[s] for s in support]
    V = fill(1, length(support))
    return AlgebraElement(sparsevec(I, V, length(b)), RG)
end
struct CyclicSubgroups{Gr, GEl}
    group::Gr
    seen::Dict{Int,Set{GEl}}
    min_order::Int
    max_order::Int

    function CyclicSubgroups(G::Group; min_order=1, max_order=order(Int, G))
        seen = Dict{Int, Set{eltype(G)}}()
        return new{typeof(G), eltype(G)}(G, seen, min_order, max_order)
    end
end

Base.eltype(citr::CyclicSubgroups) = valtype(citr.seen)
Base.IteratorSize(::CyclicSubgroups) = Base.SizeUnknown()

function Base.iterate(citr::CyclicSubgroups)
    g, state = iterate(citr.group) # g is identity here
    @assert isone(g)
    if citr.min_order ≤ 1 ≤ citr.max_order
        citr.seen[1] = Set([g])
        return Set([g]), state
    end
    return iterate(citr, state)
end

function Base.iterate(citr::CyclicSubgroups, state)
    k = iterate(citr.group, state)
    k === nothing && return nothing
    g, state = k
    ord = order(Int, g)
    if citr.min_order ≤ ord ≤ citr.max_order
        if !haskey(citr.seen, ord)
            citr.seen[ord] = Set([g])
            return Set(g^i for i in 1:ord), state
        else
            if any(g^i in citr.seen[ord] for i in 1:ord-1)
                return iterate(citr, state)
            else
                push!(citr.seen[ord], g)
                return Set(g^i for i in 1:ord), state
            end
        end
    end
    return iterate(citr, state)
end

function small_idempotents(
    RG::StarAlgebra{<:Group},
    subgroups = CyclicSubgroups(parent(RG), min_order = 2),
)
    return (algebra_elt_from_support(H, RG) // length(H) for H in subgroups)
end

function (χ::AbstractClassFunction)(α::AlgebraElement{<:StarAlgebra{<:Group}})
    @assert parent(χ) === parent(parent(α))
    iszero(α) && return zero(eltype(StarAlgebras.coeffs(α)))*zero(eltype(χ))
    return sum(α(g) * χ(g) for g in supp(α))
end

function minimal_rank_projection(χ::Character, RG::StarAlgebra{<:Group};
    idems = small_idempotents(RG), # idems are AlgebraElements over Rational{Int}
    iters=3
)
    degree(χ) == 1 && return one(RG, Rational{Int}), true

    for µ in idems
        isone(χ(µ)) && µ^2 == µ && return µ, true
    end

    for n in 2:iters
        for (i, elts) in enumerate(Iterators.product(fill(idems, n)...))
            µ = *(elts...)
            isone(χ(µ)) && µ^2 == µ && return µ, true
        end
    end
    @debug "Could not find minimal projection for $χ"
    return one(RG, Rational{Int}), false
end

function minimal_projection_system(
    chars::AbstractVector{<:AbstractClassFunction},
    RG::StarAlgebra{<:Group}
)
    res = fetch.([Threads.@spawn minimal_rank_projection(χ, RG) for χ in chars])

    r1p, simple = first.(res), last.(res) # rp1 are sparse storage

    mps = [isone(µ) ? AlgebraElement(χ, RG) : µ*AlgebraElement(χ, RG)
        for (µ,χ) in zip(r1p, chars)] # dense storage

    return mps, simple
end


text:
# Run `julia --project=@. --color=yes make.jl` in this directory to
# build the documentation. The `Documenter` must be available (i.e.
# installed globally) for this to work.

using Documenter

using DDF

makedocs(; sitename="DDF", format=Documenter.HTML(; prettyurls=false))
deploydocs(; repo="github.com/eschnett/DDF.jl")


text:
@doc doc"""
    Hyperbolic{N} <: Manifold

The hyperbolic space $\mathbb H^n$ represented by $n+1$-Tuples, i.e. in by
vectors in $\mathbb R^{n+1}$ using the Minkowsi metric, i.e.

```math
\mathbb H^n = \Bigl\{x\in\mathbb R^{n+1}
\ \Big|\ \langle x,x \rangle_{\mathrm{M}}= -x_{n+1}^2
+ \displaystyle\sum_{k=1}^n x_k^2 = -1, x_{n+1} > 0\Bigr\},
```

where $\langle\cdot,\cdot\rangle_{\mathrm{M}}$ denotes the [`minkowski_dot`](@ref)
is Minkowski inner product. The tangent space $T_x\mathbb H^n$ is given by

````math
T_x\mathbb H^n \coloneqq \bigl\{
v \in \mathbb R^{n+1} \ \bigl |\ \langle x,v\rangle_{\mathrm{M}} = 0
\bigr\}.
````
The Minkowski inner product inntroduces the [`MinkowskiMetric`](@ref), which is
a Riemannian metric on the tangent bundle $T\mathbb H^n$.

# Constructor

    Hyperbolic(n)

Generate the $\mathbb H^{n}\subset \mathbb R^{n+1}$
"""
struct Hyperbolic{N} <: Manifold end

Hyperbolic(n::Int) = Hyperbolic{n}()

@doc doc"""
    MinkowskiMetric <: LorentzMetric

The Minkowski metric is a [`LorentzMetric`](@ref) with, i.e.

````math
\langle a,b\rangle_{\mathrm{M}} = -a_{n+1}b_{n+1} +
\displaystyle\sum_{k=1}^n a_kb_k.
````
It is also the default metric e.g. for the [`Hyperbolic`](@ref) space.

!!! note
    While the `MinkowskiMetric` itself is not positive definite in the whole embedded space,
    it is positive definite when restricted to a tangent space $T_x\mathcal M$,
    $x\in\mathcal M$, of the [`Hyperbolic`](@ref) space $\mathcal M$.
"""
struct MinkowskiMetric <: LorentzMetric end

"""
    check_manifold_point(M::Hyperbolic, x; kwargs...)

Check whether `x` is a valid point on the [`Hyperbolic`](@ref) `M`, i.e. is a vector with
[`minkowski_dot`](@ref) -1. The tolerance for the last test can be set using the `kwargs...`.
"""
function check_manifold_point(M::Hyperbolic, x; kwargs...)
    if size(x) != representation_size(M)
        return DomainError(
            size(x),
            "The point $(x) does not lie on $(M), since its size is not $(representation_size(M)).",
        )
    end
    if !isapprox(minkowski_dot(x, x), -1.0; kwargs...)
        return DomainError(
            minkowski_dot(x, x),
            "The point $(x) does not lie on $(M) since its Minkowski inner product is not -1.",
        )
    end
    return nothing
end

"""
    check_tangent_vector(M::Hyperbolic, x, v; kwargs... )

Check whether `v` is a tangent vector to `x` on the [`Hyperbolic`](@ref) `M`, i.e.
after [`check_manifold_point`](@ref)`(M,x)`, `v` has to be of same dimension as `x`
and orthogonal to `x` with respect to [`minkowski_dot`](@ref).
The tolerance for the last test can be set using the `kwargs...`.
"""
function check_tangent_vector(M::Hyperbolic, x, v; kwargs...)
    perr = check_manifold_point(M, x)
    perr === nothing || return perr
    if size(v) != representation_size(M)
        return DomainError(
            size(v),
            "The vector $(v) is not a tangent to a point on $M since its size does not match $(representation_size(M)).",
        )
    end
    if !isapprox(minkowski_dot(x, v), 0.0; kwargs...)
        return DomainError(
            abs(minkowski_dot(x, v)),
            "The vector $(v) is not a tangent vector to $(x) on $(M), since it is not orthogonal (with respect to the Minkowski inner product) in the embedding.",
        )
    end
    return nothing
end

@doc doc"""
    distance(M::Hyperbolic, x, y)

Compute the distance on the [`Hyperbolic`](@ref) `M`, which reads

````math
d_{\mathbb H^n}(x,y) = \operatorname{acosh}( - \langle x, y \rangle_{\mathrm{M}}),
````

where $\langle\cdot,\cdot\rangle_{\mathrm{M}}$ denotes the [`minkowski_dot`](@ref).
"""
distance(M::Hyperbolic, x, y) = acosh(max(-minkowski_dot(x, y), 1.0))

@doc doc"""
    exp(M::Hyperbolic, x, v)

Compute the exponential map on the [`Hyperbolic`](@ref) space $\mathbb H^n$ eminating
from `x` towards `v`, which is optionally scaled by `t`. The formula reads

````math
\exp_x v = \cosh(\sqrt{\langle v,v\rangle_{\mathrm{M}}})x
+ \sinh(\sqrt{\langle v,v\rangle_{\mathrm{M}}})\frac{v}{\sqrt{\langle v,v\rangle_{\mathrm{M}}}},
````

where $\langle\cdot,\cdot\rangle_{\mathrm{M}}$ denotes the [`minkowski_dot`](@ref).
"""
exp(::Hyperbolic, ::Any...)

function exp!(M::Hyperbolic, y, x, v)
    vn = sqrt(max(minkowski_dot(v, v), 0.0))
    vn < eps(eltype(x)) && return copyto!(y, x)
    return copyto!(y, cosh(vn) * x + sinh(vn) / vn * v)
end

flat!(M::Hyperbolic, v::CoTFVector, x, w::TFVector) = copyto!(v, w)

@doc doc"""
    injectivity_radius(M::Hyperbolic[, x])

Return the injectivity radius on the [`Hyperbolic`](@ref), which is always $\infty$.
"""
injectivity_radius(H::Hyperbolic, args...) = Inf

@doc doc"""
    inner(M::Hyperbolic, x, v, w)

Compute the Riemannian inner product for two tangent vectors `v` and `w`
from $T_x\mathbb H^n$ of the [`Hyperbolic`](@ref) space $\mathbb H^n$ given by
$\langle w, v \rangle_{\mathrm{M}}$ the [`minkowski_dot`](@ref) Minkowski
inner product on $\mathbb R^{n+1}$.
"""
@inline inner(M::Hyperbolic, x, w, v) = minkowski_dot(w, v)

is_default_metric(::Hyperbolic, ::MinkowskiMetric) = Val(true)

@doc doc"""
    log(M::Hyperbolic, x, y)

Compute the logarithmic map on the [`Hyperbolic`](@ref) space $\mathbb H^n$, the tangent
vector representing the [`geodesic`](@ref) starting from `x`
reaches `y` after time 1 on the [`Hyperbolic`](@ref) space `M`.
The formula reads for $x\neq y$

```math
\log_x y = d_{\mathbb H^n}(x,y)
\frac{y-\langle x,y\rangle_{\mathrm{M}} x}{\lVert y-\langle x,y\rangle_{\mathrm{M}} x \rVert_2}
```
and is zero otherwise.
"""
log(::Hyperbolic, ::Any...)

function log!(M::Hyperbolic, v, x, y)
    scp = minkowski_dot(x, y)
    w = y + scp * x
    wn = sqrt(max(scp .^ 2 - 1, 0.0))
    wn < eps(eltype(x)) && return zero_tangent_vector!(M, v, x)
    v .= acosh(max(1.0, -scp)) / wn .* w
    return v
end

@doc doc"""
    minkowski_dot(a,b)
Compute the Minkowski inner product of two Vectors `a` and `b` of same length
`n+1`, i.e.

````math
\langle a,b\rangle_{\mathrm{M}} = -a_{n+1}b_{n+1} + \displaystyle\sum_{k=1}^n a_kb_k.
````
"""
function minkowski_dot(a::AbstractVector, b::AbstractVector)
    return -a[end] * b[end] + sum(a[1:end-1] .* b[1:end-1])
end

@doc doc"""
    manifold_dimension(H::Hyperbolic)

Return the dimension of the hyperbolic space manifold $\mathbb H^n$, i.e. $n$.
"""
manifold_dimension(::Hyperbolic{N}) where {N} = N

"""
    mean(
        M::Hyperbolic,
        x::AbstractVector,
        [w::AbstractWeights,]
        method = CyclicProximalPointEstimationMethod();
        kwargs...,
    )

Compute the Riemannian [`mean`](@ref mean(M::Manifold, args...)) of `x` on the
[`Hyperbolic`](@ref) space using [`CyclicProximalPointEstimation`](@ref).
"""
mean(::Hyperbolic, ::Any...)

function mean!(M::Hyperbolic, y, x::AbstractVector, w::AbstractVector; kwargs...)
    return mean!(M, y, x, w, CyclicProximalPointEstimation(); kwargs...)
end

@doc doc"""
    project_tangent(M::Hyperbolic, x, v)

Perform an orthogonal projection with respect to the Minkowski inner product of `v` onto
the tangent space at `x` of the [`Hyperbolic`](@ref) space `M`.

The formula reads
````math
w = v + \langle x,v\rangle_{\mathrm{M}} x,
````
where $\langle \cdot, \cdot \rangle_{\mathrm{M}}$ denotes the Minkowski inner
product in the embedding, see [`minkowski_dot`](@ref).
"""
project_tangent(::Hyperbolic, ::Any...)

project_tangent!(::Hyperbolic, w, x, v) = (w .= v .+ minkowski_dot(x, v) .* x)

@doc doc"""
    representation_size(M::Hyperbolic)

Return the representation size on the [`Hyperbolic`](@ref), i.e. for the `n`-diomensional
hyperbolic manifold the dimention of the embedding, i.e. `n+1`.
"""
@generated representation_size(::Hyperbolic{N}) where {N} = (N + 1,)

sharp!(M::Hyperbolic, v::TFVector, x, w::CoTFVector) = copyto!(v, w)

@doc doc"""
    vector_transport_to(M::Hyperbolic, x, v, y, ::ParallelTransport)

Compute the paralllel transport of the `v` from the tangent space at `x` on the
[`Hyperbolic`](@ref) space $\mathbb H^n$ to the tangent at `y` along the [`geodesic`](@ref)
connecting `x` and `y`. The formula reads

````math
P_{y\gets x}(v) = v - \frac{\langle \log_xy,v\rangle_x}{d^2_{\mathbb H^n}(x,y)}
\bigl(\log_xy + \log_yx \bigr).
````
"""
vector_transport_to(::Hyperbolic, ::Any, ::Any, ::Any, ::ParallelTransport)

function vector_transport_to!(M::Hyperbolic, vto, x, v, y, ::ParallelTransport)
    w = log(M, x, y)
    wn = norm(M, x, w)
    wn < eps(eltype(x + y)) && return copyto!(vto, v)
    return copyto!(vto, v - (inner(M, x, w, v) * (w + log(M, y, x)) / wn^2))
end

@doc doc"""
    zero_tangent_vector(M::Hyperbolic, x)

Return the zero vector from the tangent space at `x` of the [`Hyperbolic`](@ref) `M`.
"""
zero_tangent_vector(::HybridArray, ::Any...)

zero_tangent_vector!(M::Hyperbolic, v, x) = fill!(v, 0)


text:
using JLD
# using Pandas
using Base.Test
import Base.convert
using MAT
using ForwardDiff
using Polynomials
import ForwardDiff.DiffBase
using JSON
# using PyPlot

# local
# @everywhere include("/home/genis/cluster_home/model_fitting_julia/gradient_utils.jl")
# @everywhere include("/home/genis/cluster_home/model_fitting_julia/functions_julia.jl")
# path="/home/genis/cluster_home/model_fitting_julia/"

#portatil

path="/home/emma/fit_behaviour_julia/"
include(path*"functions/functions_julia.jl")
include(path*"functions/utils.jl")
include(path*"functions/analysis.jl")
include(path*"functions/parameter_functions.jl")
include(path*"functions/simulations.jl")

SUBJECTS=[24,25,36,37]
#SUBJECTS=[24]
# filename="synthetic_data_DW_Nstim10000.jld"
model_fit="StLogTimeDw"

Tframe=0.5

model_f=Dict(:c2_f=>c2_const,:c4_f=>c4_const,:c6_f=>c6_zeros,:bias_d_f=>bias_d_const,
        :hbias_d_f=>history_bias_d_const,:st_f=>st_logistic_time,:x0_f=>initial_condition_bias_hbias,
        :ro_f=>read_out_perfect_sim,:lambda=>0)

###### PARAM simulations ######

tau=1.0

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

subject=1

for subject in SUBJECTS
#subject=4
    path_load=path*"results/summary/model_fit"*model_fit*"/"
    filename_load=path_load*"rat"*string(subject)*"_05sec.jld"
    data_fit=load(filename_load)

    iminimum=find(x->x==minimum(data_fit["LL_training"]),data_fit["LL_training"])
    param_best=data_fit["param_all"][iminimum][1]
#    param_best = Dict("c2"=> [10, 0.0],"st"=> [1, 0.0, 0.0, 0.0],"sigma_a"=> [1.0],"x0"=> [0.0, 0.0],"bias_d"=> [0.0],"c4"=>[1],"hbias_d"=> [0.0])
	println("param_ebst", param_best)
    PR_fit=data_fit["PR_training"][iminimum][1]
    PR_fit=data_fit["PR_training"][iminimum][1]
    path_data=path*"scripts_rats/data/"
    println(data_fit["ci"])

    hessian = data_fit["hessian"][iminimum][1]

    # data = Dict()
    open(path_data*"processed_data_rat"*string(subject)*"_dataset2_05.json", "r") do f
        global data
        dicttxt = readstring(f)  # file information to string
        data=JSON.parse(dicttxt)  # parse and transform data
    end
    println(keys(data))
    n_stim_10 = data["n_stim_10"]
    n_stim_20 = data["n_stim_20"]

    #Nstim=100
    # stim=data["stim"][1:Nstim,:]
    # choices=data["choices"][1:Nstim]
    # i_sigma=data["i_sigma"][1:Nstim]

    stim=data["stim_10"]
    stim = permutedims(reshape(stim,(10,n_stim_10)))
    choices=data["choices_10"]
    i_sigmas = data["sigma_10"]
    Nframes=10

    tau=1.0
    Tframe=0.5
    T=Tframe*Nframes
    dt=0.04
    NT = T/dt
    # stim=data["stim"]
    # choices=data["choices"]
    # i_sigmas=data["i_sigma"]
    Nstim=50
    Ntrials=500
    PR_sim=zeros(Nstim)
    for istim in 1:Nstim
        internal_noise_sigma=param_best["sigma_a"][1]*randn(Ntrials,Int(200))
        stim_rep=repeat(transpose(stim[istim,:]),Ntrials)
        d,xfinal=simulation_Nframe_general_trace(param_best,stim_rep,internal_noise_sigma,T,tau,model_f)
        PR_sim[istim]=mean( (d+1)/2 )
    end

    save=Dict("PR_sim"=>PR_sim,"PR_fit"=>PR_fit,"hessian"=>hessian)

    path_save=path*"figures/"
    filename_save=path_save*"check_PR_fit_PR_sim"*model_fit*"_subject"*string(subject)*"_05sec.jld"
    JLD.save(filename_save,save)
    json_string = json_string = JSON.json(save)
	open(path_save*"check_PR_fit_PR_sim"*model_fit*"_subject"*string(subject)*"_05sec.json","w") do f
		write(f, json_string)
	print("done")
	end
#
# figure()
# plot(PR_sim,PR_fit[1:Nstim],".")
# title(string(subject)*" "*string(param_best["c4"]) )
# println(subject)
# println(param_best)
# show()
end
#     history_bias=[0,0]
#     Nsigma=6
#     i_sigmas_possible=0:5
#     Tframe=0.2
#
#     BINS=[]
#     XPDF=[]
#     PARAM_BEST=[]
#     for isig in i_sigmas_possible
#         i=find(x->x==isig,i_sigmas)[1]
#         stim_sig=stim[i,:]
#
#         xpdf,bins=pdf_x_time_1stim(stim_sig,Tframe,param_best,model_f,history_bias)
#         push!(BINS,bins)
#         push!(XPDF,xpdf)
#         push!(PARAM_BEST,param_best)
#     end
#
#     save=Dict("BINS"=>BINS,"XPDF"=>XPDF,"PARAM_BEST"=>PARAM_BEST)
#
#     path_save=path*"figures/"
#     filename_save=path_save*"single_trials_model_fit"*model_fit*"_subject"*string(subject)*".jld"
#     JLD.save(filename_save,save)
#
#
#     df=Pandas.DataFrame(save)
#     filename_save=path_save*"single_trials_model_fit"*model_fit*"_subject"*string(subject)*".json"
#     Pandas.to_json(df,filename_save)
# end
#hello

max_stars_repo_name: emmafdz/FinalDegreeProject

max_issues_repo_name: emmafdz/FinalDegreeProject

max_forks_repo_name: emmafdz/FinalDegreeProject


text:
using jFEMTools
import Tensors

mesh = unitSquareMesh2(RectangleCell, (3,3));

# forcing function
rhs(x::Tensors.Vec{2}) = 2*π^2*sin(π*x[1])*sin(π*x[2])

degree = 2;
dim = 2;
element = PoissonVirtualElement(dim,degree);
u = TrialFunction(VEMFunctionSpace(mesh,element))
dofs = DofHandler(mesh, [u]);
operators = VEMOperators(dofs, u;load = rhs);

K = assemble_stiffnessMat(operators);
b = assemble_load(operators);

# Boundary condition
g(x::Tensors.Vec{2}) = sin(π * x[2])*sin(π * x[1]);
dbc = Dirichlet(dofs,u,"boundary",g);
apply!(K,b,dbc);

#Solve
x = K\b

#Plot solution
using Makie
import AbstractPlotting
include("src/plot_recipes.jl")
#popdisplay(AbstractPlotting.PlotDisplay())
#AbstractPlotting.inline!(true)
scene = Scene(resolution = (400, 200))
# Plot approximation
vi = jFEMTools.vertexdofs(dofs, u)
plot!(scene, mesh, color = x[vi])

#Plot exact solution
vv = get_vertices_matrix(mesh)
xx = [g(vv) for vv in mesh.vertices];
plot!(scene, mesh, color = xx)
#display(scene)

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
module Z_laminate_examples
using FinEtools
using FinEtools.AlgoDeforLinearModule

function Z_laminate_u_ss()
    println("""
    % Three-dimensional Elasticity Solution for Uniformly Loaded Cross-ply
    % Laminates and Sandwich Plates
    % Ashraf M. Zenkour, Journal of Sandwich Structures and Materials 2007 9: 213-238
    % DOI: 10.1177/1099636207065675
    """)

    t0 = time()
    # Lamina material parameters
    E1s = 25.0e6*phun("psi")
    E2s = 1.0e6*phun("psi")
    E3s = E2s
    nu12s = nu13s = nu23s = 0.25
    G12s = 0.5e6*phun("psi")
    G13s = G12s
    G23s = 0.2e6*phun("psi")

    a = 200.0*phun("mm") # side of the square plate
    b = 600.0*phun("mm") # side of the square plate
    q0 = 1.0*phun("psi");
    # The below values come from Table 2
    # h = a/4; wc_analytical = 3.65511/(100*E3s*h^3/a^4/q0);
    # h = a/10; wc_analytical = 1.16899/(100*E3s*h^3/a^4/q0);
    # h = a/50; wc_analytical = 0.66675/(100*E3s*h^3/a^4/q0);
    h = a/100; wc_analytical = 0.65071/(100*E3s*h^3/a^4/q0);
    angles =[0,90,0];
    nLayers = length(angles);
    tolerance = 0.0001*h

    # Generate mesh
    na = 10 # number of elements along the side of the plate
    nb = 30 # number of elements along the side of the plate
    xs = collect(linspace(0.0, a, na+1))
    ys = collect(linspace(0.0, b, nb+1))
    ts = h/nLayers*ones(nLayers);# layer thicknesses
    nts= 3*ones(Int, nLayers);# number of elements per layer
    fens,fes = H8layeredplatex(xs, ys, ts, nts)
    fens,fes = H8toH20(fens,fes)

    MR = DeforModelRed3D
    laminamaterial = MatDeforElastOrtho(MR,
    0.0, E1s, E2s, E3s,
    nu12s, nu13s, nu23s,
    G12s, G13s, G23s,
    0.0, 0.0, 0.0)

    function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        rotmat3!(csmatout, angles[fe_label]/180.*pi* [0.;0.;1.]);
    end

    gr = GaussRule(3, 3)

    region = FDataDict("femm"=>FEMMDeforLinear(MR, IntegData(fes, gr), CSys(3, 3, updatecs!), laminamaterial))

    lx0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
    lxa = selectnode(fens, box=[a a -Inf Inf -Inf Inf], inflate=tolerance)
    ly0 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
    lyb = selectnode(fens, box=[-Inf Inf b b -Inf Inf], inflate=tolerance)

    ex02 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>lx0 )
    ex03 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lx0 )
    exa2 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>lxa )
    exa3 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lxa )
    ey01 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>ly0 )
    ey03 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>ly0 )
    eyb1 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lyb )
    eyb3 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lyb )

    bfes = meshboundary(fes)
    ttopl = selectelem(fens, bfes; facing=true, direction = [0.0 0.0 1.0])
    Trac = FDataDict("traction_vector"=>[0.0; 0.0; -q0],
    "femm"=>FEMMBase(IntegData(subset(bfes, ttopl), GaussRule(2, 3))))

    modeldata = FDataDict("fens"=>fens,
    "regions"=>[region],
    "essential_bcs"=>[ex02, ex03, exa2, exa3, ey01, ey03, eyb1, eyb3],
    "traction_bcs"=> [Trac]
    )
    modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

    u = modeldata["u"]
    geom = modeldata["geom"]
    lcenter = selectnode(fens, box=[a/2 a/2  b/2 b/2 -Inf Inf], inflate=tolerance)
    cdis = abs(mean(u.values[lcenter, 3]))
    println("")
    println("Normalized Center deflection: $(cdis/wc_analytical)")

    File =  "Z_laminate_u_ss.vtk"
    vtkexportmesh(File, fes.conn, geom.values, FinEtools.MeshExportModule.H20;
    scalars = [("Layer", fes.label)], vectors = [("displacement", u.values)])
    @async run(`"paraview.exe" $File`)

    println("Done")
    true

end # Z_laminate_u_ss

function allrun()
    println("#####################################################")
    println("# Z_laminate_u_ss ")
    Z_laminate_u_ss()
    return true
end # function allrun

end # module Z_laminate_examples


text:
module State

using PyCall
import ..@pyfunc

#qulacs = PyNULL()

state = PyNULL()
const state_functions = [
    :inner_product,
    :tensor_product,
    :permutate_qubit,
    :drop_qubit,
    :partial_trace,
]

#=
for func in state_functions
    @eval begin
        function $(func)(args...; kwargs...)
            qulacs.state.$(func)(args...; kwargs...)
        end
        export $(func)
    end
end
=#

for func in state_functions
    @eval begin
        @pyfunc state $(func)
        export $(func)
    end
end

function __init__()
    #copy!(qulacs, pyimport("qulacs"))
    copy!(state, pyimport("qulacs.state"))
end

end # module State


text:
using RecipesBase

export eachplot

@userplot EachPlot

@recipe function f(p::EachPlot)
    x,y = p.args
    label --> ""
    # color --> :black
    alpha --> 0.01
    m = Matrix(y)'
    x,m
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antddraggerupload

"""
    'feffery'_antddraggerupload(;kwargs...)

An AntdDraggerUpload component.

Keyword arguments:
- `id` (String; optional)
- `apiUrl` (String; optional)
- `className` (String; optional)
- `directory` (Bool; optional)
- `failedTooltipInfo` (String; optional)
- `fileListMaxLength` (Real; optional)
- `fileMaxSize` (Real; optional)
- `fileTypes` (Array of Strings; optional)
- `hint` (String; optional)
- `lastUploadTaskRecord` (optional): . lastUploadTaskRecord has the following type: lists containing elements 'fileName', 'fileSize', 'completeTimestamp', 'taskStatus', 'taskId'.
Those elements have the following types:
  - `fileName` (String; optional)
  - `fileSize` (Real; optional)
  - `completeTimestamp` (Real; optional)
  - `taskStatus` (String; optional)
  - `taskId` (String; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `locale` (a value equal to: 'zh-cn', 'en-us'; optional)
- `multiple` (Bool; optional)
- `style` (Dict; optional)
- `text` (String; optional)
- `uploadId` (String; optional)
"""
function 'feffery'_antddraggerupload(; kwargs...)
        available_props = Symbol[:id, :apiUrl, :className, :directory, :failedTooltipInfo, :fileListMaxLength, :fileMaxSize, :fileTypes, :hint, :lastUploadTaskRecord, :loading_state, :locale, :multiple, :style, :text, :uploadId]
        wild_props = Symbol[]
        return Component("'feffery'_antddraggerupload", "AntdDraggerUpload", "feffery_antd_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: ReilYoung/feffery-antd-components

max_issues_repo_name: ReilYoung/feffery-antd-components

max_forks_repo_name: ReilYoung/feffery-antd-components


text:
type USTigerProvider <: ShapefileProvider

  kind::AbstractString
  year::Int
  basename::AbstractString
  url::AbstractString

  projection::AbstractString
  polygons::Array
  records::Dict

  function USTigerProvider(; kind::AbstractString="STATE",
                             year::Int=2015,
                             basename::AbstractString="",
                             url::AbstractString="")
    year < 2010 && error("year must be 2010 or greater")
    prov = new(kind, year, basename, url)
    prov.projection = "epsg:4269"  # TIGER data is NAD83 / EPSG:4269
    _localize(prov)
    prov.polygons = _parseshapes(prov)
    prov.records = _parserecords(prov)
    return prov
  end

end

function projection(provider::USTigerProvider)
  return provider.projection
end

# Nationwide kinds
for k in ["AIANNH", "AITSN", "ANRC", "CBSA", "CNECTA", "COASTLINE", "COUNTY", "CSA",
          "FACESMIL", "METDIV", "MIL", "NECTA", "NECTADIV", "PRIMARYROADS", "RAILS",
          "STATE", "TBG", "TTRACT"]
  @eval begin
    function ($(Symbol(k)))(; year::Int=2015)
      basename = "tl_$year\_us_$(lowercase($k))"
      url = "ftp://ftp2.census.gov/geo/tiger/TIGER$year\/$(uppercase($k))/.zip"
      return USTigerProvider(kind=($k), year=year, basename=basename, url=url)
    end
  end
end

# Geoid kinds
for k in ["AREAWATER", "EDGES", "FACES", "FACESAH", "LINEARWATER", "ROADS"]
  @eval begin
    function ($(Symbol(k)))(; year::Int=2015, geoid::AbstractString="01001")
      basename = "tl_$year\_$geoid\_$(lowercase($k))"
      url = "ftp://ftp2.census.gov/geo/tiger/TIGER$year\/$(uppercase($k))/$basename\.zip"
      return USTigerProvider(kind=($k), year=year, basename=basename, url=url)
    end
  end
end

# Statefp kinds
for k in ["AREALM", "BG", "COUSUB", "FACESAL", "PLACE", "POINTLM", "PRISECROADS", "SLDL",
          "SLDU", "TABBLOCK", "TRACT", "UNSD"]
  @eval begin
    function ($(Symbol(k)))(; year::Int=2015, statefp::AbstractString="01")
      basename="tl_$year\_$statefp\_$(lowercase($k))"
      url = "ftp://ftp2.census.gov/geo/tiger/TIGER$year\/$(uppercase($k))/$basename\.zip"
      return USTigerProvider(kind=($k), year=year, basename=basename, url=url)
    end
  end
end

# Summary kinds with 3 resolutions -- these generate functions with names like
# CBSASUMMARY, COUNTYSUMMARY, etc.
for k in ["CBSA", "CD114", "COUNTY", "CSA", "DIVISION", "NATION", "NECTA", "REGION", "STATE"]
  @eval begin
    function ($(Symbol(string(k,"SUMMARY"))))(; year::Int=2014, resolution::AbstractString="20m")
      resolutions = ["500k", "5m", "20m"]
      !(resolution in resolutions) && error("resolution must be one of $resolutions")
      basename = "cb_$year\_us_$(lowercase($k))_$resolution"
      url = "http://www2.census.gov/geo/tiger/GENZ$year\/shp/$basename\.zip"
      return USTigerProvider(kind=($k), year=year, basename=basename, url=url)
    end
  end
end

# Summary kinds with 500k resolution only
for k in ["COUNTY_WITHIN_CD114", "UA10", "ZCTA510"]
  @eval begin
    function ($(Symbol(string(k,"SUMMARY"))))(; year::Int=2014)
      basename = "cb_$year\_us_$(lowercase($k))_500k"
      url = "http://www2.census.gov/geo/tiger/GENZ$year\/shp/$basename\.zip"
      return USTigerProvider(kind=($k), year=year, basename=basename, url=url)
    end
  end
end

# These (non-summary) kinds have no polygon data. Consider ways to roll these in as shapeless types.
#geoidnonmap=["ADDR", "ADDRFEAT", "ADDRFN", "FEATNAMES"]


text:
function MatrixTest(game_type)
    ## RPS
    RPS = game_type([
        (0,0) (-1,1) (1,-1);
        (1,-1) (0,0) (-1,1);
        (-1,1) (1,-1) (0,0)
    ])
    init_strategy = [0.1,0.2,0.7]
    p1 = player(RPS, 1, init_strategy); p2 = player(RPS, 2, init_strategy)
    train_both!(p1, p2, 100_000)
    # p = plot(p1, p2, lw=2);
    # title!(p, "RM RPS NE")
    # display(p)

    RPS_NE_strat = fill(1/3,3)
    @test isapprox(p1.strategy, RPS_NE_strat, atol=0.01)
    @test isapprox(p2.strategy, RPS_NE_strat, atol=0.01)

    p1 = player(RPS, 1, init_strategy)
    p2 = player(RPS, 2, init_strategy)
    train_one!(p1, p2, 10_000)
    # p = plot(p1,p2, lw=2);
    # title!(p, "RM RPS Exploitative")
    # display(p)

    @test isapprox(p1.strategy, Float64[1,0,0], atol=0.01)
    @test p2.strategy == init_strategy

    ## Prisoners Dilemma
    PD = game_type([
        (-1,-1) (-3,0);
        (0,-3) (-2,-2)
    ])
    p1 = player(PD, 1)
    p2 = player(PD, 2)

    train_both!(p1,p2, 10_000)
    # p = plot(p1,p2, lw=2);
    # title!(p, "RM Prisoner's Dilemma")
    # display(p)

    PD_NE_strat = Float64[0,1]
    @test isapprox(p1.strategy, PD_NE_strat, atol=0.01)
    @test isapprox(p2.strategy, PD_NE_strat, atol=0.01)

    # https://sites.math.northwestern.edu/~clark/364/handouts/bimatrix-mixed.pdf
    game = game_type([
        (1,1) (0,0) (0,0);
        (0,0) (0,2) (3,0);
        (0,0) (2,0) (0,3);
    ])
    p1 = player(game, 1)
    p2 = player(game, 2)
    train_both!(p1,p2, 100_000)
    NEs = [[6/11,3/11,2/11], [0,3/5,2/5], [1,0,0]]
    s1 = p1.strategy; s2 = p2.strategy
    @test begin
        ≈(s1, NEs[1], atol=0.01) ||
        ≈(s1, NEs[2], atol=0.01) ||
        ≈(s1, NEs[3], atol=0.01)
    end
    @test begin
        ≈(s2, NEs[1], atol=0.01) ||
        ≈(s2, NEs[2], atol=0.01) ||
        ≈(s2, NEs[3], atol=0.01)
    end

    v = evaluate(p1,p2)
    @test begin
        all( .≈(v,(6/5,6/5), atol=0.01)) ||
        all( .≈(v,(6/11,6/11), atol=0.01)) ||
        all( .≈(v,(1,1), atol=0.01))
    end
end

@testset "RM Matrix" begin MatrixTest(MatrixGame) end

@testset "CFR Matrix" begin MatrixTest(SimpleIIGame) end

max_stars_repo_name: WhiffleFish/HelloCFR

max_issues_repo_name: WhiffleFish/HelloCFR

max_forks_repo_name: WhiffleFish/HelloCFR


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    ProductPartition(p₁, p₂)

A method for partitioning spatial objects using the product of two
partitioning methods `p₁` and `p₂`.
"""
struct ProductPartition{P1,P2} <: PartitionMethod
  p₁::P1
  p₂::P2
end

# general case
function partition(object, method::ProductPartition)
  # individual partition results
  s₁ = indices(partition(object, method.p₁))
  s₂ = indices(partition(object, method.p₂))

  # label-based representation
  l₁ = Vector{Int}(undef, nelements(object))
  l₂ = Vector{Int}(undef, nelements(object))
  for (i, inds) in enumerate(s₁)
    l₁[inds] .= i
  end
  for (i, inds) in enumerate(s₂)
    l₂[inds] .= i
  end

  # product of labels
  counter = 0
  d = Dict{Tuple{Int,Int},Int}()
  l = Vector{Int}(undef, nelements(object))
  for i in 1:nelements(object)
    t = (l₁[i], l₂[i])
    if t ∉ keys(d)
      counter += 1
      d[t] = counter
    end
    l[i] = d[t]
  end

  # return partition using label predicate
  pred(i, j) = l[i] == l[j]
  partition(object, PredicatePartition(pred))
end

# predicate partition method
function partition(object, method::ProductPartition{P1,P2}) where {P1<:PredicatePartitionMethod,
                                                                   P2<:PredicatePartitionMethod}
  pred(i, j) = method.p₁(i, j) * method.p₂(i, j)
  partition(object, PredicatePartition(pred))
end

# spatial predicate partition method
function partition(object, method::ProductPartition{P1,P2}) where {P1<:SPredicatePartitionMethod,
                                                                   P2<:SPredicatePartitionMethod}
  pred(x, y) = method.p₁(x, y) * method.p₂(x, y)
  partition(object, SpatialPredicatePartition(pred))
end

*(p₁::PartitionMethod, p₂::PartitionMethod) = ProductPartition(p₁, p₂)


text:
using Images, FileIO, ImageIO
using QuartzImageIO, ImageMagick

filename = "128sq_p2_low"
searchdir(path, key) = filter(x -> occursin(key, x), readdir(path))
pngfiles = searchdir(pwd(), filename)

# to make gif
p = load(pwd() * "/" * pngfiles[1])
tmp = [load(pwd() * "/" * pngfiles[i]) for i in eachindex(pngfiles)]
size(tmp)
tmp2 = reduce(hcat, tmp)
tmp3 = reshape(tmp2, (size(tmp[1])...,length(tmp)))
save(filename * "_exasim.gif", tmp3)


text:
using Test, SpatialGraphs, LightGraphs, GeoData, ArchGDAL

url_base = "https://raw.githubusercontent.com/Circuitscape/datasets/main/"
# Download the NLCD tile used to create the resistance surface and load it
download(string(url_base, "data/nlcd_2016_frederick_md.tif"),
         "nlcd_2016_frederick_md.tif")

garray = GeoArray(GDALarray("nlcd_2016_frederick_md.tif", missingval = -9))

@testset "graph construction" begin
    no_data_val = -9999
    cost = [1 1 1;
            1 3 1;
            1 NaN no_data_val]

    nodemap = construct_nodemap(cost, no_data_val = no_data_val)

    # NoData entries in cost are 0 in nodemap
    @test (nodemap .== 0) == ((cost .== no_data_val) .| isnan.(cost))

    # Values in nodemap are as expected
    @test sort(collect(nodemap[nodemap .!= 0])) == collect(1:sum((cost .!= no_data_val) .& (!).(isnan.(cost))))

    ## Check that values in the resistance graph are as expected
    graph = construct_weighted_graph(cost, nodemap, no_data_val = no_data_val)
    the_edges = collect(edges(graph))

    # Test that the edges are correct and have proper weights
    for i in 1:length(the_edges)
        source_i = src(the_edges[i])
        dest_i = dst(the_edges[i])
        weight_i = the_edges[i].weight

        source_coords = findall(nodemap .== source_i)[1]
        dest_coords = findall(nodemap .== dest_i)[1]

        # Test that source row is within 1 step of dest row
        row_diff = abs(source_coords[1] - dest_coords[1])
        @test row_diff <= 1

        # Test that source column is within 1 step of dest row
        col_diff = abs(source_coords[2] - dest_coords[2])
        @test col_diff <= 1

        # Test that the weight is what it should be (assumes connect_using_avg_cost = true in graph construction)
        if (row_diff == 1 && col_diff == 1) # get diagonal average
            @test weight_i == SpatialGraphs.res_diagonal_avg(cost[source_coords], cost[dest_coords])
        else
            @test weight_i == SpatialGraphs.res_cardinal_avg(cost[source_coords], cost[dest_coords])
        end
    end

    ## Test that construct_weighted_graphs works when a vertex (ID = 1) takes up
    ## multiple pixels
    nodemap = [1 0 4;
               1 2 5;
               1 3 1]
    cost = [1 -1 -1;
            1  1  2;
            3  2  1]

    graph = construct_weighted_graph(cost, nodemap, no_data_val = -1)
    # there are multiple connections from nodes 1 to 2, want to make sure we're
    # using the minimum cost, which in the above case would be 1
    @test graph.weights[1, 2] == 1

    ## GeoArray compatibility
    old_data = deepcopy(garray.data[:, :, 1])
    nodemap = construct_nodemap(garray)

    @test size(nodemap) == size(garray[Band(Between(1, 1))])
    g = construct_weighted_graph(garray, nodemap,
                        cost_layer_is_conductance = true,
                        connect_four_neighbors_only = true,
                        connect_using_avg_cost = false)

    # Make sure none of the functions changed input
    @test old_data == float.(garray.data[:, :, 1])
end

@testset "internals" begin
    # get_lat_lon_dims, test that they are returned in the correct order
    ll_dims = SpatialGraphs.get_lat_lon_dims(garray)
    @test ll_dims[1] == dims(garray)[1]
    @test ll_dims[2] == dims(garray)[2]
end

# LightGraphs and SimpleWeightedGraphs presumably already test for correctness of
# methods, so above tests to confirm graphs are constructed properly are most
# important
@testset "paths" begin
    no_data_val = -9999
    weight = [1 1 1;
              1 3 1;
              1 1 1]

    nodemap = construct_nodemap(weight, no_data_val = no_data_val)
    g = construct_weighted_graph(weight, nodemap, no_data_val = no_data_val)

    cdist_array = cost_distance(g, nodemap, 5)
    @test size(cdist_array) == size(nodemap)
    @test all(cdist_array[nodemap[nodemap .∈ [collect(2:2:8)]]] .== 2.0)
    @test all(cdist_array[nodemap[nodemap .∈ [[1, 3, 7]]]] .==
        SpatialGraphs.res_diagonal_avg(3, 1))

    lcps, nodemap = random_lcps(weight, fill(1., size(weight)), 50, parallel = false)

    path = least_cost_path(g, 2, 8)
    path_coords = path_to_cartesian_coords(path, nodemap)
    path_array = path_to_array(path, nodemap)

    # With custom nodemap
    nodemap = [1 0 4;
               1 2 5;
               1 3 6]
    cost = [1 -1 -1;
            1  1  2;
            3  2  1]
    graph = construct_weighted_graph(cost, nodemap, no_data_val = -1)
    path = least_cost_path(graph, 1, 5)
    @test_throws BoundsError path_coords = path_to_cartesian_coords(path, nodemap)

    # GeoArray stuff
    nodemap = construct_nodemap(garray)
    g = construct_weighted_graph(garray, nodemap,
                        cost_layer_is_conductance = true,
                        connect_four_neighbors_only = false,
                        connect_using_avg_cost = false)
    path = least_cost_path(g, 2, 8)
    path_geoarray = path_to_geoarray(path, nodemap)
    path_linestring = path_to_linestring(path, nodemap)

    cd_geoarray = cost_distance(g, nodemap, [2, 3])

    # random lcps, and convert to geometries
    lcps, nodemap = random_lcps(garray, garray, 10)
    linestrings = paths_to_linestrings(lcps, nodemap)
end

rm("nlcd_2016_frederick_md.tif")

max_stars_repo_name: Circuitscape/SpatialGraphs-OLD

max_issues_repo_name: Circuitscape/SpatialGraphs-OLD

max_forks_repo_name: Circuitscape/SpatialGraphs-OLD


text:
# Hinge loss function

export HingeLoss

"""
**Hinge loss**

    HingeLoss(y, μ=1)

Returns the function
```math
f(x) = μ⋅∑_i \\max\\{0, 1 - y_i ⋅ x_i\\},
```
where `y` is an array and `μ` is a positive parameter.
"""
HingeLoss(y::T) where {T <: AbstractArray} =
    PrecomposeDiagonal(SumPositive(), -y, 1)

HingeLoss(y::T, mu::R) where {T <: AbstractArray, R <: Real} =
    Postcompose(PrecomposeDiagonal(SumPositive(), -y, 1), mu)


text:
using Test
using SparseIR

using Random, LinearAlgebra

include("__conftest.jl")

@testset "sampling.jl" begin
    @testset "alias" begin
        beta = 1
        wmax = 10
        basis = FiniteTempBasis(fermion, beta, wmax; sve_result=sve_logistic[beta * wmax])
        @test MatsubaraSampling(basis) isa MatsubaraSampling64
        @test TauSampling(basis) isa TauSampling64
    end

    @testset "decomp" begin
        Random.seed!(420)
        A = randn(49, 39)

        Ad = svd(A)
        norm_A = first(Ad.S) / last(Ad.S)
        @test all(isapprox.(A, Matrix(Ad), atol=1e-15 * norm_A, rtol=0))

        x = randn(39)
        @test A * x ≈ Ad.U * Diagonal(Ad.S) * Ad.Vt * x atol = 1e-14 * norm_A rtol = 0

        x = randn(39, 3)
        @test A * x ≈ Ad.U * Diagonal(Ad.S) * Ad.Vt * x atol = 2e-14 * norm_A rtol = 0

        y = randn(49)
        @test A \ y ≈ Ad \ y atol = 1e-14 * norm_A rtol = 0

        y = randn(49, 2)
        @test A \ y ≈ Ad \ y atol = 1e-14 * norm_A rtol = 0
    end

    @testset "fit from tau with stat = $stat, Λ = $Λ" for stat in (boson, fermion),
        Λ in (10, 42)
        basis = DimensionlessBasis(stat, Λ; sve_result=sve_logistic[Λ])
        smpl = TauSampling(basis)
        @test issorted(smpl.sampling_points)
        Random.seed!(5318008)

        shape = (2, 3, 4)
        rhol = randn(ComplexF64, (length(basis), shape...))
        originalgl = -basis.s .* rhol
        for dim in 1:ndims(rhol)
            gl = SparseIR.movedim(originalgl, 1 => dim)
            gtau = evaluate(smpl, gl; dim)
            @test size(gtau) == (
                size(gl)[1:(dim - 1)]...,
                length(smpl.sampling_points),
                size(gl)[(dim + 1):end]...,
            )

            gl_from_tau = fit(smpl, gtau; dim)
            @test gl_from_tau ≈ gl

            gl_from_tau2 = similar(gl_from_tau)
            fit!(gl_from_tau2, smpl, gtau; dim)
            @test gl_from_tau2 ≈ gl
        end
    end

    @testset "τ noise with stat = $stat, Λ = $Λ" for stat in (boson, fermion), Λ in (10, 42)
        basis = DimensionlessBasis(stat, Λ; sve_result=sve_logistic[Λ])
        smpl = TauSampling(basis)
        @test issorted(smpl.sampling_points)
        Random.seed!(5318008)

        ρℓ = basis.v([-0.999, -0.01, 0.5]) * [0.8, -0.2, 0.5]
        Gℓ = basis.s .* ρℓ
        Gℓ_magn = norm(Gℓ)
        @inferred evaluate(smpl, Gℓ)
        @inferred evaluate(smpl, Gℓ, dim=1)
        Gτ = evaluate(smpl, Gℓ)

        Gτ_inplace = similar(Gτ)
        evaluate!(Gτ_inplace, smpl, Gℓ)
        @test Gτ == Gτ_inplace

        noise = 1e-5
        Gτ_n = Gτ + noise * norm(Gτ) * randn(size(Gτ)...)
        @inferred fit(smpl, Gτ_n)
        @inferred fit(smpl, Gτ_n, dim=1)
        Gℓ_n = fit(smpl, Gτ_n)

        Gℓ_n_inplace = similar(Gℓ_n)
        fit!(Gℓ_n_inplace, smpl, Gτ_n)
        @test Gℓ_n == Gℓ_n_inplace

        @test isapprox(Gℓ, Gℓ_n, atol=12 * noise * Gℓ_magn, rtol=0)
    end

    @testset "iω noise with stat = $stat, Λ = $Λ" for stat in (boson, fermion),
        Λ in (10, 42)

        basis = DimensionlessBasis(stat, Λ; sve_result=sve_logistic[Λ])
        smpl = MatsubaraSampling(basis)
        @test issorted(smpl.sampling_points)
        Random.seed!(1312 + 161)

        ρℓ = basis.v([-0.999, -0.01, 0.5]) * [0.8, -0.2, 0.5]
        Gℓ = basis.s .* ρℓ
        Gℓ_magn = norm(Gℓ)
        @inferred evaluate(smpl, Gℓ)
        @inferred evaluate(smpl, Gℓ, dim=1)
        Giw = evaluate(smpl, Gℓ)

        Giw_inplace = similar(Giw)
        evaluate!(Giw_inplace, smpl, Gℓ)
        @test Giw == Giw_inplace

        noise = 1e-5
        Giwn_n = Giw + noise * norm(Giw) * randn(size(Giw)...)
        @inferred fit(smpl, Giwn_n)
        @inferred fit(smpl, Giwn_n, dim=1)
        Gℓ_n = fit(smpl, Giwn_n)

        Gℓ_n_inplace = similar(Gℓ_n)
        fit!(Gℓ_n_inplace, smpl, Giwn_n)
        @test Gℓ_n == Gℓ_n_inplace

        @test isapprox(Gℓ, Gℓ_n, atol=40 * noise * Gℓ_magn, rtol=0)
    end
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

function plot_histogram(h::AbstractHistogram, swap::Bool)
    if swap
      return h.weights, h.edges[1][1:end-1]
    else
        h.edges[1][1:end-1], h.weights
    end
end

# ToDo: This is type piracy, find a cleaner solution!
@recipe function f(hist::Histogram,
                    param::Integer;
                    intervals = standard_confidence_vals,
                    bins=200,
                    normalize = true,
                    colors = standard_colors,
                    intervallabels = [])

    seriestype = get(plotattributes, :seriestype, :stephist)

    orientation = get(plotattributes, :orientation, :vertical)
    (orientation != :vertical) ? swap=true : swap = false
    plotattributes[:orientation] = :vertical # without: auto-scaling of axes not correct

    dims = collect(1:ndims(hist.weights))
    if length(dims) > 1
        weights = vec(sum(hist.weights, dims=setdiff(dims, param)))
        hist = Histogram(hist.edges[param], weights, :left)
    end

    normalize ? hist=StatsBase.normalize(hist) : nothing
    weights = hist.weights

    if seriestype == :histogram

         @series begin
            seriestype := :steppost
            label --> ""
            fillrange --> 0
            fillcolor --> :dodgerblue
            linewidth --> 0
            plot_histogram(hist, swap)
        end

    elseif seriestype == :stephist || seriestype == :steps #TODO: :steps not working -> "seriestype steppost not supported"
        @series begin
            seriestype := :steppost
            label --> ""
            linecolor --> :dodgerblue
            plot_histogram(hist, swap)
        end

    elseif seriestype == :smallest_intervals
        hists, orig_hist, realintervals = split_smallest(hist, intervals)

        colors = colors[sortperm(intervals, rev=true)]

        for i in 1:length(realintervals)
            @series begin
                seriestype := :steppost
                fillcolor --> colors[i]
                linewidth --> 0
                fillrange --> 0

                if length(intervallabels) > 0
                    label := intervallabels[i]
                else
                    label := "smallest $(@sprintf("%.2f", realintervals[i]*100))% interval(s)"
                end
                plot_histogram(hists[i], swap)
            end
        end

        @series begin
            seriestype := :steppost
            linecolor --> :black
            linewidth --> 0.7
            label --> ""
            plot_histogram(orig_hist, swap)
        end

    elseif seriestype == :central_intervals
        hists, orig_hist, realintervals = split_central(hist, intervals)

        colors = colors[sortperm(intervals, rev=true)]

        for i in 1:length(realintervals)
            @series begin
                seriestype := :steppost
                fillcolor --> colors[i]
                linewidth --> 0
                fillrange --> 0
                label --> "central $(@sprintf("%.2f", realintervals[i]*100))% interval(s)"
                plot_histogram(hists[i], swap)
            end
        end

         @series begin
            seriestype := :steppost
            linecolor --> :black
            linewidth --> 0.7
            label --> ""
            plot_histogram(orig_hist, swap)
        end

    else
        error("seriestype $seriestype not supported")
    end

end


text:
global _current_param_set

ParamSet = Dict{String,String}

"""
    activate_param_set(d::Dict{String,String})

Activate the parameter set from a dictionary `d`.
"""
function activate_param_set(d::Dict{String,String})
    global _current_param_set = d
end

"""
    activate_param_set(d::String)

Activate the parameter set from a json file, whose path is given by `d`.
"""
function activate_param_set(d::String)
    global _current_param_set = JSON.parsefile(d)
end

"""
    activate_param_set(name=>formatter...)

Activate the parameter set with given `name` and `formatter` pairs.
"""
function activate_param_set(args...)
    global _current_param_set = ParamSet(args...)
end

"""
    get_param_set()

Return current active parameter set as a dictionary.
"""
function get_param_set()
    _current_param_set
end

"""
    load_config_from_json(name) -> DataEntry, Tasks

Load tasks from configuration file.
"""
function load_config_from_json(name)
    t = JSON.parsefile(name)
    if haskey(t, "default_values") && t["default_values"] != nothing
        default_vals = parse_string_to_num(t["default_values"])
    else
        default_vals = Dict()
    end
    f(x) = typeof(x[1]) == String ? [x] : x
    res = Dict{String,DataEntry}()
    for (k, v) in t["tasks"]
        if haskey(v, "folder") && (v["folder"] != nothing) && (!isempty(v["folder"]))
            folder_path = f(v["folder"])
        else
            folder_path = ()
        end
        if haskey(v, "group") && (v["group"] != nothing) && (!isempty(v["group"]))
            group_path = f(v["group"])
        else
            group_path = ()
        end
        root_path = haskey(v, "root") ? joinpath(t["root"], v["root"]) : t["root"]
        res[k] = DataEntry(root_path, folder_path, group_path)
    end
    res, default_vals
end


text:
using NFEdocs
using Test

@testset "NFEdocs.jl" begin
    # This is just for doc generation.
end


text:
##############################
# Global variables/constants #
##############################
const ADBACKEND = Ref(:forwarddiff)
setadbackend(backend_sym::Symbol) = setadbackend(Val(backend_sym))
function setadbackend(backend::Val)
    _setadbackend(backend)
    AdvancedVI.setadbackend(backend)
    Bijectors.setadbackend(backend)
end

function _setadbackend(::Val{:forwarddiff})
    CHUNKSIZE[] == 0 && setchunksize(40)
    ADBACKEND[] = :forwarddiff
end
function _setadbackend(::Val{:tracker})
    ADBACKEND[] = :tracker
end
function _setadbackend(::Val{:zygote})
    ADBACKEND[] = :zygote
end

const ADSAFE = Ref(false)
function setadsafe(switch::Bool)
    @info("[Turing]: global ADSAFE is set as $switch")
    ADSAFE[] = switch
end

const CHUNKSIZE = Ref(40) # default chunksize used by AD

function setchunksize(chunk_size::Int)
    if ~(CHUNKSIZE[] == chunk_size)
        @info("[Turing]: AD chunk size is set as $chunk_size")
        CHUNKSIZE[] = chunk_size
    end
end

abstract type ADBackend end
struct ForwardDiffAD{chunk} <: ADBackend end
getchunksize(::Type{<:ForwardDiffAD{chunk}}) where chunk = chunk
getchunksize(::Type{<:Sampler{Talg}}) where Talg = getchunksize(Talg)
getchunksize(::Type{SampleFromPrior}) = CHUNKSIZE[]

struct TrackerAD <: ADBackend end
struct ZygoteAD <: ADBackend end

ADBackend() = ADBackend(ADBACKEND[])
ADBackend(T::Symbol) = ADBackend(Val(T))

ADBackend(::Val{:forwarddiff}) = ForwardDiffAD{CHUNKSIZE[]}
ADBackend(::Val{:tracker}) = TrackerAD
ADBackend(::Val{:zygote}) = ZygoteAD
ADBackend(::Val) = error("The requested AD backend is not available. Make sure to load all required packages.")

"""
    getADbackend(alg)

Find the autodifferentiation backend of the algorithm `alg`.
"""
getADbackend(spl::Sampler) = getADbackend(spl.alg)
getADbackend(spl::SampleFromPrior) = ADBackend()()

"""
    gradient_logp(
        θ::AbstractVector{<:Real},
        vi::VarInfo,
        model::Model,
        sampler::AbstractSampler,
        ctx::DynamicPPL.AbstractContext = DynamicPPL.DefaultContext()
    )

Computes the value of the log joint of `θ` and its gradient for the model
specified by `(vi, sampler, model)` using whichever automatic differentation
tool is currently active.
"""
function gradient_logp(
    θ::AbstractVector{<:Real},
    vi::VarInfo,
    model::Model,
    sampler::AbstractSampler,
    ctx::DynamicPPL.AbstractContext = DynamicPPL.DefaultContext()
)
    return gradient_logp(getADbackend(sampler), θ, vi, model, sampler, ctx)
end

"""
gradient_logp(
    backend::ADBackend,
    θ::AbstractVector{<:Real},
    vi::VarInfo,
    model::Model,
    sampler::AbstractSampler = SampleFromPrior(),
    ctx::DynamicPPL.AbstractContext = DynamicPPL.DefaultContext()
)

Compute the value of the log joint of `θ` and its gradient for the model
specified by `(vi, sampler, model)` using `backend` for AD, e.g. `ForwardDiffAD{N}()` uses `ForwardDiff.jl` with chunk size `N`, `TrackerAD()` uses `Tracker.jl` and `ZygoteAD()` uses `Zygote.jl`.
"""
function gradient_logp(
    ::ForwardDiffAD,
    θ::AbstractVector{<:Real},
    vi::VarInfo,
    model::Model,
    sampler::AbstractSampler=SampleFromPrior(),
    ctx::DynamicPPL.AbstractContext = DynamicPPL.DefaultContext()
)
    # Define function to compute log joint.
    logp_old = getlogp(vi)
    function f(θ)
        new_vi = VarInfo(vi, sampler, θ)
        model(new_vi, sampler, ctx)
        logp = getlogp(new_vi)
        setlogp!(vi, ForwardDiff.value(logp))
        return logp
    end

    chunk_size = getchunksize(typeof(sampler))
    # Set chunk size and do ForwardMode.
    chunk = ForwardDiff.Chunk(min(length(θ), chunk_size))
    config = ForwardDiff.GradientConfig(f, θ, chunk)
    ∂l∂θ = ForwardDiff.gradient!(similar(θ), f, θ, config)
    l = getlogp(vi)
    setlogp!(vi, logp_old)

    return l, ∂l∂θ
end
function gradient_logp(
    ::TrackerAD,
    θ::AbstractVector{<:Real},
    vi::VarInfo,
    model::Model,
    sampler::AbstractSampler = SampleFromPrior(),
    ctx::DynamicPPL.AbstractContext = DynamicPPL.DefaultContext()
)
    T = typeof(getlogp(vi))

    # Specify objective function.
    function f(θ)
        new_vi = VarInfo(vi, sampler, θ)
        model(new_vi, sampler, ctx)
        return getlogp(new_vi)
    end

    # Compute forward and reverse passes.
    l_tracked, ȳ = Tracker.forward(f, θ)
    # Remove tracking info from variables in model (because mutable state).
    l::T, ∂l∂θ::typeof(θ) = Tracker.data(l_tracked), Tracker.data(ȳ(1)[1])

    return l, ∂l∂θ
end

function gradient_logp(
    backend::ZygoteAD,
    θ::AbstractVector{<:Real},
    vi::VarInfo,
    model::Model,
    sampler::AbstractSampler = SampleFromPrior(),
    context::DynamicPPL.AbstractContext = DynamicPPL.DefaultContext()
)
    T = typeof(getlogp(vi))

    # Specify objective function.
    function f(θ)
        new_vi = VarInfo(vi, sampler, θ)
        model(new_vi, sampler, context)
        return getlogp(new_vi)
    end

    # Compute forward and reverse passes.
    l::T, ȳ = ZygoteRules.pullback(f, θ)
    ∂l∂θ::typeof(θ) = ȳ(1)[1]

    return l, ∂l∂θ
end

function verifygrad(grad::AbstractVector{<:Real})
    if any(isnan, grad) || any(isinf, grad)
        @warn("Numerical error in gradients. Rejecting current proposal...")
        @warn("grad = $(grad)")
        return false
    else
        return true
    end
end

# These still seem necessary
for F in (:link, :invlink)
    @eval begin
        $F(dist::PDMatDistribution, x::Tracker.TrackedArray) = Tracker.track($F, dist, x)
        Tracker.@grad function $F(dist::PDMatDistribution, x::Tracker.TrackedArray)
            x_data = Tracker.data(x)
            T = eltype(x_data)
            y = $F(dist, x_data)
            return  y, Δ -> begin
                out = reshape((ForwardDiff.jacobian(x -> $F(dist, x), x_data)::Matrix{T})' * vec(Δ), size(Δ))
                return (nothing, out)
            end
        end
    end
end