text:
using PyCall: pyimport_conda, pycall
using Conda

function installmac()
   try
      pyimport_conda("nomkl", "nomkl")
      pyimport_conda("sklearn", "scikit-learn")
      println("mkl and scikit-learn successfully installed")
   catch
      try
         # https://github.com/JuliaPy/Conda.jl/issues/182
         # https://docs.anaconda.com/mkl-optimizations/
         Conda.add("nomkl")
         Conda.add("scikit-learn")
         Conda.rm("mkl")
         pyimport_conda("sklearn", "scikit-learn")
         pyimport_conda("sklearn.decomposition", "scikit-learn")
         pyimport_conda("mkl", "mkl")
         println("mkl and scikit-learn successfully installed")
      catch
         println("mkl and scikit-learn failed to install")
      end
   end
end

function installwinlinux()
   try
      pyimport_conda("sklearn", "scikit-learn")
      println("scikit-learn successfully installed")
   catch
      try
         # https://github.com/JuliaPy/Conda.jl/issues/182
         # https://docs.anaconda.com/mkl-optimizations/
         Conda.add("scikit-learn")
         pyimport_conda("sklearn", "scikit-learn")
         pyimport_conda("sklearn.decomposition", "scikit-learn")
         println("scikit-learn successfully installed")
      catch
         println("scikit-learn failed to install")
      end
   end
end

if Sys.isapple()
   installmac()
else
   installwinlinux()
end


text:
using Combinatorics
using DataStructures

struct PriorTracking{C,Sp,Fp} <: Tracking
  cohere::C
  tc::typeof(1.0s)
  source_prior::Sp
  freq_prior::Fp
  max_sources::Int
  min_norm::Float64
end
function checknorm(normalize,min_norm)
  if min_norm<Inf && normalize==false
    error("You must set normalize to true for `min_norm` value to be",
          " non-infinite")
  elseif min_norm == Inf && normalize==true
    error("You must specify a non-infinite value for `min_norm` if",
          " `normalize==true`.")
  end
end

Δt(x::PriorTracking) = Δt(x.cohere)
function Tracking(C,::Val{:prior};tc=1s,source_prior=nothing,
                  freq_prior=nothing, max_sources=4,
                  normalize=false,min_norm=Inf)
  checknorm(normalize,min_norm)
  PriorTracking(getmeta(C),tc,source_prior,freq_prior,max_sources,
                min_norm)
end

include("tracking_priors.jl")

# the number of ways observations A can be mapped to a set of source in B when
# it is possible that source b ∈ B is modeled as sum of elements in A
possible_groupings(n_sources,n_obs) =
  (Grouping(grouping,mapping)
   for grouping in partitions(1:n_obs)
   for sources in combinations(1:n_sources,length(grouping))
   for mapping in permutations(sources))

struct PermutedCoherence{C,P}
  cohere::C
  permuted::P
end

function prepare_coherence(C::Coherence,min_norm)
  @assert axisdim(C,Axis{:time}) == 1
  @assert axisdim(C,Axis{:scale}) == 2
  @assert axisdim(C,Axis{:freq}) == 3
  @assert axisdim(C,Axis{:component}) == 4

  # reorder for cache efficiency and then and normalize the data if necessary
  C_ = permutedims(C,[2,3,4,1])
  if min_norm < Inf
    for I in CartesianIndices((Base.axes(C_,3),Base.axes(C_,4)))
      C_[:,:,I] ./= max(min_norm,norm(vec(view(C_,:,:,I))))
    end
  end
  PermutedCoherence(C,C_)
end

function track(C::Coherence,params,args...)
  track(prepare_coherence(C,params.min_norm),params,args...)
end
function track(perm::PermutedCoherence,params::PriorTracking,progressbar=true,
               progress = track_progress(progressbar,ntimes(C),"prior"))
  C = perm.cohere
  C_ = perm.permuted
  timeax = 4

  track = TrackedSources(size(C_,1)*size(C_,2),params)

  C_out = zeros(eltype(C_),size(C_,1),size(C_,2)...,
                params.max_sources,size(C_,4))
  source_out = copy(C_out)
  sourceS_out = copy(C_out)

  window = ceil(Int,params.tc / Δt(C))
  lp_out = AxisArray(fill(0.0,ntimes(C)),AxisArrays.axes(C,1))
  buf = Array{eltype(C_)}(undef,size(C_,1),size(C_,2))
  function sumcomponents(C_t,kk)
    y = sum!(buf,view(C_t,:,:,kk,:))
    # buf ./= size(C_t,timeax)
    vec(buf)
  end

  veccomponent(C_t,k) = vec(view(C_t,:,:,k))

  groupings = CircularDeque{Grouping}(window+1)
  for t in eachindex(times(C))
    # find the MAP grouping of sources
    MAPgrouping, lp_out[t] =
      maximumby(g -> logpdf(track,view(C_,:,:,:,t:t),sumcomponents,g),
                possible_groupings(params.max_sources,ncomponents(C)))

    # arrange the sources
    for (kk,i) in iterable(MAPgrouping)
      for k in kk
        C_out[:,:,i,t] .+= C_[:,:,k,t]
      end
    end

    # update the source models
    update!(track,view(C_out,:,:,:,t),veccomponent,MAPgrouping)
    # mult!(track,decay)
    push!(groupings,MAPgrouping)
    if t > window
      downdate!(track,view(C_out,:,:,:,t-window),veccomponent,popfirst!(groupings))
    end

    next!(progress)
  end

  tracking =
    MetaArray(params,AxisArray(C_out,AxisArrays.axes(C,2),
                                     AxisArrays.axes(C,3),
                                     Axis{:component}(1:params.max_sources),
                                     AxisArrays.axes(C,1)))
  (tracking,lp_out)
end

function sort_components(x)
  order = sortperm(component_means(x),rev=true)
  x .= x[Axis{:component}(order)]
  x
end

ext: jl

max_stars_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_issues_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_forks_repo_name: JuliaTagBot/AuditoryBistabilityLE


text:
# Get the parameter values for gridsearch

function getgridparamset(learnertype, setting; Npar = 20)
    # setting[1] is the learnerparam that we vary
    if ((learnertype == TLeakyIntegrator)
        || (learnertype == GLeakyIntegrator))
        params = (setting[1],)
    elseif ((learnertype == TVarSmile)
        || (learnertype == TSmileOriginal)
        || (learnertype == GSmileOriginal)
        || (learnertype == GVarSmile))
        params = (setting[1], setting[2])#
    elseif ((learnertype == TParticleFilter)
        || (learnertype == GParticleFilter)
        || (learnertype == TMPN)
        || (learnertype == GMPN)
        || (learnertype == TSOR)
        || (learnertype == GSOR))
        params = (Npar, setting[2], setting[1])
    elseif ((learnertype == TBayesFilter)
        || (learnertype == GBayesFilter)
        || (learnertype == GNas10)
        || (learnertype == GNas12)
        || (learnertype == GNas10original)
        || (learnertype == GNas12original))
        params = (setting[2], setting[1])
    else
        println("Not sure what to do")
    end
    params
end

# If the inserted parameters are 0.,
# then load the parameter values from the provided parameterDict (optimized ones)
# within the 'setting'
function getparamset(setting)
    # Leaky: If 0s: set from Dictionary
    # VarSmile: If 0s: set m from Dictionary at position 1 and stochasticity at position 2
    # Smile: If 0s: set m from Dictionary at position 1 and stochasticity at position 2
    # PF: If 0s: set stochasticity at position 2 and changeprobability at position 3
    # Bayefilter: set stochasticity at position 1, and changeprobability at position 2
    # MPN: set stochasticity at position 1, and changeprobability at position 2
    # SOR: set stochasticity at position 1, and changeprobability at position 2
    isparamzerolist = isequal.(setting[1][1].param, 0)
    indexparam = findall(isparamzerolist)

    paramset = deepcopy(setting[1][1].param)
    if in(1, indexparam)# leak OR m OR s OR nparticles
        if ((setting[1][1].type == TLeakyIntegrator)
            || (setting[1][1].type == GLeakyIntegrator)
            || (setting[1][1].type == TVarSmile)
            || (setting[1][1].type == TSmileOriginal)
            || (setting[1][1].type == GSmileOriginal)
            || (setting[1][1].type == GVarSmile))
            if !in(:parameterDict, fieldnames(typeof(setting[1][1])))
                error("You need to insert a parameter value or a parameter dictionary")
            else
                paramvalue = setting[1][1].parameterDict[(setting[1][2], setting[1][3])]
                paramset = parameterset(paramset, paramvalue, 1)
            end
        elseif ((setting[1][1].type == TBayesFilter)
            || (setting[1][1].type == GBayesFilter)
            || (setting[1][1].type == GNas10)
            || (setting[1][1].type == GNas12)
            || (setting[1][1].type == GNas10original)
            || (setting[1][1].type == GNas12original))
            paramset = parameterset(paramset, setting[1][2], 1)
        elseif ((setting[1][1].type == TParticleFilter)
            || (setting[1][1].type == GParticleFilter)
            || (setting[1][1].type == TMPN)
            || (setting[1][1].type == GMPN)
            || (setting[1][1].type == TSOR)
            || (setting[1][1].type == GSOR))
            error("You need to insert the nr of particles")
        end
    end
    if in(2, indexparam) # pc and or
        if ((setting[1][1].type == TParticleFilter)
            || (setting[1][1].type == TVarSmile)
            || (setting[1][1].type == TSmileOriginal)
            || (setting[1][1].type == GParticleFilter)
            || (setting[1][1].type == GSmileOriginal)
            || (setting[1][1].type == GVarSmile)
            || (setting[1][1].type == TMPN)
            || (setting[1][1].type == GMPN)
            || (setting[1][1].type == TSOR)
            || (setting[1][1].type == GSOR))
            paramset = parameterset(paramset, setting[1][2], 2)
        elseif ((setting[1][1].type == GNas10)
            || (setting[1][1].type == GNas12)
            || (setting[1][1].type == GNas10original)
            || (setting[1][1].type == GNas12original)
            )
            if !in(:parameterDict, fieldnames(typeof(setting[1][1]))) #ismissing(parameterDict)
                error("You need to insert a parameter value or a parameter dictionary")
            else
                paramvalue = setting[1][1].parameterDict[(setting[1][2], setting[1][3])]
                paramset = parameterset(paramset, paramvalue, 2)
            end
        elseif ((setting[1][1].type == TBayesFilter)
            || (setting[1][1].type == GBayesFilter)
            # || (setting[1][1].type == GNas10) # Now we optimize it, so it should be imported from the parameter dictionary
            # || (setting[1][1].type == GNas12) # Now we optimize it, so it should be imported from the parameter dictionary
            # || (setting[1][1].type == GNas10original) # Now we optimize it, so it should be imported from the parameter dictionary
            # || (setting[1][1].type == GNas12original) # Now we optimize it, so it should be imported from the parameter dictionary
            )
            paramset = parameterset(paramset, setting[1][3], 2)
        end
    end
    if in(3, indexparam) # pc
        paramset = parameterset(paramset, setting[1][3], 3) # Pass changeprobability as 3rd parameter, also for other algorithms... (eg VarSmile for Shannon calculation)
    end
    paramset
