text:
module Sequencing

using ..Qubits # Includes Clifford and Waveforms
using ..Readout
using InstrumentControl
using InstrumentControl.AWG5014C
using Alazar # For readout

export sendSequence

# From qubits
export Qubit
export cosInit
export gaussInit
export generalInit

# From waveforms
export FloatWaveform
export ExactWaveform

# From clifford
export benchmark1Qubit
export benchmark2Qubit
export Pulse

# gateNames has its last use here in naming the AWG pulse - no parents need it

#============================ Basic seqeuncing ================================#
# The sequencing also has to be different depending on which type
# of waveform we are using.  In the case of FloatWaveforms, we need to insert
# idle pulses regularly, which are all the same point in the DAC's memory.
# In the case of ExactWaveforms, we need to access distinct points in the DAC's
# memory for each combination of (line, pulse).  The code allows
# for the possibility that the qubit's lines are not all on the same board.

# In preparing for 1-qubit sequencing, the qubit ensures that all pulses in its
# dictionary are up to date, and updates the 10x3 matrix of indices.
# These map the 10 basic single qubit pulses to their corresponding on-board
# waveform index.  That means that by calling q.pulseConvert[seq, n], with some
# sequence generated elsewhere in the code, we get the actual index list for
# sequencing on the nth line.
function prepForSeq(q::Qubit)
  # Communicate with DAC only if necessary
  for p in q.waveforms
    if p[2].dirty
      q.pulseConvert[p[1],1] = pushWaveform(p[2].XYI, q.lineXYI[1], gateNames[p[1]]*"I")
      q.pulseConvert[p[1],2] = pushWaveform(p[2].XYQ, q.lineXYQ[1], gateNames[p[1]]*"Q")
      if isa(q, QubitWithZ)
        q.pulseConvert[p[1],3] = pushWaveform(p[2].Z, q.lineZ[1])
      end
      p[2].dirty = false
    end
  end
end

# Return the index of the waveform in memory
function pushWaveform(wavedata::Vector{UInt16}, board::Instrument, name="")
  if length(wavedata) == 0
    # This is one of the off channels of a FloatWaveform.  Just return the
    # stored index of idle pulse
    return idles[board]
  end
  if isa(board, InsAWG5014C)
    return awgPush(wavedata, board, name)
  else
    error("Instrument not supported yet, whine at a programmer")
  end
end

# We will be using the VISA protocol and binary writing.  We will also force
# every waveform to 250 points, with front-padded zeros.  A 30-point waveform
# will then have 220 points of offsetValue and then 30 points with signal.
function awgPush(wavedata::Vector{UInt16}, ins::InsAWG5014C, name::String)
  if length(wavedata) > 250
    error("Waveforms longer than 250 points not presently supported.")
  else
    binblockwrite(ins, "WLIST:WAV:DATA "*name*", "*reinterpret(UInt8,
      [fill(offsetValue, 250 - length(wavedata)); map(htol, wavedata)]))
  end
  name # Return the name, how the AWG will refer to it later.
end

# The following functions take pulse sequences in the form produced by the code,
# and send them to the DACs based on their indices in the DAC memory.  This
# means that sendSequence(q, [Xpi2, I, Xpi2]) for example, will send XYI, XYQ,
# and Z pulses according to the prior definitions of these pulses.
function sendSequence(q::Qubit, sequence::Pulse, readout::Bool=True)
  if size(sequence, 2) == 1
    sendSequence(q, sequence[:], readout)
  else
    error("Please specify only one sequence\n")
  end
end

function sendSequence(q::Qubit, sequence::Vector{Int8}, readout::Bool=True)
  # Make sure we aren't going to send obsolete pulses
  # (no writes done if everything is current)
  prepForSeq(q)

  # TODO: More generalized version of this function.  For now, we assume the
  # qubit is controlled by the AWG, that it is a QubitNoZ object, and that
  # the other AWG lines need to idle until the sequence terminates, and then
  # perform a readout pulse.  This allows the sequencing to be done in one
  # (albeit large) command to the AWG.  If this isn't the goal, bug Brett.
  # Someday readout will be a single pulse for multiple qubits.
  len = "SEQ:LENG "*string(length(sequence))*"\n"
  istring = map((x,y)-> "SEQ:ELEM"*string(x)*":WAV"*string(q.lineXYI[2])*
                                                    " \""*string(y)*"\"\n",
                                1:length(sequence), q.pulseConvert[sequence, 1])
  qstring = map((x,y)-> "SEQ:ELEM"*string(x)*":WAV"*string(q.lineXYQ[2])*
                                                    " \""*string(y)*"\"\n",
                                1:length(sequence), q.pulseConvert[sequence, 2])
  # Get the other channels, which are IQ for readout, and instruct them to be
  # idle at this time
  a = 1:4
  a = a[(a.!=q.lineXYI[2])&(a.!=q.lineXYQ[2])] # Select channels not used by XY
  idle1 = map(x -> "SEQ:ELEM"*string(x)*":WAV"*string(a[1])*" \"I\"\n")
  idle2 = map(x -> "SEQ:ELEM"*string(x)*":WAV"*string(a[1])*" \"I\"\n")
  # Put it all together and send to the AWG
  println(*(len, istring..., qstring..., idle1..., idle2...))
end
#=
function sendSequence(q::Qubit, sequence::Vector{Int8})
  # Make sure we aren't going to send obsolete pulses
  # (no writes done if everything is current)
  prepForSeq(q)
  # Call a lower level function for each line, using indices on the board
  # that correspond to the sequence requested
  seqLowLevel(q.pulseConvert[sequence, 1], q.lineXYI)
  seqLowLevel(q.pulseConvert[sequence, 2], q.lineXYQ)
  if isa(q, QubitWithZ)
    seqLowLevel(q.pulseConvert[sequence, 3], q.lineZ)
  end
  # Turn on output where along this path?

  # How to manage readout?  Do we send a trigger pulse along with the last
  # pulse in the seqeunce, or do we calculate the length of the pulse sequence
  # and just readout based on timing?  Also i expect not all sequences to be
  # sent are necessarily to be followed by readout, so I guess readout-
  # terminated sequences get their own method.
end
=#
# Nothing here covers readout or two-qubit gates.  Both are heirarchically
# lateral from the qubits themselves - pairs of coupled qubits need to know
# waveforms for 2-qubit gates, but not the qubits themselves.  The readout
# process covers all qubits connected to the readout line.  Whatever controls
# the readout will need to be able to tell which qubit is which, but the qubits
# themselves may remain unaware of the readout process.

end # End module sequencing


text:
init(vi::VorbisInfo) = ccall((:vorbis_info_init,vorbisenc),Void,(Ref{VorbisInfo},),vi)
clear(vi::VorbisInfo) = ccall((:vorbis_info_clear,vorbisenc),Void,(Ref{VorbisInfo},),vi)

info_blocksize(vi::VorbisInfo,zo) =
    ccall((:vorbis_info_blocksize,vorbisenc),Cint,(Ref{VorbisInfo},Cint),vi,zo)

init(vc::VorbisComment) = ccall((:vorbis_comment_init,vorbisenc),Void,(Ref{VorbisComment},Cint),vc)

comment_add(vc::VorbisComment,c::ByteString) =
    ccall((:vorbis_comment_add,vorbisenc),Void,(Ref{VorbisComment},Ptr{Uint8}),vc,c)
comment_add_tag(vc::VorbisComment,t::ByteString,c::ByteString) =
    ccall((:vorbis_comment_add_tag,vorbisenc),Void,
          (Ref{VorbisComment},Ptr{Uint8},Ptr{Uint8}),vc,t,c)

comment_query(vc::VorbisComment,t::ByteString,count) =
    bytestring(ccall((:vorbis_comment_query,vorbisenc),Ptr{Uint8},
                     (Ref{VorbisComment},Ptr{Uint8},Cint),c,t,count))

comment_query_count(vc::VorbisComment,t::ByteString) =
    ccall((:vorbis_comment_query_count,vorbisenc),Cint,
          (Ref{VorbisComment},Ptr{Uint8}),c,t)

clear(vc::VorbisComment) = ccall((:vorbis_comment_clear,vorbisenc),Void,(Ref{VorbisComment},),vc)

init(v::VorbisDSPState,vb::VorbisBlock) =
    ccall((:vorbis_block_init,vorbisenc),Cint,(Ref{VorbisDSPState},Ref{VorbisBlock}),v,vb)

clear(vb::VorbisBlock) = ccall((:vorbis_block_clear,vorbisenc),Cint,(Ref{VorbisBlock},),vb)
clear(v::VorbisDSPState) = ccall((:vorbis_dsp_clear,vorbisenc),Cint,(Ref{VorbisDSPState},),v)
init(v::VorbisDSPState,vi::VorbisInfo) =
    ccall((:vorbis_analysis_init,vorbisenc),Cint,(Ref{VorbisDSPState},Ref{VorbisInfo}),v,vi)
commentheader_out(vc::VorbisComment,op::Ogg.Packet) =
    ccall((:vorbis_commentheader_out,vorbisenc),Cint,(Ref{VorbisComment},Ref{Ogg.Packet}),vc,op)
analysis_header_out(v::VorbisDSPState,vc::VorbisComment,op::Ogg.Packet) =
    ccall((:vorbis_analysis_headerout,vorbis),Cint,
          (Ref{VorbisDSPState},Ref{VorbisComment},Ref{Ogg.Packet}),v,vc,op)
analysis_buffer(v::VorbisDSPState,vals) =
    ccall((:vorbis_analysis_buffer,vorbisenc),Ptr{Ptr{Cfloat}},
          (Ref{VorbisDSPState},Cint),v,vals)
analysis_wrote(v::VorbisDSPState,vals) =
    ccall((:vorbis_analysis_wrote,vorbisenc),Cint,
          (Ref{VorbisDSPState},Cint),v,vals)

analysis_blockout(v::VorbisDSPState,vb::VorbisBlock) =
    ccall((:vorbis_analysis_blockout,vorbisenc),Cint,
          (Ref{VorbisDSPState},Ref{VorbisBlock}),v,vb)

analysis(vb::VorbisBlock,op::Ogg.Packet) =
    ccall((:vorbis_analysis,vorbisenc),Cint,
          (Ref{VorbisBlock},Ref{Ogg.Packet}),vb,op)

bitrate_addblock(vb::VorbisBlock) =
    ccall((:vorbis_bitrate_addblock,vorbisenc),Cint,
          (Ref{VorbisBlock},),vb)

bitrate_flushpacket(v::VorbisDSPState,op::Ogg.Packet) =
    ccall((:vorbis_bitrate_flushpacket,vorbisenc),Cint,
          (Ref{VorbisDSPState},Ref{Ogg.Packet}),vb,op)

ext: jl


text:
using TracyWidom
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# See https://arxiv.org/abs/0904.1581
@test TWcdf(0,beta=1) ≈ 0.83190806620295 atol=1e-14
@test TWcdf(0,beta=2) ≈ 0.96937282835526 atol=1e-14

max_stars_repo_name: AustenLamacraft/TracyWidom

max_issues_repo_name: AustenLamacraft/TracyWidom

max_forks_repo_name: AustenLamacraft/TracyWidom


text:
# ----------------------------------------------------------------
# Using immutable types
# ----------------------------------------------------------------

using SharedArrays

@everywhere struct Point{T <: Real}
    x::T
    y::T
end
A = SharedArray{Point{Float64}}(3);
A .= [Point(rand(), rand()) for in in 1:length(A)]

#=
julia> @everywhere struct Point{T <: Real}
           x::T
           y::T
       end

julia> A = SharedArray{Point{Float64}}(3);

julia> A .= [Point(rand(), rand()) for in in 1:length(A)]
3-element SharedArray{Point{Float64},1}:
 Point{Float64}(0.07859856539270682, 0.2065974967070332)
 Point{Float64}(0.7820850739788159, 0.9971987501535848)
 Point{Float64}(0.10800356384806586, 0.24558257669174854)
=#

# ----------------------------------------------------------------
# Non-bits types (such as mutable struct) are not supported
# ----------------------------------------------------------------

@everywhere mutable struct MutablePoint{T <: Real}
    x::T
    y::T
end
B = SharedArray{MutablePoint{Float64}}(3);

#=
julia> @everywhere mutable struct MutablePoint{T <: Real}
    x::T
    y::T
end

julia> B = SharedArray{MutablePoint{Float64}}(3);
ERROR: ArgumentError: type of SharedArray elements must be bits types, got MutablePoint{Float64}
=#


text:
struct CplsrAgg
    fm
    fm_da::Plsrda
    lev
    ni
end

