text:
###########################################################
##
## MAPS
##
## maps between residue rings of polynomials and in general
##
###########################################################

function domain(M::Map)
  return M.header.domain
end

function codomain(M::Map)
  return M.header.codomain
end

function image_function(f::Map)
  return f.header.image
end

function preimage_function(f::Map)
  return f.header.preimage
end

export Map, CoerceMap, ResidueRingPolyMap

type MapHeader{D, C}
  domain::D
  codomain::C
  image::Function
  preimage::Function

  function MapHeader()
    z = new{D, C}()
    return z
  end

  function MapHeader(domain::D, codomain::C)
    z = new{D, C}()
    z.domain = domain
    z.codomain = codomain
    return z
  end

  function MapHeader(domain::D, codomain::C, image::Function)
    z = new{D, C}()
    z.domain = domain
    z.codomain = codomain
    z.image = image
    return z
  end

  function MapHeader(domain::D, codomain::C, image::Function, preimage::Function)
    z = new{D, C}()
    z.domain = domain
    z.codomain = codomain
    z.image = image
    z.preimage = preimage
    return z
  end
end

function MapHeader{D, C}(domain::D, codomain::C, image::Function)
  return MapHeader{D, C}(domain, codomain, image)
end

function MapHeader{D, C}(domain::D, codomain::C, image::Function, preimage::Function)
  return MapHeader{D, C}(domain, codomain, image, preimage)
end
# this type represents a -> f(g(a))
type CompositeMap{D, C, R} <: Map{D, C}
  header::MapHeader{D, C}
  f::Map{R, C}
  g::Map{D, R}

  function CompositeMap(f::Map, g::Map)
  ##CF should be function CompositeMap(f::Map{R, C}, g::Map{D, R})
  ## but that seems to not work:
  # U, m = UnitGroup(ResidueRing(ZZ, 2^9));
  # q, mq = Hecke.quo(U, [preimage(m, codomain(m)(fmpz(-1)))])
  # z = Hecke.compose(m, inv(mq))
  # btw: m*inv(mq) also fails.

    domain(f) == codomain(g) || throw("maps not compatible")
    z = new{D, C, R}()

    image = function(x)#x::elem_type(domain(z)))
      parent(x) != domain(g) && error("Element not in domain of map")
      return f(g(x))::elem_type(codomain(z))
    end

    if isdefined(f.header, :preimage) && isdefined(g.header, :preimage)
      _preimage = function(x)#x::elem_type(codomain(z)))
        return preimage(g, preimage(f, x))::elem_type(domain(z))
      end
      z.header = MapHeader(domain(g), codomain(f), image, _preimage)
    else
      z.header = MapHeader(domain(g), codomain(f), image)
    end

    return z
  end
end

function *{D, C, R}(f::Map{R, C}, g::Map{D, R})
  return CompositeMap{D, C, R}(f, g)
end

function compose{D, C, R}(f::Map{R, C}, g::Map{D, R})
  return CompositeMap{D, C, R}(f, g)
end

function compose(f::Map, g::Map)
  return CompositeMap{typeof(domain(g)), typeof(codomain(f)), Any}(f, g)
end

type InverseMap{D, C} <: Map{D, C}
  header::MapHeader{D, C}
  origin::Map{C, D}

  function InverseMap(f::Map{C, D})
    z = new{D, C}()
    z.header = MapHeader(codomain(f), domain(f), preimage_function(f), image_function(f))
    z.origin = f
    return z
  end
end

function InverseMap{D, C}(f::Map{C, D})
  return InverseMap{D, C}(f)
end

type ResidueRingPolyMap{D, C, T} <: Map{D, C}
  header::MapHeader{D, C}
  gen_image::GenRes{T}
  coeff_map::Map # can be missing if domain and codomain have the same
                 # base_ring(base_ring())
  function ResidueRingPolyMap(domain::D, codomain::C, gen_image::GenRes{T}, coeff_map::Map)
    z = new{D, C, T}()
    z.gen_image = gen_image
    z.coeff_map = coeff_map

    image = function(a::GenRes)
      #a should be in the domain of M...
      I = codomain(0)
      for i in degree(a.data):-1:0
        I = I*z.gen_image + z.coeff_map(coeff(a.data, i))
      end
      return I::elem_type(C)
    end

    # I need to call preimage in _preimage
    _preimage = function(a::GenRes)
      R = codomain
      parent(a) == R || throw("mixed rings in preimage")
      g = gens(domain)
      im_gen = map(z, g) # apply x -> z(x) to the generatotrs ## possibly should be cached and stored
          ## need to write the elements in a matrix, solve the eqn for a
      Mt = MatrixSpace(base_ring(base_ring(R)), degree(R.modulus), degree(R.modulus))()

      for i=1:degree(R.modulus)
        elem_to_mat_row!(Mt, i, im_gen[i])
      end
      b = MatrixSpace(base_ring(base_ring(R)), 1, degree(R.modulus))()
      elem_to_mat_row!(b, 1, a)
      s = solve(Mt', b') # why, oh why is solve operating on columns????
      # This is the worst
      if isa(s, Tuple) ## again, why, oh why is solve doing things differently
                   ## over rings than fields?
        s = s[1] * inv(s[2]) # all rings here (type) are actually fields (math)
      end

      s = sum([preimage(z.coeff_map, s[i,1])*g[i] for i=1:length(g)])
      @assert parent(s) == domain
      return s::elem_type(domain)
    end

    z.header = MapHeader(domain, codomain, image, _preimage)
    return z
  end

  function ResidueRingPolyMap(domain::D, codomain::C, gen_image::GenRes{T})
    z = new{D, C, T}()
    z.gen_image = gen_image

    image = function(a::GenRes)
      I = z.codomain(0)
      for i in degree(a.data):-1:0
        I = I*z.gen_image + coeff(a.data, i)
      end
      return I::elem_type(C)
    end

    preimage = function(a::GenRes)
      R = z.codomain
      parent(a) == R || throw("mixed rings in preimage")
      g = gens(domain)
      im_gen = map(x -> z(x), g) # apply x -> z(x) to the generatotrs ## possibly should be cached and stored
          ## need to write the elements in a matrix, solve the eqn for a
      Mt = MatrixSpace(base_ring(base_ring(R)), degree(R.modulus), degree(R.modulus))()

      for i=1:degree(R.modulus)
        elem_to_mat_row!(Mt, i, im_gen[i])
      end
      b = MatrixSpace(base_ring(base_ring(R)), 1, degree(R.modulus))()
      elem_to_mat_row!(b, 1, a)
      s = solve(Mt', b') # why, oh why is solve operating on columns????
      # This is the worst
      if isa(s, Tuple) ## again, why, oh why is solve doing things differently
                   ## over rings than fields?
        s = s[1] * inv(s[2]) # all rings here (type) are actually fields (math)
      end

      s = sum([s[i,1]*g[i] for i=1:length(g)])

      @assert parent(s) == domain
      return s::elem_type(domain)
    end

    z.header = MapHeader(domain, codomain, image, preimage)
    return z
  end
end

function ResidueRingPolyMap{D, C, T}(domain::D, codomain::C, i::GenRes{T})
  return ResidueRingPolyMap{D, C, T}(domain, codomain, i)
end

function ResidueRingPolyMap{D, C, T}(domain::D, codomain::C, i::GenRes{T}, coeff_map::Map)
  return ResidueRingPolyMap{D, C, T}(domain, codomain, i, coeff_map)
end

type IdentityMap{D} <: Map{D, D}
  header::MapHeader{D, D}

  function IdentityMap(domain::D)
    z = new{D}()

    image = function(x)# (x::elem_type(D))
      return x::elem_type(D)
    end
    preimage = image

    z.header = MapHeader(domain, domain, image, preimage)
    return z
  end
end

function IdentityMap{D}(domain::D)
  return IdentityMap{D}(domain)
end

type CoerceMap{D, C} <: Map{D, C}
  header::MapHeader{D, C}

  function CoerceMap(domain::D, codomain::C)
    z = new{D, C}()
    z.header = MapHeader(domain, codomain)
    return z
  end

  function CoerceMap(domain::GenResRing{fmpz}, codomain::FqNmodFiniteField)
    z = new{GenResRing{fmpz}, FqNmodFiniteField}()

    image = function(a::GenRes{fmpz})
        parent(a) != domain && error("Element not in domain")
        return codomain(ZZ(a))
    end

    preimage = function(a::fq_nmod)
      parent(a) != codomain && error("Element not in codomain")
      a.length > 1 && throw("Element not in image")
      return domain(coeff(a, 0))::GenRes{fmpz}
    end

    z.header = MapHeader(domain, codomain, image, preimage)
    return z
  end

  function CoerceMap(domain::FqNmodFiniteField, codomain::GenResRing{fq_nmod_poly})
    z = new{FqNmodFiniteField, GenResRing{fq_nmod_poly}}()

    image = function(a::fq_nmod)
      parent(a) != domain && error("Element not in domain")
      return codomain(a)::GenRes{fq_nmod_poly}
    end

    preimage = function(a::GenRes{fq_nmod_poly})
      degree(a.data) > 0 && throw("Element not in subfield")
      return domain(coeff(a.data, 0))::fq_nmod
    end

    z.header = MapHeader(domain, codomain, image, preimage)
    return z
  end

  function CoerceMap{S, T <: PolyElem}(domain::GenResRing{S}, codomain::GenResRing{T})
    z = new{GenResRing{S}, GenResRing{T}}()

    image = function(a::GenRes)
      return codomain(a)::elem_type(codomain)
    end

    preimage = function(a::GenRes)
      while parent(a) != domain
        degree(a.data)>0 && throw("Element not in subfield")
        a = coeff(a.data, 0)
      end
      return a::elem_type(domain)
    end

    z.header = MapHeader(domain, codomain, image, preimage)
    return z
  end
end

function CoerceMap{D, C}(domain::D, codomain::C)
  return CoerceMap{D, C}(domain, codomain)
end

type FinGenGrpAbMap{D <: FinGenGrpAb, C <: FinGenGrpAb} <: Map{D, C}
  header::MapHeader{D, C}

  map::fmpz_mat
  imap::fmpz_mat
  im:: FinGenGrpAb  # if set
  ke:: FinGenGrpAb  # if set

  function FinGenGrpAbMap(From::D, To::C, M::fmpz_mat)
    r = new()
    function image(a::FinGenGrpAbElem{D})
      return FinGenGrpAbElemCreate(To, a.coeff*M)
    end

    function preimage(a::FinGenGrpAbElem{C})
      error("missing")
    end

    r.header = MapHeader(From, To, image, preimage)
    return r
  end

  function FinGenGrpAbMap(From::D, To::C, M::fmpz_mat, Mi::fmpz_mat)
    r = new()
    function image(a::FinGenGrpAbElem{D})
      return FinGenGrpAbElemCreate(To, a.coeff*M)
    end

    function preimage(a::FinGenGrpAbElem{C})
      return FinGenGrpAbElemCreate(From, a.coeff*Mi)
    end

    r.header = MapHeader(From, To, image, preimage)
    return r
  end

end

function FinGenGrpAbMap{D <: FinGenGrpAb, C <: FinGenGrpAb}(From::D, To::C, M::fmpz_mat)
  return FinGenGrpAbMap{D, C}(From, To, M)
end

function FinGenGrpAbMap{D <: FinGenGrpAb, C <: FinGenGrpAb}(From::D, To::C, M::fmpz_mat, Mi::fmpz_mat)
  return FinGenGrpAbMap{D, C}(From, To, M, Mi)
end


text:
"""
    Genser.deserialize([to::Type,] serialized, M::Module; kwargs...)

Deserialize a value with an explicit deserialization module.

If `to` is not provided the value is deserialized into its default representation.

The keyword arguments are passed to the deserialization module.
"""
function deserialize(to::Type, value, M::Module; kwargs...)
    if !isdefined(M, :deserialize_genser)
        throw(ArgumentError("$M cannot deserialize data"))
    end

    rawdata = M.deserialize_genser(value; kwargs...)
    genserdata = togenser(gensertypefor(to), rawdata)

    fromgenser(to, genserdata)
end

function deserialize(to::Type, value, ::Type{M}; kwargs...) where M <: MIME
    Mod = genser_type_for_mime(M)
    deserialize(to, value, Mod::Module; kwargs...)
end

"""
    Genser.serialize([to::Type,] serialized, mime::MIME; kwargs...)

Deserialize a value using the registered deserialization module for a mime type.
"""
@inline function deserialize(to::Type, value, m::MIME; kwargs...)
    deserialize(to, value, typeof(m); kwargs...)
end

"""
    Genser.serialize([to::Type,] serialized, mime::AbstractString; kwargs...)

Deserialize a value using the registered deserialization module for a mime type provided as a string.
"""
@inline function deserialize(to::Type, value, m::AbstractString; kwargs...)
    deserialize(to, value, MIME(m); kwargs...)
end

@inline function deserialize(value, m; kwargs...)
    deserialize(Any, value, m; kwargs...)
end

ext: jl


