text:
@par function write_spec(s::A) where {A<:@par(Simulation)}
    i = s.iteration[]
    hout = s.hspec
    vout = s.vspec

    outH = s.specH
    outV = s.specV
    tspecH = s.tspecH

    spectrum_u(hout,vout,s.u)
    writespectrum("u",i,hout,vout,outH,outV,tspecH)

    spectrum_viscosity(hout,vout,hout,vout,s)
    writespectrum("vis",i,hout,vout,outH,outV,tspecH)

    if (SPEC_FIL_D2 != 0.0)
      spectrum_u(hout,vout,s.u)
      filter_spectrum!(hout,vout)
      writespectrum("u_fil",i,hout,vout,outH,outV,tspecH)

      spectrum_viscosity(hout,vout,hout,vout,s)
      writespectrum("vis_fil",i,hout,vout,outH,outV,tspecH)
    end

    spectrum_non_linear(hout,vout,s.u,s.rhs,s.kspec)
    writespectrum("nl",i,hout,vout,outH,outV,tspecH)

    if (SPEC_FIL_D2 != 0.0)
      filter_spectrum!(hout,vout)
      writespectrum("nl_fil",i,hout,vout,outH,outV,tspecH)
    end

    spectrum_pressure(hout,vout,s)
    writespectrum("press",i,hout,vout,outH,outV,tspecH)

    if (SPEC_FIL_D2 != 0.0)
      filter_spectrum!(hout,vout)
      writespectrum("press_fil",i,hout,vout,outH,outV,tspecH)
    end

    if hasdensity(A)
        spectrum_buoyancy(vout,s.u,s.densitystratification.ρ,gravity(s))
        #write("buoyancy_v.spec3D.$i",vout)

        calculate_specHV(outH,outV,vout)

        writedlm("buoyancy_v.specH.$i",zip(KH,outH))
        writedlm("buoyancy_v.specV.$i",zip(KRZ,outV))

        outH .= (-).(outH ./ gravity(s.densitystratification).z) .* meangradient(s.densitystratification)
        outV .= (-).(outV ./ gravity(s.densitystratification).z) .* meangradient(s.densitystratification)

        writedlm("rhosource.specH.$i",zip(KH,outH))
        writedlm("rhosource.specV.$i",zip(KRZ,outV))

        spectrum_rho(hout,s.densitystratification.ρ)
        #write("rho.spec3D.$i",hout)

        calculate_specHV(outH,outV,hout)

        writedlm("rho.specH.$i",zip(KH,outH))
        writedlm("rho.specV.$i",zip(KRZ,outV))

        spectrum_rhodiss(hout,hout,s)

        #write("rhodiss.spec3D.$i",hout)

        calculate_specHV(outH,outV,hout)

        writedlm("rhodiss.specH.$i",zip(KH,outH))
        writedlm("rhodiss.specV.$i",zip(KRZ,outV))

        spectrum_rhonl(hout,s.densitystratification.ρ,s.densitystratification.flux)
        #write("rhonl.spec3D.$i",hout)

        calculate_specHV(outH,outV,hout)

        writedlm("rhonl.specH.$i",zip(KH,outH))
        writedlm("rhonl.specV.$i",zip(KRZ,outV))
        if (SPEC_FIL_D2 != 0.0)
            spectrum_buoyancy(vout,s.u,s.densitystratification.ρ,gravity(s))
            filter_spectrum!(hout,vout)
            #write("buoyancy_v.spec3D.$i",vout)

            calculate_specHV(outH,outV,vout)

            writedlm("buoyancy_fil_v.specH.$i",zip(KH,outH))
            writedlm("buoyancy_fil_v.specV.$i",zip(KRZ,outV))

            outH .= (-).(outH ./ gravity(s.densitystratification).z) .* meangradient(s.densitystratification)
            outV .= (-).(outV ./ gravity(s.densitystratification).z) .* meangradient(s.densitystratification)

            writedlm("rhosource_fil.specH.$i",zip(KH,outH))
            writedlm("rhosource_fil.specV.$i",zip(KRZ,outV))

            spectrum_rho(hout,s.densitystratification.ρ)
            filter_spectrum!(hout,vout)
            #write("rho.spec3D.$i",hout)

            calculate_specHV(outH,outV,hout)

            writedlm("rho_fil.specH.$i",zip(KH,outH))
            writedlm("rho_fil.specV.$i",zip(KRZ,outV))

            spectrum_rhodiss(hout,hout,s)
            filter_spectrum!(hout,vout)

            #write("rhodiss.spec3D.$i",hout)

            calculate_specHV(outH,outV,hout)

            writedlm("rhodiss_fil.specH.$i",zip(KH,outH))
            writedlm("rhodiss_fil.specV.$i",zip(KRZ,outV))

            spectrum_rhonl(hout,s.densitystratification.ρ,s.densitystratification.flux)
            filter_spectrum!(hout,vout)
            #write("rhonl.spec3D.$i",hout)

            calculate_specHV(outH,outV,hout)

            writedlm("rhonl_fil.specH.$i",zip(KH,outH))
            writedlm("rhonl_fil.specV.$i",zip(KRZ,outV))
        end
    end

    if hasles(A)
        spectrum_les(hout,vout,s.u,s.lesmodel.tau)
        writespectrum("les",i,hout,vout,outH,outV,tspecH)
    end

    if hashyperviscosity(A)
        if HyperViscosityType <: SpectralBarrier
            spectrum_spectral_barrier(hout,vout,s.u,s.hyperviscosity)
        else
            spectrum_hyperviscosity(hout,vout,s.u,s)
        end
        writespectrum("hvis",i,hout,vout,outH,outV,tspecH)

        if (SPEC_FIL_D2 != 0.0)
          filter_spectrum!(hout,vout)
          writespectrum("hvis_fil",i,hout,vout,outH,outV,tspecH)
        end
    end

end

@par function write_spec_forcing(s::A) where {A<:@par(Simulation)}
    i = s.iteration[]
    hout = s.hspec

    outH = s.specH
    outV = s.specV
    tspecH = s.tspecH

    fx = s.forcing.forcex
    fy = s.forcing.forcey
    ux = s.u.c.x
    uy = s.u.c.y
    dt = get_dt(s)
    spectrum_forcing(hout,ux,uy,fx,fy,dt)
    #write("forcing_h.spec3D.$i",hout)

    calculate_specHV(outH,outV,hout)

    writedlm("forcing_h.specH.$i",zip(KH,outH))
    writedlm("forcing_h.specV.$i",zip(KRZ,outV))

    return nothing
end

function writespectrum(n::String,i::Integer,hout,vout,specH,specV,tspecH)
    #write("$(n)_h.spec3D.$i",hout)
    #write("$(n)_v.spec3D.$i",vout)

    calculate_specHV(specH,specV,hout)

    copy!(tspecH,specH)

    writedlm("$(n)_h.specH.$i",zip(KH,specH))
    writedlm("$(n)_h.specV.$i",zip(KRZ,specV))

    calculate_specHV(specH,specV,vout)

    tspecH .+= specH
    writedlm("$(n).specH.$i",zip(KH,tspecH))

    writedlm("$(n)_v.specH.$i",zip(KH,specH))
    writedlm("$(n)_v.specV.$i",zip(KRZ,specV))

    return nothing
end

function writedlm2D(n,out)
    open(n,"w+") do f
        for k in RZRANGE
            for j in eachindex(KH)
                write(f,join((KZ[k],KH[j],out[j,k]),'\t'),'\n')
            end
        end
    end
end

vecouterproj(a,b) = Vec(proj(a.x,b.x),
                        proj(a.y,b.y),
                        proj(a.z,b.z))

function spectrum_non_linear(hout,vout,u,nl,ke)
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        k = ZRANGE[kk]
        j = YRANGE[jj]
            @simd for i in XRANGE
                ∇ = im*K[i,j,k]
                out = vecouterproj(nl[i,j,k] - ∇*ke[i,j,k],u[i,j,k])
                hout[i,j,k] = out.x + out.y
                vout[i,j,k] = out.z
            end
    end
end

function spectrum_les(hout,vout,u,τ)
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        k = ZRANGE[kk]
        j = YRANGE[jj]
            @inbounds @simd for i in XRANGE
                kh = K[i,j,k]
                out = vecouterproj((im*kh)⋅τ[i,j,k],u[i,j,k])
                hout[i,j,k] = out.x + out.y
                vout[i,j,k] = out.z
            end
    end
end

function spectrum_rhonl(hout,rho,flux)
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        k = ZRANGE[kk]
        j = YRANGE[jj]
            @inbounds @simd for i in XRANGE
                kh = K[i,j,k]
                out = proj((im*kh)⋅flux[i,j,k],rho[i,j,k])
                hout[i,j,k] = out
            end
    end
end

function spectrum_buoyancy(out,u,ρ,g)
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        k = ZRANGE[kk]
        j = YRANGE[jj]
            @inbounds @simd for i in XRANGE
                out[i,j,k] = vecouterproj(u[i,j,k],g*ρ[i,j,k]).z
            end
    end
end

function spectrum_forcing(out,ux,uy,fx,fy,dt)
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        k = ZRANGE[kk]
        j = YRANGE[jj]
            @inbounds @simd for i in XRANGE
                out[i,j,k] = (proj(ux[i,j,k],fx[i,j,k]) + proj(uy[i,j,k],fy[i,j,k]))/dt
            end
    end
end

function spectrum_spectral_barrier(hout,vout,u::VectorField{T},hv::SpectralBarrier{ini,cut,F}) where {T,F,ini,cut}
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        @inbounds begin
            f = hv.func
            kz2 = KZ[l]*KZ[l]
                kyz2 = muladd(KY[j], KY[j], kz2)
                for i in XRANGE
                    k2 = muladd(KX[i], KX[i], kyz2)
                    k = sqrt(k2)
                    fk = f(k)
                    fk = ifelse(fk == -Inf, 0.0,fk)
                    v = u[i,j,l]
                    out = fk*vecouterproj(v,v)
                    hout[i,j,l] = out.x + out.y
                    vout[i,j,l] = out.z
                end
        end
    end
end

function spectrum_hyperviscosity(hout,vout,u::VectorField{T},s) where {T}
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        mν = nuh(s)
        @inbounds begin
            M::Int = get_hyperviscosity_exponent(s)
            kz2 = KZ[l]*KZ[l]
                kyz2 = muladd(KY[j], KY[j], kz2)
                @simd for i in XRANGE
                    k2 = muladd(KX[i], KX[i], kyz2)
                    out = mν*(k2^M)*vecouterproj(u[i,j,l], u[i,j,l])
                    hout[i,j,l] = out.x + out.y
                    vout[i,j,l] = out.z
                end
        end
    end
end

function spectrum_u(hout,vout,u::VectorField{T}) where {T}
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
                @inbounds @simd for i in XRANGE
                    out = vecouterproj(u[i,j,l], u[i,j,l])
                    hout[i,j,l] = (out.x + out.y)/2
                    vout[i,j,l] = out.z/2
                end
    end
end

function spectrum_rho(hout,rho::ScalarField{T}) where {T}
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
                @inbounds @simd for i in XRANGE
                    out = proj(rho[i,j,l], rho[i,j,l])
                    hout[i,j,l] = out
                end
    end
end

function spectrum_viscosity(hout,vout,hin,vin,s) where {T}
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        mν = ν
            kz2 = KZ[l]*KZ[l]
                kyz2 = muladd(KY[j], KY[j], kz2)
                @inbounds @simd for i in XRANGE
                    k2 = muladd(KX[i], KX[i], kyz2)
                    hout[i,j,l] = 2*mν*k2*hin[i,j,l]
                    vout[i,j,l] = 2*mν*k2*vin[i,j,l]
                end
    end
end

function spectrum_rhodiss(hout,hin,s) where {T}
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        mν = diffusivity(s.densitystratification)
            kz2 = KZ[l]*KZ[l]
                kyz2 = muladd(KY[j], KY[j], kz2)
                @inbounds @simd for i in XRANGE
                    k2 = muladd(KX[i], KX[i], kyz2)
                    hout[i,j,l] = mν*k2*hin[i,j,l]
                end
    end
end

@par function spectrum_pressure(hout,vout,s::@par(Simulation))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        k = ZRANGE[kk]
        j = YRANGE[jj]
        spectrum_pressure(j,k,hout,vout,s)
    end
end

@par function spectrum_pressure(j::Int,k::Int,hout,vout,s::A) where {A<:@par(Simulation)}
    u = s.u.c
    rhsv = s.rhs.c
    ke = s.kspec

    if hasles(A)
        τ = s.lesmodel.tau.c
    end

    if hasdensity(A)
        ρ = s.densitystratification.ρ.field
        g = gravity(s.densitystratification)
    end

        @inbounds @msimd for i in XRANGE

            kh = K[i,j,k]
            K2 = kh⋅kh
            v = u[i,j,k]
            rhs = rhsv[i,j,k] - (im*ke[i,j,k])*kh

            if hasles(A)
                if !(!is_SandP(A) && is_FakeSmagorinsky(A))
                    rhs += (im*kh)⋅τ[i,j,k]
                end
            end

            if hasdensity(A)
                rhs += ρ[i,j,k]*g
            end

            p1 = ifelse(k==j==i==1,zero(ComplexF64),-(kh⋅rhs)/K2)
            pressure = p1*kh
            out = vecouterproj(v,pressure)

            outxy = out.x + out.y
            outz = out.z

            hout[i,j,k] = outxy
            vout[i,j,k] = outz
        end
end

function calculate_specHV(h,v,out)
    fill!(h,0.0)
    fill!(v,0.0)
    @inbounds for l in ZRANGE
        KZ2 = KZ[l]^2
        nk = round(Int, fsqrt(KZ2)/DKZ) + 1
        @inbounds for j in YRANGE
            KY2 = KY[j]^2
            nj = round(Int, fsqrt(KY2)/MAXDKH) + 1
            ee = out[1,j,l]
            h[1] += 0.5*ee
            h[nj] += 0.5*ee
            v[nk] += ee
            @msimd for i in 2:NX
                k = fsqrt(KX[i]*KX[i])

                nx = round(Int, k/MAXDKH) + 1

                ee = 2*out[i,j,l]
                h[nx] += 0.5*ee*MAXDKH
                h[nj] += 0.5*ee*MAXDKH
                v[nk] += ee*DKZ

            end
        end
    end
end