function cplsr_agg(X, Y, cla = nothing; ncla = nothing, nlv_da, nlv)
    X = ensure_mat(X)
    Y = ensure_mat(Y)
    if isnothing(cla)
        zfm = Clustering.kmeans(X', ncla; maxiter = 500, display = :none)
        cla = zfm.assignments
    end
    z = tab(cla)
    lev = z.keys
    nlev = length(lev)
    ni = collect(values(z))
    fm_da = plsrda(X, cla; nlv = nlv_da)
    fm = list(nlev)
    @inbounds for i = 1:nlev
        s = findall(cla .== lev[i])
        z = eval(Meta.parse(nlv))
        zmin = minimum(z)
        zmax = maximum(z)
        ni[i] <= zmin ? zmin = ni[i] - 1 : nothing
        ni[i] <= zmax ? zmax = ni[i] - 1 : nothing
        znlv = string(zmin:zmax)
        fm[i] = plsr_agg(X[s, :], Y[s, :]; nlv = znlv)
    end
    CplsrAgg(fm, fm_da, lev, ni)
end

function predict(object::CplsrAgg, X)
    X = ensure_mat(X)
    m = size(X, 1)
    nlev = length(object.lev)
    zp = predict(object.fm_da, X).posterior
    zp .= (mapreduce(i -> Float64.(zp[i, :] .== maximum(zp[i, :])), hcat, 1:m)')
    #zp = (mapreduce(i -> Jchemo.mweights(exp.(zp[i, :])), hcat, 1:m))'
    #zp .= (mapreduce(i -> 1 ./ (1 .+ exp.(-zp[i, :])), hcat, 1:m)')
    #zp .= (mapreduce(i -> zp[i, :] / sum(zp[i, :]), hcat, 1:m))'
    acc = zp[:, 1] .* predict(object.fm[1], X).pred
    @inbounds for i = 2:nlev
        if object.ni[i] >= 30
            acc .+= zp[:, i] .* predict(object.fm[i], X).pred
        end
    end
    (pred = acc, posterior = zp)
end


text:
# Functionality that may be loaded more than one

export drawTags!, drawTags

"""
    drawTags!

Draw axes and tag IDs on an image based on detections in `tags`.

Notes
- Will self detect tags, or requires user to pass in `tags` detections.
"""
function drawTags!(imageCol::AbstractArray{<:AbstractRGB,2},
                   K::AbstractArray{<:Real,2},
                   tags::AbstractArray{AprilTag,1} = AprilTagDetector()(imageCol),
                   drawReticle::Bool = false )
  #
  # draw the tag number for each tag
  if @isdefined(drawTagID!)
    foreach(x->drawTagID!(imageCol, x),tags)
  end
  #draw color box on tag corners
  foreach(tag->drawTagBox!(imageCol, tag, width = 2, drawReticle = drawReticle), tags)
  foreach(tag->drawTagAxes!(imageCol,tag, K), tags)
  imageCol
end

"""
    drawTags

Draw axes and tag IDs on an image based on detections in `tags`.

Notes
- Will self detect tags, or requires user to pass in `tags` detections.
"""
drawTags(imageCol::AbstractArray{<:AbstractRGB,2},K::AbstractArray{<:Real,2},tags::AbstractArray{AprilTag,1}=AprilTagDetector()(imageCol),drawReticle::Bool=false ) = drawTags!(RGB.(imageCol),K,tags,drawReticle)
# Convert image to RGB
# imageCol = RGB.(image)

#


text:
function wavenumber(ωs::AbstractVector{T}, medium::Medium{T}, species::Vector{Specie{T}};
        wavenumber_initial_guess = :auto, kws...) where T<:Number
    ks = Vector{Complex{T}}(length(ωs))
    ks[1] = wavenumber(ωs[1], medium, species; wavenumber_initial_guess = wavenumber_initial_guess, kws...)
    ks[2] = wavenumber(ωs[2], medium, species; wavenumber_initial_guess = wavenumber_initial_guess, kws...)
    map(3:length(ks)) do i
        k0 = ks[i-2] + (ωs[i] - ωs[i-2])*(ks[i-1] - ks[i-2])/(ωs[i-1] - ωs[i-2])
        ks[i] = wavenumber(ωs[i], medium, species; wavenumber_initial_guess = k0, kws...)
    end
    return ks
end

function wavenumber(ω::T, medium::Medium{T}, species::Vector{Specie{T}}; tol = 1e-8,
        wavenumber_initial_guess = :auto,
        hankel_order = maximum_hankel_order(ω, medium, species; tol=tol),
        radius_multiplier = 1.005,
        kws...) where T<:Number

    @warn("This function is now depricated as it's too unstable. Use wavenumbers(ω::T, medium::Medium{T}, Vector{Specie{T}}), which will return many effective wavenumbers.")
    k = ω/medium.c
    S = length(species)
    ho = hankel_order
    Z_l_n = Zn_matrix(ω, medium, species; hankel_order = ho)

    r = maximum(s.r for s in species)
    φ = sum(volume_fraction.(species))
    if wavenumber_initial_guess == :auto
        if abs(r*k) > 0.2 && φ < 0.4
            # use low volume fraction (valid for any ω) effective wavenumber
            k0 = wavenumber_low_volfrac(ω, medium, species)
        else
            # use low frequency (valid for any φ) effective wavenumber
            eff_medium = effective_medium(medium, species)
            k0 = ω/eff_medium.c
        end
    else
        k0 = wavenumber_initial_guess
    end
    initial_k0 = [real(k0), 0.95*imag(k0)]

    as = radius_multiplier*[(s1.r + s2.r) for s1 in species, s2 in species]
    function M(keff,j,l,m,n)
        (n==m ? 1.0:0.0)*(j==l ? 1.0:0.0) + 2.0pi*species[l].num_density*Z_l_n[l,n]*
            kernelN(n-m,k*as[j,l],keff*as[j,l])/(k^2.0-keff^2.0)
    end

    # this matrix is needed to calculate the eigenvectors
    MM(keff::Complex{T}) = reshape(
        [M(keff,j,l,m,n) for m in -ho:ho, j = 1:S, n in -ho:ho, l = 1:S]
    , ((2ho+1)*S, (2ho+1)*S))

    constraint(keff_vec::Array{T}) = ( (keff_vec[2] < zero(T)) ? one(T) : zero(T))*(-1 + exp(-T(100.0)*keff_vec[2]))
    detMM2(keff_vec::Array{T}) =  constraint(keff_vec) + map(x -> real(x*conj(x)), det(MM(keff_vec[1]+im*keff_vec[2])))

    #Note that there is not a unique effective wavenumber. The root closest to k_eff = 0.0 + 0.0im seems to be the right one, the others lead to strange transmission angles and large amplitudes As.
    result = optimize(detMM2, initial_k0;  g_tol = tol^2.0, f_tol = tol^4.0)
    if !(result.f_converged || result.x_converged || result.g_converged) # if nothing converged
        result = optimize(detMM2, [0.,0.];  g_tol = tol^2.0, f_tol = tol^4.0)
        if !(result.f_converged || result.x_converged || result.g_converged)
            @warn("Local optimisation did not converge")
        end
    end
    # Check result
    k_eff = result.minimizer[1] + im*result.minimizer[2]
    # in case wave travelling opposite direction was found
    if imag(k_eff) < zero(T) k_eff = - k_eff end
    MM_svd = svd(MM(k_eff))
    if last(MM_svd.S) > T(4)*tol
        @warn("Local optimisation was unsucessful at finding an effective wavenumber for ω = $ω and max(radius)*k = $(abs(r*k)).")
        @warn("$(last(MM_svd.S)) was the smallest eigenvalue value (should be zero) of the effective wavenumber matrix equation.")
    end

    return k_eff
end


text:
module AddCoverColumn

import SearchLight.Migrations: add_column, remove_column

function up()
  add_column(:books, :cover, :string)
end

function down()
  remove_column(:books, :cover)
end

end

max_stars_repo_name: tetsugps/julia

max_issues_repo_name: tetsugps/julia

max_forks_repo_name: tetsugps/julia


text:
"""
	Perform a move left
"""
function move!(x)
	if x[1] == 0
		# 0XXX
		if x[2] == 0
			# 00XX
			if x[3] == 0
				# 000A -> A000
				x[1], x[4] = x[4], 0
			elseif x[3] == x[4]
				# 00AA -> 2A000
				x[1], x[3], x[4] = 1+x[3], 0 , 0
			else
				# 00AB - AB00
				x[1], x[2], x[3], x[4] = x[3], x[4], 0, 0
			end
		else
			# 0AXX
			if x[3] == 0
				# 0A0X
				if x[4] == x[2]
					# 0A0A -> 2A000
					x[1], x[2], x[4] = 1+x[2], 0, 0
				else
					# 0A0B -> AB00
					x[1], x[2], x[4] = x[2], x[4], 0
				end
			else
				# 0AXX
				if x[3] == x[2]
					# 0AAX -> 2AX00
					x[1], x[2], x[3], x[4] = 1+x[2], x[4], 0, 0
				elseif x[3] == x[4]
					# 0ABB -> A2B00
					x[1], x[2], x[3], x[4] = x[2], 1+x[3], 0, 0
				else
					# 0ABC -> ABC0
					x[1], x[2], x[3], x[4] = x[2], x[3], x[4], 0
				end
			end
		end
	else
		# AXXX
		if x[2] == 0
			# A0XX
			if x[3] == 0
				if x[4] == x[1]
					# A00A -> a000
					x[1], x[4] = 1+x[1], 0
				else
					#A00B
					x[2], x[4] = x[4], 0
				end
			else
				# A0XX
				if x[3] == x[1]
					#A0AX -> # 2AX00
					x[1], x[2], x[3], x[4] = 1+x[1], x[4], 0, 0
				elseif x[3] == x[4]
					#A0BB -> A2B00
					x[2], x[3], x[4] = 1+x[3], 0, 0
				else
					#A0BC -> ABC0
					x[2], x[3], x[4] = x[3], x[4], 0
				end
			end
		else
			# AXXX
			if x[2] == x[1]
				# AAXX
				if x[3] == 0
					#AA0X -> aX00
					x[1], x[2], x[4] = 1+x[1], x[4], 0
				elseif x[3] == x[4]
					#AABB -> 2A2B00
					x[1], x[2], x[3], x[4] = 1+x[1], 1+x[3], 0, 0
				else
					#AABC -> 2ABC0
					x[1], x[2], x[3], x[4] = 1+x[1], x[3], x[4], 0
				end
			else
				# ABXX
				if x[3] == 0
					# AB0X -> ABX0
					x[3], x[4] = x[4], 0
				elseif x[3] == x[2]
					# ABBX -> A2BX0
					x[2], x[3], x[4] = 1+x[2], x[4], 0
				elseif x[3] == x[4]
					# ABCC -> ABc0
					x[3], x[4] = 1+x[4], 0
				#else
					#println("got here")
					# ABCD -> ABCD
					#x[1], x[2], x[3], x[4] = 1+x[1], x[3], x[4], 0
				end
			end
		end
	end
	x
end


text:
using PatModules
include("utils.jl")
@once module one
    import ..utils

    function say_hi()
        utils.say("hi")
    end
end


text:
δ(i,j) = i == j ? 1.0 : 0.0
Isym_func(i,j,k,l) = 0.5*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k))
Isymdev_func(i,j,k,l) = 0.5*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k)) - 1.0/3.0*δ(i,j)*δ(k,l)

function insert_σoop(σᵢⱼ::Tensor{2,3,T},σoop_guess) where T
  σᵢⱼ = Tensor{2,3,T}(SA[σᵢⱼ[1,1] σᵢⱼ[1,2] σᵢⱼ[1,3] ; σᵢⱼ[2,1] σᵢⱼ[2,2] σᵢⱼ[2,3] ; σᵢⱼ[3,1] σᵢⱼ[3,2] σoop_guess])
end
function insert_σoop(σᵢⱼ::SymmetricTensor{2,3,T},σoop_guess) where T
  σᵢⱼ = SymmetricTensor{2,3,T}(SA[σᵢⱼ[1,1] σᵢⱼ[1,2] σᵢⱼ[1,3] ; σᵢⱼ[2,1] σᵢⱼ[2,2] σᵢⱼ[2,3] ; σᵢⱼ[3,1] σᵢⱼ[3,2] σoop_guess])
end

function insert_into(tensor::SymmetricTensor{2,S},values,indices) where{S}
  return SymmetricTensor{2,S}((i,j) -> map_id_to_value(tensor,values,indices,i,j))
end

function map_id_to_value(tensor::SymmetricTensor,values,indices,i,j)
  ((i,j) == indices) | ((j,i) == indices) && return values
  if (i,j) in indices
    ind = findfirst(indices .== Ref((i,j)))
    return values[only(ind)]
  elseif (j,i) in indices
    ind = findfirst(indices .== Ref((j,i)))
    return values[only(ind)]
  else
    return tensor[i,j]
  end
end

function filter_offdiagonal(σᵢⱼ ; tol = 1e-5)
  (abs(σᵢⱼ[1,2]) <= tol) && (σᵢⱼ = insert_into(σᵢⱼ, 0.0, (1,2)))
  return σᵢⱼ
end

function set_print_flag!(p,iter,tsim)
  flags = p.flags
  if (iter==1) | print_flag_condition(p,iter,tsim)
    flags.print = true
    p.output.last_tsim_printed[] = tsim
  elseif flags.nan & (p.output.last_tsim_printed[] != tsim)
    flags.print = true
    p.output.last_tsim_printed[] = tsim
  else
    flags.print = false
  end
  return nothing
end

function set_save_flag!(p,iter,tsim)
  flags = p.flags
  if (iter==1) | save_flag_condition(p,iter,tsim)
    flags.save = true
    p.output.last_tsim_saved[] = tsim
  elseif flags.nan & (p.output.last_tsim_saved[] != tsim)
    flags.save = true
    p.output.last_tsim_saved[] = tsim
  else
    flags.save = false
  end
  return nothing
end

print_flag_condition(p::Params{TF,Nothing},iter,tsim) where {TF<:Real} = ( (tsim-p.output.last_tsim_printed[]) >= (1/p.output.print_frequency) )
print_flag_condition(p::Params{Nothing,TP},iter,tsim) where {TP<:Real} = (iter%p.output.print_period == 0)
function print_flag_condition(p::Params{TF,TP},iter,tsim) where {TF<:Real,TP<:Real}
  return ( (tsim-p.output.last_tsim_printed[]) >= (1/p.output.print_frequency) ) || (iter%p.output.print_period == 0)
end

save_flag_condition(p::Params{T1,T2,TF,Nothing},iter,tsim) where {T1,T2,TF<:Real} = ( (tsim-p.output.last_tsim_saved[]) >= (1/p.output.save_frequency) )
save_flag_condition(p::Params{T1,T2,Nothing,TP},iter,tsim) where {T1,T2,TP<:Real} = (iter%p.output.save_period == 0)
function save_flag_condition(p::Params{T1,T2,TF,TP},iter,tsim) where {T1,T2,TF<:Real,TP<:Real}
  return ( (tsim-p.output.last_tsim_saved[]) >= (1/p.output.save_frequency) ) || (iter%p.output.save_period == 0)
end

print_time_iteration(iter,tsim) = print("------","\n","time iteration $iter : $tsim","\n","------")

function print_nans_error(p::Params,Δt)
    flags = p.flags
    if flags.nan1
        flags.nan = true
        @error "integration over Δt=$(Δt) produced NaNs."
    elseif flags.nan2
        flags.nan = true
        @error "integration over Δt/2=$(Δt/2) twice produced NaNs."
    end
end


text:
# ---
# title: 539. Minimum Time Difference
# id: problem539
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/minimum-time-difference/description/>
# hidden: true
# ---
#
# Given a list of 24-hour clock time points in **" HH:MM"** format, return _the
# minimum **minutes** difference between any two time-points in the list_.
#
#
#
# **Example 1:**
#
#
#
#     Input: timePoints = ["23:59","00:00"]
#     Output: 1
#
#
# **Example 2:**
#
#
#
#     Input: timePoints = ["00:00","23:59","00:00"]
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `2 <= timePoints <= 2 * 104`
#   * `timePoints[i]` is in the format **" HH:MM"**.
#
#
## @lc code=start
using LeetCode

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


text:
using ImageView, TestImages, ImageCore, Reactive,
      GtkReactive, Gtk, IntervalSets