text:
@doc doc"""
Minimise le problème : ``min_{||s||< \delta_{k}} q_k(s) = s^{t}g + (1/2)s^{t}Hs``
                        pour la ``k^{ème}`` itération de l'algorithme des régions de confiance

# Syntaxe
```julia
sk = Gradient_Conjugue_Tronque(fk,gradfk,hessfk,option)
```

# Entrées :
   * **gradfk**           : (Array{Float,1}) le gradient de la fonction f appliqué au point xk
   * **hessfk**           : (Array{Float,2}) la Hessienne de la fonction f appliqué au point xk
   * **options**          : (Array{Float,1})
      - **delta**    : le rayon de la région de confiance
      - **max_iter** : le nombre maximal d'iterations
      - **tol**      : la tolérance pour la condition d'arrêt sur le gradient

# Sorties:
   * **s** : (Array{Float,1}) le pas s qui approche la solution du problème : ``min_{||s||< \delta_{k}} q(s)``

# Exemple d'appel:
```julia
gradf(x)=[-400*x[1]*(x[2]-x[1]^2)-2*(1-x[1]) ; 200*(x[2]-x[1]^2)]
hessf(x)=[-400*(x[2]-3*x[1]^2)+2  -400*x[1];-400*x[1]  200]
xk = [1; 0]
options = []
s = Gradient_Conjugue_Tronque(gradf(xk),hessf(xk),options)
```
"""
function Gradient_Conjugue_Tronque(gradfk,hessfk,options)

    "# Si option est vide on initialise les 3 paramètres par défaut"
    if options == []
        deltak = 2
        max_iter = 100
        tol = 1e-6
    else
        deltak = options[1]
        max_iter = options[2]
        tol = options[3]
    end

    n = length(gradfk)
    s = zeros(n)
    q(s) = gradfk' * s + (1/2)*s'*hessfk*s

    sj = zeros(n)
    gj = gradfk
    pj = -gradfk

    for j = 1:max_iter
        kj = pj' * hessfk * pj
        if kj<0
            "on écrit ||sj+rho*pj||=deltak sous forme de ax^2 +bx+c=0"
            a = norm(pj)^2
            b = 2*sj'*pj
            c = norm(sj)^2 - deltak^2
            d = sqrt(b^2 - 4a*c)
            rho1, rho2 = (-b-d)/(2a), (-b+d)/(2a)

            if q(sj+rho1*pj)>q(sj+rho2*pj)
                rho = rho2
            else
                rho = rho1
            end # if
            s = sj+rho*pj
            break
        end # if

        alphaj = (gj' * gj)/kj

        "Cas où : ||sj+alpha*pj||>deltak"
        if norm(sj+alphaj*pj)>=deltak
            "on écrit ||sj+alpha*pj||=deltak sous forme de ax^2 +bx+c=0"
            a = norm(pj)^2
            b = 2*sj'*pj
            c = norm(sj)^2 - deltak^2

            d = sqrt(b^2 - 4a*c)
            alpha1, alpha2 = (-b-d)/(2a), (-b+d)/(2a)

            if alpha1>0
                alpha = alpha1
            elseif alpha2>0
                alpha = alpha2
                s = sj+alpha*pj
                break
            end # if
        end # if

        sj = sj+alphaj*pj
        gj_1 = gj+alphaj*hessfk*pj
        betaj = (gj_1'*gj_1)/(gj'*gj)
        pj = -gj_1 + betaj*pj

        if norm(gj_1)<=tol*norm(gradfk)
            s = sj
            break
        end # if
        gj = gj_1
    end # for

   return s
end

max_stars_repo_path: 2A/S7/UE Optimisation et RO/Optimisation/Optinum/src/Gradient_Conjugue_Tronque.jl

max_issues_repo_name: Sajid-Badr/ENSEEIHT

max_forks_repo_name: Sajid-Badr/ENSEEIHT


text:
data = readlines()

function is_nice(s)
    count(∈("aeiou"), s) ≥ 3 && any(c == cc for (c, cc) in zip(s[1:end-1], s[2:end])) && !any(occursin.(split("ab cd pq xy"), Ref(s)))
end

function is_nice_2(s)
    !(isnothing(match(r"(..).*\1", s)) || isnothing(match(r"(.).\1", s)))
end

function solve_1()
    is_nice.(data) |> count
end

function solve_2()
    is_nice_2.(data) |> count
end

solve_1() |> println
solve_2() |> println

max_stars_repo_name: pLOPeGG/JuliaOfCode

max_issues_repo_name: pLOPeGG/JuliaOfCode

max_forks_repo_name: pLOPeGG/JuliaOfCode


text:
function linear_step(hmm, observations, obs_lengths)
    O,T = size(observations);
    a = transpose(log.(hmm.a)); A = log.(hmm.A) #less expensive to transpose transmatrix in the backwards_sweep than to transpose here and take a ranged view ie view(A,m:m,:) is more expensive than transpose(view(A,m,:))
    N = length(hmm.B); Γ = length(hmm.B[1].support);
    mask=observations.!=0
    #INITIALIZATION
    βoi_T = zeros(O,N); βoi_t = zeros(O,N) #log betas at T initialised as zeros
    Eoγim_T = fill(-Inf,O,Γ,N,N); Eoγim_t = fill(-Inf,O,Γ,N,N)
    for m in 1:N, i in 1:N, γ in 1:Γ, o in 1:O
        observations[o, obs_lengths[o]] == γ && m == i && (Eoγim_T[o, γ, i, m] = 0)
    end
    Tijm_T = fill(-Inf,O,N,N,N); Tijm_t = fill(-Inf,O,N,N,N) #Ti,j(T,m) = 0 for all m; in logspace

    #RECURRENCE
    βoi_T,Tijm_T,Eoγim_T=backwards_sweep!(hmm,A,N,Γ,βoi_T,βoi_t,Tijm_T,Tijm_t,Eoγim_T,Eoγim_t,observations,mask,obs_lengths)

    #TERMINATION
    lls = c_llhs(hmm,observations[:,1])
    α1om = lps.(lls,a) #first position forward msgs
    log_pobs = [c_lse(lps.(α1om[o,:], βoi_T[o,:])) for o in 1:O]

    Toij = [c_lse([lps(Tijm_T[o,i,j,m], α1om[o,m], -log_pobs[o]) for m in 1:N]) for o in 1:O, i in 1:N, j in 1:N] #terminate Tijs with forward messages
    Eoiγ = [c_lse([lps(Eoγim_T[o,γ,i,m], α1om[o,m], -log_pobs[o]) for m in 1:N]) for o in 1:O, i in 1:N, γ in 1:Γ] #terminate Eids with forward messages

    #INTEGRATE ACROSS OBSERVATIONS AND SOLVE FOR NEW BHMM PARAMS
    #INITIAL STATE DIST
    a_o=α1om.+βoi_T.-c_lse.(eachrow(α1om.+βoi_T)) #estimate a for each o
    obs_penalty=log(O) #broadcast subtraction to normalise log prob vals by obs number
    new_a=c_lse.(eachcol(a_o)).-obs_penalty #sum over obs and normalise by number of obs
    #TRANSITION MATRIX
    new_A = [c_lse(view(Toij,:,i,j)) for i in 1:N, j in 1:N].-c_lse.(eachcol(c_lse.([view(Toij,o,i,:) for o in 1:O, i in 1:N])))
    #EMISSION MATRIX
    new_b=[c_lse(view(Eoiγ,:,i,γ)) for i in 1:N, γ in 1:Γ].-c_lse.(eachcol(c_lse.([view(Eoiγ,o,i,:) for o in 1:O, i in 1:N])))
    new_B=[Categorical(exp.(new_b[i,:])) for i in 1:N]

    return BHMM(exp.(new_a), exp.(new_A), new_B, hmm.partition), lps(log_pobs)
end
                #LINEAR_STEP SUBFUNCS
                function backwards_sweep!(hmm, A, N, Γ, βoi_T, βoi_t, Tijm_T, Tijm_t, Eoγim_T, Eoγim_t, observations, mask, obs_lengths)
                    for t in maximum(obs_lengths)-1:-1:1
                        last_β=copy(βoi_T)
                        lls = c_llhs(hmm,observations[:,t+1])
                        omask = findall(mask[:,t+1])
                        βoi_T[omask,:] .+= view(lls,omask,:)
                        Threads.@threads for m in 1:N
                            βoi_t[omask,m] = c_lse.(eachrow(view(βoi_T,omask,:).+transpose(view(A,m,:))))
                            for j in 1:N, i in 1:N
                                Tijm_t[omask, i, j, m] .= c_lse.(eachrow(lps.(view(Tijm_T,omask,i,j,:), view(lls,omask,:), transpose(view(A,m,:)))))
                                i==m && (Tijm_t[omask, i, j, m] .= logaddexp.(Tijm_t[omask, i, j, m], lps.(last_β[omask,j], A[m,j],lls[omask,j])))
                            end
                            for i in 1:N, γ in 1:Γ
                                Eoγim_t[omask, γ, i, m] .= c_lse.(eachrow(lps.(view(Eoγim_T,omask,γ,i,:),view(lls,omask,:),transpose(view(A,m,:)))))
                                if i==m
                                    symmask = findall(observations[:,t].==γ)
                                    Eoγim_t[symmask, γ, i, m] .= logaddexp.(Eoγim_t[symmask, γ, i, m], βoi_t[symmask,m])
                                end
                            end
                        end
                        βoi_T=copy(βoi_t); Tijm_T=copy(Tijm_t); Eoγim_T = copy(Eoγim_t);
                    end
                    return βoi_T, Tijm_T, Eoγim_T
                end

                #logsumexp with single dispatch for speed
                function c_lse(X::AbstractArray{T}; dims=:) where {T<:Real}
                    u=maximum(X)
                    u isa AbstractArray || isfinite(u) || return float(u)
                    return u + log(sum(x -> exp(x-u), X))
                end

                #log likelihoods
                function c_llhs(hmm, observation)
                    lls = zeros(length(observation),length(hmm.B))
                    Threads.@threads for d in 1:length(hmm.B)
                        lls[:,d] = logpdf.(hmm.B[d], observation)
                    end
                    return lls
                end


text:
################################################################################
# setting up
################################################################################
using SortingLab, SortingAlgorithms, BenchmarkTools, Base.Threads
rangelen = 1_000_000;
minval = 1;
srand(1);
a = rand(UInt32(minval):UInt32(rangelen+minval-1), 100_000_000);
lo = 1;
hi = length(a);
RADIX_SIZE = 10;
include("src/fsortperm_dev_code.jl");

################################################################################
# Hybrid LSD and counting sort
# I have shown below that coutning sorts beats radix sort up to 2^11= 2048 unique values
# so it makes sense to create a better one
################################################################################

a[fsortperm_msd_hybrid(a, rangelen, minval, 10, sortcutoff = 2048)] |> issorted
a[fsortperm_msd_hybrid(a, rangelen, minval, 10, sortcutoff = -1)] |> issorted

@benchmark fsortperm_msd_hybrid($a, rangelen, minval, 10, sortcutoff = -1) samples = 5 seconds = 120
@benchmark fsortperm_msd_hybrid($a, rangelen, minval, 10, sortcutoff = 2048) samples = 5 seconds = 120

################################################################################
# Straight counting sort
################################################################################

if false
    function counting_vs_radix(n)
        a = rand(1:n, 100_000_000)
        (@belapsed(sortperm_int_range2($a, $n, 1)), @belapsed(sortperm_int_range_p1($a, $n, 1, 10)))
    end

    tic()
    res = [counting_vs_radix(n) for n in 2.^(2:20)]
    toc()

    using Plots
    plot([r[1] for r in res], title="Counting-sortperm vs Radix-sortperm (# of elements = 100m)", label = "Counting-sortperm",
        ylabel = "seconds", xlabel = "log_2(n) unique values", ylim = (0, 15))
    plot!([r[2] for r in res], label = "Radix-sortperm")
    plot!(xticks = (2:20))
    savefig("Counting-sortperm vs Radix-sortperm ( of elements = 100m).png")

    a = rand(1:1_000_000, 100_000_000)
    @time sortperm_int_range2(a, 1_000_000, 1)
    @code_warntype sortperm_int_range2(a, 1_000_000, 1)
    # @code_lowered sortperm_int_range2(a, 1_000_000, 1)
    # @code_typed sortperm_int_range2(a, 1_000_000, 1)
    # @code_llvm sortperm_int_range2(a, 1_000_000, 1)
    # @code_native sortperm_int_range2(a, 1_000_000, 1)

    @benchmark sortperm_int_range2($a, $rangelen, $minval)
end

################################################################################
# MSD sortperm
################################################################################
function fsortperm_msd(a, rangelen, minval, RADIX_SIZE;  sortcutoff = -1)
    # println(RADIX_SIZE," ", now())
    @assert 2^32 > rangelen
    @assert 2^32 >= length(a)
    vs = Vector{VT}(length(a))
    @inbounds for i in eachindex(a)
        vs[i] = Pair(UInt32(i),UInt32(a[i]))
    end
    # @time @inbounds vs = [Pair(UInt32(i),UInt32(ai)) for (i, ai) in enumerate(a)]

    iters = Int(ceil(log(2, rangelen)/RADIX_SIZE))
    @assert iters <= ceil(Integer, sizeof(typeof(vs[1].second))*8/RADIX_SIZE)

    _fsortperm_msd(vs, 1, length(vs), rangelen, minval, RADIX_SIZE, iters, sortcutoff = sortcutoff)
end

function _fsortperm_msd(vs, lo, hi, rangelen, minval, RADIX_SIZE, iters, ts = similar(vs); sortcutoff = -1)
    len = hi-lo+1
    if len  == 1
        return [Int(vs[lo].first)]
    elseif len <= sortcutoff
        vsv = vs[lo:hi]
        sort!(vsv, by = x->x.second)
        return [Int(vsi.first) for vsi in vsv]
    end

    RADIX_MASK = UInt32(1<<RADIX_SIZE-1)

    # Init
    bin = zeros(UInt32, 2^RADIX_SIZE)

    # Histogram for each element, radix
    @inbounds for i = lo:hi
        v = vs[i].second
        idx = Int((v >> (iters-1)*RADIX_SIZE) & RADIX_MASK) + 1
        @inbounds bin[idx] += 1
    end

    # Sort!

    # Unroll first data iteration, check for degenerate case
    v = vs[hi].second
    idx = Int((v >> (iters-1)*RADIX_SIZE) & RADIX_MASK) + 1

    # are all values the same at this radix?
    # if bin[idx,iters] == len;  continue;  end
    if bin[idx] == len
        if iters == 1
            return [Int(vsi.first) for vsi in vs[lo:hi]]
        else
            return _fsortperm_msd(vs, lo, hi, rangelen, minval, RADIX_SIZE, iters-1, ts)
        end
    end

    # cbin = cumsum(bin[:])
    cbin  = copy(bin)
    cumsum!(bin, bin)
    cumsumbin = copy(bin)

    # Finish the loop...
    @inbounds for i in lo:hi
        ts[i] = vs[i]
    end

    @inbounds for i in lo:hi
        v = ts[i].second
        idx = Int((v >> (iters-1)*RADIX_SIZE) & RADIX_MASK) + 1
        ci = bin[idx]
        vs[lo - 1 + ci] = ts[i]
        bin[idx] -= 1
    end

    if iters == 1
        return [Int(vsi.first) for vsi in vs[lo:hi]]
    end

    # now that it's finished sorting sort each chunk recursive now
    res = Vector{Int}(length(vs))
    start_now = 1
    for i = 1:length(cumsumbin)
        if cbin[i] != 0
            res[start_now:cumsumbin[i]] = _fsortperm_msd(vs, start_now, cumsumbin[i], rangelen, minval, RADIX_SIZE, iters-1, ts)
            start_now = cumsumbin[i] + 1
        end
    end
    res
end

@time fsortperm_msd(a, rangelen, minval, 10)
@time fsortperm_msd(a, rangelen, minval, 10, sortcutoff = 64)

using RCall
@rput a
r = R"""
memory.limit(2^31-1)
system.time(order(a))
"""

srand(1);
a = rand((minval):(rangelen+minval-1), 1_000_000_000);
@benchmark fsortperm_msd($a, rangelen, minval, 10) samples = 5 seconds = 120
@benchmark fsortperm_msd($a, rangelen, minval, 10, sortcutoff = 64) samples = 5 seconds = 120
@benchmark fsortperm_msd($a, rangelen, minval, 10, sortcutoff = 2048)

srand(1)
a = rand(UInt32(minval):UInt32(rangelen+minval-1), 2^16);
vs = [Pair(UInt32(i),UInt32(ai)) for (i,ai) in enumerate(a)];
@btime fsortperm_msd($a, rangelen, minval, 10);
@btime sort($vs, by=x->x.second);

function test(n)
    srand(1);
    c = rand(UInt32(1):UInt32(1_000_000), n)
    vs = Vector{VT}(length(c))
    @inbounds for i in eachindex(c)
        vs[i] = Pair(UInt32(i),UInt32(c[i]))
    end

    res = [
    @belapsed(fsortperm_msd($c, 1_000_000, 1, 10, sortcutoff = -1)),
    @belapsed(sort($vs, by=x->x.second, alg = InsertionSort)),
    @belapsed(sort($vs, by=x->x.second))]
    res
end

hehe = test.(2.^(1:10))

using Plots
plot([h[1] for h in hehe])
plot!([h[2] for h in hehe])
plot!([h[3] for h in hehe])

@benchmark fsortperm_msd($a, rangelen, minval, 10) samples = 5 seconds = 30
@benchmark fsortperm_msd($a, rangelen, minval, 10, sortcutoff = 64) samples = 5 seconds = 30
@benchmark fsortperm_msd($a, rangelen, minval, 10, sortcutoff = 128) samples = 5 seconds = 30
@benchmark fsortperm_msd($a, rangelen, minval, 10, sortcutoff = 256) samples = 5 seconds = 30
@benchmark fsortperm_msd($a, rangelen, minval, 10, sortcutoff = 16) samples = 5 seconds = 30
@benchmark fsortperm_msd($a, rangelen, minval, 10, sortcutoff = 512) samples = 5 seconds = 30
@benchmark fsortperm_msd($a, rangelen, minval, 10, sortcutoff = 1024) samples = 5 seconds = 30
@benchmark fsortperm_msd($a, rangelen, minval, 10, sortcutoff = 2048) samples = 5 seconds = 30

@benchmark fsortperm_msd($a, rangelen, minval, 10) samples = 5 seconds = 30
@benchmark fsortperm_msd($a, rangelen, minval, 11) samples = 5 seconds = 30

b = a[as];
issorted(b)

################################################################################
# sorting algorithms from discourse
# https://discourse.julialang.org/t/ironic-observation-about-sort-and-sortperm-speed-for-small-intergers-vs-r/8715/20?u=xiaodai
################################################################################
# a = rand(1:1_000_000, 100_000_000)

@time sortperm_int_range_p1(a[1:1], rangelen, 1, 11);
gc()
@time sortperm_int_range_p1(a, rangelen, 1, 10);
gc()
@time sortperm_int_range_p1(a, rangelen, 1, 11);

@benchmark sortperm_int_range_p1(a, rangelen, 1, 10)

aint = Int.(a)
@benchmark sortperm_int_range_p1(aint, rangelen, 1, 10)
issorted(a[aa])

using RCall
@rput a

R"""
system.time(order(a))
"""

function bathtub(rangelen, minval, label, ontop)
    a = rand(UInt32(minval):UInt32(rangelen+minval-1), 100_000_000)
    res = [@elapsed sortperm_int_range_p1(a, rangelen, minval, i) for i = 1:20]
    if ontop
        plot!(res, label=label)
    else
        plot(res, label=label)
    end
    res
end

using Plots
gc()
res1m = bathtub(1_000_000, 1, "1m", false)
res10k = bathtub(10_000, 1, "10k", true)
res2k = bathtub(2_000, 1, "2k", true)

df = DataFrame(res1m = res1m, res10k = res10k, res2k = res2k)
using CSV
# CSV.write("sortperm_int_range_p1 1 billion uniques (2k-1m).csv", df)

plot(res1m, label="1m", title = "sortperm_int_range_p1 1 billion uniques (2k-1m)")
plot!(res10k, label="10k")
plot!(res2k, label="2k")
# savefig("sortperm_int_range_p1 1 billion uniques (2k-1m).png")

################################################################################
# MSD sortperm in place
################################################################################
function fsortperm_msd_in_place(a, rangelen, minval, RADIX_SIZE)
    # println(RADIX_SIZE," ", now())
    @assert 2^32 > rangelen
    @assert 2^32 >= length(a)
    vs = Vector{VT}(length(a))
    @inbounds for i in eachindex(a)
        vs[i] = Pair(UInt32(i),UInt32(a[i]))
    end
    # @time @inbounds vs = [Pair(UInt32(i),UInt32(ai)) for (i, ai) in enumerate(a)]
    _fsortperm_msd_in_place(vs, rangelen, minval, RADIX_SIZE)
end

function _fsortperm_msd_in_place(vs, rangelen, minval, RADIX_SIZE)
    RADIX_MASK = UInt32(1<<RADIX_SIZE-1)

    # Init
    lo = 1
    hi = length(vs)
    iters = Int(ceil(log(2, rangelen)/RADIX_SIZE))

    # println("iters: $iters")
    @assert iters <= ceil(Integer, sizeof(typeof(vs[1].second))*8/RADIX_SIZE)
    bin = zeros(UInt32, 2^RADIX_SIZE, iters)

    # Histogram for each element, radix
    @inbounds for i = lo:hi
        v = vs[i].second
        for j = 1:iters
            idx = Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK) + 1
            @inbounds bin[idx,j] += 1
        end
    end

    # Sort!
    len = hi-lo+1
    done = BitArray(len)
    @inbounds for j = iters:iters
        done .= false
        # Unroll first data iteration, check for degenerate case
        v = vs[1].second
        idx = Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK) + 1

        # are all values the same at this radix?
        if bin[idx,j] == len;  continue;  end

        cbin = cumsum(bin[:,j])

        # ci = cbin[idx]
        # ts[ci] = vs[hi]
        # cbin[idx] -= 1

        # Finish the loop...
        while !all(done)
            @inbounds for i in lo:hi
                if !done[i]
                    v = vs[i].second
                    idx = Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK) + 1
                    ci = cbin[idx]
                    vs[ci], vs[i] = vs[i], vs[ci]
                    # after the above swap vs[ci] is where is wants to be
                    done[ci] = true
                    cbin[idx] -= 1
                end
            end
            # println(sum(done))
        end
    end

    res = Vector{Int}(length(vs))
    @inbounds for i in eachindex(vs)
        res[i] = Int(vs[i].first)
    end
    res
    # [Int(vs1.first) for vs1 in vs]