function filter_spectrum!(hout,vout)
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        @inbounds begin
            l = ZRANGE[kk]
            j = YRANGE[jj]
            KZ2 = KZ[l]^2
                KY2 = KY[j]^2
                KYZ2 = KZ2 + KY2
                @simd for i in XRANGE
                    KX2 = KX[i]*KX[i]
                    K2 = KYZ2 + KX2
                    #G2 = Gaussfilter(SPEC_FIL_D2,K2) * ((KX2 <= SPEC_FIL_CUTOFF) & (KY2 <= SPEC_FIL_CUTOFF) & (KZ2 <= SPEC_FIL_CUTOFF))
                    G2 = Gaussfilter(SPEC_FIL_D2,K2)
                    G2 *=G2
                    hout[i,j,l] = G2*hout[i,j,l]
                    vout[i,j,l] = G2*vout[i,j,l]
                end
        end
    end
end


text:
#=
Tools for working with Markov Chains

@author : Spencer Lyon, Zac Cranko

@date: 07/10/2014

References
----------

http://quant-econ.net/jl/finite_markov.html
=#
import LightGraphs: DiGraph, period, attracting_components,
                    strongly_connected_components, is_strongly_connected

@inline check_stochastic_matrix(P) = maxabs(sum(P, 2) - 1) < 1e-15 ? true : false

"""
Finite-state discrete-time Markov chain.

It stores useful information such as the stationary distributions, and
communication, recurrent, and cyclic classes, and allows simulation of state
transitions.

##### Fields

- `p::Matrix` The transition matrix. Must be square, all elements must be
positive, and all rows must sum to unity
"""
type MarkovChain{T, TM<:AbstractMatrix, TV<:AbstractVector}
    p::TM # valid stochastic matrix
    state_values::TV

    function MarkovChain(p::AbstractMatrix, state_values)
        n, m = size(p)

        eltype(p) != T &&
            throw(ArgumentError("Types must be consistent with given matrix"))

        n != m &&
            throw(DimensionMismatch("stochastic matrix must be square"))

        any(p.<0) &&
            throw(ArgumentError("stochastic matrix must have nonnegative elements"))

        !check_stochastic_matrix(p) &&
            throw(ArgumentError("stochastic matrix rows must sum to 1"))

        length(state_values) != n &&
            throw(DimensionMismatch("state_values should have $n elements"))

        return new(p, state_values)
    end
end

# Provide constructor that infers T from eltype of matrix
MarkovChain(p::AbstractMatrix, state_values=1:size(p, 1)) =
    MarkovChain{eltype(p), typeof(p), typeof(state_values)}(p, state_values)

"Number of states in the markov chain `mc`"
n_states(mc::MarkovChain) = size(mc.p, 1)

function Base.show{T,TM}(io::IO, mc::MarkovChain{T,TM})
    println(io, "Discrete Markov Chain")
    println(io, "stochastic matrix of type $TM:")
    println(io, mc.p)
end

"""
solve x(P-I)=0 using an algorithm presented by Grassmann-Taksar-Heyman (GTH)

##### Arguments

- `p::Matrix` : valid stochastic matrix

##### Returns

- `x::Matrix`: A matrix whose columns contain stationary vectors of `p`

##### References

The following references were consulted for the GTH algorithm

- W. K. Grassmann, M. I. Taksar and D. P. Heyman, "Regenerative Analysis and
Steady State Distributions for Markov Chains, " Operations Research (1985),
1107-1116.
- W. J. Stewart, Probability, Markov Chains, Queues, and Simulation, Princeton
University Press, 2009.

"""
gth_solve{T<:Integer}(a::AbstractMatrix{T}) =
    gth_solve(convert(AbstractArray{Rational, 2}, a))

# solve x(P-I)=0 by the GTH method
function gth_solve{T<:Real}(original::AbstractMatrix{T})
    # NOTE: the full here will convert a sparse matrix to dense before the gth
    #       algorithm begins. Given the nature of the algorithm this will
    #       almost certainly be more efficient than operating on the sparse
    #       matrix itself.
    a = copy(full(original))
    n = size(a, 1)
    x = zeros(T, n)

    @inbounds for k in 1:n-1
        scale = sum(a[k, k+1:n])
        if scale <= zero(T)
            # There is one (and only one) recurrent class contained in
            # {1, ..., k};
            # compute the solution associated with that recurrent class.
            n = k
            break
        end
        a[k+1:n, k] /= scale

        for j in k+1:n, i in k+1:n
            a[i, j] += a[i, k] * a[k, j]
        end
    end

    # backsubstitution
    x[n] = 1
    @inbounds for k in n-1:-1:1, i in k+1:n
        x[k] += x[i] * a[i, k]
    end

    # normalisation
    x / sum(x)
end

"""
Find the recurrent classes of the `MarkovChain`

##### Arguments

- `mc::MarkovChain` : MarkovChain instance containing a valid stochastic matrix

##### Returns

- `x::Vector{Vector}`: A `Vector` containing `Vector{Int}`s that describe the
recurrent classes of the transition matrix for p

"""
recurrent_classes(mc::MarkovChain) = attracting_components(DiGraph(mc.p))

communication_classes(mc::MarkovChain) = strongly_connected_components(DiGraph(mc.p))

is_irreducible(mc::MarkovChain) =  is_strongly_connected(DiGraph(mc.p))

is_aperiodic(mc::MarkovChain) = period(mc) == 1

function period(mc::MarkovChain)
    g = DiGraph(mc.p)
    recurrent = attracting_components(g)
    periods   = Int[]

    for r in recurrent
        push!(periods, period(g[r]))
    end
    return lcm(periods)
end
"""
calculate the stationary distributions associated with a N-state markov chain

##### Arguments

- `mc::MarkovChain` : MarkovChain instance containing a valid stochastic matrix

##### Returns

- `dists::Matrix{Float64}`: N x M matrix where each column is a stationary
distribution of `mc.p`

"""
function mc_compute_stationary{T}(mc::MarkovChain{T})
    recurrent = recurrent_classes(mc)

    # unique recurrent class
    length(recurrent) == 1 && return gth_solve(mc.p)

    # more than one recurrent classes
    stationary_dists = Array(T, n_states(mc), length(recurrent))
    for (i, class) in enumerate(recurrent)
        dist        = zeros(T, n_states(mc))
        dist[class] = gth_solve(mc.p[class, class]) # recast to correct dimension
        stationary_dists[:, i] = dist
    end
    return stationary_dists
end

"""
Simulate time series of state transitions of the Markov chain `mc`.

The sample path from the `j`-th repetition of the simulation with initial state
`init[i]` is stored in the `(j-1)*num_reps+i`-th column of the matrix X.

##### Arguments

- `mc::MarkovChain` : MarkovChain instance.
- `ts_length::Int` : Length of each simulation.
- `init::Vector{Int}` : Vector containing initial states.
- `;num_reps::Int(1)` : Number of repetitions of simulation for each element
of `init`

##### Returns

- `X::Matrix{Int}` : Array containing the sample paths as columns, of shape
(ts_length, k), where k = length(init)* num_reps

"""
function simulate(mc::MarkovChain, ts_length::Int, init::Vector{Int};
                  num_reps::Int=1)
    k = length(init)*num_reps
    X = Array(Int, ts_length, k)
    X[1, :] = repmat(init, num_reps)

    simulate!(mc, X)
    return X
end

"""
Simulate time series of state transitions of the Markov chain `mc`.

##### Arguments

- `mc::MarkovChain` : MarkovChain instance.
- `ts_length::Int` : Length of each simulation.
- `init::Int` : Initial state.
- `;num_reps::Int(1)` : Number of repetitions of simulation

##### Returns

- `X::Matrix{Int}` : Array containing the sample paths as columns, of shape
(ts_length, k), where k = num_reps

"""
simulate(mc::MarkovChain, ts_length::Int, init::Int; num_reps::Int=1) =
    simulate(mc, ts_length, [init], num_reps=num_reps)

"""
Simulate time series of state transitions of the Markov chain `mc`.

##### Arguments

- `mc::MarkovChain` : MarkovChain instance.
- `ts_length::Int` : Length of each simulation.
- `;num_reps::Union{Int, Void}(nothing)` : Number of repetitions of simulation.

##### Returns

- `X::Matrix{Int}` : Array containing the sample paths as columns, of shape
(ts_length, k), where k = num_reps

"""
function simulate(mc::MarkovChain, ts_length::Int; num_reps::Int=1)
    init = rand(1:size(mc.p, 1), num_reps)
    return simulate(mc, ts_length, init)
end

"""
Fill `X` with sample paths of the Markov chain `mc` as columns.

##### Arguments

- `mc::MarkovChain` : MarkovChain instance.
- `X::Matrix{Int}` : Preallocated matrix of integers to be filled with sample
paths of the markov chain `mc`. The elements in `X[1, :]` will be used as the
initial states.

"""
function simulate!(mc::MarkovChain, X::Matrix{Int})
    n = size(mc.p, 1)

    # NOTE: ensure dense array and transpose before slicing the array. Then
    #       when passing to DiscreteRV use `sub` to avoid allocating again
    p = full(mc.p)'
    P_dist = [DiscreteRV(sub(p, :, i)) for i in 1:n]

    ts_length, k = size(X)

    for i in 1:k
        for t in 1:ts_length-1
            X[t+1, i] = draw(P_dist[X[t]])
        end
    end
    X
end

"""
Simulate time series of state transitions of the Markov chain `mc`.

##### Arguments

- `mc::MarkovChain` : MarkovChain instance.
- `ts_length::Int` : Length of each simulation.
- `init_state::Int(rand(1:n_states(mc)))` : Initial state.

##### Returns

- `x::Vector`: A vector of transition indices for a single simulation
"""
function simulation(mc::MarkovChain, ts_length::Int,
                    init_state::Int=rand(1:n_states(mc)))
    simulate(mc, ts_length, init_state; num_reps=1)[:, 1]
end

# ----------------------- #
# simulate_values methods #
# ----------------------- #

function simulate_values(mc::MarkovChain,
                         ts_length::Int,
                         init_state::Vector{Int};
                         num_reps::Int=1)
    k = length(init_state)*num_reps
    X = Array(eltype(mc.state_values), ts_length, k)
    init_state = repmat(init_state, num_reps)
    X[1, :] = mc.state_values[init_state]

    simulate_values!(mc, init_state, X)
    return X
end

function simulate_values(mc::MarkovChain, ts_length::Int, init_state::Int;
                         num_reps::Int=1)
    simulate_values(mc, ts_length, [init_state], num_reps=num_reps)
end

function simulate_values(mc::MarkovChain, ts_length::Int; num_reps::Int=1)
    init_state = rand(1:size(mc.p, 1), num_reps)
    simulate_values(mc, ts_length, init_state)
end

function simulate_values!(mc::MarkovChain, init_state::Vector{Int}, X::Matrix)
    n = size(mc.p, 1)

    # NOTE: ensure dense array and transpose before slicing the array. Then
    #       when passing to DiscreteRV use `sub` to avoid allocating again
    p = full(mc.p)'
    P_dist = [DiscreteRV(sub(p, :, i)) for i in 1:n]

    ts_length, k = size(X)

    for i in 1:k
        i_state = init_state[i]
        for t in 1:ts_length-1
            i_state = draw(P_dist[i_state])
            X[t+1, i] = mc.state_values[i_state]
        end
    end
    X
end

@doc """ Like `simulate(::MarkovChain, args...; kwargs...)`, but instead of
returning integers specifying the state indices, this routine returns the
values of the `mc.state_values` at each of those indices. See docstring
for `simulate` for more information
""" simulate_values, simulate_values!

"""
Simulate time series of state transitions of the Markov chain `mc`.

##### Arguments

- `mc::MarkovChain` : MarkovChain instance.
- `ts_length::Int` : Length of each simulation.
- `init_state::Int(rand(1:n_states(mc)))` : Initial state.

##### Returns

- `x::Vector`: A vector of state values along a simulated path
"""
function value_simulation(mc::MarkovChain, ts_length::Int,
                          init_state::Int=rand(1:n_states(mc)))
    simulate_values(mc, ts_length, init_state; num_reps=1)[:, 1]
end

ext: jl


text:
@doc """ read data from Homo_sapiens.gene_info
""" ->
function read_data()
    if !isdir(data_dir)
        mkdir(data_dir)
    end
    if !isfile(hs_fl)
        cur_dir =  pwd()
        cd(data_dir)
        cmd = `wget ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/GENE_INFO/Mammalia/Homo_sapiens.gene_info.gz`
        run(cmd)
        #download("ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/GENE_INFO/Mammalia/Homo_sapiens.gene_info.gz")
        run(`gunzip $(string(hs_fl,".gz"))`)
        cd(cur_dir)
    end
    @assert isfile(hs_fl)
    data = readdlm(hs_fl,ASCIIString)::Array{ASCIIString,2}
end

@doc """ get gene name and its synonym from data,
         save it to file, and return it
""" ->
function genename_synonym(data::Array{ASCIIString,2})
    num_gene = size(data,1)
    genenames = data[:,3]
    synonyms  = map(x->map(y->convert(ASCIIString,y), split(x,"|")), data[:,5])

    id_genes = Dict{Int64, Set{ASCIIString}}()
    for i = 1:num_gene
        id_genes[i] = filter(x->x!="-", Set(vcat(genenames[i], synonyms[i])))
    end

    save(id_genes_fl, "id_genes_dict", id_genes)

    genes_id_pair = map(x->(collect(x[1]),x[2]), map(reverse, collect(id_genes)))
    gene_id = Dict{ASCIIString,Int64}()
    for i = 1:num_gene
        samplegenes = genes_id_pair[i][1]
        id          = genes_id_pair[i][2]
        num_samplegene = length(samplegenes)
        for j = 1:num_samplegene
            genename = samplegenes[j]
            gene_id[genename] = id
        end
    end
    save(gene_id_fl, "gene_id_dict", gene_id)

    nothing
end

function build_index_synonym()
    data = read_data()
    genename_synonym(data)
end


text:
"""
Finite set {0,1} × ... × {0,1}
"""
struct BooleanCube{S} <: AbstractProduct
    itr :: UnitRange{S}
    n :: Int
end
const BoolCube = BooleanCube
const BCube = BooleanCube

"""
construct a boolean cube with given number of dimensions
"""
BCube(n :: Int) = BCube(0:(2^n - 1), n)

# INTERFACE
iterator(bc :: BCube) = bc
ndims(bc :: BCube) = bc.n
length(bc :: BCube) = 2^ndims(bc)
factors(bc :: BCube) = tuple(fill(BCube(1), ndims(bc))...)

# Iteration utilities
postprocess(:: Nothing, n) = nothing
postprocess(p, n) = (BoolVec(p[1], n), p[2])