using Test
import AxisArrays
using AxisArrays: Axis

@testset "1d" begin
    img = rand(N0f8, 5)
    guidict = imshow_now(img)
    win = guidict["gui"]["window"]
    destroy(win)
end

@testset "Aspect ratio" begin
    img = rand(N0f8, 20, 20)
    guidict = imshow_now(img)
    win, frame = guidict["gui"]["window"], guidict["gui"]["frame"]
    @test isa(frame, Gtk.GtkAspectFrameLeaf)
    zr = guidict["roi"]["zoomregion"]
    Reactive.run_till_now()
    sleep(1.0)  # give compilation a chance to catch up

    @test get_gtk_property(frame, :ratio, Float32) == 1.0
    push!(zr, (1:20, 8:10))  # The first one sometimes gets dropped. Reactive bug?
    Reactive.run_till_now()
    push!(zr, (1:20, 9:10))
    Reactive.run_till_now()
    sleep(1.0)
    @test value(zr).currentview.x == 9..10
    @test get_gtk_property(frame, :ratio, Float32) ≈ 0.1
    push!(zr, (9:10, 1:20))
    Reactive.run_till_now()
    Gtk.showall(win)
    sleep(0.1)
    @test get_gtk_property(frame, :ratio, Float32) ≈ 10.0

    destroy(win)

    guidict = imshow_now(img, aspect=:none)
    win, frame = guidict["gui"]["window"], guidict["gui"]["frame"]
    @test isa(frame, Gtk.GtkFrameLeaf)
    destroy(win)
end

# image display
@testset "Image display" begin
    img_n0f8 = rand(N0f8, 3,3)
    imsd = imshow_now(img_n0f8; name="N0f8")
    @test get_gtk_property(imsd["gui"]["window"], :title, String) == "N0f8"

    img_n0f16 = rand(N0f16, 3,3)
    imshow_now(img_n0f16; name="N0f16")

    img_rgb = rand(RGB{N0f8}, 3, 3)
    imshow_now(img_rgb; name="RGB{N0f8}")

    img_int = rand(Int, 3,3)
    imshow_now(img_int; name="Int")

    img_float16 = rand(Float16, 3,3)
    imshow_now(img_float16; name="Float16")

    img_float32 = rand(Float32, 3,3)
    img_float32[1,1] = NaN
    img_float32[2,1] = Inf
    img_float32[3,1] = -5
    imshow_now(img_float32; name="Float32")

    img_float64 = rand(Float64, 3,3)
    imshow_now(img_float64; name="Float64")

    img_nan = fill(NaN, (3,3))
    imshow_now(img_nan; name="NaN")

    img_rgbfloat = rand(RGB{Float32}, 3, 3)
    imshow_now(img_rgbfloat; name="RGB{Float32}")

    img = testimage("lighthouse")
    hlh = imshow_now(img, name="Lighthouse")

    # a large image
    img = testimage("earth")
    hbig = imshow_now(img, name="Earth")
    win = hbig["gui"]["window"]
    w, h = size(win)
    ws, hs = screen_size(win)
    !Sys.iswindows() && @test w <= ws && h <= hs
end

@testset "imshow!" begin
    img = testimage("mri")
    guidict = imshow(img[:,:,1])
    c = guidict["gui"]["canvas"]
    ImageView.imshow!(c, img[:,:,2])

    imgsig = Signal(img[:,:,1])
    imshow(c, imgsig)
    push!(imgsig, img[:,:,8])
end

@testset "Orientation" begin
    img = [1 2; 3 4]
    guidict = imshow_now(img)
    @test parent(value(guidict["roi"]["image roi"])) == [1 2; 3 4]
    guidict = imshow_now(img, flipy=true)
    @test parent(value(guidict["roi"]["image roi"])) == [3 4; 1 2]
    guidict = imshow_now(img, flipx=true)
    @test parent(value(guidict["roi"]["image roi"])) == [2 1; 4 3]
    guidict = imshow_now(img, flipx=true, flipy=true)
    @test parent(value(guidict["roi"]["image roi"])) == [4 3; 2 1]
end

if Gtk.libgtk_version >= v"3.10"
    # These tests use the player widget
    @testset "Multidimensional" begin
        # Test that we can use positional or named axes with AxisArrays
        img = AxisArray(rand(3, 5, 2), :x, :y, :z)
        guin = imshow_now(img; name="AxisArray Named")
        @test isa(guin["roi"]["slicedata"].axs[1], Axis{:z})
        guip = imshow_now(img; axes=(1,2), name="AxisArray Positional")
        @test isa(guip["roi"]["slicedata"].axs[1], Axis{3})
        guip2 = imshow_now(img; axes=(1,3), name="AxisArray Positional")
        @test isa(guip2["roi"]["slicedata"].axs[1], Axis{2})

        ## 3d images
        img = testimage("mri")
        hmri = imshow_now(img; name="P,R view")
        @test isa(hmri["roi"]["slicedata"].axs[1], Axis{:S})

        # Use a custom CLim here because the first slice is not representative of the intensities
        hmrip = imshow_now(img, Signal(CLim(0.0, 1.0)), axes=(:S, :P), name="S,P view")
        @test isa(hmrip["roi"]["slicedata"].axs[1], Axis{:R})
        push!(hmrip["roi"]["slicedata"].signals[1], 84)

        ## Two coupled images
        mriseg = RGB.(img)
        mriseg[img .> 0.5] .= colorant"red"
        # version 1
        guidata = imshow_now(img, axes=(1,2))
        zr = guidata["roi"]["zoomregion"]
        slicedata = guidata["roi"]["slicedata"]
        guidata2 = imshow_now(mriseg, nothing, zr, slicedata)
        @test guidata2["roi"]["zoomregion"] === zr

        # version 2
        zr, slicedata = roi(img, (1,2))
        gd = imshow_gui((200, 200), (1,2); slicedata=slicedata)
        guidata1 = imshow(gd["frame"][1,1], gd["canvas"][1,1], img, nothing, zr, slicedata)
        guidata2 = imshow(gd["frame"][1,2], gd["canvas"][1,2], mriseg, nothing, zr, slicedata)
        Gtk.showall(gd["window"])
        sleep(0.01)
        @test guidata1["zoomregion"] === guidata2["zoomregion"] === zr

        # imlink
        gd = imlink(img, mriseg)
        Gtk.showall(gd["window"])
        sleep(0.01)
        @test gd["guidata"][1]["zoomregion"] === gd["guidata"][2]["zoomregion"]
    end

    @testset "Non-AbstractArrays" begin
        include("cone.jl")
    end
end

nothing


text:
using
  GeophysicalFlows,
  Statistics,
  Random,
  Test

import # use 'import' rather than 'using' for submodules to keep namespace clean
  GeophysicalFlows.TwoDNavierStokes,
  GeophysicalFlows.SingleLayerQG,
  GeophysicalFlows.BarotropicQGQL,
  GeophysicalFlows.MultiLayerQG,
  GeophysicalFlows.SurfaceQG

using FourierFlows: parsevalsum
using GeophysicalFlows: lambdipole, peakedisotropicspectrum

# the devices on which tests will run
devices = CUDA.has_cuda() ? (CPU(), GPU()) : (CPU(),)

const rtol_lambdipole = 1e-2 # tolerance for lamb dipole tests
const rtol_twodnavierstokes = 1e-13 # tolerance for twodnavierstokes forcing tests
const rtol_singlelayerqg = 1e-13 # tolerance for singlelayerqg forcing tests
const rtol_multilayerqg = 1e-13 # tolerance for multilayerqg forcing tests
const rtol_surfaceqg = 1e-13 # tolerance for surfaceqg forcing tests

# Run tests
testtime = @elapsed begin
for dev in devices

  @info "testing on " * string(typeof(dev))

  @testset "Utils" begin
    include("test_utils.jl")

    @test testpeakedisotropicspectrum(dev)
    @test_throws ErrorException("the domain is not square") testpeakedisotropicspectrum_rectangledomain()
  end

  @testset "TwoDNavierStokes" begin
    include("test_twodnavierstokes.jl")

    @test test_twodnavierstokes_advection(0.0005, "ForwardEuler", dev)
    @test test_twodnavierstokes_lambdipole(256, 1e-3, dev)
    @test test_twodnavierstokes_deterministicforcing_energybudget(dev)
    @test test_twodnavierstokes_stochasticforcing_energybudget(dev)
    @test test_twodnavierstokes_deterministicforcing_enstrophybudget(dev)
    @test test_twodnavierstokes_stochasticforcing_enstrophybudget(dev)
    @test test_twodnavierstokes_energyenstrophy(dev)
    @test test_twodnavierstokes_problemtype(dev, Float32)
    @test TwoDNavierStokes.nothingfunction() == nothing
  end

  @testset "SingleLayerQG" begin
    include("test_singlelayerqg.jl")

    for deformation_radius in [Inf, 1.23]
      @test test_1layerqg_rossbywave("ETDRK4", 1e-2, 20, dev, deformation_radius=deformation_radius)
      @test test_1layerqg_rossbywave("FilteredETDRK4", 1e-2, 20, dev, deformation_radius=deformation_radius)
      @test test_1layerqg_rossbywave("RK4", 1e-2, 20, dev, deformation_radius=deformation_radius)
      @test test_1layerqg_rossbywave("FilteredRK4", 1e-2, 20, dev, deformation_radius=deformation_radius)
      @test test_1layerqg_rossbywave("AB3", 1e-3, 200, dev, deformation_radius=deformation_radius)
      @test test_1layerqg_rossbywave("FilteredAB3", 1e-3, 200, dev, deformation_radius=deformation_radius)
      @test test_1layerqg_rossbywave("ForwardEuler", 1e-4, 2000, dev, deformation_radius=deformation_radius)
      @test test_1layerqg_rossbywave("FilteredForwardEuler", 1e-4, 2000, dev, deformation_radius=deformation_radius)
      @test test_1layerqg_problemtype(dev, Float32, deformation_radius=deformation_radius)
    end
    @test test_1layerqg_advection(0.0005, "ForwardEuler", dev)
    @test test_streamfunctionfrompv(dev; deformation_radius=1.23)
    @test test_1layerqg_energies_EquivalentBarotropicQG(dev; deformation_radius=1.23)
    @test test_1layerqg_energyenstrophy_BarotropicQG(dev)
    @test test_1layerqg_deterministicforcing_energybudget(dev)
    @test test_1layerqg_stochasticforcing_energybudget(dev)
    @test test_1layerqg_deterministicforcing_enstrophybudget(dev)
    @test test_1layerqg_stochasticforcing_enstrophybudget(dev)
    @test SingleLayerQG.nothingfunction() == nothing
    @test_throws ErrorException("not implemented for finite deformation radius") test_1layerqg_energy_dissipation(dev; deformation_radius=2.23)
    @test_throws ErrorException("not implemented for finite deformation radius") test_1layerqg_enstrophy_dissipation(dev; deformation_radius=2.23)
    @test_throws ErrorException("not implemented for finite deformation radius") test_1layerqg_energy_work(dev; deformation_radius=2.23)
    @test_throws ErrorException("not implemented for finite deformation radius") test_1layerqg_enstrophy_work(dev; deformation_radius=2.23)
    @test_throws ErrorException("not implemented for finite deformation radius") test_1layerqg_energy_drag(dev; deformation_radius=2.23)
    @test_throws ErrorException("not implemented for finite deformation radius") test_1layerqg_enstrophy_drag(dev; deformation_radius=2.23)
  end

  @testset "BarotropicQGQL" begin
    include("test_barotropicqgql.jl")

    @test test_bqgql_rossbywave("ETDRK4", 1e-2, 20, dev)
    @test test_bqgql_rossbywave("FilteredETDRK4", 1e-2, 20, dev)
    @test test_bqgql_rossbywave("RK4", 1e-2, 20, dev)
    @test test_bqgql_rossbywave("FilteredRK4", 1e-2, 20, dev)
    @test test_bqgql_rossbywave("AB3", 1e-3, 200, dev)
    @test test_bqgql_rossbywave("FilteredAB3", 1e-3, 200, dev)
    @test test_bqgql_rossbywave("ForwardEuler", 1e-4, 2000, dev)
    @test test_bqgql_rossbywave("FilteredForwardEuler", 1e-4, 2000, dev)
    @test test_bqgql_deterministicforcingbudgets(dev)
    @test test_bqgql_stochasticforcingbudgets(dev)
    @test test_bqgql_advection(0.0005, "ForwardEuler", dev)
    @test test_bqgql_energyenstrophy(dev)
    @test test_bqgql_problemtype(dev, Float32)
    @test BarotropicQGQL.nothingfunction() == nothing
  end

  @testset "SurfaceQG" begin
    include("test_surfaceqg.jl")

    @test test_sqg_kineticenergy_buoyancyvariance(dev)
    @test test_sqg_advection(0.0005, "ForwardEuler", dev)
    @test test_sqg_deterministicforcing_buoyancy_variance_budget(dev)
    @test test_sqg_stochasticforcing_buoyancy_variance_budget(dev)
    @test test_sqg_stochasticforcedproblemconstructor(dev)
    @test test_sqg_problemtype(dev, Float32)
    @test test_sqg_paramsconstructor(dev)
    @test test_sqg_noforcing(dev)
    @test SurfaceQG.nothingfunction() == nothing
  end

  @testset "MultiLayerQG" begin
    include("test_multilayerqg.jl")

    @test test_pvtofromstreamfunction_2layer(dev)
    @test test_pvtofromstreamfunction_3layer(dev)
    @test test_mqg_rossbywave("RK4", 1e-2, 20, dev)
    @test test_mqg_nonlinearadvection(0.005, "ForwardEuler", dev)
    @test test_mqg_linearadvection(0.005, "ForwardEuler", dev)
    @test test_mqg_energies(dev)
    @test test_mqg_energysinglelayer(dev)
    @test test_mqg_fluxes(dev)
    @test test_mqg_fluxessinglelayer(dev)
    @test test_mqg_setqsetψ(dev)
    @test test_mqg_paramsconstructor(dev)
    @test test_mqg_stochasticforcedproblemconstructor(dev)
    @test test_mqg_problemtype(dev, Float32)
    @test MultiLayerQG.nothingfunction() == nothing
  end
end
end # time

println("Total test time: ", testtime)


text:
using CellFreeModelGenerationKit
using Logging

# path to test metabolism file -
path_to_test_metabolism_file = "./test/data/Test.vff"
path_to_output_dir = "./test/generated_code/model"

# setup logging -
io = open("log.txt", "w+")
logger = SimpleLogger(io)