end

@time fsortperm_msd_in_place(a, rangelen, minval, 11)

@code_warntype fsortperm_msd_in_place(a, rangelen, minval, 10)

################################################################################
# Staged counting sort
# this is a modified counting sort
################################################################################
function sortperm_int_range_p8(a, rangelen, minval, RADIX_SIZE)
    @assert 2^32 > rangelen
    @assert 2^32 >= length(a)
    vs::VTV = Vector{VT}(length(a))
    # vs = Vector{VT}(length(a))
    @inbounds for i in eachindex(a)
        vs[i] = Pair(UInt32(i),UInt32(a[i]))
    end
    _sortperm_int_range_p8(vs, rangelen, minval, RADIX_SIZE)
end

function _sortperm_int_range_p8(vs, rangelen, minval, RADIX_SIZE)
    cnts = fill(0, rangelen)

    # create a count first
    @inbounds for ai in vs
        cnts[ai.second - minval + 1] +=  1
    end

    # create cumulative sum
    cumsum!(cnts, cnts)

    # note: this is for performance reasons only
    # do not change unless it's been benchmarked
    # creat a smaller sum
    iter, rem_iter = divrem(log(2,rangelen), RADIX_SIZE)
    step = Int(ceil(2^(rem_iter+RADIX_SIZE*(iter-1))))
    la = length(vs)
    res = Vector{Int}(la)

    for iter1 in 1:iter
        range = step:step:rangelen
        cnts_smaller = Vector{UInt32}(1 << RADIX_SIZE)
        for (i, rng) in enumerate(range)
            cnts_smaller[i] = cnts[rng]
        end
        last_index = Int(ceil(rangelen/step))
        cnts_smaller[last_index] = cnts[rangelen]

        # now distribute the cnts smaller according to the below cnts_smaller
        @inbounds for i in la:-1:1
            ai = Int(ceil((vs[i].second - minval + 1)/step))
            c = cnts_smaller[ai]
            cnts_smaller[ai] -= 1
            res[c] = i
        end
    end

    res
end

@time a1 = sortperm_int_range_p8(a, rangelen, minval, RADIX_SIZE)

(c1 |> diff .|> abs |> mean) |> round |> Int # 33 330 788 #33 333 154

# const RADIX_SIZE = 11
# if RADIX_SIZE < 9
#     const RADIX_MASK = UInt32(2<<(RADIX_SIZE-1)-1)
# elseif RADIX_SIZE < 17
#     const RADIX_MASK = UInt32(2<<(RADIX_SIZE-1)-1)
# else
#     const RADIX_MASK = UInt32(2<<(RADIX_SIZE-1)-1)
# end

################################################################################
# Straight counting sort mean travel distance
################################################################################
rangelen = 2000
minval = 1
a = rand(UInt32(minval):UInt32(rangelen+minval-1), 100_000_000)

function sortperm_int_range2c(a, rangelen, minval)
    cnts = fill(0, rangelen)

    # create a count first
    @inbounds for ai in a
        cnts[ai - minval + 1] +=  1
    end

    # create cumulative sum
    cumsum!(cnts, cnts)

    la = length(a)
    res = zeros(Int, la)
    crecord = Vector{Int}(la)
    @inbounds for i in la:-1:1
        ai = a[i] - minval + 1
        c = cnts[ai]
        cnts[ai] -= 1
        crecord[la-i+1] = c
        res[c] = i
    end

    res, crecord
end

@time a2, c2 = sortperm_int_range2c(a, rangelen, minval)

(c2 |> diff .|> abs |> mean) |> round |> Int # 33 328 745 # 33 333 180

################################################################################
# plotting the results
################################################################################
fsortperm_bench = @benchmark fsortperm(a)
using Plots
bar(
    ["Base.sortperm_int_range", "fsortperm", "R order(...)"],
    [17.218, fsortperm_bench.times./1e9 |> mean,  2.17 ],
    title = "`sortperm_int_range` perf (values = 1:1m, vector len = 100m)",
    ylabel="seconds"
)
savefig("julia_vs_r_sortperm_integer.png")

################################################################################
# proved that the res[c] is the problem statement in straight counting sort
# yes it's proven as I just removed that line but still sum counts and it completed
# very quickly
################################################################################
function prove_scs(a, rangelen, minval)
    cnts = fill(0, rangelen)

    # create a count first
    @inbounds for ai in a
        cnts[ai - minval + 1] +=  1
    end

    # create cumulative sum
    cumsum!(cnts, cnts)

    la = length(a)
    # res = Vector{Int}(la)
    sumc = 0
    @inbounds for i in la:-1:1
        ai = a[i] - minval + 1
        c = cnts[ai]
        cnts[ai] -= 1
        # res[c] = i
        sumc += c
    end

    sumc
