text:
module SpinPackage

using Printf: @printf,@sprintf
using ..Spatials: PID
using ..DegreesOfFreedom: IID,Internal,Index,FilteredAttributes,OID,Operator,LaTeX,latexformat
using ..Terms: wildcard,constant,Subscript,Subscripts,subscriptexpr,Coupling,Couplings,couplingcenters,Term,TermCouplings,TermAmplitude,TermModulate
using ...Interfaces: rank,kind
using ...Prerequisites: Float,decimaltostr,delta
using ...Mathematics.VectorSpaces: VectorSpace,IsMultiIndexable,MultiIndexOrderStyle
using ...Mathematics.AlgebraOverFields: SimpleID,ID

import ..DegreesOfFreedom: script,otype,isHermitian
import ..Terms: couplingcenter,statistics,abbr
import ...Interfaces: dims,inds,expand,matrix,permute

export SID,Spin,SIndex
export usualspinindextotuple
export SOperator,soptdefaultlatex
export SCID,SpinCoupling
export Heisenberg,Ising,Gamma,DM,Sˣ,Sʸ,Sᶻ
export @sc_str,@heisenberg_str,@ising_str,@gamma_str,@dm_str,@sˣ_str,@sʸ_str,@sᶻ_str
export SpinTerm

const sidtagmap=Dict(1=>'x',2=>'y',3=>'z',4=>'+',5=>'-')
const sidseqmap=Dict(v=>k for (k,v) in sidtagmap)
const sidajointmap=Dict('x'=>'x','y'=>'y','z'=>'z','+'=>'-','-'=>'+')

"""
    SID <: IID

The spin id.
"""
struct SID <: IID
    orbital::Int
    spin::Float
    tag::Char
    function SID(orbital::Int,spin::Real,tag::Char)
        @assert tag in ('x','y','z','+','-') "SID error: not supported tag($tag)."
        new(orbital,spin,tag)
    end
end
Base.fieldnames(::Type{<:SID})=(:orbital,:spin,:tag)

"""
    SID(;orbital::Int=1,spin::Real=0.5,tag::Char='z')

Create a spin id.
"""
SID(;orbital::Int=1,spin::Real=0.5,tag::Char='z')=SID(orbital,spin,tag)

"""
    adjoint(sid::SID) -> SID

Get the adjoint of a spin id.
"""
Base.adjoint(sid::SID)=SID(sid.orbital,sid.spin,sidajointmap[sid.tag])

"""
    matrix(sid::SID,dtype::Type{<:Number}=Complex{Float}) -> Matrix{dtype}

Get the matrix representation of a sid.
"""
function matrix(sid::SID,dtype::Type{<:Number}=Complex{Float})
    N=Int(2*sid.spin+1)
    result=zeros(dtype,(N,N))
    for i=1:N,j=1:N
        row,col=N+1-i,N+1-j
        m,n=sid.spin+1-i,sid.spin+1-j
        result[row,col]=sid.tag=='x' ? (delta(i+1,j)+delta(i,j+1))*sqrt(sid.spin*(sid.spin+1)-m*n)/2 :
                        sid.tag=='y' ? (delta(i+1,j)-delta(i,j+1))*sqrt(sid.spin*(sid.spin+1)-m*n)/2im :
                        sid.tag=='z' ? delta(i,j)*m :
                        sid.tag=='+' ? delta(i+1,j)*sqrt(sid.spin*(sid.spin+1)-m*n) :
                                       delta(i,j+1)*sqrt(sid.spin*(sid.spin+1)-m*n)
    end
    return result
end

"""
    Spin <: Internal{SID}

The spin interanl degrees of freedom.
"""
struct Spin <: Internal{SID}
    atom::Int
    norbital::Int
    spin::Float
end
IsMultiIndexable(::Type{Spin})=IsMultiIndexable(true)
MultiIndexOrderStyle(::Type{Spin})=MultiIndexOrderStyle('C')
dims(sp::Spin)=(sp.norbital,length(sidtagmap))
inds(sid::SID,::Spin)=(sid.orbital,sidseqmap[sid.tag])
SID(inds::NTuple{2,Int},sp::Spin)=SID(inds[1],sp.spin,sidtagmap[inds[2]])

"""
    Spin(;atom::Int=1,norbital::Int=1,spin::Real=0.5)

Construct a spin degrees of freedom.
"""
Spin(;atom::Int=1,norbital::Int=1,spin::Real=0.5)=Spin(atom,norbital,spin)

"""
    SIndex{S} <: Index{PID{S},SID}

The spin index.
"""
struct SIndex{S} <: Index{PID{S},SID}
    scope::S
    site::Int
    orbital::Int
    spin::Float
    tag::Char
end
Base.fieldnames(::Type{<:SIndex})=(:scope,:site,:orbital,:spin,:tag)

"""
    union(::Type{P},::Type{SID}) where {P<:PID}

Get the union type of `PID` and `SID`.
"""
Base.union(::Type{P},::Type{SID}) where {P<:PID}=SIndex{fieldtype(P,:scope)}

"""
    script(oid::OID{<:SIndex},::Val{:site}) -> Int
    script(oid::OID{<:SIndex},::Val{:orbital}) -> Int
    script(oid::OID{<:SIndex},::Val{:spin}) -> Float
    script(oid::OID{<:SIndex},::Val{:tag}) -> Char

Get the required script of a spin oid.
"""
script(oid::OID{<:SIndex},::Val{:site})=oid.index.site
script(oid::OID{<:SIndex},::Val{:orbital})=oid.index.orbital
script(oid::OID{<:SIndex},::Val{:spin})=oid.index.spin
script(oid::OID{<:SIndex},::Val{:tag})=oid.index.tag

"""
    usualspinindextotuple

Indicate that the filtered attributes are `(:scope,:site,:orbital)` when converting a spin index to tuple.
"""
const usualspinindextotuple=FilteredAttributes(:scope,:site,:orbital)

"""
    SOperator(value::Number,id::ID{OID}=ID())

Spin operator.
"""
struct SOperator{V<:Number,I<:ID{OID}} <: Operator{V,I}
    value::V
    id::I
    SOperator(value::Number)=new{typeof(vlalue),Tuple{}}(value,())
    SOperator(value::Number,id::ID{OID})=new{typeof(value),typeof(id)}(value,id)
end

"""
    soptdefaultlatex

The default LaTeX format of the oids of a spin operator.
"""
const soptdefaultlatex=LaTeX{(:tag,),(:site,:orbital)}('S')
latexformat(SOperator,soptdefaultlatex)

"""
    statistics(opt::SOperator) -> Char
    statistics(::Type{<:SOperator}) -> Char

Get the statistics of SOperator.
"""
statistics(opt::SOperator)=opt|>typeof|>statistics
statistics(::Type{<:SOperator})='B'

"""
    permute(::Type{<:SOperator},id1::OID{<:SIndex},id2::OID{<:SIndex},table) -> Tuple{Vararg{SOperator}}

Permute two fermionic oid and get the result.
"""
function permute(::Type{<:SOperator},id1::OID{<:SIndex},id2::OID{<:SIndex},table)
    @assert id1.index≠id2.index || id1.rcoord≠id2.rcoord || id1.icoord≠id2.icoord "permute error: permuted ids should not be equal to each other."
    if usualspinindextotuple(id1.index)==usualspinindextotuple(id2.index) && id1.rcoord==id2.rcoord && id1.icoord==id2.icoord
        @assert id1.index.spin==id2.index.spin "permute error: noncommutable ids should have the same spin attribute."
        if id1.index.tag=='x'
            id2.index.tag=='y' && return (SOperator(+1im,ID(permutesoid(id1,'z',table))),SOperator(1,ID(id2,id1)))
            id2.index.tag=='z' && return (SOperator(-1im,ID(permutesoid(id1,'y',table))),SOperator(1,ID(id2,id1)))
            id2.index.tag=='+' && return (SOperator(-1,ID(permutesoid(id1,'z',table))),SOperator(1,ID(id2,id1)))
            id2.index.tag=='-' && return (SOperator(+1,ID(permutesoid(id1,'z',table))),SOperator(1,ID(id2,id1)))
        elseif id1.index.tag=='y'
            id2.index.tag=='x' && return (SOperator(-1im,ID(permutesoid(id1,'z',table))),SOperator(1,ID(id2,id1)))
            id2.index.tag=='z' && return (SOperator(+1im,ID(permutesoid(id1,'x',table))),SOperator(1,ID(id2,id1)))
            id2.index.tag=='+' && return (SOperator(-1im,ID(permutesoid(id1,'z',table))),SOperator(1,ID(id2,id1)))
            id2.index.tag=='-' && return (SOperator(-1im,ID(permutesoid(id1,'z',table))),SOperator(1,ID(id2,id1)))
        elseif id1.index.tag=='z'
            id2.index.tag=='x' && return (SOperator(+1im,ID(permutesoid(id1,'y',table))),SOperator(1,ID(id2,id1)))
            id2.index.tag=='y' && return (SOperator(-1im,ID(permutesoid(id1,'x',table))),SOperator(1,ID(id2,id1)))
            id2.index.tag=='+' && return (SOperator(+1,ID(id2)),SOperator(1,ID(id2,id1)))
            id2.index.tag=='-' && return (SOperator(-1,ID(id2)),SOperator(1,ID(id2,id1)))
        elseif id1.index.tag=='+'
            id2.index.tag=='x' && return (SOperator(+1,ID(permutesoid(id1,'z',table))),SOperator(1,ID(id2,id1)))
            id2.index.tag=='y' && return (SOperator(+1im,ID(permutesoid(id1,'z',table))),SOperator(1,ID(id2,id1)))
            id2.index.tag=='z' && return (SOperator(-1,ID(id1)),SOperator(1,ID(id2,id1)))
            id2.index.tag=='-' && return (SOperator(+2,ID(permutesoid(id1,'z',table))),SOperator(1,ID(id2,id1)))
        elseif id1.index.tag=='-'
            id2.index.tag=='x' && return (SOperator(-1,ID(permutesoid(id1,'z',table))),SOperator(1,ID(id2,id1)))
            id2.index.tag=='y' && return (SOperator(1im,ID(permutesoid(id1,'z',table))),SOperator(1,ID(id2,id1)))
            id2.index.tag=='z' && return (SOperator(+1,ID(id1)),SOperator(1,ID(id2,id1)))
            id2.index.tag=='+' && return (SOperator(-2,ID(permutesoid(id1,'z',table))),SOperator(1,ID(id2,id1)))
        end
    else
        return (SOperator(1,ID(id2,id1)),)
    end