# need to build a model object -
build_result = build_julia_model_object(path_to_test_metabolism_file, path_to_output_dir)
if (isa(build_result.value, Exception) == true)
    @show build_result.value
end
julia_model_object = build_result.value

# generate code -
generate(julia_model_object; logger=logger)
flush(io)
close(io)

max_stars_repo_name: varnerlab/CellFreeModelGenerationKit

max_issues_repo_name: varnerlab/CellFreeModelGenerationKit

max_forks_repo_name: varnerlab/CellFreeModelGenerationKit


text:
using PyCall
using SparseArrays
using LinearAlgebra
using Flux.Optimise
import Flux.Tracker: data, grad, back!, update!

const scipy_sparse_find = pyimport("scipy.sparse")["find"]
function jlsparse(Apy::PyObject)
    IA, JA, SA = scipy_sparse_find(Apy)
    return sparse(Int[i+1 for i in IA], Int[i+1 for i in JA], SA)
end

function jlload_data(dataset_str)
    # Loads input data from gcn/data directory

    py"""
import numpy as np
import scipy.sparse as sp
# import torch

def encode_onehot(labels):
    classes = set(labels)
    classes_dict = {c: np.identity(len(classes))[i, :] for i, c in
                    enumerate(classes)}
    labels_onehot = np.array(list(map(classes_dict.get, labels)),
                             dtype=np.int32)
    return labels_onehot

def load_data(path="./data/cora/", dataset="cora"):

    print('Loading {} dataset...'.format(dataset))

    idx_features_labels = np.genfromtxt("{}{}.content".format(path, dataset),
                                        dtype=np.dtype(str))
    features = sp.csr_matrix(idx_features_labels[:, 1:-1], dtype=np.float32)
    labels = encode_onehot(idx_features_labels[:, -1])

    # build graph
    idx = np.array(idx_features_labels[:, 0], dtype=np.int32)
    idx_map = {j: i for i, j in enumerate(idx)}
    edges_unordered = np.genfromtxt("{}{}.cites".format(path, dataset),
                                    dtype=np.int32)
    edges = np.array(list(map(idx_map.get, edges_unordered.flatten())),
                     dtype=np.int32).reshape(edges_unordered.shape)
    adj = sp.coo_matrix((np.ones(edges.shape[0]), (edges[:, 0], edges[:, 1])),
                        shape=(labels.shape[0], labels.shape[0]),
                        dtype=np.float32)

    # build symmetric adjacency matrix
    adj = adj + adj.T.multiply(adj.T > adj) - adj.multiply(adj.T > adj)

    features = normalize(features)
    adj = normalize(adj + sp.eye(adj.shape[0]))

    # idx_train = range(140)
    # idx_val = range(200, 500)
    # idx_test = range(500, 1500)

    return adj, features, labels

def normalize(mx):
    rowsum = np.array(mx.sum(1))
    r_inv = np.power(rowsum, -1).flatten()
    r_inv[np.isinf(r_inv)] = 0.
    r_mat_inv = sp.diags(r_inv)
    mx = r_mat_inv.dot(mx)
    return mx

def accuracy(output, labels):
    preds = output.max(1)[1].type_as(labels)
    correct = preds.eq(labels).double()
    correct = correct.sum()
    return correct / len(labels)

    """

    adj, features, labels = py"load_data"()

    features = jlsparse(features)
    adj = jlsparse(adj)
    idx_train = 1:140
    idx_val = 141:500
    idx_test = 501:size(features, 2)
    return adj, features, labels, idx_train, idx_val, idx_test
end

function jlaccuracy(preds, labels)
    # accuracy
    correct_pred = equal(argmax(preds, 1), argmax(labels, 1))
    accuracy_all = cast(correct_pred, Float32)
    return reduce_mean(accuracy_all)
end

macro interrupts(ex)
  :(try $(esc(ex))
    catch e
      e isa InterruptException || rethrow()
      throw(e)
    end)
end

function update!(opt, x, x̄)
  update!(x, apply!(opt, x, copy(data(x̄))))
end

function _update_params!(opt, xs)
  for x in xs
    Δ = Optimise.apply!(opt, x.data, x.grad)
    x.data .-= Δ
    Δ .= 0
  end
end

max_stars_repo_name: mcbaron/fluxGCN

max_issues_repo_name: mcbaron/fluxGCN

max_forks_repo_name: mcbaron/fluxGCN


text:
using Epidemia
using Test

@testset "Epidemia.jl" begin
    # Write your tests here.
end


text:
"""
    points(A::AbstractRaster; dims=(YDim, XDim), ignore_missing) => Array{Tuple}

Returns a generator of the points in `A` for dimensions in `dims`,
where points are a tuple of the values in each specified dimension
index.

# Keywords

- `dims` the dimensions to return points from. The first slice of other
    layers will be used.
- `ignore_missing`: wether to ignore missing values in the array when considering
    points. If `true`, all points in the dimensions will be returned, if `false`
    only the points that are not `=== missingval(A)` will be returned.

The order of `dims` determines the order of the points.

$EXPERIMENTAL
"""
function points(A::AbstractRaster; ignore_missing=false, order=(XDim, YDim, ZDim))
    ignore_missing ? _points(A; order) : _points_missing(A; order)
end
function points(dims::DimTuple; order=(XDim, YDim, ZDim))
    return DimPoints(dims; order=DD.dims(dims, order))
end

_points(A::AbstractRaster; kw...) = points(dims(A); kw...)
function _points_missing(A::AbstractRaster; order)
    points = DimPoints(A; order=dims(A, order))
    function ordered_point_or_missing(I)
        isequal(A[I], missingval(A)) ? missing : points[I]
    end
    return (ordered_point_or_missing(I) for I in CartesianIndices(A))
end


text:
struct SLk0k1lMixingRule <: SLMixingRule
    components::Vector{String}
    k0::PairParam{Float64}
    k1::PairParam{Float64}
    l::PairParam{Float64}
end

@registermodel SLk0k1lMixingRule

function SLk0k1lMixingRule(components; userlocations=String[], verbose=false)
    params = getparams(components, ["LatticeFluid/SanchezLacombe/mixing/k0k1l_unlike.csv"]; userlocations=userlocations, verbose=verbose)
    k0 = params["k0"]
    k1 = params["k1"]
    l = params["l"]
    model = SLk0k1lMixingRule(components,k0,k1,l)
    return model
end

function sl_mix(unmixed_vol,unmixed_epsilon,mixmodel::SLk0k1lMixingRule)
    #dont mind the function names, it performs the correct mixing
    premixed_vol= epsilon_LorentzBerthelot(unmixed_vol,mixmodel.l)
    premixed_epsilon = sigma_LorentzBerthelot(unmixed_epsilon)
    return premixed_vol,premixed_epsilon
end

function mix_vε(model::SanchezLacombe,V,T,z,mix::SLk0k1lMixingRule,r̄,Σz)
    r =  model.params.segment.values
    v = model.params.vol.values
    k0 = mix.k0.values
    k1 = mix.k1.values
    ε = model.params.epsilon.values
    r̄inv = one(r̄)/r̄
     ϕ = @. r* z* r̄inv/Σz
    v_r = zero(V+T+first(z))
    ε_r = v_r
    Σz2 = 1/(Σz*Σz)
    for i in @comps
        for j in @comps
            ϕi = ϕ[i]
            ϕj = ϕ[j]
            ϕiϕj = ϕi*ϕj
            v_r += ϕiϕj*v[i,j]
            δij = (i===j)
            kmi = view(k1,:,i)
            kmj = view(k1,:,j)
            kij = k0[i,j] + (1-δij)* (dot(ϕ,kmi) + dot(ϕ,kmj))

            εij = ε[i,j]*(1-kij)
            ε_r += ϕiϕj*εij
        end
    end
    return v_r,ε_r
end


text:
# Control script for calculating CPCP for multiple timesteps.
# Using Channels for parallel tasks.
#
# To run this on a server without x-server, you must switch to the headless
# version of Paraview.
#
# Hongyang Zhou, hyzhou@umich.edu 03/19/2020

## Import packages
using Distributed, Glob
@everywhere using DelimitedFiles, ParallelDataTransfer
if !isdefined(Main, :ReconnectionRate)
   @everywhere include("utility.jl")
   @everywhere using .ReconnectionRate
end

## Function and global definitions

# The first parameter of a channel is the Type (e.g.: Int) and the second one,
# the maximum number of elements allowed in the channel.
const jobs    = RemoteChannel(()->Channel{Int}(32))
const results = RemoteChannel(()->Channel{Tuple}(32))

@everywhere function process_file(i)
   file = files[i]
   @info "processing file $file ..."
   # Create Python script with the correct input filename for tracing
   txt = read(filePy_orig[1], String)
   tracePy = filePy_orig[1][1:end-3]*string(i)*".py"
   fileTrace = "streamline"*lpad(i,4,"0")*".txt"
   fileStream = joinpath(fileDir, fileTrace)
   fileBoundaryXYZ = "boundary"*lpad(i,4,"0")*".txt"

   if isfile(fileDir*fileBoundaryXYZ)
      return # skip the finished snapshots
   end

   open(tracePy, "w") do f
      txt = replace(txt, "test.vtu"=>file)
      txt = replace(txt, "streamline.txt"=>fileTrace)
      write(f, txt)
   end

   # Step 2: field line tracing (in Paraview)
   run(`$run_pvbatch $tracePy`)
   rm(tracePy, force=true)

   # Step 3: find magnetosphere boundary
   x, y, z = findBoundary(fileStream, method=1, nSmooth=100, DoPlot=false)
   saveBoundary(x, y, z, joinpath(fileDir,fileBoundaryXYZ))
   rm(fileDir*fileTrace, force=true)

   # Create Python script with the correct input filename for interpolation
   txt = read(filePy_orig[2], String)
   interpPy = filePy_orig[2][1:end-3]*string(i)*".py"
   fileBoundary = "boundary_value"*lpad(i,4,"0")*".txt"
   open(interpPy, "w") do f
      txt = replace(txt, "test.vtu" => file)
      txt = replace(txt, "boundary.txt" => fileBoundaryXYZ)
      txt = replace(txt, "boundary_value.txt" => fileBoundary)
      write(f, txt)
   end

   # Step 4: interpolate onto the boundary points (in Paraview)
   run(`$run_pvbatch $interpPy`)
   rm(interpPy)
end

function write_CPCP(filenames, fileDir="./")

   ϕ = Vector{Float32}(undef,length(filenames))

   for (i,file) in enumerate(filenames)
      fileBoundary = "boundary_value"*lpad(i,4,"0")*".txt"
      # Step 5: 2D line integral
      file_integration = joinpath(fileDir, fileBoundary)
      ϕ[i] = integrate_along_boundary_hall(file_integration, false)

      open(joinpath(fileDir,"CPCP.txt"), "a") do io
         i_end  = findfirst("_n",file)[1] - 1
         second = parse(Int32, file[i_end-1:i_end])
         minute = parse(Int32, file[i_end-3:i_end-2])
         timestep = 60*minute + second
         writedlm(io, [timestep round(ϕ[i],digits=4)])
      end
   end

end

@everywhere function do_work(jobs, results) # Define work function everywhere.
   while true
      job_id = take!(jobs)
      exec_time = @elapsed process_file(job_id)
      put!(results, (job_id, exec_time, myid()))
   end
end

function make_jobs(n)
   for i in 1:n
      put!(jobs, i)
   end
end

## Main
if Sys.isapple() # Sys.MACHINE might be better
   run_pvbatch = "pvbatch"
else
   run_pvbatch = "run_pvbatch"
end
# Step 1: create seeding points (same for all timesteps)
@everywhere dataDir = "../../GM/"
@everywhere fileDir = "data/"

!isdir(fileDir) && mkdir(fileDir)

createSeeds("seeds.txt", fileDir)
# Loop over multiple snapshots
files = basename.(glob("cut*.vtu", dataDir))

sendto(workers(), files=files)

@everywhere filePy_orig = ["streamtrace_save.py", "Interpolate_On_Boundary_batch.py"]

nFile = length(files)

@async make_jobs(nFile) # Feed the jobs channel with "n" jobs.

for p in workers() # Start tasks on the workers to process requests in parallel.
   @async remote_do(do_work, p, jobs, results) # Similar to remotecall.
end