end

@benchmark prove_scs(a, rangelen, minval)
# BenchmarkTools.Trial:
#   memory estimate:  7.63 MiB
#   allocs estimate:  3
#   --------------
#   minimum time:     1.802 s (0.00% GC)
#   median time:      1.819 s (0.00% GC)
#   mean time:        1.879 s (0.03% GC)
#   maximum time:     2.018 s (0.00% GC)
#   --------------
#   samples:          3
#   evals/sample:     1

################################################################################
# yet to be classified
################################################################################
function sortperm_int_range7(a::Vector{T}, rangelen, minval) where T
    cnts = fill(0, rangelen)

    # create a count first
    @inbounds for ai in a
        cnts[ai - minval + 1] +=  1
    end

    # create cumulative sum
    cumsum!(cnts, cnts)

    la = length(a)
    res = Vector{Int}(la)

    upto = zero(UInt8)
    cache = zeros(Int, 256)

    jumptimes = 0
    for i in la:-1:1
        ai = a[i] - minval + 1
        c = cnts[ai]
        cnts[ai] -= 1

        upto += UInt8(1)
        # print(upto)
        if upto != 0
            cache[upto] = c
        else
            cache[256] = c
            res[cache] .= collect(la-jumptimes*256:-1:la-(jumptimes+1)*256+1)
        end
    end

    res
end
@time sortperm_int_range7(a, rangelen, minval)

function sortperm_int_range6(a::Vector{T}, rangelen, minval) where T
    cnts = fill(0, rangelen)

    # create a count first
    @inbounds for ai in a
        cnts[ai - minval + 1] +=  1
    end

    # create cumulative sum
    cumsum!(cnts, cnts)

    la = length(a)
    res = Vector{Int}(la)

    upto = zero(UInt8)
    cache = zeros(Int, 256, 2)

    for i in la:-1:1
        ai = a[i] - minval + 1
        c = cnts[ai]
        cnts[ai] -= 1

        upto += UInt8(1)
        # print(upto)
        if upto != 0
            cache[upto, 1], cache[upto, 2] = c, i
        else
            cache[256, 1], cache[256,2] = c,i
            res[cache[:,1]] .= cache[:,2]
        end
    end

    res
end
@time sortperm_int_range6(a, rangelen, minval)

function sortperm_int_range5(a, rangelen, minval)
    cnts = fill(0, rangelen)

    # create a count first
    @inbounds for ai in a
        cnts[ai - minval + 1] +=  1
    end

    # create cumulative sum
    cumsum!(cnts, cnts)

    la = length(a)
    res = Vector{Int}(la)

    upto = zeros(UInt8, rangelen)

    indexes = Array{Int, 2}(256, rangelen)

    @inbounds for i in 1:la
        ai = a[i] - minval + 1
        # c = cnts[ai]
        # cnts[ai] -= 1
        u = upto[ai]
        upto[ai] += 1
        indexes[u, ai] = i
    end

    indexes
end

@time sortperm_int_range5(a, rangelen, minval)

@time sortperm(a);
@time Base.Sort.sortperm_int_range(a, 999_999, 1);

function sortperm_int_range2(a, rangelen, minval)
    cnts = fill(0, rangelen)

    # create a count first
    @inbounds for ai in a
        cnts[ai - minval + 1] +=  1
    end

    # create cumulative sum
    cumsum!(cnts, cnts)

    la = length(a)
    res = Vector{Int}(la)
    @inbounds for i in la:-1:1
        ai = a[i] - minval + 1
        c = cnts[ai]
        cnts[ai] -= 1
        res[c] = i
    end

    res
end

rangelen = 1_000_000
minval = 1
@time x = sortperm_int_range2(a, 1_000_000, 1);
@code_warntype sortperm_int_range2(a, 1_000_000, 1);

reinterpret(Vector{UInt16}, 10)

a[x]
issorted(a[x])
using SortingLab
i = collect(1:length(a))
@time SortingLab.sorttwo!(a, i)

function sortperm_int_range3(a)
    ca = copy(a)
    i = collect(1:length(ca))
    SortingLab.sorttwo!(ca, i)
    i
end

@time ii = sortperm_int_range3(a)
issorted(a[ii])

using BenchmarkTools
@benchmark sortperm_int_range2(a, 1_000_000, 1) samples = 5 seconds = 120
# BenchmarkTools.Trial:
#   memory estimate:  770.57 MiB
#   allocs estimate:  4
#   --------------
#   minimum time:     15.931 s (0.00% GC)
#   median time:      16.970 s (0.00% GC)
#   mean time:        17.218 s (0.00% GC)
#   maximum time:     18.552 s (0.00% GC)
#   --------------  samples:          5
#   evals/sample:     1
@benchmark sortperm_int_range3(a) samples = 5 seconds = 120
# BenchmarkTools.Trial:
#   memory estimate:  2.99 GiB
#   allocs estimate:  557
#   --------------
#   minimum time:     10.916 s (0.00% GC)
#   median time:      13.304 s (3.02% GC)
#   mean time:        12.914 s (2.89% GC)
#   maximum time:     14.062 s (6.70% GC)
#   --------------
#   samples:          5
#   evals/sample:     1

function sortperm_int_range4(a, minval)
    ula = UInt32(length(a))
    lz = leading_zeros(ula)

    ca = Vector{UInt64}(ula)

    tz = trailing_zeros(UInt64(minval))

    # ltz = lz - tz
    for i = UInt32(1):ula
        ca[i] = (UInt64(a[i]) << 32) | i
    end

    SortingLab.sort32!(ca, skipbits = 32+tz)

    Int.(ca .& (UInt64(2)^(32-lz)-1))
end

using SortingAlgorithms, BenchmarkTools, SortingLab

const a = rand(1:1_000_000, 100_000_000)
function test1(a)
    ca = (UInt32.(a) .<< 32) .| collect(UInt(1):UInt(length(a)))
    SortingLab.sort32!(ca)
end
@time test1(a)

@time x=sortperm_int_range4(a, 1);
@code_warntype test1(a);
issorted(a[x])

@benchmark sortperm_int_range4($a,1) samples = 5 seconds = 120
# BenchmarkTools.Trial:
#   memory estimate:  2.24 GiB
#   allocs estimate:  20
#   --------------
#   minimum time:     6.899 s (1.16% GC)
#   median time:      7.528 s (4.39% GC)
#   mean time:        7.508 s (3.58% GC)
#   maximum time:     8.004 s (4.13% GC)
#   --------------
#   samples:          5
#   evals/sample:     1

using StatPlots

bar(
    ["sortperm_int_range", "my sortperm_int_range4", "R order(...)"],
    [17.218, 9.577, 2.17 ],
    title = "`sortperm_int_range` perf (values = 1:1m, vector len = 100m)",
    ylabel="seconds"
)
savefig("julia_vs_r_sortperm_integer.png")

@time aa = sortperm_int_range4(a, 1_000_000, 1)
issorted(a[aa])

a = rand(1:1_000_000, 100_000_000)

js = @elapsed sort(a)
jsp = @elapsed sortperm(a)

using RCall
@rput a
r = R"""
list(system.time(sort(a)),  system.time(order(a)))
"""

using StatPlots

groupedbar(
    ["sort", "sort", "sortperm", "sortperm"],
    [js, r[1][3], jsp, r[2][3]],
    group = ["Julia","R","Julia","R"],
    title = "Julia vs R sort and perm (values = 1:1m, vector len = 100m)",
    ylabel="seconds"
)
savefig("julia_vs_r_sortandperm_integer.png")

Got it!!! Almost there. Thanks to @sdanisch for this discussion. I have simply repurposed the code. Basically, I am using sort to create fsortperm which on par with R now! I am sure my code is pretty bad and that’s why it doesn’t beat R!!

@tshort, in a sense you are right, using sort seems to be the key, but basically, I have used Julia underlying counting sort to get this sort of performance. I still don’t understand why my original version is so slow as it’s just an application of counting sort (one pass radix sort).

julia_vs_r_sortperm_integerjulia_vs_r_sortperm_integer.png3600x2400 66.9 KB
See code below

# function fsortperm(a::Integer)
#     A2 = transpose(hcat(a, collect(1:length(a))))
#     Atup = reinterpret(NTuple{2, Int}, A2, (size(A2, 2),));
#     sort!(Atup)
#     [atup[2] for atup in Atup]
# end

# @time i = fsortperm(a)
# issorted(a[i]) # true

# fsortperm_bench = @benchmark fsortperm(a)
# using Plots
# bar(
#     ["Base.sortperm_int_range", "fsortperm", "R order(...)"],
#     [17.218, fsortperm_bench.times./1e9 |> mean,  2.17 ],
#     title = "`sortperm_int_range` perf (values = 1:1m, vector len = 100m)",
#     ylabel="seconds"
# )
# savefig("julia_vs_r_sortperm_integer.png")

@benchmark sort($a)
a32 = Int32.(a)

@benchmark sort($a32)

A2 = transpose(hcat(a, collect(1:length(a))))

ai = hcat(a, collect(1:length(a)))
function tshort1(ai, rangelen, minval)
    cnts = fill(0, rangelen)

    # create a count first
    @inbounds for i in 1:size(ai,1)
        cnts[ai[i,1] - minval + 1] +=  1
    end

    # create cumulative sum
    cumsum!(cnts, cnts)

    la = size(ai,1)
    res = Vector{Int}(la)
    aic = similar(ai)
    @inbounds for i in la:-1:1
        ai1 = ai[i,1] - minval + 1
        c = cnts[ai1]
        cnts[ai1] -= 1
        aic[c,:] .= ai[i,:]
    end

    aic
end

################################################################################
# failed ideas
################################################################################

@time tshort1(ai, 1_000_000, 1)
@code_warntype tshort1(ai, 1_000_000, 1)

Atup = reinterpret(NTuple{2, Int}, A2, (size(A2, 2),));

@time Atup_sorted = sort(Atup, by = x->x[1]);
@time Atup_sorted = sort!(Atup, by = x->x[1], alg=RadixSort)

sort(Atup; alg=QuickSort)

@benchmark sort(Atup)
@benchmark sort(Atup; alg=QuickSort)
@benchmark sort(Atup; alg=RadixSort)

# https://discourse.julialang.org/t/faster-sorting-with-gpu/8409/22?u=xiaodai
using SortingAlgorithms
function fsortperm(a::Vector{Int})
    @time A2 = transpose(hcat(a, collect(1:length(a))))
    @time Atup = reinterpret(NTuple{2, Int}, A2, (size(A2, 2),));
    @time sort!(Atup, by = x->x[1], alg=RadixSort)
    @time [atup[2] for atup in Atup]
end

@time fsortperm(a)

a = rand(1:1_000_000, 100_000_000)

using DataFrames
df = DataFrame(a = a, i = collect(1:length(a)))
@time abc(df) = sort!(df, cols = :a)
@time abc(df)

@time i = fsortperm(a)
issorted(a[i])[1]

ra = rand(Int, 100_000_000);
@time ira = fsortperm(ra);
issorted(ra[ira]) # 29.0015

@time sortperm(ra); # 51


text:
using IntervalOptimisation
using Test

include("sorted_vector.jl")
include("heaped_vector.jl")
include("optimise.jl")


text:
#=
string:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-08
=#

module String
    # valuables
    str1 = "Hello Julia"
    str2 = "こんにちは、Julia"

    # methods
    function get_str1()
        return str1
    end

    function get_str1_char(idx)
        return str1[idx]
    end

    function connect_str(str1, str2)
        println("How to connect str1:$str1 and str2:$str2")
        println("Using string(): ", string(str1, " ", str2))
        println("Using *: ", str1 * " " * str2)
        println("Using interpolation: ", "$str1 $str2")
    end

    function get_str2()
        return str2
    end

    function get_str2_unicode(idx)
        return str2[idx]
    end

    function main()
        # normal string and char
        println("str1 = ", get_str1())
        println("str1[1] = ", get_str1_char(1))
        println("typeof(str1[1]) = ", typeof(get_str1_char(1)))
        println("str1[end] = ", get_str1()[end])
        println("str1[1:5] = ", get_str1_char(1:5))
        # unicode string
        println("str2 = ", get_str2())
        # almost Hiragana and Kanji are represented as 3 bytes.
        println("str2[1] = ", get_str2_unicode(1))
        println("str2[4] = ", get_str2_unicode(4))
        # get next index
        println("nextind(str2, 1) = ", nextind(get_str2(), 1))
        println("nextind(str2, 4) = ", nextind(get_str2(), 4))
        # convert string array
        println("Vector{Char}(str2) = ", Vector{Char}(get_str2()))
        println("Vector{Char}(str2)[1] = ", Vector{Char}(get_str2())[1])
        println("Vector{Char}(str2)[2] = ", Vector{Char}(get_str2())[2])
        # string function
        println("length(Julia) = ", length("Julia"))
        println("repeat(Julia, 3) = ", repeat("Julia", 3))
        println("replace(Python is the best!!, Python => Julia) = ", replace("Python is the best!!", "Python" => "Julia"))
        println("split(Julia-lang, -) = ", split("Julia-lang", "-"))
        println("startswith(JuliaLang, Julia) = ", startswith("JuliaLang", "Julia"))
        println("endswith(JuliaLang, Julia) = ", endswith("JuliaLang", "Julia"))
        println("join([Julia, Lang], -) = ", join(["Julia", "Lang"], "-"))
        println("findfirst(Julia, JuliaLang) = ", findfirst("Julia", "JuliaLang"))
        # regular expression
        println("regex = ", r"J.*g")
        println("match(regex, JuliaLang is the best) = ", match(r"J.*g", "JuliaLang is the best").match)
        println("offset(regex, JuliaLang is the best) = ", match(r"J.*g", "JuliaLang is the best").offset)
        connect_str("Goodbye", "Julia")
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .String
    String.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