end


text:
using Images, GeometryTypes, GLVisualize, Reactive, GLWindow, Colors
using GLAbstraction, GLFW
import GLAbstraction: singlepressed, imagespace
import GLVisualize: moving_average

if !isdefined(:runtests)
    window = glscreen()
end

description = """
Drawing lines and then adding symmetrie to them.
"""

window.color = RGBA(0f0, 0f0, 0f0, 1f0)

const linebuffer = Signal(fill(Point2f0(NaN), 4))
const middle = Point2f0(widths(window)) / 2
const symmetries = 18
const unsmoothed = copy(value(linebuffer))

function angle_between(a, b)
    a = normalize(a)
    b = normalize(b)
    atan2(cross(a, b), dot(a, b))
end

symmetry_lines = foldp(Point2f0[], linebuffer) do v0, lines
    resize!(v0, symmetries * length(lines))
    i = 1;
    for fi = linspace(0, 2pi, symmetries)
        for point in lines
            diff = middle - point
            radius = norm(diff)
            _fi = (fi + angle_between(diff, Point2f0(0, 1))) % 2pi
            v0[i] = radius * Point2f0(
                sin(_fi + pi),
                cos(_fi + pi)
            ) + middle
            i += 1
        end
    end
    v0
end

const lineobj = visualize(
    symmetry_lines, :lines,
    color = RGBA(0.8f0, 0.8f0, 0.8f0, 1f0),
    thickness = 1.5f0,
    boundingbox = AABB{Float32}(value(window.area)) # boundingbox for center!
)
_view(lineobj, window, camera = :orthographic_pixel)

@materialize mouseposition, mouse_buttons_pressed, mouseinside, buttons_pressed = window.inputs

camera = window.cameras[:orthographic_pixel]
const history = Point2f0[]
s = map(mouseposition, mouse_buttons_pressed, init = nothing) do mp, mbp
    l0 = value(linebuffer)
    if singlepressed(mbp, GLFW.MOUSE_BUTTON_LEFT) && value(mouseinside) && isempty(value(buttons_pressed))
        p = imagespace(mp, camera)
        keep, p = moving_average(p, 1.5f0, history)
        if keep
            push!(linebuffer, push!(l0, p))
        end
    else
        if !isnan(last(l0)) # only push one NaN to seperate
            empty!(history) # reset
            push!(linebuffer, push!(l0, Point2f0(NaN)))
        end
    end
    nothing
end
# preserve signals, so that it doesn't get garbage collected.
preserve(s)
# we need to define init, because otherwise it will be initialised by calling
# the function one time, which

if !isdefined(:runtests)
    renderloop(window)
end


text:
using ShipMMG
using Test
using Distributions

@testset "ShipMMG.jl" begin

    include("test_data.jl")
    include("test_kt.jl")
    include("test_mmg.jl")
    include("test_draw.jl")

end

max_stars_repo_name: mitsuyukiLab/JuliaShipMMG

max_issues_repo_name: mitsuyukiLab/JuliaShipMMG


text:
import Random

export CubaAlg

struct CubaAlg{A,K}
    algorithm::A
    kw::K
end

function CubaAlg(a; kw...)
    CubaAlg(a, kw)
end

struct CubaFun{T,F,D<:Domain}
    inner::F
    domain::D
    ndim::Int
    ncomp::Int
end

function CubaFun(f, dom::Domain)
    inner = f
    ndim   = ndims(dom)
    rng = Random.GLOBAL_RNG
    x = draw(rng, dom).position
    y = f(x)
    ncomp  = length(y)
    T = typeof(float.(y))
    F = typeof(f)
    D = typeof(dom)
    CubaFun{T,F,D}(inner, dom, ndim, ncomp)
end

function (f::CubaFun)(pt0, out)
    dom = f.domain
    vol = volume(dom)
    lo = dom.lower
    up = dom.upper
    pt = lo .+ pt0 .* (up .- lo)
    y = vol * f.inner(pt)
    copyto!(out, y)
end

function from_vector(::Type{T}, v) where {T <: Number}
    @assert length(v) == 1
    T(first(v))
end

function from_vector(::Type{T}, v) where {T}
    T(v)
end

canonicalize(f, dom::NTuple{2,Number}, alg::CubaAlg) = canonicalize_cuba(f∘first, dom, alg)
canonicalize(f, dom, alg::CubaAlg) = canonicalize_cuba(f, dom, alg)

function canonicalize_cuba(f, dom, alg::CubaAlg)
    dom2 = Domain(dom)
    CubaFun(f, dom2), dom2, alg
end

@noinline function integral_kernel(f::CubaFun{T}, dom::Domain, cuba::CubaAlg) where {T}
    raw = cuba.algorithm(f, f.ndim, f.ncomp; cuba.kw...)

    (value=from_vector(T, raw.integral),
        std=from_vector(T, raw.error),
        neval=raw.neval, raw=raw)
end


text:
"""
    SparseKKTSystem{T, MT} <: AbstractReducedKKTSystem{T, MT}

Implement the [`AbstractReducedKKTSystem`](@ref) in sparse COO format.

"""
struct SparseKKTSystem{T, MT} <: AbstractReducedKKTSystem{T, MT}
    hess::StrideOneVector{T}
    jac::StrideOneVector{T}
    pr_diag::StrideOneVector{T}
    du_diag::StrideOneVector{T}
    # Augmented system
    aug_raw::SparseMatrixCOO{T,Int32}
    aug_com::MT
    aug_csc_map::Union{Nothing, Vector{Int}}
    # Jacobian
    jac_raw::SparseMatrixCOO{T,Int32}
    jac_com::MT
    jac_csc_map::Union{Nothing, Vector{Int}}
    # Info
    ind_ineq::Vector{Int}
    ind_fixed::Vector{Int}
    ind_aug_fixed::Vector{Int}
    jacobian_scaling::Vector{T}
end

"""
    SparseUnreducedKKTSystem{T, MT} <: AbstractUnreducedKKTSystem{T, MT}

Implement the [`AbstractUnreducedKKTSystem`](@ref) in sparse COO format.

"""
struct SparseUnreducedKKTSystem{T, MT} <: AbstractUnreducedKKTSystem{T, MT}
    hess::StrideOneVector{T}
    jac::StrideOneVector{T}
    pr_diag::StrideOneVector{T}
    du_diag::StrideOneVector{T}

    l_diag::StrideOneVector{T}
    u_diag::StrideOneVector{T}
    l_lower::StrideOneVector{T}
    u_lower::StrideOneVector{T}

    aug_raw::SparseMatrixCOO{T,Int32}
    aug_com::MT
    aug_csc_map::Union{Nothing, Vector{Int}}

    jac_raw::SparseMatrixCOO{T,Int32}
    jac_com::MT
    jac_csc_map::Union{Nothing, Vector{Int}}
    ind_ineq::Vector{Int}
    ind_fixed::Vector{Int}
    ind_aug_fixed::Vector{Int}
    jacobian_scaling::Vector{T}
end

# Template to dispatch on sparse representation
const AbstractSparseKKTSystem{T, MT} = Union{SparseKKTSystem{T, MT}, SparseUnreducedKKTSystem{T, MT}}

#=
    Generic sparse methods
=#

function mul!(y::AbstractVector, kkt::AbstractSparseKKTSystem, x::AbstractVector)
    mul!(y, Symmetric(kkt.aug_com, :L), x)
end
function mul!(y::AbstractKKTVector, kkt::AbstractSparseKKTSystem, x::AbstractKKTVector)
    mul!(full(y), Symmetric(kkt.aug_com, :L), full(x))
end