@elapsed while nFile > 0 # Print out results.
   job_id, exec_time, wid = take!(results)
   println("$job_id finished $(round(exec_time; digits=2)) seconds "*"
   on worker $wid")
   global nFile -= 1
end

write_CPCP(files, fileDir)


text:
module FESpacesTests1

using Test

@testset "ConformingFESpaces" begin include("ConformingFESpacesTests.jl") end

@testset "FESpacesInterfaces" begin include("FESpacesInterfacesTests.jl") end

@testset "SingleFieldFESpaces" begin include("SingleFieldFESpacesTests.jl") end

@testset "TrialFESpaces" begin include("TrialFESpacesTests.jl") end

@testset "Assemblers" begin include("AssemblersTests.jl") end

@testset "SparseMatrixAssemblers" begin include("SparseMatrixAssemblersTests.jl") end

@testset "FEOperators" begin include("FEOperatorsTests.jl") end

@testset "AffineFEOperators" begin include("AffineFEOperatorsTests.jl") end

@testset "FEOperatorsFromWeakForm" begin include("FEOperatorsFromWeakFormTests.jl") end

@testset "FESolvers" begin include("FESolversTests.jl") end

@testset "DiscontinuousFESpaces" begin include("DiscontinuousFESpacesTests.jl") end

@testset "DivConformingFESpaces" begin include("DivConformingFESpacesTests.jl") end

@testset "CurlConformingFESpaces" begin include("CurlConformingFESpacesTests.jl") end

end # module


text:
export log_pdf, gradlog_pdf, hesslog_pdf

"""
    log_pdf(x)

Compute the log of the density of the univariate standard Gaussian distribution (zero mean and unitary standard deviation).
"""
log_pdf(x) = -0.5*(log(2*π) + x^2)

"""
    gradlog_pdf(x)

Compute the gradient of the log of the density of the univariate standard Gaussian distribution (zero mean and unitary standard deviation).
"""
gradlog_pdf(x) = -x

"""
    hesslog_pdf(x)

Compute the hessian of the log of the density of the univariate standard Gaussian distribution (zero mean and unitary standard deviation).
"""
hesslog_pdf(x) = -1.0


text:
# Linear Regression with L1 penalty
PlotResults=false

useproc("GPU") # GPU is significantly faster for larger systems
#useproc("CPU")

StartCode()
X=ADnode()
w=ADvariable()
Y=ADnode()
loss=meanSquareLoss(X*w,Y) # mixed regularisation term
net=EndCode() # defines the graph

# instantiate parameter nodes and inputs:

D=100 # dimension of weight vector
N=100 # number of datapoints
truew=randn(D,1).*(rand(D,1).>0.8)
net.value[w]=truew
net.value[X]=randn(N,D)
net.value[Y]=net.value[X]*net.value[w] # make a realisable problem
net.value[w]=randn(D,1) # start with a random w

net=compile(net) # compile and preallocate memory

@gpu CUDArt.init([0]) # let the user do device management

@gpu net=convert(net,PROC,Float32)
#@gpu net=convert(net,PROC,Float64)

if 1==1

gradcheck(net)

# Training:
    tic()
println("Training: using $(net.gpu==true? "GPU" : "CPU") ")
ParsToUpdate=Parameters(net)
nupdates=500
LearningRate=0.01
error=Array(Float64,0)
velo=NesterovInit(net) # Nesterov velocity
for i=1:nupdates
    ADforward!(net)
    ADbackward!(net)
    push!(error,extract(net.value[net.FunctionNode])[1])
    printover("iteration $i: training loss = $(error[i])")
    for par in ParsToUpdate
        #GradientDescentUpdate!(net.value[par],net.gradient[par],LearningRate)
        NesterovGradientDescentUpdate!(net.value[par],net.gradient[par],velo[par],LearningRate,i)
    end
end
@gpu CUDArt.device_synchronize()
println()
    toc()

#println("\ntrue weights   ",round(100*truew')/100)
#println("learned weights",round(100*extract(net.value[w])')/100)

if PlotResults
    figure(1);plot(truew);title("true weights")
    figure(2);plot(extract(net.value[w])); title("learned weights")
end

end

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
function toggle(::UIkit; value=false, label="", bind="checked", valueUpdate="change", kwargs...)
    ip = input(NativeHTML(), value; typ="checkbox", wdgtyp="toggle", bind=bind, valueUpdate="change", kwargs...)
    InteractBase.scope(ip).dom = dom"div.field"(
        Node(:label, className = "switch-light switch-material", attributes = Dict("onclick"=>""))(
            InteractBase.scope(ip).dom,
            label,
            Node(:span, Node(:span, "Off"), Node(:span, "On"), Node(:a))
        )
    )
    Widget{:toggle}(ip)
end


text:
module Exp
include("./observable.jl")

experiments = Array{Dict{String,Array{Float64,1}},1}(undef, length(observables))
error_bars = Array{Dict{String,Array{Float64,1}},1}(undef, length(observables))

experiments[observables_index("Phosphorylated_MEKc")] = Dict(
    "EGF" => [0.000,0.773,0.439,0.252,0.130,0.087,0.080,0.066],
    "HRG" => [0.000,0.865,1.000,0.837,0.884,0.920,0.875,0.789],
)
error_bars[observables_index("Phosphorylated_MEKc")] = Dict(
    "EGF" => [0.000,0.030,0.048,0.009,0.009,0.017,0.012,0.008] ./ sqrt(3),
    "HRG" => [0.000,0.041,0.000,0.051,0.058,0.097,0.157,0.136] ./ sqrt(3),
)

experiments[observables_index("Phosphorylated_ERKc")] = Dict(
    "EGF" => [0.000,0.867,0.799,0.494,0.313,0.266,0.200,0.194],
    "HRG" => [0.000,0.848,1.000,0.971,0.950,0.812,0.747,0.595],
)
error_bars[observables_index("Phosphorylated_ERKc")] = Dict(
    "EGF" => [0.000,0.137,0.188,0.126,0.096,0.087,0.056,0.012] ./ sqrt(3),
    "HRG" => [0.000,0.120,0.000,0.037,0.088,0.019,0.093,0.075] ./ sqrt(3),
)

experiments[observables_index("Phosphorylated_RSKw")] = Dict(
    "EGF" => [0,0.814,0.812,0.450,0.151,0.059,0.038,0.030],
    "HRG" => [0,0.953,1.000,0.844,0.935,0.868,0.779,0.558],
)
error_bars[observables_index("Phosphorylated_RSKw")] = Dict(
    "EGF" => [0,0.064,0.194,0.030,0.027,0.031,0.043,0.051] ./ sqrt(3),
    "HRG" => [0,0.230,0.118,0.058,0.041,0.076,0.090,0.077] ./ sqrt(3),
)

experiments[observables_index("Phosphorylated_cFos")] = Dict(
    "EGF" => [0,0.060,0.109,0.083,0.068,0.049,0.027,0.017],
    "HRG" => [0,0.145,0.177,0.158,0.598,1.000,0.852,0.431],
)
error_bars[observables_index("Phosphorylated_cFos")] = Dict(
    "EGF" => [0,0.003,0.021,0.013,0.016,0.007,0.003,0.002] ./ sqrt(3),
    "HRG" => [0,0.010,0.013,0.001,0.014,0.000,0.077,0.047] ./ sqrt(3),
)

experiments[observables_index("Phosphorylated_CREBw")] = Dict(
    "EGF" => [0,0.446,0.030,0.000,0.000],
    "HRG" => [0,1.000,0.668,0.460,0.340],
)
error_bars[observables_index("Phosphorylated_CREBw")] = Dict(
    "EGF" => [0,0.0,0.0,0.0,0.0] ./ sqrt(3),
    "HRG" => [0,0.0,0.0,0.0,0.0] ./ sqrt(3),
)

experiments[observables_index("cfos_mRNA")] = Dict(
    "EGF" => [0,0.181,0.476,0.518,0.174,0.026,0.000],
    "HRG" => [0,0.353,0.861,1.000,0.637,0.300,0.059],
)
error_bars[observables_index("cfos_mRNA")] = Dict(
    "EGF" => [0.017,0.004,0.044,0.004,0.023,0.007,0.008] ./ sqrt(3),
    "HRG" => [0.017,0.006,0.065,0.044,0.087,0.023,0.001] ./ sqrt(3),
)

experiments[observables_index("cFos_Protein")] = Dict(
    "EGF" => [0,0.078,0.216,0.240,0.320,0.235],
    "HRG" => [0,0.089,0.552,0.861,1.000,0.698],
)
error_bars[observables_index("cFos_Protein")] = Dict(
    "EGF" => [0,0.036,0.028,0.056,0.071,0.048] ./ sqrt(3),
    "HRG" => [0,0.021,0.042,0.063,0.000,0.047] ./ sqrt(3),
)

experiments[observables_index("dusp_mRNA")] = Dict(
    "EGF" => [0.000,0.177,0.331,0.214,0.177,0.231],
    "HRG" => [0.000,0.221,0.750,1.000,0.960,0.934],
)
error_bars[observables_index("dusp_mRNA")] = Dict(
    "EGF" => [0.033,0.060,0.061,0.032,0.068,0.050] ./ sqrt(3),
    "HRG" => [0.027,0.059,0.094,0.124,0.113,0.108] ./ sqrt(3),
)

function get_timepoint(obs_name::String)::Vector{Float64}
    if obs_name in ["Phosphorylated_MEKc", "Phosphorylated_ERKc",
                    "Phosphorylated_RSKw", "Phosphorylated_cFos"]
        return [0., 300., 600., 900., 1800., 2700., 3600., 5400.]

    elseif obs_name == "Phosphorylated_CREBw"
        return [0., 600., 1800., 3600., 5400.]

    elseif obs_name == "cfos_mRNA"
        return [0.,600.,1200.,1800.,2700.,3600.,5400.]

    elseif obs_name in ["cFos_Protein", "dusp_mRNA"]
        return [0.,900.,1800.,2700.,3600.,5400.]
    end
end
end # module

max_stars_repo_name: himoto/paramEstim

max_issues_repo_name: himoto/paramEstim


text:
function loss(sim::SimResult, measurements::Vector)
  loss = 0.0
  for dp in measurements
    loss += loss_point(sim, dp)
  end
  return loss
end

function loss_point(sim::SimResult, dp::NormalMeasurementPoint{MU,S}) where {MU,S}
  sim_val = _param_value(sim, dp.μ, dp)
  measurements_val = dp.val
  sigma = _param_value(sim, dp.σ, dp)
  sigma_sq = (sigma)^2

  log(2π) + log(sigma_sq) + (sim_val - measurements_val)^2 / sigma_sq
end

function loss_point(sim::SimResult, dp::LogNormalMeasurementPoint{MU,S}) where {MU,S}
  sim_val = _param_value(sim, dp.μ, dp)
  measurements_val = dp.val
  sigma = _param_value(sim, dp.σ, dp)
  sigma_sq = (sigma)^2

  log(2π) + log(sigma_sq) + log((measurements_val)^2) + (log(sim_val) - log(measurements_val))^2 / sigma_sq
end

_param_value(sim, p::Float64, dp) = p

function _param_value(sim, p::Symbol, dp)
  if p ∈ keys(parameters(sim))
    val = parameters(sim)[p]
  elseif p ∈ observables(sim)
    val = sim(dp.t, p, dp.scope)
  else
    error("$p not found in simulated results.")
  end
  return val
end


text:
# Autogenerated wrapper script for SQLCipher_jll for armv7l-linux-gnueabihf
export libsqlcipher, sqlcipher

using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("SQLCipher")
JLLWrappers.@declare_library_product(libsqlcipher, "libsqlcipher.so.0")
JLLWrappers.@declare_executable_product(sqlcipher)
function __init__()
    JLLWrappers.@generate_init_header(OpenSSL_jll)
    JLLWrappers.@init_library_product(
        libsqlcipher,
        "lib/libsqlcipher.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        sqlcipher,
        "bin/sqlcipher",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using DiffEqBase: numargs
using StaticArrays: SVector, push
using ForwardDiff

using .Codim1: Codim1Problem

"""
Fixed point bifurcation problem.

See also: [`AbstractContinuationProblem`](@ref)

# Fields
* `homotopy::Function`:
  A function to compute ``H(x, t)`` where ``H`` is a homotopy
  ``H: \\mathbb R^N \\times \\mathbb R \\to \\mathbb R^N``.
  Function `homotopy` must be callable in one of the following form:
  `homotopy(x, p, t) ↦ H` (return `H`) for mutable state type or
  `homotopy(H, x, p, t)` (mutate `H`) for immutable state type.
* `homotopy_jacobian::Union{Function, Nothing}`:
  A function to compute ``H(x, t)`` and its Jacobian
  ``J = \\partial H / \\partial (x, t) \\in \\mathbb R^{N \\times (N+1)}``.
  Function `homotopy_jacobian` must be callable in one of the following form:
  `homotopy_jacobian(x, p, t) ↦ (H, J)` (return `(H, J)`) or
  `homotopy_jacobian(H, J, x, p, t)` (mutate `H` and `J`).
* `u0::Union{AbstractArray, Real}`: Initial state.
* `t0::Real`: Initial parameter.
* `t_domain::Tuple{<:Real, <:Real}`: Range of the parameter.
* `phase_space::Tuple{typeof(u0), typeof(u0)}`: A pair of lower and
  upper bound of the phase space.  Default is unbounded.
* `p`: Model parameter (constants).
"""
struct FixedPointBifurcationProblem{skind <: StateKind,
                                    tkind <: TimeKind,
                                    HJ, H, U, T, P,
                                    } <: Codim1Problem{skind, tkind}
    statekind::skind
    timekind::tkind
    homotopy_jacobian::HJ
    homotopy::H
    # TODO: rename u0 -> x0 (consistently use x for dynamical system state)
    u0::U
    t0::T
    t_domain::Tuple{T, T}
    phase_space::Tuple{U, U}
    p::P

    # TODO: Use Domains.jl?
end

"""
    unbounded_phase_space(u0::AbstractVector)
    unbounded_phase_space(u0::Real)
"""
function unbounded_phase_space(x0::X) where {X <: AbstractVector}
    return (X(fill(-Inf, length(x0))),
            X(fill(+Inf, length(x0))))
end
unbounded_phase_space(::X) where {X <: Real} = (X(-Inf), X(Inf))

function FixedPointBifurcationProblem(
        statekind::StateKind, timekind::TimeKind,
        homotopy::H, u0::U, t0::Real, t_domain::Tuple,
        p::P = nothing;
        homotopy_jacobian::HJ = nothing,
        phase_space = unbounded_phase_space(u0),
        ) where{HJ, H, U, P}
    T = promote_type(typeof(t0), map(typeof, t_domain)...)
    return FixedPointBifurcationProblem(
        statekind, timekind,
        homotopy_jacobian, homotopy,
        u0, t0, t_domain, phase_space, p)
end

struct HasJac end
struct NoJac end
hasjac(::FixedPointBifurcationProblem{<:Any, <:Any, <:Function}) = HasJac()
hasjac(::FixedPointBifurcationProblem{<:Any, <:Any, Nothing}) = NoJac()

function FixedPointBifurcationProblem(tkind::TimeKind,
                                      homotopy, args...; kwargs...)
    skind = numargs(homotopy) == 4 ? MutableState() : ImmutableState()
    return FixedPointBifurcationProblem(
        skind, tkind,
        homotopy, args...; kwargs...)
end

as_immutable_state(x::Tuple) = SVector(x)
as_immutable_state(x::Number) = x

function FixedPointBifurcationProblem(tkind::TimeKind,
                                      homotopy,
                                      u0::Union{Tuple, Number},
                                      args...; kwargs...)
    return FixedPointBifurcationProblem{false, tkind}(
        ImmutableState(), tkind,
        homotopy,
        as_immutable_state(u0),
        args...; kwargs...)
end

# TODO: trait
BifurcationsBase.contkind(::FixedPointBifurcationProblem) =
    BifurcationsBase.FixedPointCont()

struct FixedPointBifurcationCache{P, C, F} <: AbstractProblemCache{P}
    prob::P
    cfg::C
    residual::F
end

FixedPointBifurcationCache(prob::FixedPointBifurcationProblem) =
    _FixedPointBifurcationCache(statekind(prob), hasjac(prob), prob)

_FixedPointBifurcationCache(::Any, ::HasJac, prob) =
    FixedPointBifurcationCache(prob, nothing, nothing)

function _FixedPointBifurcationCache(::MutableState, ::NoJac, prob)
    x = get_u0(prob)
    y = similar(x, length(x) - 1)
    residual = (y, x) -> _residual!(y, x, prob)
    cfg = ForwardDiff.JacobianConfig(residual, y, x)
    return FixedPointBifurcationCache(prob, cfg, residual)
end

function _FixedPointBifurcationCache(::ImmutableState, ::NoJac, prob)
    x = get_u0(prob)
    residual = (x) -> _residual!(nothing, x, prob)
    cfg = ForwardDiff.JacobianConfig(residual, x)
    return FixedPointBifurcationCache(prob, cfg, residual)
end

TimeKind(::Type{<: FixedPointBifurcationCache{P}}) where P = TimeKind(P)

get_prob_cache(prob::FixedPointBifurcationProblem) =
    FixedPointBifurcationCache(prob)

get_u0(prob::FixedPointBifurcationProblem) = _get_u0(prob, prob.u0)

function _get_u0(prob::FixedPointBifurcationProblem, ::AbstractVector)
    u0 = similar(prob.u0, length(prob.u0) + 1)
    u0[1:end-1] = prob.u0
    u0[end] = prob.t0
    return u0
end

function _get_u0(prob::FixedPointBifurcationProblem, ::Real)
    return SVector(prob.u0, prob.t0)
end

function _get_u0(prob::FixedPointBifurcationProblem, ::SVector)
    return push(prob.u0, prob.t0)
end

residual!(H, u, cache::FixedPointBifurcationCache) =
    _residual!(H, u, cache.prob)

_residual!(H, u, prob) = _residual!(H, u, prob, statekind(prob), prob.u0)

residual_jacobian!(H, J, u, cache::FixedPointBifurcationCache) =
    _residual_jacobian!(H, J, u, cache,
                        statekind(cache.prob),
                        hasjac(cache.prob))

function isindomain(u, cache::FixedPointBifurcationCache)
    tmin, tmax = cache.prob.t_domain
    xmin, xmax = cache.prob.phase_space
    umin = vcat(xmin, SVector(tmin))
    umax = vcat(xmax, SVector(tmax))
    return all(umin .<= u .<= umax)
end

# ----------------------------------------------------------- inplace interface

function _residual!(H, u, prob::FixedPointBifurcationProblem,
                    ::MutableState, ::Any)
    x = u[1:end-1]
    t = u[end]
    prob.homotopy(H, x, prob.p, t)
    return H
end

function _residual_jacobian!(H, J, u, cache::FixedPointBifurcationCache,
                             ::MutableState, ::HasJac)
    prob = cache.prob
    x = u[1:end-1]
    t = u[end]
    prob.homotopy_jacobian(H, J, x, prob.p, t)
    return (H, J)
end

function _residual_jacobian!(H, J, u, cache::FixedPointBifurcationCache,
                             ::MutableState, ::NoJac)
    ForwardDiff.jacobian!(
        J,
        cache.residual,
        H,  # y
        u,  # x
        cache.cfg,
    )
    return (H, J)
end

# ------------------------------------------------------ out-of-place interface

ds_state(u::AbstractArray) = @view u[1:end-1]

@generated function ds_state(u::SVector{S, T}) where {S, T}
    values = [:(u[$i]) for i in 1:S-1]
    quote
        SVector{$(length(values)), $T}($(values...))
    end
end

function _residual!(::Any, u, prob::FixedPointBifurcationProblem,
                    ::ImmutableState, ::Real)
    x = u[1]
    t = u[end]
    return SVector(prob.homotopy(x, prob.p, t))
end

function _residual!(::Any, u, prob::FixedPointBifurcationProblem,
                    ::ImmutableState, ::AbstractArray)
    x = ds_state(u)
    t = u[end]
    return prob.homotopy(x, prob.p, t)
end

function _residual_jacobian!(_H, _J, u, cache::FixedPointBifurcationCache,
                             ::ImmutableState, ::HasJac)
    prob = cache.prob
    x = ds_state(u)
    t = u[end]
    return prob.homotopy_jacobian(x, prob.p, t)
end

function _residual_jacobian!(_H, _, u, cache::FixedPointBifurcationCache,
                             ::ImmutableState, ::NoJac)
    # TODO: Can I compute H and J in one go?  Or is it already
    # maximally efficient?
    H = cache.residual(u)
    J = ForwardDiff.jacobian(
        cache.residual,
        u,  # x
        cache.cfg,
    )
    return (H, J)
end


text:
using YaoBase
export KronBlock, kron

"""
    KronBlock{N, T, MT<:AbstractBlock} <: CompositeBlock{N, T}

composite block that combine blocks by kronecker product.
"""
struct KronBlock{N, MT<:AbstractBlock} <: CompositeBlock{N}
    slots::Vector{Int}
    locs::Vector{Int}
    blocks::Vector{MT}
end

KronBlock{N}(slots::Vector{Int}, locs::Vector{Int}, blocks::Vector{MT}) where {N, MT <: AbstractBlock} =
    KronBlock{N, MT}(slots, locs, blocks)

function KronBlock{N}(locs::Vector{Int}, blocks::Vector{MT}) where {N, MT<:AbstractBlock}
    perm = sortperm(locs)
    permute!(locs, perm)
    permute!(blocks, perm)
    @assert_locs_safe N collect(i:i+nqubits(b)-1 for (i, b) in zip(locs, blocks))

    slots = zeros(Int, N)
    for (i, each) in enumerate(locs)
        slots[each] = i
    end
    return KronBlock{N, MT}(slots, locs, blocks)
end

function KronBlock{N}(itr::Pair{Int,<:AbstractBlock}...) where N
    blocks = AbstractBlock[]
    locs = Int[]

    for (addr, block) in itr
        push!(locs, addr)
        push!(blocks, block)
    end
    return KronBlock{N}(locs, blocks)
end

function KronBlock(itr::AbstractBlock...)
    N = sum(nqubits, itr)
    locs = Int[]
    count = 1
    for each in itr
        push!(locs, count)
        count += nqubits(each)
    end
    return KronBlock{N}(locs, collect(itr))
end

KronBlock(blk::KronBlock) = copy(blk)

"""
    kron(n, blocks::Pair{Int, <:AbstractBlock}...)

Return a [`KronBlock`](@ref), with total number of qubits `n` and pairs of blocks.

# Example

Use `kron` to construct a `KronBlock`, it will put an `X` gate on the `1`st qubit,
and a `Y` gate on the `3`rd qubit.

```jldoctest
julia> kron(4, 1=>X, 3=>Y)
nqubits: 4
kron
├─ 1=>X gate
└─ 3=>Y gate

```
"""
Base.kron(total::Int, blocks::Pair{Int, <:AbstractBlock}...) = KronBlock{total}(blocks...)

"""
    kron(blocks::AbstractBlock...)
    kron(n, itr)

Return a [`KronBlock`](@ref), with total number of qubits `n`, and `blocks` should use all
the locations on `n` wires in quantum circuits.

# Example

You can use kronecker product to composite small blocks to a large blocks.

```jldoctest
julia> kron(X, Y, Z, Z)
nqubits: 4
kron
├─ 1=>X gate
├─ 2=>Y gate
├─ 3=>Z gate
└─ 4=>Z gate

```
"""
Base.kron(blocks::AbstractBlock...) = KronBlock(blocks...)
Base.kron(fs::Union{Function, AbstractBlock}...) = @λ(n->kron(n, fs...))

function Base.kron(total::Int, blocks::AbstractBlock...)
    sum(nqubits, blocks) == total || error("total number of qubits mismatch")
    return kron(blocks...)
end

Base.kron(total::Int, blocks::Union{AbstractBlock, Pair}...) =
    error("location of sparse distributed blocks must be explicit declared with pair (e.g 2=>X)")

Base.kron(total::Int, blocks::Base.Generator) = kron(total, blocks...)

"""
    kron(blocks...) -> f(n)
    kron(itr) -> f(n)

Return a lambda, which will take the total number of qubits as input.

# Example

If you don't know the number of qubit yet, or you are just too lazy, it is fine.

```jldoctest
julia> kron(put(1=>X) for _ in 1:2)
(n -> kron(n, (n  ->  put(n, 1 => X gate)), (n  ->  put(n, 1 => X gate))))

julia> kron(X for _ in 1:2)
nqubits: 2
kron
├─ 1=>X gate
└─ 2=>X gate

julia> kron(1=>X, 3=>Y)
(n -> kron(n, 1 => X gate, 3 => Y gate))
```
"""
Base.kron(blocks::Pair{Int, <:AbstractBlock}...,) = @λ(n->kron(n, blocks...))
Base.kron(blocks::Base.Generator) = kron(blocks...)

occupied_locs(k::KronBlock) = Iterators.flatten(map(x-> x + i - 1, occupied_locs(b)) for (i, b) in zip(k.locs, subblocks(k)))
subblocks(x::KronBlock) = x.blocks
chsubblocks(pb::KronBlock{N}, it) where N = KronBlock{N}(pb.locs, collect(it))
cache_key(x::KronBlock) = [cache_key(each) for each in x.blocks]
color(::Type{T}) where {T <: KronBlock} = :cyan

function _prepair_kronmat(k::KronBlock{N}) where N
    sizes = map(nqubits, subblocks(k))
    start_locs = @. N - $(k.locs) - sizes + 1

    order = sortperm(start_locs)
    sorted_start_locs = start_locs[order]
    num_bit_list = vcat(diff(push!(sorted_start_locs, N)) .- sizes[order])
    return order, num_bit_list, sorted_start_locs
end

function mat(::Type{T}, k::KronBlock{N}) where {T, N}
    order, num_bit_list, sorted_start_locs = _prepair_kronmat(k)
    blocks = subblocks(k)[order]
    return reduce(zip(blocks, num_bit_list), init=IMatrix{1 << sorted_start_locs[1], T}()) do x, y
        kron(x, mat(T, y[1]), IMatrix(1<<y[2]))
    end
end

function apply!(r::ArrayReg, k::KronBlock)
    _check_size(r, k)
    for (locs, block) in zip(k.locs, k.blocks)
        _instruct!(state(r), block, Tuple(locs:locs+nqubits(block)-1))
    end
    return r
end

_instruct!(state::AbstractArray{T}, block::AbstractBlock, locs) where T = instruct!(state, mat(T, block), locs)

# specialization
for G in [:X, :Y, :Z, :T, :S, :Sdag, :Tdag]
    GT = Expr(:(.), :ConstGate, QuoteNode(Symbol(G, :Gate)))
    @eval _instruct!(state::AbstractArray, block::$GT, locs) = instruct!(state, Val($(QuoteNode(G))), locs)
end

function Base.copy(k::KronBlock{N}) where N
    slots = copy(k.slots)
    locs = copy(k.locs)
    blocks = copy(k.blocks)
    return KronBlock{N}(slots, locs, blocks)
end

function Base.similar(k::KronBlock{N}) where N
    slots = zeros(Int, N)
    locs = empty!(similar(k.locs))
    blocks = empty!(similar(k.blocks))
    return KronBlock{N}(slots, locs, blocks)
end

function Base.getindex(k::KronBlock, addr)
    index = k.slots[addr]
    index == 0 && throw(KeyError(addr))
    return k.blocks[index]
end

function Base.setindex!(k::KronBlock, val, addr)
    index = k.slots[addr]
    index == 0 && return _insert_new!(k, val, addr)
    k.blocks[index] = val
    return k
end

function _insert_new!(k::KronBlock, val, addr)
    push!(k.locs, addr)
    push!(k.blocks, val)
    k.slots[addr] = lastindex(k.locs)
    return k
end

function Base.iterate(k::KronBlock, st = 1)
    if st > length(k)
        return nothing
    else
        return (k.locs[st], k.blocks[st]), st + 1
    end
end

Base.eltype(k::KronBlock) = Tuple{Int, AbstractBlock}
Base.length(k::KronBlock) = length(k.blocks)
Base.eachindex(k::KronBlock) = k.locs

function Base.:(==)(lhs::KronBlock{N}, rhs::KronBlock{N}) where N
    return all(lhs.locs .== rhs.locs) && all(lhs.blocks .== rhs.blocks)
end

Base.adjoint(blk::KronBlock{N}) where N = KronBlock{N}(blk.slots, blk.locs, map(adjoint, blk.blocks))

YaoBase.ishermitian(k::KronBlock) = all(ishermitian, k.blocks) || ishermitian(mat(k))
YaoBase.isunitary(k::KronBlock) = all(isunitary, k.blocks) || isunitary(mat(k))
YaoBase.isreflexive(k::KronBlock) = all(isreflexive, k.blocks) || isreflexive(mat(k))


text:
##
# Generic base overloads
Base.extrema(primitive::GeometryPrimitive) = (minimum(primitive), maximum(primitive))
function widths(x::AbstractRange)
    mini, maxi = Float32.(extrema(x))
    return maxi - mini
end

##
# conversion & decompose
convert_simplex(::Type{T}, x::T) where T = (x,)

function convert_simplex(NFT::Type{NgonFace{N, T1}}, f::Union{NgonFace{N, T2}}) where {T1, T2, N}
    return (convert(NFT, f),)
end

convert_simplex(NFT::Type{NgonFace{3,T}}, f::NgonFace{3,T2}) where {T, T2} = (convert(NFT, f),)

"""
    convert_simplex(::Type{Face{3}}, f::Face{N})

Triangulate an N-Face into a tuple of triangular faces.
"""
@generated function convert_simplex(::Type{TriangleFace{T}}, f::Union{SimplexFace{N}, NgonFace{N}}) where {T, N}
    3 <= N || error("decompose not implemented for N <= 3 yet. N: $N")# other wise degenerate
    v = Expr(:tuple)
    for i = 3:N
        push!(v.args, :(TriangleFace{T}(f[1], f[$(i-1)], f[$i])))
    end
    return v
end

"""
    convert_simplex(::Type{Face{2}}, f::Face{N})

Extract all line segments in a Face.
"""
@generated function convert_simplex(::Type{LineFace{T}}, f::Union{SimplexFace{N}, NgonFace{N}}) where {T, N}
    2 <= N || error("decompose not implented for N <= 2 yet. N: $N")# other wise degenerate

    v = Expr(:tuple)
    for i = 1:N-1
        push!(v.args, :(LineFace{$T}(f[$i], f[$(i+1)])))
    end
    # connect vertices N and 1
    push!(v.args, :(LineFace{$T}(f[$N], f[1])))
    return v
end

to_pointn(::Type{T}, x) where T<:Point = convert_simplex(T, x)[1]

# disambiguation method overlords
convert_simplex(::Type{Point}, x::Point) = (x,)
convert_simplex(::Type{Point{N,T}}, p::Point{N,T}) where {N, T} = (p,)
function convert_simplex(::Type{Point{N, T}}, x) where {N, T}
    N2 = length(x)
    return (Point{N, T}(ntuple(i-> i <= N2 ? T(x[i]) : T(0), N)),)
end

function convert_simplex(::Type{Vec{N, T}}, x) where {N, T}
    N2 = length(x)
    return (Vec{N, T}(ntuple(i-> i <= N2 ? T(x[i]) : T(0), N)),)
end

collect_with_eltype(::Type{T}, vec::Vector{T}) where T = vec
collect_with_eltype(::Type{T}, vec::AbstractVector{T}) where T = collect(vec)

function collect_with_eltype(::Type{T}, iter) where T
    # TODO we could be super smart about allocating the right length
    # but its kinda annoying, since e.g. T == Triangle and first(iter) isa Quad
    # will need double the length etc - but could all be figured out ;)
    result = T[]
    for element in iter
        # convert_simplex always returns a tuple,
        # so that e.g. convert(Triangle, quad) can return 2 elements
        for telement in convert_simplex(T, element)
            push!(result, telement)
        end
    end
    return result
end

"""
The unnormalized normal of three vertices.
"""
function orthogonal_vector(v1, v2, v3)
    a = v2 - v1
    b = v3 - v1
    return cross(a, b)
end

"""
```
normals{VT,FD,FT,FO}(vertices::Vector{Point{3, VT}},
                    faces::Vector{Face{FD,FT,FO}},
                    NT = Normal{3, VT})
```
Compute all vertex normals.
"""
function normals(vertices::AbstractVector{<: AbstractPoint{3, T}},
                 faces::AbstractVector{F};
                 normaltype=Vec{3, T}) where {T, F <: NgonFace}
    normals_result = zeros(normaltype, length(vertices)) # initilize with same type as verts but with 0
    for face in faces
        v = metafree.(vertices[face])
        # we can get away with two edges since faces are planar.
        n = orthogonal_vector(v[1], v[2], v[3])
        for i =1:length(F)
            fi = face[i]
            normals_result[fi] = normals_result[fi] + n
        end
    end
    normals_result .= normalize.(normals_result)
    return normals_result
end

##
# Some more primitive types

"""
    HyperSphere{N, T}

A `HyperSphere` is a generalization of a sphere into N-dimensions.
A `center` and radius, `r`, must be specified.
"""
struct HyperSphere{N, T} <: GeometryPrimitive{N, T}
    center::Point{N, T}
    r::T
end
"""
    Circle{T}

An alias for a HyperSphere of dimension 2. (i.e. `HyperSphere{2, T}`)
"""
const Circle{T} = HyperSphere{2, T}

"""
    Sphere{T}

An alias for a HyperSphere of dimension 3. (i.e. `HyperSphere{3, T}`)
"""
const Sphere{T} = HyperSphere{3, T}

"""
    Quad{T}

A rectangle in 3D space.
"""
struct Quad{T} <: GeometryPrimitive{3, T}
    downleft::Vec{3, T}
    width   ::Vec{3, T}
    height  ::Vec{3, T}
end

struct Pyramid{T} <: GeometryPrimitive{3, T}
    middle::Point{3, T}
    length::T
    width ::T
end

struct Particle{N, T} <: GeometryPrimitive{N, T}
    position::Point{N, T}
    velocity::Vec{N, T}
end

"""
    Cylinder{N, T}

A `Cylinder` is a 2D rectangle or a 3D cylinder defined by its origin point,
its extremity and a radius. `origin`, `extremity` and `r`, must be specified.
"""
struct Cylinder{N, T} <: GeometryPrimitive{N, T}
    origin::Point{N,T}
    extremity::Point{N,T}
    r::T
end

"""
    Cylinder2{T}
    Cylinder3{T}

A `Cylinder2` or `Cylinder3` is a 2D/3D cylinder defined by its origin point,
its extremity and a radius. `origin`, `extremity` and `r`, must be specified.
"""
const Cylinder2{T} = Cylinder{2, T}
const Cylinder3{T} = Cylinder{3, T}

origin(c::Cylinder{N, T}) where {N, T} = c.origin
extremity(c::Cylinder{N, T}) where {N, T} = c.extremity
radius(c::Cylinder{N, T}) where {N, T} = c.r
height(c::Cylinder{N, T}) where {N, T} = norm(c.extremity - c.origin)
direction(c::Cylinder{N, T}) where {N, T} = (c.extremity .- c.origin) ./ height(c)

function rotation(c::Cylinder{2, T}) where T
    d2 = direction(c); u = @SVector [d2[1], d2[2], T(0)]
    v = @MVector [u[2], -u[1], T(0)]
    normalize!(v)
    return hcat(v, u, @SVector T[0, 0, 1])
end

function rotation(c::Cylinder{3, T}) where T
    d3 = direction(c); u = @SVector [d3[1], d3[2], d3[3]]
    if abs(u[1]) > 0 || abs(u[2]) > 0
        v = @MVector [u[2], -u[1], T(0)]
    else
        v = @MVector [T(0), -u[3], u[2]]
    end
    normalize!(v)
    w = @SVector [u[2] * v[3] - u[3] * v[2], -u[1] * v[3] + u[3] * v[1], u[1] * v[2] - u[2] * v[1]]
    return hcat(v, w, u)
end

function coordinates(c::Cylinder{2, T}, nvertices=(2, 2)) where T
    r = Rect(c.origin[1] - c.r/2, c.origin[2], c.r, height(c))
    M = rotation(c)
    points = coordinates(r, nvertices)
    vo = to_pointn(Point3{T}, origin(c))
    return (M * (to_pointn(Point3{T}, point) .- vo) .+ vo for point in points)
end

function faces(sphere::Cylinder{2}, nvertices=(2, 2))
    return faces(Rect(0, 0, 1, 1), nvertices)
end

function coordinates(c::Cylinder{3, T}, nvertices=30) where T
    if isodd(nvertices)
        nvertices = 2 * (nvertices ÷ 2)
    end
    nvertices = max(8, nvertices);
    nbv = nvertices ÷ 2

    M = rotation(c)
    h = height(c)
    range = 1:(2 * nbv + 2)
    function inner(i)
        if i == length(range)
            return c.extremity
        elseif i == length(range) - 1
            return origin(c)
        else
            phi = T((2π * (((i + 1) ÷ 2) - 1)) / nbv)
            up = ifelse(isodd(i), 0, h)
            return (M * Point(c.r * cos(phi), c.r * sin(phi), up)) .+ c.origin
        end
    end

    return (inner(i) for i in range)
end

function faces(c::Cylinder{3}, facets=30)
    isodd(facets) ? facets = 2 * div(facets, 2) : nothing
    facets < 8 ? facets = 8 : nothing; nbv = Int(facets / 2)
    indexes = Vector{TriangleFace{Int}}(undef, facets)
    index = 1
    for j = 1:(nbv-1)
        indexes[index] = (index + 2, index + 1, index)
        indexes[index + 1] = ( index + 3, index + 1, index + 2)
        index += 2
    end
    indexes[index] = (1, index + 1, index)
    indexes[index + 1] = (2, index + 1, 1)

    for i = 1:length(indexes)
        i%2 == 1 ? push!(indexes, (indexes[i][1], indexes[i][3], 2*nbv+1)) : push!(indexes,(indexes[i][2], indexes[i][1], 2*nbv+2))
    end
    return indexes
end

##
# Sphere

HyperSphere{N}(p::Point{N, T}, number) where {N, T} = HyperSphere{N, T}(p, convert(T, number))

widths(c::HyperSphere{N, T}) where {N, T} = Vec{N, T}(radius(c)*2)
radius(c::HyperSphere) = c.r
origin(c::HyperSphere) = c.center

Base.minimum(c::HyperSphere{N, T}) where {N, T} = Vec{N, T}(origin(c)) - Vec{N, T}(radius(c))
Base.maximum(c::HyperSphere{N, T}) where {N, T} = Vec{N, T}(origin(c)) + Vec{N, T}(radius(c))

function Base.in(x::AbstractPoint{2}, c::Circle)
    @inbounds ox, oy = origin(c)
    xD = abs(ox - x)
    yD = abs(oy - y)
    return xD <= c.r && yD <= c.r
end

centered(S::Type{HyperSphere{N, T}}) where {N, T} = S(Vec{N,T}(0), T(0.5))
centered(::Type{T}) where {T <: HyperSphere} = centered(HyperSphere{ndims_or(T, 3), eltype_or(T, Float32)})

function coordinates(s::Circle, nvertices=64)
    rad = radius(s)
    inner(fi) = Point(rad*sin(fi + pi), rad*cos(fi + pi)) .+ origin(s)
    return (inner(fi) for fi in LinRange(0, 2pi, nvertices))
end

function texturecoordinates(s::Circle, nvertices=64)
    return coordinates(Circle(Point2f0(0.5), 0.5f0), nvertices)
end

function coordinates(s::Sphere, nvertices=24)
    θ = LinRange(0, pi, nvertices); φ = LinRange(0, 2pi, nvertices)
    inner(θ, φ) = Point(cos(φ)*sin(θ), sin(φ)*sin(θ), cos(θ)) .* s.r .+ s.center
    return ivec((inner(θ, φ) for θ in θ, φ in φ))
end

function texturecoordinates(s::Sphere, nvertices=24)
    ux = LinRange(0, 1, nvertices)
    return ivec(((φ, θ) for θ in reverse(ux), φ in ux))
end

function faces(sphere::Sphere, nvertices=24)
    return faces(Rect(0, 0, 1, 1), (nvertices, nvertices))
end

function normals(s::Sphere{T}, nvertices=24) where {T}
    return coordinates(Sphere(Point{3, T}(0), 1), nvertices)
end

function coordinates(p::Pyramid{T}, nvertices=nothing) where {T}
    leftup = Point{3, T}(-p.width , p.width, 0) / 2
    leftdown = Point(-p.width, -p.width, 0) / 2
    tip = Point{3, T}(p.middle + Point{3, T}(0, 0, p.length))
    lu = Point{3, T}(p.middle + leftup)
    ld = Point{3, T}(p.middle + leftdown)
    ru = Point{3, T}(p.middle - leftdown)
    rd = Point{3, T}(p.middle - leftup)
    return Point{3, T}[
        tip, rd, ru,
        tip, ru, lu,
        tip, lu, ld,
        tip, ld, rd,
        rd,  ru, lu,
        lu,  ld, rd
    ]
end

function faces(r::Pyramid, nvertices=nothing) where FT
    return (TriangleFace(triangle) for triangle in TupleView{3}(1:18))
end


text:
export pai

"""
    pai(x, y, d, τ; w = 0, correspondence_measure = Statistics.cor) → Float64
    pai(x, y, d, τ, bootstrap_method::Symbol; w = 0, correspondence_measure = Statistics.cor,
        method = :segment, L = ceil(Int, (length(x)-d*τ)*0.2), nreps = 100) → Vector{Float64}

Compute the pairwise asymmetric inference (PAI; McCracken, 2014)[^McCracken2014] between `x` and `y`.
Returns the correspondence between original and cross mapped values (the default is
`ρ = correspondence_measure(y(t), ỹ(t) | M_xy)`).

PAI is a modification to Sugihara et al. (2012)'s CCM algorithm[^Sugihara2012], where instead of
using completely out-of-sample prediction when trying to predict ``y(t)``, values about *both* variables
are included in the embedding used to make predictions. Specifically, PAI computes the
correspondence between the values ``y(t)`` and the cross-map estimated values ``ỹ(t) | M_xy``,
where the ``\\tilde{y}(t)`` are the values estimated using the embedding ``M_{xy} = \\{ ( x_t, x_{t-\\tau}, x_{t-2\\tau}, \\ldots, x_{t-(d - 1)\\tau} ) \\}``.
*Note: a `d+1`-dimensional embedding is used, rather than the `d`-dimensional embedding used for CCM.

Like for the CCM algorithm, the Theiler window `r` indicates how many temporal neighbors of the predictee is to be excluded
during the nearest neighbors search (the default `r = 0` excludes only the predictee itself, while
`r = 2` excludes the point itself plus its two nearest neighbors in time).

If `bootstrap_method` is specified, then `nreps` different bootstrapped estimates of
`correspondence_measure(y(t), ỹ(t) | M_x)` are returned. The following bootstrap methods are available:

- `bootstrap_method = :random` selects training sets of length `L` consisting of randomly selected
    points from the embedding ``M_x``  (time ordering does not matter). This is method 3 from Luo
    et al. (2015)[^Luo2015], which critiqued the original Sugihara et al. methodology.
- `bootstrap_method = :segment` selects training sets consisting of time-contiguous segments
    (each of lenght `L`) of embedding vectors in ``M_x`` (time ordering matters). This is
    method 2 from Luo et al. (2015)[^Luo2015].

[^McCracken2014]: McCracken, James M., and Robert S. Weigel. "Convergent cross-mapping and pairwise asymmetric inference." Physical Review E 90.6 (2014): 062903.
[^Sugihara2012]: Sugihara, George, et al. "Detecting causality in complex ecosystems." Science (2012): 1227079.[http://science.sciencemag.org/content/early/2012/09/19/science.1227079](http://science.sciencemag.org/content/early/2012/09/19/science.1227079)
[^Luo2015]: "Questionable causality: Cosmic rays to temperature." Proceedings of the National Academy of Sciences Aug 2015, 112 (34) E4638-E4639; DOI: 10.1073/pnas.1510571112 Ming Luo, Holger Kantz, Ngar-Cheung Lau, Wenwen Huang, Yu Zhou
"""
function pai(x, y, d, τ; correspondence_measure = Statistics.cor, r = 0)
    # Embed `x` and `y`, using only current values for `y` and current+ historical values for `x`.
    Mₓy = crossmapembed(x, y, d, τ, PAIEmbedding())

    # Compute `ỹ | Mₓy` and return `correspondence_measure(ỹ | Mₓy, y).
    crossmap_basic(Mₓy, y, d, τ; correspondence_measure = correspondence_measure, r = r)
end

function pai(x, y, d, τ, bootstrap_method::Symbol;
        L = ceil(Int, (length(x) - d * τ) * 0.2), nreps = 100,
        r = 0, correspondence_measure = Statistics.cor)
    # Embed `x` and `y`, using only current values for `y` and current+ historical values for `x`.
    @assert length(x) == length(y)
    @assert length(x) - d * τ > 0
    Mₓy = crossmapembed(x, y, d, τ, PAIEmbedding())

    # Compute `ỹ | Mₓy` and return `nreps` values for `correspondence_measure(ỹ | Mₓy, y),
    # each value being computed on a random library (sequential or not) consisting of `L` points.
    return crossmap_bootstrap(Mₓy, y, d, τ, bootstrap_method;
        correspondence_measure = correspondence_measure, r = r)
end


text:
# # Sequential Simulations with [PowerSimulations.jl](https://github.com/NREL-SIIP/PowerSimulations.jl)

# **Originally Contributed by**: Clayton Barrows

# ## Introduction

# PowerSimulations.jl supports simulations that consist of sequential optimization problems
# where results from previous problems inform subsequent problems in a variety of ways. This
# example demonstrates some of these capabilities to represent electricity market clearing.

# ## Dependencies
# Since the `OperatiotnsProblem` is the fundamental building block of a sequential
# simulation in PowerSimulations, we will build on the [OperationsProblem example](https://nbviewer.jupyter.org/github/NREL-SIIP/SIIPExamples.jl/blob/master/notebook/3_PowerSimulations_examples/01_operations_problems.ipynb)
# by sourcing it as a dependency.
using SIIPExamples
pkgpath = dirname(dirname(pathof(SIIPExamples)))
include(
    joinpath(pkgpath, "test", "3_PowerSimulations_examples", "01_operations_problems.jl"),
)

# ### 5-Minute system
# We had already created a `sys::System` from hourly RTS data in the OperationsProblem example.
# The RTS data also includes 5-minute resolution time series data. So, we can create another
# `System`:
sys_RT = build_system(PSITestSystems, "modified_RTS_GMLC_RT_sys")

# ## `OperationsProblemTemplate`s define `Stage`s
# Sequential simulations in PowerSimulations are created by defining `OperationsProblems`
# that represent `Stages`, and how information flows between executions of a `Stage` and
# between different `Stage`s.
#
# Let's start by defining a two stage simulation that might look like a typical day-Ahead
# and real-time electricity market clearing process.

# ### We've already defined the reference model for the day-ahead unit commitment
#set_device_model!(template_ed, GenericBattery, BookKeeping)
template_uc

# ### Define the reference model for the real-time economic dispatch
# In addition to the manual specification process demonstrated in the OperationsProblem
# example, PSI also provides pre-specified templates for some standard problems:
template_ed = template_economic_dispatch()

# ### Define the `SimulaitonProblems`
# `OperationsProblem`s define models. The actual problem will change as the stage gets updated to represent
# different time periods, but the formulations applied to the components is constant within
# a stage. In this case, we want to define two stages with the `OperationsProblemTemplate`s
# and the `System`s that we've already created.
problems = SimulationProblems(
    UC = OperationsProblem(template_uc, sys, optimizer = solver),
    ED = OperationsProblem(
        template_ed,
        sys_RT,
        optimizer = solver,
        balance_slack_variables = true,
    ),
)
# Note that the "ED" problem has a `balance_slack_variables = true` argument. This adds slack
# variables with a default penalty of 1e6 to the nodal energy balance constraint and helps
# ensure feasibility with some performance impacts.

# ### `SimulationSequence`
# Similar to an `OperationsProblemTemplate`, the `SimulationSequence` provides a template of
# how to execute a sequential set of operations problems.

#nb # print_struct(SimulationSequence)

# Let's review some of the `SimulationSequence` arguments.

# ### Chronologies
# In PowerSimulations, chronologies define where information is flowing. There are two types
# of chronologies.
#  - inter-stage chronologies: Define how information flows between stages. e.g. day-ahead
# solutions are used to inform economic dispatch problems
#  - intra-stage chronologies: Define how information flows between multiple executions of a
# single stage. e.g. the dispatch setpoints of the first period of an economic dispatch problem
# are constrained by the ramping limits from setpoints in the final period of the previous problem.
#

# Let's define an inter-stage chronology that synchronizes information from 24 periods of
# the first stage with a set of executions of the second stage:

feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24))