struct JointToLinkPosition{Rn,R3,S,JC} <: TaskMap{Rn,R3,S}
    state_cache::StateCache{S,JC}
    link::RigidBody{S}
end
struct JointToLinkPositionT{Rn,R3,S,JC} <: TaskMapT{Rn,R3,S}
    base_map::JointToLinkPosition{Rn,R3,S,JC}
end
JointToLinkPosition{ℝ{n},R3}(s::StateCache{S,JC}, l::RigidBody{S}) where {n,R3,S,JC} =
    JointToLinkPosition{ℝ{n},R3,S,JC}(s,l)

function task_map_emb(::EmbRep, ::EmbRep, xme, task_map::JointToLinkPosition{ℝ{n},R3}) where n
    state = task_map.state_cache[eltype(xme)]
    set_configuration!(state, xme)
    translation(transform_to_root(state, task_map.link))
end


text:
#for running hand-constructed test cases

using MetaGen
using JSON

#try to make objects_observed::Array{Array{Array{Array{Detection2D}}}} of observed objects.
#outer array is for scenes, then frames, the receptive fields, then last is an array of detections

################################################################################
params = Video_Params()

receptive_fields = make_receptive_fields()
num_receptive_fields = length(receptive_fields)

num_videos = 30
num_frames = 10

objects_observed = Matrix{Array{Detection2D}}(undef, num_videos, num_frames)
camera_trajectories = Matrix{Camera_Params}(undef, num_videos, num_frames)
for v=1:num_videos
	for f = 1

		labels = [2, 5]
		xs = [50.1, 225.5]
		ys = [60.2, 173.1]
		temp = Array{Detection2D}(undef, length(labels))
		for i = 1:length(labels)
			label = labels[i]
			x = xs[i]
			y = ys[i]
			temp[i] = (x, y, label) #so now I have an array of detections
		end
		#turn that array of detections into an array of an array of detections sorted by receptive_field
		#temp_sorted_into_rfs = map(rf -> filter(p -> within(p, rf), temp), receptive_fields)
		objects_observed[v, f] = temp

		c = Camera_Params(camera_location = Coordinate(0.01,0.02,0.001), camera_focus = Coordinate(1.0,1.0,1.0))
		camera_trajectories[v, f] = c
	end

	#################################
	#frame 2
	for f = 2:10

		labels = [2]
		xs = [50.1]
		ys = [60.2]
		temp = Array{Detection2D}(undef, length(labels))
		for i = 1:length(labels)
			label = labels[i]
			x = xs[i]
			y = ys[i]
			temp[i] = (x, y, label) #so now I have an array of detections
		end
		#turn that array of detections into an array of an array of detections sorted by receptive_field
		#temp_sorted_into_rfs = map(rf -> filter(p -> within(p, rf), temp), receptive_fields)
		#objects_observed[v, f] = temp_sorted_into_rfs
		objects_observed[v, f] = temp

		c = Camera_Params(camera_location = Coordinate(0.01,0.02,0.001), camera_focus = Coordinate(1.0,1.0,1.0))
		camera_trajectories[v, f] = c
	end
end

#count_observations(objects_observed)

outfile = string("test_case.csv")
file = open(outfile, "w")

#file header
print(file, "gt_V & ")
for v=1:29
	print(file, "avg V ", v, " & ")
	print(file, "dictionary realities PF for scene ", v, " & ")
	print(file, "mode realities PF for scene ", v, " & ")
end
print(file, "avg V 30 &")
print(file, "dictionary realities PF for scene 30 & ")
print(file, "mode realities PF for scene 30")

print(file, "\n")

gt_v = [0.5 0.0; 0.0 0.0; 0.5 0.0; 0.5 0.5; 0.5 0.1]
print(file, gt_v, " & ")
num_particles = 1

traces = unfold_particle_filter(num_particles, objects_observed, camera_trajectories, file)
println("done")
#visualize_observations(objects_observed, 1, 1, receptive_fields)
#visualize_trace(traces, 1, camera_trajectories, 1, 1, params)

# for v=1:4
# 	print_Vs_and_Rs_to_file(file, traces, num_particles, params, v)
# end
# print_Vs_and_Rs_to_file(file, traces, num_particles, params, 5, true)
#print_Vs_and_Rs_to_file(file, traces, num_particles, params, 1, true)
close(file)

max_stars_repo_name: zhangir-azerbayev/MetaGen

max_issues_repo_name: zhangir-azerbayev/MetaGen

max_forks_repo_name: zhangir-azerbayev/MetaGen


text:
"""
    cnse_boundary_state!(::Union{NumericalFluxFirstOrder, NumericalFluxGradient}, ::Insulating, ::HBModel)

apply insulating boundary condition for temperature
sets transmissive ghost point
"""
function cnse_boundary_state!(
    ::Union{NumericalFluxFirstOrder, NumericalFluxGradient},
    ::Insulating,
    ::CNSE3D,
    state⁺,
    aux⁺,
    n⁻,
    state⁻,
    aux⁻,
    t,
)
    state⁺.ρθ = state⁻.ρθ

    return nothing
end

"""
    cnse_boundary_state!(::NumericalFluxSecondOrder, ::Insulating, ::HBModel)

apply insulating boundary condition for velocity
sets ghost point to have no numerical flux on the boundary for κ∇θ
"""
@inline function cnse_boundary_state!(
    ::NumericalFluxSecondOrder,
    ::Insulating,
    ::CNSE3D,
    state⁺,
    gradflux⁺,
    aux⁺,
    n⁻,
    state⁻,
    gradflux⁻,
    aux⁻,
    t,
)
    state⁺.ρθ = state⁻.ρθ
    gradflux⁺.κ∇θ = n⁻ * -0

    return nothing
end

"""
    cnse_boundary_state!(::Union{NumericalFluxFirstOrder, NumericalFluxGradient}, ::TemperatureFlux, ::HBModel)

apply temperature flux boundary condition for velocity
applies insulating conditions for first-order and gradient fluxes
"""
function cnse_boundary_state!(
    nf::Union{NumericalFluxFirstOrder, NumericalFluxGradient},
    ::TemperatureFlux,
    model::CNSE3D,
    args...,
)
    return cnse_boundary_state!(nf, Insulating(), model, args...)
end

"""
    cnse_boundary_state!(::NumericalFluxSecondOrder, ::TemperatureFlux, ::HBModel)

apply insulating boundary condition for velocity
sets ghost point to have specified flux on the boundary for κ∇θ
"""
@inline function cnse_boundary_state!(
    ::NumericalFluxSecondOrder,
    bc::TemperatureFlux,
    ::CNSE3D,
    state⁺,
    gradflux⁺,
    aux⁺,
    n⁻,
    state⁻,
    gradflux⁻,
    aux⁻,
    t,
)
    state⁺.ρθ = state⁻.ρθ
    gradflux⁺.κ∇θ = n⁻ * bc(state⁻, aux⁻, t)

    return nothing
end