end
function permutesoid(id::OID{<:SIndex},tag::Char,table)
    index=replace(id.index,tag=tag)
    fieldtype(id|>typeof,:seq)<:Nothing && return replace(id,index=index)
    seq=get(table,index,nothing)
    seq===nothing && error("permutesoid error: incomplete index-sequence table with $index not found.")
    return replace(id,index=index,seq=seq)
end

"""
    SCID(;center=wildcard,atom=wildcard,orbital=wildcard,tag='z',subscript=wildcard)=SCID(center,atom,orbital,tag,subscript)

The id of a spin coupling.
"""
struct SCID{C,A,O,S} <: SimpleID
    center::C
    atom::A
    orbital::O
    tag::Char
    subscript::S
    function SCID(center,atom,orbital,tag::Char,subscript)
        @assert tag in ('x','y','z','+','-') "SCID error: not supported tag($tag)."
        new{typeof(center),typeof(atom),typeof(orbital),typeof(subscript)}(center,atom,orbital,tag,subscript)
    end
end
SCID(;center=wildcard,atom=wildcard,orbital=wildcard,tag='z',subscript=wildcard)=SCID(center,atom,orbital,tag,subscript)
Base.fieldnames(::Type{<:SCID})=(:center,:atom,:orbital,:tag,:subscript)

"""
    SpinCoupling(value::Number,id::ID{SCID},subscripts::Subscripts)
    SpinCoupling{N}(    value::Number=1;
                        tags::NTuple{N,Char},
                        centers::Union{NTuple{N,Int},Nothing}=nothing,
                        atoms::Union{NTuple{N,Int},Nothing}=nothing,
                        orbitals::Union{NTuple{N,Int},Subscript,Nothing}=nothing
                        ) where N

Spin coupling.
"""
struct SpinCoupling{V,I<:ID{SCID},S<:Subscripts} <: Coupling{V,I}
    value::V
    id::I
    subscripts::S
    function SpinCoupling(value::Number,id::ID{SCID},subscripts::Subscripts)
        new{typeof(value),typeof(id),typeof(subscripts)}(value,id,subscripts)
    end
end
function SpinCoupling{N}(   value::Number=1;
                            tags::NTuple{N,Char},
                            centers::Union{NTuple{N,Int},Nothing}=nothing,
                            atoms::Union{NTuple{N,Int},Nothing}=nothing,
                            orbitals::Union{NTuple{N,Int},Subscript,Nothing}=nothing
                            ) where N
    centers===nothing && (centers=ntuple(i->wildcard,N))
    atoms===nothing && (atoms=ntuple(i->wildcard,N))
    subscript=isa(orbitals,NTuple{N,Int}) ? Subscript(orbitals) : isa(orbitals,Nothing) ? Subscript{N}() : orbitals
    orbitals=subscript.opattern
    subscripts=ntuple(i->subscript.identifier,N)
    return SpinCoupling(value,ID(SCID,centers,atoms,orbitals,tags,subscripts),Subscripts(subscript))
end

"""
    show(io::IO,sc::SpinCoupling)

Show a spin coupling.
"""
function Base.show(io::IO,sc::SpinCoupling)
    @printf io "SpinCoupling(value=%s" decimaltostr(sc.value)
    cache=[]
    for attrname in (:centers,:atoms,:orbitals,:tags)
        attrvalue=getproperty(sc.id,attrname)
        all(attrvalue.≠wildcard) && @printf io ",%s=(%s)" attrname join(attrvalue,",")
    end
    all((subscripts=sc.id.subscripts).==wildcard) || all((subscripts.==constant)) || @printf io ",subscripts=(%s)" join(subscripts,",")
    @printf io ")"
end

"""
    repr(sc::SpinCoupling) -> String

Get the repr representation of a spin coupling.
"""
function Base.repr(sc::SpinCoupling)
    cache=[]
    for (attrname,abbr) in zip((:atoms,:orbitals),("sl","ob"))
        attrvalue=getproperty(sc.id,attrname)
        all(attrvalue.≠wildcard) && push!(cache,@sprintf "%s(%s)" abbr join(attrvalue,','))
    end
    result=@sprintf "%s S%s" decimaltostr(sc.value) join(sc.id.tags,'S')
    length(cache)>0 && (result=@sprintf "%s %s" result join(cache,"⊗"))
    any((centers=sc.id.centers).≠wildcard) && (result=@sprintf "%s%s@(%s)" result (length(cache)>0 ? "" : " ") join(centers,','))
    subscripts=sc.id.subscripts
    all(subscripts.==wildcard) || all(subscripts.==constant) || (result=@sprintf "%s with %s" result join(subscripts," && "))
    return result
end

"""
    *(sc1::SpinCoupling,sc2::SpinCoupling) -> SpinCoupling

Get the multiplication between two spin couplings.
"""
Base.:*(sc1::SpinCoupling,sc2::SpinCoupling)=SpinCoupling(sc1.value*sc2.value,sc1.id*sc2.id,sc1.subscripts*sc2.subscripts)

"""
    expand(sc::SpinCoupling,pid::PID,spin::Spin,kind::Union{Val{K},Nothing}=nothing) where K -> Union{SCExpand,Tuple{}}
    expand(sc::SpinCoupling,pids::NTuple{N,PID},spins::NTuple{N,Spin},kind::Union{Val{K},Nothing}=nothing) where {N,K} -> Union{SCExpand,Tuple{}}

Expand a spin coupling with the given set of point ids and spin degrees of freedom.
"""
expand(sc::SpinCoupling,pid::PID,spin::Spin,kind::Union{Val{K},Nothing}=nothing) where K=expand(sc,(pid,),(spin,),kind)
function expand(sc::SpinCoupling,pids::NTuple{N,PID},spins::NTuple{N,Spin},kind::Union{Val{K},Nothing}=nothing) where {N,K}
    centers=couplingcenters(typeof(sc),sc.id.centers,Val(N))
    rpids=NTuple{rank(sc),eltype(pids)}(pids[centers[i]] for i=1:rank(sc))
    rspins=NTuple{rank(sc),eltype(spins)}(spins[centers[i]] for i=1:rank(sc))
    sps=NTuple{rank(sc),Float}(rspins[i].spin for i=1:rank(sc))
    for (i,atom) in enumerate(sc.id.atoms)
        isa(atom,Int) && (atom≠rspins[i].atom) && return ()
    end
    sbexpands=expand(sc.subscripts,NTuple{rank(sc),Int}(rspins[i].norbital for i=1:rank(sc)))|>collect
    return SCExpand(sc.value,rpids,sbexpands,sps,sc.id.tags)
end
couplingcenter(::Type{<:SpinCoupling},i::Int,n::Int,::Val{2})=i%2==1 ? 1 : 2
struct SCExpand{V,N,S} <: VectorSpace{Tuple{V,NTuple{N,SIndex{S}}}}
    value::V
    pids::NTuple{N,PID{S}}
    sbexpands::Vector{NTuple{N,Int}}
    spins::NTuple{N,Float}
    tags::NTuple{N,Char}
end
IsMultiIndexable(::Type{<:SCExpand})=IsMultiIndexable(true)
MultiIndexOrderStyle(::Type{<:SCExpand})=MultiIndexOrderStyle('C')
dims(sce::SCExpand)=(length(sce.sbexpands),)
@generated function Tuple(index::Tuple{Int},sce::SCExpand{V,N}) where {V,N}
    exprs=[:(SIndex(sce.pids[$i],SID(sce.sbexpands[index[1]][$i],sce.spins[$i],sce.tags[$i]))) for i=1:N]
    return Expr(:tuple,:(sce.value),Expr(:tuple,exprs...))
end