# ### `FeedForward` and `Cache`
# The definition of exactly what information is passed using the defined chronologies is
# accomplished with `FeedForward` and `Cache` objects. Specifically, `FeedForward` is used
# to define what to do with information being passed with an inter-stage chronology. Let's
# define a `FeedForward` that affects the semi-continuous range constraints of thermal generators
# in the economic dispatch problems based on the value of the unit-commitment variables.

feedforward = Dict(
    ("ED", :devices, :ThermalStandard) => SemiContinuousFF(
        binary_source_problem = PSI.ON,
        affected_variables = [PSI.ACTIVE_POWER],
    ),
)

# ### Sequencing
# The stage problem length, look-ahead, and other details surrounding the temporal Sequencing
# of stages are controlled using the `intervals` argument and the structure of the `Forecast`
# data in the `System` of each problem.
#  - intervals::Dict(String, Dates.Period) : defines the interval with which stage problems
# advance after each execution. e.g. day-ahead problems have an interval of 24-hours
#
# So, to define a typical day-ahead - real-time sequence, we can define the following:
#  - Day ahead problems should represent 48 hours, advancing 24 hours after each execution (24-hour look-ahead)
#  - Real time problems should represent 1 hour (12 5-minute periods), advancing 15 min after each execution (15 min look-ahead)