text:
export iauAtciq
"""
Quick ICRS, epoch J2000.0, to CIRS transformation, given precomputed
star-independent astrometry parameters.

Use of this function is appropriate when efficiency is important and
where many star positions are to be transformed for one date.  The
star-independent parameters can be obtained by calling one of the
functions iauApci[13], iauApcg[13], iauApco[13] or iauApcs[13].

If the parallax and proper motions are zero the iauAtciqz function
can be used instead.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   rc,dc  double     ICRS RA,Dec at J2000.0 (radians)
   pr     double     RA proper motion (radians/year; Note 3)
   pd     double     Dec proper motion (radians/year)
   px     double     parallax (arcsec)
   rv     double     radial velocity (km/s, +ve if receding)
   astrom iauASTROM* star-independent astrometry parameters:
      pmt    double       PM time interval (SSB, Julian years)
      eb     double[3]    SSB to observer (vector, au)
      eh     double[3]    Sun to observer (unit vector)
      em     double       distance from Sun to observer (au)
      v      double[3]    barycentric observer velocity (vector, c)
      bm1    double       sqrt(1-|v|^2): reciprocal of Lorenz factor
      bpn    double[3][3] bias-precession-nutation matrix
      along  double       longitude + s' (radians)
      xpl    double       polar motion xp wrt local meridian (radians)
      ypl    double       polar motion yp wrt local meridian (radians)
      sphi   double       sine of geodetic latitude
      cphi   double       cosine of geodetic latitude
      diurab double       magnitude of diurnal aberration vector
      eral   double       "local" Earth rotation angle (radians)
      refa   double       refraction constant A (radians)
      refb   double       refraction constant B (radians)

Returned:
   ri,di   double    CIRS RA,Dec (radians)

Notes:

   1. All the vectors are with respect to BCRS axes.

   2. Star data for an epoch other than J2000.0 (for example from the
   Hipparcos catalog, which has an epoch of J1991.25) will require a
   preliminary call to iauPmsafe before use.

   3. The proper motion in RA is dRA/dt rather than cos(Dec)*dRA/dt.

Called:
   iauPmpx      proper motion and parallax
   iauLdsun     light deflection by the Sun
   iauAb        stellar aberration
   iauRxp       product of r-matrix and pv-vector
   iauC2s       p-vector to spherical
   iauAnp       normalize angle into range 0 to 2pi

This revision:   2013 October 9

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

function iauAtciq(rc::Real, dc::Real,
                   pr::Real, pd::Real,
                   px::Real, rv::Real,
                   astrom::iauASTROM)

   # Allocate return value
   ref_astrom = Ref{iauASTROM}(astrom)
   ref_ri     = Ref{Float64}(0.0)
   ref_di     = Ref{Float64}(0.0)

   ccall((:iauAtciq, libsofa_c), Cvoid,
            (Cdouble, Cdouble,
            Cdouble, Cdouble,
            Cdouble, Cdouble, Ref{iauASTROM},
            Ref{Cdouble}, Ref{Cdouble}),
            convert(Float64, rc),
            convert(Float64, dc),
            convert(Float64, pr),
            convert(Float64, pd),
            convert(Float64, px),
            convert(Float64, rv),
            ref_astrom, ref_ri, ref_di)

   return ref_ri[], ref_di[]
end


text:
using JLD
using Distances
using ProgressMeter

if !isdefined(:loaded) || !loaded
    #jldopen("90.jld") do file
    jldopen("85.jld") do file
        global X, vocab, vocab2idx, idx2vocab
        X = read(file, "X")
        vocab = read(file, "vocab")
        vocab2idx = read(file, "vocab2idx")
        idx2vocab = read(file, "idx2vocab")
        global loaded = true
    end
    # @show length(vocab)
end

combined = readcsv("combined.csv")[2:end,:]

for idx in 1:size(combined,1)
    word1, word2, score = combined[idx,:]

    D = -1
    try
        v1 = X[vocab2idx[lowercase(word1)],:]
        v2 = X[vocab2idx[lowercase(word2)],:]
        D = dot(v1, v2) / norm(v1) / norm(v2)
        isnan(D) && (D = 0)
    catch e
        if isa(e, KeyError)
            D = 01
        else
            rethrow(e)
        end
    end

    println("$word1,$word2,$score,$D")
end


text:
using Documenter, DocumenterMarkdown
using Arrow

makedocs(
    format = Markdown(),
    modules = [Arrow],
    pages = [
        "Home" => "index.md",
        "User Manual" => "manual.md",
        "API Reference" => "reference.md"
    ]
)

deploydocs(repo = "github.com/JuliaData/Arrow.jl.git")

max_stars_repo_name: sweb/arrow

max_issues_repo_name: sweb/arrow

max_forks_repo_name: sweb/arrow


text:
using WAVI
function MISMIP_PLUS_Ice1r_2km()
    #Grid and boundary conditions
    nx = 320
    ny = 40
    nσ = 4
    x0 = 0.0
    y0 = -40000.0
    dx = 2000.0
    dy = 2000.0
    h_mask=trues(nx,ny)
    u_iszero = falses(nx+1,ny); u_iszero[1,:].=true
    v_iszero=falses(nx,ny+1); v_iszero[:,1].=true; v_iszero[:,end].=true
    grid = Grid(nx = nx,
                ny = ny,
                nσ = nσ,
                x0 = x0,
                y0 = y0,
                dx = dx,
                dy = dy,
                h_mask = h_mask,
                u_iszero = u_iszero,
                v_iszero = v_iszero)

    #Bed
    bed = WAVI.mismip_plus_bed #function definition of bed

    # Inputing thickness profile, reading from binary file
    fname = "examples\\Favier2019_melt_params\\data\\MISMIP_ice0_2km_SteadyThickness.bin";
    # fname = joinpath(dirname(@__FILE__), "data",  "MISMIP_ice0_2km_SteadyThickness.bin")
    h = Array{Float64,2}(undef, nx, ny)
    read!(fname, h)
    h = ntoh.(h)

    #solver parameters
    maxiter_picard = 30
    solver_params = SolverParams(maxiter_picard = maxiter_picard)

    # melt rate Model
    pme = PlumeEmulator(α=1.5)

    #Physical parameters
    accumulation_rate = 0.3
    params = Params(accumulation_rate = accumulation_rate)

    #make the model
    model = Model(grid = grid,
                     bed_elevation = bed,
                     params = params,
                     melt_rate = pme,
                     solver_params = solver_params,
                     initial_conditions = initial_conditions)

    #timestepping parameters
    niter0 = 0
    dt = 0.02
    end_time = 100.
    chkpt_freq = 1.
    pchkpt_freq = 1.
    timestepping_params = TimesteppingParams(niter0 = niter0,
                                            dt = dt,
                                            end_time = end_time,
                                            chkpt_freq = chkpt_freq,
                                            pchkpt_freq = pchkpt_freq)

    #output parameters
    folder = "outputs"
    isdir(folder) && rm(folder, force = true, recursive = true)
    mkdir(folder) #make a clean folder for outputs
    outputs = (h   = model.fields.gh.h,
                u  = model.fields.gh.u,
                v  = model.fields.gh.v,
                melt = model.fields.gh.basal_melt,
                grfrac = model.fields.gh.grounded_fraction)
    output_freq = 1.
    output_params = OutputParams(outputs = outputs,
                            output_freq = output_freq,
                            output_format = "mat",
                            output_path = folder,
                            zip_format = "nc")

    simulation = Simulation(model = model,
                        timestepping_params = timestepping_params,
                        output_params = output_params)

    #perform the simulation
    run_simulation!(simulation)
    return simulation
end

@time simulation = Favier2019_PME_retreat();


text:
using Unitful, UnitfulAngles
import Dates
using Test

units = (u"doubleTurn", u"turn", u"halfTurn", u"quadrant", u"sextant", u"octant", u"clockPosition", u"hourAngle", u"compassPoint", u"hexacontade", u"brad", u"°", u"grad", u"arcminute", u"arcsecond", u"rad", u"diameterPart")
quantities = (0.5, 1, 2, 4, 6, 8, 12, 24, 32, 60, 256, 360, 400, 21600, 1296000, 2π, 120π)

@test all(1u"turn" ≈ q*u for (q, u) in zip(quantities, units))
for _f in (:sin, :cos, :tan, :sec, :csc, :cot), (q, u) in zip(quantities, units), a in 13:17
    @test @eval $_f(2π/$a) ≈ $_f($q*$u/$a)
end

for (_f, _x) in zip((:asin, :acos, :atan, :asec, :acsc, :acot), (.5, √3/2, √3/3, 2/√3, 2, √3))
    @test @eval $_f(u"°", $_x) ≈ 30u"°"
end

@test atan(u"°", 1,1) == 45u"°"

@test convert(Dates.Time, 45u"°") == Dates.Time(3,0,0)

@test convert(u"°", Dates.Time(12,0,0)) == 180u"°"

####################
# for _f in (:sin, :cos), (q, u) in zip(qrats, rats)
    # @test @eval $(Symbol("$(_f)pi"))(2) === $_f($q*$u)
# end


text:
@testset "Repositories" begin
    @test SL.repository_ordering_is_valid()
    @test SL.prettify_loc(1000) == "1,000"

    test_repo(; kwargs...) = Repo("https://example.com", "test/test"; kwargs...)

    p = Pattern(1, "Title", ["Julia", 1.3], "Description", "toml")
    r = test_repo(; predicates=[])
    @test SL.filter_patterns(r, [p]) == [p]
    r = test_repo()
    @test SL.filter_patterns(r, [p]) == []
    r = test_repo(; predicates=[SL.version(1.5)])
    @test SL.filter_patterns(r, [p]) == [p]
end

max_stars_repo_name: rikhuijzer/display-lint

max_issues_repo_name: rikhuijzer/display-lint

max_forks_repo_name: rikhuijzer/ShowLint


text:
module Pseudea

# package code goes here
include("lattice.jl") # definition of lattice types

end # module


text:
using FactCheck

include("fred.jl")
include("yahoo.jl")
include("colors.jl")
include("const.jl")
include("plot.jl")
include("kplot.jl")
include("uplot.jl")
include("uhist.jl")
include("first.jl")
include("last.jl")

exitstatus()


text:
mutable struct SwitchRegenerateState{T}
    weight::Float64
    score::Float64
    noise::Float64
    prev_trace::Trace
    trace::Trace
    index::Int
    retdiff::Diff
    SwitchRegenerateState{T}(weight::Float64, score::Float64, noise::Float64, prev_trace::Trace) where T = new{T}(weight, score, noise, prev_trace)
end

function process!(gen_fn::Switch{C, N, K, T},
                  index::Int,
                  index_argdiff::Diff,
                  args::Tuple,
                  kernel_argdiffs::Tuple,
                  selection::Selection,
                  state::SwitchRegenerateState{T}) where {C, N, K, T}
    branch_fn = getfield(gen_fn.branches, index)
    merged = get_selected(get_choices(state.prev_trace), complement(selection))
    new_trace, weight = generate(branch_fn, args, merged)
    retdiff = UnknownChange()
    weight -= project(state.prev_trace, complement(selection))
    weight += (project(new_trace, selection) - project(state.prev_trace, selection))
    state.index = index
    state.weight = weight
    state.noise = project(new_trace, EmptySelection()) - project(state.prev_trace, EmptySelection())
    state.score = get_score(new_trace)
    state.trace = new_trace
    state.retdiff = retdiff
end

function process!(gen_fn::Switch{C, N, K, T},
                  index::Int,
                  index_argdiff::NoChange,
                  args::Tuple,
                  kernel_argdiffs::Tuple,
                  selection::Selection,
                  state::SwitchRegenerateState{T}) where {C, N, K, T}
    new_trace, weight, retdiff = regenerate(getfield(state.prev_trace, :branch), args, kernel_argdiffs, selection)
    state.index = index
    state.weight = weight
    state.noise = project(new_trace, EmptySelection()) - project(state.prev_trace, EmptySelection())
    state.score = get_score(new_trace)
    state.trace = new_trace
    state.retdiff = retdiff
end

@inline process!(gen_fn::Switch{C, N, K, T}, index::C, index_argdiff::Diff, args::Tuple, kernel_argdiffs::Tuple, selection::Selection, state::SwitchRegenerateState{T}) where {C, N, K, T} = process!(gen_fn, getindex(gen_fn.cases, index), index_argdiff, args, kernel_argdiffs, selection, state)

function regenerate(trace::SwitchTrace{T},
                    args::Tuple,
                    argdiffs::Tuple,
                    selection::Selection) where T
    gen_fn = trace.gen_fn
    index, index_argdiff = args[1], argdiffs[1]
    state = SwitchRegenerateState{T}(0.0, 0.0, 0.0, trace)
    process!(gen_fn, index, index_argdiff, args[2 : end], argdiffs[2 : end], selection, state)
    return SwitchTrace(gen_fn, state.index, state.trace, get_retval(state.trace), args, state.score, state.noise), state.weight, state.retdiff
end

max_stars_repo_name: georgematheos/Gen

max_issues_repo_name: georgematheos/Gen

max_forks_repo_name: georgematheos/Gen


text:
using StatsBase
using Random
using Memoize

"""
    weights(state::Vector{ComplexF64})

extended StatsBase weights to allow complex values
"""
function StatsBase.weights(state::Vector{ComplexF64})
    weights(norm.(state).^2)
end

"""
    measure(state:Vector{ComplexF64}, n::Int)

sample 'n' possible quantum state after measurement.
"""
function measure(state::Vector{ComplexF64}, shots::Int=1)

    (shots > 0) || error("Shots has to be a natural integer above 0.")
    ispow2(length(state)) || error("Statevector is not a viable quantum state (length $(length(state) )")
    norm(state) ≈ 1 || error("Statevector is not a viable statevector (norm is $(norm(state))")
    N = Int(log(2, length(state)))

    # generate samples and count number of measurements of each basis state
    s = sample(Random.GLOBAL_RNG, [0:(2^N)-1 ... ], weights(state), shots)
    u=unique(s)
    d=Dict([(i,count(x->x==i,s)) for i in u])

    return d
end


text:
# [
#  [ MPIStateArray Name, Field Name, Maximum, Minimum, Mean, Standard Deviation ],
#  [         :                :          :        :      :          :           ],
# ]

parr = [
    ["∫u", "∫x[1]", 12, 12, 0, 12],
    ["∫u", "∫x[2]", 15, 15, 15, 15],
    ["U", "η", 12, 12, 0, 12],
    ["U", "U[1]", 12, 12, 0, 12],
    ["U", "U[2]", 15, 15, 15, 15],
]

initial = [
    [
        "∫u",
        "∫x[1]",
        -6.36104748927323058183e+01,
        6.36104748927324763486e+01,
        3.59432306140661235747e-14,
        3.16776758596447187699e+01,
    ],
    [
        "∫u",
        "∫x[2]",
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
    ],
    [
        "U",
        "η",
        -1.00000000000000000000e+00,
        1.00000000000000000000e+00,
        3.73034936274052574533e-18,
        7.07673146340372483110e-01,
    ],
    [
        "U",
        "U[1]",
        -9.95282537582876658533e-01,
        9.95282537582876547511e-01,
        -8.82958196287834196438e-18,
        7.07673146340372150043e-01,
    ],
    [
        "U",
        "U[2]",
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
    ],
]

day = [
    [
        "∫u",
        "∫x[1]",
        -6.40437964821475844701e+01,
        6.40437964821476981570e+01,
        2.21189111471176159638e-14,
        2.60887556680598677872e+01,
    ],
    [
        "∫u",
        "∫x[2]",
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
    ],
    [
        "U",
        "η",
        -8.52761706341222169847e-01,
        8.52761706341222169847e-01,
        -1.86517468137026310378e-17,
        6.03476559805077306109e-01,
    ],
    [
        "U",
        "U[1]",
        -3.15397076561132330141e+01,
        3.15397076561132294614e+01,
        -2.82332564752296143530e-15,
        2.24255960582435314166e+01,
    ],
    [
        "U",
        "U[2]",
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
    ],
]

month = [
    [
        "∫u",
        "∫x[1]",
        -3.51564130159214158766e+01,
        3.51564130159213874549e+01,
        -1.16415321826934820032e-14,
        1.41614248818105643579e+01,
    ],
    [
        "∫u",
        "∫x[2]",
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
    ],
    [
        "U",
        "η",
        -5.30601964901526002016e-01,
        5.30601964901526002016e-01,
        -1.24344978758017535116e-17,
        3.75492761956246756672e-01,
    ],
    [
        "U",
        "U[1]",
        -3.51564130159208829696e+01,
        3.51564130159208758641e+01,
        -2.47942593560761348802e-15,
        2.49971726354604832920e+01,
    ],
    [
        "U",
        "U[2]",
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
    ],
]

year = [
    [
        "∫u",
        "∫x[1]",
        -4.25110151351051179791e-01,
        4.25110151351050846724e-01,
        -1.36424205265939223736e-16,
        1.74619680058462317662e-01,
    ],
    [
        "∫u",
        "∫x[2]",
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
    ],
    [
        "U",
        "η",
        -4.39687965160425187072e-02,
        4.39687965160425187072e-02,
        -6.10622663543836061796e-19,
        3.11155365713074068268e-02,
    ],
    [
        "U",
        "U[1]",
        -4.25110151351044573964e-01,
        4.25110151351044518453e-01,
        -1.99959679826973105365e-17,
        3.02264961945818255717e-01,
    ],
    [
        "U",
        "U[2]",
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
    ],
]

decade = [
    [
        "∫u",
        "∫x[1]",
        -1.88298126764002322427e-12,
        1.88298126764002483985e-12,
        5.62482816536058859320e-28,
        7.73459738532225749582e-13,
    ],
    [
        "∫u",
        "∫x[2]",
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
    ],
    [
        "U",
        "η",
        -3.38352596570510824820e-15,
        3.38352596570510824820e-15,
        6.31088724176809474572e-34,
        2.39443046587488032416e-15,
    ],
    [
        "U",
        "U[1]",
        -1.88298126763999575929e-12,
        1.88298126763999575929e-12,
        3.78999175038214983828e-29,
        1.33885125866565214937e-12,
    ],
    [
        "U",
        "U[2]",
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
    ],
]

century = [
    [
        "∫u",
        "∫x[1]",
        -3.97994692621354800600e-134,
        3.97994692621354409568e-134,
        -1.20124985710986044942e-149,
        1.63481642745144852479e-134,
    ],
    [
        "∫u",
        "∫x[2]",
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
    ],
    [
        "U",
        "η",
        -2.07116984864437321907e-136,
        2.07116984864437321907e-136,
        -5.62108290883927252001e-153,
        1.46571128339547670108e-136,
    ],
    [
        "U",
        "U[1]",
        -3.97994692621348592969e-134,
        3.97994692621348544091e-134,
        -4.15669815960076969677e-150,
        2.82985128060348631090e-134,
    ],
    [
        "U",
        "U[2]",
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
        0.00000000000000000000e+00,
    ],
]

refVals = (
    intitial = (initial, parr),
    day = (day, parr),
    month = (month, parr),
    year = (year, parr),
    decade = (decade, parr),
    century = (century, parr),
)


text:
function xmetrics(x)
  2x
end


text:
__precompile__(true)

module TreeViews

"""
    hastreeview(x::T)::Bool