function iterate(bc :: BCube)
    p = iterate(bc.itr)
    return postprocess(p, ndims(bc))
end

function iterate(bc :: BCube, state)
    p = iterate(bc.itr, state)
    return postprocess(p, ndims(bc))
end

# direct product of two bool cubes
direct(bc1 :: BCube, bc2 :: BCube) = BCube(ndims(bc1) + ndims(bc2))

# show utilities
Base.show(io :: IO, bc :: BCube) = print(io, "{0,1}^$(ndims(bc))")

#### Boolean vectors

struct BooleanVector
    k :: Int
    n :: Int
end

const BoolVec = BooleanVector
const BooleanVec = BooleanVector

BoolVec(vs :: Union{NTuple, Vector}) = BoolVec(itr2bin(vs), length(vs))
BoolVec(v :: BoolVec) = v
value(v :: BoolVec) = v.k
value(x) = x   # fallback

function getindex(v :: BoolVec, i :: Int)
    n = ndims(v)
    return (value(v) >> (n - i)) & 1
end

ndims(v :: BoolVec) = v.n
length(v :: BoolVec) = ndims(v)
iterate(v :: BoolVec) = (v[1], 2)
iterate(v :: BoolVec, i :: Int) = i > ndims(v) ? nothing : (v[i], i+1)

dot(v1 :: BoolVec, v2 :: BoolVec) = sum(collect(v1) .* collect(v2)) % 2
# v1 <= v2 in lexicographic order?
isless(v1 :: BoolVec, v2 :: BoolVec) = (value(v1) & value(v2) == value(v1))

function Base.show(io :: IO, v :: BoolVec)
    n = ndims(v) - 1
    s = bitstring(value(v))[end - n : end]
    print(io, "[", s, ']')
end

# change representation to BoolVec
# (bc :: BCube)(x :: Union{NTuple, Array}) = BoolVec(x)
# (bc :: BCube)(x :: Int) = x

max_stars_repo_name: kirtsar/DiscreteFunctions

max_issues_repo_name: kirtsar/DiscreteFunctions


text:
# Define abstract interface
"""
$(TYPEDEF)

Abstract supertype for receive coils.
"""
abstract type ReceiveCoils{T,N} end

"""
$(SIGNATURES)

Initialize an receive coil setup with `samplingRate` [Hz] and multiple receive
coils specified by their coil sensitivity [ToDo].
`excitationFields`.
"""
function ReceiveCoils(samplingRate::Real, sensitivities::U) where {T,U<:Tuple{Vararg{W} where W<:MagneticFields{T}}}
    error("ExcitationField is not implemented for selectionField::$(typeof(selectionField)) and excitationFields::$(typeof(excitationFields)).")
end

function ReceiveCoils(samplingRate::Real, sensitivities...)
    return ReceiveCoils(samplingRate,tuple(sensitivities...))
end

"""
$(SIGNATURES)

Return a receive coil sensitivities [??] as SMatrix{N,3,T} for a given location
r [m] and time t [s].
"""
function receiveCoilSensitivities(receiveCoil::ReceiveCoils{T,N},r::SVector{3,T},t::T) where {T,N}
    error("receiveCoilSensitivities is not implemented for receiveCoil::$(typeof(receiveCoil)).")
end

# define specific subtypes
## N-dimensional receive coil setup all running with the same syncronized clock
struct StaticReceiveCoils{T,N,V<:Tuple{Vararg{W} where W<:StaticMagneticFields{T}},U} <: ReceiveCoils{T,N}
    samplingRate::U
    sensitivities::V
end

function receiveCoilSensitivities(receiveCoil::StaticReceiveCoils{T,N},r::SVector{3,T},t::T) where {T,N}
    m = MMatrix{N,3,T}(undef)
    @inbounds for i =1:N
        m[i,:] = magneticFieldStrength(receiveCoil.sensitivities[i],r,t)
    end
    return SMatrix{N,3,T}(m)
end

function ReceiveCoils(samplingRate::Real,sensitivities::V) where {T,V<:Tuple{Vararg{W} where W<:StaticMagneticFields{T}}}
    N = length(sensitivities)
    return StaticReceiveCoils{T,N,V,typeof(samplingRate)}(samplingRate,sensitivities)
end


text:
using DelimitedFiles

function sum_sliding_window(depths, i)
  depths[i] + depths[i+1] + depths[i+2]
end

no_of_depth_increases = 0

depths = readdlm("input.txt", '\n', Int)
  for i = 2:length(depths)-2 # to avoid index error
    if sum_sliding_window(depths, i) > sum_sliding_window(depths, i-1)
      global no_of_depth_increases += 1
    end
  end
return no_of_depth_increases

print("Total Number of depth increases is: ", no_of_depth_increases)


text:
import Base: readbytes!, close
let
    global default_streamer, set_default_streamer, readbytes!, close
    _default_streamer::Tuple{Function, String}                = (recorder, DEFAULT_RECORDER_ID)
    default_streamer()::Union{Base.Process,PyObject,IOBuffer} = _default_streamer[1](_default_streamer[2])
    set_default_streamer(streamerkind::Function, id::String)  = (_default_streamer = (streamerkind, id); return)

    function readbytes!(stream::PyObject, b::AbstractVector{UInt8}, nb=length(b))
        nb_frames = Int(nb/sizeof(AUDIO_ELTYPE))
        if stream.__class__ == Wave.Wave_read
            tmp = Vector{UInt8}(stream.readframes(nb_frames))
            b[1:length(tmp)] .= tmp
            return length(tmp)
        elseif stream.__class__ == Sounddevice.RawInputStream
            b .= stream.read(nb_frames)[1]
            return nb
        else
            @APIUsageError("PyObject of unknown class.")
        end
    end

    close(stream::PyObject) = stream.close()
end


text:
using ConstraintDomains
using QUBOConstraints
using Test

@testset "QUBOConstraints.jl" begin
    include("encoding.jl")
    include("learn.jl")
end


text:
"""
Using Grunwald Letnikov discretization operator to discrete fractional operator.
"""
struct GLDiff <: FractionalDiffEqAlgorithm end

function solve(α, d, rightfun, M, N, initial_condition, left_boundry, right_boundry, ::GLDiff)
    ue=zeros(M+1, N+1)
    u=copy(ue)
    D=zeros(M-1, 1)
    a=copy(D)
    #

    #exact(x, t) = exp(-t).*x.^3
    #=
    for k=1:N+1
        ue[1:end, k]=exact(x[1:end], t[k])
    end
    =#
    u[1:end, 1]=initial_condition.(x)
    u[1, 1:end]=left_boundry.(t)
    u[end, 1:end]=right_boundry.(t)
    A=zeros(M-1, M-1)
    for i=1:M-1
        D[i, 1]=d(x[i+1])
    end
    a=tau*D/(2*h^α)
    gg = g(M, α)
    for i=1:M-1
        for k=1:N-1
            if k <= i-1
                A[i, k]=a[i, 1]*gg[i-k+2, 1]
            elseif k==i
                A[i, k]=a[i, 1]*gg[2, 1]
            elseif k==i+1
                A[i, k]=a[i, 1]*gg[1, 1]
            else
                A[i, k]=0
            end
        end
    end

    for k=1:N
        b = (zeros(M-1, M-1)+I+A)*u[2:end-1, k].+tau.*rightfun.(x[2:end-1], t[k]+tau/2) .+a.*(gg[3:end].*(u[1, k+1]+u[1, k]).*ones(M-1, 1))
        b[end, 1]=b[end, 1]+a[end, 1]*gg[1]*(u[end, k+1]+u[end, k])
        u[2:end-1, k+1]=(zeros(M-1, M-1)+I-A)\b
    end
    return u
end

function g(M, α)
    gg_alpha = zeros(M+1, 1)
    gg_alpha[1, 1]=1
    for i=1:M
        gg_alpha[i+1, 1]=gamma(i-α)/(gamma(-α)*gamma(i+1))
    end
    return gg_alpha
end


text:
import Base: isapprox