intervals = Dict("UC" => (Hour(24), Consecutive()), "ED" => (Minute(15), Consecutive()))

# Finally, we can put it all together:

DA_RT_sequence = SimulationSequence(
    problems = problems,
    intervals = intervals,
    ini_cond_chronology = InterProblemChronology(),
    feedforward_chronologies = feedforward_chronologies,
    feedforward = feedforward,
)

# ## `Simulation`
# Now, we can build and execute a simulation using the `SimulationSequence` and `Stage`s
# that we've defined.
sim = Simulation(
    name = "rts-test",
    steps = 1,
    problems = problems,
    sequence = DA_RT_sequence,
    simulation_folder = dirname(dirname(pathof(SIIPExamples))),
)

# ### Build simulation

build!(sim)

# ### Execute simulation
# the following command returns the status of the simulation (0: is proper execution) and
# stores the results in a set of HDF5 files on disk.
execute!(sim, enable_progress_bar = false)

# ## Results
# To access the results, we need to load the simulation result metadata and then make
# requests to the specific data of interest. This allows you to efficiently access the
# results of interest without overloading resources.
results = SimulationResults(sim);
uc_results = get_problem_results(results, "UC"); # UC stage result metadata
ed_results = get_problem_results(results, "ED"); # ED stage result metadata