"""
    Heisenberg( mode::String="+-z";
                centers::Union{NTuple{2,Int},Nothing}=nothing,
                atoms::Union{NTuple{2,Int},Nothing}=nothing,
                orbitals::Union{NTuple{2,Int},Subscript,Nothing}=nothing
                ) -> Couplings{ID{SCID,2},SpinCoupling{Rational{Int}/Int,<:ID{SCID,2}}}

The Heisenberg couplings.
"""
function Heisenberg(mode::String="+-z";
                    centers::Union{NTuple{2,Int},Nothing}=nothing,
                    atoms::Union{NTuple{2,Int},Nothing}=nothing,
                    orbitals::Union{NTuple{2,Int},Subscript,Nothing}=nothing
                    )
    @assert mode=="+-z" || mode=="xyz" "Heisenberg error: not supported mode($mode)."
    if mode=="+-z"
        sc1=SpinCoupling{2}(1//2,centers=centers,atoms=atoms,orbitals=orbitals,tags=('+','-'))
        sc2=SpinCoupling{2}(1//2,centers=centers,atoms=atoms,orbitals=orbitals,tags=('-','+'))
        sc3=SpinCoupling{2}(1//1,centers=centers,atoms=atoms,orbitals=orbitals,tags=('z','z'))
    else
        sc1=SpinCoupling{2}(1,centers=centers,atoms=atoms,orbitals=orbitals,tags=('x','x'))
        sc2=SpinCoupling{2}(1,centers=centers,atoms=atoms,orbitals=orbitals,tags=('y','y'))
        sc3=SpinCoupling{2}(1,centers=centers,atoms=atoms,orbitals=orbitals,tags=('z','z'))
    end
    return Couplings(sc1,sc2,sc3)
end

"""
    Ising(  tag::Char;
            centers::Union{NTuple{2,Int},Nothing}=nothing,
            atoms::Union{NTuple{2,Int},Nothing}=nothing,
            orbitals::Union{NTuple{2,Int},Subscript,Nothing}=nothing
            ) -> Couplings{ID{SCID,2},SpinCoupling{Rational{Int}/Int,<:ID{SCID,2}}}

The Ising couplings.
"""
function Ising(tag::Char;centers::Union{NTuple{2,Int},Nothing}=nothing,atoms::Union{NTuple{2,Int},Nothing}=nothing,orbitals::Union{NTuple{2,Int},Subscript,Nothing}=nothing)
    @assert tag in ('x','y','z') "Ising error: not supported input tag($tag)."
    return Couplings(SpinCoupling{2}(1,centers=centers,atoms=atoms,orbitals=orbitals,tags=(tag,tag)))
end

"""
    Gamma(  tag::Char;
            centers::Union{NTuple{2,Int},Nothing}=nothing,
            atoms::Union{NTuple{2,Int},Nothing}=nothing,
            orbitals::Union{NTuple{2,Int},Subscript,Nothing}=nothing
            ) -> Couplings{ID{SCID,2},SpinCoupling{Rational{Int}/Int,<:ID{SCID,2}}}

The Gamma couplings.
"""
function Gamma(tag::Char;centers::Union{NTuple{2,Int},Nothing}=nothing,atoms::Union{NTuple{2,Int},Nothing}=nothing,orbitals::Union{NTuple{2,Int},Subscript,Nothing}=nothing)
    @assert tag in ('x','y','z') "Gamma error: not supported input tag($tag)."
    t1,t2=tag=='x' ? ('y','z') : tag=='y' ? ('z','x') : ('x','y')
    sc1=SpinCoupling{2}(1,centers=centers,atoms=atoms,orbitals=orbitals,tags=(t1,t2))
    sc2=SpinCoupling{2}(1,centers=centers,atoms=atoms,orbitals=orbitals,tags=(t2,t1))
    return Couplings(sc1,sc2)
end

"""
    DM( tag::Char;
        centers::Union{NTuple{2,Int},Nothing}=nothing,
        atoms::Union{NTuple{2,Int},Nothing}=nothing,
        orbitals::Union{NTuple{2,Int},Subscript,Nothing}=nothing
        ) -> Couplings{ID{SCID,2},SpinCoupling{Rational{Int}/Int,<:ID{SCID,2}}}
"""
function DM(tag::Char;centers::Union{NTuple{2,Int},Nothing}=nothing,atoms::Union{NTuple{2,Int},Nothing}=nothing,orbitals::Union{NTuple{2,Int},Subscript,Nothing}=nothing)
    @assert tag in ('x','y','z') "DM error: not supported input tag($tag)."
    t1,t2=tag=='x' ? ('y','z') : tag=='y' ? ('z','x') : ('x','y')
    sc1=SpinCoupling{2}(1,centers=centers,atoms=atoms,orbitals=orbitals,tags=(t1,t2))
    sc2=SpinCoupling{2}(-1,centers=centers,atoms=atoms,orbitals=orbitals,tags=(t2,t1))
    return Couplings(sc1,sc2)
end

scsinglewrapper(::Nothing)=nothing
scsinglewrapper(value::Int)=(value,)

"""
    Sˣ(;atom::Union{Int,Nothing}=nothing,orbital::Union{Int,Nothing}=nothing) -> Couplings{ID{SCID,1},SpinCoupling{Rational{Int}/Int,<:ID{SCID,1}}}

The single Sˣ coupling.
"""
function Sˣ(;atom::Union{Int,Nothing}=nothing,orbital::Union{Int,Nothing}=nothing)
    Couplings(SpinCoupling{1}(1,tags=('x',),atoms=scsinglewrapper(atom),orbitals=scsinglewrapper(orbital)))
end

"""
    Sʸ(;atom::Union{Int,Nothing}=nothing,orbital::Union{Int,Nothing}=nothing) -> Couplings{ID{SCID,1},SpinCoupling{Rational{Int}/Int,<:ID{SCID,1}}}

The single Sʸ coupling.
"""
function Sʸ(;atom::Union{Int,Nothing}=nothing,orbital::Union{Int,Nothing}=nothing)
    Couplings(SpinCoupling{1}(1,tags=('y',),atoms=scsinglewrapper(atom),orbitals=scsinglewrapper(orbital)))
end

"""
    Sᶻ(;atom::Union{Int,Nothing}=nothing,orbital::Union{Int,Nothing}=nothing) -> Couplings{ID{SCID,1},SpinCoupling{Rational{Int}/Int,<:ID{SCID,1}}}

The single Sᶻ coupling.
"""
function Sᶻ(;atom::Union{Int,Nothing}=nothing,orbital::Union{Int,Nothing}=nothing)
    Couplings(SpinCoupling{1}(1,tags=('z',),atoms=scsinglewrapper(atom),orbitals=scsinglewrapper(orbital)))
end

"""
    sc"..." -> SpinCoupling

Construct a SpinCoupling from a literal string.
"""
macro sc_str(str::String)
    ps=split(str," with ")
    condition=length(ps)==2 ? Meta.parse(ps[2]) : :nothing
    ps=split(ps[1],' ')
    coeff=eval(Meta.parse(ps[1]))
    tags=Tuple(replace(ps[2],"S"=>""))
    attrpairs=Any[:tags=>tags]
    if length(ps)==3
        components=split(ps[3],'@')
        centers=length(components)==2 ? couplingcenters(components[2]) : nothing
        centers===nothing || @assert length(centers)==length(tags) "@sc_str error: dismatched ranks."
        push!(attrpairs,:centers=>centers)
        if length(components[1])>0
            for component in split(components[1],'⊗')
                attrname,attrvalue=sccomponent(component)
                if isa(attrvalue,Expr)
                    @assert attrname==:orbitals "@sc_str error: wrong input pattern."
                    @assert length(tags)==length(attrvalue.args) "@sc_str error: dismatched ranks."
                    push!(attrpairs,attrname=>eval(subscriptexpr(attrvalue,condition)))
                else
                    @assert length(tags)==length(attrvalue) "@sc_str error: dismatched ranks."
                    push!(attrpairs,attrname=>attrvalue)
                end
            end
        end
    end
    return SpinCoupling{length(tags)}(coeff;attrpairs...)
end
function sccomponent(str::AbstractString)
    @assert str[3]=='(' && str[end]==')' "sccomponent error: wrong input pattern."
    attrname=str[1:2]=="sl" ? :atoms : str[1:2]=="ob" ? :orbitals : error("sccomponent error: wrong input pattern.")
    expr=Meta.parse(str[3:end])
    attrvalue=isa(expr,Expr) ? (all(isa(arg,Int) for arg in expr.args) ? Tuple(expr.args) : expr) : (expr,)
    return attrname=>attrvalue
end

function scpairs(str::AbstractString,::Val{R}) where R
    attrpairs=Pair{Symbol,NTuple{R,Int}}[]
    if length(str)>0
        components=split(str,'@')
        length(components)==2 && push!(attrpairs,:centers=>couplingcenters(components[2]))
        if length(components[1])>0
            for component in split(components[1],'⊗')
                attrname,attrvalue=sccomponent(component)
                @assert isa(attrvalue,NTuple{R,Int}) "scpairs error: wrong input pattern."
                push!(attrpairs,attrname=>attrvalue)
            end
        end
    end
    return attrpairs
end

"""
    heisenberg"sl(a₁,a₂)⊗ob(o₁,o₂)@(c₁,c₂)" -> Couplings

The Heisenberg couplings.
"""
macro heisenberg_str(str::String)
    ps=split(str," ")
    @assert length(ps)∈(1,2) "heisenberg_str error: not supported pattern."
    mode,str=length(ps)==2 ? (ps[1],ps[2]) : ps[1]∈("+-z","xyz") ?  (ps[1],"") : ("+-z",ps[1])
    @assert mode=="+-z" || mode=="xyz" "heisenberg_str error: not supported mode($mode)."
    attrpairs=scpairs(str,Val(2))
    if mode=="+-z"
        sc1=SpinCoupling{2}(1//2;:tags=>('+','-'),attrpairs...)
        sc2=SpinCoupling{2}(1//2;:tags=>('-','+'),attrpairs...)
        sc3=SpinCoupling{2}(1//1;:tags=>('z','z'),attrpairs...)
    else
        sc1=SpinCoupling{2}(1;:tags=>('x','x'),attrpairs...)
        sc2=SpinCoupling{2}(1;:tags=>('y','y'),attrpairs...)
        sc3=SpinCoupling{2}(1;:tags=>('z','z'),attrpairs...)
    end
    return Couplings(sc1,sc2,sc3)
end

"""
    ising"x sl(a₁,a₂)⊗ob(o₁,o₂)@(c₁,c₂)" -> Couplings
    ising"y sl(a₁,a₂)⊗ob(o₁,o₂)@(c₁,c₂)" -> Couplings
    ising"z sl(a₁,a₂)⊗ob(o₁,o₂)@(c₁,c₂)" -> Couplings

The Ising couplings.
"""
macro ising_str(str::String)
    @assert str[1] in ('x','y','z') "@ising_str error: wrong input pattern."
    attrpairs=length(str)>1 ? (@assert str[2]==' ' "@ising_str error: wrong input pattern."; scpairs(str[3:end],Val(2))) : Pair{Symbol,NTuple{2,Int}}[]
    return Couplings(SpinCoupling{2}(1;:tags=>(str[1],str[1]),attrpairs...))
end

"""
    gamma"x sl(a₁,a₂)⊗ob(o₁,o₂)@(c₁,c₂)" -> Couplings
    gamma"y sl(a₁,a₂)⊗ob(o₁,o₂)@(c₁,c₂)" -> Couplings
    gamma"z sl(a₁,a₂)⊗ob(o₁,o₂)@(c₁,c₂)" -> Couplings

The Gamma couplings.
"""
macro gamma_str(str::String)
    @assert str[1] in ('x','y','z') "@gamma_str error: wrong input pattern."
    t1,t2=str[1]=='x' ? ('y','z') : str[1]=='y' ? ('z','x') : ('x','y')
    attrpairs=length(str)>1 ? (@assert str[2]==' ' "@gamma_str error: wrong input pattern."; scpairs(str[3:end],Val(2))) : Pair{Symbol,NTuple{2,Int}}[]
    sc1=SpinCoupling{2}(1;:tags=>(t1,t2),attrpairs...)
    sc2=SpinCoupling{2}(1;:tags=>(t2,t1),attrpairs...)
    return Couplings(sc1,sc2)
end

"""
    dm"x sl(a₁,a₂)⊗ob(o₁,o₂)@(c₁,c₂)" -> Couplings
    dm"y sl(a₁,a₂)⊗ob(o₁,o₂)@(c₁,c₂)" -> Couplings
    dm"z sl(a₁,a₂)⊗ob(o₁,o₂)@(c₁,c₂)" -> Couplings

The DM couplings.
"""
macro dm_str(str::String)
    @assert str[1] in ('x','y','z') "@dm_str error: wrong input pattern."
    t1,t2=str[1]=='x' ? ('y','z') : str[1]=='y' ? ('z','x') : ('x','y')
    attrpairs=length(str)>1 ? (@assert str[2]==' ' "@dm_str error: wrong input pattern."; scpairs(str[3:end],Val(2))) : Pair{Symbol,NTuple{2,Int}}[]
    sc1=SpinCoupling{2}(1;:tags=>(t1,t2),attrpairs...)
    sc2=SpinCoupling{2}(-1;:tags=>(t2,t1),attrpairs...)
    return Couplings(sc1,sc2)
end

"""
    sˣ"sl(a)⊗ob(o)" -> Couplings
    sʸ"sl(a)⊗ob(o)" -> Couplings
    sᶻ"sl(a)⊗ob(o)" -> Couplings

The single Sˣ/Sʸ/Sᶻ coupling.
"""
macro sˣ_str(str::String) Couplings(SpinCoupling{1}(1;:tags=>('x',),scpairs(str,Val(1))...)) end
macro sʸ_str(str::String) Couplings(SpinCoupling{1}(1;:tags=>('y',),scpairs(str,Val(1))...)) end
macro sᶻ_str(str::String) Couplings(SpinCoupling{1}(1;:tags=>('z',),scpairs(str,Val(1))...)) end

"""
    SpinTerm{R}(id::Symbol,value::Any,bondkind::Any;
                couplings::Union{Function,Coupling,Couplings},
                amplitude::Union{Function,Nothing}=nothing,
                modulate::Union{Function,Bool}=false,
                ) where R

Spin term.

Type alias for `Term{'B',:SpinTerm,R,id,V,<:Any,<:TermCouplings,<:TermAmplitude,<:Union{TermModulate,Nothing}}`.
"""
const SpinTerm{R,id,V,B<:Any,C<:TermCouplings,A<:TermAmplitude,M<:Union{TermModulate,Nothing}}=Term{'B',:SpinTerm,R,id,V,B,C,A,M}
function SpinTerm{R}(id::Symbol,value::Any,bondkind::Any;
                    couplings::Union{Function,Coupling,Couplings},
                    amplitude::Union{Function,Nothing}=nothing,
                    modulate::Union{Function,Bool}=false,
                    ) where R
    isa(couplings,TermCouplings) || (couplings=TermCouplings(couplings))
    isa(amplitude,TermAmplitude) || (amplitude=TermAmplitude(amplitude))
    isa(modulate,TermModulate) || (modulate=modulate===false ? nothing : TermModulate(id,modulate))
    Term{'B',:SpinTerm,R,id}(value,bondkind,couplings,amplitude,modulate,1)
end
abbr(::Type{<:SpinTerm})=:sp
isHermitian(::Type{<:SpinTerm})=true

"""
    otype(T::Type{<:SpinTerm},I::Type{<:OID})

Get the operator type of a spin term.
"""
otype(T::Type{<:SpinTerm},I::Type{<:OID})=SOperator{T|>valtype,ID{I,T|>rank}}

end # module


text:
fbroadcast(f::Function, a::NamedTuple{K}) where K = NamedTuple{K}(f.(values(a)))
fbroadcast(f::Function, a::OrderedDict) = OrderedDict(k => f(v) for (k, v) in a)
fbroadcast(f::Function, a) = f.(a)

treemap(f::Function, t; kwargs...) = treemap((t, chs) -> f(only(t), chs), (t,); kwargs...)
treemap(f::Function, ts...; kwargs...) = treemap(f, ts; kwargs...)
function treemap(f::Function, ts::Tuple; complete::Bool=false, order::AbstractOrder=PostOrder())
    _treemap(f, ts, complete, order)
end

# TODO finish preorder, think about how both maps should behave in case of complete traversals
# maybe annotate about named tuples? With traits?
#
# select_keys(x, y) = error("Inconsistent types of children of mapped vertex and original vertices")
# select_keys(x::Tuple, y::Tuple) = x[eachindex(y)]
# function select_keys(x::NamedTuple, y::NamedTuple)
#     NamedTuple{keys(y)}(tuple((x[k] for k in keys(y))...))
# end

function _treemap(f::Function, ts::Tuple, complete::Bool, order::PreOrder)
    error("Treemaps using PreOrder() are not supported yet")
end

# function _treemap(f::Function, ts::Tuple, complete::Bool, order::PreOrder)
#     n = f(ts, _children_sorted.(ts))
#     isleaf(n) && return n
#     nchs = children(n)
#     chs = _children_pairs_keys(ts, complete)
#     @show keys(nchs), keys(chs)
#     @assert keys(nchs) ⊆ keys(chs)
#     mchs = fbroadcast(ts -> _treemap(f, ts, complete, order), select_keys(chs, nchs))
#     set_children(n, mchs)
# end

function _treemap(f::Function, ts::Tuple, complete::Bool, order::PostOrder)
    chs = fbroadcast(ts -> _treemap(f, ts, complete, order), _children_pairs(ts, complete))
    return f(ts, chs)
end

function _treemap(f::Function, ts::Tuple, complete::Bool, order::LevelOrder)
    error("Treemaps using LevelOrder() are not supported yet")
end

function treemap!(f::Function, ts...; complete::Bool=false, order::AbstractOrder=PostOrder())
    foreach(f, NodeIterator(ts...; order=order, complete=complete))
end
function leafmap!(f::Function, ts...; complete::Bool=false, order::AbstractOrder=PostOrder())
    foreach(f, LeafIterator(ts...; order=order, complete=complete))
end
function typemap!(f::Function, t::Union{Type, Tuple{Vararg{Type}}}, ts...; complete::Bool=false,
                  order::AbstractOrder=PostOrder())
    foreach(f, TypeIterator(t, ts...; order=order, complete=complete))
end
function predicatemap!(f::Function, pred::Function, ts...; complete::Bool=false,
                  order::AbstractOrder=PostOrder())
    foreach(f, PredicateIterator(pred, ts...; order=order, complete=complete))
end

ext: jl


text:
export async

import Base: close
import Base: show

"""
    async(size::Int = typemax(Int))

Creates an async operator, which sends items from the source Observable asynchronously.

# Arguments
- `size`: Asynchronous messages buffer size, default is a `typemax(Int)`

# Producing

Stream of type `<: Subscribable{L}` where `L` refers to type of source stream

See also: [`AbstractOperator`](@ref), [`InferableOperator`](@ref), [`ProxyObservable`](@ref)
"""
async(size::Int = typemax(Int)) = schedule_on(AsyncScheduler(size))


text:
@testset "objective_coefficients.jl" begin

    @testset "set_dual_model_sense" begin
        # ERROR: FEASIBILITY_SENSE is not supported
        # @test_throws ErrorException Dualization.set_dual_model_sense(lp11_test(), lp11_test())
        # obj = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(0.0, VI(1))], 0.0)
        # @test_throws ErrorException Dualization.PrimalObjective{Float64}(obj)

        model = lp1_test()
        @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MIN_SENSE
        Dualization.set_dual_model_sense(model, model)
        @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE

        model = lp4_test()
        @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
        Dualization.set_dual_model_sense(model, model)
        @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MIN_SENSE
    end

    @testset "get_primal_objective" begin
        model = lp1_test()
        primal_objective = Dualization.get_primal_objective(model)

        @test Dualization.get_affine_terms(primal_objective)[1] == MOI.ScalarAffineTerm{Float64}(-4.0, MOI.VariableIndex(2))
        @test MOI.constant(Dualization.get_raw_obj(primal_objective)) == -1.0

        model = lp10_test()

        @test model.objective == MOI.SingleVariable(MOI.VariableIndex(1))
        primal_objective = Dualization.get_primal_objective(model)
        @test Dualization.get_affine_terms(Dualization.get_primal_objective(model))[1] == MOI.ScalarAffineTerm{Float64}(1.0, MOI.VariableIndex(1))
        @test MOI.constant(Dualization.get_raw_obj(primal_objective)) == 0.0
    end

end


text:
function processexpr!(expr) end # fallback

function processexpr!(expr::Expr)
    head = expr.head
    args = expr.args

    if head === :call && (args[1] === :|> || args[1] === :.|>)
        # process inside `|>` operators first
        processexpr!(args[2])

        if args[3] isa Expr
            if args[3].head === :call
                # inject into the first argument
                args[3] = makechainfunc(args[3], 1)
            elseif args[3].head === :tuple
                chainarg, fcall = args[3].args
                args[3] = makechainfunc(fcall, chainarg)
            end
        end
    end
end

@doc """
    makechainfunc

manipulate chanining expression
"""
function makechainfunc end

function makechainfunc(fcall::Symbol, chainarg::Int)
    f = fcall
    :( chained -> $f(chained) )
end

function makechainfunc(fcall::Symbol, chainarg::QuoteNode)
    f = fcall
    :( chained -> $f(; $(chainarg.value) = chained) )
end

function makechainfunc(fcall::Expr, chainarg::Int)
    offset = haskeywordarg(fcall) ? 2 : 1
    insert!(fcall.args, chainarg + offset, :chained)

    :( chained -> $fcall )
end

function makechainfunc(fcall::Expr, chainarg::QuoteNode)
    keyarg = Expr(:kw, chainarg.value, :chained)
    if haskeywordarg(fcall)
        push!(fcall.args[2].args, keyarg)
    else
        insert!(fcall.args, 2, Expr(:parameters, keyarg))
    end

    :( chained -> $fcall )
end

function haskeywordarg(fcall)
    for arg ∈ fcall.args
        arg isa Expr && arg.head === :parameters && return true
    end
    return false
end

# @TODO: absolutely better to manipulate ASTs directly rather than parsing from string
function concatexpr(exprs, dotfuse = false)
    concated = ["@>"]
    push!(concated, stringify(exprs[1]))
    for expr ∈ exprs[2:end]
        if dotfuse
            push!(concated, ".|>")
        else
            push!(concated, "|>")
        end
        push!(concated, stringify(expr))
    end

    join(concated, " ") |> Meta.parse
end

stringify(expr) = string(expr)
stringify(expr::String) = "\"$(expr)\""


text:
abstract type AbstractOdeSolver
end

abstract type AbstractAdaptiveExplicitSolver <: AbstractOdeSolver
end

abstract type AbstractRosenbrockSolver <: AbstractOdeSolver
end

function order(::S) where {S <: AbstractOdeSolver}
    error("To be implemented")
end
function n_intermediates(::S) where {S <: AbstractOdeSolver}
    error("To be implemented")
end
function stepper!(ode_solver::S, ode_fcn, stepper_options, stepper_temp::StepperTemporaries) where {S <: AbstractOdeSolver}
    error("To be implemented")
end

# some errors with those: probably wrong formulas
# include("solver_heuneuler.jl")
# include("solver_ros34.jl")

include("solver_bsh23.jl")
include("solver_dp45.jl")
include("solver_ts45.jl")
include("solver_ros23.jl")

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
include("IOM/src/share/LogSystem.jl")
include("IOM/src/share/PolelikeCoordinate.jl")
include("IOM/src/models/EMOM/ENGINE_EMOM.jl")
include("IOM/src/driver/driver_working.jl")
include("IOM/src/share/CyclicData.jl")

using MPI
using CFTime, Dates
using ArgParse

using .PolelikeCoordinate
using .LogSystem
using .CyclicData
function parse_commandline()

    s = ArgParseSettings()

    @add_arg_table s begin

        "--config-file"
            help = "Configuration file."
            arg_type = String
            required = true

        "--forcing-file"
            help = "Annual forcing file. It should contain: TAUX, TAUY, SWFLX, NSWFLX, VSFLX"
            arg_type = String
            required = true

        "--nudging-timescale-day"
            help = "Nudging timescale used to replace weak restoring timescale when finding Q-flux in days. Default is 5 days."
            arg_type = Float64
            default = 5.0

        "--stop-n"
            help = "Core of the model."
            arg_type = Int64

        "--time-unit"
            help = "."
            arg_type = String

    end

    return parse_args(s)
end

parsed = parse_commandline()

include(parsed["config-file"])

time_unit = Dict(
    "year" => Dates.Year,
    "month" => Dates.Month,
    "day"   => Dates.Day,
)[parsed["time-unit"]]

t_simulation = time_unit(parsed["stop-n"])

Δt = Dates.Second(86400)
read_restart = false

cfgmc = config[:MODEL_CORE]
cfgmm = config[:MODEL_MISC]

cfgmc[:weak_restoring] = :on
cfgmc[:τwk_TEMP] = parsed["nudging-timescale-day"] * 86400.0
cfgmc[:τwk_SALT] = parsed["nudging-timescale-day"] * 86400.0

t_start = DateTimeNoLeap(1, 1, 1, 0, 0, 0)
t_end = t_start + t_simulation

MPI.Init()
comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)
is_master = rank == 0

coupler_funcs = (

    master_before_model_init = function()
        global cdatam = CyclicDataManager(;
            timetype     = getproperty(CFTime, Symbol(cfgmm[:timetype])),
            filename     = parsed["forcing-file"],
            varnames     = ["TAUX", "TAUY", "SWFLX", "NSWFLX", "VSFLX"],
            beg_time     = DateTimeNoLeap(1, 1, 1),
            end_time     = DateTimeNoLeap(2, 1, 1),
            align_time   = DateTimeNoLeap(1, 1, 1),
        )

        global datastream = makeDataContainer(cdatam)

        return read_restart, t_start
    end,

    master_after_model_init! = function(OMMODULE, OMDATA)
            # setup forcing

    end,

    master_before_model_run! = function(OMMODULE, OMDATA)

        global datastream

        t_end_reached = OMDATA.clock.time >= t_end

        if ! t_end_reached

            interpData!(cdatam, OMDATA.clock.time, datastream)
            OMDATA.x2o["SWFLX"]       .= datastream["SWFLX"]
            OMDATA.x2o["NSWFLX"]      .= datastream["NSWFLX"]
            OMDATA.x2o["VSFLX"]       .= datastream["VSFLX"]
            OMDATA.x2o["TAUX_east"]   .= datastream["TAUX"]
            OMDATA.x2o["TAUY_north"]  .= datastream["TAUY"]

            return_values = ( :RUN,  Δt, t_end_reached )
        else
            return_values = ( :END, 0.0, t_end_reached  )
        end

        return return_values
    end,

    master_after_model_run! = function(OMMODULE, OMDATA)
    end,

    master_finalize! = function(OMMODULE, OMDATA)
        writeLog("[Coupler] Finalize")
    end,
)

runModel(
    ENGINE_EMOM,
    coupler_funcs,
    t_start,
    t_simulation,
    read_restart,
    config,
)

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
export BatchExplorer

"""
    BatchExplorer(explorer::AbstractExplorer)
    BatchExplorer(explorers::AbstractExplorer...)
"""
struct BatchExplorer{E} <: AbstractExplorer
    explorer::E
end

BatchExplorer(explorers::AbstractExplorer...) = BatchExplorer(explorers)

(x::BatchExplorer)(values::AbstractMatrix) = [x.explorer(v) for v in eachcol(values)]
(x::BatchExplorer{<:Tuple})(values::AbstractMatrix) =
    [explorer(v) for (explorer, v) in zip(x.explorer, eachcol(values))]


text:
#using PyPlot
#pygui(true)
using Random
using LinearAlgebra
using Plots

function mySliceSampler(pdf_log, x_0, w, m, N)
    # REQUIRES MODULES:
    # Random
    #
    # INPUT:
    # pdf_log --> log-pdf of the target distribution (a function)
    # x_init --> inital vector of x (dimension 1xD), where D dimension of pdf
    # w --> typical slice size, for a window (dimension 1xD)
    # m --> integer limiting the slice size to m*w
    # N --> number of sample points
    #
    # OUTPUT:
    # x_s --> sampled x values (NxD), only keep values after full permutations along dimensions
    # pdflog_x_s --> log-pdf of sampled values

    D = length(x_0)  # the dimension of the distribution (nr of coordinates)
    x_s = Array{Float64}(undef,N,D)   # Samples (NxD), preallocate array, with undefined values
    pdflog_x_s = Array{Float64}(undef,N)   # (unnormalised) log-prob of sampled points (Nx1)
                                           # for reducing evaluations
    pdflog_x_1 = 0.0;

    for ii in 1:N
        # Update the new x progressively for each dimension:
        L   = 1*x_0  # R and L need to be of same dimension, because we apply the pdf upon those!
        R   = 1*x_0
        x_1 = 1*x_0  # "1*" needed, otherwise x_0 suddendly takes mysterious values
        for d in randperm(D)
            ### 1) Make the Slice
            # random "vertical" position
            if (ii == 1)
                #vertical = rand() * pdf_log(x_0) # in normal space
                z = pdf_log(x_0) + log(rand()) # in log-space
            else
                # saving evaluations, pdflog_x_1 was updated for current permutation in the end
                z = pdflog_x_1 + log(rand())
            end

            ### 2) Make the window, "Stepping-out" alogrithm
            # Randomly place a window, containing x_0:
            L[d] = x_0[d] - rand() * w[d]
            R[d] = L[d] + w[d]

            # Randomly share the max window size among left/right:
            J = floor(m*rand());
            K = (m-1) - J;

            # Extend window to the left, until outside slice or allowance seizes:
            while ((J > 0) && (z < pdf_log(L)))
                #println("Lefting")
                L[d] -= w[d]
                J -= 1
            end

            # Extend window to the right, until outside slice or allowance seizes:
            while ((K > 0) && (z < pdf_log(R)))
                #println("Righting")
                R[d] += w[d]
                K -= 1
            end

            ### 3) Sample from window
            # finding an allowable point:
            while true
                x_1[d] = L[d] + rand() * (R[d] - L[d])

                # this + breaking out of loop reduces the amount of log-pdf evaluations:
                pdflog_x_1 = pdf_log(x_1)

                if (pdflog_x_1 >= z)     # new value found
                    x_0[d] = x_1[d]     # update value for this dimension
                    break
                # Value was not within slice, shrink the interval:
                elseif (x_1[d] < x_0[d])
                    L[d] = x_1[d]
                elseif (x_1[d] > x_0[d])
                    R[d] = x_1[d]
                else
                    throw(ErrorException("Error during shrinking the interval"))
                end
            end

        end # end of permutating through dimensions

        ### 4) Update the chain:
        x_s[ii,:] = x_1    # sample point, after one round of permutations
        pdflog_x_s[ii] = pdflog_x_1 # now contains log_pdf of all new coordinates of new point, since has been updated at each permutation
        # note that x_0 has been updated during the permutations

    end

    return x_s, pdflog_x_s
end

#%%

using LinearAlgebra

function mylogNormal(x)
    x = x[1]
    m = 1;
    s = 1;
    y = log(1/sqrt(2*pi*s^2)) + (-(x-m)^2 / (2*s^2))
    return y
end

function mylogNormal_bivariate(x)
    # Requires LinearAlgebra
    X = [0.0; 0.0]
    X[1] = x[1]
    X[2] = x[2]

    M = [1;
         1]
    S = [1 0;
         0 1]

    y = log(1 / sqrt(2*pi)) +
        log(1 / sqrt(det(S))) +
        (-1/2) * (X-M)' * inv(S) * (X-M)

    return y
end

function myCrazyDistribution(x)
    if x < 1
        y = 0
    elseif 1 <= x <= 2
        y = 1
    elseif 2 < x < 3
        y = 0
    elseif 3 <= x <= 4
        y = 1
    else
        y = 0
    end
    return y
end

function log_square(x)
    x=x[1]
    if 0 < x < 4
        y = 2*log(x)
    else
        y = -Inf
    end
    return y
end

function log_uniform(x)
    x = x[1]
    if 0 < x < 1
        y = 0
    else
        y = -Inf
    end
end

#%%
N = 1000000
#x_start = [0.5];
x_start = [0.5 0.5];
w = [.1 0.1]; # typical window width
m = 100;

using Distributions
bivarLogNormal(X) = log(pdf(MvNormal([1.0;1.0],[1. 0.; 0. 1.]), [X[1];X[2]]))

bivarLogNormal([1 1])

#%%
#x_s, log_vals = mySliceSampler(log_uniform, x_start, w, m, N)
#x_s, log_vals = mySliceSampler(log_square, x_start, w, m, N)
#x_s, log_vals = mySliceSampler(mylogNormal, x_start, w, m, N)
#x_s, log_vals = mySliceSampler(bivarLogNormal, x_start, w, m, N)
x_s, log_vals = mySliceSampler(mylogNormal_bivariate, x_start, w, m, N)

mean = sum(x_s,dims=1)/N

#%%
#figure()
#hist(x_s[:,2],bins=100)
histogram2d(x_s[:,1],x_s[:,2])#bins=100)

#%%

x = -2:.1:4
y = -2:.1:4
f(x,y) = exp(mylogNormal_bivariate([x;y]))
pyplot()
plot(x,y,f,st=:surface)

#%%


text:
#
# Power spectrum
#

immutable PowerSpectrum <: Statistic; end

Base.eltype{T<:Real}(::PowerSpectrum, X::AbstractVecOrMat{Complex{T}}) = T
allocwork{T<:Complex}(::PowerSpectrum, X::AbstractVecOrMat{Complex{T}}) = nothing
allocoutput{T<:Real}(::PowerSpectrum, X::AbstractVecOrMat{Complex{T}}) =
    Array(T, 1, nchannels(X))

# Single input matrix
computestat!{T<:Real}(::PowerSpectrum, out::AbstractMatrix{T}, ::Void,
                      X::AbstractVecOrMat{Complex{T}}) =
    scale!(sumabs2!(out, X), 1/ntrials(X))

#
# Cross spectrum
#

immutable CrossSpectrum <: PairwiseStatistic; end
Base.eltype{T<:Real}(::CrossSpectrum, X::AbstractArray{Complex{T}}) = Complex{T}

# Single input matrix
allocwork{T<:Complex}(::CrossSpectrum, X::AbstractVecOrMat{T}) = nothing
computestat!{T<:Complex}(::CrossSpectrum, out::AbstractMatrix{T}, ::Void,
                         X::AbstractVecOrMat{T}) =
    scale!(Ac_mul_A!(out, X), 1/ntrials(X))

# Two input matrices
allocwork{T<:Complex}(::CrossSpectrum, X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T}) = nothing
computestat!{T<:Complex}(::CrossSpectrum, out::AbstractMatrix{T}, ::Void,
                         X::AbstractVecOrMat{T}, Y::AbstractVecOrMat{T}) =
    scale!(Ac_mul_B!(out, X, Y), 1/ntrials(X))

accumulator{T<:Real}(::Type{CrossSpectrum}, ::Type{T}) = zero(Complex{T})
@inline accumulate{T<:Real}(::Type{CrossSpectrum}, x::Complex{T},
                            v1::Complex{T}, v2::Complex{T}) = (x + conj(v1)*v2)
@inline accumulate{T<:Real}(::Type{CrossSpectrum}, x::Complex{T},
                            v1::Complex{T}, v2::Complex{T}, weight::Real) = (x + conj(v1)*v2*weight)
finish(::Type{CrossSpectrum}, x::Complex, n::Int) = x/n


text:
using Distributed
using JLD2, FileIO

# addprocs(Sys.CPU_THREADS-1)
# addprocs(Int(Sys.CPU_THREADS/2)-1)

@everywhere using Pkg
@everywhere Pkg.activate("./Module1")
@everywhere using Module1

# for i = 1:10
#     gapsdict10xi = primegapsdict(10^i)
#     @save "gapsdict10x$(i).jld2" gapsdict10xi
# end

# rmprocs(workers())

# gapsdict10x12 = primegapsdict(10^12)
# @save "gapsdict10x12.jld2" gapsdict10x12

max_stars_repo_name: chortbauer/primegaps

max_issues_repo_name: chortbauer/primegaps

max_forks_repo_name: chortbauer/primegaps


text:
using DistributionVectors
using Test, Distributions, LinearAlgebra, Missings
using BenchmarkTools: @btime
using Plots, StatsPlots

function benchmarkSums()
    nrep = 30
    mu = log.(rand(Normal(120, 10), nrep));
    sigma = log.(rand(Normal(1.2, 0.1), nrep));
    acf1 = [0.8,0.3,0.1];
    dv = SimpleDistributionVector(LogNormal{eltype(mu)}, mu, sigma)
    @btime dsum_v = sum($dv, $acf1, SkipMissing(); method = Val(:vector) )
    @btime dsum_m = sum($dv, $acf1, SkipMissing(); method = Val(:bandedmatrix) )
    # the Bandematrix based slower after optimizing the vector for stability in missings
    #
    # repeat with missings
    mum = allowmissing(copy(mu)); mum[1] = missing
    dv = SimpleDistributionVector(LogNormal{nonmissingtype(eltype(mum))}, mum, sigma)
    @btime dsum_v = sum($dv, $acf1, SkipMissing(); method = Val(:vector) )
    @btime dsum_m = sum($dv, $acf1, SkipMissing(); method = Val(:bandedmatrix) )
    #
    #S = similar(mum);
    @code_warntype sum_lognormals(S, dv, acf1, SkipMissing())
    f() = sum_lognormals(S, dv, acf1, SkipMissing())
    @time f(); @time f()

    # # try allocating instead of view (replace line of view_nonmissing)
    # # allocating is faster
    # function sum_lognormals2!(S, dv, corr::AbstractMatrix,
    # ms::MissingStrategy=PassMissing())
    #     parms = params(dv)
    #     μ = @view parms[1,:]
    #     σ = @view parms[2,:]
    #     # S = allowmissing(similar(μ))
    #     @. S = exp(μ + abs2(σ)/2)
    #     nmissing = count(ismissing, S)
    #     anymissing = nmissing != 0
    #     !(isa(ms, HandleMissingStrategy)) && anymissing && error(
    #          "Found missing values. Use argument 'SkipMmissing()' to sum over nonmissing.")
    #     Ssum::nonmissingtype(eltype(S)) = sum(skipmissing(S))
    #     @. S = σ * S  # do only after Ssum
    #     # setting S to zero results in summing zero for missing records
    #     # which is the same as filtering both S and corr
    #     anymissing && replace!(S, missing => 0.0)
    #     #s = transpose(disallowmissing(S)) * corr * disallowmissing(S)
    #     #Spure = view_nonmissing(S) # non-allocating
    #     Spure = disallowmissing(S) # allocating
    #     s = transpose(Spure) * corr * Spure
    #     σ2eff = s/abs2(Ssum)
    #     μ_sum = log(Ssum) - σ2eff/2
    #     #@show Ssum, s, length(S) - nmissing
    #     LogNormal(μ_sum, √σ2eff)
    # end
    # @btime sum_lognormals($storage, $dv, $corMa, SkipMissing())
    # @btime sum_lognormals2!($storage, $dv, $corMa, SkipMissing())
end

function bootstrap_sums_lognormal()
    #using StatsBase
    nObs = 5
    xTrue = fill(10.0, nObs)
    sigmaStar = fill(1.5, nObs) # multiplicative stddev of 1.5
    dv = SimpleDistributionVector(fit.(LogNormal, xTrue, Σstar.(sigmaStar))...)
    dsum = sum(dv)
    nboot = 100_000
    x = rand(dv, nboot);
    sums = vec(sum(x, dims = 2));
    #stds = std.(Ref(x));
    cdf_sums = ecdf(sums);
    @test cdf_sums(median(dsum)) ≈ 0.5 rtol = 0.02
    @test cdf_sums(quantile(dsum, [0.025,0.975])) ≈ [0.025, 0.975] rtol = 0.02
end

function fplot(dsum, sums)
    plot(dsum, lab="computed", xlabel="sum", ylabel="density")
    density!(sums, lab="random sample")
    vline!([mean(dsum)], lab="mean computed")
    vline!([mean(sums)], lab="mean random")
    vline!(quantile(dsum, [0.025, 0.975]), lab="cf computed")
    vline!(quantile(sums, [0.025, 0.975]), lab="cf random")
end

function bootstrap_sums_lognormal_acf()
    #using StatsBase
    mu = log.([110,100,80,120,160.0])
    sigma = log.([1.2,1.5,1.1,1.3,1.1])
    acf1 = [0.4,0.1]
    dv = SimpleDistributionVector(LogNormal{eltype(mu)}, mu, sigma);
    #ps = 0.025:0.025:0.975; qs = quantile.(dv[1], ps); plot(qs, pdf(dv[1],qs))
    drsum = boot_dvsums_acf(dv, acf1)
    dsum = sum(dv, acf1)
    @test dsum ≈ drsum rtol = 0.02
    # see plotting function in bootstrap_sums_lognormal above
end


text:
function test_sparsity(x, target_sparsity; atol = 0.05)
    number_of_zeros = count(ele -> (ele == 0.0), x)
    sparsity = number_of_zeros / sum(length, x)
    @test isapprox(target_sparsity, sparsity; atol)
end

@testset "MCDense GPU" begin
    dropout_rate = 0.35
    # Test MC Dense layer
    a = gpu(rand(Float32, 8, 32))
    layer = gpu(MCDense(8, 16, dropout_rate))
    output = layer(a)
    # Test if it's CuArray
    @test output isa CuArray
    @test isequal(size(output), (16, 32))
    test_sparsity(output, dropout_rate)
    # Test MC dense dropout toggle
    output = layer(a, dropout = false)
    test_sparsity(output, 0)
end

@testset "MCConv GPU" begin
    dropout_rate = 0.4
    # Test MC conv layer
    a = gpu(rand(Float32, 32, 32, 3, 32))
    layer = MCConv((5, 5), 3 => 6, dropout_rate) |> gpu
    output = layer(a)
    # Test if it's CuArray
    @test output isa CuArray
    # Test the output shape
    @test isequal(size(output), (28, 28, 6, 32))
    # Test the sparsity percentage in the array
    test_sparsity(output, dropout_rate)
    # Test MC conv dropout toggle
    output = layer(a, dropout = false)
    test_sparsity(output, 0)
end


text:
#=

  Euler #42 in Julia.

  Problem 42:
  """
  The nth term of the sequence of triangle numbers is given by,
      tn = 1/2*n*(n+1);
  so the first ten triangle numbers are:

  1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...

  By converting each letter in a word to a number corresponding to its
  alphabetical position and adding these values we form a word value. For example,
  the word value for SKY is 19 + 11 + 25 = 55 = t10. If the word value
  is a triangle number then we shall call the word a triangle word.

  Using words.txt (right click and 'Save Link/Target As...'), a 16K text file
  containing nearly two-thousand common English words, how many
  are triangle words?
  """

  This Julia model was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

# n'th triangle number
function triangle_number(n)
    return floor(Int,(n*(n-1)) / 2)
end

# get the score for a name
function get_score(name)
    return ([name[i] for i in 1:length(name)]
           .|>c->(Int(c)-64))|>sum
end

# 0.09908583s
function euler42a()
    t20 = Set((1:20).|>i->triangle_number(i))
    words = split( (readline("euler42_words.txt").|>word->replace(word,"\""=>"")),",")
    return (words
           .|>word-> get_score(word) in t20
           )|>sum
end

run_euler(euler42a)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
function _show_acronym_component_html(io::IO, c::String, n::Int)
    a, b = splitword(c, n)
    print(io, "", uppercase(a), "", b)
end

function Base.show(io::IO, ::MIME"text/html", a::Acronym)
    p = 1
    for (i, (comp, n)) in enumerate(zip(a.components, a.nvals))
        if isnothing(comp)
            comp2 = a.word[p:(p+n-1)]
            _show_acronym_component_html(io, comp2, n)
            print(io, "???")
        else
            _show_acronym_component_html(io, comp, n)
        end
        i != length(a.components) && print(io, " ")

        p += n
    end
end

function Base.show(io::IO, ::MIME"text/html", af::AcronymFamily)
    println(io, "
")

    println(io, "
")
    print(io, "
")
    for letter in af.word
        print(io, "", uppercase(letter), "	")
    end
    print(io, "
")
    println(io, "
")

    println(io, "
")
    for (i, row) in enumerate(eachrow(af.components))
        println(io, "
")

        for words in row
            println(io, "")

            isempty(words) && print(io, "?")

            for word in words
                _show_acronym_component_html(io, word, i)
                print(io, "
 ")
            end

            println(io, "	")
        end

        println(io, "
")
    end
    println(io, "
")

    println(io, "
")
end

function Base.show(io::IO, ::MIME"text/html", saf::ScoredAcronymFamily)
    af = saf.family
    l = size(af.components, 2)

    println(io, "
")

    println(io, "
")

    print(io, "
")
    print(io, "", uppercase(af.word), "	")
    println(io, "", saf.best_score, "	")
    println(io, "
")

    print(io, "
")
    for letter in af.word
        print(io, "", uppercase(letter), "	")
    end
    println(io, "
")

    println(io, "
")

    println(io, "
")
    for (i, row) in enumerate(eachrow(af.components))
        println(io, "
")

        for words in row
            println(io, "")

            isempty(words) && println(io, "?")

            for word in sort(words, by=word -> score_word(word, saf.params), rev=true)
                _show_acronym_component_html(io, word, i)
                print(io, "
 ")
            end

            println(io, "	")
        end

        println(io, "
")
    end
    println(io, "
")

    println(io, "
")
end


text:
## Turing Complete Interpreter
### A Very Tiny Turing Complete Programming Language defined with denotational semantics

# semantica dalle dispense degano
using Metatheory

import Base.ImmutableDict
Mem = Dict{Symbol,Union{Bool,Int}}

read_mem = @theory v σ begin
	(v::Symbol, σ::Mem) => if v == :skip σ else σ[v] end
end

@testset "Reading Memory" begin
	ex = :((x), $(Mem(:x => 2)))
	@test true == areequal(read_mem, ex, 2)
end

arithm_rules = @theory a b σ begin
	(a + b, σ::Mem) 	--> (a, σ) + (b, σ)
	(a * b, σ::Mem) 	--> (a, σ) * (b, σ)
	(a - b, σ::Mem) 	--> (a, σ) - (b, σ)
	(a::Int, σ::Mem) 	--> a
	(a::Int + b::Int) 	=> a + b
	(a::Int * b::Int) 	=> a * b
	(a::Int - b::Int) 	=> a - b
end

@testset "Arithmetic" begin
	@test areequal(read_mem ∪ arithm_rules, :((2 + 3), $(Mem())), 5)
end

# don't need to access memory
bool_rules = @theory a b σ begin
	(a < b, σ::Mem) --> (a, σ) < (b, σ)
	(a ∨ b, σ::Mem) --> (a, σ) ∨ (b, σ)
	(a ∧ b, σ::Mem) --> (a, σ) ∧ (b, σ)
	(¬(a), σ::Mem) 	--> ¬((a, σ))

	(a::Bool, σ::Mem) 		=> a
	(¬a::Bool) 				=> !a
	(a::Bool ∨ b::Bool) 	=> (a || b)
	(a::Bool ∧ b::Bool) 	=> (a && b)
	(a::Int < b::Int) 		=> (a < b)
end

t = read_mem ∪ arithm_rules ∪ bool_rules

@testset "Booleans" begin
	@test areequal(t, :((false ∨ false), $(Mem())), false)

	exx = :((false ∨ false) ∨ ¬(false ∨ false), $(Mem(:x => 2)))
	g = EGraph(exx)
	saturate!(g, t)
	ex = extract!(g, astsize)
	@test ex == true
	params = SaturationParams(timeout=12)
	@test areequal(t, exx, true; params=params)

	@test areequal(t, :((2 < 3) ∧ (3 < 4), $(Mem(:x => 2))), true)
	@test areequal(t, :((2 < x) ∨ ¬(3 < 4), $(Mem(:x => 2))), false)
	@test  areequal(t, :((2 < x) ∨ ¬(3 < 4), $(Mem(:x => 4))), true)
end

if_rules = @theory guard t f σ begin
	(if guard; t end) --> (if guard; t else :skip end)
	(if guard; t else f end, σ::Mem) --> (if (guard, σ); t else f end, σ)
	(if true; t else f end, σ::Mem) --> (t, σ)
	(if false; t else f end, σ::Mem) --> (f, σ)
end

if_language = read_mem ∪ arithm_rules ∪ bool_rules ∪ if_rules

@testset "If Semantics" begin
	@test areequal(if_language, 2, :(if true x else 0 end, 	$(Mem(:x => 2))))
	@test areequal(if_language, 0, :(if false x else 0 end, 	$(Mem(:x => 2))))
	@test areequal(if_language, 2, :(if ¬(false) x else 0 end, $(Mem(:x => 2))))
	params = SaturationParams(timeout=10)
	@test areequal(if_language, 0, :(if ¬(2 < x) x else 0 end, $(Mem(:x => 3))); params=params)
end

while_rules = @theory a b σ begin
	(:skip, σ::Mem) 	--> σ
	((a; b), σ::Mem) 	--> ((a, σ); b)
	(a::Int; b) 		--> b
	(a::Bool; b) 		--> b
	(σ::Mem; b) 		--> (b, σ)
	(while a b end, σ::Mem) -->
		(if a; (b; while a b end) else :skip end, σ)
end

write_mem = @theory sym val σ begin
	(sym::Symbol = val, σ::Mem) --> (sym = (val, σ), σ)
	(sym::Symbol = val::Int, σ::Mem) => merge(σ, Dict(sym => val))
end

while_language = if_language ∪ write_mem ∪ while_rules;

@testset "While Semantics" begin
	exx = :((x = 3), $(Mem(:x => 2)))
	g = EGraph(exx)
	saturate!(g, while_language)
	ex = extract!(g, astsize)

	@test areequal(while_language, Mem(:x => 3), exx)

	exx = :((x = 4; x = x + 1), $(Mem(:x => 3)))
	g = EGraph(exx)
	saturate!(g, while_language)
	ex = extract!(g, astsize)

	params = SaturationParams(timeout=10)
	@test areequal(while_language, Mem(:x => 5), exx; params=params)

	params = SaturationParams(timeout=14)
	exx = :((if x < 10 x = x + 1 else skip end), $(Mem(:x => 3)))
	@test areequal(while_language, Mem(:x => 4), exx; params=params)

	exx = :(
		(while x < 10
			x = x + 1
		end; x), $(Mem(:x => 3)))
	g = EGraph(exx)
	params = SaturationParams(timeout=100, scheduler=Schedulers.SimpleScheduler)
	saturate!(g, while_language, params)
	@test 10 == extract!(g, astsize)
end


text:
# ---
# title: 976. Largest Perimeter Triangle
# id: problem976
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math, Sort
# link: <https://leetcode.com/problems/largest-perimeter-triangle/description/>
# hidden: true
# ---
#
# Given an array `A` of positive lengths, return the largest perimeter of a
# triangle with **non-zero area** , formed from 3 of these lengths.
#
# If it is impossible to form any triangle of non-zero area, return `0`.
#
#
#
# **Example 1:**
#
#
#
#     Input: [2,1,2]
#     Output: 5
#
#
# **Example 2:**
#
#
#
#     Input: [1,2,1]
#     Output: 0
#
#
# **Example 3:**
#
#
#
#     Input: [3,2,3,4]
#     Output: 10
#
#
# **Example 4:**
#
#
#
#     Input: [3,6,2,3]
#     Output: 8
#
#
#
#
# **Note:**
#
#   1. `3 <= A.length <= 10000`
#   2. `1 <= A[i] <= 10^6`
#
#
## @lc code=start
using LeetCode

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


text:
# Extension operator: joExtend

## helper module
module joExtend_etc
    using JOLI: jo_convert, LocalVector, LocalMatrix
    ### zeros
    function zeros_fwd(v::LocalVector{VT},n::Integer,pad_upper::Integer,pad_lower::Integer,rdt=DataType) where {VT<:Number}
        w = [zeros(VT,pad_lower); v; zeros(VT,pad_upper)]
        return jo_convert(rdt,w,false)
    end
    function zeros_fwd(v::LocalMatrix{VT},n::Integer,pad_upper::Integer,pad_lower::Integer,rdt=DataType) where {VT<:Number}
        w = [zeros(VT,pad_lower,size(v,2)); v; zeros(VT,pad_upper,size(v,2))]
        return jo_convert(rdt,w,false)
    end
    function zeros_tran(v::LocalVector{VT},n::Integer,pad_upper::Integer,pad_lower::Integer,rdt=DataType) where {VT<:Number}
        w = v[pad_lower+1:n+pad_lower]
        return jo_convert(rdt,w,false)
    end
    function zeros_tran(v::LocalMatrix{VT},n::Integer,pad_upper::Integer,pad_lower::Integer,rdt=DataType) where {VT<:Number}
        w = v[pad_lower+1:n+pad_lower,:]
        return jo_convert(rdt,w,false)
    end
    ### border
    function border_fwd(v::LocalVector{VT},n::Integer,pad_upper::Integer,pad_lower::Integer,rdt=DataType) where {VT<:Number}
        w = [repeat([v[1]],pad_lower); v; repeat([v[end]],pad_upper)]
        return jo_convert(rdt,w,false)
    end
    function border_fwd(v::LocalMatrix{VT},n::Integer,pad_upper::Integer,pad_lower::Integer,rdt=DataType) where {VT<:Number}
        w = [repeat(v[1:1,:],pad_lower,1); v; repeat(v[end:end,:],pad_upper,1)]
        return jo_convert(rdt,w,false)
    end
    function border_tran(v::LocalVector{VT},n::Integer,pad_upper::Integer,pad_lower::Integer,rdt=DataType) where {VT<:Number}
        w = v[pad_lower+1:end-pad_upper]
        w[1] += sum(v[1:pad_lower])
        w[end] += sum(v[end-pad_upper+1:end])
        return jo_convert(rdt,w,false)
    end
    function border_tran(v::LocalMatrix{VT},n::Integer,pad_upper::Integer,pad_lower::Integer,rdt=DataType) where {VT<:Number}
        w = v[pad_lower+1:end-pad_upper,:]
        w[1:1,:] += sum(v[1:pad_lower,:],dims=1)
        w[end:end,:] += sum(v[end-pad_upper+1:end,:],dims=1)
        return jo_convert(rdt,w,false)
    end
    ### mirror
    function mirror_fwd(v::LocalVector{VT},n::Integer,pad_upper::Integer,pad_lower::Integer,rdt=DataType) where {VT<:Number}
        w = [reverse(v[1:pad_lower]); v; reverse(v[end-pad_upper+1:end])]
        return jo_convert(rdt,w,false)
    end
    function mirror_fwd(v::LocalMatrix{VT},n::Integer,pad_upper::Integer,pad_lower::Integer,rdt=DataType) where {VT<:Number}
        w = [flipdim(v[1:pad_lower,:],1); v; flipdim(v[end-pad_upper+1:end,:],1)]
        return jo_convert(rdt,w,false)
    end
    function mirror_tran(v::LocalVector{VT},n::Integer,pad_upper::Integer,pad_lower::Integer,rdt=DataType) where {VT<:Number}
        w = v[pad_lower+1:end-pad_upper]
        w[1:pad_lower] += reverse(v[1:pad_lower])
        w[end-pad_upper+1:end] += reverse(v[end-pad_upper+1:end])
        return jo_convert(rdt,w,false)
    end
    function mirror_tran(v::LocalMatrix{VT},n::Integer,pad_upper::Integer,pad_lower::Integer,rdt=DataType) where {VT<:Number}
        w = v[pad_lower+1:end-pad_upper,:]
        w[1:pad_lower,:] += flipdim(v[1:pad_lower,:],1)
        w[end-pad_upper+1:end,:] += flipdim(v[end-pad_upper+1:end,:],1)
        return jo_convert(rdt,w,false)
    end
    ### periodic
    function periodic_fwd(v::LocalVector{VT},n::Integer,pad_upper::Integer,pad_lower::Integer,rdt=DataType) where {VT<:Number}
        w = [v[end-pad_lower+1:end]; v; v[1:pad_upper]]
        return jo_convert(rdt,w,false)
    end
    function periodic_fwd(v::LocalMatrix{VT},n::Integer,pad_upper::Integer,pad_lower::Integer,rdt=DataType) where {VT<:Number}
        w = [v[end-pad_lower+1:end,:]; v; v[1:pad_upper,:]]
        return jo_convert(rdt,w,false)
    end
    function periodic_tran(v::LocalVector{VT},n::Integer,pad_upper::Integer,pad_lower::Integer,rdt=DataType) where {VT<:Number}
        w = v[pad_lower+1:end-pad_upper]
        w[1:pad_upper] += v[end-pad_upper+1:end]
        w[end-pad_lower+1:end] += v[1:pad_lower]
        return jo_convert(rdt,w,false)
    end
    function periodic_tran(v::LocalMatrix{VT},n::Integer,pad_upper::Integer,pad_lower::Integer,rdt=DataType) where {VT<:Number}
        w = v[pad_lower+1:end-pad_upper,:]
        w[1:pad_upper,:] += v[end-pad_upper+1:end,:]
        w[end-pad_lower+1:end,:] += v[1:pad_lower,:]
        return jo_convert(rdt,w,false)
    end
end
using .joExtend_etc

export joExtend
"""
    julia> op = joExtend(n,pad_type;[pad_lower=...,][pad_upper=...,][DDT=...,][RDT=...,][name=...])

1D extension operator

# Signature

    function joExtend(n::Integer,pad_type::Symbol;
        pad_upper::Integer=0,pad_lower::Integer=0,
        DDT::DataType=joFloat,RDT::DataType=DDT,
        name="joExtend")

# Arguments

- `n` : size of input vector
- `pad_type` : one of the symbols
   - `:zeros` : pad signal with zeros
   - `:border` : pad signal with values at the edge of the domain
   - `:mirror` : mirror extension of the signal
   - `:periodic` : periodic extension of the signal
- keywords
    - `pad_lower` : number of points to pad on the lower index end
    - `pad_upper` : number of points to pad on the upper index end
    - `DDT`: domain data type
    - `RDT`: range data type
    - `name`: custom name

# Examples

extend a n-length vector with 10 zeros on either side

    joExtend(n,:zeros; pad_lower=10,pad_upper=10)

append, to a n-length vector, so that x[n+1:n+10] = x[n]

    joExtend(n,:border; pad_upper=10)

prepend, to n-length vector, its mirror extension: y=[reverse(x[1:10]);x]

    joExtend(n,:mirror; pad_lower=10)

append, to n-length vector, its periodic extension: y=[x;x[1:10]]

    joExtend(n,:periodic; pad_upper=10)

examples with DDT/RDT

    joExtend(n,:mirror; pad_lower=10,DDT=Float32)
    joExtend(n,:periodic; pad_upper=10,DDT=Float32,RDT=Float64)

"""
function joExtend(n::Integer,pad_type::Symbol;
    pad_upper::Integer=0,pad_lower::Integer=0,
    DDT::DataType=joFloat,RDT::DataType=DDT,
    name="joExtend")

    (pad_upper>=0) || throw(joLinearFunctionException("joExtend: invalid pad_upper size; should be >=0"))
    (pad_lower>=0) || throw(joLinearFunctionException("joExtend: invalid pad_lower size; should be >=0"))
    if pad_type==:zeros
        return joLinearFunctionFwd(n+pad_lower+pad_upper,n,
                                v1->joExtend_etc.zeros_fwd(v1,n,pad_upper,pad_lower,RDT),
                                v2->joExtend_etc.zeros_tran(v2,n,pad_upper,pad_lower,DDT),
                                v2->joExtend_etc.zeros_tran(v2,n,pad_upper,pad_lower,DDT),
                                v1->joExtend_etc.zeros_fwd(v1,n,pad_upper,pad_lower,RDT),
                                DDT,RDT;fMVok=true,
                                name=name*"_zeros")
    elseif pad_type==:border
        return joLinearFunctionFwd(n+pad_lower+pad_upper,n,
                                v1->joExtend_etc.border_fwd(v1,n,pad_upper,pad_lower,RDT),
                                v2->joExtend_etc.border_tran(v2,n,pad_upper,pad_lower,DDT),
                                v2->joExtend_etc.border_tran(v2,n,pad_upper,pad_lower,DDT),
                                v1->joExtend_etc.border_fwd(v1,n,pad_upper,pad_lower,RDT),
                                DDT,RDT;fMVok=true,
                                name=name*"_border")
    elseif pad_type==:mirror
        (pad_upper<=n) || throw(joLinearFunctionException("joExtend: invalid pad_upper size; should be <=n"))
        (pad_lower<=n) || throw(joLinearFunctionException("joExtend: invalid pad_lower size; should be <=n"))
        return joLinearFunctionFwd(n+pad_lower+pad_upper,n,
                                v1->joExtend_etc.mirror_fwd(v1,n,pad_upper,pad_lower,RDT),
                                v2->joExtend_etc.mirror_tran(v2,n,pad_upper,pad_lower,DDT),
                                v2->joExtend_etc.mirror_tran(v2,n,pad_upper,pad_lower,DDT),
                                v1->joExtend_etc.mirror_fwd(v1,n,pad_upper,pad_lower,RDT),
                                DDT,RDT;fMVok=true,
                                name=name*"_mirror")
    elseif pad_type==:periodic
        (pad_upper<=n) || throw(joLinearFunctionException("joExtend: invalid pad_upper size; should be <=n"))
        (pad_lower<=n) || throw(joLinearFunctionException("joExtend: invalid pad_lower size; should be <=n"))
        return joLinearFunctionFwd(n+pad_lower+pad_upper,n,
                                v1->joExtend_etc.periodic_fwd(v1,n,pad_upper,pad_lower,RDT),
                                v2->joExtend_etc.periodic_tran(v2,n,pad_upper,pad_lower,DDT),
                                v2->joExtend_etc.periodic_tran(v2,n,pad_upper,pad_lower,DDT),
                                v1->joExtend_etc.periodic_fwd(v1,n,pad_upper,pad_lower,RDT),
                                DDT,RDT;fMVok=true,
                                name=name*"_periodic")
    else
        throw(joLinearFunctionException("joExtend: unknown extension type."))
    end
    return nothing
end