Called by a frontend to decide whether a tree view for type `T` is available.
Defaults to `false`.
"""
hastreeview(x) = false

"""
    numberofnodes(x)

Number of direct descendents.
Defaults to `fieldcount(typeof(x))`.
"""
numberofnodes(x::T) where {T} = fieldcount(T)

"""
    treelabel(io::IO, x, mime = MIME"text/plain"())

Prints `x`'s tree header to `io`.
Like with `Base.show` there are also methods with `mime::AbstractString` and no `mime` argument at all (which
falls back to `MIME"text/plain"()`). Please only overload the `treelabel(io::IO, x, mime::MIME)` form.
"""
treelabel(io::IO, x::T, mime::MIME"text/plain") where {T} = show(io, mime, T)
treelabel(io::IO, x::T, mime::AbstractString)  where {T} = treelabel(io, x, MIME(mime))
treelabel(io::IO, x::T)  where {T} = treelabel(io, x, MIME"text/plain"())

"""
    nodelabel(io::IO, x::T, i::Integer, mime::MIME"text/plain" = MIME"text/plain"())

Prints the label of `x`'s `i`-th child to `io`.
Like with `Base.show` there are also methods with `mime::AbstractString` and no `mime` argument at all (which
falls back to `MIME"text/plain"()`). Please only overload the `nodelabel(io::IO, x, i::Integer, mime::MIME)` form.
"""
function nodelabel(io::IO, x::T, i::Integer, mime::MIME"text/plain") where {T}
  show(io, mime, Text(String(fieldname(T, i))))
end
nodelabel(io::IO, x::T, i::Integer, mime::AbstractString)  where {T} = nodelabel(io, x, i, MIME(mime))
nodelabel(io::IO, x::T, i::Integer)  where {T} = nodelabel(io, x, i, MIME"text/plain"())

"""
    treenode(x::T, i::Integer)

Returns the `i`-th node of `x`, which is usually printed by the display frontend next to
the corresponding `treelabel`.
"""
treenode(x::T, i::Integer) where {T} = getfield(x, fieldname(T, i))

@deprecate(treelabel(io::IO, x::T, i::Integer, mime) where {T},
           nodelabel(io, x, i, mime))
@deprecate(treelabel(io::IO, x::T, i::Integer) where {T}, nodelabel(io, x, i))

end # module


text:
module AlgebraOfGraphics

using Base: front, tail
using Dates
using Tables: rows, columns, getcolumn, columnnames
using StructArrays: StructArrays, components, uniquesorted, GroupPerm, StructArray
using GeometryBasics: AbstractGeometry, Polygon, MultiPolygon
using GeoInterface: coordinates, geotype
using Colors: RGB, RGBA, red, green, blue, Color
using PlotUtils: optimize_datetime_ticks
using Makie
using Makie: current_default_theme, to_value, automatic, Automatic, PlotFunc, ATTRIBUTES
import Makie.MakieLayout: hidexdecorations!,
                          hideydecorations!,
                          hidedecorations!,
                          linkaxes!,
                          linkxaxes!,
                          linkyaxes!
using GridLayoutBase: determinedirsize, Col, Row
using PooledArrays: PooledArray
using KernelDensity: kde, pdf
using StatsBase: fit, histrange, Histogram, normalize, weights, AbstractWeights, sturges

import GLM, Loess
import FileIO
import RelocatableFolders

export hideinnerdecorations!, deleteemptyaxes!
export Entry, AxisEntries
export renamer, sorter, nonnumeric, verbatim
export density, histogram, linear, smooth, expectation, frequency
export visual, data, geodata, dims, mapping
export datetimeticks
export draw, draw!
export facet!, colorbar!, legend!
export set_aog_theme!

include("theme.jl")
include("helpers.jl")
include("scales.jl")
include("entries.jl")
include("facet.jl")
include("geometry.jl")
include("algebra/layer.jl")
include("algebra/layers.jl")
include("algebra/select.jl")
include("algebra/processing.jl")
include("recipes/linesfill.jl")
include("transformations/visual.jl")
include("transformations/linear.jl")
include("transformations/smooth.jl")
include("transformations/density.jl")
include("transformations/histogram.jl")
include("transformations/groupreduce.jl")
include("transformations/frequency.jl")
include("transformations/expectation.jl")
include("guides/legendelements.jl")
include("guides/legend.jl")

end


text:
#################### File I/O ####################
"""
    Base.read(name::AbstractString, ::Type{T}) where {T<:AbstractChains}

Read a chain from an external file.

Returns an `AbstractChains` subtype read from an external file.

* `name` : file to read or write. Recommended convention is for the file name to be specified with a `.jls` extension.

* `T` : chain type to read.

* `c` : chain to write.
"""
function Base.read(name::AbstractString, ::Type{T}) where {T<:AbstractChains}
    c = open(deserialize, name, "r")
    isa(c, T) || throw(TypeError(:open, "read(\"$name\", $T)", T, c))
    c
end
"""
    Base.write(name::AbstractString, c::AbstractChains)

Write a chain to an external file.

Returns a written external file containing a subtype.

* `name` : file to read or write. Recommended convention is for the file name to be specified with a `.jls` extension.

* `T` : chain type to read.

* `c` : chain to write.
"""
function Base.write(name::AbstractString, c::AbstractChains)
    open(file -> serialize(file, c), name, "w")
end

"""
    readcoda(output::AbstractString, index::AbstractString)

Read MCMC sampler output generated in the CODA format by OpenBUGS. The function only retains those sampler iterations at which all model parameters were monitored.

Returns a `Chains` object containing the read sampler output.

* `output` : text file containing the iteration numbers and sampled values for the model parameters.

* `index` : text file containing the names of the parameters, followed by the first and last rows in which their output can be found in the `output` file.

"""
function readcoda(output::AbstractString, index::AbstractString)
    out = readdlm(output, Any)
    ind = readdlm(index, Any)

    firstind = ind[:, 2]
    firstiter = out[firstind, 1]
    lastind = ind[:, 3]
    lastiter = out[lastind, 1]

    thin = Int((lastiter[1] - firstiter[1]) / (lastind[1] - firstind[1]))
    window = maximum(firstiter):thin:minimum(lastiter)
    startind = firstind .+ (first(window) .- firstiter) ./ step(window)
    stopind = lastind .- (lastiter .- last(window)) ./ step(window)

    names = AbstractString[ind[:, 1]...]

    value = Array{Float64}(undef, length(window), length(names))
    for i = 1:size(value, 2)
        inds = Int(startind[i]):Int(stopind[i])
        value[:, i] = out[inds, 2]
    end

    Chains(value, start = first(window), thin = step(window), names = names)
end

"""
    to_file(model::ModelChains, outpath::AbstractString)