# Now we can read the specific results of interest for a specific problem, time window (optional),
# and set of variables, duals, or parameters (optional)

read_variables(uc_results, names = [:P__ThermalStandard, :P__RenewableDispatch])

# Or if we want the result of just one variable, parameter, or dual (must be defined in the
# problem definition), we can use:

read_parameter(
    ed_results,
    :P__max_active_power__RenewableFix,
    initial_time = DateTime("2020-01-01T06:00:00"),
    count = 5,
)

# * note that this returns the results of each execution step in a separate dataframe *
# If you want the realized results (without lookahead periods), you can call `read_realized_*`:

read_realized_variables(uc_results, names = [:P__ThermalStandard, :P__RenewableDispatch])

# ## Plotting
# Take a look at the [plotting examples.](https://nbviewer.jupyter.org/github/NREL-SIIP/SIIPExamples.jl/blob/master/notebook/3_PowerSimulations_examples/04_bar_stack_plots.ipynb)


text:
# ---
# title: 1078. Occurrences After Bigram
# id: problem1078
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Hash Table
# link: <https://leetcode.com/problems/occurrences-after-bigram/description/>
# hidden: true
# ---
#
# Given words `first` and `second`, consider occurrences in some `text` of the
# form "`first second third`", where `second` comes immediately after `first`,
# and `third` comes immediately after `second`.
#
# For each such occurrence, add "`third`" to the answer, and return the answer.
#
#
#
# **Example 1:**
#
#
#
#     Input: text = "alice is a good girl she is a good student", first = "a", second = "good"
#     Output: ["girl","student"]
#
#
# **Example 2:**
#
#
#
#     Input: text = "we will we will rock you", first = "we", second = "will"
#     Output: ["we","rock"]
#
#
#
#
# **Note:**
#
#   1. `1 <= text.length <= 1000`
#   2. `text` consists of space separated words, where each word consists of lowercase English letters.
#   3. `1 <= first.length, second.length <= 10`
#   4. `first` and `second` consist of lowercase English letters.
#
#
## @lc code=start
using LeetCode

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


text:
immutable TDist <: ContinuousUnivariateDistribution
    df::Float64 # non-integer degrees of freedom allowed
    function TDist(d::Real)
    	d > zero(d) || error("df must be positive")
        new(float64(d))
    end
end

@_jl_dist_1p TDist t

function entropy(d::TDist)
    ((d.df + 1.0) / 2.0) *
        (digamma((d.df + 1.0) / 2.0) - digamma((d.df) / 2.0)) +
        (1.0 / 2.0) * log(d.df) + lbeta(d.df + 1.0, 1.0 / 2.0)
end

insupport(::TDist, x::Real) = isfinite(x)
insupport(::Type{TDist}, x::Real) = isfinite(x)

mean(d::TDist) = d.df > 1.0 ? 0.0 : NaN

median(d::TDist) = 0.0

mode(d::TDist) = 0.0
modes(d::TDist) = [0.0]

function pdf(d::TDist, x::Real)
    1.0 / (sqrt(d.df) * beta(0.5, 0.5 * d.df)) *
        (1.0 + x^2 / d.df)^(-0.5 * (d.df + 1.0))
end

function var(d::TDist)
    d.df > 2.0 && return d.df / (d.df - 2.0)
    d.df > 1.0 && return Inf
    NaN
end


text:
# # Robust approximate fitting
# Section 6.4.2
# Boyd & Vandenberghe "Convex Optimization"
# Original by Lieven Vandenberghe
# Adapted for Convex by Joelle Skaf - 10/03/05
#
# Adapted for Convex.jl by Karanveer Mohan and David Zeng - 26/05/14
# Original cvx code and plots here:
# <http://web.cvxr.com/cvx/examples/cvxbook/Ch06_approx_fitting/html/fig6_15.html>
#
# Consider the least-squares problem:
#       minimize $\|(A + tB)x - b\|_2$
# where $t$ is an uncertain parameter in [-1,1]
# Three approximate solutions are found:
#
#   1. nominal optimal (i.e. letting t=0)
#   2. stochastic robust approximation:
#           minimize $\mathbb{E}\|(A+tB)x - b\|_2$
#      assuming $u$ is uniformly distributed on [-1,1].
#      (reduces to minimizing $\mathbb{E} \|(A+tB)x-b\|^2 = \|A*x-b\|^2  + x^TPx$
#        where $P = \mathbb{E}(t^2) B^TB = (1/3) B^TB$ )
#   3. worst-case robust approximation:
#           minimize $\mathrm{sup}_{-1\leq u\leq 1} \|(A+tB)x - b\|_2$
#      (reduces to minimizing $\max\{\|(A-B)x - b\|_2, \|(A+B)x - b\|_2\}$ ).
#
using Convex, LinearAlgebra, SCS
if VERSION < v"1.2.0-DEV.0"
     LinearAlgebra.diagm(v::AbstractVector) = diagm(0 => v)
end
# Input Data
m = 20;
n = 10;
A = randn(m,n);
(U,S,V) = svd(A);
S = diagm(exp10.(range(-1, stop=1, length=n)));
A = U[:, 1:n] * S * V';

B = randn(m, n);
B = B / norm(B);

b = randn(m, 1);
x = Variable(n)

# Case 1: Nominal optimal solution
p = minimize(norm(A * x - b, 2))
solve!(p, () -> SCS.Optimizer(verbose=0))
x_nom = evaluate(x)

# Case 2: Stochastic robust approximation
P = 1 / 3 * B' * B;
p = minimize(square(pos(norm(A * x - b))) + quadform(x, Symmetric(P)))
solve!(p, () -> SCS.Optimizer(verbose=0))
x_stoch = evaluate(x)

# Case 3: Worst-case robust approximation
p = minimize(max(norm((A - B) * x - b), norm((A + B) * x - b)))
solve!(p, () -> SCS.Optimizer(verbose=0))
x_wc = evaluate(x)

# Plot residuals:
parvals = range(-2, stop=2, length=100);

errvals(x) = [ norm((A + parvals[k] * B) * x - b) for k = eachindex(parvals)]
errvals_ls = errvals(x_nom)
errvals_stoch = errvals(x_stoch)
errvals_wc = errvals(x_wc)

using Plots
plot(parvals, errvals_ls, label="Nominal problem")
plot!(parvals, errvals_stoch, label="Stochastic Robust Approximation")
plot!(parvals, errvals_wc, label="Worst-Case Robust Approximation")
plot!(title="Residual r(u) vs a parameter u for three approximate solutions", xlabel="u", ylabel="r(u) = ||A(u)x-b||_2")


text:
function test_is_same_object(this_object, that_object; except=[])
  cleaned_except = Symbol.(except)

  @test typeof(this_object) == typeof(this_object)
  all_same = _test_is_same_object(this_object, that_object, cleaned_except)

  @test all_same
end

function test_is_not_same_object(this_object, that_object; except=[])
  cleaned_except = Symbol.(except)

  is_same_object = typeof(this_object) == typeof(this_object)

  if is_same_object
    is_same_object = _test_is_same_object(this_object, that_object, cleaned_except, true)
  end

  @test !is_same_object
end

function _test_is_same_object(this_object, that_object, cleaned_except, is_silent=false)
  all_same = true

  for cur_field in fieldnames(typeof(this_object))
    ( cur_field in cleaned_except ) && continue

    this_value, that_value = getfield.([this_object, that_object], cur_field)
    ( this_value == that_value ) && continue

    if isa(this_value, AbstractDict) && isa(that_value, AbstractDict)
      this_dict, that_dict = map(cur_flat -> cur_flat.it, Base.Iterators.flatten.([this_value, that_value]))

      if sort(collect(keys(this_dict))) == sort(collect(keys(that_dict)))
        cur_keys = sort(collect(keys(this_dict)))

        actually_good = true
        for cur_key in cur_keys
          ( Symbol(cur_key) in cleaned_except ) && continue
          ( this_dict[cur_key] == that_dict[cur_key] ) && continue

          actually_good = false
          println([cur_key, this_dict[cur_key], that_dict[cur_key]])
          break
        end

        actually_good && continue
      else
        println([cur_field, this_value, that_value])
      end
    else
      println([cur_field, this_value, that_value])
    end

    all_same = false
    break
  end

  all_same
end