text:
function computus(year; servois=false)

    # Year's position on the 19 year metonic cycle
    a = mod(year, 19)

    # Century index
    k = floor(Int, year/100)

    # Shift of metonic cycle, add a day offset every 300 years
    p = floor(Int, (13 + 8 * k) / 25)

    # Correction for non-observed leap days
    q = floor(Int, k / 4)

    # Correction to starting point of calculation each century
    M = mod(15 - p + k - q, 30)

    # Number of days from March 21st until the full moon
    d = mod(19 * a + M, 30)

    # Returning if user wants value for Servois' table
    if servois
        return mod(21+d,31)
    end

    # Finding the next Sunday
    # Century-based offset in weekly calculation
    N = mod(4 + k - q, 7)

    # Correction for leap days
    b = mod(year, 4)
    c = mod(year, 7)

    # Days from d to next Sunday
    e = mod(2 * b + 4 * c + 6 * d + N, 7)

    # Historical corrections for April 26 and 25
    if (d == 29 && e == 6) || (d == 28 && e == 6 && a > 10)
        e = -1
    end

    # Determination of the correct month for Easter
    if(22+ d + e > 31)
        return "April "*string(d + e - 9)
    else
        return "March "*string(22+d+e)
    end
end

# Here, we will output the date of the Paschal full moon
# (using Servois notation), and Easter for 2020-2030

a = [i for i = 2020:2030]
servois_numbers = computus.(a;servois=true)
easter_dates = computus.(a)

println("The following are the dates of the Paschal full moon (using Servois "*
        "notation) and the date of Easter for 2020-2030 AD:")
println("Year	Servois number	Easter")
for i = 1:length(a)
    println(string(a[i])*"\t"*string(servois_numbers[i])*"\t\t"*
            string(easter_dates[i]))
end

max_stars_repo_name: mfarberbrodsky/algorithm-archive

max_issues_repo_name: mfarberbrodsky/algorithm-archive

max_forks_repo_name: mfarberbrodsky/algorithm-archive


text:
#export: degree_relative, random_elem, one_root, norm_equation

degree_relative(L::Hecke.LocalField, K::Union{FlintQadicField, Hecke.LocalField}) = divexact(absolute_degree(L), absolute_degree(K))

degree_relative(L::Hecke.RelFinField, K::Union{FqNmodFiniteField, Hecke.RelFinField}) = divexact(absolute_degree(L), absolute_degree(K))

degree_relative(L::FqNmodFiniteField, K::Hecke.Nemo.GaloisField) = divexact(absolute_degree(L), absolute_degree(K))

##############################################
#random element with small coefficients
##############################################

function random_elem(L::Union{FlintQadicField, Hecke.LocalField})
   b = basis(L)
   n = degree(L)
   r = [rand(1:5*n) for i in 1:n]   # Choose small coefficients
   return sum( [r[i]*b[i] for i in 1:n])
end

################ one_root computes a single root in the finite field extensions#############################
function one_root(f::Union{gfp_poly, fq_nmod_poly}, F::Union{FqNmodFiniteField, Hecke.RelFinField})
   g = polynomial(F, [coeff(f,i) for i = 0:degree(f) ] )
   fac = factor(g)
   if length(fac) == 1
      error("no roots")
   end
   r = first(fac)[1]
   @assert degree(r) == 1
   return -coeff(r,0)
end

function roots(f::Union{gfp_poly, fq_nmod_poly}, F::Union{FqNmodFiniteField, Hecke.RelFinField})
   g = polynomial(F, [coeff(f,i) for i = 0:degree(f) ] )
   fac = factor(g)
   return [-coeff(r,0) for (r, _) in fac if degree(r) == 1]
end

function one_root(f::Hecke.AbstractAlgebra.Generic.Poly, F::Hecke.RelFinField)
   g = polynomial(F, [coeff(f,i) for i = 0:degree(f)])
   fac = factor(g)
   if length(fac) == 1
      error("no roots")
   end
   r = first(fac)[1]
   @assert degree(r) == 1
   return -coeff(r,0)
end

######################### norm equation over finite fields ##############

function norm_equation(F::Union{FqNmodFiniteField, Hecke.RelFinField}, b::Union{gfp_elem, fq_nmod})
   if iszero(b)
      return zero(F)
   end
   k = parent(b)
   n = degree_relative(F,k)
   f = polynomial(k,vcat([b],[rand(k) for i = 1:n-1],[1]))
   while !is_irreducible(f)
      f = polynomial(k,vcat([b],[rand(k) for i = 1:n-1],[1]))
   end
   return (-1)^(n)*one_root(f,F)
end

################ norm equation over local field extensions###########################