This function writes the results of the MCMC runs into files. The destination of
the files is specified using `outpath`. It will create a files for each chain. A
`params_x.log` file storing each parameter sample. In this case `x` specifies the
index of the chain. The file is compatible with MCMC analysis tools like `Tracer`
(http://tree.bio.ed.ac.uk/software/tracer/). If in addition trees are sampled,
they are stored in newick format in a file called `trees_x.nwk`, where `x` again
specifies the index of the respective chain.
"""
function to_file(model::ModelChains, outpath::AbstractString)
    for run in 1:size(model.value,3)
        thin = model.range.step
        if isassigned(model.trees, 1)
            to_file(model.value[:, :, run], model.trees[:,:,run], outpath, string(run), thin, model.names, model.tree_names)
        else
            to_file(model.value[:, :, run], outpath, string(run), thin, model.names)
        end
    end

end

function to_file(df::Array, outpath::String, run::String, thin::Int64, names::Array)
    io = open(string(outpath, "params_"*run*".log"),"w")
    write(io, "it\t")
    join(io, names, "\t")
    write(io, "\n")
    for (ind,row) in enumerate(eachrow(df))
        write(io, "$(ind*thin)\t")
        join(io, row, "\t")
        write(io, "\n")
    end
    close(io)
end

function to_file(df::Array, tdf::Array, outpath::AbstractString,
                run::AbstractString, thin::Int64, names::Array, tree_names::Array{A}) where A <: AbstractString

    to_file(df, outpath, run, thin, names)
    for (ind, n_name) in enumerate(tree_names)
        io = open(string(outpath, "trees_"*n_name*"_"*run*".nwk"), "w")
        for x = 1:length(tdf[:,ind])
            write(io, tdf[x,ind])
            write(io, "\n")
        end
        close(io)
    end
end

max_stars_repo_name: erathorn/MCPhylo

max_issues_repo_name: erathorn/MCPhylo

max_forks_repo_name: erathorn/MCPhylo


text:
export AdaptiveForwardSensFspSparse, ForwardSensFspInitialConditionSparse
export forwardsens_initial_condition, get_states, get_probability, get_sensitivity

"""
    $(TYPEDEF)

Type to store algorithmic options for solving a Foward Sensitivity CME problem.

# Fields

$(FIELDS)

# See also
[`ForwardSensRStepAdapter`](@ref)
"""
Base.@kwdef struct AdaptiveForwardSensFspSparse <: ForwardSensCmeAlgorithm
    """
    Method to adapt state space when the current solution error exceeds tolerance.
    """
    space_adapter::AbstractForwardSensSpaceAdapterSparse
    """
    An ODE solver from `DifferentialEquations.jl` to integrate the FSP-truncated problem.
    """
    ode_method::Union{Nothing,AbstractODEAlgorithm}
end

"""
    $(TYPEDEF)

Initial condition for Foward Sensitivity FSP in sparse format.

# Fields

$(FIELDS)

"""
struct ForwardSensFspInitialConditionSparse{NS,IntT<:Integer,RealT<:AbstractFloat}
    """
    List of states with nonzero initial probabilities.
    """
    states::Vector{MVector{NS,IntT}}
    """
    List of nonzero probabilities, with `p[i]` being the probability of state `states[i]`.
    """
    p::Vector{RealT}
    """
    List of partial derivatives, with `S[j][i]` being the `j`-th partial derivative of `p[i]`.
    """
    S::Vector{Vector{RealT}}
end
get_states(ic::ForwardSensFspInitialConditionSparse) = ic.states
get_probability(ic::ForwardSensFspInitialConditionSparse) = ic.p
get_sensitivity(ic::ForwardSensFspInitialConditionSparse) = ic.S

"""
    $(TYPEDSIGNATURES)

Return a `ForwardSensFspInitialConditionSparse` instance from a list of initial states `states`, initial probabilities `probabilities` and initial sensitivities `sensitivity`.
"""
function forwardsens_initial_condition(states::Vector{<:AbstractVector{IntT}}, probabilities::Vector{RealT}, sensitivity::Vector{Vector{RealT}}) where {IntT<:Integer, RealT<:AbstractFloat}
    isempty(states) && throw(ArgumentError("Empty state list in input."))
    species_count = length(states[1])

    _states = Vector{MVector{species_count, IntT}}()
    for eachstate in states
        push!(_states, eachstate)
    end
    _p = copy(probabilities)
    _S = deepcopy(sensitivity)
    return ForwardSensFspInitialConditionSparse{species_count, IntT, RealT}(
        _states,
        _p,
        _S
    )
end

"""
$(TYPEDSIGNATURES)

Foward sensitivity Finite State Projection to compute the time-dependent probability distributions of a stochastic reaction network along with their partial derivatives with respect to model parameters.

# Arguments
- `model`: Stochastic reaction network model with propensity sensitivities.
- `initial_condition`: Initial condition.
- `tspan`: Time span.
- `sensfspalgorithm`: Forward Sesitivity FSP method, a combination of space adaptation strategy and the ODE method to integrate the FSP-truncated problem.
- `saveat`: timepoints to save the intermediate solutions. If left empty, all timesteps computed by the ODE solver is saved.
- `fsptol`: FSP error tolerance. The solver dynamically adapts the projection space so that the sum of probability mass in the truncated solution is above `1.0 - fsptol`. Default: 1.0E-6.
- `odeatol`: Absolute tolerance for the ODE solver. Default: 1.0E-10.
- `odertol`: Relative tolerance for the ODE solver. Default: 1.0E-4.
- `verbose`: Whether to print status when adapting state space. Default: false.

# Returns
A `ForwardSensFspOutputSparse` instance.

# See also
[`FspVectorSparse`](@ref), [`ForwardSensFspOutputSparse`](@ref), [`AdaptiveForwardSensFspSparse`](@ref).
"""
function solve(model::CmeModelWithSensitivity,
    initial_condition::ForwardSensFspInitialConditionSparse{NS, IntT, RealT},
    tspan::Tuple{AbstractFloat,AbstractFloat},
    sensfspalgorithm::AdaptiveForwardSensFspSparse;
    saveat = [],
    fsptol::AbstractFloat = 1.0E-6,
    odeatol::AbstractFloat = 1.0E-10,
    odertol::AbstractFloat = 1.0E-4,
    verbose::Bool=false) where {NS, IntT<:Integer, RealT<:AbstractFloat}

    tstart = min(tspan...)
    tend = max(tspan...)
    parameter_count = get_parameter_count(model)

    (length(get_sensitivity(initial_condition)) ≠ parameter_count) && throw(ArgumentError("Initial condition does not match CME model. Initial condition must contain `np` sensitivity vectors where `np` is the number of CME model parameters."))

    adapter = sensfspalgorithm.space_adapter
    statespace = StateSpaceSparse(get_stoich_matrix(model), get_states(initial_condition))
    sink_count = get_sink_count(statespace)
    p0 = copy(get_probability(initial_condition))
    S0 = deepcopy(get_sensitivity(initial_condition))
    init!(statespace, adapter, p0, S0, tstart, fsptol)

    output = ForwardSensFspOutputSparse{NS,IntT,RealT}(
        t = Vector{RealT}(),
        p = Vector{FspVectorSparse{NS,IntT,RealT}}(),
        sinks = Vector{Vector{RealT}}(),
        S = Vector{Vector{FspVectorSparse{NS,IntT,RealT}}}(),
        dsinks = Vector{Vector{Vector{RealT}}}()
    )

    tnow = tstart
    n = length(p0)+sink_count
    unow = zeros(RealT, n*(parameter_count+1))
    copyto!(unow, 1, p0, 1, length(p0))
    for i in 1:parameter_count
        copyto!(unow, i*n+1, S0[i], 1, length(S0[i]))
    end
    while tnow < tend
        SA = ForwardSensFspMatrixSparse{RealT}(model, statespace)
        n = get_rowcount(SA.fspmatrix)
        function sensfsprhs!(du, u, θ, t)
            matvec!(du, t, SA, u)
        end
        function affect!(integrator)
            DE.terminate!(integrator)
        end
        # Error constraint for the intermediate solutions. The intermediate sinks must not grow beyond a linear function of time that reaches `fsptol` at the end time
        function fsp_error_constraint(u, t, integrator)
            sinks = u[n-sink_count+1:n]
            return sum(sinks) - fsptol * t / tend
        end

        fsp_cb = DE.ContinuousCallback(
            fsp_error_constraint,
            affect!,
            save_positions = (false, false),
            interp_points = 100,
            abstol = eps()
        )

        fsensfspprob = DE.ODEProblem(sensfsprhs!, unow, (tnow, tend), p = get_parameters(model))
        integrator = DE.init(fsensfspprob, sensfspalgorithm.ode_method, atol = odeatol, rtol = odertol, callback = fsp_cb, saveat = saveat, progress=true)
        DE.step!(integrator, tend - tnow, true)

        for (t, u) in zip(integrator.sol.t, integrator.sol.u)
            push!(output.t, t)
            push!(output.p, FspVectorSparse(statespace.states, u[1:n-sink_count]))
            push!(output.sinks, u[n-sink_count+1:n])
            push!(output.S, [FspVectorSparse(get_states(statespace), u[ip*n+1:(ip+1)*n-sink_count]) for ip in 1:parameter_count])
            push!(output.dsinks, [u[(ip+1)*n-sink_count+1:(ip+1)*n] for ip in 1:parameter_count])
        end

        tnow = integrator.t
        # If the ODE integrator exits before reaching `tend`, it means the intermediate solution with current state space has reached error tolerance threshold
        if tnow < tend
            n = get_rowcount(SA.fspmatrix)
            p = integrator.u[1:n-sink_count]
            sinks = integrator.u[n-sink_count+1:n]
            S = [integrator.u[ip*n+1:(ip+1)*n-sink_count] for ip in 1:parameter_count]
            dsinks = [integrator.u[(ip+1)*n-sink_count+1:(ip+1)*n] for ip in 1:parameter_count]

            # Call adapter to modify the state space and shrink/expand the solution vector as appropriate
            adapt!(statespace, adapter, p, S, sinks, dsinks, tnow, tend, fsptol, integrator=integrator)

            # This numerical trick ensures that DifferentialEquations.jl's event detection does not miss the event the updated state space is insufficient
            if sum(sinks) >= tnow * fsptol / tend
                sinks .-= eps()
            end

            # Extended solution vector to start the ODE integrator again
            n = length(p)+sink_count
            unow = zeros(RealT, n*(parameter_count+1))
            copyto!(unow, 1, p, 1, length(p))
            copyto!(unow, n - sink_count+1, sinks, 1, sink_count)
            for ip in 1:parameter_count
                copyto!(unow, ip*n+1, S[ip], 1, length(S[ip]))
                copyto!(unow, (ip+1)*n - sink_count + 1, dsinks[ip], 1, sink_count)
            end
            if verbose
                println("At t = $(round(tnow, digits=2)): update sate space. New size: $(get_state_count(statespace)) states.")
            end
        else # Otherwise, add the final slice to the output
            u = integrator.u
            push!(output.t, tnow)
            push!(output.p, FspVectorSparse(statespace.states, u[1:n-sink_count]))
            push!(output.sinks, u[n-sink_count+1:n])
            push!(output.S, [FspVectorSparse(get_states(statespace), u[ip*n+1:ip*n+n-sink_count]) for ip in 1:parameter_count])
            push!(output.dsinks, [u[ip*n+n-sink_count+1:ip*n+n] for ip in 1:parameter_count])
        end
    end
    return output
end


text:
function listencodeblocks(encodings, blocks)
    filledblocks = Any[blocks]
    changedblocks = Any[tuplemap(_ -> false, blocks)]
    for encoding in encodings
        push!(changedblocks, tuplemap(x -> !isnothing(x), encodedblock(encoding, blocks)))
        blocks = encodedblockfilled(encoding, blocks)
        push!(filledblocks, blocks)
    end
    return filledblocks, changedblocks
end

function listdecodeblocks(encodings, blocks)
    filledblocks = Any[blocks]
    changedblocks = Any[tuplemap(_ -> false, blocks)]
    for encoding in Iterators.reverse(encodings)
        push!(changedblocks, tuplemap(x -> !isnothing(x), decodedblock(encoding, blocks)))
        blocks = decodedblockfilled(encoding, blocks)
        push!(filledblocks, blocks)
    end
    return filledblocks, changedblocks
end

tuplemap(f, args...) = f(args...)
tuplemap(f, args::Vararg{Tuple}) = map((as...) -> tuplemap(f, as...), args...)

function blockcolumn(encodings, block; decode = false)
    blocks, changed = decode ? listdecodeblocks(encodings, block) : listencodeblocks(encodings, block)
    blockscol = [tuplemap((b, c) -> c ? "**`$(typeof(b))`**" : "`$(typeof(b))`", bs, ch) for (bs, ch) in zip(blocks, changed)]
    if block isa Tuple
        blockscol = [join(row, ", ") for row in blockscol]
    end
    return reshape(blockscol, :, 1)
end

encodingscolumn(encodings) = reshape(
    ["", ["`$(typeof(enc).name.name)`" for enc in encodings]...], :, 1)

function describeencodings(
        encodings,
        blocks::Tuple;
        inname="Input",
        outname="Output",
        blocknames=repeat([""], length(blocks)),
        decode=false,
        markdown=false,
        tf=tf_markdown)
    namescol = reshape([inname, ["" for _ in 2:length(encodings)]..., outname], :, 1)

    data = hcat(
        encodingscolumn(decode ? reverse(encodings) : encodings),
        namescol,
        [blockcolumn(encodings, block; decode=decode) for block in blocks]...,
    )

    s = pretty_table(
        String,
        data,
        header=[decode ? "Decoding" : "Encoding", "Name", blocknames...],
        alignment=[:r, :r, [:l for _ in 1:length(blocknames)]...],
        tf=tf,
        )
    return markdown ? Markdown.parse(s) : s
end

function describemethod(method::BlockMethod)
    xblock = encodedblock(method.encodings, method.blocks[1])
    s = """
    #### `LearningMethod` summary

    - Task: `$(typeof(method.blocks[1])) -> $(typeof(method.blocks[2]))`
    - Model blocks: `$(typeof(xblock)) -> $(typeof(method.outputblock))`

    Encoding a sample (`encode(method, context, sample)`)

    $(describeencodings(
        method.encodings,
        method.blocks,
        blocknames=["`method.blocks[1]`", "`method.blocks[2]`"],
        inname="`(input, target)`", outname="`(x, y)`"))

    Decoding a model output (`decode(method, context, ŷ)`)

    $(describeencodings(
        method.encodings,
        (method.outputblock,),
        blocknames=["`method.outputblock`"],
        inname="`ŷ`", outname="`target_pred`", decode=true))
    """

    return Markdown.parse(s)
end


text:
struct FNN <: AbstractApproximator
    n::Int
    m::Int
    NN::Flux.Chain
    function FNN(n::Int, m::Int, h_array::Vector{Int}, act)
        node_array = [n+m, h_array..., 1]
        new(n, m, construct_layer_array(node_array, act))
    end
end

"""
    (nn::FNN)(x, u)

Infer using two input arguments `x` and `u`.
# Notes
size(x) = (n, d)
size(u) = (m, d)
"""
function (nn::FNN)(x, u)
    res = nn.NN(vcat(x, u))
end

Flux.params(approximator::FNN) = Flux.params(approximator.NN)


text:
function sledge_down(trees)
    is_tree(i, j) = trees[i, j] == '#'

    x, y = 1, 1
    cnt = 0
    nr_rows, nr_cols = size(trees)
    while x <= nr_rows
        if is_tree(x, y)
            cnt += 1
        end
        x += 1
        y = ((y + 2) % nr_cols) + 1
    end
    cnt
end

function main()
    lines = readlines("../resources/input_03.txt")
    trees = permutedims(hcat(map(collect, lines)...))
    result = sledge_down(trees)
    println(result)
end

main()


text:
using Documenter
using MixedModelsSim

makedocs(
    modules = [MixedModelsSim],
    root = joinpath(dirname(pathof(MixedModelsSim)), "..", "docs"),
    sitename = "MixedModelsSim.jl",
    doctest = true,
    authors = "Phillip Alday, Douglas Bates, Lisa DeBruine, Reinhold Kliegl",
    pages = ["Home" => "index.md"
             "Full Example" => "simulation.md"],
)

deploydocs(; repo = "github.com/RePsychLing/MixedModelsSim.jl", push_preview = true)


text:
using YaoPlots
using Compose
using Test
using Yao

@testset "gate styles" begin
    c = YaoPlots.CircuitGrid(1)
	@test YaoPlots.get_brush_texts(c, X)[2] == "X"
	@test YaoPlots.get_brush_texts(c, Rx(0.5))[2] == "Rx(0.5)"
	@test YaoPlots.get_brush_texts(c, shift(0.5))[2] == "ϕ(0.5)"
	@test YaoPlots.get_brush_texts(c, YaoBlocks.phase(0.5))[2] == "    0.5\n"
end

@testset "circuit canvas" begin
	c = CircuitGrid(5)
	@test YaoPlots.nline(c) == 5
	@test YaoPlots.frontier(c, 2, 3) == 0
	@test YaoPlots.depth(c) == 0
	circuit_canvas(5) do c
		YaoPlots.draw!(c, put(5, 3=>X), 1:5, [])
		@test YaoPlots.frontier(c, 1, 2) == 0
		@test YaoPlots.frontier(c, 3, 5) == 1
		@test YaoPlots.depth(c) == 1
	end

	gg = circuit_canvas(5) do c
		put(3=>X) >> c
		put(3=>matblock(rand_unitary(2); tag="xxx")) >> c
		put(3=>igate(1)) >> c
		time_evolve(put(3,3=>igate(1)), 0.1) >> c
		control(2, 3=>X) >> c
		control(2, 3=>shift(0.5)) >> c
		chain(5, control(2, 3=>X), put(1=>X)) >> c
		@test YaoPlots.depth(c) >= 6
	end
	@test gg isa Context

	g = put(5, (3, 4)=>SWAP) |> vizcircuit(; scale=0.7, w_line=0.8, w_depth=0.9)
	@test g isa Context
	@test vizcircuit(put(5, (3,4)=>kron(X, Y)); scale=0.7, w_line=0.8, w_depth=0.9) isa Context

	@test vizcircuit(control(10, (2, -3), 6=>X)) isa Context
	@test vizcircuit(control(10, (2, -3), 6=>im*X)) isa Context
	@test plot(put(7, (2,3)=>matblock(randn(4,4)))) isa Context
end

@testset "fix #3" begin
	@test (control(4, 2, (1,3)=>kron(X, X)) |> vizcircuit) isa Context
end

@testset "pretty_angle" begin
	@test YaoPlots.pretty_angle(π) == "π"
	@test YaoPlots.pretty_angle(π*1.0) == "π"
	@test YaoPlots.pretty_angle(-π*1.0) == "-π"
	@test YaoPlots.pretty_angle(-π*0.0) == "0"
	@test YaoPlots.pretty_angle(-π*0.5) == "-π/2"
	@test YaoPlots.pretty_angle(1.411110) == "1.41"
end

@testset "readme" begin
    YaoPlots.CircuitStyles.linecolor[] = "pink"
    YaoPlots.CircuitStyles.gate_bgcolor[] = "yellow"
    YaoPlots.CircuitStyles.textcolor[] = "#000080" # the navy blue color
    YaoPlots.CircuitStyles.fontfamily[] = "JuliaMono"
    YaoPlots.CircuitStyles.lw[] = 2.5pt
    YaoPlots.CircuitStyles.textsize[] = 13pt
    YaoPlots.CircuitStyles.paramtextsize[] = 8pt

    @test plot(chain(3, put(1=>X), repeat(3, H), put(2=>Y), repeat(3, Rx(π/2)))) isa Compose.Context
    darktheme!()
    @test plot(chain(3, put(1=>X), repeat(3, H), put(2=>Y), repeat(3, Rx(π/2)))) isa Compose.Context
    lighttheme!()
    @test plot(chain(3, put(1=>X), repeat(3, H), put(2=>Y), repeat(3, Rx(π/2)))) isa Compose.Context
end


text:
using Test, InteractiveChaos, Observables

include("od_test.jl")
include("psos_test.jl")