"""
    PauliBasis(num_qubits::Int)

Basis for an N-qubit space where `num_qubits` specifies the number of qubits.
The dimension of the basis is 2²ᴺ.
"""
struct PauliBasis{S,B} <: Basis
    shape::S
    bases::B
    function PauliBasis(num_qubits::T) where {T<:Integer}
        shape = [2 for _ in 1:num_qubits]
        bases = Tuple(SpinBasis(1//2) for _ in 1:num_qubits)
        return new{typeof(shape),typeof(bases)}(shape, bases)
    end
end
==(pb1::PauliBasis, pb2::PauliBasis) = length(pb1.bases) == length(pb2.bases)

"""
    Base class for Pauli transfer matrix classes.
"""
abstract type PauliTransferMatrix{B1, B2} end

"""
    DensePauliTransferMatrix(B1, B2, data)

DensePauliTransferMatrix stored as a dense matrix.
"""
mutable struct DensePauliTransferMatrix{B1,B2,T<:Matrix} <: PauliTransferMatrix{B1, B2}
    basis_l::B1
    basis_r::B2
    data::T
    function DensePauliTransferMatrix(basis_l::BL, basis_r::BR, data::T) where {BL,
                                                                                BR,
                                                                                T<:Matrix}
        if length(basis_l[1])*length(basis_l[2]) != size(data, 1) ||
           length(basis_r[1])*length(basis_r[2]) != size(data, 2)
            throw(DimensionMismatch())
        end
        new{BL, BR, T}(basis_l, basis_r, data)
    end
end

PauliTransferMatrix(ptm::DensePauliTransferMatrix) = ptm

function *(ptm0::DensePauliTransferMatrix{B, B},ptm1::DensePauliTransferMatrix{B, B}) where B
    return DensePauliTransferMatrix(ptm0.basis_l, ptm1.basis_r, ptm0.data*ptm1.data)
end

"""
    Base class for χ (process) matrix classes.
"""
abstract type ChiMatrix{B1, B2} end

"""
    DenseChiMatrix(b, b, data)

DenseChiMatrix stored as a dense matrix.
"""
mutable struct DenseChiMatrix{B1,B2,T<:Matrix} <: PauliTransferMatrix{B1, B2}
    basis_l::B1
    basis_r::B2
    data::T
    function DenseChiMatrix(basis_l::BL, basis_r::BR, data::T) where {BL,BR,T<:Matrix}
        if length(basis_l[1])*length(basis_l[2]) != size(data, 1) ||
           length(basis_r[1])*length(basis_r[2]) != size(data, 2)
            throw(DimensionMismatch())
        end
        new{BL, BR, T}(basis_l, basis_r, data)
    end
end

ChiMatrix(chi_matrix::DenseChiMatrix) = chi_matrix

"""
A dictionary that represents the Pauli algebra - for a pair of Pauli operators
σᵢσⱼ information about their product is given under the key "ij". The first
element of the dictionary value is the Pauli operator, and the second is the
scalar multiplier. For example, σ₀σ₁ = σ₁, and `"01" => ("1", 1)`.
"""
const pauli_multiplication_dict = Dict(
  "00" => ("0", 1.0+0.0im),
  "23" => ("1", 0.0+1.0im),
  "30" => ("3", 1.0+0.0im),
  "22" => ("0", 1.0+0.0im),
  "21" => ("3", -0.0-1.0im),
  "10" => ("1", 1.0+0.0im),
  "31" => ("2", 0.0+1.0im),
  "20" => ("2", 1.0+0.0im),
  "01" => ("1", 1.0+0.0im),
  "33" => ("0", 1.0+0.0im),
  "13" => ("2", -0.0-1.0im),
  "32" => ("1", -0.0-1.0im),
  "11" => ("0", 1.0+0.0im),
  "03" => ("3", 1.0+0.0im),
  "12" => ("3", 0.0+1.0im),
  "02" => ("2", 1.0+0.0im),
)

"""
    multiply_pauli_matirices(i4::String, j4::String)

A function to algebraically determine result of multiplying two
(N-qubit) Pauli matrices. Each Pauli matrix is represented by a string
in base 4. For example, σ₃⊗σ₀⊗σ₂ would be "302". The product of any pair of
Pauli matrices will itself be a Pauli matrix multiplied by any of the 1/4 roots
of 1.
"""
cache_multiply_pauli_matrices() = begin
    local pauli_multiplication_cache = Dict()
    function _multiply_pauli_matirices(i4, j4)
        if (i4, j4) ∉ keys(pauli_multiplication_cache)
            pauli_multiplication_cache[(i4, j4)] = mapreduce(x -> pauli_multiplication_dict[prod(x)],
                                                             (x,y) -> (x[1] * y[1], x[2] * y[2]),
                                                             zip(i4, j4))
        end
        return pauli_multiplication_cache[(i4, j4)]
    end
end
multiply_pauli_matirices = cache_multiply_pauli_matrices()

function *(chi_matrix0::DenseChiMatrix{B, B},chi_matrix1::DenseChiMatrix{B, B}) where B

    num_qubits = length(chi_matrix0.basis_l[1].shape)
    sop_dim = 2 ^ prod(chi_matrix0.basis_l[1].shape)
    ret = zeros(ComplexF64, (sop_dim, sop_dim))

    for ijkl in Iterators.product(0:(sop_dim-1),
                                  0:(sop_dim-1),
                                  0:(sop_dim-1),
                                  0:(sop_dim-1))
        i, j, k, l = ijkl
        if (chi_matrix0.data[i+1, j+1] != 0.0) & (chi_matrix1.data[k+1, l+1] != 0.0)
            i4, j4, k4, l4 = map(x -> string(x, base=4, pad=2), ijkl)

            pauli_product_ik = multiply_pauli_matirices(i4, k4)
            pauli_product_lj = multiply_pauli_matirices(l4, j4)

            ret[parse(Int, pauli_product_ik[1], base=4)+1,
                parse(Int, pauli_product_lj[1], base=4)+1] += (pauli_product_ik[2] * pauli_product_lj[2] * chi_matrix0.data[i+1, j+1] * chi_matrix1.data[k+1, l+1])
        end
    end
    return DenseChiMatrix(chi_matrix0.basis_l, chi_matrix0.basis_r, ret / 2^num_qubits)
end

# TODO MAKE A GENERATOR FUNCTION
"""
    pauli_operators(num_qubits::Integer)

Generate a list of N-qubit Pauli operators.
"""
function pauli_operators(num_qubits::Integer)
    pauli_funcs = (identityoperator, sigmax, sigmay, sigmaz)
    po = []
    for paulis in Iterators.product((pauli_funcs for _ in 1:num_qubits)...)
        basis_vector = reduce(⊗, f(SpinBasis(1//2)) for f in paulis)
        push!(po, basis_vector)
    end
    return po
end

"""
    pauli_basis_vectors(num_qubits::Integer)

Generate a matrix of basis vectors in the Pauli representation given a number
of qubits.
"""
function pauli_basis_vectors(num_qubits::Integer)
    po = pauli_operators(num_qubits)
    sop_dim = 4 ^ num_qubits
    return mapreduce(x -> sparse(reshape(x.data, sop_dim)), (x, y) -> [x y], po)
end

"""
    PauliTransferMatrix(sop::DenseSuperOpType)

Convert a superoperator to its representation as a Pauli transfer matrix.
"""
function PauliTransferMatrix(sop::DenseSuperOpType)
    num_qubits = length(sop.basis_l[1].bases)
    pbv = pauli_basis_vectors(num_qubits)
    sop_dim = 4 ^ num_qubits
    data = real.(pbv' * sop.data * pbv / √sop_dim)
    return DensePauliTransferMatrix(sop.basis_l, sop.basis_r, data)
end

SuperOperator(unitary::DenseOpType) = spre(unitary) * spost(unitary')
SuperOperator(sop::DenseSuperOpType) = sop

"""
    SuperOperator(ptm::DensePauliTransferMatrix)

Convert a Pauli transfer matrix to its representation as a superoperator.
"""
function SuperOperator(ptm::DensePauliTransferMatrix)
    num_qubits = length(ptm.basis_l[1].bases)
    pbv = pauli_basis_vectors(num_qubits)
    sop_dim = 4 ^ num_qubits
    data = pbv * ptm.data * pbv' / √sop_dim
    return DenseSuperOperator(ptm.basis_l, ptm.basis_r, data)
end

"""
    PauliTransferMatrix(unitary::DenseOpType)

Convert an operator, presumably a unitary operator, to its representation as a
Pauli transfer matrix.
"""
PauliTransferMatrix(unitary::DenseOpType) = PauliTransferMatrix(SuperOperator(unitary))

"""
    ChiMatrix(unitary::DenseOpType)

Convert an operator, presumably a unitary operator, to its representation as a χ matrix.
"""
function ChiMatrix(unitary::DenseOpType)
    num_qubits = length(unitary.basis_l.bases)
    pbv = pauli_basis_vectors(num_qubits)
    aj = pbv' * reshape(unitary.data, 4 ^ num_qubits)
    return DenseChiMatrix((unitary.basis_l, unitary.basis_l), (unitary.basis_r, unitary.basis_r), aj * aj' / (2 ^ num_qubits))
end

"""
    ChiMatrix(sop::DenseSuperOpType)

Convert a superoperator to its representation as a Chi matrix.
"""
function ChiMatrix(sop::DenseSuperOpType{B, B, T}) where {B, T}
    num_qubits = length(sop.basis_l)
    sop_dim = 4 ^ num_qubits
    po = pauli_operators(num_qubits)
    data = Matrix{eltype(T)}(undef, (sop_dim, sop_dim))
    for (idx, jdx) in Iterators.product(1:sop_dim, 1:sop_dim)
        data[idx, jdx] = tr((spre(po[idx]) * spost(po[jdx])).data' * sop.data) / √sop_dim
    end
    return DenseChiMatrix(sop.basis_l, sop.basis_r, data)
end

"""
    PauliTransferMatrix(chi_matrix::DenseChiMatrix)

Convert a χ matrix to its representation as a Pauli transfer matrix.
"""
function PauliTransferMatrix(chi_matrix::DenseChiMatrix{B, B, T}) where {B, T}
    num_qubits = length(chi_matrix.basis_l)
    sop_dim = 4 ^ num_qubits
    po = pauli_operators(num_qubits)
    data = Matrix{real(eltype(T))}(undef, (sop_dim, sop_dim))
    for (idx, jdx) in Iterators.product(1:sop_dim, 1:sop_dim)
        data[idx, jdx] = tr(mapreduce(x -> po[idx] * po[x[1]] * po[jdx] * po[x[2]] * chi_matrix.data[x[1], x[2]],
                                      +,
                                      Iterators.product(1:16, 1:16)).data) / sop_dim |> real
    end
    return DensePauliTransferMatrix(chi_matrix.basis_l, chi_matrix.basis_r, data)
end

"""
    SuperOperator(chi_matrix::DenseChiMatrix)

Convert a χ matrix to its representation as a superoperator.
"""
SuperOperator(chi_matrix::DenseChiMatrix) = SuperOperator(PauliTransferMatrix(chi_matrix))

"""
    ChiMatrix(ptm::DensePauliTransferMatrix)

Convert a Pauli transfer matrix to its representation as a χ matrix.
"""
ChiMatrix(ptm::DensePauliTransferMatrix) = ChiMatrix(SuperOperator(ptm))

"""Equality for all varieties of superoperators."""
==(sop1::T, sop2::T) where T<:Union{DensePauliTransferMatrix, DenseSuperOpType, DenseChiMatrix} = sop1.data == sop2.data
==(sop1::Union{DensePauliTransferMatrix, DenseSuperOpType, DenseChiMatrix}, sop2::Union{DensePauliTransferMatrix, DenseSuperOpType, DenseChiMatrix}) = false

"""Approximate equality for all varieties of superoperators."""
function isapprox(sop1::T, sop2::T; kwargs...) where T<:Union{DensePauliTransferMatrix, DenseSuperOpType, DenseChiMatrix}
    return isapprox(sop1.data, sop2.data; kwargs...)
end


text:
"""
---
Apply WOLA-OFDM demodulator to input signal and returns the T/F QAM matrix
# ---
Syntax
qamRx	= wolaSigDecode(sigId,nFFT,nCP,allocSubcarriers);
# ---  Input parameters
- sigId	  : OFDM signal in time domain [Array{Complex{Float64},nbEch}]
- nbEch	: Number of samples: nbSymb*(nFFT+nCP)
- nFFT	  : FFT size [Int]
- nCP	  : Cyclic prefix size (in samples) [Int]
- allocSubcarriers : Vector of index of allocated subcarriers [Array{Int,nbSubcarriers}]
- winLengthTx	: Window size @Tx side (FFT rotation)
- doTailBiting	: Do tail biting approach (by default 1) [Int]
- winFunc : Type of window. Can be a string of supported window or directly the window taps.
- winLength : Length of window. This parameter is not used if winFunc is an array of the window tap.
# ---  Output parameters
- qamMat  : Time frequency matrix : [Array{Complex{Float64},nbSubcarriers,nbSymb}]
# ---
# Supported window
-	"Triangle"		: Triangle window
-	"srrc"			: Square Root Raised Cosine
-	"Meyer"			: Meyer window (See [1])
# ---
# References
- [1] R. Zayani, Y. Medjahdi, H. Shaiek and D. Roviras, "WOLA-OFDM: A Potential Candidate for Asynchronous 5G," 2016.
- [2] Y. Medjahdi and al, "On the road to 5G: Comparative study of Physical layer in MTC context", 2017.
- [3] R. Gerzaguet and al, "Comparison of Promising Candidate Waveforms for 5G: WOLA-OFDM Versus BF-OFDM", 2017.
# ---
# v 1.0 - Robin Gerzaguet.
"""
function wolaSigDecode(sigChan,nFFT,nCP,allocatedSubcarriers,winLengthTx,doTailBiting,winFunc::String,winLength=0,window=[])
	# Set 4 core for FFT computation
	#FFTW.set_num_threads(4)
	# --- WOLA signal parameters
	sizeSymb		= nFFT + nCP;
	# ----------------------------------------------------
	# --- Getting window
	# ----------------------------------------------------
	if isempty(window)
		# --- We create a window with a dedicated function
		# No control on window and parameters, it will be done in the function
		objWindow = getWolaWindow(winFunc,nFFT,0,winLength);
		window	  = objWindow.window;
	end
	# --- Rotation applied after FFT
	# As we have overlapping between symbol, FFT position is set @ middle (and not at the end
	# as for CP-OFDM)
	# --> We need to calculate this position (function of Tx side)
	# --> We need to apply a rotation factor
	#delta           = (nFFT+nCP)÷2 - (nFFT + winLength)÷2;
	delta            = winLengthTx÷2 + nCP + nFFT -nFFT - winLength;	# End of apod.
	rotFactor       = -winLengthTx÷2 - nCP + delta +winLength÷2 +1;
	# --- Getting number of symbol to decode
	nbSymbDecode   = Int(floor((length(sigChan)-winLength)/sizeSymb));
	# --- Initiate Rx T/F matrix
	qamRx		   = zeros(Complex{Float64},length(allocatedSubcarriers),nbSymbDecode);
	sigW		   = zeros(Complex{Float64},nFFT + winLength);
	# --- OFDM demodulator
	for n = 1  : 1 : nbSymbDecode
		# --- Getting current extnded symbol
		currSymb  = sigChan[ 1+delta + (n-1)*sizeSymb .+ (1:nFFT+winLength)];
		# --- Rx windowing
		currSymb = currSymb .* window;
		# --- Tail biting
		if Bool(doTailBiting)
			# --- Append end to begining
			sigW[ 1:winLength ] = currSymb[ end+1-winLength:end ];
			# --- Append beginning to end
			sigW[ end+1-winLength:end ] = currSymb[ 1:winLength ];
			# --- Add
			sigTB	  = currSymb + sigW;
		else
			sigTB     = currSymb;
		end
		# -- Removing CP
		currSymb  = sigTB[ 1+winLength÷2:end-winLength÷2 ]
		# --- FFT and rotation
		fftEl	  = fft(currSymb,1) .* exp.(-2im*pi*rotFactor/nFFT.*(0:nFFT-1));
		# --- Mapping in output matrix
		qamRx[:,n]		= fftEl[allocatedSubcarriers];
	end
	return qamRx;
end

function wolaSigDecode(sigChan,wola::StrucWOLA,doTailBiting=1);
	return wolaSigDecode(sigChan,wola.nFFT,wola.nCP,wola.allocatedSubcarriers,wola.windowTx.winLength,doTailBiting,wola.windowRx.winFunc,wola.windowRx.winLength,wola.windowRx.window);
end


text:
using Test

using FluidSynth

settings = Settings()
synth = Synth(settings)

sfont_id = sfload(synth, "violin_sample.sf2")
program_select(synth, Int32(0), sfont_id, Int32(0), Int32(0))
noteon(synth, Int32(0), Int32(60), Int32(70))
wav = write_s16_stereo(synth, Int32(1024))
@test length(wav) == 2048


text:
using CSV, DataFrames, Statistics, Plots, LaTeXStrings; pyplot()

include("simParams.jl")

function statsSummary(;doPlot = true)
    cd(@__DIR__)
    results = DataFrame(CSV.File("results.csv"))

    plotDict = Dict()
    resultsDict = Dict()

    for (j,γ) in enumerate(γVals)
        for n in nVals[j]
            tempResults =  filter([:γ, :n] => (γdf, ndf) -> γdf == γ  && ndf == n, results)
            tempResultsCG = filter(:gappingPolicy => gp -> gp,tempResults)
            tempResultsCO = filter(:gappingPolicy => gp -> !gp,tempResults)
            diffs = zeros(length(cFracVals))
            for (i,cf) in enumerate(cFracVals)
                rwdsCG = filter(:cFrac => cfv -> cfv == cf,tempResultsCG).reward
                avgCG = mean(rwdsCG)
                rwdsCO = filter(:cFrac => cfv -> cfv == cf,tempResultsCO).reward
                avgCO = mean(rwdsCO)
                diffs[i] = avgCO - avgCG
            end
            push!(plotDict,(γ,n) => diffs)
            md, i = findmax(diffs)
            cMax = cFracVals[i]
            diffsAtMax = tempResultsCO[tempResultsCO.cFrac .== cMax,:reward] - tempResultsCG[tempResultsCG.cFrac .== cMax,:reward]
            reps = length(diffsAtMax)
            serr = std(diffsAtMax)/sqrt(reps)
            # @show n, γ,  md, md-1.96serr,md+1.96serr, cMax*γ^2
            push!(resultsDict,(γ,n) => (round(md,digits=5),round(1.96serr,digits =5),round(cMax*γ^2,digits=5)))
        end
    end
    plotDict,resultsDict
end

function doPlot()
    plts = []
    doLegend = true
    for (j,γ) in enumerate(γVals)
        cVals = cFracVals*γ^2
        nnVals = nVals[j]
        p = plot(cVals,plotDict[(γ,nnVals[1])],label="n=$(nnVals[1])")
        for i in 2:length(nnVals)
            p = plot!(cVals,plotDict[(γ,nnVals[i])],
                label="n=15",title = L"\gamma = "*string(γ),
                xlabel="c",ylabel="Policy Diff",legend=doLegend)
        end
        doLegend = false
        push!(plts,p)
    end
    plot(plts...)
end

function formatTable(rd)
    gg(g,n) = rd[(g,n)][1]
    ee(g,n) = rd[(g,n)][2]
    cc(g,n) = rd[(g,n)][3]

    tableString = """
    \\begin{table}[h!]
    \\begin{center}
      \\caption{Evaluating the difference between the cool-off and call-gapping policy.}
      \\label{tab:table1}
      \\begin{tabular}{|c|c|l|l|}
        \\hline
        \\textbf{System} & \\textbf{Number of channels} & \\textbf{Maximal Gap}& \\textbf{Worst cost} \\\\
        \\hline
        \\multirow{18}{*}{ \$\\gamma = 0.2 \$}
        &  \$n=3 \$ &  \$$(gg(0.2,3)) \\pm $(ee(0.2,3)) \$ &  \$c = $(cc(0.2,3)) \$\\\\
        &  \$n=4 \$ &  \$$(gg(0.2,4)) \\pm $(ee(0.2,4)) \$ &  \$c = $(cc(0.2,4)) \$ \\\\
        &  \$n=5 \$ &  \$$(gg(0.2,5)) \\pm $(ee(0.2,5)) \$ &  \$c = $(cc(0.2,5)) \$ \\\\
        &  \$n=6 \$ &  \$$(gg(0.2,6)) \\pm $(ee(0.2,6)) \$ &  \$c = $(cc(0.2,6)) \$ \\\\
        &  \$n=7 \$ &  \$$(gg(0.2,7)) \\pm $(ee(0.2,7)) \$ &  \$c = $(cc(0.2,7)) \$ \\\\
        &  \$n=8 \$ &  \${ $(gg(0.2,8))} \\pm $(ee(0.2,8)) \$ &  \$c = $(cc(0.2,8)) \$ \\\\
        &  \$n=9 \$ &  \${  $(gg(0.2,9))} \\pm $(ee(0.2,9)) \$ &  \$c = $(cc(0.2,9)) \$ \\\\
        &  \$n=10 \$ &  \${ $(gg(0.2,10))} \\pm $(ee(0.2,10)) \$ &  \$c = $(cc(0.2,10)) \$ \\\\
        &  \$n=11 \$ &  \${ $(gg(0.2,11))} \\pm $(ee(0.2,11)) \$ &  \$c = $(cc(0.2,11)) \$ \\\\
        &  \$n=12 \$ &  \${ $(gg(0.2,12))} \\pm $(ee(0.2,12)) \$ &  \$c = $(cc(0.2,12)) \$ \\\\
        &  \$n=13 \$ &  \${ $(gg(0.2,13))} \\pm $(ee(0.2,13)) \$ &  \$c = $(cc(0.2,13)) \$ \\\\
        &   \${\\bf n=14} \$ &  \${\\bf $(gg(0.2,14)) \\pm $(ee(0.2,14))} \$ &  \$c = $(cc(0.2,14)) \$ \\\\
        &  \$ {\\bf n=15} \$ &  \${\\bf $(gg(0.2,15)) \\pm $(ee(0.2,15))} \$ &  \$c = $(cc(0.2,15)) \$ \\\\
        &  \$ {\\bf n=16} \$ &  \${\\bf $(gg(0.2,16)) \\pm $(ee(0.2,16))} \$ &  \$c = $(cc(0.2,16)) \$ \\\\
        &  \$n=17 \$ &  \${ $(gg(0.2,17))} \\pm $(ee(0.2,17)) \$ &  \$c = $(cc(0.2,17)) \$ \\\\
        &  \$n=18 \$ &  \${ $(gg(0.2,18))} \\pm $(ee(0.2,18)) \$ &  \$c = $(cc(0.2,18)) \$ \\\\
        &  \$n=19 \$ &  \${ $(gg(0.2,19))} \\pm $(ee(0.2,19)) \$ &  \$c = $(cc(0.2,19)) \$ \\\\
        &  \$n=20 \$ &  \${ $(gg(0.2,20))} \\pm $(ee(0.2,20)) \$ &  \$c = $(cc(0.2,20)) \$ \\\\
        \\hline
        \\multirow{13}{*}{ \$\\gamma = 0.3\$}
        &  \$n=3 \$ &  \$$(gg(0.3,3)) \\pm $(ee(0.3,3)) \$ &  \$c = $(cc(0.3,3)) \$\\\\
        &  \$n=4 \$ &  \$$(gg(0.3,4)) \\pm $(ee(0.3,4)) \$ &  \$c = $(cc(0.3,4)) \$ \\\\
        &  \$n=5 \$ &  \${ $(gg(0.3,5))} \\pm $(ee(0.3,5)) \$ &  \$c = $(cc(0.3,5)) \$ \\\\
        &  \$n=6 \$ &  \$$(gg(0.3,6)) \\pm $(ee(0.3,6)) \$ &  \$c = $(cc(0.3,6)) \$ \\\\
        &  \$n=7 \$ &  \$$(gg(0.3,7)) \\pm $(ee(0.3,7)) \$ &  \$c = $(cc(0.3,7)) \$ \\\\
        &  \$n=8 \$ &  \${ $(gg(0.3,8))} \\pm $(ee(0.3,8)) \$ &  \$c = $(cc(0.3,8)) \$ \\\\
        &  \${\\bf n=9} \$ &  \${\\bf $(gg(0.3,9)) \\pm $(ee(0.3,9))} \$ &  \$c = $(cc(0.3,9)) \$ \\\\
        &  \$n=10 \$ &  \$$(gg(0.3,10)) \\pm $(ee(0.3,10)) \$ &  \$c = $(cc(0.3,10)) \$ \\\\
        &  \$n=11 \$ &  \$$(gg(0.3,11)) \\pm $(ee(0.3,11)) \$ &  \$c = $(cc(0.3,11)) \$ \\\\
        &  \$n=12 \$ &  \$$(gg(0.3,12)) \\pm $(ee(0.3,12)) \$ &  \$c = $(cc(0.3,12)) \$ \\\\
        &  \$n=13 \$ &  \$$(gg(0.3,13)) \\pm $(ee(0.3,13)) \$ &  \$c = $(cc(0.3,13)) \$ \\\\
        &  \$n=14 \$ &  \$$(gg(0.3,14)) \\pm $(ee(0.3,14)) \$ &  \$c = $(cc(0.3,14)) \$ \\\\
        &  \$n=15 \$ &  \$$(gg(0.3,15)) \\pm $(ee(0.3,15)) \$ &  \$c = $(cc(0.3,15)) \$ \\\\
        \\hline
        \\multirow{8}{*}{ \$\\gamma = 0.4\$}
        &  \$n=3 \$ &  \$$(gg(0.4,3)) \\pm $(ee(0.4,3)) \$ &  \$c = $(cc(0.4,3)) \$\\\\
        &  \$n=4 \$ &  \$$(gg(0.4,4)) \\pm $(ee(0.4,4)) \$ &  \$c = $(cc(0.4,4)) \$ \\\\
        &  \$n=5 \$ &  \${ $(gg(0.4,5))} \\pm $(ee(0.4,5)) \$ &  \$c = $(cc(0.4,5)) \$ \\\\
        &  \$n=6 \$ &  \${ $(gg(0.4,6))} \\pm $(ee(0.4,6)) \$ &  \$c = $(cc(0.4,6)) \$ \\\\
        &  \${\\bf n=7} \$ &  \${\\bf $(gg(0.4,7)) \\pm $(ee(0.4,7))} \$ &  \$c = $(cc(0.4,7)) \$ \\\\
        &  \$n=8 \$ &  \$$(gg(0.4,8)) \\pm $(ee(0.4,8)) \$ &  \$c = $(cc(0.4,8)) \$ \\\\
        &  \$n=9 \$ &  \$$(gg(0.4,9)) \\pm $(ee(0.4,9)) \$ &  \$c = $(cc(0.4,9)) \$ \\\\
        &  \$n=10 \$ &  \$$(gg(0.4,10)) \\pm $(ee(0.4,10)) \$ &  \$c = $(cc(0.4,10)) \$ \\\\

        \\hline
        \\multirow{3}{*}{ \$\\gamma = 0.6 \$}
        &  \$n=3 \$ &  \$$(gg(0.6,3)) \\pm $(ee(0.6,3)) \$ &  \$c = $(cc(0.6,3)) \$\\\\
        &  \${\\bf n=4} \$ &  \${ \\bf $(gg(0.6,4)) \\pm $(ee(0.6,4))} \$ &  \$c = $(cc(0.6,4)) \$ \\\\
        &  \$n=5 \$ &  \$$(gg(0.6,5)) \\pm $(ee(0.6,5)) \$ &  \$c = $(cc(0.6,5)) \$ \\\\
        \\hline
        \\multirow{3}{*}{ \$\\gamma = 0.8 \$}
        &  \${\\bf n=3 }\$ &  \${\\bf $(gg(0.8,3)) \\pm $(ee(0.8,3))} \$ &  \$c = $(cc(0.8,3)) \$\\\\
        &  \$n=4 \$ &  \$$(gg(0.8,4)) \\pm $(ee(0.8,4)) \$ &  \$c = $(cc(0.8,4)) \$ \\\\
        &  \$n=5 \$ &  \$ { $(gg(0.8,5))} \\pm $(ee(0.8,5)) \$ &  \$c = $(cc(0.8,5)) \$ \\\\
        \\hline
      \\end{tabular}
    \\end{center}
  \\end{table}"""
  println(tableString)
  nothing
end

plotDict, resultsDict = statsSummary(doPlot=true)
formatTable(resultsDict)
doPlot()


text:
# Mon 31 Dec 2012 01:39:13 PM EST
#
# NeHe Tut 12 - Rotate a textured cube
#
# Q - quit
# Up/Down - increase/decrease rotation of each cube about it's own x-axis
# Left/Right - increase/decrease rotation of each cube about it's own y-axis

# load necessary GLUT/OpenGL routines

global OpenGLver="1.0"
using OpenGL
using GLUT

# initialize variables

global window

global xrot = 0
global yrot = 0

global box  = 0
global top  = 0

boxcol      = [1.0 0.0 0.0;
               1.0 0.5 0.0;
               1.0 1.0 0.0;
               0.0 1.0 0.0;
               0.0 1.0 1.0]

topcol      = [0.5 0.0  0.0;
               0.5 0.25 0.0;
               0.5 0.5  0.0;
               0.0 0.5  0.0;
               0.0 0.5  0.5]

global tex  = Array(Uint32,1) # generating 1 texture

width       = 640
height      = 480

# load textures from images

function LoadGLTextures()
    global tex

    img, w, h = glimread(expanduser("~/.julia/GLUT/Examples/NeHe/tut12/cube.bmp"))

    glGenTextures(1,tex)
    glBindTexture(GL_TEXTURE_2D,tex[1])
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR_MIPMAP_NEAREST)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR)
    glTexImage2D(GL_TEXTURE_2D, 0, 3, w, h, 0, GL_RGB, GL_UNSIGNED_BYTE, img)
    gluBuild2DMipmaps(GL_TEXTURE_2D, 3, w, h, GL_RGB, GL_UNSIGNED_BYTE, img)
end

# function to init OpenGL context

function initGL(w::Integer,h::Integer)
    global box
    global top

    glViewport(0,0,w,h)
    LoadGLTextures()
    glClearColor(0.0, 0.0, 0.0, 0.5)
    glClearDepth(1.0)
    glDepthFunc(GL_LEQUAL)
    glEnable(GL_DEPTH_TEST)
    glShadeModel(GL_SMOOTH)
    glHint(GL_PERSPECTIVE_CORRECTION_HINT, GL_NICEST)

    # enable texture mapping
    glEnable(GL_TEXTURE_2D)

    #enable simple lighting
    glEnable(GL_LIGHT0)
    glEnable(GL_LIGHTING)
    glEnable(GL_COLOR_MATERIAL)

    glMatrixMode(GL_PROJECTION)
    glLoadIdentity()

    gluPerspective(45.0,w/h,0.1,100.0)

    glMatrixMode(GL_MODELVIEW)

    # build the display lists

    box = glGenLists(2)

    glNewList(box, GL_COMPILE)
        glBegin(GL_QUADS)
            # Bottom Face
            glTexCoord(0.0, 1.0)
            glVertex(-1.0, -1.0, -1.0)
            glTexCoord(1.0, 1.0)
            glVertex(1.0, -1.0, -1.0)
            glTexCoord(1.0, 0.0)
            glVertex(1.0, -1.0,  1.0)
            glTexCoord(0.0, 0.0)
            glVertex(-1.0, -1.0,  1.0)

            # Front Face
            glTexCoord(1.0, 0.0)
            glVertex(-1.0, -1.0,  1.0)
            glTexCoord(0.0, 0.0)
            glVertex(1.0, -1.0,  1.0)
            glTexCoord(0.0, 1.0)
            glVertex(1.0,  1.0,  1.0)
            glTexCoord(1.0, 1.0)
            glVertex(-1.0,  1.0,  1.0)

            # Back Face
            glTexCoord(0.0, 0.0)
            glVertex(-1.0, -1.0, -1.0)
            glTexCoord(0.0, 1.0)
            glVertex(-1.0,  1.0, -1.0)
            glTexCoord(1.0, 1.0)
            glVertex(1.0,  1.0, -1.0)
            glTexCoord(1.0, 0.0)
            glVertex(1.0, -1.0, -1.0)

            # Right Face
            glTexCoord(0.0, 0.0)
            glVertex(1.0, -1.0, -1.0)
            glTexCoord(0.0, 1.0)
            glVertex(1.0,  1.0, -1.0)
            glTexCoord(1.0, 1.0)
            glVertex(1.0,  1.0,  1.0)
            glTexCoord(1.0, 0.0)
            glVertex(1.0, -1.0,  1.0)

            # Left Face
            glTexCoord(1.0, 0.0)
            glVertex(-1.0, -1.0, -1.0)
            glTexCoord(0.0, 0.0)
            glVertex(-1.0, -1.0,  1.0)
            glTexCoord(0.0, 1.0)
            glVertex(-1.0,  1.0,  1.0)
            glTexCoord(1.0, 1.0)
            glVertex(-1.0,  1.0, -1.0)
        glEnd()
    glEndList()

    top = uint32(box+1)

    glNewList(top, GL_COMPILE)
        glBegin(GL_QUADS)
            # Top Face
            glTexCoord(1.0, 1.0)
            glVertex(-1.0, 1.0, -1.0)
            glTexCoord(1.0, 0.0)
            glVertex(-1.0, 1.0,  1.0)
            glTexCoord(0.0, 0.0)
            glVertex(1.0, 1.0,  1.0)
            glTexCoord(0.0, 1.0)
            glVertex(1.0, 1.0, -1.0)
        glEnd()
    glEndList()
end

# prepare Julia equivalents of C callbacks that are typically used in GLUT code

function ReSizeGLScene(w::Int32,h::Int32)
    if h == 0
        h = 1
    end

    glViewport(0,0,w,h)

    glMatrixMode(GL_PROJECTION)
    glLoadIdentity()

    gluPerspective(45.0,w/h,0.1,100.0)

    glMatrixMode(GL_MODELVIEW)

    return nothing
end

_ReSizeGLScene = cfunction(ReSizeGLScene, Void, (Int32, Int32))

function DrawGLScene()
    global tex
    global xrot
    global yrot
    global box
    global top

    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)

    glBindTexture(GL_TEXTURE_2D, tex[1])

    for yloop = 1:5
	      for xloop = 1:yloop
            glLoadIdentity()

            glTranslate(1.4+2.8xloop-1.4yloop, ((6.0-yloop)*2.4)-7.0, -20.0)

            glRotate(45.0-(2.0yloop)+xrot, 1.0, 0.0, 0.0)
            glRotate(45.0+yrot, 0.0, 1.0, 0.0)

            glColor(boxcol[yloop,:])
            glCallList(box)

            glColor(topcol[yloop,:])
            glCallList(top)
        end
    end

    glutSwapBuffers()

    return nothing
end

_DrawGLScene = cfunction(DrawGLScene, Void, ())

function keyPressed(the_key::Char,x::Int32,y::Int32)
    if the_key == int('q')
        glutDestroyWindow(window)
    end

    return nothing # keyPressed returns "void" in C. this is a workaround for Julia's "automatically return the value of the last expression in a function" behavior.
end

_keyPressed = cfunction(keyPressed, Void, (Char, Int32, Int32))

function specialKeyPressed(the_key::Int32,x::Int32,y::Int32)
    global xrot
    global yrot

    if the_key == GLUT_KEY_UP
        xrot -=0.2
    elseif the_key == GLUT_KEY_DOWN
        xrot +=0.2
    elseif the_key == GLUT_KEY_LEFT
        yrot -=0.2
    elseif the_key == GLUT_KEY_RIGHT
        yrot +=0.2
    end

    return nothing # specialKeyPressed returns "void" in C-GLUT. this is a workaround for Julia's "automatically return the value of the last expression in a function" behavior.
end

_specialKeyPressed = cfunction(specialKeyPressed, Void, (Int32, Int32, Int32))

# run GLUT routines

glutInit()
glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_ALPHA | GLUT_DEPTH)
glutInitWindowSize(width, height)
glutInitWindowPosition(0, 0)

window = glutCreateWindow("NeHe Tut 12")

glutDisplayFunc(_DrawGLScene)
glutFullScreen()

glutIdleFunc(_DrawGLScene)
glutReshapeFunc(_ReSizeGLScene)
glutKeyboardFunc(_keyPressed)
glutSpecialFunc(_specialKeyPressed)

initGL(width, height)

glutMainLoop()


text:
mel = 1.6021766208000002e-22
A = 1e-10

a = 4.3*A
Ei = 4.196*mel
theta = 45.75
D = 5*mel
h = 0.1
xi = 0.38/A

m_He3 = 5.0082343984338715e-27

ki = sqrt(Ei/HASlib.h2m(m_He3))
grid = HASlib.HexGrid(a)
b1,b2 = HASlib.rotate_reciprocal(grid.l,[1,0])
channels = HASlib.get_channels(b1,b2,ki,theta*pi/180,5)
channels = HASlib.sort_and_cut(channels,5,0)

@test_nowarn cc(channels,a,h,D,xi,theta,m_He3,(-1*A,6*A))


text:
# Kmer Composition
# ================
#
# K-mer composition counter of biological sequences.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

"""
    DNAKmerCounts{k}(seq::Sequence[, step=1])
    RNAKmerCounts{k}(seq::Sequence[, step=1])

Count ocurrences of short (≤ 32) `k`-mers in a sequence.

This method uses a dense table to store counts, requiring O(4^k) memory, so is
not recommended for large k-mer sizes.

# Arguments
* `seq`: A nucleotide sequence.
* `step=1`: K-mers counted are separated by this many nucleotides.
"""
immutable KmerCounts{T,K}
    data::Vector{UInt32}

    function KmerCounts(seq::Sequence, step::Integer=1)
        data = zeros(UInt32, 4^K)
        @inbounds for (_, x) in each(Kmer{T,K}, seq, step)
            data[convert(UInt64, x) + 1] += 1
        end
        return new(data)
    end
end

typealias DNAKmerCounts{K} KmerCounts{DNANucleotide,K}
typealias RNAKmerCounts{K} KmerCounts{RNANucleotide,K}

function Base.getindex{T,K}(counts::KmerCounts{T,K}, x::Kmer{T,K})
    @inbounds c = counts.data[convert(UInt64, x) + 1]
    return c
end

function Base.show{T,K}(io::IO, counts::KmerCounts{T,K})
    println(io, (T == DNANucleotide ? "DNA" : "RNA"), "KmerCounts{", K, "}:")
    for x in UInt64(1):UInt64(4^K)
        s = string(Kmer{T,K}(x - 1))
        println(io, "  ", s, " => ", counts.data[x])
    end
end


text:
const CompactSARTSATrajectory =
    Union{CircularCompactSARTSATrajectory,VectorialCompactSARTSATrajectory}

function get_trace(b::CompactSARTSATrajectory, s::Symbol)
    if s == :state || s == :action
        select_last_dim(b[s], 1:(nframes(b[s]) > 1 ? nframes(b[s]) - 1 : nframes(b[s])))
    elseif s == :reward || s == :terminal
        b[s]
    elseif s == :next_state
        select_last_dim(b[:state], 2:nframes(b[:state]))
    elseif s == :next_action
        select_last_dim(b[:action], 2:nframes(b[:action]))
    else
        throw(ArgumentError("unknown trace name: $s"))
    end
end

Base.length(b::CompactSARTSATrajectory) = length(b[:terminal])
Base.isempty(b::CompactSARTSATrajectory) = all(isempty(b[s]) for s in RTSA)

function Base.getindex(b::CompactSARTSATrajectory, i::Int)
    (
        state = select_last_dim(b[:state], i),
        action = select_last_dim(b[:action], i),
        reward = select_last_dim(b[:reward], i),
        terminal = select_last_dim(b[:terminal], i),
        next_state = select_last_dim(b[:state], i + 1),
        next_action = select_last_dim(b[:action], i + 1),
    )
end

function Base.empty!(b::CompactSARTSATrajectory)
    for s in RTSA
        empty!(b[s])
    end
    b
end

function Base.push!(b::CompactSARTSATrajectory, kv::Pair{Symbol})
    k, v = kv
    if k == :state || k == :next_state
        push!(b[:state], v)
    elseif k == :action || k == :next_action
        push!(b[:action], v)
    elseif k == :reward || k == :terminal
        push!(b[k], v)
    else
        throw(ArgumentError("unknown trace name: $k"))
    end
    b
end

function Base.pop!(t::CompactSARTSATrajectory, s::Symbol)
    if s == :state || s == :next_state
        pop!(t[:state])
    elseif s == :action || s == :next_action
        pop!(t[:action])
    elseif s == :reward || s == :terminal
        pop!(t[s])
    else
        throw(ArgumentError("unknown trace name: $s"))
    end
    t
end

function Base.pop!(t::CompactSARTSATrajectory)
    if length(t) <= 0
        throw(ArgumentError("can not pop! from an empty trajectory"))
    else
        NamedTuple{RTSA}((
            pop!(t, :reward),
            pop!(t, :terminal),
            pop!(t, :state),
            pop!(t, :action),
        ))
    end
    t
end


text:
#---------------------------------------------------------
# # [SPECTrecon rotation](@id 2-rotate)
#---------------------------------------------------------

# This page explains the image rotation portion of the Julia package
# [`SPECTrecon.jl`](https://github.com/JeffFessler/SPECTrecon.jl).

# ### Setup

# Packages needed here.

using SPECTrecon
using MIRTjim: jim, prompt
using Plots: scatter, scatter!, plot!, default
default(markerstrokecolor=:auto, markersize=3)

# The following line is helpful when running this example.jl file as a script;
# this way it will prompt user to hit a key after each figure is displayed.

isinteractive() ? jim(:prompt, true) : prompt(:draw);

# ### Overview

# The first step in SPECT image forward projection
# is to rotate each slice of a 3D image volume
# to the appropriate view angle.
# In principle
# one could use any of numerous candidate interpolation methods
# for this task.
# However, because emission images are nonnegative
# and maximum likelihood methods
# for SPECT image reconstruction
# exploit that nonnegativity,
# it is desirable to use interpolators
# that preserve nonnegativity.
# This constraint rules out quadratic and higher B-splines,
# including the otherwise attractive cubic B-spline methods.
# On the other hand, nearest-neighbor interpolation
# (equivalent to 0th-order B-splines)
# does not provide adequate image quality.
# This leaves 1st-order interpolation methods
# as the most viable options.

# This package supports two 1st-order linear interpolators:
# * 2D bilinear interpolation,
# * a 3-pass rotation method based on 1D linear interpolation.

# Because image rotation is done repeatedly
# (for every slice of both the emission image and the attenuation map,
# for both projection and back-projection,
# and for multiple iterations),
# it is important for efficiency
# to use mutating methods
# rather than to repeatedly make heap allocations.

# Following other libraries like
# [FFTW.jl](https://github.com/JuliaMath/FFTW.jl),
# the rotation operations herein start with a `plan`
# where work arrays are preallocated
# for subsequent use.
# The `plan` is a `Vector` of `PlanRotate` objects:
# one for each thread.
# (Parallelism is across slices for a 3D image volume.)
# The number of threads defaults to `Threads.nthreads()`.

# ### Example

# Start with a 3D image volume (just 2 slices here for simplicity).

T = Float32 # work with single precision to save memory
image = zeros(T, 64, 64, 2)
image[30:50,20:30,1] .= 1
image[25:28,20:40,2] .= 1
jim(image, "Original image")

# Now plan the rotation
# by specifying
# * the image size `nx` (it must be square, so `ny=nx` implicitly)
# * the `DataType` used for the work arrays.

plan2 = plan_rotate(size(image, 1); T)

# Here are the internals for the plan for the first thread:

plan2[1]

# With this `plan` preallocated, now we can rotate the image volume,
# specifying the rotation angle in radians:

result2 = similar(image) # allocate memory for the result
imrotate!(result2, image, π/6, plan2) # mutates the first argument
jim(result2, "Rotated image by π/6 (2D bilinear)")

# The default, shown above, uses 2D bilinear interpolation for rotation.
# That default is the recommended approach because it is faster.

# Here is the 3-pass 1D interpolation approach,
# included mainly for checking consistency
# with the historical ASPIRE approach used in Matlab version of MIRT.

plan1 = plan_rotate(size(image, 1); T, method=:one)

# Here are the plan internals for the first thread:

plan1[1]

# The results of rotation using 3-pass 1D interpolation look quite similar:

result1 = similar(image)
imrotate!(result1, image, π/6, plan1)
jim(result1, "Rotated image by π/6 (3-pass 1D)")

# Here are the difference images for comparison.

jim(result1 - result2, "Difference images")

# ### Adjoint

# To ensure adjoint consistency between SPECT forward- and back-projection,
# there is also an adjoint routine:

adj2 = similar(result2)
imrotate_adj!(adj2, result2, π/6, plan2)
jim(adj2, "Adjoint image rotation (2D)")

adj1 = similar(result1)
imrotate_adj!(adj1, result1, π/6, plan1)
jim(adj1, "Adjoint image rotation (3-pass 1D)")

# The adjoint is *not* the same as the inverse
# so one does not expect the output here to match the original image!

# ### LinearMap

# One can form a linear map corresponding to image rotation using `LinearMapAA`.
# An operator like this may be useful
# as part of a motion-compensated image reconstruction method.

using LinearMapsAA: LinearMapAA

nx = 20 # small size for illustration
r1 = plan_rotate(nx; T, nthread = 1, method=:two)[1]
r2 = plan_rotate(nx; T, nthread = 1, method=:one)[1]
idim = (nx,nx)
odim = (nx,nx)
forw! = (y,x) -> imrotate!(y, x, π/6, r1)
back! = (x,y) -> imrotate_adj!(x, y, π/6, r1)
A1 = LinearMapAA(forw!, back!, (prod(odim),prod(idim)); T, odim, idim)
forw! = (y,x) -> imrotate!(y, x, π/6, r2)
back! = (x,y) -> imrotate_adj!(x, y, π/6, r2)
A2 = LinearMapAA(forw!, back!, (prod(odim),prod(idim)); T, odim, idim)

Afull1 = Matrix(A1)
Aadj1 = Matrix(A1')
Afull2 = Matrix(A2)
Aadj2 = Matrix(A2')
jim(cat(dims=3, Afull1', Aadj1', Afull2', Aadj2'), "Linear map for 2D rotation and its adjoint")

# The following verify adjoint consistency:
@assert Afull1' ≈ Aadj1
@assert Afull2' ≈ Aadj2

# Applying this linear map to a 2D or 3D image performs rotation:

image2 = zeros(nx,nx); image2[4:6, 5:13] .= 1
jim(cat(dims=3, image2, A2 * image2), "Rotation via linear map: 2D")

# Here is 3D too.
# The `A2 * image3` here uses the advanced "operator" feature of
# [LinearMapsAA.jl](https://github.com/JeffFessler/LinearMapsAA.jl).

image3 = cat(dims=3, image2, image2')
jim(cat(dims=4, image3, A2 * image3), "Rotation via linear map: 3D")

# Examine row and column sums of linear map

scatter(xlabel="pixel index", ylabel="row or col sum")
scatter!(vec(sum(Afull1, dims=1)), label="dim1 sum1", marker=:x)
scatter!(vec(sum(Afull1, dims=2)), label="dim2 sum1", marker=:square)
scatter!(vec(sum(Afull2, dims=1)), label="dim1 sum2")
scatter!(vec(sum(Afull2, dims=2)), label="dim2 sum2")


text:
@inline function DRA(Cell)
    """
    Function for Discrete Realisation Algorithm.

    DRA(Cell,s,f)

    """

    # Method Check
    if Cell.RA.Fs==(1/Cell.RA.SamplingT)
        println("CI-DRA Utilised")
    else
        println("DRA Utilised")
    end

    #Additional Pulse Setup
    tfft = (1/Cell.RA.Fs)*Cell.RA.f
    OrgT = Cell.RA.SamplingT*(0:floor(tfft[end]/Cell.RA.SamplingT))

    #Initialise Loop Variables
    A = Matrix{Float64}(I,Cell.RA.M+1,Cell.RA.M+1)
    B = Vector{Float64}(undef,Cell.RA.M+1)
    C = Matrix{Float64}(undef,Cell.RA.Outs,Cell.RA.M+1)
    D = Vector{Float64}(undef,Cell.RA.Outs)
    puls = Array{Float64}(undef,Cell.RA.Outs,size(OrgT,1)-1)
    C_Aug = Vector{Float64}(undef,Cell.RA.Outs)
    i=Int(1)
    l=Int(1)
    u=Int(0)

    for run in Cell.Transfer.tfs
        tf = Array{ComplexF64}(undef,size(Cell.Transfer.Locs[i],1),size(Cell.RA.s,2))
        Di = Vector{Float64}(undef,size(Cell.Transfer.Locs[i],1))::Vector{Float64}
        res0 = Vector{Float64}(undef,size(Cell.Transfer.Locs[i],1))::Vector{Float64}
        samplingtf = Array{Float64}(undef,size(Cell.Transfer.Locs[i],1),length(OrgT))

        if Cell.Transfer.Elec[i] == "Pos"
            run(Cell,Cell.RA.s,Cell.Transfer.Locs[i],"Pos",tf,Di,res0)
        elseif Cell.Transfer.Elec[i] == "Neg"
            run(Cell,Cell.RA.s,Cell.Transfer.Locs[i],"Neg",tf,Di,res0)
        else
            run(Cell,Cell.RA.s,Cell.Transfer.Locs[i],tf,Di,res0)
        end

        # if Cell.RA.Fs==(1/Cell.RA.SamplingT)
        #     f = one(Int)
        #     λ = Int(Cell.RA.Fs*Cell.RA.SamplingT)
        #     for i in 1:λ:size(tfft,2), j in 1:λ:size(tfft,1)
        #         Org[j,f] = tf[j,i]
        #         f+=1
        #     end
        # end

        u += Int(size(Cell.Transfer.Locs[i],1))::Int

        if Cell.RA.Fs==(1/Cell.RA.SamplingT)
            puls[l:u,:] .= real(ifft(tf,2))[:,2:end]
        else
            stpsum = cumsum(real(ifft(tf,2)), dims=2) #cumulative sum of tf response * sample time
            #Interpolate H(s) to obtain h_s(s) to obtain discrete-time impulse response
            for k in 1:size(stpsum,1)
                spl1 = CubicSplineInterpolation(tfft, stpsum[k,:]; bc=Line(OnGrid()), extrapolation_bc=Throw())
                samplingtf[k,:] = spl1(OrgT)
            end
            puls[l:u,:] .= diff(samplingtf, dims=2)
        end

        #puls[l:u,:] .= diff(stpsum, dims=2)
        #puls[l:u,:] .= real(ifft(tf,2))[:,2:end]
        #puls[l:u,:] .= diff(samplingtf, dims=2)
        D[l:u,:] .= Di
        C_Aug[l:u,:] .= res0
        l = u+one(u)::Int
        i += one(i)

    end

    #Scale Transfer Functions in Pulse Response
    SFactor = sqrt.(sum(puls.^2,dims=2))
    puls .= puls./SFactor

    #Pre-Allocation for Hankel & SVD
    Puls_L = size(puls,1)
    𝐇 = Array{Float64}(undef,length(Cell.RA.H1)*Puls_L,length(Cell.RA.H2))
    U,S,V = fh!(𝐇,Cell.RA.H1,Cell.RA.H2,puls,Cell.RA.M,Puls_L)

    # for i in 1:size(puls,1)
    #     if U[i,1] < 0
    #         U[i,:] .= -U[i,:]
    #         V[:,i] .= -V[:,i]
    #     end
    # end

    # Create Observibility and Control Matrices -> Create A, B, and C
    ʃ = sqrt(Diagonal(S))
    #Observibility = U*ʃ
    #Control = ʃ*V'

    U .= U*ʃ
    V .= ʃ*V

    #A[2:end,2:end] .= (Observibility\𝐇)/Control
    A[2:end,2:end] .= (U\𝐇)/V

    #Performance check
    if any(i -> i>1., real(eigvals(A)))
        println("Oscilating System: A has indices of values > 1")
    end

    if any(i -> i<0., real(eigvals(A)))
        println("Unstable System: A has indices of negative values")
    end

    #B .= [Cell.RA.SamplingT; Control[:,Cell.RA.N]]
    #C .= [C_Aug SFactor.*Observibility[1:size(puls,1),:]]
    B .= [Cell.RA.SamplingT; V[:,1:Cell.RA.N]]
    C .= [C_Aug SFactor.*U[1:Puls_L,:]]

    #  Final State-Space Form
    # d, Sᵘ = eigen(A,sortby=nothing)
    # d = mag!(d) # taking the magnitude of S and maintaining the sign of the real values
    # S = mag!(S)
    # A = Diagonal(d)
    # C = C*Sᵘ.*(inv(Sᵘ)*B)'
    # B = ones(size(B))

return mag!(A), mag!(B), mag!(C), D

end

max_stars_repo_name: BradyPlanden/LiBRA

max_forks_repo_name: BradyPlanden/LiBRA


text:
#cd("Dropbox/PhD/JonAustin/Model")
using DelimitedFiles

# Include my functions.....
include("func_Gnp.jl")

#p = 0.1 # Probability of forming simplex
N = 100 # Number of nodes
m = 3 # simplex size
nvar_max = 100

filename=string("Outputs/jl_3simplex","_N",N,"_m",m,".txt")

# Fill simplices with probability p
# -------------

# Keep track of simplices

for nvar=1:nvar_max

    p = 1/(nvar*nvar)
    println("doing run\t",nvar,"\tout of \t",nvar_max)

    Triples = Dict()
    Pairs = Dict()

    for a=0:N-m
        for b=a+1:N-(m-1)
            for c=b+1:N-(m-2)
                simplex = [a,b,c]
                random_float = rand(Float64)
                if (random_float<p)
                    Triples[simplex]=1
                    Pairs = fill_pairs(Pairs,simplex)
                end
            end
        end
    end

    # Evaluate prevalence of patterns
    # -------------
    results = Dict()
    for a=0:N-m
        for b=a+1:N-(m-1)
            for c=b+1:N-(m-2)
                simplex = [a,b,c]
                results = find_pattern(Triples,Pairs,results,simplex)

            end
        end
    end

    results = normalize_results(results)
    results["p"] = p
    write_results(filename,results)
    #println(results)
end

max_stars_repo_name: jonassjuul/team-formation

max_issues_repo_name: jonassjuul/team-formation

max_forks_repo_name: jonassjuul/team-formation


text:
estimate_move_time(fex::nGraph.FluxExecutable, frame::Frame{ILPHolder{Static}}) = zero(Float64)
estimate_move_time(fex::nGraph.FluxExecutable, frame::Frame) =
    estimate_move_time(fex.ex.ngraph_function, frame)

function estimate_move_time(f::nGraph.NFunction, frame::Frame)
    move_time = zero(Float64)
    for _node in f
        node = NodeDescriptor(_node)

        # If this is a move, determine which direction data is being moved and add the move
        # time estimate to the rolling counter.
        if description(node) == "Move"
            tensor = first(outputs(node))
            if nGraph.is_persistent(tensor)
                move_time += sizeof(tensor) / wb(frame.modeltype)
            else
                move_time += sizeof(tensor) / rb(frame.modeltype)
            end
        end
    end
    return move_time
end

# your moves are weak
function profile_moves(fex)
    timing_data = read_timing_data(fex.ex.ngraph_function)
    computed_stats = Dict{String, NamedTuple}()
    for node_unwrapped in fex.ex.ngraph_function
        node = NodeDescriptor(node_unwrapped)
        ismove(node) || continue

        time = timing_data[findfirst(x -> x["name"] == nGraph.name(node), timing_data)]["dur"]
        # Convert bytes to GB, time from μs to s
        bytes = sizeof(first(inputs(node)))
        bandwidth = (bytes / 1E9) / (time / 1E6)
        computed_stats[nGraph.name(node)] = (
            bytes = bytes,
            bandwidth = bandwidth,
            write_to_pmem = !nGraph.is_persistent(first(inputs(node))),
        )
    end

    # Summarize read and write bandwidth
    println("Read Bandwidths")
    for f in (ismoveasync, !ismoveasync)
        s = 0
        count = 0
        for (node_name, stats) in computed_stats
            f(node_name) || continue
            if stats.write_to_pmem == false
                println("$node_name => $(stats.bandwidth) GB/s")
                println("    size: $(stats.bytes) B")

                if !isinf(stats.bandwidth)
                    s += stats.bandwidth
                    count += 1
                end
            end
        end

        println()
        println("Average Bandwidth: ", s / count)
        println()
    end
    println()
    println("Write Bandwidths")
    for f in (ismoveasync, !ismoveasync)
        s = 0
        count = 0
        for (node_name, stats) in computed_stats
            f(node_name) || continue
            if stats.write_to_pmem == true
                println("$node_name => $(stats.bandwidth) GB/s")
                println("    size: $(stats.bytes) B")

                if !isinf(stats.bandwidth)
                    s += stats.bandwidth
                    count += 1
                end
            end
        end
        println()
        println("Average Bandwidth: ", s / count)
        println()
    end

    return computed_stats
end

# Inspection routines for the post optimized ILP model
function list_overlaps(frame::Frame)
    model = frame.model
    node_times = model[:node_times]

    total_node_execution_time = 0.0
    total_move_time = 0.0
    for node in nodes(frame.profile_data)
        node_name = nGraph.name(node)
        haskey(node_times, node_name) || continue
        _async = get(model[:tensor_async], node_name, nothing)
        if !isnothing(_async) && !iszero(JuMP.value(_async))
            # Get the values of the asynchronous move times
            async_move_time = JuMP.value(_async)
            node_execution_time = JuMP.value(node_times[node_name])

            total_move_time += async_move_time

            @info """
            Overlap times for $node_name
            Node Execution Time: $node_execution_time
            Async Move Time : $async_move_time
            """
        end
        total_node_execution_time += JuMP.value(node_times[node_name])
    end

    println()
    @info """
    Total Node Execution Time: $total_node_execution_time
    Total Move Time: $total_move_time
    """
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
abstract type CollisionBox2 end
collide(a::CollisionBox2,b::CollisionBox2)=collide(b,a)
struct EmptyCollisionBox<:CollisionBox2 end
collide(::CollisionBox2,::EmptyCollisionBox)=false
struct BCollisionBox{T}<:CollisionBox2 where T<:Number # 横平竖直的矩形
	lx::T
	ly::T
	rx::T
	ry::T
end
struct SegCollisionBox{T}<:CollisionBox2 where T<:Number # 线段,保证lx<=rx
	lx::T
	ly::T
	rx::T
	ry::T
end
struct CirCollisionBox{T}<:CollisionBox2 where T<:Number # 圆
	x::T
	y::T
	r::T
end
function collide(s::SegCollisionBox,b::BCollisionBox)
	lxin=(b.lx<s.lx<b.rx)
	rxin=(b.lx<s.rx<b.rx)
	lyin=(b.ly<s.ly<b.ry)
	ryin=(b.ly<s.ry<b.ry)
	if  lxin&&rxin&&lyin&&ryin # 线段在矩形内
		return true
	end
	if s.ly==s.ry # 线段水平
		if lyin && (lxin||rxin)
			return true
		end
		return false
	end
	if s.lx==s.rx # 线段竖直
		if lxin && (lyin||ryin)
			return true
		end
		return false
	end
	# s: y=ax+k
	a=(s.ry-s.ly)/(s.rx-s.lx)
	k=(s.rx*s.ly-s.lx*s.ry)/(s.rx-s.lx)
	if (s.lx<b.lx<s.rx&&b.ly<a*b.lx+k<b.ry) || # 交x=b.lx
		(s.lx<b.rx<s.rx&&b.ly<a*b.rx+k<b.ry) # 交x=b.rx
		return true
	end
	if s.ly<s.ry
		sl=s.ly
		sr=s.ry
	else
		sl=s.ry
		sr=s.ly
	end
	if (sl<b.ly<sr&&b.lx<(b.ly-k)/a<b.rx) || # 交y=b.ly
		(sl<b.ry<sr&&b.lx<(b.ry-k)/a<b.rx) # 交y=b.ry
		return true
	end
	return false
end
function collide_with(ci::CirCollisionBox,a::Number,b::Number,c::Number) # ax+by+c=0
	return abs(a*ci.x+b*ci.y+c)<hypot(a,b)*ci.r
end
function collide(c::CirCollisionBox,b::BCollisionBox)
	# 圆角矩形，分成2个矩形4个圆
	return (b.lx<c.x<b.ly&&b.ly-c.r<c.y<b.ry+c.r) ||
	(b.ly<c.y<b.ry&&b.lx-c.r<c.x<b.rx+c.r) ||
	(hypot(c.x-b.lx,c.y-b.ly)<c.r) ||
	(hypot(c.x-b.lx,c.y-b.ry)<c.r) ||
	(hypot(c.x-b.rx,c.y-b.ly)<c.r) ||
	(hypot(c.x-b.rx,c.y-b.ry)<c.r)
end
function collide(c::CirCollisionBox,bs::Dimension)
	lx=floor(c.x-c.r);rx=floor(c.x+c.r)
	ly=floor(c.y-c.r);ry=floor(c.y+c.r)
	for i in lx:ly
		for j in rx:ry
			if collide(c,colbox(blk(bs,i,j),i,j))
				return true
			end
		end
	end
	return false
end
function collide_move(c::CirCollisionBox,bs::Dimension,vx::Number,vy::Number,max_depth::UInt8=0x8)
	lx=0;ly=0
	rx=vx;ry=vy
	anx=0;any=0
	while true # 二分，原则：保证不碰撞
		mx=(lx+rx)/2
		my=(ly+ry)/2
		t=CirCollisionBox(c.x+mx,c.y+my,c.r)
		if collide(t,bs)
			if max_depth==0x0
				return Pair(anx,any)
			end
			rx=mx
			ry=my
		else
			if max_depth==0x0
				return Pair(mx,my)
			end
			anx=lx=mx
			any=ly=my
		end
		max_depth-=0x1
	end
end


text:
module codejam_io

export parse_int, parse_bigint, rep
export read_straight, read_simple, read_simple_2, read_N, read_f, read_f_multi, read_fixed
export write_simple, write_slst, write_slst_h
export getFinFout, getFinFout_A, getFinFout_B, getFinFout_C, getFinFout_D, getFinFout_E, getFinFout_F
export nestedToMatrix

function nestedToMatrix{T<:Any}(x::Array{T,1})
    C, maxC = extrema(length.(x))
    if C != maxC
        error("Data is not actually a matrix")
    end
    R = length(x)
    return [x[i][j] for i in 1:R, j in 1:C]
end

function parse_int(x)
	parse(Int, x)
end

function parse_bigint(x)
	parse(BigInt, x)
end

function rep(x, n)
	[x for i in 1:n]
end

# Read a file straight into a list of strings
function read_straight(fname::AbstractString)
	open(fname, "r") do f
        strip.(readlines(f))
    end
end

# Quick read if the data is one-per-line integers
function read_simple(fname::AbstractString, fMap=parse_int)
	open(fname, "r") do f
		[fMap.(split(strip(x))) for x in readlines(f)[2:end]]
	end
end

# One-per-twoline with first line removable
function read_simple_2(fname::AbstractString, fMap=parse_int)
	open(fname, "r") do f
		[fMap.(split(strip(x))) for x in readlines(f)[3:2:end]]
	end
end

# Quick read T groups of N lines each (N given a line before as first element)
function read_N(fname::AbstractString; N_pos = 1, int_parser=parse_int, fMap=string, fOfN = (x->x), storeHead = false)
	L = open(fname, "r") do f
		split.(strip.(readlines(f)))
	end
	T = int_parser(L[1][1])
	output = []
	ln = 2
	for i in 1:T
		N = fOfN(int_parser(L[ln][N_pos]))
		if storeHead
			push!(output,[[fMap.(L[ln])], [fmap.(x) for x in L[ln+1:ln+N]]])
		else
			push!(output, [fMap.(x) for x in L[ln+1:ln+N]])
		end
		ln += (N+1)
	end
	return output
end

# Quick read T groups of lines given by some function of the first line among them
# Exact same default behaviour as read_N; but more flexible
function read_f(fname::AbstractString, fMapFn=x->rep(string, x[1]); int_parser=parse_int, storeHead = false)
	L = open(fname, "r") do f
		split.(strip.(readlines(f)))
	end
	T = int_parser(L[1][1])
	output = []
	ln = 2
	for i in 1:T
		headerRow = int_parser.(L[ln])
		fMap = fMapFn(headerRow)
		N = length(fMap)
		if storeHead
			push!(output, [headerRow, [fMap[k].(x) for (k,x) in enumerate(L[ln+1:ln+N])]])
		else
			push!(output, [fMap[k].(x) for (k,x) in enumerate(L[ln+1:ln+N])])
		end
		ln += (N+1)
	end
	return output
end
# e.g. L = read_f('C-test.in', fMapFn=x->vcat(rep(string, x[1]), rep(parse_int,x[1]), rep(parse_int, x[2])))

# READ_F with multi-line headers (also saves the headers if storeHead=true)
# This is the most general read function and can mimic (almost) the behaviour of all the others
function read_f_multi(fname::AbstractString,fMapFn=x->rep(string,x[1][1]); nlines_head=1,int_parser=parse_int, storeHead=true)
	L = open(fname, "r") do f
		split.(strip.(readlines(f)))
	end
	T = int_parser(L[1][1])
	output = []
	ln = 2
	for i in 1:T
		head = [int_parser.(x) for x in L[ln:ln+nlines_head-1]]
		fMap = fMapFn(head)
		N = length(fMap)
		ln += nlines_head
		if storeHead
			push!(output, [head,[fMap[k].(x) for (k,x) in enumerate(L[ln:ln+N-1])]])
		else
			push!(output, [fMap[k].(x) for (k,x) in enumerate(L[ln:ln+N-1])])
		end
		ln += N
	end
	return output
end

# Read a sequence of inputs each one of which is a fixed number of lines
function read_fixed(fname::AbstractString, fMaps; int_parser=parse_int)
	L = open(fname, "r") do f
		split.(strip.(readlines(f)))
	end
	T = int_parser(L[1][1])
	ln = 2
	N = length(fMaps)
	output = []
	for i in 1:T
		push!(output, [fMaps[k].(L[ln+k-1]) for k in 1:N])
		ln+=N
	end
	return output
end

# Quick one-item-per-line outputs
function write_simple(fname, lst)
	f = open(fname, "w")
	for (i,elt) in enumerate(lst)
		println(f, "Case #$i: $(string(elt))")
	end
	close(f)
end

# Write a list of lists of strings
function write_slst(fname, strlst)
	f = open(fname, "w")
	for (i, L) in enumerate(strlst)
		println(f,"Case #$i:")
		for s in L
			println(f,join(string.(s), " "))
		end
	end
	close(f)
end

# Write a list of lists of strings, but include a word/number on the first line with 'CASE'
function write_slst_h(fname, strlst)
	f = open(fname, "w")
	for (i, L) in enumerate(strlst)
		println(f,"Case #$i: $(L[1])")
		for s in L[2:end]
			println(f,join(string.(s), " "))
		end
	end
	close(f)
end

function getFinFout(v, default_in)
    if length(v) == 0
        fin = default_in
    else
        fin = v[1]
    end
    if length(v) < 2
        fpos = search(fin, '.')
        if fpos == 0
            fin = default_in
        end
        fout = fin[1:fpos]*"out"
    else
        fout = v[2]
    end

    return fin, fout
end

function getFinFout_A(v, default_in = "A-sample.in")
	return getFinFout(v, default_in)
end

function getFinFout_B(v, default_in = "B-sample.in")
	return getFinFout(v, default_in)
end

function getFinFout_C(v, default_in = "C-sample.in")
	return getFinFout(v, default_in)
end

function getFinFout_D(v, default_in = "D-sample.in")
	return getFinFout(v, default_in)
end

function getFinFout_E(v, default_in = "E-sample.in")
	return getFinFout(v, default_in)
end

function getFinFout_F(v, default_in = "F-sample.in")
	return getFinFout(v, default_in)
end

end

max_stars_repo_name: DarioFanucchi/CompetitionCode

max_issues_repo_name: DarioFanucchi/CompetitionCode

max_forks_repo_name: DarioFanucchi/CompetitionCode


text:
using Test, LinearMaps, LinearAlgebra

@testset "transpose/adjoint" begin
    A = 2 * rand(ComplexF64, (20, 10)) .- 1
    v = rand(ComplexF64, 10)
    w = rand(ComplexF64, 20)
    V = rand(ComplexF64, 10, 3)
    W = rand(ComplexF64, 20, 3)
    Av = A * v
    AV = A * V
    M = @inferred LinearMap(A)
    N = @inferred LinearMap(M)

    @test @inferred M' * w == A' * w
    @test @inferred mul!(copy(V), adjoint(M), W) ≈ A' * W
    @test @inferred transpose(M) * w == transpose(A) * w
    @test @inferred transpose(LinearMap(transpose(M))) * v == A * v
    @test @inferred LinearMap(M')' * v == A * v
    @test @inferred transpose(transpose(M)) == M
    @test (M')' == M
    Mherm = @inferred LinearMap(A'A)
    @test @inferred ishermitian(Mherm)
    @test @inferred !issymmetric(Mherm)
    @test @inferred !issymmetric(transpose(Mherm))
    @test @inferred ishermitian(transpose(Mherm))
    @test @inferred ishermitian(Mherm')
    @test @inferred isposdef(Mherm)
    @test @inferred isposdef(transpose(Mherm))
    @test @inferred isposdef(adjoint(Mherm))
    @test @inferred !(transpose(M) == adjoint(M))
    @test @inferred !(adjoint(M) == transpose(M))
    @test @inferred transpose(M') * v ≈ transpose(A') * v
    @test @inferred transpose(LinearMap(M')) * v ≈ transpose(A') * v
    @test @inferred LinearMap(transpose(M))' * v ≈ transpose(A)' * v
    @test @inferred transpose(LinearMap(transpose(M))) * v ≈ Av
    @test @inferred adjoint(LinearMap(adjoint(M))) * v ≈ Av

    @test @inferred mul!(copy(w), transpose(LinearMap(M')), v) ≈ transpose(A') * v
    @test @inferred mul!(copy(w), LinearMap(transpose(M))', v) ≈ transpose(A)' * v
    @test @inferred mul!(copy(w), transpose(LinearMap(transpose(M))), v) ≈ Av
    @test @inferred mul!(copy(w), adjoint(LinearMap(adjoint(M))), v) ≈ Av
    @test @inferred mul!(copy(V), transpose(M), W) ≈ transpose(A) * W
    @test @inferred mul!(copy(V), adjoint(M), W) ≈ A' * W

    B = @inferred LinearMap(Symmetric(rand(10, 10)))
    @test transpose(B) == B
    @test B == transpose(B)

    B = @inferred LinearMap(Hermitian(rand(ComplexF64, 10, 10)))
    @test adjoint(B) == B
    @test B == B'

    CS = @inferred LinearMap{ComplexF64}(cumsum, x -> reverse(cumsum(reverse(x))), 10; ismutating=false)
    for transform in (adjoint, transpose)
        @test transform(transform(CS)) == CS
    end
    @test transpose(CS) != adjoint(CS)
    @test adjoint(CS) != transpose(CS)
    M = Matrix(CS)
    x = rand(ComplexF64, 10)
    for transform1 in (adjoint, transpose), transform2 in (adjoint, transpose)
        @test transform2(LinearMap(transform1(CS))) * x ≈ transform2(transform1(M))*x
    end

    id = @inferred LinearMap(identity, identity, 10; issymmetric=true, ishermitian=true, isposdef=true)
    for transform in (adjoint, transpose)
        @test transform(id) == id
        for prop in (issymmetric, ishermitian, isposdef)
            @test prop(transform(id))
        end
    end
end


text:
import StructArrays
using StructArrays: StructVector
"""
chainvec(x::T, n=1)

Return a vector used to store the type T in a SampleChain. Satisfies the law
chainvec(x,n)[1] == x

If `n>1`, the remaining entries will be `undef`. Note that `x` should correspond
to information stored in a single sample.
"""
function chainvec end

function chainvec(nt::NamedTuple, n = 1)
    tv = TupleVector(undef, nt, n)
    @inbounds tv[1] = nt
    return tv
end

function chainvec(x::T, n = 1) where {T<:Real}
    ev = ElasticVector{T}(undef, n)
    @inbounds ev[1] = x
    return ev
end

function chainvec(x::T, n = 1) where {T}
    if isstructtype(T)
        v = StructArrays.replace_storage(ElasticVector, StructVector{T}(undef, n))
        @inbounds v[1] = x
        return v
    end
    @inbounds ev[1] = x
    return ev
end

function chainvec(x::DenseArray{T,N}, n = 1) where {T,N}
    nv = nestedview(ElasticArray{T,N + 1}(undef, size(x)..., n), N)
    @inbounds nv[1] = x
    return nv
end


text:
using Flux, Flux.Tracker, CuArrays, Test
using Flux: gpu

@info "Testing GPU Support"

@testset "CuArrays" begin

CuArrays.allowscalar(false)

x = param(randn(5, 5))
cx = gpu(x)
@test cx isa TrackedArray && cx.data isa CuArray

@test Flux.onecold(param(gpu([1.,2.,3.]))) == 3

x = Flux.onehotbatch([1, 2, 3], 1:3)
cx = gpu(x)
@test cx isa Flux.OneHotMatrix && cx.data isa CuArray
@test (cx .+ 1) isa CuArray

m = Chain(Dense(10, 5, tanh), Dense(5, 2), softmax)
cm = gpu(m)

@test all(p isa TrackedArray && p.data isa CuArray for p in params(cm))
@test cm(gpu(rand(10, 10))) isa TrackedArray{Float32,2,CuArray{Float32,2}}

x = [1,2,3]
cx = gpu(x)
@test Flux.crossentropy(x,x) ≈ Flux.crossentropy(cx,cx)

xs = param(rand(5,5))
ys = Flux.onehotbatch(1:5,1:5)
@test collect(cu(xs) .+ cu(ys)) ≈ collect(xs .+ ys)

c = gpu(Conv((2,2),3=>4))
l = c(gpu(rand(10,10,3,2)))
Flux.back!(sum(l))

c = gpu(CrossCor((2,2),3=>4))
l = c(gpu(rand(10,10,3,2)))
Flux.back!(sum(l))

end

@testset "onecold gpu" begin
  y = Flux.onehotbatch(ones(3), 1:10) |> gpu;
  @test Flux.onecold(y) isa CuArray
  @test y[3,:] isa CuArray
end

if CuArrays.libcudnn != nothing
    @info "Testing Flux/CUDNN"
    include("cudnn.jl")
    if !haskey(ENV, "CI_DISABLE_CURNN_TEST")
      include("curnn.jl")
    end
end


text:
using WAVI, Plots
"""
Run retreat experiment (100 years from MISMIP steady state)
***Options***
    - PICO_nboxP_zQ    : PICO with P boxes, ambient temperature taken at Qm depth for P = 10, 8, 5, 2 and Q = 500, 700
                        (e.g. "PICO_nbox2_z700")
    - PME              : Plume model emulator (using Lazeroms 2018 algorithm for grounding line and basal slope -- doi:10.5194/tc-12-49-2018)
    - MISMIP_1r        : Melt rate according to the MISMIP+ ice 1r experiment (doi: 10.5194/tc-14-2283-2020) scaled to match the mean melt rate
"""

melt_rate_model = "PME"

function Favier2019_4km_retreat(melt_model)
# Grid and boundary conditions
nx = 160
ny = 20
nσ = 4
x0 = 0.0
y0 = -40000.0
dx = 4000.0
dy = 4000.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

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

# make the model
initial_conditions = InitialConditions(initial_thickness=h) # set thickness
model = Model(grid=grid,
            bed_elevation=bed,
            melt_rate = melt_model,
            initial_conditions=initial_conditions,
            solver_params=SolverParams(maxiter_picard=10) ,
            params=Params(accumulation_rate=0.3)
            ) ;

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

#output parameters
folder =  joinpath(dirname(@__FILE__),"retreat_output")
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 = 0.5
output_params = OutputParams(outputs = outputs,
                        output_freq = output_freq,
                        output_format = "jld2",
                        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

#Need these for PICO
nx = 160; ny = 20;
ice_front_mask = zeros(nx,ny);
ice_front_mask[end,:] .= 1;

if melt_rate_model == "PICO_nbox10_z700"
    melt_model = PICO(ice_front_mask = ice_front_mask,
                T0 =1.2,
                S0 = 34.6,
                use_box_mean_depth = true,
                γT = 0.94e-5,
                nbox = 10);

elseif melt_rate_model == "PICO_nbox8_z700"
    melt_model = PICO(ice_front_mask = ice_front_mask,
                T0 =1.2,
                S0 = 34.6,
                use_box_mean_depth = true,
                γT = 0.91e-5,
                nbox = 8);

elseif melt_rate_model == "PICO_nbox5_z700"
    melt_model = PICO(ice_front_mask = ice_front_mask,
                T0 = 1.2,
                S0 = 34.6,
                use_box_mean_depth = true,
                γT = 0.87e-5,
                nbox = 5);

elseif melt_rate_model == "PICO_nbox2_z700"
    melt_model = PICO(ice_front_mask = ice_front_mask,
                T0 = 1.2,
                S0 = 34.6,
                use_box_mean_depth = true,
                γT = 0.85e-5,
                nbox = 2);

elseif melt_rate_model == "PME"
    melt_model = PlumeEmulator(α=1.49)

elseif melt_rate_model == "MISMIP_1r"
    melt_model = MISMIPMeltRateOne(α = 0.32)
else
    throw(ArgumentError("Specified melt rate model not found"))
end

@time simulation = Favier2019_4km_retreat(melt_model);

#m = deepcopy(model.fields.gh.basal_melt);
#zb = model.fields.gh.b .* (model.fields.gh.grounded_fraction .== 1) + -  918.0 / 1028.0  .* model.fields.gh.h .* (model.fields.gh.grounded_fraction .< 1)
#idx = (model.fields.gh.grounded_fraction .== 0)
#mean_melt = sum(m[idx])./length(m[idx])
#println("mean melt rate for shelf points is $mean_melt")