function norm_equation(R:: Hecke.LocalField, b::Union{qadic,Hecke.LocalFieldElem})
   K = parent(b)
   prec_b = precision(b)
   @assert base_field(R) == K  #"since trace(a,_) is not defined"
   f,mf = ResidueField(K)
   F,mF = ResidueField(R)
   ee = absolute_ramification_index(K)
   if degree(R) == ramification_index(R) && mf(b) !=f(1)
      error("To be implemented")
   end
   if mf(b) == f(1)
      f_nm = R(1)
   else
      f_nm = norm_equation(F,mf(b))
      f_nm = mF\(f_nm)
   end
   b = b//norm(f_nm)
   b = setprecision(b,prec_b)
   p = prime(R)
   if valuation(b-1) < ee//(p-1)+1//ee
      error("To be implemented or try norm_equation_unramified")
   end
   r = random_elem(R)
   while valuation(trace(r)) != 0 || valuation(r//R(trace(r))) != 0
      r = random_elem(R)
   end
   s = r*R(trace(r)^-1)*R(log(b))
   return exp(s)*f_nm
end

###########################################################################################################
#   The following "norm_equation_unramified" solves the norm equations only in unramified extensions
###########################################################################################################

function norm_equation_unramified(L::Hecke.LocalField, b::Hecke.LocalFieldElem)
   K = parent(b)
   @assert degree(L) == inertia_degree(L)
   prec_b = precision(b)
   piK = uniformizer(K)
   piL = uniformizer(L)
   f,mf = ResidueField(K)
   F,mF = ResidueField(L)
   ee = absolute_ramification_index(K)
   if mf(b) == f(1)
      f_nm = L(1)
   else
      f_nm = norm_equation(F,mf(b))
      f_nm = mF\(f_nm)
   end
   b = b//norm(f_nm)
   b = setprecision(b,prec_b)
   c = L(1)
   C = [L(1)]
   n = ee*valuation((b//norm(C[1]))-1)
   r = random_elem(L)
   while valuation(trace(r)) != 0 || valuation(r//L(trace(r))) != 0
      r = random_elem(L)
   end
   z = ((b//norm(c))-1)//piK^ZZ(n)
   setprecision!(z, prec_b)
   push!(C, 1+ piL^ZZ(n)* (L(z)*r//L(trace(r))))
   c = prod(C)
   nc = norm(c)
   n = ee*valuation((b//nc)-1)
   while prec_b >= n+1 #  "solve till precision n-1"
      z = ((b//nc)-1)*piK^-ZZ(n)
      setprecision!(z, prec_b)
      push!(C, 1+ piL^ZZ(n)*(L(z)*r//L(trace(r))))
      c = prod(C)
      setprecision!(c, precision(L))
      nc = norm(c)
      chk = (nc*b^-1)-1
      if iszero(chk) == true
         n = prec_b
      else
         n = ee*valuation((b//nc)-1)
      end
   end
   return c*f_nm
end

ext: jl


text:
export computeHinvTRec

function computeHinvTRec(pMis::Array{RemoteChannel})
HinvTRec = Array{Array{ComplexF64}}(undef,length(pMis))
@sync begin
	for k=1:length(pMis)
		@async begin
			HinvTRec[k] = remotecall_fetch(computeHinvTRec,pMis[k].where,pMis[k]);
		end
	end
end
return HinvTRec;
end

function calculateShiftedH(pMis, m)
	pFor = pMis.pFor;
	M       	= pFor.Mesh
	omega   	= pFor.omega
	wavelet 	= pFor.WaveletCoef;
	gamma   	= pFor.gamma
	Q       	= pFor.Sources
	P       	= pFor.Receivers
	Ainv    	= pFor.ForwardSolver;
	batchSize 	= pFor.forwardSolveBatchSize;
	select  	= pFor.sourceSelection;
	An2cc = getNodalAverageMatrix(M);

	m, = pMis.modelfun(m)
	m = interpGlobalToLocal(m,pMis.gloc.PForInv,pMis.gloc.sigmaBackground);

	m = An2cc'*m;
	gamma = An2cc'*gamma;

	H = GetHelmholtzOperator(M,m,omega, gamma, true,useSommerfeldBC);
	Ainv.helmParam = HelmholtzParam(M,gamma,m,omega,true,useSommerfeldBC);
	H = H + GetHelmholtzShiftOP(m, omega,Ainv.shift[1]);
	H = sparse(H');

	return H;
end

function computeHinvTRec(pMisRF::RemoteChannel)
pMis = fetch(pMisRF)
HinvTP, = solveLinearSystem(spzeros(ComplexF64,0,0),complex(Matrix(pMis.pFor.Receivers)),pMis.pFor.ForwardSolver,1);
return HinvTP;
end

function computeHinvTRec(pMisRF::RemoteChannel, m)
pMis = fetch(pMisRF)
H = calculateShiftedH(pMis, m)
HinvTP, = solveLinearSystem(H,complex(Matrix(pMis.pFor.Receivers)),pMis.pFor.ForwardSolver,1);
return HinvTP;
end

function computeHinvTRec(pMisRF::RemoteChannel, x, m, doTranspose)
pMis = fetch(pMisRF)
if doTranspose == 0
	HinvTP, = solveLinearSystem(spzeros(ComplexF64,0,0),x,pMis.pFor.ForwardSolver,doTranspose);
	return complex(pMis.pFor.Receivers' * HinvTP);
else
	HinvTP, = solveLinearSystem(spzeros(ComplexF64,0,0),complex(pMis.pFor.Receivers * x),pMis.pFor.ForwardSolver,doTranspose);
	return  HinvTP;
end
end

function computeHinvTRecX(pMis::Array{RemoteChannel}, x, m, doTranspose=1)
HinvTRec = Array{Array{ComplexF64}}(undef,length(pMis))
@sync begin
		for k=1:length(pMis)
	@async begin
			HinvTRec[k] = remotecall_fetch(computeHinvTRec,pMis[k].where,pMis[k],x, m,doTranspose);
		end
	end
end
return HinvTRec;
end

function computeHinvTRecXarr(pMis::Array{RemoteChannel}, x, m, doTranspose=1)
HinvTRec = Array{Array{ComplexF64}}(undef,length(pMis))
@sync begin
	for k=1:length(pMis)
	@async begin
			HinvTRec[k] = remotecall_fetch(computeHinvTRec,pMis[k].where,pMis[k],x[k], m,doTranspose);
		end
	end
end
return HinvTRec;
end


text:
## #__precompile__()
#"""
## tsZeroCross Module
#
# Zero crossing analysis routines including many outputs: heights, periods, peaks, troughs, etc
#
#---
#
### Functions
#
#    o=tsCalcStats(z)
#    o=tsCalcZeroCrossDiscrete(z)
#    o=tsCalcZeroCrossFull(t,z)
#
### Types
#
#    tsStats -  Type containing basic statistics of the input z-signal
#    zCross - Type containing zero crossing minima, maxima and height data
#    zCrossDiscrete - Type containing zero crossing low level data
#    zCrossUpDown - Type containing height and period data and statistics
#    ZeroCrossFull - Type containing full zero crossing data
#
### Description
#
# An OM_JPL Function.
#
# Zero crossing analysis routine including many outputs: heights, periods, peaks, troughs, etc
#
# Determines zero crossing parameters from the input time series. Both up crossing
#  and downcrossing is peformed and reported.
#
# If the start and end of the time series is not at zero then these are not included in the
#  statistics.  Only full zero-crossing wave cycles are used and reported.
#
# Values of z which are Inf or NaN are treated as missing values. If a zero-crossing occurs between
#  these times linear interpolation is used to find the crossing point.
#
### Examples
#
#```jldoctest
#julia> f=1./[2;4; 8; 10; 15; 20]
#julia> a= 0.5 * [1; 2; 3; 6; 2; 1]
#julia> p=rand(6)*2*pi
#julia> t=[xi for xi=0:0.1:110 ]
#julia> z=0*t
#julia> for i=1:size(a,1)
#julia> z.=z.+a[i].*sin.(2.*pi.*f[i].*t+p[i])
#julia> end
#julia> z=z-mean(z)
#julia> oStats=tsCalcStats(z)
#julia> oDis=tsCalcZeroCrossDiscrete(z)
#julia> oFull=tsCalcZeroCrossFull(t,z)
#```
#
### History
#
# Created by Jason McConochie
#   Revision 1, 20Dec2016 , Jason McConochie
#      based upon Revision 3, of matlab version
#"""

#using PyPlot
########################################################################################
# Functions and types for basic time series data statistics
#export tsCalcStats, tsStats
"""
# tsStats
Type containing basic statistics of the input z-signal
  mean median max imax min imin std rms var
"""
type tsStats
   mean::Float64
   median::Float64
   max::Float64
   imax::Float64
   min::Float64
   imin::Float64
   std::Float64
   rms::Float64
   var::Float64
   function tsStats()
     return new()
   end
end

########################################################################################
# Functions and types for discrete zero crossing analysis using z only, no t
#export tsCalcZeroCrossDiscrete, zCrossDiscrete, zCross
# Type to contain more detail but usually less relevant information
#  this is made by the zerocross_discrete function
"""
# zCross
Type containing zero crossing minima, maxima and height data
"""
type zCross
   zmax::Array{Float64,1}
   zmin::Array{Float64,1}
   izmax::Array{Int64,1}
   izmin::Array{Int64,1}
   height::Array{Float64,1}
   function zCross()
     return new()
   end
end
"""
# zCrossDiscrete

Type containing zero crossing low level data
     .izc       indicies of the points surrounding all zero crossings selected
     .zctype      value indicating the type of crossing
     .zctypemeta  ' 1 is an upcrossing, -1 is a downcrossing'
     .ilcrest   indicies of all local crests identified
     .iltrough  indicies of all local troughs identified
     .ipslope   boolean indicating postively sloping z
     .i0slope   boolean indicating sloping z
     .inslope   boolean indicating sloping z
     .uc 	contains zCross structure foe up-crossings
     .dz   	contains zCross structure foe down-crossings
"""
type zCrossDiscrete
   # low-level zero crossing data
   izc::Array{Int64,2}
   zctype::Array{Int64,1}
   zctypemeta::AbstractString
   ilcrest::Array{Int64,1}
   iltrough::Array{Int64,1}
   # slopes
   ipslope::Array{Bool,1}
   i0slope::Array{Bool,1}
   inslope::Array{Bool,1}
   # crossing data
   uc::zCross
   dc::zCross
   function zCrossDiscrete()
     t=new()
     t.uc=zCross()
     t.dc=zCross()
     return t
   end
end

########################################################################################
# Functions and types for full zero crossing analysis using t,z
#export tsCalcZeroCrossFull, ZeroCrossFull
"""
# zCrossUpDown

Type containing height and period data and statistics

       .H         vector of individual wave heights
       .T         vector of individual wave periods

       .Hmax      maximum individual wave height
       .Hm        mean wave height
       .H1on10        mean of the highest wave exceeding a 1/10 probability of exceedence
       .H1on3     mean of the highest wave exceeding a 1/3 probability of exceedence
       .H1on2     mean of the highest wave exceeding a 1/2 probability of exceedence

       .THmax     period of Hmax
       .Tm        mean wave period
       .T1on10    mean wave period of all waves with heights exceeding a 1/10 probability of exceedence
       .T1on3     mean wave period of all waves with heights exceeding a 1/3 probability of exceedence
       .T1on2     mean wave period of all waves with heights exceeding a 1/2 probability of exceedence

       .zcTimes   vector of times of zero crossings (either up or down)

"""
type zCrossUpDown
   H::Array{Float64,1}
   T::Array{Float64,1}
   # height statistics
   Hmax::Float64
   Hm::Float64
   H1on10::Float64
   H1on3::Float64
   H1on2::Float64
   # period statistics
   THmax::Float64
   Tm::Float64
   T1on10::Float64
   T1on3::Float64
   T1on2::Float64
   # time of zero crossings
   zcTimes::Array{Float64,1}
   function zCrossUpDown()
     return new()
   end
end
"""
# ZeroCrossFull
 Type containing full zero crossing data
"""
type ZeroCrossFull
   # time of zero crossings
   zcTimes::Array{Float64,1}
   upCross::zCrossUpDown
   downCross::zCrossUpDown
   zCrossDiscrete::zCrossDiscrete
   Hsigma::Float64
   function ZeroCrossFull()
     t=new()
     t.upCross=zCrossUpDown()
     t.downCross=zCrossUpDown()
     t.zCrossDiscrete=zCrossDiscrete()
     return t
   end
end

########################################################################################
# Basic Statistics - tsCalcStats
########################################################################################
"""
# tsCalcStats

Calculate basic statistics min, max, std, var, mean, median

"""
function tsCalcStats(z::Array{Float64,1})::tsStats
  # 0. Make new output structure
  o=tsStats()
  # A. Find typical statistics
  o.mean=mean(z)
  o.median=median(z)
  o.max=maximum(z)
  o.imax=indmax(z)
  o.min=minimum(z)
  o.imin=indmin(z)
  o.std=std(z)
  a=0
  for i=1:length(z)
    a=a+z[i].^2
  end
  o.rms=(a/length(z))^0.5
  o.var=var(z)
  return o
end

########################################################################################
# Discrete Zero Crossing - tsZeroCrossDiscrete
########################################################################################
"""

# tsCalcZeroCrossDiscrete

 Function for discrete time series analysis

## INPUT

  z - Vector of time series values

"""
function tsCalcZeroCrossDiscrete(z::Array{Float64,1})::zCrossDiscrete

  # 0. Input handling
  o=zCrossDiscrete()

  # A. Return empty data if there is only none or point or all zeros
  nz=size(z,1)
  allZero=true
  for i=1:size(z,1)
    if z[i]!=0
      allZero=false
      break
    end
  end
  if allZero
    return o
  end

  # C. Find zero-crossing indicies
  #sz=convert(Array{Int64,1},sign(z)) # -1 -1 0 0 1 1
  sz=sign(z) # -1 -1 0 0 1 1
  sz[1]==0 ? sz[1]=sz[2] : 0
  # ..Find patterns of -1 n0 -1 or 1 n0 1 in sz - remove these from
  foundzero=false
  izero=0
  ncross=0
  for i=2:nz
    if sz[i]==0 && !foundzero
      foundzero=true
      izero=i
    elseif (sz[i]!=0) && foundzero
      sz[izero:(i-1)]=sz[i]
      foundzero=false
    end
    abs(sz[i]-sz[i-1])>0 ? ncross+=1 : 0
  end

  icross::Array{Int64,1}=zeros(Int64,ncross)
  kcross=0
  for i=2:nz
    if abs(sz[i]-sz[i-1])>0
      kcross+=1
      icross[kcross]=i-1
    end
  end

  # D. Plotting check
  #println("Ncrossings=$(ncross)")
  #PyPlot.plot(1:1:length(z),z,color="black")
  #hold(true)
  #PyPlot.plot(icross,icross.*0,"rx")
  #PyPlot.plot([1;length(z)],[0;0])

  # E. Put into icross array of indicies which bound the zero-crossing
  #	 start1,end1
  #	 start2,end2 etc
  o.izc=zeros(Int64,(ncross,2))
  for i=1:ncross
     o.izc[i,1]=icross[i]
     o.izc[i,2]=icross[i]+1
  end

  # F1. Check if the first point is zero-touching - and include
  if z[1]==0
    o.izc=[1 1;o.izc];tcross1=convert(Int64,sign(z[2]))
  end
  # F2. Check if last point is zero touching
  if z[end]==0
    o.izc[end,:]=[length(z) length(z)];tcrossN=-1*convert(Int64,sign(z[end-1]))
  end

  # G. Got all zero-crossing indicies now mark which is upcrossing and which is
  #  down crossing
  # 0  - touching zero
  # 1 - up-crossing
  # -1	- down-crossing

  o.zctype=Array(Int64,size(o.izc,1))
  for i=1:size(o.izc,1)
    o.zctype[i]=-1*sign(z[o.izc[i,1]])
  end

  # H1. Assign zero-touches to up or down crossing
  z[1]==0 ? o.zctype[1]=tcross1 : 0
  z[end]==0 ? o.zctype[end]=tcrossN : 0
  # H2. Check if we are still okay
  if !all(abs(o.zctype[2:end]-o.zctype[1:end-1]).==2)
   println("**** Serious problem with zero crosssing routine")
   return
  end

  o.zctypemeta=" 1 is an upcrossing, -1 is a downcrossing"

  # I1. Find crests and troughs
  ra1::Array{Bool,1}=falses(nz)
  ra2::Array{Bool,1}=falses(nz)

  # I2.1 Find all 121 style crests
  for i=1:nz-2
    ra1[i]=(z[i]<z[i+1]) & (z[i+1]>z[i+2])
  end
  # I2.2 Find all 1221 style crests
  for i=1:nz-3
    ra2[i]=(z[i]<z[i+1]) & (z[i+1]==z[i+2]) & (z[i+2]>z[i+3])
  end
  o.ilcrest=find(ra1 | ra2)+1

  # I3.1 Find all 212 style troughs
  for i=1:nz-2
    ra1[i]=(z[i]>z[i+1]) & (z[i+1]<z[i+2])
  end
  # I3.2 Find all 2112 style crests
  for i=1:nz-3
    ra2[i]=(z[i]>z[i+1]) & (z[i+1]==z[i+2]) & (z[i+2]<z[i+3])
  end
  o.iltrough=find(ra1 | ra2)+1

  # J. Categorise slope
  o.ipslope=falses(nz)
  o.i0slope=falses(nz)
  o.inslope=falses(nz)
  for i=1:nz-2
    # .. find all points with neutral gradient
    o.ipslope[i]=z[i] < z[i+1]
    # .. find all points with neutral gradient
    o.i0slope[i]=z[i] == z[i+1]
    # .. find all points with negative gradient
    o.inslope[i]=z[i] > z[i+1]
  end

  # K. Find all zero-crossing heights
  #    construct array of waves
  o.uc=zCross()
  o.dc=zCross()

  # K1. upcrossing heights
  #    ..from first up cross to last upcross
  # .. first upcross
  iuc=find(o.zctype.==1)  #upcrossings
  nst=size(iuc,1)-1
  o.uc.zmax=zeros(Float64,nst)
  o.uc.zmin=zeros(Float64,nst)
  o.uc.izmax=zeros(Int64,nst)
  o.uc.izmin=zeros(Int64,nst)
  for i=1:nst
    st=o.izc[iuc[i],2]   #first point after upcross
    ed=o.izc[iuc[i+1],1] #last point in upcross wave
    mz=z[st:ed]
    o.uc.zmax[i]=maximum(mz)
    o.uc.zmin[i]=minimum(mz)
    o.uc.izmax[i]=st+indmax(mz)-1
    o.uc.izmin[i]=st+indmin(mz)-1
  end
  o.uc.height=o.uc.zmax-o.uc.zmin

  # K2. downcrossing heights
  #   ..from first up cross to last upcross
  # .. first upcross
  iuc=find(o.zctype.==-1)  #down-crossings
  nst=size(iuc,1)-1
  o.dc.zmax=zeros(Float64,nst)
  o.dc.zmin=zeros(Float64,nst)
  o.dc.izmax=zeros(Int64,nst)
  o.dc.izmin=zeros(Int64,nst)
  for i=1:nst
    st=o.izc[iuc[i],2]        #first point after upcross
    ed=o.izc[iuc[i+1],1]      #last point in upcross wave
    mz=z[st:ed]
    o.dc.zmax[i]=maximum(mz)
    o.dc.zmin[i]=minimum(mz)
    o.dc.izmax[i]=st+indmax(mz)-1
    o.dc.izmin[i]=st+indmin(mz)-1
  end
  o.dc.height=o.dc.zmax-o.dc.zmin

  #plot(z)
  #hold("true")
  #plot(o.izc[:,1],o.izc[:,1]*0,"rx")

  #println("Done $(o)")

  return o

end # end function dzca

########################################################################################
# Full Zero Crossing - tsZeroCross
########################################################################################
"""

# tsCalcZeroCrossFull

 Function for full zero crossing analysis of time series

## INPUT

  t - Vector of times as Float64
  z - Vector of time series values

"""
function tsCalcZeroCrossFull(t::Array{Float64,1},z::Array{Float64,1})::ZeroCrossFull

     # 0. Do some checking
     t=t[:]
     z=z[:]
     if size(t,1)!=size(z,1)
         return NaN
         #error('tsZeroCross: length of t (%d) does not match z (%d)\n',length(t),length(z));
     end
     # .. check for unique and monotonically increasing values of t
     bad=any(diff(t).<=0.)
     if bad
         return NaN
         # error('tsZeroCross: t is not monotonically increasing\n');
     end
     # .. check for NaN's/Inf's
     bad1=any(isnan(z))
     bad2=any(isinf(z))
     if bad1 | bad2
         println("tsZeroCross: NaNs or Infs detected in z, removing them and treating as missing values")
         m=isnan(z) | isinf(z)
         ind=1:1:size(m,1)
         deleteat!(z,ind(m))
         deleteat!(t,ind(m))
     end

     # A. Perform discrete part of time series analysis
     # A1. Make new output strcuture
     om=ZeroCrossFull()
     # A2. Run the discrete zero crossing
     # .. this only requires the z component
     om.zCrossDiscrete=tsCalcZeroCrossDiscrete(z)
     # B. Find zero crossing times and up and down crossing periods
     if size(om.zCrossDiscrete.izc,1)>1
         # B1. Find zero crossing times
         x1=z[om.zCrossDiscrete.izc[:,1]]
         x2=z[om.zCrossDiscrete.izc[:,2]]
         y1=t[om.zCrossDiscrete.izc[:,1]]
         y2=t[om.zCrossDiscrete.izc[:,2]]
         # ..find slope of known data points
         Slope = (y2 - y1) ./ (x2 - x1)
         # ..find intercept
         Intercept = y2 - Slope .* x2
         inan=isnan(Intercept)
         Intercept[inan]=0
         Intercept=inan.*y1+Intercept
         om.zcTimes=Intercept

         # B2. Upcrossing Periods
         ic=find(om.zCrossDiscrete.zctype.==1)  # upcrossings
         if !isempty(ic) # must have at least one upcrossing
             om.upCross.zcTimes=om.zcTimes[ic]
             om.upCross.H=om.zCrossDiscrete.uc.height
         end
         if isdefined(om.upCross,:zcTimes)
             om.upCross.T=diff(om.upCross.zcTimes)
         end

         # B3. Downcrossing Periods
         ic=find(om.zCrossDiscrete.zctype.==-1)  # downcrossings
         if !isempty(ic) # must have at least one downcrossing
             om.downCross.zcTimes=om.zcTimes[ic]
             om.downCross.H=om.zCrossDiscrete.dc.height
         end
         if isdefined(om.downCross,:zcTimes)
             om.downCross.T=diff(om.downCross.zcTimes)
         end
     end

     # C. Find periods statistics
     # C1. Upcrossing period statistics
     if !isdefined(om,:upCross)
         om.upCross.Tm=mean(om.upCross.T)
         isp=sortperm(om.upCross.H)
         N=size(isp,1)
         ist=ceil(N-N/10+0.5)
         om.upCross.T1on10=mean(om.upCross.T[ist:N])
         ist=ceil(N-N/3+0.5)
         om.upCross.T1on3=mean(om.upCross.T[ist:N])
         ist=ceil(N-N/2+0.5)
         om.upCross.T1on2=mean(om.upCross.T[ist:N])
     end
     # C2. Downcrossing period statistics
     if !isdefined(om,:downCross)
         om.downCross.Tm=mean(om.downCross.T)
         isp=sortperm(om.downCross.H)
         N=size(isp,1)
         ist=ceil(N-N/10+0.5)
         om.downCross.T1on10=mean(om.downCross.T[ist:N])
         ist=ceil(N-N/3+0.5)
         om.downCross.T1on3=mean(om.downCross.T[ist:N])
         ist=ceil(N-N/2+0.5)
         om.downCross.T1on2=mean(om.downCross.T[ist:N])
     end

     # D. Find height statistics
     # D1. UpCrossing Height statistics
     if !isdefined(om,:upCross)
         om.upCross.Hm=mean(om.upCross.H)
         sh=sort(om.upCross.H)
         N=size(sh,1)
         ist=ceil(N-N/10+0.5)
         om.upCross.H1on10=mean(sh[ist:N])
         ist=ceil(N-N/3+0.5)
         om.upCross.H1on3=mean(sh[ist:N])
         ist=ceil(N-N/2+0.5)
         om.upCross.H1on2=mean(sh[ist:N])
     end
     # D2. DownCrossing Height statistics
     if !isdefined(om,:downCross)
         om.downCross.Hm=mean(om.downCross.H)
         sh=sort(om.downCross.H)
         N=size(sh,1)
         ist=ceil(N-N/10+0.5)
         om.downCross.H1on10=mean(sh[ist:N])
         ist=ceil(N-N/3+0.5)
         om.downCross.H1on3=mean(sh[ist:N])
         ist=ceil(N-N/2+0.5)
         om.downCross.H1on2=mean(sh[ist:N])
     end

     # E. Find data for the maximum wave
     # E1. Upcrossing
     if !isdefined(om,:upCross)
         om.upCross.Hmax=max(om.upCross.H)
         om.upCross.iHmax=indmax(om.upCross.H)
         om.upCross.THmax=om.upCross.T[om.upCross.iHmax]
     end
     # E2. Downcrossing
     if !isdefined(om,:downCross)
         om.downCross.Hmax=max(om.downCross.H)
         om.downCross.iHmax=indmax(om.downCross.H)
         om.downCross.THmax=om.downCross.T[om.downCross.iHmax]
     end

     # F. Calculate significant wave height
     om.Hsigma=4*sqrt(var(z))

     # G. If requested make a plot of all of the information
     #if exist('doplot','var')==1
     #    figure
         #ind=1:1:length(z);
     #    plot(t,z,'k-x');hold on
         #plot(t(dtsa.ZeroUpCrossingIndicies_1st),z(dtsa.ZeroUpCrossingIndicies_1st),'ro');
         #plot(t(dtsa.ZeroDownCrossingIndicies_1st),z(dtsa.ZeroDownCrossingIndicies_1st),'b+');
     #    plot(t(o.dzc.ilcrest),z(o.dzc.ilcrest),'g^');
     #    plot(t(o.dzc.iltrough),z(o.dzc.iltrough),'cv');

     #    plot(o.uc.tzc,zeros(size(o.uc.tzc)),'+');
         #h=text(tsa.ZeroUpCrossingTimes,zeros(size(tsa.ZeroUpCrossingTimes)),char(tsa.Periods));
         #set(h,'verticalalignment','bottom','horizontalalignment','left');
     #    h=line(t,zeros(size(t)));
     #    set(h,'color',[0 0 0]);

     #end
     return om

 end

max_stars_repo_name: jadamc/MetoceanTools

max_issues_repo_name: jadamc/MetoceanTools

max_forks_repo_name: jadamc/MetoceanTools


text:
import DSP
import Multirate

sampleRate    = 48000
interpolation = 147
decimation    = 160
ratio         = interpolation  // decimation
numTaps       = 24*interpolation
x             = rand( 1_000_000 )
h             = Multirate.firdes( numTaps, 0.5/interpolation, Multirate.kaiser, beta = 7.8562  )

function naiveresampler{T}( h::Vector{T}, x::Vector{T}, ratio::Rational{Int} )
    upfactor   = numerator( ratio )
    downfactor = denominator( ratio )
    xStuffed   = zeros( T, length(x) * upfactor )

    for n = 0:length(x)-1;
        xStuffed[ n*upfactor+1 ] = x[ n+1 ]
    end

    yInterpolated = DSP.filt( h, xStuffed )
    y = [ yInterpolated[n] for n = 1:downfactor:length( yInterpolated ) ]
end

@time y = naiveresampler( h, x, 147//160 );
@time y = Multirate.filt( h, x, ratio );


text:
lenn_file   = string(path, "/SampleFiles/ASCII/0215162000.c00")
geocsv_file = string(path, "/SampleFiles/ASCII/geo-tspair.csv")
slist_file  = string(path, "/SampleFiles/ASCII/2m-62.5hz.slist")
lenn_pat    = string(path, "/SampleFiles/ASCII/021516*c00")
geocsv_pat  = string(path, "/SampleFiles/ASCII/geo-tspair.*")
slist_pat   = string(path, "/SampleFiles/ASCII/*.slist")

printstyled("  Lennartz ASCII\n", color=:light_green)
C = read_data("lennartz", lenn_file)[1]
@test ≈(C.fs, 62.5)

printstyled("    wildcard support\n", color=:light_green)
S = read_data("lennartz", lenn_pat)

printstyled("  GeoCSV timeseries\n", color=:light_green)
S = read_data("geocsv", geocsv_file)
@test S.n == 8
if Sys.iswindows()
  println("    IDs = ", join(S.id, ","))
end
i = findid("CC.JRO..BHZ", S.id)
if i > 0
  @test ≈(S.fs[i], 50.0)
  @test ==(S.loc[i], GeoLoc(lat=46.275269, lon=-122.218262, el=1219.0, inc=180.0))
end

printstyled("    wildcard support\n", color=:light_green)
S = read_data("geocsv", geocsv_pat)

printstyled("  slist\n", color=:light_green)
S = read_data("slist", slist_file)
nx = length(S.x[1])
@test S.id[1] == "YY.ERTA..EHZ"
@test ≈(S.fs[1], 62.5)
@test isapprox(C.x[1:nx], S.x[1])

printstyled("    wildcard support\n", color=:light_green)
S = read_data("slist", slist_pat)
@test S.id[1] == "YY.ERTA..EHZ"
@test ≈(S.fs[1], 62.5)
@test isapprox(C.x[1:nx], S.x[1])


text:
@testset "Spelled pitches" begin
    @testset "types ok and compact" begin
        @test isbitstype(SpelledInterval)
        @test isbitstype(SpelledIC)
        @test isbitstype(Pitch{SpelledInterval})
        @test isbitstype(Pitch{SpelledIC})
    end

    @testset "constructors" begin
        @test spelled(-3, 0) == spelled(-3, 0)
        @test spelledp(0, 4) == spelledp(0, 4)
        @test sic(-3) == sic(-3)
        @test spc(3) == spc(3)
    end

    @testset "string macros" begin
        @test i"M3:1" == spelled(4, -1)
        @test i"-M3:0" == spelled(-4, 2)
        @test p"C♭4" == spelledp(-7, 8)
        @test p"Cb4" == spelledp(-7, 8)
        @test i"m3" == sic(-3)
        @test i"-m3" == sic(3)
        @test p"C♯" == spc(7)
        @test p"C#" == spc(7)
        # TODO: test more complex parsing cases
    end

    @testset "accessors" begin
        @test octaves(i"M3:1") == 1
        @test internalocts(i"M3:1") == -1
        @test fifths(i"M3:1") == 4
        @test degree(i"M3:1") == 2
        @test generic(i"M3:1") == 2
        @test diasteps(i"M3:1") == 9
        @test alteration(i"M3:1") == 0

        @test octaves(i"-M3:1") == -2
        @test internalocts(i"-M3:1") == 1
        @test fifths(i"-M3:1") == -4
        @test degree(i"-M3:1") == 5
        @test generic(i"-M3:1") == -2
        @test diasteps(i"-M3:1") == -9
        @test alteration(i"-M3:1") == 0

        @test octaves(i"a5") == 0
        @test internalocts(i"a5") == 0
        @test fifths(i"a5") == 8
        @test degree(i"a5") == 4
        @test generic(i"a5") == 4
        @test diasteps(i"a5") == 4
        @test alteration(i"a5") == 1

        @test octaves(p"Ebb5") == 5
        @test fifths(p"Ebb5") == -10
        @test degree(p"Ebb5") == 2
        @test alteration(p"Ebb5") == -2
        @test letter(p"Ebb5") == 'E'

        @test octaves(p"F#") == 0
        @test fifths(p"F#") == 6
        @test degree(p"F#") == 3
        @test alteration(p"F#") == 1
        @test letter(p"F#") == 'F'

        # edge cases
        @test alteration(i"P4") == 0
        @test alteration(i"M7") == 0
        @test alteration(i"-P4:0") == 0
        @test alteration(i"-M7:0") == 0

        @test alteration(i"a4:0") == 1
        @test alteration(i"m7:0") == -1
        @test alteration(i"-a4:0") == 1
        @test alteration(i"-m7:0") == -1

        @test alteration(p"F") == 0
        @test alteration(p"B") == 0
    end

    @testset "printing" begin
        @test string(i"m3:1") == "m3:1"
        @test string(p"Eb4") == "E♭4"
        @test string(i"m3") == "m3"
        @test string(p"E#") == "E♯"
    end

    @testset "hashing" begin
        @test hash(spelled(-3,1)) == hash(spelled(-3,1))
        @test hash(spelledp(3,0)) == hash(spelledp(3,0))
        @test hash(sic(3)) == hash(sic(3))
        @test hash(spc(3)) == hash(spc(3))
    end

    @testset "interval interface" begin
        @test i"m3:0" + i"M3:0" == i"P5:0"
        @test i"m3:0" + i"M7:0" == i"M2:1"
        @test i"P5:0" + i"P5:0" == i"M2:1"
        @test i"-m3:0" + i"M3:0" == i"a1:0"
        @test i"m3:0" + i"-M3:0" == i"-a1:0"
        @test i"m3:0" - i"M3:0" == i"-a1:0"
        @test i"m3:0" - i"M6:0" == i"-a4:0"
        @test -i"P4:0" == i"-P4:0"
        @test -i"P4:0" == i"P5:-1"
        @test -i"P5:0" == zero(SpelledInterval) - i"P5:0"

        @test zero(SpelledInterval) == i"P1:0"
        @test zero(i"m3:0") == i"P1:0"

        @test i"P5:0" * 2 == i"M2:1"
        @test i"M2:0" * 4 == i"a5:0"
        @test i"-m3:0" * 4 == i"-d2:1"
        @test i"M3:0" * -3 == i"-a7:0"
        @test 4 * i"M2:0" == i"a5:0"
        @test 4 * i"-M3:0" == i"-aa2:1"
        @test 5 * i"M3:0" == i"aaa4:1"

        @test tomidi(i"aaa4:1") == midi(20)
        @test tomidi(i"-P5:0") == midi(-7)
        @test octave(SpelledInterval) == i"P1:1"
        @test sign(i"m2:0") == 1
        @test sign(i"P1:0") == 0
        @test sign(i"d1:0") == 0
        @test sign(i"a1:0") == 0
        @test sign(i"-m3:0") == -1
        @test abs(i"-m3:0") == i"m3:0"

        @test ic(i"M3:3") == i"M3"
        @test ic(i"-M3:1") == i"m6"
        @test embed(i"M3:3") == i"M3:3"
        @test intervaltype(SpelledInterval) == SpelledInterval
        @test intervalclasstype(SpelledInterval) == SpelledIC

        @test isstep(i"d1:0")
        @test isstep(i"P1:0")
        @test isstep(i"a1:0")
        @test isstep(i"d2:0")
        @test isstep(i"m2:0")
        @test isstep(i"M2:0")
        @test isstep(i"a2:0")
        @test isstep(i"-d2:0")
        @test isstep(i"-m2:0")
        @test isstep(i"-M2:0")
        @test isstep(i"-a2:0")

        @test !isstep(i"d3:0")
        @test !isstep(i"-d3:0")
        @test !isstep(i"M7:0")
        @test !isstep(i"-M7:0")
        @test !isstep(i"P1:1")
        @test !isstep(i"-P1:1")
        @test !isstep(i"m2:1")
        @test !isstep(i"-m2:1")

        @test chromsemi(SpelledInterval) == i"a1:0"
    end

    @testset "interval class interface" begin
        @test i"m3" + i"M3" == i"P5"
        @test i"m3" + i"M7" == i"M2"
        @test i"P5" + i"P5" == i"M2"
        @test i"-m3" + i"M3" == i"a1"
        @test i"m3" + i"-M3" == i"-a1"
        @test i"m3" - i"M3" == i"-a1"
        @test i"m3" - i"M6" == i"-a4"
        @test -i"P4" == i"-P4"
        @test -i"P4" == i"P5"
        @test -i"P5" == zero(SpelledIC) - i"P5"

        @test zero(SpelledIC) == i"P1"
        @test zero(i"m3") == i"P1"

        @test i"P5" * 2 == i"M2"
        @test i"M2" * 4 == i"a5"
        @test i"-m3" * 4 == i"-d2"
        @test i"M3" * -3 == i"-a7"
        @test 4 * i"M2" == i"a5"
        @test 4 * i"-M3" == i"-aa2"
        @test 5 * i"M3" == i"aaa4"

        @test tomidi(i"aaa4") == midic(8)
        @test tomidi(i"-P5") == midic(-7)
        @test octave(SpelledIC) == i"P1"
        @test sign(i"m2") == 1
        @test sign(i"P1") == 0
        @test sign(i"d1") == 0
        @test sign(i"a1") == 0
        @test sign(i"-m3") == -1
        @test abs(i"-m3") == i"m3"

        @test ic(i"-M3") == i"m6"
        @test embed(i"M3") == i"M3:0"
        @test embed(i"M7") == i"M7:0"
        @test embed(i"P4") == i"P4:0"
        @test embed(i"a1") == i"a1:0"
        @test intervaltype(SpelledIC) == SpelledInterval
        @test intervalclasstype(SpelledIC) == SpelledIC

        @test isstep(i"d1")
        @test isstep(i"P1")
        @test isstep(i"a1")
        @test isstep(i"d2")
        @test isstep(i"m2")
        @test isstep(i"M2")
        @test isstep(i"a2")
        @test isstep(i"-d2")
        @test isstep(i"-m2")
        @test isstep(i"-M2")
        @test isstep(i"-a2")

        @test !isstep(i"d3")
        @test !isstep(i"-d3")

        @test chromsemi(SpelledIC) == i"a1"
    end

    @testset "pitch interface" begin
        @test topitch(i"m3:4") == p"Eb4"
        @test tointerval(p"C#3") == i"a1:3"

        @test p"Eb4" + i"P5:0"  == p"Bb4"
        @test p"Eb4" + i"-m3:0" == p"C4"
        @test p"Eb4" - i"P5:0"  == p"Ab3"
        @test p"G4"  - p"C#4"   == i"d5:0"

        @test alteration(p"Ab-1") == -1
        @test alteration(p"A#-1") == 1

        @test pc(p"Eb4") == p"Eb"
        @test embed(p"Eb4") == p"Eb4"

        @test tomidi(p"C#3") == midip(49)
        @test tomidi(p"Db3") == midip(49)
    end

    @testset "pitch class interface" begin
        @test topitch(i"m3") == p"Eb"
        @test tointerval(p"E") == i"M3"

        @test p"Eb" + i"P5"  == p"Bb"
        @test p"Eb" + i"-m3" == p"C"
        @test p"Eb" - i"P5"  == p"Ab"
        @test p"G"  - p"C#"   == i"d5"

        @test pc(p"Eb") == p"Eb"
        @test embed(p"Eb") == p"Eb0"
        @test embed(p"Eb", 4) == p"Eb4"

        @test tomidi(p"C#") == midipc(1)
        @test tomidi(p"Db") == midipc(1)
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "SpglibBuilder"
version = v"0.2.0"

# Collection of sources required to build SpglibBuilder
sources = [
    "https://github.com/atztogo/spglib/archive/v1.14.1.tar.gz" =>
    "9803b0648d9c2d99377f3e1c4cecf712320488403cd674192ec5cbe956bb3c78",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd spglib-1.14.1/
if [[ ${target} == *-mingw32 ]]; then
    sed -i -e 's/LIBRARY/RUNTIME/' CMakeLists.txt
fi
mkdir _build
cd _build/
cmake -DCMAKE_INSTALL_PREFIX=$prefix \
      -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain \
      ../
make
make install VERBOSE=1
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    MacOS(:x86_64),
    FreeBSD(:x86_64),
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Windows(:i686),
    Windows(:x86_64),
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libsymspg", :libsymspg)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
import Base: +, -, *, /, //, %, &, |, ^, <<, >>, ⊻

# multiply
Base.:*(a::T, b::Number) where {T<:AbstractJaxArray} =
    convert(T, PyObject(b) * a.o)
Base.:*(a::Number, b::T) where {T<:AbstractJaxArray} =
    convert(T,  b.o * PyObject(a))
Base.:*(a::AbstractJaxMatrix, b::AbstractJaxMatrix) =
    jax.numpy.matmul( b.o, a.o)
Base.:*(a::AbstractJaxMatrix, b::AbstractJaxVector) =
    jax.numpy.dot( b.o, a.o)
Base.:*(a::AbstractJaxVector, b::AbstractJaxMatrix) =
    jax.numpy.dot( b.o, a.o)

const overridenbfuncs = :[+, -, *, /].args

# Julia base dispatches vec + vec to broadcast.
# But since we override the broadcast machinery to un-broadcast everything,
# this will result in an endless loop.
# So we override this and dispatch vec + vec to vec _+ vec which when broadcasted
# will call the correct jax function.
for op in overridenbfuncs
    isdefined(Base, op) || continue
    fname = Symbol("_", op)

    @eval begin
        $fname(a::T, b::T) where {T<:AbstractJaxArray} = convert(T, $op(a.o, b.o))
        $fname(a::T, b::Number) where {T<:AbstractJaxArray} = convert(PyAny, $op(a.o, (b)))
        $fname(a::Number, b::T) where {T<:AbstractJaxArray} = convert(PyAny, $op((a), b.o))

        # If a tracer is applied to an JaxArray (a constant) then return a tracer.
        $fname(a::T, b::JaxArray) where {T<:AbstractJaxArray} = convert(T, $op(a.o, b.o))
        $fname(a::JaxArray, b::T) where {T<:AbstractJaxArray} = convert(T, $op(a.o, b.o))
        $fname(a::JaxArray, b::JaxArray) = convert(JaxArray, $op(a.o, b.o))

        $fname(a::AbstractJaxArray, b::AbstractJaxArray) = convert(PyAny, $op(a.o, b.o))

        # for many arguments like *(1,2,3)
        ($fname)(a, b, c, xs...) = Base.afoldl($fname, ($fname)(($fname)(a,b),c), xs...)

        jaxfunc(::typeof(Base.$op)) = $fname

        PyCall.pycall(::typeof($fname), ::Type, arg1, arg2) = $fname(arg1, arg2)
    end
end

@eval _pmul = $(Symbol("_", :*))

@eval $(Symbol("_", :-))(a::T) where {T<:AbstractJaxArray} = convert(T, -(a.o))
@eval $(Symbol("_", :+))(a::T) where {T<:AbstractJaxArray} = convert(T, +(a.o))


text:
function generate_logo()
    # Load original image (from https://github.com/JuliaImages/juliaimages.github.io/blob/source/docs/src/assets/logo.png)
    img = load(joinpath(@__DIR__, "logo_original.png"))

    # Define some colors
    c_red, c_green, c_blue, c_purple = RGBA.([Colors.JULIA_LOGO_COLORS...])
    gray_dark = RGBA{N0f8}(0.1,0.1,0.1)
    gray_light = RGBA{N0f8}(0.9,0.9,0.9)

    # Change base color black to gray
    img = img*N0f8(0.5) .+ RGBA{N0f8}(0.5,0.5,0.5,0)
    img_dark = copy(img)
    img_light = copy(img)

    # Define points on the image
    p_green = Point(87,35)
    p_red = Point(75,55)
    p_purple = Point(99,55)
    p_1 = Point(13,96)
    p_2 = Point(52,50)
    p_3 = Point(73,73)
    p_4 = Point(84,62)
    p_5 = Point(116,97)

    # Define path
    path_dark = joinpath(@__DIR__, "src", "assets", "logo-dark.png")
    path_light = joinpath(@__DIR__, "src", "assets", "logo.png")

    for (img,gray,name) in ((img_dark,gray_light,path_dark), (img_light,gray_dark,path_light))
        draw!(img, LineTwoPoints(p_1,p_2), gray)
        draw!(img, LineTwoPoints(p_2,p_3), gray)
        draw!(img, LineTwoPoints(p_3,p_4), gray)
        draw!(img, LineTwoPoints(p_4,p_5), gray)

        draw!(img, Cross(p_1,8), gray)
        draw!(img, Cross(p_2,8), gray)
        draw!(img, Cross(p_3,8), gray)
        draw!(img, Cross(p_4,8), gray)
        draw!(img, Cross(p_5,8), gray)

        draw!(img, Ellipse(CirclePointRadius(p_green,10)), c_green)
        draw!(img, Ellipse(CirclePointRadius(p_red,10)), c_red)
        draw!(img, Ellipse(CirclePointRadius(p_purple,10)), c_purple)

        save(name,img)
    end
end


text:
struct OmniscapeFlags
    calc_flow_potential::Bool
    calc_normalized_current::Bool
    compute_flow_potential::Bool #bad name, but need a variable that stores the "or" of the previous two, internal use only
    write_raw_currmap::Bool
    parallelize::Bool
    correct_artifacts::Bool
    source_from_resistance::Bool
    conditional::Bool
    mask_nodata ::Bool
    resistance_is_conductance::Bool
    write_as_tif::Bool
    allow_different_projections::Bool
    reclassify::Bool
    write_reclassified_resistance::Bool
end

struct Target
    x_coord::Int64
    y_coord::Int64
    amps::Float64
end

struct Conditions
    comparison1::String
    comparison2::String
    condition1_lower::Number
    condition1_upper::Number
    condition2_lower::Number
    condition2_upper::Number
end

struct ConditionLayers
    condition1_present::Array{Union{Missing, Number}, 2} # where T <: Number
    condition1_future::Array{Union{Missing, Number}, 2} # where U <: Number
    condition2_present::Array{Union{Missing, Number}, 2} # where V <: Number
    condition2_future::Array{Union{Missing, Number}, 2} # where W <: Number
end


text:
module RcTemperature
using Parameters
export ExpTemp, ExpTempResults

@with_kw struct ExpTemp
    notebook_pdf::String # path to pdf-scan of notebook containing the measurements
end

struct ExpTempResults
end

end


text:
module ICNBenchmarks

# usings
using BenchmarkTools
using CompositionalNetworks
using ConstraintDomains
using Constraints
using CSV
using DataFrames
using DataVoyager
using Dictionaries
using Distributed
using DrWatson
using JSON
using Statistics
using StatsBase
using Tables
#using SharedArrays

# imports
import Constraints: make_error

# constants
export ALL_PARAMETERS
export BENCHED_CONSTRAINTS

# others
export analyse_composition
export analyze_icn
export compositions_benchmark
export icn_benchmark
export search_space
export visualize_compositions
export visualize_icn

# includes
include("constants.jl")
include("search_space.jl")
include("extra_constraints.jl")
include("icn.jl")
include("composition.jl")
include("analyze.jl")

end


text:
function Balances(dx,x,problem,t)

    # alias -
    mRNA = x[1]
    G = x[2]
    σ70 = x[3]

    # calculate kinetics -
    transcription_kinetics = calculate_transcription_kinetics(t,x,problem)
    degradation_kinetics = calculate_mRNA_degradation_kinetics(t,x,problem)

    # compute the control -
    u_variable = calculate_transcriptional_control_array(t,x,problem)

    # compute the material balances -
    dx[1] = transcription_kinetics*u_variable - degradation_kinetics    # mRNA balances
    dx[2] = 0   # G
    dx[3] = 0   # σ70
end


text:
module mt008
using Random
using StaticArrays
using MeshCore: VecAttrib
using MeshCore: P1, L2, T3, ShapeColl, manifdim, nfacets, facetdesc, nshapes
using MeshCore: Q4ShapeDesc, shapedesc, n1storderv, nridges, nshifts, nvertices
using MeshCore: IncRel, attribute
using MeshCore: ir_skeleton, ir_transpose, ir_bbyfacets
using MeshSteward: T3block
using MeshSteward: vtkwrite
using MeshSteward: import_ABAQUS, vtkwrite, export_MESH, import_MESH
using ShellStructureTopo: orient_surface_mesh, _normal, make_topo_faces
using Test
function test()
    connectivities = import_ABAQUS(joinpath("../models", "extrusion.inp"))
    @test length(connectivities) == 1
    t2v = connectivities[1]
    for k in 1:nshapes(t2v.left)
        t2v._v[k] = shuffle!(Vector(t2v._v[k]))
    end
    e2v = ir_skeleton(t2v)
    # println(summary(e2v))
    t2e = ir_bbyfacets(t2v, e2v)
    t2v, orientable = orient_surface_mesh(t2v)
    @test orientable == true
    @test length(unique(attribute(t2v.left, "surfid"))) == 3
    geom = t2v.right.attributes["geom"]
    t2v.right.attributes["normals"] = VecAttrib([_normal(t2v[i], geom) for i in 1:nshapes(t2v.left)])
    vtkwrite("mt008-normals", t2v, [(name = "normals", )])
    # try rm("mt008-consistent.vtu"); catch end
    t2v = make_topo_faces(t2v)
    vtkwrite("mt008-classified", t2v, [(name = "tfid", )])
    true
end
end
using .mt008
mt008.test()


text:
module TrajOpt

    using FixedSizeArrays
    using Polynomials

    export HermiteBasis, QuadratureRule, BCTypes
    export compute_interpolation_matrices
    export assemble_stiffness_matrix
    export apply_bc!

    immutable Block{T,N} <: MutableFixedMatrix{T,N,N}

    end

    """
        hb = HermiteBasis(degree=N)

    Return left and right Hermite interpolation polynomials of specified degree.

    # Description
    Hermite basis function are used to interpolate a function over the interval
    τ = [0, 1] when the value of the function and it's first k derviatives are
    known at endpoints of the interval. The basis functions are simple τ-poly-
    nomials with degree N = 2k-1. To compute the value of the interpolant at a
    point within the interval, the basis functions are used to construct a
    weighted sum of the function data as follows (note: this is pseudo code):

        x_left  = [ x(0.0), x'(0.0) ... ]
        x_right = [ x(1.0), x'(1.0) ... ]
        x(τ) = sum( h.left(τ).*x_left ) + sum( h.right(τ).*x_right )

    # Fields
    * left::Vector{Poly}   Basis functions for left-endpoint function data
    * right::Vector{Poly}  Basis functions for right-endpoint function data
    """
    immutable HermiteBasis

        left  :: Vector{Poly{Rational{Int64}}}
        right :: Vector{Poly{Rational{Int64}}}

        function HermiteBasis(; degree=3)
            if degree == 1
                left  = [ Poly([ 1,  -1 ]//1, :τ) ]
                right = [ Poly([ 0,   1 ]//1, :τ) ]
            elseif degree == 3
                left  = [
                    Poly([ 1,   0,  -3,   2 ]//1, :τ)
                    Poly([ 0,   1,  -2,   1 ]//1, :τ)
                ]
                right = [
                    Poly([ 0,   0,   3,  -2 ]//1, :τ)
                    Poly([ 0,   0,  -1,   1 ]//1, :τ)
                ]
            elseif degree == 5
                left  = [
                    Poly([ 1,   0,   0, -10,  15,  -6 ]//1, :τ)
                    Poly([ 0,   1,   0,  -6,   8,  -3 ]//1, :τ)
                    Poly([ 0,   0,   1,  -3,   3,  -1 ]//2, :τ)
                ]
                right = [
                    Poly([ 0,   0,   0,  10, -15,   6 ]//1, :τ)
                    Poly([ 0,   0,   0,  -4,   7,  -3 ]//1, :τ)
                    Poly([ 0,   0,   0,   1,  -2,   1 ]//2, :τ)
                ]
            else
                throw(ArgumentError("degree=$degree: must be 1, 3, or 5."))
            end
            return new(left, right)
        end

    end

    """
        qr = QuadratureRule(num_points=N)

    Returns the N-point Guassian quadrature rule.

    # Description
    Numerical quadrature approximates the integral of a continuous function
    by sampling the integrand at discrete points in the domiain of integration
    and then taking a weighted sum of the sample values. Guassian quadrature
    rules are optimal for a 1D domain of integration: an N-point rule can
    exactly integrate a polynomial of degree 2*N-1. This type returns the
    sample locations and weights for a Gaussian Quadrature Rule with a
    specified number of sample locations.

    # Fields
    * locations::Vector{Real}  Quad point locations on bi-unit interval
    * weights::Vector{Real}    Summation weights for function samples

    """
    immutable QuadratureRule

        locations :: Vector{Float64}
        weights   :: Vector{Float64}

        function QuadratureRule(; num_points=2)
            if num_points == 1
                locations = [  0.0 ]
                weights   = [  2.0 ]
            elseif num_points == 2
                locations = [ -1.0,  1.0 ] * sqrt(1/3)
                weights   = [  1.0,  1.0 ]
            elseif num_points == 3
                locations = [ -1.0,  0.0,  1.0 ] * sqrt(3/5)
                weights   = [  5.0,  8.0,  5.0 ] / 9.0
            else
                throw(ArgumentError(
                    "num_points=$num_points: must be integer between 1 and 3."
                ))
            end
            return new(locations, weights)
        end
    end

    function compute_interpolation_matrices(hb, qr, time_points)

        # Preliminaries
        n_interval = length(time_points)-1
        n_quad = length(qr.locations)
        n_dof = length(hb.left)

        # Allocate interpolation matrices
        zero_block = MutableFixedMatrix{n_dof, n_quad, Float64}(0.0)
        HL = repmat([zero_block], n_interval)
        HR = repmat([zero_block], n_interval)

        # Map quadrature locations to unit interval
        τ = (qr.locations + 1)/2  # Maps [-1,1] => [0,1]
        τ_scaling = 1/2

        # Compute unit-interval interp matrices
        HL_unit = τ_scaling * polyval.( hb.left,  τ' )
        HR_unit = τ_scaling * polyval.( hb.right, τ' )

        # Apply interval scaling to interp matrices
        Δt = diff(time_points)
        scale_factors(Δt) = [ Δt^(q+1) for q = 0:n_dof-1 ]
        for i = 1:n_interval
            S = scale_factors(Δt[i]);
            HL[i] = S .* HL_unit
            HR[i] = S .* HR_unit
        end

        return (HL, HR)

    end

    """
        K = assemble_stiffness_matrix(hb, time)

    Compute stiffness matrix for time integration of 2nd order ODEs.

    # Description
    The Continuous Galerkin method can be used to discretize and solve systems
    of differential equations of the form ``\ddot{x} = f(t,x,\dot{x})``. The
    discretization of this system using Hermite basis functions results in a
    system of simultaneous non-linear equations ``Kx + f(x) = 0``, where ``K``
    is the stiffness matrix. This matrix is a block tri-diagonal matrix that
    is independent of the solution, x, and depends only on the distribution of
    time points used to discretize the solution domain. See the Jupyter note-
    book 'Theory.ipynb' for details.

    # Inputs
    * hb::HermiteBasis  Basis functions used to represent solution
    * time::AbstractVector{Float64}  Time points discretizing solution domain

    # Outputs
    * K::Tridiagonal{FixedSizeArrays.Mat}  Block tridiagonal stiffness matrix
    """
    function assemble_stiffness_matrix(hb, time_points)

        # Preliminaries
        n_point = length(time_points)
        n_interval = n_point-1
        n_dof = length(hb.left)

        # Error checking
        @assert n_dof>1   "Must use HermiteBasis with order >= 3"
        @assert n_point>2 "Must be at least 3 time points"

        # Allocate stiffness matrix
        zero_block = MutableFixedMatrix{Float64, n_dof, n_dof}(0.0)
        K = Tridiagonal(
            repmat( [zero_block], n_point-1 ), # Lower Diagonal
            repmat( [zero_block], n_point   ), # Main  Diagonal
            repmat( [zero_block], n_point-1 ), # Upper Diagonal
        )

        # Compute unit-interval stiffness matrices
        dhb_left  = polyder(hb.left)
        dhb_right = polyder(hb.right)
        integrate(x) = polyval( polyint(x), 1.0 )
        KLL = [ integrate(bi*bj) for bi=dhb_left,  bj=dhb_left  ]
        KLR = [ integrate(bi*bj) for bi=dhb_left,  bj=dhb_right ]
        KRL = [ integrate(bi*bj) for bi=dhb_right, bj=dhb_left  ]
        KRR = [ integrate(bi*bj) for bi=dhb_right, bj=dhb_right ]

        # Assemble global matrix from scaled unit-interval matrices
        Δt = diff(time_points)
        scale_factors(Δt) = [ Δt^-(p*q+1) for p=0:n_dof-1, q=0:n_dof-1 ]
        for i = 1:n_interval
            S = scale_factors(Δt[i])
            K[ i , i ] += S.*KLL
            K[ i ,i+1] += S.*KLR
            K[i+1, i ] += S.*KRL
            K[i+1,i+1] += S.*KRR
        end

        # Add boundary operators
        BL = [ polyval(bi*bj, 0.0) for bi=hb.left,  bj=dhb_left  ]
        BR = [ polyval(bi*bj, 1.0) for bi=hb.right, bj=dhb_right ]
        K[ 1 , 1 ] += scale_factors(Δt[ 1 ]) .* BL
        K[end,end] -= scale_factors(Δt[end]) .* BR

        return K

    end

    """
        F = assemble_load_vector(HL, HR, accel)

    Constructs the discrete load vector (RHS) for a uniform acceleration.

    # Description
    This function is a specicialized version of the assemble_load_vector that
    only works for constant accelerations that are independent of the trajectory
    solution (e.g. a constant gravitational acceleration). This is mainly used
    to verify correctness and accuracy of the Continuous Galerkin stiffness
    matrix and the associated boundary conditions.

    # Inputs
    HL::Vector{Mat}  Left-hand Hermite interpolation matrices
    HR::Vector{Mat}  Right-hand Hermite interpolation matrices
    accel::Float64   Constant acceleration applied to trajectory

    # Outputs
    F::Vector{Vec}   Load vector for discrete system of equations
    """
    function assemble_load_vector(HL, HR, weights, accel)

        # Preliminaries
        n_interval = length(HL)
        n_point = n_interval+1
        n_dof = size(HL[1],1)

        # Error checking
        @assert n_dof>1   "Must use HermiteBasis with order >= 3"
        @assert n_point>2 "Must be at least 3 time points"

        # Allocate load vector
        F = repmat(MutableFixedVector{Float64, n_dof}(0.0), n_point)

        # Evaluate load integrals
        aw = accel * weights
        for i = 1:n_interval
            F[ i ] += HL[i] * aw
            F[i+1] += HR[i] * aw
        end

    end

    """
    Defines named constants for different boundary condition types
    """
    module BCTypes
        @enum BCType dirichlet=1 neumann=2
    end

    """
        apply_bc!(K, bc, i)

    Modifies stiffness matrix to enforce Dirichlet/Neumman boundary conditions.

    # Description:
    Modifies stiffness matrix to enforce Dirichlet/Neumman boundary conditions
    at a specified time-point in the solution. Note that this doesn't actually
    define what the *value* at the boundary will be; the value is specified by
    calling apply_bc! with the load vector as the first argument, e.g.

        apply_bc!(f, bc, i, value)

    This version of apply_bc! simply modifies the stiffness matrix so that the
    required degree of freedom is set equal to the value specified in the load
    vector.
    """
    function apply_bc!{T <: MutableFixedMatrix}(K::Tridiagonal{T}, bc, i)
        # How this works:
        # To embed a boundary condition in the K matrix, we need to replace the
        # existing discrete equation for the "p"-th degree of freedom such that
        #   (K*x)[p] = x[p] = f[p]
        # Since K block tri-diagonal, this boils down to zero-ing out the
        # "q"-th row in of the three blocks coupled to the i-th node and then
        # setting the diagonal element to unity.
        #
        # Since we are using Hermite shape functions, the value of "q" depends
        # on the type of BC. If we want a Dirichlet BC we modify the 1st nodal
        # degree of freedom; if we want a Neumann BC, we modifiy the 2nd nodal
        # degree of freedom. This mapping from BC type to row index is codified
        # using the BCType enumeration.
        N = size(K,2)
        ibc = Int(bc)
        K[i,i][ibc,:]  = 0.0
        K[i,i][ibc,ibc] = 1.0
        if i > 1;  K[i,i-1][ibc,:] = 0.0;  end
        if i < N;  K[i,i+1][ibc,:] = 0.0;  end
        return nothing
    end

    """
        apply_bc!(F, bc, i, value)

    Modifies load vector to apply a Dirichlet/Neumman boundary conditions.

    # Description:
    Modifies assembled load vector to apply a Dirichlet or Neumman boundary
    conditions at a specified time-point in the solution. Note that for this
    to be effective, the stiffness matrix must also be modified at the specified
    node to accept a boundary condition. This is achieved by calling apply_bc!
    with the stiffness matrix as the first argument.
    """
    function apply_bc!{T <: MutableFixedVector}(F::Vector{T}, bc, i, value)
        F[i][Int(bc)] = value
        return nothing
    end

end


text:
"""
    stability_2ph(eos, c, [K])

Determine if mixture is single-phase stable under conditions `c`.

This is done using a version of Michelsen's stability test.

Reference: [The isothermal flash problem. Part I. Stability](https://doi.org/10.1016/0378-3812(82)85001-2)
"""
function stability_2ph(eos, c, K = initial_guess_K(eos, c); kwarg...)
    storage = flash_storage(eos, c)
    stability_2ph!(storage, K, eos, c)
end

function stability_2ph!(storage, K, eos, c; verbose = false, kwarg...)
    forces = storage.forces
    f_z = storage.buffer1
    f_xy = storage.buffer2
    x, y = storage.x, storage.y
    z, p, T = c.z, c.p, c.T
    liquid = (p = p, T = T, z = x)
    vapor = (p = p, T = T, z = y)
    # Update fugacities for current conditions used in both tests
    mixture_fugacities!(f_z, eos, c, forces)
    props = eos.mixture.properties

    wilson_estimate!(K, props, p, T)
    stable_vapor, i_v = michelsen_test!(vapor, f_z, f_xy, vapor.z, z, K, eos, c, forces, Val(true); kwarg...)
    wilson_estimate!(K, props, p, T)
    stable_liquid, i_l = michelsen_test!(liquid, f_z, f_xy, liquid.z, z, K, eos, c, forces, Val(false); kwarg...)
    stable = stable_vapor && stable_liquid
    if !stable
        @. K = y/x
    end
    if verbose
        @info "Stability done. Iterations:\nV: $i_v\nL: $i_l" stable_vapor stable_liquid stable
    end
    return stable
end

f_ratio(f_z, f_xy, ::Val{true}) = f_z/f_xy
f_ratio(f_z, f_xy, ::Val{false}) = f_xy/f_z
xy_value(z, K, ::Val{true}) = z*K
xy_value(z, K, ::Val{false}) = z/K

"""
    stability_2ph!(storage, eos, c, [K])

In-place version of [`stability_2ph`](@ref). `storage` should be allocated by `flash_storage`.
"""
function michelsen_test!(c_inside, f_z, f_xy, xy, z, K, eos, cond, forces, inside_is_vapor;
    tol_equil = 1e-10, tol_trivial = 1e-5, maxiter = 1000)
    trivial = false
    S = 1.0
    iter = 0
    done = false
    while !done
        iter += 1
        S = 0.0
        @inbounds for c in eachindex(xy)
            xy_i = xy_value(z[c], K[c], inside_is_vapor)
            xy[c] = xy_i
            S += xy_i
        end
        @. xy /= S
        mixture_fugacities!(f_xy, eos, c_inside, forces)

        R_norm = 0.0
        K_norm = 0.0
        @inbounds for c in eachindex(K)
            R = f_ratio(f_z[c], S*f_xy[c], inside_is_vapor)
            K[c] *= R

            R_norm += (R-1)^2
            K_norm += log(K[c])^2
        end
        # Two convergence criteria:
        # - Approaching trivial solution (K-values are all 1)
        # - Equilibrium for a small amount of the "other" phase,
        #   the single-phase conditions are not stable.
        trivial = K_norm < tol_trivial
        converged = R_norm < tol_equil

        # Termination of loop
        ok = trivial || converged
        done = ok || iter == maxiter
        if done && !ok
            trivial = true
            @warn "Stability test failed to converge in $maxiter iterations. Assuming stability." cond xy K_norm R_norm K
        end
    end
    stable = trivial || S <= 1 + tol_trivial
    return (stable, iter)
end


text:
abstract type AbstractReducer{T} end
struct Reducer{T} <: AbstractReducer{T} end
struct ReducerWrapper{T,V} <: AbstractReducer{T}
    data::V
end
@inline ReducerWrapper{T}(data::V) where {T,V} = ReducerWrapper{T,V}(data)

Base.size(::AbstractReducer) = ()

@inline Base.:*(::Reducer{:row}, A::LinearAlgebra.Diagonal{<:Real,<:AbstractFixedSizeVector}) = A.diag'
@inline Base.:*(::Reducer{:row}, A::Number) = A
@inline Base.:*(::Reducer{:row}, A::LinearAlgebra.Adjoint{<:Any,<:AbstractFixedSizeVector}) = A

@generated function Base.:*(
    sp::StackPointer,
    A::PaddedMatrices.AbstractFixedSizeMatrix{M,N,T,P},
    ::Reducer{:row}
) where {M,N,T,P}
    quote
        reduction = PtrVector{$N,$T,$N}(pointer(sp,$T))
        @inbounds for n ∈ 0:(N-1)
            sₙ = zero($T)
            @vvectorize $T for m ∈ 1:$M
                sₙ += A[m + $P*n]
            end
            reduction[n+1] = sₙ
        end
        sp + $(VectorizationBase.align(N*sizeof(T))), reduction'
    end
end
@generated function RESERVED_INCREMENT_SEED_RESERVED!(
    C::AbstractMutableFixedSizeVector{N,T,PC},
    ::Reducer{:row},
    A::AbstractFixedSizeMatrix{M,N,T,PA}
) where {M,N,T,PA,PC}
    quote
        @inbounds for n ∈ 0:$(N-1)
            sₙ = zero(T)
            @vvectorize $T for m ∈ 1:$M
                sₙ += A[m + $PA*n]
            end
            C[n+1] = Base.FastMath.add_fast(C[n+1], sₙ)
        end
        nothing
    end
end
@generated function RESERVED_INCREMENT_SEED_RESERVED!(
    C::UninitializedVector{N,T,PC},
    ::Reducer{:row},
    A::AbstractFixedSizeMatrix{M,N,T,PA}
) where {M,N,T,PA,PC}
    quote
        @inbounds for n ∈ 0:$(N-1)
            sₙ = zero(T)
            @vvectorize $T for m ∈ 1:$M
                sₙ += A[m + $PA*n]
            end
            C[n+1] = sₙ
        end
        nothing
    end
end
@generated function Base.:*(sp::StackPointer, a::LinearAlgebra.Adjoint{T,<:AbstractVector{T}}, ::Reducer{S}) where {T,S}
    S == true && return quote
        $(Expr(:meta,:inline))
        (sp, sum(a.parent))
    end
    # @assert sum(S) == M
    N = length(S)
    q = quote
        out = PtrVector{$N,$T,$N}(pointer(sp,$T))
    end
    outtup = Expr(:tuple,)
    ind = 0
    for (j,s) ∈ enumerate(S)
        ind += 1
        accumulate_sym = gensym()
        push!(q.args, :($accumulate_sym = a[$ind]))
        for i ∈ 2:s
            ind += 1
            push!(q.args, :($accumulate_sym += a[$ind]))
        end
        push!(q.args, :(out[$j] = $accumulate_sym))
    end
    quote
        @fastmath @inbounds begin
            $q
        end
        sp + $(VectorizationBase.align(sizeof(T)*N)), out'
    end
end
@generated function Base.:*(sp::StackPointer, A::AbstractMatrix{T}, ::Reducer{S}) where {T,S}
    S == true && return :(sum(A))
    # y = Vector{T}(undef, 0)
    quote
        (sp, y) = PtrVector{$(sum(S)),$T}(sp)
        # resize!(y, size(A,1))
        *(sp, sum!(y, A)', Reducer{$S}())
    end
end
@generated function Base.:*(
    sp::StackPointer,
    a::LinearAlgebra.Adjoint{T, <: AbstractVector{T}},
    b::ReducerWrapper{S, <: AbstractVector{T}}
) where {T,S}
    # @assert sum(S) == M
    N = length(S)
    q = quote
        (sp, out) = PtrVector{$N,$T}(sp)
    end
    ind = 0
    ind2 = 0
    for (j,s) ∈ enumerate(S)
        ind += 1
        accumulate_sym = gensym()
        push!(q.args, :($accumulate_sym = a[$ind] * b[$ind]))
        for i ∈ 2:s
            ind += 1
            push!(q.args, :($accumulate_sym += a[$ind] * b[$ind]))
        end
        ind2 += 1
        push!(q.args, :(out[$ind2] = accumulate_sym))
    end
#    P = PaddedMatrices.calc_padding(N, T)
#    for p ∈ N+1:P
#        push!(outtup.args, zero(T))
#    end
    quote
        @fastmath @inbounds begin
            $q
 #           ConstantFixedSizeVector{$N,$T,$P}($outtup)'
        end
        sp, out'
    end
end
@generated function Base.:*(sp::StackPointer, A::AbstractMatrix{T}, rw::ReducerWrapper{S}) where {T,S}
    # y = Vector{T}(undef, 0)
    quote
        (sp, y) = PtrVector{$(sum(S)),$T}(sp)

        # resize!(y, size(A,1))
        *(sp, sum!($y, A)', rw)
    end
end

@generated function RESERVED_INCREMENT_SEED_RESERVED!(
    seedout::AbstractMutableFixedSizeVector{M,T,L},
    jac::UniformScaling{T},
    seedin::AbstractFixedSizeVector{M,T,L}
) where {M,T,L}
    op = seedout <: UninitializedVector ? :(=) : :(+=)
    quote
        $(Expr(:meta,:inline))
        j = jac.λ
        @vvectorize $T for m ∈ 1:$M
            $(Expr(op, :(seedout[m]), :(j * seedin[m])))
        end
    end
end

#@inline ∂mul(x, y, ::Val{(true,true)}) = ((@show typeof(x), typeof(y)); return y, x)
@inline ∂mul(x, y, ::Val{(true,true)}) = y', x'
@inline ∂mul(x, y, ::Val{(true,false)}) = y'
@inline ∂mul(x, y, ::Val{(false,true)}) = x'
@inline ∂mul(x::Real, y::Real, ::Val{(true,true)}) = y, x
@inline ∂mul(x::Real, y::Real, ::Val{(true,false)}) = y
@inline ∂mul(x::Real, y::Real, ::Val{(false,true)}) = x

@inline ∂mul(x::AbstractVector, y::Real, ::Val{(true,true)}) = y, Diagonal(x)
@inline ∂mul(x::AbstractVector, y::Real, ::Val{(true,false)}) = y
@inline ∂mul(x::AbstractVector, y::Real, ::Val{(false,true)}) = Diagonl(x)

@inline ∂mul(x::Real, y::AbstractVector, ::Val{(true,true)}) = y, UniformScaling(x)
@inline ∂mul(x::Real, y::AbstractVector, ::Val{(true,false)}) = y
@inline ∂mul(x::Real, y::AbstractVector, ::Val{(false,true)}) = UniformScaling(x)

@inline function ∂mul(
    D::LinearAlgebra.Diagonal{T,<:AbstractFixedSizeVector{M,T,P}},
    L::StructuredMatrices.AbstractLowerTriangularMatrix{M,T,N}, ::Val{(true,true)}
) where {M,N,P,T}
    StructuredMatrices.∂DiagLowerTri∂Diag(L), StructuredMatrices.∂DiagLowerTri∂LowerTri(D)
end
@inline function ∂mul(
    D::LinearAlgebra.Diagonal{T,<:AbstractFixedSizeVector{M,T,P}},
    L::StructuredMatrices.AbstractLowerTriangularMatrix{M,T,N}, ::Val{(true,false)}
) where {M,N,P,T}
    StructuredMatrices.∂DiagLowerTri∂Diag(L)
end
@inline function ∂mul(
    D::LinearAlgebra.Diagonal{T,<:AbstractFixedSizeVector{M,T,P}},
    L::StructuredMatrices.AbstractLowerTriangularMatrix{M,T,N}, ::Val{(false,true)}
) where {M,N,P,T}
    StructuredMatrices.∂DiagLowerTri∂LowerTri(D)
end


text:
using .Turing, Random, Test

# Import utility functions and reused models.
include("staging.jl")
include("numerical_tests.jl")
include("ad_utils.jl")
include("models.jl")
include("random_measure_utils.jl")
include("testing_functions.jl")


text:
# Construct n-grams using single space concatenation
function ngramize(lang::S, words::Vector{T}, n::Int=DEFAULT_NGRAM_COMPLEXITY
                 ) where {S<:Language, T<:AbstractString}
    if n == 1
        return onegramize(lang, words)
    end
    n_words = length(words)
    tokens = Dict{T, Int}()
    @inbounds for m in 1:n
        for index in 1:(n_words - m + 1)
            token = join(words[index:(index + m - 1)], " ")
            tokens[token] = get(tokens, token, 0) + 1
        end
    end
    return tokens
end

function onegramize(lang::S, words::Vector{T}) where {S<:Language, T<:AbstractString}
    n_words = length(words)
    tokens = Dict{T, Int}()
    @inbounds for word in words
        tokens[word] = get(tokens, word, 0) + 1
    end
    return tokens
end


text:
module ConstLab

using Tensors, Parameters
import NLsolve, ForwardDiff, DiffResults

# Stress and strain control
abstract type ControlType end
struct StrainControl <: ControlType
    f # f(t) -> strain
end
struct StressControl <: ControlType
    # TODO: Make this just a MixedControl instead,
    # no gain in having a separate type
    f # f(t) -> stress
end
struct MixedControl <: ControlType
    # TODO: improve this masking by e.g. making masks 4th order tensors
    strain
    strain_mask::Vector{Bool}
    stress
    stress_mask::Vector{Bool}
    function MixedControl(strain, strain_mask, stress, stress_mask)
        strain_mask = collect(strain_mask.data)
        stress_mask = collect(stress_mask.data)
        for (a, b) in zip(strain_mask, stress_mask)
            a ⊻ b || error("no can do")
        end
        return new(strain, strain_mask, stress, stress_mask)
    end
end

###############################################
###############################################
###############################################
###############################################
abstract type Material end

abstract type MaterialState end

####################
## Linear Elastic ##
####################
struct LinearElastic <: Material
    G::Float64
    K::Float64
end
struct LinearElasticState{D,N} <: MaterialState
    ε::SymmetricTensor{2,D,Float64,N}
    σ::SymmetricTensor{2,D,Float64,N}
end
function LinearElasticState{D}() where D
    z = zero(SymmetricTensor{2,D})
    return LinearElasticState(z, z)
end

function constitutive_driver(model::LinearElastic, ε::SymmetricTensor, state::LinearElasticState)
    G, K = model.G, model.K
    σ = 2G * dev(ε) + 3K * vol(ε)
    return LinearElasticState(ε, σ)
end

#############################
## Plastic Prototype Model ##
#############################
struct Plastic <: Material
    G::Float64
    K::Float64
    σ_y::Float64 # Yield stress
    H::Float64
    κ∞::Float64
    α∞::Float64
    r::Float64
end
struct PlasticState{D,N} <: MaterialState
    ε::SymmetricTensor{2,D,Float64,N}
    σ::SymmetricTensor{2,D,Float64,N}
    εₚ::SymmetricTensor{2,D,Float64,N}
    κ::Float64
    α::SymmetricTensor{2,D,Float64,N}
    μ::Float64
end
# Initial state
function PlasticState{D}() where D
    z = zero(SymmetricTensor{2,D})
    return PlasticState(z, z, z, 0.0, z, 0.0)
end

function elastic_tangent(model::Plastic)
    G, K = model.G, model.K
    # Use Lamé parameters
    λ = K - 2G/3
    μ = G
    δ(i,j) = i == j ? 1.0 : 0.0
    E = SymmetricTensor{4,3}() do i, j, k, l
        λ * δ(i,j) * δ(k,l) + μ * (δ(i,k) * δ(j,l) + δ(i,l) * δ(j,k))
    end
    return E
end

mise(s) = √(3/2) * norm(dev(s))

function constitutive_driver(model::Plastic, ε::SymmetricTensor, state::PlasticState;
        solver_params=Dict{Symbol,Any}())
    G, K, σ_y = model.G, model.K, model.σ_y

    # Assume elastic step
    Δε = ε - state.ε
    σ_tr = state.σ + 2G * dev(Δε) + K * tr(Δε) * one(Δε)
    Φ = mise(σ_tr - state.α) - σ_y - state.κ
    if Φ < 0
        # elastic
        σ = σ_tr
        εₚ = state.εₚ
        κ = state.κ
        α = state.α
        μ = state.μ
        dσdε = elastic_tangent(model)
    else
        # plastic; need to iterate
        σ, εₚ, κ, α, μ, dσdε = solve_local_problem(model, Δε, state; solver_params=solver_params)
    end
    return σ, dσdε, PlasticState(ε, σ, εₚ, κ, α, μ)
end

function pack_variables(σ, κ, α, μ::T) where T
    n = length(σ.data) + length(κ) + length(α.data) + length(μ)
    @assert n == 14 # hehe
    X = Vector{T}(undef, n)
    return pack_variables!(X, σ, κ, α, μ)
end

function pack_variables!(X, σ, κ, α, μ)
    tomandel!(X, σ)
    X[7] = κ
    tomandel!(X, α; offset=6+1)
    X[end] = μ
    return X
end

function extract_variables(X)
    N = 6
    σ = frommandel(SymmetricTensor{2,3}, X)
    κ = X[N+1]
    α = frommandel(SymmetricTensor{2,3}, X; offset=N+1)
    μ = X[end]
    return σ, κ, α, μ
end

function solve_local_problem(model::Plastic, Δε::SymmetricTensor, state::PlasticState{D,N}; solver_params=Dict{Symbol,Any}()) where {D,N}
    # extract material parameters
    @unpack G, K, σ_y, κ∞, α∞, r, H = model

    σₙ = state.σ
    εₚₙ = state.εₚ
    κₙ = state.κ
    αₙ = state.α
    μₙ = state.μ * 0 # ??
    σ_tr = σₙ + 2G * dev(Δε) + K * tr(Δε) * one(Δε)

    # Initial guess
    X0 = pack_variables(σ_tr, κₙ, αₙ, μₙ)

    # Residual function
    function residual!(R, X)
        # fill!(R, 0) # pack_variables! should overwrite
        # Unpack variables
        σ, κ, α, μ = extract_variables(X)
        σ_red_dev = dev(σ - α)

        # Compute residuals
        Rσ = σ - σ_tr + 3G*μ/mise(σ_red_dev) * σ_red_dev
        Rκ = κ - κₙ - r * H * μ * (1 - κ/κ∞)
        Rα = α - αₙ - (1-r) * H * μ * (σ_red_dev / mise(σ_red_dev) - α/α∞)
        RΦ = mise(σ_red_dev) - σ_y - κ

        # Scale stresses with σ_y to make the residual dimensionless
        Rσ /= σ_y
        Rκ /= σ_y
        Rα /= σ_y
        RΦ /= σ_y

        # Populate residual vector R
        pack_variables!(R, Rσ, Rκ, Rα, RΦ)
        return R
    end

    # Use NLsolve to solve non-linear system
    params = Dict{Symbol,Any}(:ftol=>1e-14, :show_trace=>false, :method => :newton)
    merge!(params, solver_params)
    R = similar(X0) # residual cache for constructing OnceDifferentiable (not used :( )
    cache = NLsolve.OnceDifferentiable(residual!, X0, R; autodiff=:forward)
    res = NLsolve.nlsolve(cache, X0; params...)
    NLsolve.converged(res) || error("local problem did not converge: $res")
    X = res.zero # Solution to the problem
    J = cache.DF # Extract Jacobian from last step

    # # Custom implementation of Newton's algorithm
    # X = X0 # zeros(14)
    # R = zeros(14)
    # cfg = ForwardDiff.JacobianConfig(residual!, R, X)
    # out = DiffResults.JacobianResult(X)
    # iters = 0
    # local J
    # while true; iters += 1
    #     # residual!(R, X)
    #     ForwardDiff.jacobian!(out, residual!, R, X, cfg)
    #     g = DiffResults.value(out)
    #     J = DiffResults.jacobian(out)
    #     if norm(g) < 1e-6
    #         @info "local problem converged" iters norm(g)
    #         break
    #     elseif iters > 20
    #         error("did not converge")
    #     end
    #     ΔX = J \ g
    #     X .-= ΔX
    # end

    # Extract variables from the solution X
    σ, κ, α, μ = extract_variables(X)
    # Extract ATS tensor from final Jacobian
    dRdε = elastic_tangent(model) # For fixed X
    Jinv = inv(J)
    Jinv_t = frommandel(SymmetricTensor{4,3}, Jinv; offset_i=0, offset_j=0)
    dσdε = (Jinv_t / σ_y) ⊡ dRdε # scale Jacobian with σ_y since residual scaled with σ_y

    # Compute plastic strain; dεₚdt = λ ν -> εₚ = εₚₙ + Δt λ ν = εₚₙ + μ ν
    ν = 3/2 * dev(σ - α) / mise(dev(σ - α))
    εₚ = εₚₙ + μ * ν

    return σ, εₚ, κ, α, μ, dσdε
end

###############################################
###############################################
###############################################
###############################################

"""
    integrate(model, time, s0)

Integrate the material response.
"""
function integrate(model::Material, ctrl::StrainControl, time, state::T;
        solver_params=Dict{Symbol,Any}()) where {T <: MaterialState}
    states = T[]
    for (i, t) in enumerate(time)
        ε = ctrl.f(t)
        σ, dσdε, state′ = constitutive_driver(model, ε, state; solver_params=solver_params)
        push!(states, state′)
        state = state′
    end
    return states
end

function integrate(model::Material, ctrl::MixedControl, time, state::T;
        solver_params=Dict{Symbol,Any}()) where {T <: MaterialState}
    states = T[]
    for (i, t) in enumerate(time)
        # Compute controlled strain
        ε_c = ctrl.strain(t)
        # Compute controlled stress
        σ_c = ctrl.stress(t)

        # Iterate to find the correct strain
        local state′
        iters = 0
        Δε = zero(state.ε) # initial guess
        # println("----- strain iterations -----")
        while true; iters += 1
            # Current guess
            ε = SymmetricTensor{2,3}(ε_c.data .* ctrl.strain_mask .+ (state.ε.data .+ Δε.data) .* ctrl.stress_mask)

            σ′, dσdε, state′ = constitutive_driver(model, ε, state; solver_params=solver_params)

            # Compute stress residual
            R = tovoigt(σ′ - σ_c)[ctrl.stress_mask]
            # println("norm(R) = ", norm(R))
            if norm(R)/model.σ_y < 1e-6 # TOL
                # println("strain iterations converged in $iters iterations")
                break
            elseif iters > 20
                error("strain iterations did not converge")
            end
            # Take a step in strain
            J = tovoigt(dσdε)[ctrl.stress_mask, ctrl.stress_mask]
            ΔΔε = - J \ R
            ΔΔε_a = zeros(6)
            ΔΔε_a[ctrl.stress_mask] .= ΔΔε
            ΔΔε_t = fromvoigt(SymmetricTensor{2,3}, ΔΔε_a)
            Δε += ΔΔε_t
        end
        # Done, store stuff
        push!(states, state′)
        state = state′
    end
    return states
end

end # module


text:
module Io

import Base.write

using Avro.Common
using Avro.Schemas

export Encoder,
       Decoder,
       BinaryEncoder,
       BinaryDecoder,
       encode_null,
       encode_boolean,
       encode_int,
       encode_long,
       encode_float,
       encode_double,
       encode_byte,
       encode_bytes,
       encode_fixed,
       encode_string,
       decode_null,
       decode_boolean,
       decode_int,
       decode_long,
       decode_float,
       decode_double,
       decode_byte,
       decode_bytes,
       decode_fixed,
       decode_string

abstract type Encoder end;
abstract type Decoder end;

"""
Writes binary data of builtin and primtive types to an output stream.
"""
struct BinaryEncoder <: Encoder
    stream::IO
end

"""
Reads data of simple and primitive types from an input stream.
"""
struct BinaryDecoder <: Decoder
    stream::IO
end

# Encoders

encode_null(::BinaryEncoder, ::Nothing) = 0
encode_boolean(encoder::BinaryEncoder, value::Bool) = write(encoder.stream, value)

encode_int(encoder::BinaryEncoder, value::Int32) = _encode_varint(encoder.stream, _encode_zigzag(value))
encode_long(encoder::BinaryEncoder, value::Int64) = _encode_varint(encoder.stream, _encode_zigzag(value))

function encode_float(encoder::BinaryEncoder, value::Float32)
    isnan(value) ? write(encoder.stream, NaN32) : write(encoder.stream, value)
end

function encode_double(encoder::BinaryEncoder, value::Float64)
    isnan(value) ? write(encoder.stream, NaN64) : write(encoder.stream, value)
end

encode_byte(encoder::BinaryEncoder, value::UInt8) = write(encoder.stream, value)

function encode_bytes(encoder::BinaryEncoder, value::Vector{UInt8})
    encode_long(encoder, length(value))
    write(encoder.stream, value)
end

encode_fixed(encoder::BinaryEncoder, value::Vector{UInt8}) = write(encoder.stream, value)

function encode_string(encoder::BinaryEncoder, value::String)
    encode_long(encoder, sizeof(value)) + write(encoder.stream, value)
end

function _encode_zigzag(n::T) where T <: Integer
	num_bits = sizeof(T) * 8
	(n << 1) ⊻ (n >> (num_bits - 1))
end

function _encode_varint(stream::IO, n::T) where T <: Integer
    max_bytes = sizeof(T) + ceil(Int, sizeof(T) / 8)
    bytes_written = 0
    while n > 0x7f && bytes_written < max_bytes
        bytes_written += write(stream, ((n | 0x80) & 0xff) % UInt8)
        n >>>= 7
    end
    if n > 0x7f
        throw(Exception("Invalid variable-length integer encoding"))
    end
    bytes_written += write(stream, n % UInt8)
end

# Decoders

decode_null(::Decoder) = nothing
decode_boolean(decoder::Decoder) = read(decoder.stream, Bool)

decode_int(decoder::Decoder) = _decode_zigzag(_decode_varint(decoder.stream, Int32))
decode_long(decoder::Decoder) = _decode_zigzag(_decode_varint(decoder.stream, Int64))

decode_float(decoder::Decoder) = read(decoder.stream, Float32)
decode_double(decoder::Decoder) = read(decoder.stream, Float64)

decode_byte(decoder::Decoder) = read(decoder.stream, UInt8)

function decode_bytes(decoder::Decoder)
    nb = decode_long(decoder)
    read(decoder.stream, nb)
end

decode_fixed(decoder::Decoder, nb::Integer) = read(decoder.stream, nb)

function decode_string(decoder::Decoder)
    nb = decode_long(decoder)
    String(decode_fixed(decoder, nb))
end

_decode_zigzag(n::Integer) = (n >>> 1) ⊻ -(n & 1)

function _decode_varint(stream::IO, ::Type{T}) where T <: Integer
    max_bytes = sizeof(T) + ceil(Int, sizeof(T) / 8)
    b = read(stream, UInt8) % T
    n = b & 0x7f
    bytes_read = 1
    while b > 0x7f && bytes_read < max_bytes
        b = read(stream, UInt8) % T
        n ⊻= (b & 0x7f) << (7 * bytes_read)
        bytes_read += 1
    end
    if b > 0x7f
        throw(Exception("Invalid variable-length integer encoding"))
    end
    n
end

end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

if VERSION >= v"0.5-dev+7720"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

export @test, @testset, @test_throws


text:
"""
    freeze_condition(u,t,integ)

Take the current step u, the current time step t, the previous step integ.uprev
and the previous time step integ.tprev, and compute the spatial and temporal difference.
If the spatial difference is less than 1Å, or if the temporal difference
is less than one femtosecond, return true. Otherwise return false.

Inputs:
    - u = the next step suggested by the integrator. u is an array containing the elements
    specified in GuidingCenterOrbits.jl/orbit.jl/integrate(). The elements of u fully characterize
    the guiding-centre particle. The elements of u are [r, phi, z, p_para, mu] where r,phi,z are the
    cylindrical coordinates, p_para is the momentum parallel to the magnetic field and mu is the
    magnetic moment.
    - t = the current time step of the integrator
    - integ = the integrator. Read more at the 'Integrator interface' page of the DifferentialEquations.jl package website.
"""
function freeze_condition(u,t,integ)
    Δu = abs.(u .- integ.uprev) # The difference between the current and previous time-step (cocos_vec,p_para,mu)
    Δspace = norm(Δu)
    #println(Δspace)
    Δtime = abs(t - integ.tprev)
    #println(Δtime)
    (Δspace < 1.0e-10) ||  (Δtime < 1.0e-15) # If the integration has stagnated in space or time (less than an Ångström or femtosecond)...
end
function freeze_affect!(integ)
    #println("I did it!")
    terminate!(integ) # If the freeze_condition returs true, terminate the integration!
end
# Read more about Callbacks at the DifferentialEquations.jl 'Event handling and Callback Functions' page
brr_cb = DiscreteCallback(freeze_condition,freeze_affect!,save_positions=(false,false)) # save_positions = (false,false). If freeze_condition, do not save position before or after suggested time step.

"""
    phi_condition(u,t,integ)

Check whether the guiding-centre particle has gone around the tokamak toroidally more than
maxphi in phi angle, without completing even one poloidal orbit. If so, return true.

Inputs:
    Same as for freeze_condition.
"""
function phi_condition(maxphi, u, t, integ)
    stat = integ.f.f.stat # The GCStatus mutable struct defined in orbit.jl
    ir, iphi, iz = cylindrical_cocos_indices(cocos(integ.f.f.M))
    !stat.poloidal_complete && (abs(u[iphi]-stat.ri[iphi]) > maxphi) # If it has gone around the tokamak more than (maxphi/(2*pi)) times toroidally, but not once poloidally...
end
function phi_affect!(integ)
    integ.p && terminate!(integ) # If only one poloidal transit desired, terminate the integration
end
phi_callback(maxphi) = DiscreteCallback((u,t,integ)->phi_condition(maxphi,u,t,integ),phi_affect!,save_positions=(false,false))
# Check callback syntax at https://diffeq.sciml.ai/stable/features/callback_functions/ for further info

"""
    r_condition(u,t,integ)

Return the ratio of the r component of the velocity, to the total speed of the particle.

Inputs:
    Same as freeze_condition.
"""
function r_condition(u,t,integ)
    v = integ.f(u,integ.p,t)
    ir, iphi, iz = cylindrical_cocos_indices(cocos(integ.f.f.M))
    v[ir]/sqrt(v[ir]^2 +v[iz]^2)
end
function r_affect!(integ)
    stat = integ.f.f.stat # Get the stat struct. Please see GuidingCenterOrbits.jl/orbit.jl/GCStatus{} for fields.
    ir, iphi, iz = cylindrical_cocos_indices(cocos(integ.f.f.M))
    if !stat.poloidal_complete && stat.nr < 20 # prevent infinite loop
        stat.nr += 1
        if (integ.u[ir] > stat.rm) # If the guiding-centre particles has gone outside its previous rm coordinate
            stat.rm = integ.u[ir] # Update the rm coordinate with the current r coordinate
            stat.zm = integ.u[iz] # Update the zm coordinate with the current z coordinate
            stat.tm = integ.t # Update the poloidal transit time with the current time
            M = integ.f.f.M # Get the axisymmetric equilibrium
            gcp = integ.f.f.gcp # Get the guiding-centre particle itself
            stat.pm = get_pitch(M,gcp,integ.u[4],integ.u[5],integ.u[ir],integ.u[iz]) # Update the pm coordinate with the current pitch
        end
    else
        integ.p && terminate!(integ)
    end
end
r_cb = ContinuousCallback(r_condition,r_affect!,abstol=1e-6)

"""
    poloidal_condition(u,t,integ)

If it's not the first time step of the integration, return the difference between
the current speed and the initial speed, both evaluated in the initial velocity
direction of the guiding-centre particle.

Inputs:
    Same as freeze_condition.
"""
function poloidal_condition(u,t,integ)
    stat = integ.f.f.stat
    i = stat.initial_dir
    (t != 0.0)*(u[i] - stat.ri[i])
end
function poloidal_affect!(integ)
    stat = integ.f.f.stat # Get the stat struct. Please see GuidingCenterOrbits.jl/orbit.jl/GCStatus{} for fields.
    ir, iphi, iz = cylindrical_cocos_indices(cocos(integ.f.f.M))
    vi = stat.vi # Get the initial velocity vector
    ri = stat.ri # Get the initial guiding-centre elements vector
    vc = integ.f(integ.u,integ.p,integ.t) # Get the current velocity vector
    dp = dot(vi,vc)/(norm(vi)*norm(vc)) # Compute the normalized dot product between the current velocity vector and the initial velocity vector
    vi_rz = SVector(vi[ir],vi[iz]) # Create a 2D vector which is the initial velocity vector projected onto the cross-section of the tokamak
    vc_rz = SVector(vc[ir],vc[iz]) # Create a 2D vector which is the current velocity vector projected onto the cross-section of the tokamak
    dprz = dot(vi_rz,vc_rz)/(norm(vi_rz)*norm(vc_rz)) # Compute the normalized dot product between the 2D projected initial and current velocities
    if !stat.poloidal_complete && (stat.nr >= 2) &&
        (dp > 0.99 && dprz > 0.99) && # If the dot products are almost 1...
        (abs(integ.u[ir]-ri[ir]) < 0.01) # ... and the particle is within 1 cm of the initial r coordinate

        stat.poloidal_complete=true # Deem the particle as poloidally complete
        stat.tau_p = integ.t # Set the current time to be the poloidal transit time
        stat.tau_t = 2pi*stat.tau_p/abs(integ.u[iphi] - integ.sol.u[1][iphi]) # Compute the toroidal transit time
        stat.class = :unknown # We don't yet know what class the orbit is (potato, banana... etc)
        integ.p && terminate!(integ) # Terminate the integration, if only one poloidal transit is desired
    end
end
pol_cb = ContinuousCallback(poloidal_condition,poloidal_affect!,abstol=1e-6)

"""
    out_of_bounds_condition(u,t,integ)

Check whether the guiding-centre particle is still inside the boundaries provided
by the axisymmetric equilibrium M. If that is not the case, return true.

Inputs:
    Same as freeze_condition.
"""
function out_of_bounds_condition(u,t,integ)
    M = integ.f.f.M
    ir, iphi, iz = cylindrical_cocos_indices(cocos(M))
    rlims, zlims = limits(M)
    !((rlims[1] < u[ir] < rlims[end]) && (zlims[1] < u[iz] < zlims[end]))
end
function out_of_bounds_affect!(integ)
    integ.f.f.stat.hits_boundary=true # The particle hit the boundary
    integ.f.f.stat.class = :lost # Classify it as lost
    terminate!(integ) # Terminate the integration
end
oob_cb = DiscreteCallback(out_of_bounds_condition, out_of_bounds_affect!,save_positions=(false,false))

"""
    wall_condition(u,t,integ)

Check whether the guiding-centre particle has hit the tokamak wall. If that is the case,
return true.

Inputs:
    Same as freeze_condition.
"""
function wall_condition(wall, u, t, integ)
    !in_vessel(wall,(u[1],u[3]))
end
function wall_affect!(integ)
    integ.f.f.stat.hits_boundary=true
    integ.f.f.stat.class = :lost
    terminate!(integ)
end
wall_callback(wall) = DiscreteCallback((u,t,integ)->wall_condition(wall,u,t,integ),wall_affect!,save_positions=(false,false))

transit_callback = CallbackSet(r_cb, pol_cb, oob_cb, brr_cb) # A collection of callbacks


text:
using Revise
using Distributed
addprocs(3)
@everywhere using Distributed, Jets, JournaledJets, LinearAlgebra, Test

@everywhere foo(iblock) = prod(iblock.I) * ones(2)
@everywhere bar(iblock) = LinearIndices((3,4))[iblock[1],iblock[2]] * ones(2,3)

@everywhere JopBar_f!(d,m;x) = d .= x[1] .* m.^2
@everywhere JopBar_df!(δd,δm;mₒ,x,kwargs...) = δd .= 2 .* x[1] .* mₒ .* δm
@everywhere function JopBar(n,x=[1.0]) spc = JetSpace(Float64, n)
    JopNl(f! = JopBar_f!, df! = JopBar_df!, df′! = JopBar_df!, dom = spc, rng = spc, s=(x=x,))
end

@testset "JArray construct, 1D" begin
    x = JArray(foo, (3,))
    @test typeof(x) == JArray{Float64, 1, Vector{Float64}}
end

@testset "JArray construct, 2D" begin
    x = JArray(bar, (3,2))
    @test typeof(x) == JArray{Float64, 2, Array{Float64,2}}
end

@testset "JArray size, 1D" begin
    x = JArray(foo, (3,))
    @test size(x) == (6,)
    @test length(x) == 6
end

@testset "JArray size, 2D" begin
    x = JArray(bar, (3,4))
    @test size(x) == (6,12)
end

@testset "JArray getblock, 1D" begin
    x = JArray(foo, (3,))
    getblock(x,1)
    for i = 1:3
        @test getblock(x, i) ≈ i*ones(2)
    end
end

@testset "JArray getblock, 2D" begin
    x = JArray(bar, (3,4))
    for jblock = 1:4, iblock = 1:3
        xᵢⱼ = getblock(x, iblock, jblock)
        @test xᵢⱼ ≈ LinearIndices((3,4))[iblock,jblock]*ones(2,3)
    end
end

@testset "JArray setblock!, 1D" begin
    x = JArray(foo, (3,))
    for iblock = 1:3
        setblock!(x, 2*iblock*ones(2), iblock)
        xᵢ = getblock(x, iblock)
        @test xᵢ ≈ 2*iblock*ones(2)
    end
end

@testset "JArray setblock!, 2D" begin
    x = JArray(bar, (3,4))
    for jblock = 1:4, iblock = 1:3
        setblock!(x, 2*LinearIndices((3,4))[iblock,jblock]*ones(2,3), iblock, jblock)
        xᵢⱼ = getblock(x, iblock, jblock)
        @test xᵢⱼ ≈ 2*LinearIndices((3,4))[iblock,jblock]*ones(2,3)
    end
end

@testset "JArray getindex, 1D" begin
    x = JArray(foo, (3,))
    for i = 1:length(x)
        if i <= 2
            @test x[i] ≈ 1
        elseif 2 < i <= 4
            @test x[i] ≈ 2
        elseif 4 < i <= 6
            @test x[i] ≈ 3
        end
    end
end

@testset "JArray getindex, 2D" begin
    x = JArray(bar, (3,2))
    for j = 1:size(x,2), i = 1:size(x,1)
        if 1 <= i <= 2 && 1 <= j <= 3
            @test x[i,j] ≈ 1
        elseif 3 <= i <= 4 && 1 <= j <= 3
            @test x[i,j] ≈ 2
        elseif 5 <= i <= 6 && 1 <= j <= 3
            @test x[i,j] ≈ 3
        elseif 1 <= i <= 2 && 4 <= j <= 6
            @test x[i,j] ≈ 4
        elseif 3 <= i <= 4 && 4 <= j <= 6
            @test x[i,j] ≈ 5
        elseif 5 <= i <= 6 && 4 <= j <= 6
            @test x[i,j] ≈ 6
        end
    end
end

@testset "JArray similar, 1D" begin
    x = JArray(i->rand(2), (3,))
    y = similar(x)
    z = similar(x, Float32)
    @test length(x) == length(y)
    @test length(z) == length(y)
    @test eltype(z) == Float32
    @test z.indices == x.indices
    @test y.indices == x.indices
end

@testset "JArray similar, 2D" begin
    x = JArray(i->rand(2,3), (3,2))
    y = similar(x)
    z = similar(x, Float32)
    @test length(x) == length(y)
    @test length(z) == length(y)
    @test size(x) == size(y)
    @test size(z) == size(y)
    @test eltype(z) == Float32
    @test z.indices == x.indices
    @test y.indices == x.indices
end

@testset "JArray procs" begin
    x = JArray(i->rand(2), (3,))
    @test procs(x) == workers()
end

@testset "JArray collect/convert, 1D" begin
    x = JArray(i->rand(2), (3,))
    workers()
    x.blockmap
    y = collect(x)
    workers()

    z = convert(Jets.BlockArray, x)
    u = convert(Array, x)

    @test isa(y, Jets.BlockArray)
    @test isa(z, Jets.BlockArray)
    @test isa(u, Array)
    for i = 1:6
        @test x[i] ≈ y[i]
        @test x[i] ≈ z[i]
        @test x[i] ≈ u[i]
    end
end

@testset "JArray norm, 1D" begin
    x = rand(90)
    y = JArray(i->x[(i[1]-1)*30+1:i[1]*30], (3,))
    z = collect(y)
    @test x ≈ z

    @test norm(x) ≈ norm(y)
    for p in (0,1,2,Inf,-Inf)
        @test norm(x,p) ≈ norm(y,p)
    end
end

@testset "JArray dot, 1D" begin
    x₁ = rand(90)
    x₂ = rand(90)
    y₁ = JArray(i->x₁[(i[1]-1)*30+1:i[1]*30], (3,))
    y₂ = JArray(i->x₂[(i[1]-1)*30+1:i[1]*30], (3,))
    z₁ = collect(y₁)
    z₂ = collect(y₂)
    @test x₁ ≈ z₁
    @test x₂ ≈ z₂
    @test dot(x₁,x₂) ≈ dot(y₁,y₂)
end

@testset "JArray extrema, 1D" begin
    x = rand(90)
    y = JArray(i->x[(i[1]-1)*30+1:i[1]*30], (3,))
    z = collect(y)
    mnmx = extrema(y)
    _mnmx = extrema(x)
    @test mnmx[1] ≈ _mnmx[1]
    @test mnmx[2] ≈ _mnmx[2]
end

@testset "JArray fill!, 1D" begin
    y = JArray(i->zeros(30), (3,))
    y .= 3
    z = collect(y)
    for i = 1:90
        @test z[i] ≈ 3
    end
end

@testset "JArray broadcasting" begin
    u = JArray(i->rand(10), (4,))
    v = JArray(i->rand(10), (4,))
    w = JArray(i->rand(10), (4,))
    a = rand(Float64)
    b = rand(Float64)
    c = rand(Float64)
    x = a*u .+ b*v .+ c*w

    _u = collect(u)
    _v = collect(v)
    _w = collect(w)
    @test a*_u .+ b*_v .+ c*_w ≈ collect(x)

    y = JArray(i->zeros(10), (4,))
    y .= x
    @test typeof(y) == JournaledJets.JArray{Float64,1,Array{Float64,1}}
    @test a*_u .+ b*_v .+ c*_w ≈ collect(y)

    z = zeros(40)
    z .= x
    @test a*_u .+ b*_v .+ c*_w ≈ z

    x = rand(Int32, 40)
    y = JArray(i->rand(10), (4,))
    z = x .* y
    @test isa(z, JournaledJets.JArray)

    @test convert(Array, z) ≈ x .* convert(Array, y)
    x = similar(y, Int)
    @test isa(x, JArray{Int})
end

@testset "JetJSpace construct" begin
    R = JetJSpace(JArray(_->[JetSpace(Float64,10)], (5,)))
    @test isa(R, JetJSpace)
end

@testset "JetJSpace size" begin
    R = JetJSpace(JArray(_->[JetSpace(Float64,10)], (5,)))
    @test size(R) == (50,)
end

@testset "JetJSpace eltype" begin
    R = JetJSpace(JArray(_->[JetSpace(Float64,10)], (5,)))
    @test eltype(R) == Float64
end

@testset "JetJSpace indices" begin
    R = JetJSpace(JArray(_->[JetSpace(Float64,10)], (5,)))
    @test indices(R) == [1:10, 11:20, 21:30, 31:40, 41:50]
end

@testset "JetJSpace space" begin
    R = JetJSpace(JArray(_->[JetSpace(Float64,10)], (5,)))
    @test space(R, 2) == JetSpace(Float64,10)
end

@testset "JetJSpace nblocks" begin
    R = JetJSpace(JArray(_->[JetSpace(Float64,10)], (5,)))
    @test nblocks(R) == 5
end

@testset "JetJSpace Array" begin
    R = JetJSpace(JArray(_->[JetSpace(Float64,10)], (5,)))
    x = Array(R)
    @test isa(x, JArray)
    @test size(x) == (50,)
end

@testset "JetJSpace ones" begin
    R = JetJSpace(JArray(_->[JetSpace(Float64,10)], (5,)))
    x = ones(R)
    @test isa(x, JArray)
    @test size(x) == (50,)
    @test collect(x) ≈ ones(50)
end

@testset "JetJSpace zeros" begin
    R = JetJSpace(JArray(_->[JetSpace(Float64,10)], (5,)))
    x = zeros(R)
    @test isa(x, JArray)
    @test size(x) == (50,)
    @test collect(x) ≈ zeros(50)
end

@testset "JetJSpace rand" begin
    R = JetJSpace(JArray(_->[JetSpace(Float64,10)], (5,)))
    x = rand(R)
    @test isa(x, JArray)
    @test size(x) == (50,)
end

@testset "block operator, tall and skinny" begin
    A = @blockop JArray(i->[JopBar(10)], (4,1))
    @test isa(A, Jop)

    _F  = JArray(_->[JopBar(10)], (4,1))
    F = @blockop _F

    _G = [_F[i,1] for i in 1:4, j in 1:1]
    G = @blockop _G

    m = rand(domain(F))
    @test collect(F*m) ≈ G*m

    J = jacobian!(F, m)
    _J = jacobian!(G, m)

    δm = rand(domain(J))
    @test collect(J*δm) ≈ _J*δm
    δd = rand(range(J))
    @test J'*δd ≈ _J'*collect(δd)
end

@testset "scale-up" begin
    x = JArray(_->rand(100), (10,1))
    y = JArray(_->rand(100), (10,1))
    a = dot(x,y)
    _x = convert(Array, x)

    addprocs(3)
    @everywhere using JournaledJets

    @everywhere _getpids_fromid(id) = JournaledJets.registry[id].pids
    @everywhere function getpids(x)
        pids = Dict{Any,Any}()
        for pid in x.pids
            pids[pid] = remotecall_fetch(_getpids_fromid, pid, x.id)
        end
        pids
    end

    @everywhere _getblkmap_fromid(id) = JournaledJets.registry[id].blockmap
    @everywhere function getblkmaps(x)
        bm = [1:10;]
        for pid in x.pids
            bm = [bm remotecall_fetch(_getblkmap_fromid, pid, x.id)]
        end
        bm
    end

    _a = dot(x,y)
    @test a ≈ _a

    p = getpids(x)
    bm = getblkmaps(x)

    for key in keys(p)
        @test p[key] ≈ p[first(keys(p))]
    end

    for i = 2:size(bm,2)
        @test bm[:,i] ≈ bm[:,2]
    end

    @test _x ≈ convert(Array, x)

    close(x)
    close(y)
    rmprocs(workers()[4:6])
end

#@testset "journal, broadcasting" begin # TODO - does not work inside a test-set
######
######
___x = @journal JArray(i->i.I[1]*ones(2), (2,))
___y = @journal JArray(i->2*i.I[1]*ones(2), (2,))
@journal ___y ___y .*= 2
@journal ___y ___y .*= 2*___x .+ ___y
@journal ___y ___y .*= 2*___x .+ 3*___y

____y = copy(___y)
___y .= 0
@replay ___y 1
@test getblock(___y, 1) ≈ getblock(____y, 1)
@test (getblock(___y, 2) ≈ getblock(____y, 1)) == false
@replay ___y 2
@test getblock(___y, 2) ≈  getblock(____y, 2)

___y .= @journal :reset ___y 2*___x
@test length(___y.journal) == 2

____y = copy(___y)
___y .= 0
@replay ___y 1
@test getblock(___y, 1) ≈ getblock(____y, 1)
@test (getblock(___y, 2) ≈ getblock(____y, 1)) == false
@replay ___y 2
@test getblock(___y, 2) ≈  getblock(____y, 2)
######
######
#end

#@testset "journal, single block" begin  # TODO - does not work inside a test-set
######
######
x = @journal JArray(i->ones(2), (2,))

@journal x begin
    _x = getblock(x,1)
    _x[1] = π
end

x .= 0
@replay x 1
@test getblock(x, 1) ≈ [π, 1]

@replay x 2
@test getblock(x, 2) ≈ [1, 1]
######
######
#end

#@testset "journal, mul!, single operator"
######
######
A = @blockop @journal JArray(_->[JopBar(2)], (2,1))
x = rand(domain(A))

y = @journal zeros(range(A))
@journal y mul!(y, A, x)

_y = copy(y)
y .= 0

@test (y ≈ _y) == false
@replay y 1
@replay y 2
@test y ≈ _y

A₁₁ = getblock(A, 1, 1)
state(A₁₁).x[1] = 2.0
mul!(y, A, x)
@test (y ≈ _y) == false

@replay A 1

mul!(y, A, x)
@test y ≈ _y
######
######
#end

#@testset "journal, mul!, composite operator" begin
######
######
A = @blockop @journal JArray(_->[JopBar(2) ∘ JopBar(2)], (2,1))
x = rand(domain(A))

y = @journal zeros(range(A))
y = @journal y mul!(y, A, x)

_y = copy(y)
y .= 0

@replay y 1
@replay y 2
@test y ≈ _y

A₁₁ = getblock(A, 1, 1)
A₁₁₁ = state(A₁₁).ops[1]
state(A₁₁₁).x[1] = 2.0
mul!(y, A, x)
@test (y ≈ _y) == false

@replay A 1

mul!(y, A, x)
@test y ≈ _y
######
######
#end

rmprocs(workers())