function jtprod!(y::AbstractVector, kkt::AbstractSparseKKTSystem, x::AbstractVector)
    mul!(y, kkt.jac_com', x)
end

nnz_jacobian(kkt::AbstractSparseKKTSystem) = nnz(kkt.jac_raw)

function compress_jacobian!(kkt::AbstractSparseKKTSystem{T, MT}) where {T, MT<:SparseMatrixCSC{T, Int32}}
    ns = length(kkt.ind_ineq)
    kkt.jac[end-ns+1:end] .= -1.0
    kkt.jac .*= kkt.jacobian_scaling # scaling
    transfer!(kkt.jac_com, kkt.jac_raw, kkt.jac_csc_map)
end

function compress_jacobian!(kkt::AbstractSparseKKTSystem{T, MT}) where {T, MT<:Matrix{T}}
    ns = length(kkt.ind_ineq)
    kkt.jac[end-ns+1:end] .= -1.0
    kkt.jac .*= kkt.jacobian_scaling # scaling
    copyto!(kkt.jac_com, kkt.jac_raw)
end

function set_jacobian_scaling!(kkt::AbstractSparseKKTSystem{T, MT}, constraint_scaling::AbstractVector) where {T, MT}
    nnzJ = length(kkt.jac)::Int
    @inbounds for i in 1:nnzJ
        index = kkt.jac_raw.I[i]
        kkt.jacobian_scaling[i] = constraint_scaling[index]
    end
end

#=
    SparseKKTSystem
=#

function SparseKKTSystem{T, MT}(
    n::Int, m::Int, ind_ineq::Vector{Int}, ind_fixed::Vector{Int},
    hess_sparsity_I, hess_sparsity_J,
    jac_sparsity_I, jac_sparsity_J,
) where {T, MT}
    n_jac = length(jac_sparsity_I)
    n_hess = length(hess_sparsity_I)

    aug_vec_length = n+m
    aug_mat_length = n+m+n_hess+n_jac

    I = Vector{Int32}(undef, aug_mat_length)
    J = Vector{Int32}(undef, aug_mat_length)
    V = Vector{T}(undef, aug_mat_length)
    fill!(V, 0.0)  # Need to initiate V to avoid NaN

    offset = n+n_jac+n_hess+m

    I[1:n] .= 1:n
    I[n+1:n+n_hess] = hess_sparsity_I
    I[n+n_hess+1:n+n_hess+n_jac] .= (jac_sparsity_I.+n)
    I[n+n_hess+n_jac+1:offset] .= (n+1:n+m)

    J[1:n] .= 1:n
    J[n+1:n+n_hess] = hess_sparsity_J
    J[n+n_hess+1:n+n_hess+n_jac] .= jac_sparsity_J
    J[n+n_hess+n_jac+1:offset] .= (n+1:n+m)

    pr_diag = view(V, 1:n)
    du_diag = view(V, n_jac+n_hess+n+1:n_jac+n_hess+n+m)

    hess = view(V, n+1:n+n_hess)
    jac = view(V, n_hess+n+1:n_hess+n+n_jac)

    aug_raw = SparseMatrixCOO(aug_vec_length,aug_vec_length,I,J,V)
    jac_raw = SparseMatrixCOO(m,n,jac_sparsity_I,jac_sparsity_J,jac)

    aug_com = MT(aug_raw)
    jac_com = MT(jac_raw)

    aug_csc_map = get_mapping(aug_com, aug_raw)
    jac_csc_map = get_mapping(jac_com, jac_raw)

    ind_aug_fixed = if isa(aug_com, SparseMatrixCSC)
        _get_fixed_variable_index(aug_com, ind_fixed)
    else
        zeros(Int, 0)
    end
    jac_scaling = ones(T, n_jac)

    return SparseKKTSystem{T, MT}(
        hess, jac, pr_diag, du_diag,
        aug_raw, aug_com, aug_csc_map,
        jac_raw, jac_com, jac_csc_map,
        ind_ineq, ind_fixed, ind_aug_fixed, jac_scaling,
    )
end

# Build KKT system directly from AbstractNLPModel
function SparseKKTSystem{T, MT}(nlp::AbstractNLPModel, ind_cons=get_index_constraints(nlp)) where {T, MT}
    n_slack = length(ind_cons.ind_ineq)
    # Deduce KKT size.
    n = get_nvar(nlp) + n_slack
    m = get_ncon(nlp)
    # Evaluate sparsity pattern
    jac_I = Vector{Int32}(undef, get_nnzj(nlp))
    jac_J = Vector{Int32}(undef, get_nnzj(nlp))
    jac_structure!(nlp,jac_I, jac_J)

    hess_I = Vector{Int32}(undef, get_nnzh(nlp))
    hess_J = Vector{Int32}(undef, get_nnzh(nlp))
    hess_structure!(nlp,hess_I,hess_J)

    force_lower_triangular!(hess_I,hess_J)
    # Incorporate slack's sparsity pattern
    append!(jac_I, ind_cons.ind_ineq)
    append!(jac_J, get_nvar(nlp)+1:get_nvar(nlp)+n_slack)

    return SparseKKTSystem{T, MT}(
        n, m, ind_cons.ind_ineq, ind_cons.ind_fixed,
        hess_I, hess_J, jac_I, jac_J,
    )
end

is_reduced(::SparseKKTSystem) = true
num_variables(kkt::SparseKKTSystem) = length(kkt.pr_diag)

#=
    SparseUnreducedKKTSystem
=#

function SparseUnreducedKKTSystem{T, MT}(
    n::Int, m::Int, nlb::Int, nub::Int, ind_ineq, ind_fixed,
    hess_sparsity_I, hess_sparsity_J,
    jac_sparsity_I, jac_sparsity_J,
    ind_lb, ind_ub,
) where {T, MT}
    n_jac = length(jac_sparsity_I)
    n_hess = length(hess_sparsity_I)

    aug_mat_length = n + m + n_hess + n_jac + 2*nlb + 2*nub
    aug_vec_length = n+m+nlb+nub

    I = Vector{Int32}(undef, aug_mat_length)
    J = Vector{Int32}(undef, aug_mat_length)
    V = zeros(aug_mat_length)

    offset = n+n_jac+n_hess+m

    I[1:n] .= 1:n
    I[n+1:n+n_hess] = hess_sparsity_I
    I[n+n_hess+1:n+n_hess+n_jac].=(jac_sparsity_I.+n)
    I[n+n_hess+n_jac+1:offset].=(n+1:n+m)

    J[1:n] .= 1:n
    J[n+1:n+n_hess] = hess_sparsity_J
    J[n+n_hess+1:n+n_hess+n_jac].=jac_sparsity_J
    J[n+n_hess+n_jac+1:offset].=(n+1:n+m)

    I[offset+1:offset+nlb] .= (1:nlb).+(n+m)
    I[offset+nlb+1:offset+2nlb] .= (1:nlb).+(n+m)
    I[offset+2nlb+1:offset+2nlb+nub] .= (1:nub).+(n+m+nlb)
    I[offset+2nlb+nub+1:offset+2nlb+2nub] .= (1:nub).+(n+m+nlb)
    J[offset+1:offset+nlb] .= (1:nlb).+(n+m)
    J[offset+nlb+1:offset+2nlb] .= ind_lb
    J[offset+2nlb+1:offset+2nlb+nub] .= (1:nub).+(n+m+nlb)
    J[offset+2nlb+nub+1:offset+2nlb+2nub] .= ind_ub

    pr_diag = view(V,1:n)
    du_diag = view(V,n_jac+n_hess+n+1:n_jac+n_hess+n+m)

    l_diag = view(V,offset+1:offset+nlb)
    l_lower= view(V,offset+nlb+1:offset+2nlb)
    u_diag = view(V,offset+2nlb+1:offset+2nlb+nub)
    u_lower= view(V,offset+2nlb+nub+1:offset+2nlb+2nub)

    hess = view(V, n+1:n+n_hess)
    jac = view(V, n_hess+n+1:n_hess+n+n_jac)

    aug_raw = SparseMatrixCOO(aug_vec_length,aug_vec_length,I,J,V)
    jac_raw = SparseMatrixCOO(m,n,jac_sparsity_I,jac_sparsity_J,jac)

    aug_com = MT(aug_raw)
    jac_com = MT(jac_raw)

    aug_csc_map = get_mapping(aug_com, aug_raw)
    jac_csc_map = get_mapping(jac_com, jac_raw)

    jac_scaling = ones(T, n_jac)

    ind_aug_fixed = if isa(aug_com, SparseMatrixCSC)
        _get_fixed_variable_index(aug_com, ind_fixed)
    else
        zeros(Int, 0)
    end

    return SparseUnreducedKKTSystem{T, MT}(
        hess, jac, pr_diag, du_diag,
        l_diag, u_diag, l_lower, u_lower,
        aug_raw, aug_com, aug_csc_map,
        jac_raw, jac_com, jac_csc_map,
        ind_ineq, ind_fixed, ind_aug_fixed, jac_scaling,
    )
end

function SparseUnreducedKKTSystem{T, MT}(nlp::AbstractNLPModel, ind_cons=get_index_constraints(nlp)) where {T, MT}
    n_slack = length(ind_cons.ind_ineq)
    nlb = length(ind_cons.ind_lb)
    nub = length(ind_cons.ind_ub)
    # Deduce KKT size.
    n = get_nvar(nlp) + n_slack
    m = get_ncon(nlp)
    # Evaluate sparsity pattern
    jac_I = Vector{Int32}(undef, get_nnzj(nlp))
    jac_J = Vector{Int32}(undef, get_nnzj(nlp))
    jac_structure!(nlp,jac_I, jac_J)

    hess_I = Vector{Int32}(undef, get_nnzh(nlp))
    hess_J = Vector{Int32}(undef, get_nnzh(nlp))
    hess_structure!(nlp,hess_I,hess_J)

    force_lower_triangular!(hess_I,hess_J)
    # Incorporate slack's sparsity pattern
    append!(jac_I, ind_cons.ind_ineq)
    append!(jac_J, get_nvar(nlp)+1:get_nvar(nlp)+n_slack)

    return SparseUnreducedKKTSystem{T, MT}(
        n, m, nlb, nub, ind_cons.ind_ineq, ind_cons.ind_fixed,
        hess_I, hess_J, jac_I, jac_J, ind_cons.ind_lb, ind_cons.ind_ub,
    )
end

function initialize!(kkt::SparseUnreducedKKTSystem)
    kkt.pr_diag.=1
    kkt.du_diag.=0
    kkt.hess.=0
    kkt.l_lower.=0
    kkt.u_lower.=0
    kkt.l_diag.=1
    kkt.u_diag.=1
end

is_reduced(::SparseUnreducedKKTSystem) = false
num_variables(kkt::SparseUnreducedKKTSystem) = length(kkt.pr_diag)


text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2020 Varnerlab
# Robert Frederick Smith School of Chemical and Biomolecular Engineering
# Cornell University, Ithaca NY 14850
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #
#
# ----------------------------------------------------------------------------------- #
# Function: DataDictionary
# Description: Holds simulation and model parameters as key => value pairs in a Julia Dict()
# Generated on: 2020-02-05T11:38:24.145
#
# Input arguments:
# time_start::Float64 => Simulation start time value (scalar)
# time_stop::Float64 => Simulation stop time value (scalar)
# time_step::Float64 => Simulation time step (scalar)
#
# Output arguments:
# data_dictionary::Dict{String,Any} => Dictionary holding model and simulation parameters as key => value pairs
# ----------------------------------------------------------------------------------- #
function build_data_dictionary(time_span::Tuple{Float64,Float64,Float64}, path_to_biophysical_constants_file::String = "./Default.json", host_type::Symbol = :bacteria)::Dict{String,Any}

	# load the biophysical_constants dictionary
	biophysical_constants_dictionary = build_biophysical_dictionary(path_to_biophysical_constants_file, host_type)

	# stoichiometric_matrix and dilution_matrix -
	stoichiometric_matrix = readdlm("./Network.dat")

	# number of states, and rates -
	(number_of_states,number_of_rates) = size(stoichiometric_matrix)

	# array of species types -
	species_symbol_type_array = [
		:gene	;	# 1	deGFP
		:gene	;	# 2	sigma_70
		:mrna	;	# 3	mRNA_deGFP
		:mrna	;	# 4	mRNA_sigma_70
		:protein	;	# 5	protein_deGFP
		:protein	;	# 6	protein_sigma_70
	]

	# we need to store the species symbol array for later -
	biophysical_constants_dictionary["species_symbol_type_array"] = species_symbol_type_array

	# array of gene lengths -
	gene_coding_length_array = [
		1000.0	;	# 1	deGFP
		1000.0	;	# 2	sigma_70
	]

	# array of mRNA coding lengths -
	mRNA_coding_length_array = [
		gene_coding_length_array[1]	;	# 3	1	mRNA_deGFP
		gene_coding_length_array[2]	;	# 4	2	mRNA_sigma_70
	]

	# array of mRNA coding lengths -
	protein_coding_length_array = [
		round((0.33)*mRNA_coding_length_array[1])	;	# 5	1	protein_deGFP
		round((0.33)*mRNA_coding_length_array[2])	;	# 6	2	protein_sigma_70
	]

	# array of gene concentrations -
	gene_abundance_array = [
		5.0	;	# (nM) 1	deGFP
		5.0	;	# (nM) 2	sigma_70
	]

	# initial condition array -
	initial_condition_array = [
		gene_abundance_array[1]	;	# 1	deGFP
		gene_abundance_array[2]	;	# 2	sigma_70
		0.0	;	# 3	mRNA_deGFP
		0.0	;	# 4	mRNA_sigma_70
		0.0	;	# 5	protein_deGFP
		0.0	;	# 6	protein_sigma_70

		# translation capacity -
		100.0	;	# translation capacity
	]

	binding_parameter_dictionary = Dict{String,Float64}()
	binding_parameter_dictionary["n_deGFP_sigma_70"] = 1.0
	binding_parameter_dictionary["K_deGFP_sigma_70"] = 0.05

	# Alias the control function parameters -
	control_parameter_dictionary = Dict{String,Float64}()
	control_parameter_dictionary["W_deGFP_RNAP"] = 0.001
	control_parameter_dictionary["W_deGFP_sigma_70"] = 1.0
	control_parameter_dictionary["W_sigma_70_RNAP"] = 0.001

	# degradation modifiers -
	degradation_modifier_array = [
		0.0	;	# 1	deGFP
		0.0	;	# 2	sigma_70
		1.0	;	# 3	mRNA_deGFP
		1.0	;	# 4	mRNA_sigma_70
		1.0	;	# 5	protein_deGFP
		1.0	;	# 6	protein_sigma_70
	]

	# time constant modifiers -
	time_constant_modifier_array = [
		0.0	;	# 1	deGFP
		0.0	;	# 2	sigma_70
		1.0	;	# 3	mRNA_deGFP
		1.0	;	# 4	mRNA_sigma_70
		1.0	;	# 5	protein_deGFP
		1.0	;	# 6	protein_sigma_70
	]

	# Dilution degrdation matrix -
	dilution_degradation_matrix = build_dilution_degradation_matrix(biophysical_constants_dictionary,species_symbol_type_array,degradation_modifier_array)

	# Precompute the translation parameters -
	translation_parameter_array = precompute_translation_parameter_array(biophysical_constants_dictionary, protein_coding_length_array, time_constant_modifier_array, host_type)

	# Precompute the kinetic limit of transcription -
	transcription_kinetic_limit_array = precompute_transcription_kinetic_limit_array(biophysical_constants_dictionary, gene_coding_length_array, gene_abundance_array, time_constant_modifier_array, host_type)

	# Parameter name index array -
	parameter_name_mapping_array = [
		"n_deGFP_sigma_70"	;	# 1
		"K_deGFP_sigma_70"	;	# 2
		"W_deGFP_RNAP"	;	# 3
		"W_deGFP_sigma_70"	;	# 4
		"W_sigma_70_RNAP"	;	# 5
		"rnapII_concentration"	;	# 6
		"ribosome_concentration"	;	# 7
		"degradation_constant_mRNA"	;	# 8
		"degradation_constant_protein"	;	# 9
		"kcat_transcription"	;	# 10
		"kcat_translation"	;	# 11
		"maximum_specific_growth_rate"	;	# 12
		"saturation_constant_transcription"	;	# 13
		"saturation_constant_translation"	;	# 14
	]

	# =============================== DO NOT EDIT BELOW THIS LINE ============================== #
	data_dictionary = Dict{String,Any}()
	data_dictionary["number_of_states"] = number_of_states
	data_dictionary["species_symbol_type_array"] = species_symbol_type_array
	data_dictionary["initial_condition_array"] = initial_condition_array
	data_dictionary["gene_coding_length_array"] = gene_coding_length_array
	data_dictionary["mRNA_coding_length_array"] = mRNA_coding_length_array
	data_dictionary["protein_coding_length_array"] = protein_coding_length_array
	data_dictionary["stoichiometric_matrix"] = stoichiometric_matrix
	data_dictionary["dilution_degradation_matrix"] = dilution_degradation_matrix
	data_dictionary["binding_parameter_dictionary"] = binding_parameter_dictionary
	data_dictionary["control_parameter_dictionary"] = control_parameter_dictionary
	data_dictionary["parameter_name_mapping_array"] = parameter_name_mapping_array
	data_dictionary["transcription_kinetic_limit_array"] = transcription_kinetic_limit_array
	data_dictionary["translation_parameter_array"] = translation_parameter_array
	data_dictionary["degradation_modifier_array"] = degradation_modifier_array
	data_dictionary["time_constant_modifier_array"] = time_constant_modifier_array
	data_dictionary["biophysical_constants_dictionary"] = biophysical_constants_dictionary

	# extra stuff -
	data_dictionary["R"] = 8.314 			# J mol^-1 K^-1
	data_dictionary["T_K"] = 273.15 + 29.0 	# K
	data_dictionary["half_life_translation_capacity"] = 8.0	 # hr
	data_dictionary["base_line_weight"] = (1.0/1.0)
	# =============================== DO NOT EDIT ABOVE THIS LINE ============================== #
	return data_dictionary
end

function customize_data_dictionary(default_data_dictionary::Dict{String,Any}, host_type::Symbol)::Dict{String,Any}

	# make a deepcopy -
	customized_data_dictionary = deepcopy(default_data_dictionary)

	# set the gene lengths -
	gene_length_array = customized_data_dictionary["gene_coding_length_array"]
	gene_length_array[1] = 683.0    # deGFP
	gene_length_array[2] = 720.0    # S70 => doesn't matter, gene expression = 0
	customized_data_dictionary["gene_coding_length_array"] = gene_length_array

	# set the protein lengths -
	protein_coding_length_array = customized_data_dictionary["protein_coding_length_array"]
	protein_coding_length_array[1] = 239.0  # deGFP
	protein_coding_length_array[2] = 240.0  # S70 => doesn't matter, translation = 0
	customized_data_dictionary["protein_coding_length_array"] = protein_coding_length_array

	# set ICs -
	initial_condition_array = customized_data_dictionary["initial_condition_array"]
	initial_condition_array[1] = 0.005  # deGFP gene
	initial_condition_array[2] = 0.0    # no S70 gene, just S70 protein in the extract
	initial_condition_array[6] = 0.035 	# S70 protein in muM
	customized_data_dictionary["initial_condition_array"] = initial_condition_array

	# setup the W's -
	control_parameter_dictionary = customized_data_dictionary["control_parameter_dictionary"]
  	control_parameter_dictionary["W_deGFP_RNAP"] = 0.000014
  	control_parameter_dictionary["W_deGFP_sigma_70"] = 10.0
	control_parameter_dictionary["W_sigma_70_RNAP"] = 0.0
	customized_data_dictionary["control_parameter_dictionary"] = control_parameter_dictionary

	# setup binding parameters -
	binding_parameter_dictionary = customized_data_dictionary["binding_parameter_dictionary"]
	binding_parameter_dictionary["n_deGFP_sigma_70"] = 1.0
	binding_parameter_dictionary["K_deGFP_sigma_70"] = 0.05
	customized_data_dictionary["binding_parameter_dictionary"] = binding_parameter_dictionary

	# setup degradation_modifier_array -
	degradation_modifier_array = [
		0.0		;	# 1	deGFP
		0.0		;	# 2	sigma_70
		0.05	;	# 3	mRNA_deGFP
		0.05	;	# 4	mRNA_sigma_70
		1.0		;	# 5	protein_deGFP
		1.0		;	# 6	protein_sigma_70
	]
	customized_data_dictionary["degradation_modifier_array"] = degradation_modifier_array
	biophysical_constants_dictionary = customized_data_dictionary["biophysical_constants_dictionary"]
	species_symbol_type_array = customized_data_dictionary["species_symbol_type_array"]
	protein_coding_length_array = customized_data_dictionary["protein_coding_length_array"]
	gene_coding_length_array = customized_data_dictionary["gene_coding_length_array"]

	# get gene IC -
	idx_gene = findall(x->x==:gene,species_symbol_type_array)
	gene_abundance_array = initial_condition_array[idx_gene]

	# no sigma_70 gene -
	gene_abundance_array[2] = 0.0

	# Dilution degrdation matrix -
  	dilution_degradation_matrix = build_dilution_degradation_matrix(biophysical_constants_dictionary,species_symbol_type_array,degradation_modifier_array)
	customized_data_dictionary["dilution_degradation_matrix"] = dilution_degradation_matrix

	# time constant modifiers -
	time_constant_modifier_array = [
		0.0	;	# 1	deGFP
		0.0	;	# 2	sigma_70
		1.0	;	# 3	mRNA_deGFP
		1.0	;	# 4	mRNA_sigma_70
		3.0	;	# 5	protein_deGFP
		1.0	;	# 6	protein_sigma_70
	]
	customized_data_dictionary["time_constant_modifier_array"] = time_constant_modifier_array

	# Precompute the translation parameters -
  	translation_parameter_array = precompute_translation_parameter_array(biophysical_constants_dictionary, protein_coding_length_array, time_constant_modifier_array,host_type)
	customized_data_dictionary["translation_parameter_array"] = translation_parameter_array

  	# Precompute the kinetic limit of transcription -
  	transcription_kinetic_limit_array = precompute_transcription_kinetic_limit_array(biophysical_constants_dictionary, gene_coding_length_array, gene_abundance_array, time_constant_modifier_array, host_type)
	customized_data_dictionary["transcription_kinetic_limit_array"] = transcription_kinetic_limit_array

	# return -
	return customized_data_dictionary
  end

max_stars_repo_name: varnerlab/Biophysical-TXTL-Model-Code

max_issues_repo_name: varnerlab/Biophysical-TXTL-Model-Code

max_forks_repo_name: varnerlab/Biophysical-TXTL-Model-Code


text:
# Preliminary (and partial) support for conversion of images from Images.jl to OpenCV Mat
# Images in Images.jl are encoded as 2d arrays of RGB4{Ufixed8}, Float64, Ufixed8, or Ufixed12

# The following function converts both grayscale and color images to Mat from Image.jl arrays
#   TO_DO: colorspace is NOT correctly transformed
#   TO_DO:  improve algorithm(currently, takes 300 ms for a 512x512 image
#   Use * pointer indexing instead of the "at" method implemented in pixset

function convertToMat(image)
    img = Images.separate(image);
    cd = Images.colordim(img);
    if (typeof(img[1,1,1].i) == UInt8)
       if (cd < 3); mat = Mat(Base.size(img,2), Base.size(img,1), CV_8UC1); end
       if (cd == 3); mat = Mat(Base.size(img,2), Base.size(img,1), CV_8UC3); end
    elseif (typeof(img[1,1,1].i) == Float32)
       if (cd < 3); mat = Mat(Base.size(img,2), Base.size(img,1), CV_32FC1); end
       if (cd == 3); mat = Mat(Base.size(img,2), Base.size(img,1), CV_32FC3); end
    elseif (typeof(img[1,1,1].i) == Float64)
       if (cd < 3); mat = Mat(Base.size(img,2), Base.size(img,1), CV_64FC1); end
       if (cd == 3); mat = Mat(Base.size(img,2), Base.size(img,1), CV_64FC3); end
    else
       throw(ArgumentError("Pixel format not supported!"))
    end

    if (cd < 3)   # grayscale or binary
        for j = 1:Base.size(img,2)     # index row first (Mat is row-major order)
            for k =1:Base.size(img,1)  # index column second
                # slow algorithm  - will try to use pointer method (C++)!
                pixset(mat, k, j, float(img[k,j,1].i))
            end
        end
    end

   if (cd == 3)   # color (RGB) image
        for j = 1:Base.size(img,2)     # index row first (Mat is row-major order)
            for k =1:Base.size(img,1)  # index column second
                colorvec = tostdvec([float(img[k,j,1].i),float(img[k,j,2].i),float(img[k,j,3].i)])
                pixset(mat, k-1, j-1, colorvec)   # -1 to have 0-indexing per C++
            end
        end
    end

    return(mat)
end


text:
abstract type AbstractControlSystem <: AbstractEnv end

# BacksteppingControl_StaticAllocator_ControlSystem
struct BacksteppingControl_StaticAllocator_ControlSystem <: AbstractControlSystem
    controller::BacksteppingPositionController
    allocator::StaticAllocator
end

function State(cs::BacksteppingControl_StaticAllocator_ControlSystem)
    @unpack controller, allocator = cs
    return function (pos0, m, g)
        x_controller = State(controller)(pos0, m, g)
        ComponentArray(controller=x_controller)
    end
end

function Command(cs::BacksteppingControl_StaticAllocator_ControlSystem)
    @unpack controller, allocator = cs
    return function (p, v, R, ω, xd, vd, ad, ȧd, äd, Td, m, J, g, Λ̂)
        νd, Ṫd, _... = Command(controller)(p, v, R, ω, xd, vd, ad, ȧd, äd, Td, m, J, g)
        u_cmd = allocator(νd, Λ̂)
        ComponentArray(νd=νd, Ṫd=Ṫd, u_cmd=u_cmd)
    end
end

function Dynamics!(cs::BacksteppingControl_StaticAllocator_ControlSystem)
    @unpack controller = cs
    @Loggable function dynamics!(dx, x, p, t; pos_cmd=nothing, Ṫd)
        @nested_log :controller Dynamics!(controller)(dx.controller, x.controller, p, t; pos_cmd=pos_cmd, Ṫd=Ṫd)
    end
end

# BacksteppingControl_AdaptiveAllocator_ControlSystem
struct BacksteppingControl_AdaptiveAllocator_ControlSystem <: AbstractControlSystem
    controller::BacksteppingPositionController
    allocator::AdaptiveAllocator
end

function State(cs::BacksteppingControl_AdaptiveAllocator_ControlSystem)
    @unpack controller, allocator = cs
    return function (pos0, m, g; Θ̂=zeros(6, 4))
        x_controller = State(controller)(pos0, m, g)
        x_allocator = State(allocator)(Θ̂)
        ComponentArray(controller=x_controller, allocator=x_allocator)
    end
end

function Command(cs::BacksteppingControl_AdaptiveAllocator_ControlSystem)
    @unpack controller, allocator = cs
    @unpack Ap, Bp, P, Kp, Kt, Kω = controller
    return function (p, v, R, ω, xd, vd, ad, ȧd, äd, Td, m, J, g, Λ̂, Θ̂)
        νd, Ṫd, e, zB, T = Command(controller)(p, v, R, ω, xd, vd, ad, ȧd, äd, Td, m, J, g)
        Θ̂̇ = Command(allocator)(νd, e, zB, T, Θ̂, R, J,
                               Ap, Bp, P, Kp, Kt, Kω)
        u_cmd = allocator(νd, Θ̂)
        ComponentArray(νd=νd, Ṫd=Ṫd, u_cmd=u_cmd, Θ̂̇=Θ̂̇)
    end
end

function Dynamics!(cs::BacksteppingControl_AdaptiveAllocator_ControlSystem)
    @unpack controller, allocator = cs
    @Loggable function dynamics!(dx, x, p, t; pos_cmd=nothing, Ṫd, Θ̂̇)
        @nested_log :controller Dynamics!(controller)(dx.controller, x.controller, p, t; pos_cmd=pos_cmd, Ṫd=Ṫd)
        @nested_log :allocator Dynamics!(allocator)(dx.allocator, x.allocator, p, t; Θ̂̇=Θ̂̇)
    end
end


text:
#
# This file is part of the Actors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

"""
    @msg [Msg] A B C

Define empty structs as message types. `Msg` is an
existing abstract datatype.

To call `@msg Msg A B C` is equivalent to
```
struct A <: Msg end
struct B <: Msg end
struct C <: Msg end
```

To call `@msg D E F` is equivalent to
```
struct D end
struct E end
struct F end
```
"""
macro msg(syms...)
    isempty(syms) && throw(ArgumentError("no arguments given for msg"))

    dt = syms[1]
    subtype = false
    if isa(dt, Symbol)
        try
            dt = Core.eval(__module__, dt)
            if isabstracttype(dt)
                subtype = true
                if length(syms) > 1
                    syms = syms[2:end]
                else
                    throw(ArgumentError("no arguments given for \"$dt\""))
                end
            end
        catch exc
            exc isa UndefVarError || rethrow()
        end
    else
        throw(ArgumentError("\"$dt\" is not a symbol"))
    end
    n = Symbol[]
    for s in syms
        isa(s, Symbol) || throw(ArgumentError("\"$s\" is not a symbol"))
        Base.isidentifier(s) || throw(ArgumentError("\"$s\" is not a valid identifier"))
        s in n && throw(ArgumentError("name \"$s\" is not unique"))
        push!(n, s)
    end
    for m in n
        subtype ?
            Core.eval(__module__, :(struct $m <: $dt end)) :
            Core.eval(__module__, :(struct $m end))
    end
end

"""
    tid(n::Int)

Return a threadid where an `n`-th calculation could run.

This is for emulating a `@threads for` loop with actors.
"""
tid(n::Int) = t = n % nthreads() == 0 ? nthreads() : t


text:
Pkg.init()
#Pkg.add("SQLite3")
#using SQLite
Pkg.add("ODBC")
using ODBC
println("So it begins")

ODBC.connect("tsp.db")
println("Ya me conecte!!")
total_cities = 1092 #10
decreasingTRate = 0.85
batchSize = 400
epsilonTemp = 0.0001
epsilonP = 0.007
constantN = 100
epsilonT = 0.00001
epsilonA = 0.01
initInitialTemp = 0.8
acceptancePercentageParameter = 0.87
cities = [6,9,10,12,13,33,34,35,39,40,43,44,46,50,52,62,63,69,73,74,79,83,86,90,93,94,95,101,104,107,109,111,113,114,115,117,119,121,123,133,137,138,141,142,145,149,151,156,158,159,177,180,181,185,189,193,194,203,212,213,221,222,223,231,233,234,235,242,243,245,246,255,256,257,258,261,262,263,266,269,272,275,278,279,280,281,285,289,290,297,298,299,302,303,305,308,310,313,317,319,322,323,324,326,332,333,334,337,339,341,342,345,346,347,358,359,360,365,366,369,371,372,382,395,398,399,400,404,406,407,408,410,412,414,415,419,422,423,430,431,433,437,438,441,442,447,449,452,453,454,456,462,463,464,469,470,473,475,477,478,480,484,492,494,497,498,499,503,508,515,519,520,526,529,534,535,541,543,545,546,552,562,564,565,566,568,574,575,576,582,583,591,593,598,601,606,611,615,620,621,625,626,627,628,630,637,638,639,644,650,654,664,669,673,675,678,679,681,690,694,695,697,700,701,710,711,717,718,719,722,735,740,741,743,757,759,761,766,768,770,772,773,777,780,781,782,785,792,796,797,803,812,814,817,822,823,829,830,832,837,839,840,844,848,853,854,858,860,863,865,871,872,876,877,880,881,890,892,893,895,897,899,900,906,910,915,916,920,922,923,924,935,936,940,943,955,957,960,964,970,973,978,980,981,991,992,995,996,999,1001,1002,1006,1008,1009,1010,1011,1012,1016,1022,1025,1031,1034,1037,1043,1046,1048,1057,1060,1062,1063,1067,1069,1070,1072,1077,1079,1082,1083,1085,1087,1090,1091] #[22,74,109,113,117,137,178,180,200,216,272,299,345,447,492,493,498,505,521,572,607,627,642,679,710,717,747,774,786,829,830,839,853,857,893,921,935,986,1032,1073] #[2,3,5,6,7,9]
cities_length = length(cities)
seed = 450000#1597
conections = rand(0:seed,total_cities, total_cities)
F = 50
final_solution = []

#db = SQLite.DB("tsp.db")
#println("ME CONECTE!!!!")

#query = "SELECT * FROM cities WHERE id = 12;"
#query(db, query)

#db = SQLite("passwd.sqlite")
#create(db, "passwd", readdlm("/etc/passwd", ':'), ["username", "password", "UID", "GID", "comment", "homedir", "shell"])
#query(db, "SELECT username, homedir FROM passwd LIMIT 10;")
#close(db,

function randpath(n)
    path = 1:n |> collect |> shuffle
    push!(path, path[1]) # loop
    return path
end

function solucion_inicial(s)
    n = rand(1:seed)
    for counter = 0:n
        s = vecino(s)
    end
    return s
end

#used to returs sum of weights, av weigth and punish
function weight_path(s)
    counter = 0
    stop = length(s)-1
    #caminito = zeros(stop)
    for i = 1:length(s) - 1
        peso = conections[s[i],s[i+1]]
        counter = counter + peso
       #caminito[i] = peso
    end
    return counter
end

function av_weight(s)
    return weight_path(s)/length(s) - 1
end

function castigo(s)
    stop = length(s)-1
    caminito = zeros(stop)
    for i = 1:length(s) - 1
        peso = conections[s[i],s[i+1]]
        caminito[i] = peso
    end
    return F * caminito[indmax(caminito)]
end

function overweight(s)
    p = castigo(s)
    stop = length(s)-1
    caminito = zeros(stop)
    counter = 0
    for i = 1:length(s) - 1
        peso = conections[s[i],s[i+1]]
        if peso == 0
            caminito[i] = p
            counter = counter + p
        else
            caminito[i] = peso
            counter = counter + peso
        end
    end
    return counter
end

function calculaLote(T, s)
    c = 0
    r = 0.0
    while c < batchSize
        s1 = vecino(s)
        f_s1 = f_costo(s1)
        f_s = f_costo(s)
        if  f_s1 <= f_s + T
            s = s1
            c = c + 1
            r = r + f_s1
            if f_s1 < f_s
                final_solution = s1
            end
        end
    end
    return r/batchSize, s
end

function vecino(s)
    i,j = rand_positions()
    return swap(s,i,j)
end

function swap(arr, i, j)
    temp = arr[i]
    arr[i] = arr[j]
    arr[j] = temp
    return arr
end

function rand_positions()
    i = j = 0
    while i == j
        i = rand(1:cities_length)
        j = rand(1:cities_length)
    end
    return i,j
end

function f_costo(s)
    suma = overweight(s)
    return suma / av_weight(s)
    return
end

function aceptacionPorUmbrales(T, s)
    p = 0
    q = -1
    while T > epsilonT
        p1 = q
        while p <= p1 || q == -1
            p1 = p
            p, s = calculaLote(T, s)
            q = p
            println(string("nueva solucion: ", s))
            println(string("funcion de costo: ", p))
        end
        T = decreasingTRate * T
        println(string("nueva temperatura : ",T))
    end
    return
end

function simulated_annealing{T <: Real}(distmat::Matrix{T}; steps = 50*length(distmat),
                                        num_starts = 1,
                                        init_temp = exp(8), final_temp = exp(-6.5),
                                        init_path::Nullable{Vector{Int}} = Nullable{Vector{Int}}())

    cool_rate = (final_temp / init_temp)^(1 / (steps - 1))

    function sahelper()
        temp = init_temp / cool_rate
        n = size(distmat, 1)
        path = isnull(init_path) ? randpath(n) : copy(get(init_path))
        cost_cur = pathcost(distmat, path)

        for i in 1:steps
            temp *= cool_rate

            first, last = rand(2:n), rand(2:n)
            if first > last
                first, last = last, first
            end
            cost_other = pathcost_rev(distmat, path, first, last)
            @fastmath accept = cost_other < cost_cur ? true : rand() < exp((cost_cur - cost_other) / temp)
            if accept
                reverse!(path, first, last)
                cost_cur = cost_other
            end
        end

        return path, cost_cur
    end

    path, cost = sahelper()
    for _ in 2:num_starts
        otherpath, othercost = sahelper()
        if othercost < cost
            cost = othercost
            path = otherpath
        end
    end

    return path, cost
end

#r = [1,2,3,4,5,6]
#print(swap(r, 1, 4))
#print(rand_positions())
#path = randpath(cities_length)
#println(cities)
#println(conections)
#println(castigo(cities))
#SOLUTION = solucion_inicial(cities)
#final_solution = SOLUTION
#println(SOLUTION)
#println(weight_path(SOLUTION))
#println(overweight(SOLUTION))
#println(av_weight(SOLUTION))
#println(f_costo(SOLUTION))
#aceptacionPorUmbrales(initInitialTemp, SOLUTION)
#println(string("solucion final: ", final_solution))
#println(string("costo final: ", f_costo(final_solution)))


text:
module Klio

end


text:
using Pkg # hideall
Pkg.activate("_literate/D0-dataframe/Project.toml")
Pkg.update()

# This tutorial is loosely adapted from [this pandas tutorial](https://pandas.pydata.org/pandas-docs/stable/getting_started/10min.html) as well as the [DataFrames.jl documentation](http://juliadata.github.io/DataFrames.jl/latest/man/getting_started/).
# It is by no means meant to be a complete introduction, rather, it focuses on some key functionalities that are particularly useful in a classical machine learning context.
#
# ## Basics
#
# To start with, we will use the `Boston` dataset which is very simple.

using RDatasets
using DataFrames

boston = dataset("MASS", "Boston");

# The `dataset` function returns a `DataFrame` object:

typeof(boston)

# ### Accessing data
#
# Intuitively a DataFrame is just a wrapper around a number of columns, each of which is a `Vector` of some type with a name"

names(boston)

# You can view the first few rows using `first` and specifying a number of rows:

first(boston, 4)

# You can access one of those columns easily using `.colname`, this returns a vector that you can access like any Julia vector:

boston.Crim[1:5]

# You can also just access the dataframe as you would a big matrix:

boston[3, 5]

# or specifying a range of rows/columns:

boston[1:5, [:Crim, :Zn]]

# or, similarly,

boston[1:5, 1:2]

# The `select` function is very convenient to get sub dataframes of interest:

b1 = select(boston, [:Crim, :Zn, :Indus])
first(b1, 2)

# The `Not` syntax is  also very  useful:

b2 = select(boston, Not(:NOx))
first(b2, 2)

# Finally, if you would like to drop columns, you can use `select!` which will mutate the dataframe in place:

select!(b1, Not(:Crim))
first(b1, 2)

# ### Describing the data
#
# `StatsBase` offers a convenient `describe` function which you can use on a DataFrame to get an overview of the data:

using StatsBase
describe(boston, :min, :max, :mean, :median, :std)

# You can pass a number of symbols to the `describe` function to indicate which statistics to compute for each feature:
#
# * `mean`, `std`, `min`, `max`, `median`, `first`, `last` are all fairly self explanatory
# * `q25`, `q75` are respectively for the 25th and 75th percentile,
# * `eltype`, `nunique`, `nmissing` can also be used
#
# You can also  pass your custom function with a pair `name => function` for instance:

foo(x) = sum(abs.(x)) / length(x)
d = describe(boston, :mean, :median, foo => :foo)
first(d, 3)

# The `describe` function returns a derived object with one row per feature and one column per required statistic.
#
# Further to `StatsBase`, `Statistics` offers a range of useful functions for data analysis.

using Statistics

#
# ### Converting the data
#
# If you want to get the content of the dataframe as one big matrix, use `convert`:

mat = Matrix(boston)
mat[1:3, 1:3]

# ### Adding columns
#
# Adding a column to a dataframe is very easy:

boston.Crim_x_Zn = boston.Crim .* boston.Zn;

# that's it! Remember also that you can drop columns or make subselections with `select` and `select!`.

# ### Missing values
#
# Let's load a dataset with missing values

mao = dataset("gap", "mao")
describe(mao, :nmissing)

# Lots of missing values...
# If  you wanted to compute simple functions on columns, they  may just return `missing`:

std(mao.Age)

# The `skipmissing` function can help counter this  easily:

std(skipmissing(mao.Age))

#
# ## Split-Apply-Combine
#
# This is a shorter version of the [DataFrames.jl tutorial](http://juliadata.github.io/DataFrames.jl/latest/man/split_apply_combine/).

iris = dataset("datasets", "iris")
first(iris, 3)

# ### `groupby`
#
# The `groupby` function allows to form "sub-dataframes" corresponding to groups of rows.
# This can be very convenient to run specific analyses for specific groups without copying the data.
#
# The basic usage is `groupby(df, cols)` where `cols` specifies one or several columns to use for the grouping.
#
# Consider a simple example: in `iris` there is a `Species` column with 3 species:

unique(iris.Species)

# We can form views for each of these:

gdf = groupby(iris, :Species);

# The `gdf` object now corresponds to **views** of the original dataframe for each of the 3 species; the first species is `"setosa"` with:

subdf_setosa = gdf[1]
describe(subdf_setosa, :min, :mean, :max)

# Note that `subdf_setosa` is a `SubDataFrame` meaning that it is just a view of the parent dataframe `iris`; if you modify that parent dataframe then the sub dataframe is also  modified.
#
# See `?groupby` for more information.
#
# ### `combine`
#
# The `combine` function allows to derive a new dataframe out of transformations of an existing one.
# Here's an example taken from the official doc (see `?combine`):

df = DataFrame(a=1:3, b=4:6)
combine(df, :a => sum, nrow)

# what happened here is that the derived DataFrame has two columns obtained respectively by (1) computing the sum of the first column and (2) applying the `nrow` function on the `df`.
#
# The transformation can produce one or several values, `combine` will try to concatenate these columns as it can, for instance:

foo(v) = v[1:2]
combine(df, :a => maximum, :b => foo)

# here the maximum value of `a` is copied twice so that the two columns have the same number of rows.
#

bar(v) = v[end-1:end]
combine(df, :a => foo, :b => bar)

# ### `combine` with `groupby`
#
# Combining `groupby` with `combine` is very useful.
# For instance you might want to compute statistics across groups for different variables:

combine(groupby(iris, :Species), :PetalLength => mean)

# let's decompose that:
#
# 1. the `groupby(iris, :Species)` creates groups using the `:Species` column (which has values `setosa`, `versicolor`, `virginica`)
# 2. the `combine` creates a derived dataframe by applying the `mean` function to the `:PetalLength` column
# 3. since there are three groups, we get one column (mean of `PetalLength`) and three rows (one per group).
#
#
# You can do this for several columns/statistics at the time and give new column names to the results:

gdf = groupby(iris, :Species)
combine(gdf, :PetalLength => mean => :MPL, :PetalLength => std => :SPL)

# so here we assign the names `:MPL` and `:SPL` to the derived columns.
# If you want to apply something on all columns apart from the grouping one, using `names` and `Not` comes in handy:

combine(gdf, names(iris, Not(:Species)) .=> std)

#
# where

names(iris, Not(:Species))

#
# and note the use of `.` in `.=>` to indicate that we broadcast the function over each column.

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
module RandomExtensions

export make, Uniform, Normal, Exponential, CloseOpen, OpenClose, OpenOpen, CloseClose, Rand, Bernoulli

# re-exports from Random, which don't overlap with new functionality and not from misc.jl
export rand!, AbstractRNG, MersenneTwister, RandomDevice

import Random: Sampler, rand, rand!, gentype

using Random
using Random: GLOBAL_RNG, SamplerTrivial, SamplerSimple, SamplerTag, SamplerType, Repetition

using SparseArrays: sprand, sprandn, AbstractSparseArray, SparseVector, SparseMatrixCSC

## a dummy container type to take advangage of SamplerTag constructor

struct Cont{T} end

Base.eltype(::Type{Cont{T}}) where {T} = T

## some helper functions, not to be overloaded, except default_sampling

default_sampling(::Type{X}) where {X} = error("default_sampling($X) not defined")
default_sampling(::X)       where {X} = default_sampling(X)

default_gentype(::Type{T}) where {T} = val_gentype(default_sampling(T))
default_gentype(::X)       where {X} = default_gentype(X)

val_gentype(X)                   = gentype(X)
val_gentype(::Type{X}) where {X} = X

## includes

include("distributions.jl")
include("sampling.jl")
include("containers.jl")
include("iteration.jl")

## updated rand docstring (TODO: replace Base's one)

"""
    rand([rng=GLOBAL_RNG], [S], [C...]) # RandomExtensions

Pick a random element or collection of random elements from the set of values specified by `S`;
`S` can be

* an indexable collection (for example `1:n` or `['x','y','z']`),
* an `AbstractDict` or `AbstractSet` object,
* a string (considered as a collection of characters), or
* a type: the set of values to pick from is then equivalent to `typemin(S):typemax(S)` for
  integers (this is not applicable to [`BigInt`](@ref)), and to ``[0, 1)`` for floating
  point numbers;
* a `Distribution` object, e.g. `Normal()` for a normal distribution (like `randn()`),
  or `CloseOpen(10.0, 20.0)` for uniform `Float64` numbers in the range ``[10.0, 20.0)``;
* a `make` object, which can be e.g. `make(Pair, S1, S2)` or `make(Complex, S1, S2)`,
  where `S1` and `S2` are one of the specifications above; `Pair` or `Complex` can optionally be
  given as concrete types, e.g. `make(ComplexF64, 1:3, Int)` to generate `ComplexF64` instead
  of `Complex{Int}`.

`S` usually defaults to [`Float64`](@ref).

If `C...` is not specified, `rand` produces a scalar. Otherwise, `C` can be:

* a set of integers, or a tuple of `Int`, which specify the dimensions of an `Array` to generate;
* `(Array, dims...)`: same as above, but with `Array` specified explicitely;
* `(p::AbstractFloat, m::Integer, [n::Integer])`, which produces a sparse array of dimensions `(m, n)`,
  in which the probability of any element being nonzero is independently given by `p`;
* `(String, [n=8])`, which produces a random `String` of length `n`; the generated string consists of `Char`
  taken from a predefined set like `randstring`, and can be specified with the `S` parameter;
* `(Dict, n)`, which produces a `Dict` of length `n`; if `Dict` is given without type parameters,
  then `S` must be specified;
* `(Set, n)` or `(BitSet, n)`, which produces a set of length `n`;
* `(BitArray, dims...)`, which produces a `BitArray` with the specified dimensions.

For `Array`, `Dict` and `Set`, a less abstract type can be specified, e.g. `Set{Float64}`, to force
the type of the result regardless of the `S` parameter. In particular, in the absence of `S`, the
type parameter(s) of the container play the role of `S`; for example, `rand(Dict{Int,Float64}, n)`
is equivalent to `rand(make(Pair, Int, Float64), Dict, n)`.

# Examples
```julia-repl
julia> rand(Int, 2)
2-element Array{Int64,1}:
 1339893410598768192
 1575814717733606317

julia> rand(MersenneTwister(0), Dict(1=>2, 3=>4))
1=>2

julia> rand("abc", String, 12)
"bccaacaabaac"

julia> rand(1:10, Set, 3)
Set([3, 8, 6])

julia> rand(1:10, Set{Float64}, 3)
Set([10.0, 5.0, 6.0])
```

!!! note
    The complexity of `rand(rng, s::Union{AbstractDict,AbstractSet})`
    is linear in the length of `s`, unless an optimized method with
    constant complexity is available, which is the case for `Dict`,
    `Set` and `BitSet`. For more than a few calls, use `rand(rng,
    collect(s))` instead, or either `rand(rng, Dict(s))` or `rand(rng,
    Set(s))` as appropriate.
"""
rand

end # module


text:
println("pid: $(getpid())")

using PyCall
using BSON

# pyimports
plt = pyimport("matplotlib.pyplot")
cytofpy = pyimport(:cytofpy)
torch = pyimport(:torch)
Gamma = torch.distributions.Gamma
Beta = torch.distributions.Beta

SIMDAT_PATH = "../../../sims/sim_study/simdata/kills-flowsom/N500/K5/90/simdat.bson"

simdat = BSON.load(SIMDAT_PATH)[:simdat]
y = [torch.tensor(yi) for yi in simdat[:y]]
priors = cytofpy.model.default_priors(y, K=30, L=[5, 3],
                                      y_quantiles=[0.02, .25, .48],
                                      p_bounds=[.05, .8, .05])
priors["sig2"] = Gamma(.1, 1)
priors["alpha"] = Gamma(.1, .1)
priors["delta0"] = Gamma(1, 1)
priors["delta1"] = Gamma(1, 1)
priors["noisy_var"] = 10.0
priors["eps"] = Beta(100, 9900)

every(max_iter, nsave) = round(Int, max_iter / nsave)
max_iter = 10000
nsave = 30
out = cytofpy.model.fit(y, max_iter=max_iter, lr=1e-2, print_freq=10, eps_conv=0,
                        priors=priors, minibatch_size=500, tau=0.005,
                        trace_every=50, backup_every=every(max_iter, nsave),
                        verbose=0, seed=1, use_stick_break=false)

out = Dict(Symbol(k) => v for (k, v) in out)
# BSON.bson("results/out.bson", out)

#= Plots
import Distributions
out = BSON.load("results/out.bson")
mp = out[:mp]
v = mp["v"]
H = mp["H"]
Z = [v.dist().rsample().reshape(1, priors["K"]) > H.dist().rsample() for i in 1:100]
Z_mean = Distributions.mean(Z)
eps_samps = [mp["eps"].transform(mp["eps"].dist().rsample()) for i in 1:100]

plt.plot(out[:elbo]); plt.show()
plt.imshow(Z_mean.numpy()); plt.show()
=#


text:
type TwoFactorShortRateTree{S <: ShortRateDynamics} <: ShortRateTree
  tree1::TrinomialTree
  tree2::TrinomialTree
  dynamics::S
  treeLattice::TreeLattice2D

  function TwoFactorShortRateTree{S}(tree1::TrinomialTree, tree2::TrinomialTree, dyn::S)
    twoFactorTree = new(tree1, tree2, dyn)
    twoFactorTree.treeLattice = TreeLattice2D(tree1, tree2, dyn.correlation, twoFactorTree)

    return twoFactorTree
  end
end

function tree(model::TwoFactorModel, grid::TimeGrid)
  #stuff
  dyn = get_dynamics(model)

  tree1 = TrinomialTree(dyn.xProcess, grid)
  tree2 = TrinomialTree(dyn.yProcess, grid)

  return TwoFactorShortRateTree{typeof(dyn)}(tree1, tree2, dyn)
end

get_size(tr::TwoFactorShortRateTree, i::Int) = get_size(tr.treeLattice, i)

descendant(tr::TwoFactorShortRateTree, i::Int, idx::Int, branch::Int) = descendant(tr.treeLattice, i, idx, branch)
probability(tr::TwoFactorShortRateTree, i::Int, idx::Int, branch::Int) = probability(tr.treeLattice, i, idx, branch)

function discount(tr::TwoFactorShortRateTree, i::Int, idx::Int)
  modulo = get_size(tr.tree1, i)
  new_idx = idx - 1

  index1 = round(Int, floor(new_idx % modulo)) + 1
  index2 = round(Int, floor(new_idx / modulo)) + 1

  x = get_underlying(tr.tree1, i, index1)
  y = get_underlying(tr.tree2, i, index2)

  r = short_rate(tr.dynamics, tr.treeLattice.tg.times[i], x, y)

  # println("$i $idx $r $x $y $modulo $index1 $index2")

  return exp(-r * tr.treeLattice.tg.dt[i])
end

type G2Dynamics{P <: Parameter} <: ShortRateDynamics
  fitting::P
  xProcess::OrnsteinUhlenbeckProcess
  yProcess::OrnsteinUhlenbeckProcess
  correlation::Float64
end

G2Dynamics{P <: Parameter}(fitting::P, a::Float64, sigma::Float64, b::Float64, eta::Float64, rho::Float64) =
          G2Dynamics{P}(fitting, OrnsteinUhlenbeckProcess(a, sigma), OrnsteinUhlenbeckProcess(b, eta), rho)

short_rate(dyn::G2Dynamics, t::Float64, x::Float64, y::Float64) = dyn.fitting(t) + x + y

type G2{AffineModelType, T <: TermStructure} <: TwoFactorModel{AffineModelType}
  modT::AffineModelType
  a::ConstantParameter
  sigma::ConstantParameter
  b::ConstantParameter
  eta::ConstantParameter
  rho::ConstantParameter
  phi::G2FittingParameter
  ts::T
  privateConstraint::PrivateConstraint
  common::ModelCommon
end

function G2{T <: TermStructure}(ts::T, a::Float64 = 0.1, sigma::Float64 = 0.01, b::Float64 = 0.1, eta::Float64 = 0.01, rho::Float64 = -0.75)
  a_const = ConstantParameter([a], PositiveConstraint())
  sigma_const = ConstantParameter([sigma], PositiveConstraint())
  b_const = ConstantParameter([b], PositiveConstraint())
  eta_const = ConstantParameter([eta], PositiveConstraint())
  rho_const = ConstantParameter([rho], BoundaryConstraint(-1.0, 1.0))

  privateConstraint = PrivateConstraint(ConstantParameter[a_const, sigma_const, b_const, eta_const, rho_const])

  phi = G2FittingParameter(a, sigma, b, eta, rho, ts)

  return G2{AffineModelType, T}(AffineModelType(), a_const, sigma_const, b_const, eta_const, rho_const, phi, ts, privateConstraint, ModelCommon())
end

get_eta(m::G2) = m.eta.data[1]
get_rho(m::G2) = m.rho.data[1]

get_params(m::G2) = Float64[get_a(m), get_sigma(m), get_b(m), get_eta(m), get_rho(m)]

generate_arguments!(m::G2) = m.phi = G2FittingParameter(get_a(m), get_sigma(m), get_b(m), get_eta(m), get_rho(m), m.ts)

get_dynamics(m::G2) = G2Dynamics(m.phi, get_a(m), get_sigma(m), get_b(m), get_eta(m), get_rho(m))

function V(m::G2, t::Float64)
  expat = exp(-get_a(m) * t)
  expbt = exp(-get_b(m) * t)
  cx = get_sigma(m) / get_a(m)
  cy = get_eta(m) / get_b(m)
  valuex = cx * cx * (t + (2.0 * expat - 0.5 * expat * expat - 1.5) / get_a(m))
  valuey = cy * cy * (t + (2.0 * expbt - 0.5 * expbt * expbt - 1.5) / get_b(m))
  value = 2.0 * get_rho(m) * cx * cy * (t + (expat - 1.0) / get_a(m) + (expbt - 1.0) / get_b(m) - (expat * expbt - 1.0) / (get_a(m) + get_b(m)))

  return valuex + valuey + value
end

A(m::G2, t::Float64, T::Float64) = discount(m.ts, T) / discount(m.ts, t) * exp(0.5 * (V(m, T - t) - V(m, T) + V(m, t)))
B(m::G2, x::Float64, t::Float64) = (1.0 - exp(-x * t)) / x

type G2SwaptionPricingFunction
  a::Float64
  sigma::Float64
  b::Float64
  eta::Float64
  rho::Float64
  w::Float64
  startTime::Float64
  fixedRate::Float64
  sigmax::Float64
  sigmay::Float64
  rhoxy::Float64
  mux::Float64
  muy::Float64
  payTimes::Vector{Float64}
  A::Vector{Float64}
  Ba::Vector{Float64}
  Bb::Vector{Float64}
end

function G2SwaptionPricingFunction(model::G2, w::Float64, startTime::Float64, payTimes::Vector{Float64}, fixedRate::Float64)
  a = get_a(model)
  sigma = get_sigma(model)
  b = get_b(model)
  eta = get_eta(model)
  rho = get_rho(model)

  sigmax = sigma * sqrt(0.5 * (-expm1(-2.0 * a * startTime)) / a)
  sigmay = eta * sqrt(0.5 * (-expm1(-2.0 * b * startTime)) / b)
  rhoxy = rho * eta * sigma * (-expm1(-(a + b) * startTime)) / ((a + b) * sigmax * sigmay)

  temp = sigma * sigma / (a * a)
  mux = -((temp + rho * sigma * eta / (a * b)) * (-expm1(-a * startTime)) - 0.5 * temp * (-expm1(-2.0 * a * startTime)) -
        rho * sigma * eta / (b * (a + b)) * (-expm1(-(b + a) * startTime)))

  temp = eta * eta / (b * b)
  muy = -((temp + rho * sigma * eta / (a * b)) * (-expm1(-b * startTime)) - 0.5 * temp * (-expm1(-2.0 * b * startTime)) -
        rho * sigma * eta / (a * (a + b)) * (-expm1(-(b + a) * startTime)))

  n = length(payTimes)
  A_ = zeros(n)
  Ba = zeros(n)
  Bb = zeros(n)
  for i = 1:n
    A_[i] = A(model, startTime, payTimes[i])
    Ba[i] = B(model, a, payTimes[i] - startTime)
    Bb[i] = B(model, b, payTimes[i] - startTime)
  end

  return G2SwaptionPricingFunction(a, sigma, b, eta, rho, w, startTime, fixedRate, sigmax, sigmay, rhoxy, mux, muy, payTimes, A_, Ba, Bb)
end

function operator(pricingFunc::G2SwaptionPricingFunction)
  phi = Normal()
  n = length(pricingFunc.payTimes)
  function _inner(x::Float64)
    temp = (x - pricingFunc.mux) / pricingFunc.sigmax
    txy = sqrt(1.0 - pricingFunc.rhoxy * pricingFunc.rhoxy)
    lambda = zeros(n)

    for i = 1:n
      tau = i == 1 ? pricingFunc.payTimes[1] - pricingFunc.startTime : pricingFunc.payTimes[i] - pricingFunc.payTimes[i - 1]
      c = i == n ? (1.0 + pricingFunc.fixedRate * tau) : pricingFunc.fixedRate * tau
      lambda[i] = c * pricingFunc.A[i] * exp(-pricingFunc.Ba[i] * x)
    end

    func = SolvingFunction(lambda, pricingFunc.Bb)
    s1d = BrentSolver(1000)

    yb = solve(s1d, operator(func), 1e-6, 0.00, -100.0, 100.0)

    h1 = (yb - pricingFunc.muy) / (pricingFunc.sigmay * txy) - pricingFunc.rhoxy * (x - pricingFunc.mux) / (pricingFunc.sigmax * txy)
    val = cdf(phi, -pricingFunc.w * h1)

    for i = 1:n
      h2 = h1 + pricingFunc.Bb[i] * pricingFunc.sigmay * sqrt(1.0 - pricingFunc.rhoxy * pricingFunc.rhoxy)
      kappa = -pricingFunc.Bb[i] * (pricingFunc.muy - 0.5 * txy * txy * pricingFunc.sigmay * pricingFunc.sigmay * pricingFunc.Bb[i] +
              pricingFunc.rhoxy * pricingFunc.sigmay * (x - pricingFunc.mux) / pricingFunc.sigmax)

      val -= lambda[i] * exp(kappa) * cdf(phi, -pricingFunc.w * h2)
    end

    blah = exp(-0.5 * temp * temp) * val / (pricingFunc.sigmax * sqrt(2.0 * pi))
    return exp(-0.5 * temp * temp) * val / (pricingFunc.sigmax * sqrt(2.0 * pi))
  end

  return _inner
end

function gen_swaption(model::G2, swaption::Swaption, fixedRate::Float64, range::Float64, intervals::Int)
  settlement = reference_date(model.ts)
  dc = model.ts.dc
  startTime = year_fraction(dc, settlement, swaption.swap.args.floatingResetDates[1])
  w = swaption.swap.payer[2]

  fixedPayTimes = zeros(length(swaption.swap.args.fixedPayDates))
  for i = 1:length(fixedPayTimes)
    fixedPayTimes[i] = year_fraction(dc, settlement, swaption.swap.args.fixedPayDates[i])
  end
  func = G2SwaptionPricingFunction(model, w, startTime, fixedPayTimes, fixedRate)

  upper = func.mux + range * func.sigmax
  lower = func.mux - range * func.sigmax

  integrator = SegmentIntegral(intervals)
  return swaption.swap.nominal * w * discount(model.ts, startTime) * QuantLib.Math.operator(integrator, QuantLib.operator(func), lower, upper)
end

discount_bond(model::G2, t::Float64, T::Float64, x::Float64, y::Float64) = A(model, t, T) * exp(-B(model, get_a(model), (T - t)) * x - B(model, get_b(model), (T - t)) * y)
discount_bond(model::G2, now::Float64, maturity::Float64, factors::Vector{Float64}) = discount_bond(model, now, maturity, factors[1], factors[2])


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(value), 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)

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, Val(N)))
    (atoms === nothing) && (atoms = ntuple(i->wildcard, Val(N)))
    subscript = isa(orbitals, NTuple{N, Int}) ? Subscript(orbitals) : isa(orbitals, Nothing) ? Subscript{N}() : orbitals
    orbitals = subscript.opattern
    subscripts = ntuple(i->subscript.identifier, Val(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
    subscripts = sc.id.subscripts
    all(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, "⊗"))
    centers = sc.id.centers
    any(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
    fpos = findfirst(r"S[+-xyz0]", ps[1]).start
    lpos = 1 + lastindex(ps[1]) - findfirst(r"[+-xyz0]S", ps[1]|>reverse).start
    coeff = eval(Meta.parse(ps[1][1:fpos-1]))
    tags = Tuple(replace(ps[1][fpos:lpos], "S"=>""))
    attrpairs = Any[:tags=>tags]
    temp = replace(ps[1][lpos+1:end], " "=>"")
    if length(temp) > 0
        components = split(temp, '@')
        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
    str = replace(str, " "=>"")
    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)
    slice = findfirst(r"xyz|\+\-z", str)
    mode = isnothing(slice) ? "+-z" : str[slice]
    isnothing(slice) || (str = str[slice.stop+1:end])
    @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] ∈ ('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