text:
# function convolution_kernel(
#         state,
#         imgSrc::AbstractArray{T},
#         kernelValues,
#         kernel_width,
#         imgConvolved,
#         ::Val{BLOCK_SIZE},
#         ::Val{LOCAL_WIDTH}
#     ) where {T, BLOCK_SIZE, LOCAL_WIDTH}
#     1 = 1; 0 = 0
#     w = kernel_width
#     wBy2 = w >> 1 #w divided by 2
#     #Goes up to 15x15 filters
#     ptr = LocalMemory(state, T, LOCAL_WIDTH) # local width need to be static, so calculating it from block size won't cut it
#     P = CLArrays.LocalArray{T, 2}(ptr, (LOCAL_WIDTH, LOCAL_WIDTH))
#
#     i = blockidx_x(state)
#     j = blockidx_y(state) #Identification of work-item
#     idX = threadidx_x(state)
#     idY = threadidx_y(state)
#
#     ii = i*BLOCK_SIZE + idX; # == get_global_id(0);
#     jj = j*BLOCK_SIZE + idY; # == get_global_id(1);
#     #Reads pixels
#     P[idX, idY] = imgSrc[ii, jj]
#     #Needs to read extra elements for the filter in the borders
#     if (idX < w)
#         P[idX + BLOCK_SIZE, idY] = imgSrc[ii + BLOCK_SIZE, jj]
#     end
#     if (idY < w)
#         P[idX, idY + BLOCK_SIZE] = imgSrc[ii, jj + BLOCK_SIZE]
#     end
#     synchronize_threads(state)
#     ##############
#     convPix = zero(T);
#     for ix = 0:(w - 1)
#         for jy = 0:(w - 1)
#             temp = P[ix, jy]
#             convPix += temp * kernelValues[ix + w*jy]
#         end
#     end
#     ##############
#     synchronize_threads(state)
#     imgConvolved[ii + wBy2, jj + wBy2] = P[idX + wBy2, idY + wBy2]
#     return
# end

function convolution_kernel(state, A::AbstractArray{T}, out, K, Asize, Ksize) where T
    ilin = linear_index(state)
    idx = gpu_ind2sub(Asize, ilin)
    if idx[1] >= Asize[1] - Ksize[1] || idx[2] >= Asize[2] - Ksize[2]
        return
    end
    accum = zero(T)
    kw, kh = Ksize[1], Ksize[2]
    for ix = 0:(kw - 1)
        for jy = 0:(kh - 1)
            temp = A[gpu_sub2ind(Asize, idx .+ (ix, jy))]
            accum += temp * K[ix + kw*jy + 1]
        end
    end
    out[ilin] = accum
    return
end

function convolution!(a, out, k)
    gpu_call(convolution_kernel, a, (a, out, k, size(a), size(k)))
    synchronize(out)
    out
end

struct FFTKernel{T}
    kernel::T
    irfftplan
    rfftplan
end

function fftkernel(A, kernel)
    plan_rfft!(A)

end

function convolution_fft!(a, out, k)
    irfft(rfft(A) .* conj(rfft(krn)), length(axes(A, 1)))
    out
end


text:
#Find Sqrt of a number using simple math

function sqrtnum(num)
	return num^0.5
end

print("Enter the Number : ")
num = parse(Int64,readline())
print(sqrtnum(num))

ext: jl

max_stars_repo_name: cashewcodes/Basic-Julia-Projects

max_issues_repo_name: cashewcodes/Basic-Julia-Projects

max_forks_repo_name: cashewcodes/Basic-Julia-Projects


text:
# Sze, S. M., & Irvin, J. C. (1968). Resistivity, mobility and impurity levels in GaAs, Ge, and Si at 300°K. Solid-State Electronics, 11(6), 599–602.
# https://doi.org/10.1016/0038-1101(68)90012-9

const pSi_data = [141241312439906.34  99.9999999999999
194499370591877.62  68.46461639456666
295571204689629.06  43.45358317903116
460381216646110.75  27.579412410858527
634034977180747.5  19.860447554400164
851955412933557.2  14.66769850372177
1089774332370062.4  11.393909366437962
1615876286771859.2  7.606239889282181
2225377681726411.2  5.477394737670293
4018372494062520  3.142373253889367
8007263534267784  1.6712242736564364
12473787575660818  1.1441972880517628
21445649112993600  0.7638325682336568
32598724792618628  0.5641195037265176
56048119498669190  0.3862222542332937
85204322253373860  0.3000185311225828
142935426217454880  0.21333771623188874
217295259830406370  0.16782740849919364
284918519217386880  0.1371233823953912
346988588390428600  0.1208556366261384
412285586887165300  0.10651782831751473
514617602216427500  0.09153955041992637
674829761472440400  0.07866748152341368
929519246185815700  0.06188572649236628
1312362103166414300  0.04992920623643243
1638099075619741200  0.04290828271566711
2429349658378808000  0.03168940440349858
3346142217080559600  0.024616418285271273
4723791645850668000  0.018645186320546278
6668636915521611000  0.014122402735408967
8743552375008310000  0.011250926818107056
13959649881815632000  0.0073234654360613645
18759330743571850000  0.005688888505961023
23996982286868873000  0.00453217978952392
33048537964912603000  0.0032637069157250395
42277677767648020000  0.00266661170206652
56811251490196600000  0.0020197687353780474
76340954998056850000  0.0015298311866137819
116016856048370400000  0.0009957993682923774
144799950023552320000  0.00081361786366439
189845222834170920000  0.000632020640126649
225550331083675270000  0.0005296010185684199
342788745840653200000  0.0003535461151468624
407258681371465600000  0.0002962535885776857
483853789148391500000  0.00024824537729314016
589156034787413000000  0.00019777021128675618
699993018629596900000  0.000169960245280993
852372589147036700000  0.0001388660166661159]


text:
using Weave

weave("guiding_center_dynamics_splitting.jmd",
  out_path="build",
  doctype = "github")

weave("guiding_center_dynamics_vprk.jmd",
  out_path="build",
  doctype = "github")


text:
"""
Normalize brightness for each time point
Version Comment
0.1     hf, initial version
0.2     hf, don't use fitting
0.2.1   hf, medianfilter with repeat padding
"""

using Images
using Statistics
using LsqFit
using Clustering

"""
Extract feature from image
"""
function extract_feature(imgs)
    z_depth, t_len = size(imgs)[3:4]
    low_line = zeros(t_len)
    mean_line  = zeros(t_len)
    @inbounds Threads.@threads for t in 1:t_len
		nozeros_pixel = imgs[:,:,:,t][imgs[:,:,:,t].≠0]
		if length(nozeros_pixel) == 0
			#println("no cell in time point $t")
			continue # skip empty stack
		end
        low_line[t] = minimum(nozeros_pixel)
        mean_line[t] = mean(nozeros_pixel)
    end
    low_line, mean_line
end

medianfilter1(v, ws) = [median(v[i:(i+ws-1)]) for i=1:(length(v)-ws+1) ]
#medianfilter(v, ws) = medianfilter1(vcat(0,v,0),ws)
medianfilter(v, ws) = medianfilter1(vcat(v[1], v, v[end]),ws) # repeat

"""
Use k-means to remove outliers
Input raw data series, return index without outliers
"""
function kmeans_rm_outliers( raw_line )
	# TODO: if only one cluster?
	t_len = length(raw_line)
    X = zeros(2, t_len)
    X[1, :] = fill(0.02, t_len) # just choose 0.02 in random.
    X[2, :] = raw_line; # forgive me just add one more dimension for k-means
    result = kmeans(X, 2; init=:kmpp)
    if result.counts[1]> result.counts[2] # chose the one have more points
        index = 1
    else
        index = 2
    end
    t_nooutliers = [ i  for i in 1:t_len  if result.assignments[i] == index ]
end

"""
Removing outliers and smoothing feature
"""
function correct_feature( low_line, mean_line, fitted::Bool=false )
    t_len = length(low_line)
    ws = 3 # use small window to remove spark
    low_line = medianfilter(low_line, ws)
    mean_line = medianfilter(mean_line, ws)
	# TODO: if only one cluster? Add additional judge by std(mean_line[1:end-1].-mean_line[2:end])
	if fitted
		t_nooutliers = kmeans_rm_outliers(mean_line)
	    mean_line_nooutliers = mean_line[t_nooutliers]
		# assume mean_line and low_line share same shape
	    low_line_nooutliers = low_line[t_nooutliers]
    	@. model(x, p) = p[1]*exp(-x*p[2]) # fit with expontional function
	    fit_mean = curve_fit(model, t_nooutliers, mean_line_nooutliers, [0.003, 0.0003])
    	fit_low = curve_fit(model, t_nooutliers, low_line_nooutliers, [0.002, 0.0003])
	    corrected_mean_line = model(1:t_len, fit_mean.param)
    	corrected_low_line = model(1:t_len, fit_low.param)
		# TODO: what if mean.(fit, no_outliers) ?
	else
		corrected_mean_line = mean_line
		corrected_low_line = low_line
	end

    corrected_low_line, corrected_mean_line
end

"Normalize nucleus to uniform intensity"
function normalize(imgs)
    #z_depth = 20
    #t_len = size(imgs)[3] ÷ z_depth
    #d1_len, d2_len  = size(imgs)[1],  size(imgs)[2]
    d1_len, d2_len, z_depth, t_len = size(imgs)
    low_line, mean_line = extract_feature(imgs)
    corrected_low_line, corrected_mean_line = correct_feature(low_line, mean_line)
    imgs_norm = zeros(N0f16, size(imgs))
	imgs_len = length(imgs_norm[:, :, :, 1])
    @inbounds Threads.@threads for t in 1:t_len
		"""
		z = (t-1)*20+1 : 20*t
        imgs_norm[:,:,z] = ((imgs[:,:,z] .- corrected_low_line[t]) ./ (corrected_mean_line[t]-corrected_low_line[t]).*0.005 .+ 0.02) .* ( imgs[:, :, z] .>0)
        #print("$t ")
		"""
		tmp = view(imgs, :, :, :, t)
		tmp_norm = view(imgs_norm, :, :, :, t)
		for i in 1:imgs_len
			if tmp[i] > 0
				tmp_norm_value = (tmp[i] - corrected_low_line[t]) /
					(corrected_mean_line[t] - corrected_low_line[t])*0.005 + 0.02
				if tmp_norm_value < 0
					tmp_norm[i] = 0
				elseif tmp_norm_value > 1
					tmp_norm[i] = 1
				else
					tmp_norm[i]= tmp_norm_value
				end
			end
		end
    end

	imgs_norm, Dict("low_line"=>low_line, "mean_line"=>mean_line,
					"corrected_low_line"=>corrected_low_line, "corrected_mean_line"=>corrected_mean_line)
end

#@time s14c1 = load(File(format"TIFF", "$data_dir/$(cell_name[14])"));
#@time imgs_norm, low_line, mean_line = normalize(s14c1);
#plot(low_line)
#plot!(low_line_raw, marker=(:star))
#plot!(mean_line)
#plot!(mean_line_raw, marker=(:dot))
#plot!(ylim=(0,0.004))


text:
"""
    abstract type FETerm <: GridapType end

A `FETerm` is a lazy representation of a summand of a finite element problem. It is not assembled.
See also [FEOperator](@ref).
"""
abstract type FETerm <: GridapType end

"""
Returns an object representing
the contribution to the Jacobian of the given term. Returns nothing if the term
has not contribution to the Jacobian (typically for source terms)
"""
function get_cell_jacobian(t::FETerm,uh,du,v)
  @abstractmethod
end

"""
Returns an object representing the contribution to the
residual of the given term. Returns always something.
"""
function get_cell_residual(t::FETerm,uh,v)
  @abstractmethod
end

"""
"""
function get_cell_id(t::FETerm)
  @abstractmethod
end

"""
"""
function get_cell_jacobian_and_residual(t::FETerm,uh,du,v)
  @assert is_a_fe_function(uh)
  @assert is_a_fe_cell_basis(v)
  @assert is_a_fe_cell_basis(du)
  celljac = get_cell_jacobian(t,uh,du,v)
  cellres = get_cell_residual(t,uh,v)
  _setup_jac_and_res(celljac,cellres)
end

function _setup_jac_and_res(celljac,cellres)
  celljacres = pair_arrays(celljac,cellres)
  (celljacres, nothing, nothing)
end

function _setup_jac_and_res(celljac::Nothing,cellres)
  (nothing, nothing, cellres)
end

function _setup_jac_and_res(celljac,cellres::Nothing)
  (nothing, celljac, nothing)
end

function _setup_jac_and_res(celljac::Nothing,cellres::Nothing)
  (nothing, nothing, nothing)
end

"""
See also [FETerm](@ref).
"""
abstract type AffineFETerm <: FETerm end

"""
Returns an object representing the contribution to the
system matrix of the given term. Returns nothing if the term has not
contribution (typically for source terms)
"""
function get_cell_matrix(t::AffineFETerm,u,v)
  @abstractmethod
end

"""
Returns an object (e.g. a CellVector) representing the contribution to the
system rhs of the given term (with Dirichlet bcs included). Returns nothing if the term has not
contribution (typically for linear terms)
"""
function get_cell_vector(t::AffineFETerm,uhd,v)
  @abstractmethod
end

"""
Returns an object (e.g. a CellVector) representing the contribution to the
system rhs of the given term (without Dirichlet bcs). Returns nothing if the term has not
contribution (typically for linear terms)
"""
function get_cell_vector(t::AffineFETerm,v)
  @abstractmethod
end

function get_cell_jacobian(t::AffineFETerm,uh,du,v)
  @assert is_a_fe_function(uh)
  @assert is_a_fe_cell_basis(v)
  @assert is_a_fe_cell_basis(du)
  get_cell_matrix(t,du,v)
end

"""
"""
function get_cell_matrix_and_vector(t::AffineFETerm,uhd,u,v)
  @assert is_a_fe_function(uhd)
  @assert is_a_fe_cell_basis(v)
  @assert is_a_fe_cell_basis(u)
  cellmat = get_cell_matrix(t,u,v)
  cellvec = get_cell_vector(t,v)
  cellvals = get_cell_values(uhd,get_cell_id(t))
  _setup_cell_matrix_and_vector(cellmat,cellvec,cellvals)
end

function  _setup_cell_matrix_and_vector(cellmat,cellvec,cellvals)
  cellmatvec = pair_arrays(cellmat,cellvec)
  cellmatvec_with_diri = attach_dirichlet(cellmatvec,cellvals)
  (cellmatvec_with_diri, nothing, nothing)
end

function _setup_cell_matrix_and_vector(cellmat,cellvec::Nothing,cellvals)
  cellmatvec_with_diri = attach_dirichlet(cellmat,cellvals)
  (cellmatvec_with_diri, nothing, nothing)
end

function _setup_cell_matrix_and_vector(cellmat::Nothing,cellvec,cellvals)
  (nothing, nothing, cellvec)
end

function _setup_cell_matrix_and_vector(cellmat::Nothing,cellvec::Nothing,cellvals)
  (nothing, nothing, nothing)
end

"""
See also [FETerm](@ref).
"""
abstract type LinearFETerm <: AffineFETerm end

"""
"""
abstract type FESource <: AffineFETerm end

function get_cell_matrix(t::FESource,u,v)
  nothing
end

# Working with a collection of FETerms

"""
"""
function collect_cell_jacobian(uh,du,v,terms)
  @assert is_a_fe_function(uh)
  @assert is_a_fe_cell_basis(v)
  @assert is_a_fe_cell_basis(du)
  w = []
  r = []
  c = []
  for term in terms
    cellvals = get_cell_jacobian(term,uh,du,v)
    cellids = get_cell_id(term)
    _push_matrix_contribution!(w,r,c,cellvals,cellids)
  end
  (w,r,c)
end

"""
"""
function collect_cell_matrix(u,v,terms)
  @assert is_a_fe_cell_basis(v)
  @assert is_a_fe_cell_basis(u)
  w = []
  r = []
  c = []
  for term in terms
    cellvals = get_cell_matrix(term,u,v)
    cellids = get_cell_id(term)
    _push_matrix_contribution!(w,r,c,cellvals,cellids)
  end
  (w,r,c)
end

"""
"""
function collect_cell_residual(uh,v,terms)
  @assert is_a_fe_function(uh)
  @assert is_a_fe_cell_basis(v)
  w = []
  r = []
  for term in terms
    cellvals = get_cell_residual(term,uh,v)
    cellids = get_cell_id(term)
    _push_vector_contribution!(w,r,cellvals,cellids)
  end
  (w,r)
end

"""
"""
function collect_cell_vector(uhd,v,terms)
  @assert is_a_fe_cell_basis(v)
  @assert is_a_fe_function(uhd)
  w = []
  r = []
  for term in terms
    cellvals = get_cell_vector(term,uhd,v)
    cellids = get_cell_id(term)
    _push_vector_contribution!(w,r,cellvals,cellids)
  end
  (w,r)
end

"""
"""
function collect_cell_matrix_and_vector(uhd,u,v,terms)
  @assert is_a_fe_cell_basis(v)
  @assert is_a_fe_cell_basis(u)
  @assert is_a_fe_function(uhd)

  matvecdata = ([],[],[])
  matdata = ([],[],[])
  vecdata = ([],[])

  for term in terms
    cellmatvec, cellmat, cellvec = get_cell_matrix_and_vector(term,uhd,u,v)
    cellids = get_cell_id(term)
    _push_matrix_contribution!(matvecdata...,cellmatvec,cellids)
    _push_matrix_contribution!(matdata...,cellmat,cellids)
    _push_vector_contribution!(vecdata...,cellvec,cellids)
  end

  (matvecdata, matdata, vecdata)
end

"""
"""
function collect_cell_jacobian_and_residual(uh,du,v,terms)
  @assert is_a_fe_function(uh)
  @assert is_a_fe_cell_basis(v)
  @assert is_a_fe_cell_basis(du)

  matvecdata = ([],[],[])
  matdata = ([],[],[])
  vecdata = ([],[])

  for term in terms
    cellmatvec, cellmat, cellvec = get_cell_jacobian_and_residual(term,uh,du,v)
    cellids = get_cell_id(term)
    _push_matrix_contribution!(matvecdata...,cellmatvec,cellids)
    _push_matrix_contribution!(matdata...,cellmat,cellids)
    _push_vector_contribution!(vecdata...,cellvec,cellids)
  end

  (matvecdata, matdata, vecdata)
end

function _push_matrix_contribution!(w,r,c,cellvals,cellids)
  push!(w,cellvals)
  push!(r,cellids)
  push!(c,cellids)
  nothing
end

function _push_matrix_contribution!(w,r,c,cellvals::Nothing,cellids)
  nothing
end

function _push_vector_contribution!(v,r,cellvals,cellids)
  push!(v,cellvals)
  push!(r,cellids)
  nothing
end

function _push_vector_contribution!(v,r,cellvals::Nothing,cellids)
  nothing
end

# Concrete implementations

struct AffineFETermFromIntegration <: AffineFETerm
  biform::Function
  liform::Function
  trian::Triangulation
  quad::CellQuadrature
end

"""
"""
function AffineFETerm(
  biform::Function,liform::Function,trian::Triangulation,quad::CellQuadrature)
  AffineFETermFromIntegration(biform,liform,trian,quad)
end

function get_cell_matrix(t::AffineFETermFromIntegration,u,v)
  @assert is_a_fe_cell_basis(v)
  @assert is_a_fe_cell_basis(u)
  _v = restrict(v,t.trian)
  _u = restrict(u,t.trian)
  integrate(t.biform(_u,_v),t.trian,t.quad)
end

function get_cell_vector(t::AffineFETermFromIntegration,uhd,v)
  @assert is_a_fe_function(uhd)
  @assert is_a_fe_cell_basis(v)
  _v = restrict(v,t.trian)
  _uhd = restrict(uhd,t.trian)
  integrate(t.liform(_v)-t.biform(_uhd,_v),t.trian,t.quad)
end

function get_cell_vector(t::AffineFETermFromIntegration,v)
  @assert is_a_fe_cell_basis(v)
  _v = restrict(v,t.trian)
  integrate(t.liform(_v),t.trian,t.quad)
end

function get_cell_residual(t::AffineFETermFromIntegration,uh,v)
  @assert is_a_fe_function(uh)
  @assert is_a_fe_cell_basis(v)
  _v = restrict(v,t.trian)
  _uh = restrict(uh,t.trian)
  integrate(t.biform(_uh,_v)-t.liform(_v),t.trian,t.quad)
end

function get_cell_id(t::AffineFETermFromIntegration)
  get_cell_id(t.trian)
end

struct FESourceFromIntegration <: FESource
  liform::Function
  trian::Triangulation
  quad::CellQuadrature
end

"""
"""
function FESource(
  liform::Function,
  trian::Triangulation,
  quad::CellQuadrature)
  FESourceFromIntegration(liform,trian,quad)
end

function get_cell_matrix(t::FESourceFromIntegration,u,v)
  nothing
end

function get_cell_vector(t::FESourceFromIntegration,uhd,v)
  @assert is_a_fe_function(uhd)
  @assert is_a_fe_cell_basis(v)
  _v = restrict(v,t.trian)
  integrate(t.liform(_v),t.trian,t.quad)
end

function get_cell_vector(t::FESourceFromIntegration,v)
  @assert is_a_fe_cell_basis(v)
  _v = restrict(v,t.trian)
  integrate(t.liform(_v),t.trian,t.quad)
end

function get_cell_residual(t::FESourceFromIntegration,uh,v)
  @assert is_a_fe_function(uh)
  @assert is_a_fe_cell_basis(v)
  _v = restrict(v,t.trian)
  integrate(-t.liform(_v),t.trian,t.quad)
end

function get_cell_id(t::FESourceFromIntegration)
  get_cell_id(t.trian)
end

struct LinearFETermFromIntegration <: LinearFETerm
  biform::Function
  trian::Triangulation
  quad::CellQuadrature
end

"""
"""
function LinearFETerm(
  biform::Function,trian::Triangulation,quad::CellQuadrature)
  LinearFETermFromIntegration(biform,trian,quad)
end

function get_cell_matrix(t::LinearFETermFromIntegration,u,v)
  @assert is_a_fe_cell_basis(v)
  @assert is_a_fe_cell_basis(u)
  _v = restrict(v,t.trian)
  _u = restrict(u,t.trian)
  integrate(t.biform(_u,_v),t.trian,t.quad)
end

function get_cell_vector(t::LinearFETermFromIntegration,uhd,v)
  @assert is_a_fe_function(uhd)
  @assert is_a_fe_cell_basis(v)
  _v = restrict(v,t.trian)
  _uhd = restrict(uhd,t.trian)
  integrate(-t.biform(_uhd,_v),t.trian,t.quad)
end

function get_cell_vector(t::LinearFETermFromIntegration,v)
  @assert is_a_fe_cell_basis(v)
  nothing
end

function get_cell_residual(t::LinearFETermFromIntegration,uh,v)
  @assert is_a_fe_function(uh)
  @assert is_a_fe_cell_basis(v)
  _v = restrict(v,t.trian)
  _uh = restrict(uh,t.trian)
  integrate(t.biform(_uh,_v),t.trian,t.quad)
end

function get_cell_id(t::LinearFETermFromIntegration)
  get_cell_id(t.trian)
end

struct NonlinearFETerm <: FETerm
  res::Function
  jac::Function
  trian::Triangulation
  quad::CellQuadrature
end

"""
"""
function FETerm(
  res::Function, jac::Function, trian::Triangulation, quad::CellQuadrature)
  NonlinearFETerm(res,jac,trian,quad)
end

function get_cell_jacobian(t::NonlinearFETerm,uh,du,v)
  @assert is_a_fe_function(uh)
  @assert is_a_fe_cell_basis(v)
  @assert is_a_fe_cell_basis(du)
  _v = restrict(v,t.trian)
  _uh = restrict(uh,t.trian)
  _du = restrict(du,t.trian)
  integrate(t.jac(_uh,_du,_v),t.trian,t.quad)
end

function get_cell_residual(t::NonlinearFETerm,uh,v)
  @assert is_a_fe_function(uh)
  @assert is_a_fe_cell_basis(v)
  _v = restrict(v,t.trian)
  _uh = restrict(uh,t.trian)
  integrate(t.res(_uh,_v),t.trian,t.quad)
end

function get_cell_id(t::NonlinearFETerm)
  get_cell_id(t.trian)
end

"""
"""
struct AffineFETermFromCellMatVec <: AffineFETerm
  matvecfun::Function
  trian::Triangulation
end

function get_cell_matrix(t::AffineFETermFromCellMatVec,u,v)
  @assert is_a_fe_cell_basis(v)
  @assert is_a_fe_cell_basis(u)
  _v = restrict(v,t.trian)
  _u = restrict(u,t.trian)
  cellmatvec = t.matvecfun(_u,_v)
  cellmat, _ = unpair_arrays(cellmatvec)
  cellmat
end

function get_cell_vector(t::AffineFETermFromCellMatVec,uhd,v)
  @assert is_a_fe_function(uhd)
  @assert is_a_fe_cell_basis(v)
  trial = TrialFESpace(get_fe_space(uhd))
  u = get_cell_basis(trial)
  _v = restrict(v,t.trian)
  _u = restrict(u,t.trian)
  _cellvals = get_cell_values(uhd,get_cell_id(t.trian))
  cellmatvec = t.matvecfun(_u,_v)
  cellmatvec_with_diri = attach_dirichlet(cellmatvec,_cellvals)
  _, cellvec_with_diri = unpair_arrays(cellmatvec_with_diri)
  cellvec_with_diri
end

function get_cell_vector(t::AffineFETermFromCellMatVec,v)
  @assert is_a_fe_cell_basis(v)
  @unreachable "This function cannot be implemented for $(typeof(t)) objects."
end

function get_cell_id(t::AffineFETermFromCellMatVec)
  get_cell_id(t.trian)
end

function get_cell_matrix_and_vector(t::AffineFETermFromCellMatVec,uhd,u,v)
  @assert is_a_fe_function(uhd)
  @assert is_a_fe_cell_basis(v)
  @assert is_a_fe_cell_basis(u)
  _v = restrict(v,t.trian)
  _u = restrict(u,t.trian)
  _cellvals = get_cell_values(uhd,get_cell_id(t.trian))
  cellmatvec = t.matvecfun(_u,_v)
  cellmatvec_with_diri = attach_dirichlet(cellmatvec,_cellvals)
  (cellmatvec_with_diri, nothing, nothing)
end

"""
"""
struct FETermFromCellJacRes <: FETerm
  jacresfun::Function
  trian::Triangulation
end

function get_cell_jacobian(t::FETermFromCellJacRes,uh,du,v)
  @assert is_a_fe_function(uh)
  @assert is_a_fe_cell_basis(v)
  @assert is_a_fe_cell_basis(du)
  _v = restrict(v,t.trian)
  _du = restrict(du,t.trian)
  _uh = restrict(uh,t.trian)
  celljacres = t.jacresfun(_uh,_du,_v)
  celljac, _ = unpair_arrays(celljacres)
  celljac
end

function get_cell_residual(t::FETermFromCellJacRes,uh,v)
  @assert is_a_fe_function(uh)
  @assert is_a_fe_cell_basis(v)
  trial = TrialFESpace(get_fe_space(uh))
  du = get_cell_basis(trial)
  _v = restrict(v,t.trian)
  _du = restrict(du,t.trian)
  _uh = restrict(uh,t.trian)
  celljacres = t.jacresfun(_uh,_du,_v)
  _, cellres = unpair_arrays(celljacres)
  cellres
end

function get_cell_id(t::FETermFromCellJacRes)
  get_cell_id(t.trian)
end

function get_cell_jacobian_and_residual(t::FETermFromCellJacRes,uh,du,v)
  @assert is_a_fe_function(uh)
  @assert is_a_fe_cell_basis(v)
  @assert is_a_fe_cell_basis(du)
  _v = restrict(v,t.trian)
  _du = restrict(du,t.trian)
  _uh = restrict(uh,t.trian)
  celljacres = t.jacresfun(_uh,_du,_v)
  (celljacres, nothing, nothing)
end


text:
#! /usr/bin/julia

# Rosetta Code, Show the epoch

println("Time zero (the epoch) is ", strftime("%c", 0), ".")

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
# examples/cancellation.jl shows effects of these settings
const POWERSERIES_ZMAX = 1.5
const POWERSERIES_NMAX = 11
const POWERSERIES_ZMAX_PRIME = 2.0
const POWERSERIES_NMAX_PRIME = 13

# GSL implementation faster than SpecialFunctions
import GSL
# This survives large arguments better
besselK0(z) = GSL.sf_bessel_K0_scaled(z)*exp(-z)
besselK1(z) = GSL.sf_bessel_K1_scaled(z)*exp(-z)
besselI0(z) = GSL.sf_bessel_I0_scaled(z)*exp(z)
besselI1(z) = GSL.sf_bessel_I1_scaled(z)*exp(z)

# besselK0(z) = besselk(0, z)
# besselK1(z) = besselk(1, z)
# besselI0(z) = besseli(0, z)
# besselI1(z) = besseli(1, z)

# ========== STOKESLET
"Computes S_{ij} f_j"
function stokeslet(r, f, alpha)
    rsq = r[1]*r[1] + r[2]*r[2]
    zsq = rsq*alpha*alpha
    z = sqrt(zsq)
    K0 = besselK0(z)
    K1 = besselK1(z)
    S1 = (zsq*K0 + z*K1 - 1) / (2*pi*zsq)
    S2 = -(zsq*K0 + 2*z*K1 - 2) / (2*pi*zsq*zsq)
    rdotf = (r[1]*f[1] + r[2]*f[2])
    return f.*S1 .+ r.*(alpha*alpha*S2*rdotf)
end

function stokeslet_kernel_noarrays(r1, r2, alpha)
    rsq = r1*r1 + r2*r2
    zsq = rsq*alpha*alpha
    z = sqrt(zsq)
    K0 = besselK0(z)
    K1 = besselK1(z)
    S1 = (zsq*K0 + z*K1 - 1) / (2*pi*zsq)
    S2 = -(zsq*K0 + 2*z*K1 - 2) / (2*pi*zsq*zsq)
    A11 = S1 + alpha^2*S2*r1*r1
    A12 = alpha^2*S2*r1*r2
    A21 = A12
    A22 = S1 + alpha^2*S2*r2*r2
    return A11, A12, A21, A22
end

# ========== PRESSURE(LET)
"Computes p_i f_i"
function pressure(r, f, alpha)
    rsq = r[1]*r[1] + r[2]*r[2]
    rdotf = (r[1]*f[1] + r[2]*f[2])
    return rdotf/(2*pi*rsq)
end

# ========= STRESSLET
"Computes T_{ijk} f_i n_k"
function stresslet(r::Array{Float64}, f::Array{Float64},
                   n::Array{Float64}, alpha::Float64)
    rsq = r[1]*r[1] + r[2]*r[2]
    rnorm = sqrt(rsq)
    z = rnorm*alpha
    (T1, T2, T3) = Ti(z)

    rdotf = r[1]*f[1] + r[2]*f[2]
    rdotn = r[1]*n[1] + r[2]*n[2]
    fdotn = f[1]*n[1] + f[2]*n[2]

    alpha2 = alpha*alpha
    alpha4 = alpha2*alpha2

    u1 = alpha2*T1*(r*fdotn + f*rdotn + n*rdotf)
    u2 = alpha4*T2*rdotf*rdotn*r
    u3 = alpha2*T3*r*fdotn;

    return u1 + u2 + u3
end

"Computes K_{ji} = T_{ijk} n_k"
function dblkernel(r::Array{Float64},
                   n::Array{Float64},
                   alpha::Float64 )
    K = Array{Float64}(undef, 2,2)
    dblkernel!(r, n, alpha, K)
    return K
end

function dblkernel!(r::Array{Float64},
                    n::Array{Float64},
                    alpha::Float64,
                    K::Array{Float64})
    K11, K12, K21, K22 = dblkernel_noarrays(r[1], r[2], n[1], n[2], alpha)
    K[1,1] = K11
    K[1,2] = K12
    K[2,1] = K21
    K[2,2] = K22
end

# Compute kernel without passing around arrays
function dblkernel_noarrays(r1, r2, n1, n2, alpha)
    rsq = r1*r1 + r2*r2
    rnorm = sqrt(rsq)
    z = rnorm*alpha
    (T1, T2, T3) = Ti(z)
    alpha2 = alpha*alpha
    alpha4 = alpha2*alpha2
    T1 *= alpha2
    T2 *= alpha4
    T3 *= alpha2
    rdotn = r1*n1 + r2*n2
    K11 = T1*( n1*r1 + n1*r1 + rdotn ) +
        T2*r1*r1*rdotn +
        T3*n1*r1
    K21 = T1*( n2*r1 + n1*r2  ) +
        T2*r1*r2*rdotn +
        T3*n1*r2
    K12 = T1*( n1*r2 + n2*r1  ) +
        T2*r2*r1*rdotn +
        T3*n2*r1
    K22 = T1*( n2*r2 + n2*r2 + rdotn ) +
        T2*r2*r2*rdotn +
        T3*n2*r2
    K11, K12, K21, K22
end

"""
Computes K_^L{ji}=T_^L{ijk} n_k
# Returns: KL
"""
function dblkernel_log(r::Array{Float64},
                       n::Array{Float64},
                       alpha::Float64 )
    KL = Array{Float64}(undef, 2,2)
    dblkernel_log!(r, n, alpha, KL)
    return KL
end

function dblkernel_log!(r::Array{Float64},
                        n::Array{Float64},
                        alpha::Float64,
                        KL::Array{Float64,2} )
    rsq = r[1]*r[1] + r[2]*r[2]
    rnorm = sqrt(rsq)
    z = rnorm*alpha
    (T1S, T1L, T2S, T2L, T3S, T3L) = Ti_split(z)
    alpha2 = alpha*alpha
    alpha4 = alpha2*alpha2
    T1L *= alpha2
    T2L *= alpha4
    T3L *= alpha2
    rdotn = r[1]*n[1] + r[2]*n[2]
    for i=1:2
        for j=1:2
            KL[j + 2*(i-1)] =
                T1L*( n[j]*r[i] + n[i]*r[j] + (i==j)*rdotn ) +
                T2L*r[i]*r[j]*rdotn +
                T3L*n[i]*r[j]
        end
    end
end

"""
Computes T_{ijk} f_i n_k, but with kernel-splitting
# Returns: (uS, uL, uC, uQ)
such that
u = uS + uL*log(r) + uC*r_k*n_k/r^2 + uQ*f_i*r_i*r_j*r_k*n_k/r^4
"""
function stresslet_split(r, f, n, alpha)
    rnorm = sqrt(r[1]*r[1] + r[2]*r[2])
    z = rnorm*alpha
    (T1S, T1L, T2S, T2L, T3S, T3L) = Ti_split(z)
    rdotf = r[1]*f[1] + r[2]*f[2]
    rdotn = r[1]*n[1] + r[2]*n[2]
    fdotn = f[1]*n[1] + f[2]*n[2]

    vec1 = r*fdotn + f*rdotn + n*rdotf
    vec2 = r*(rdotf*rdotn)
    vec3 = r*fdotn

    alpha2 = alpha*alpha
    alpha4 = alpha2*alpha2

    uS = alpha2*(T1S*vec1 + alpha2*T2S*vec2 + T3S*vec3)
    uL = alpha2*(T1L*vec1 + alpha2*T2L*vec2 + T3L*vec3)
    uS = uS + log(alpha)*uL # Split on log(r)
    uC = alpha2/(8*pi)*r*rdotf
    uQ = -1/pi

    return (uS, uL, uC, uQ)
end

function stresslet_split_noarrays(r1, r2, f1, f2, n1, n2, alpha)
    rnorm = sqrt(r1*r1 + r2*r2)
    z = rnorm*alpha
    (T1S, T1L, T2S, T2L, T3S, T3L) = Ti_split(z)
    rdotf = r1*f1 + r2*f2
    rdotn = r1*n1 + r2*n2
    fdotn = f1*n1 + f2*n2

    vecA1 = r1*fdotn + f1*rdotn + n1*rdotf
    vecA2 = r2*fdotn + f2*rdotn + n2*rdotf
    vecB1 = r1*(rdotf*rdotn)
    vecB2 = r2*(rdotf*rdotn)
    vecC1 = r1*fdotn
    vecC2 = r2*fdotn

    alpha2 = alpha*alpha
    alpha4 = alpha2*alpha2

    uS1 = alpha2*(T1S*vecA1 + alpha2*T2S*vecB1 + T3S*vecC1)
    uS2 = alpha2*(T1S*vecA2 + alpha2*T2S*vecB2 + T3S*vecC2)
    uL1 = alpha2*(T1L*vecA1 + alpha2*T2L*vecB1 + T3L*vecC1)
    uL2 = alpha2*(T1L*vecA2 + alpha2*T2L*vecB2 + T3L*vecC2)
    uS1 += log(alpha)*uL1 # Split on log(r)
    uS2 += log(alpha)*uL2 # Split on log(r)
    uC1 = alpha2/(8*pi)*r1*rdotf
    uC2 = alpha2/(8*pi)*r2*rdotf
    uQ = -1/pi

    return (uS1, uS2, uL1, uL2, uC1, uC2, uQ)
end

"""
Stresslet functions
# Returns: (T1, T2, T3)
"""
function Ti(z::Float64)
    if z < POWERSERIES_ZMAX
        (T1S, T1L, T2S, T2L, T3S, T3L) = Ti_split_pow(z)
        T1 = T1S + T1L*log(z)
        T2 = T2S + T2L*log(z) + 1/(8*pi*z^2) - 1/(pi*z^4)
        T3 = T3S + T3L*log(z)
        return (T1, T2, T3)
    else
        return Ti_direct(z)
    end
end

function Ti_direct(z::Float64)
    zsq = z*z
    z4 = zsq*zsq
    z6 = z4*zsq
    K0 = besselK0(z)
    K1 = besselK1(z)
    T1 = -(2*zsq*K0 + (zsq+4)*z*K1 - 4) / (2*pi*z4)
    T2 = (4*zsq*K0 + (zsq+8)*z*K1 - 8) / (pi*z6)
    T3 = (z*K1 - 1) / (2*pi*zsq)
    return (T1, T2, T3)
end

function Ti_prime(z::Float64)
    if z < POWERSERIES_ZMAX_PRIME
        (T1S, T1L, T2S, T2L, T3S, T3L) = Ti_split_pow(z, false)
        (DT1S, DT1L, DT2S, DT2L, DT3S, DT3L) = Ti_split_pow(z, true)
        T1p = DT1S + DT1L*log(z) + T1L / z
        T2p = DT2S + DT2L*log(z) + T2L/z - 2/(8*pi*z^3) + 4/(pi*z^5)
        T3p = DT3S + DT3L*log(z) + T3L/z
        return (T1p, T2p, T3p)
    else
        return Ti_direct_prime(z)
    end
end

function Ti_direct_prime(z::Float64)
    zsq = z*z
    z4 = zsq*zsq
    z6 = z4*zsq
    K0 = besselK0(z)
    K1 = besselK1(z)
    T1p = ((2*z^2 + 8)*K1)/(z^4*pi) - 8/(z^5*pi) + ((z^2 + 8)*K0)/(2*z^3*pi)
    T2p = 48/(z^7*pi) - ((8*z^2 + 48)*K1)/(z^6*pi) - ((z^2 + 24)*K0)/(z^5*pi)
    T3p = 1/(z^3*pi) - K0/(2*z*pi) - K1/(z^2*pi)
    return (T1p, T2p, T3p)
end

"""
Split of stresslet functions
# Returns: (T1S, T1L, T2S, T2L, T3S, T3L)
"""
function Ti_split(z::Float64)
    if z < POWERSERIES_ZMAX
        return Ti_split_pow(z)
    else
        return Ti_split_direct(z)
    end
end

function Ti_split_prime(z::Float64)
    if z < POWERSERIES_ZMAX_PRIME
        return Ti_split_pow(z,true)
    else
        return Ti_split_direct(z,true)
    end
end

function Ti_split_direct(z::Float64,ifprime::Bool=false)

    if !ifprime

        K0 = besselK0(z)
        K1 = besselK1(z)
        I0 = besselI0(z)
        I1 = besselI1(z)

        K0S = K0 + I0*log(z)
        K1S = K1 - I1*log(z) - 1/z

        T1S = -(2*z*K0S + (z^2+4)*K1S + z) / (2*pi*z^3)
        T1L = (2*z*I0 - (z^2+4)*I1) / (2*pi*z^3)

        T2S = (32*z*K0S + 8*(z^2+8)*K1S  - z*(z^2-16)) / (8*pi*z^5)
        T2L = ((z^2+8)*I1 - 4*z*I0) / (pi*z^5)

        T3S = K1S/(2*pi*z)
        T3L = I1/(2*pi*z)
        return (T1S, T1L, T2S, T2L, T3S, T3L)

    else

        logz = log(z)
        z2 = z*z
        z3 = z*z2
        z4 = z*z3
        z5 = z*z4
        z6 = z*z5
        z7 = z*z6

        K0 = besselK0(z)
        K1 = besselK1(z)
        I0 = besselI0(z)
        I1 = besselI1(z)

        DT1S = (-16 +
                z*(z*(8+z2)*K0+4*(4+z2)*K1-(4+z2)*I1*(4*logz-1)+
                   z*I0*((8+z2)*logz-2)))/(2*pi*z5)
        DT2S = (192-16*z2+z4+
                4*z*(-z*(24+z2)*K0-8*(6+z2)*K1+
                     I1*(-8-z2+8*(6+z2)*logz)-
                     z*I0*(-4+(24+z2)*logz)))/(4*pi*z7)
        DT3S = (2-I1*(z-2*z*logz)-
                z2*(2*K1/z+K0+I0*logz))/(2*pi*z3)

        DT1L = (-2*z*(8 + z2)*I0 + 8*(4 + z2)*I1)/(4*pi*z4)
        DT2L = (z*(24 + z2)*I0 - 8*(6 + z2)*I1)/(pi*z6)
        DT3L = (I0-2*I1/z)/(2*pi*z)

        return (DT1S, DT1L, DT2S, DT2L, DT3S, DT3L)

    end

end

"""
Direct power series for split stresslet functions
# Returns: (T1S, T1L, T2S, T2L, T3S, T3L)
"""
function Ti_split_pow(z::Float64,ifprime::Bool=false)
    if !ifprime
        # Outputs
        T1S = 0.0
        T1L = 0.0
        T2S = 0.0
        T2L = 0.0
        T3S = 0.0
        T3L = 0.0
        # Add n=0 term
        log2 = log(2)
        k1L0 = 0.5
        k1S0 = (-2*log2-(-2*eulergamma+1.0)) * 0.25
        T3S += k1S0
        T3L += k1L0
        # Recursion variables, init at n=1
        n=1
        z2 = z*z
        z2nm2 = 1; # z^(2*n-2)
        z2n = z2; # z^(2*n)
        digamma_np1 = -eulergamma + 1.0 # digamma(n+1)
        digamma_np2 = -eulergamma + 1.5
        ifact_n = 1.0 # 1/factorial(n)
        ifact_np1 = 1.0/2
        i4n = 1.0/4 # 1/4^n
        i4np1 = 1.0/16
        k1Snm1 = k1S0 # k1S(n-1)
        k1Lnm1 = k1L0 # k1L(n-1)
        ## First iteration unrolled
        # k0, k1
        k0n_denom = i4n*ifact_n*ifact_n
        k1n_denom = i4np1*ifact_n*ifact_np1
        k0Ln = (-1)*k0n_denom
        k0Sn = (digamma_np1 + log2)*k0n_denom
        k1Ln = (2)*k1n_denom
        k1Sn = (-2*log2-digamma_np1-digamma_np2)*k1n_denom
        # T1
        t1Sn = 2*k0Sn + 4*k1Sn + k1Snm1
        t1Ln = 2*k0Ln + 4*k1Ln + k1Lnm1
        T1S += t1Sn*z2nm2
        T1L += t1Ln*z2nm2
        # T3
        T3S += k1Sn*z2n
        T3L += k1Ln*z2n
        # Update recursions
        n += 1
        inp1 = 1/(n+1)
        z2n *= z2
        z2nm2 = z2
        z2nm4 = 1
        (digamma_np1, digamma_np2) = (digamma_np2, digamma_np2+inp1)
        (ifact_n, ifact_np1) = (ifact_np1, ifact_np1*inp1)
        (i4n, i4np1) = (i4np1, i4np1*0.25)
        k1Snm1 = k1Sn
        k1Lnm1 = k1Ln
        ## Loop
        while n<POWERSERIES_NMAX
            # k0, k1
            k0n_denom = i4n*ifact_n*ifact_n
            k1n_denom = i4np1*ifact_n*ifact_np1
            k0Ln = (-1)*k0n_denom
            k0Sn = (digamma_np1 + log2)*k0n_denom
            k1Ln = (2)*k1n_denom
            k1Sn = (-2*log2-digamma_np1-digamma_np2)*k1n_denom
            # T1
            t1Sn = 2*k0Sn + 4*k1Sn + k1Snm1
            t1Ln = 2*k0Ln + 4*k1Ln + k1Lnm1
            T1S += t1Sn*z2nm2
            T1L += t1Ln*z2nm2
            # T2
            t2Sn = 4*k0Sn + 8*k1Sn + k1Snm1
            t2Ln = 4*k0Ln + 8*k1Ln + k1Lnm1
            T2S += t2Sn*z2nm4
            T2L += t2Ln*z2nm4
            # T3
            T3S += k1Sn*z2n
            T3L += k1Ln*z2n
            # Update recursions
            n += 1
            inp1 = 1/(n+1)
            (z2nm4, z2nm2, z2n) = (z2nm2, z2n, z2n*z2)
            (digamma_np1, digamma_np2) = (digamma_np2, digamma_np2+inp1)
            (ifact_n, ifact_np1) = (ifact_np1, ifact_np1*inp1)
            (i4n, i4np1) = (i4np1, i4np1*0.25)
            k1Snm1 = k1Sn
            k1Lnm1 = k1Ln
        end
        # Multiply with constant factors
        T1S *= -1/(2*pi)
        T1L *= -1/(2*pi)
        T2S *= 1/pi
        T2L *= 1/pi
        T3S *= 1/(2*pi)
        T3L *= 1/(2*pi)
        return (T1S, T1L, T2S, T2L, T3S, T3L)
    else
        # Outputs
        T1S = 0.0
        T1L = 0.0
        T2S = 0.0
        T2L = 0.0
        T3S = 0.0
        T3L = 0.0
        # Add n=0 term (derivative is zero)
        log2 = log(2)
        k1L0 = 0.5
        k1S0 = (-2*log2-(-2*eulergamma+1.0)) * 0.25
        #T3S += k1S0
        #T3L += k1L0
        # Recursion variables, init at n=1
        n=1
        z2 = z*z
        z2nm3 = 0; # z^(2*n-3)*(2n-2)
        z2nm1 = z*2; # z^(2*n-1)*2n
        digamma_np1 = -eulergamma + 1.0 # digamma(n+1)
        digamma_np2 = -eulergamma + 1.5
        ifact_n = 1.0 # 1/factorial(n)
        ifact_np1 = 1.0/2
        i4n = 1.0/4 # 1/4^n
        i4np1 = 1.0/16
        k1Snm1 = k1S0 # k1S(n-1)
        k1Lnm1 = k1L0 # k1L(n-1)
        ## First iteration unrolled
        # k0, k1
        k0n_denom = i4n*ifact_n*ifact_n
        k1n_denom = i4np1*ifact_n*ifact_np1
        k0Ln = (-1)*k0n_denom
        k0Sn = (digamma_np1 + log2)*k0n_denom
        k1Ln = (2)*k1n_denom
        k1Sn = (-2*log2-digamma_np1-digamma_np2)*k1n_denom
        # T1 (constant ... derivative is zero)
        #t1Sn = 2*k0Sn + 4*k1Sn + k1Snm1
        #t1Ln = 2*k0Ln + 4*k1Ln + k1Lnm1
        #T1S += t1Sn*z2nm3
        #T1L += t1Ln*z2nm3
        # T3
        T3S += k1Sn*z2nm1
        T3L += k1Ln*z2nm1
        # Update recursions
        n += 1
        inp1 = 1/(n+1)
        z2nm1 *= z2*(n/(n-1.0))
        z2nm3 = z*2
        z2nm5 = 0
        (digamma_np1, digamma_np2) = (digamma_np2, digamma_np2+inp1)
        (ifact_n, ifact_np1) = (ifact_np1, ifact_np1*inp1)
        (i4n, i4np1) = (i4np1, i4np1*0.25)
        k1Snm1 = k1Sn
        k1Lnm1 = k1Ln
        ## Loop
        while n<POWERSERIES_NMAX_PRIME
            # k0, k1
            k0n_denom = i4n*ifact_n*ifact_n
            k1n_denom = i4np1*ifact_n*ifact_np1
            k0Ln = (-1)*k0n_denom
            k0Sn = (digamma_np1 + log2)*k0n_denom
            k1Ln = (2)*k1n_denom
            k1Sn = (-2*log2-digamma_np1-digamma_np2)*k1n_denom
            # T1
            t1Sn = 2*k0Sn + 4*k1Sn + k1Snm1
            t1Ln = 2*k0Ln + 4*k1Ln + k1Lnm1
            T1S += t1Sn*z2nm3
            T1L += t1Ln*z2nm3
            # T2
            t2Sn = 4*k0Sn + 8*k1Sn + k1Snm1
            t2Ln = 4*k0Ln + 8*k1Ln + k1Lnm1
            T2S += t2Sn*z2nm5
            T2L += t2Ln*z2nm5
            # T3
            T3S += k1Sn*z2nm1
            T3L += k1Ln*z2nm1
            # Update recursions
            n += 1
            inp1 = 1/(n+1)
            (z2nm5, z2nm3, z2nm1) = (z2nm3, z2nm1, z2nm1*z2*(n/(n-1.0)))
            (digamma_np1, digamma_np2) = (digamma_np2, digamma_np2+inp1)
            (ifact_n, ifact_np1) = (ifact_np1, ifact_np1*inp1)
            (i4n, i4np1) = (i4np1, i4np1*0.25)
            k1Snm1 = k1Sn
            k1Lnm1 = k1Ln
        end
# Multiply with constant factors
        T1S *= -1/(2*pi)
        T1L *= -1/(2*pi)
        T2S *= 1/pi
        T2L *= 1/pi
        T3S *= 1/(2*pi)
        T3L *= 1/(2*pi)
        return (T1S, T1L, T2S, T2L, T3S, T3L)
    end
end


text:
using SolverLogging
using Crayons
using Test

##
lg = SolverLogging.Logger()
SolverLogging.resetlogger!(lg)
setentry(lg, "iter", Int64, width=5)
setentry(lg, "cost")
setentry(lg, "ΔJ")
setentry(lg, "α", fmt="%6.4f")
setentry(lg, "info", String)

a = 1+2
itr = 1
@log lg "iter" itr
@log lg "cost" 10a
@log lg "ΔJ" 1e-3
cost = 12.0
@log lg cost
printheader(lg)
SolverLogging.formrow(lg)
printrow(lg)
@test lg.opts._count == 1
printrow(lg)
@test lg.opts._count == 2

iter = 2
@log lg iter
lg.data
printlog(lg)

# Test the add operation
@log lg "iter" 1 :add
@test parse(Int,lg.data[1]) == 3

setentry(lg, "info", String, width=25)
@log lg "info" "Some info"

SolverLogging.resetcount!(lg)
for iter = 1:12
    @log lg iter
    iter == 12 && @log lg "info" "last iter" :append
    printlog(lg)
end
@test occursin("Some info. last iter", SolverLogging.formrow(lg))
@test SolverLogging.isenabled(lg)

println("\nNothing should print between this line...")
SolverLogging.disable(lg)
SolverLogging.resetcount!(lg)
for iter = 1:12
    @log lg iter
    printlog(lg)
end
println("...and this line")
SolverLogging.enable(lg)
@test lg.opts.enable == true

SolverLogging.resetlogger!(lg)
@test length(lg.data) == 0

## Test Info field
setentry(lg, "info", String, width=20)
@log lg "info" "Cost Increase"
occursin("Cost Increase", lg.data[1])

# Should insert a period
@log lg "info" "Max Iters." :append
occursin("Cost Increase. Max Iters", lg.data[1])

# Try inserting spaces
@log lg "info" "   "
@log lg "info" "New message"
@test lg.data[1][1:3] == "New"

ENV["JULIA_DEBUG"] = "SolverLogging"
@test_logs (:debug,r"Rejecting") @log lg "new" 1.0

## Test output to a file
filename = joinpath(@__DIR__, "log.out")
lg = SolverLogging.Logger(filename)
setentry(lg, "iter", Int)
setentry(lg, "cost", Float64)
setentry(lg, "tol", Float64, level=2)
lg.opts.freq = 5
lg.opts.linechar = 0
iter = 1
for i = 1:10
    local iter = i
    @log lg iter
    @log lg "cost" log(10*i)
    @log lg "tol" exp(-i)
    printlog(lg)
end
SolverLogging.resetcount!(lg)
lg.opts.linechar = '*'
for i = 1:10
    local iter = i
    @log lg iter
    @log lg "cost" log(10*i)
    @log lg "tol" exp(-i)
    printlog(lg)
end
@test lg.io isa IOStream
lines = readlines(filename)
@test length(lines) == 26
for i in (1,7)
    @test occursin("iter", lines[i])
    @test occursin("cost", lines[i])
    @test !occursin("tol", lines[i])
end
for i = 2:6
    @test occursin("$(i-1)", lines[i])
end
@test length(filter(x->occursin("***", x), lines)) == 2
rm(filename)

## Test Condition formatting
ccrayon_tol = ConditionalCrayon(1e-6,Inf, reverse=false)

goodctrl = x->abs(x) < 1
badctrl = x->abs(x) > 10
defcolor = crayon"208"  # the ANSI color for a dark orange.
ccrayon_control = ConditionalCrayon(badctrl, goodctrl, defaultcolor=defcolor)

logger = SolverLogging.Logger()
setentry(logger, "iter", Int, width=5)
setentry(logger, "tol", Float64, ccrayon=ccrayon_tol)
setentry(logger, "ctrl", Float64, fmt="%.1f", ccrayon=ccrayon_control)

logger.fmt["iter"].ccrayon(10)
for iter = 1:10
    tol = exp10(-iter)
    ctrl = 0.1*(iter-1)*iter^2
    @log logger iter
    @log logger tol
    @log logger ctrl
    printlog(logger)
end


text:
import FLOWFarm; const ff = FLOWFarm

rotor_diameter = 0.57
hub_height = 0.7
yaw = 0.0
ct = 0.81
constcp = 0.8
generator_efficiency = 0.944
ai = 1.0/3.0
wind_speed = 8.1
air_density = 1.1716  # kg/m^3
ambient_ti = 0.1
turbine_x = [0.0]
turbine_y = [0.0]
turbine_z = [0.0]
turbine_yaw = [0.0]
turbine_ct = [ct]
turbine_ai = [ai]
winddirections = [270.0*pi/180.0]
windspeeds = [wind_speed]
windprobabilities = [1.0]
measurementheights = [hub_height]
shearexponent = 0.15
turbine_inflow_velcities = [wind_speed]
ambient_tis = [ambient_ti]
nturbines = 1

ct_model = ff.ThrustModelConstantCt(ct)
power_model = ff.PowerModelConstantCp(constcp)
power_models = Vector{typeof(power_model)}(undef, nturbines)
for i = 1:nturbines
    power_models[i] = power_model
end
wind_shear_model = [ff.PowerLawWindShear(shearexponent)]

cut_in_speed = 0.0
cut_out_speed = 25.0
rated_speed = 12.0
rated_power = 1.0176371581904552e6

turbine1 = ff.TurbineDefinition(1, [rotor_diameter], [hub_height], [cut_in_speed], [rated_speed], [cut_out_speed], [rated_power], [generator_efficiency], ct_model, power_model)
turbine_definitions = [turbine1]
turbine_definition_ids = [1]
sorted_turbine_index = [1]

windfarm = ff.WindFarm(turbine_x, turbine_y, turbine_z, turbine_definition_ids, turbine_definitions)
windfarmstate = ff.SingleWindFarmState(1, turbine_x, turbine_y, turbine_z, turbine_yaw, turbine_ct, turbine_ai, turbine_inflow_velcities, [0.0], [ambient_ti],sorted_turbine_index)
windresource = ff.DiscretizedWindResource(winddirections, windspeeds, windprobabilities, measurementheights, air_density, ambient_tis, wind_shear_model)

loc = [7.0*rotor_diameter, 0.0, hub_height]
alpha = 0.1
wakedeficitmodel = ff.JensenTopHat(alpha)
horizontal_spread_rate = alpha
wakedeflectionmodel = ff.JiminezYawDeflection(horizontal_spread_rate)
wakecombinationmodel = ff.SumOfSquaresFreestreamSuperposition()
localtimodel = ff.LocalTIModelNoLocalTI()

ms1 = ff.WindFarmModelSet(wakedeficitmodel, wakedeflectionmodel, wakecombinationmodel, localtimodel)
# pd1 = ff.WindFarmProblemDescription(windfarm, windresource, [windfarmstate])


text:
@testset "227.basic-calculator-ii.jl" begin
    @test calculate_ii("3+2*2") == 7
    @test calculate_ii("14-3/2") == 13
    @test calculate_ii("3/2") == 1
    @test calculate_ii("3+ 5 / 2") == 5
end


text:
module Sentiment

using ZipFile
using ..Data: deps

function load()
  isfile(deps("sentiment.zip")) || return
  info("Downloading sentiment treebank dataset")
  download("https://cache.julialang.org/https://nlp.stanford.edu/sentiment/trainDevTestTrees_PTB.zip",
           deps("sentiment.zip"))
end

getfile(r, name) = r.files[findfirst(x -> x.name == name, r.files)]

function getfile(name)
  r = ZipFile.Reader(deps("sentiment.zip"))
  text = readstring(getfile(r, "trees/$name"))
  close(r)
  return text
end

using ..Data: Tree

totree_(n, w) = Tree{Any}((parse(Int, n), w))
totree_(n, a, b) = Tree{Any}((parse(Int, n), nothing), totree(a), totree(b))
totree(t::Expr) = totree_(t.args...)

function parsetree(s)
  s = replace(s, r"\$", s -> "\\\$")
  s = replace(s, r"[^\s\(\)]+", s -> "\"$s\"")
  s = replace(s, " ", ", ")
  return totree(parse(s))
end

function gettrees(name)
  load()
  ss = split(getfile("$name.txt"), '\n', keep = false)
  return parsetree.(ss)
end

train() = gettrees("train")
test() = gettrees("test")
dev() = gettrees("dev")

end


text:
import Base.length

abstract type Cone end

"""
# Description
Type representing a finitely generated cone. This is the positive
hull of a finite collection of vectors a₁, … , aₙ:

{ ∑ λᵢ xᵢ : λ ≥ 0 }

# Arguments
`A::Matrix{Float64}`: Matrix of cone generators, where each column corresponds to a generator
"""
struct FiniteCone <: Cone
    A::Matrix{Float64}    # Matrix where each column is a cone generator
    AtA::Matrix{Float64}  # Cross-product of cone generator matrix
    num_gen::Int64
    FiniteCone(A::Matrix{T}) where T<:Real = new(float(A), float(A'A), size(A,2))
end

function project(cone::FiniteCone, p::Vector{Float64})
    w,z = lcp_solve(cone.AtA, -cone.A'p)
    cone.A * z
end

# Dimension of ambient space
length(cone::FiniteCone) = size(cone.A, 1)

"""
# Description
Type representing a Polyhedra cone. This is the
intersection of a finite collection of half-spaces

{x : aᵢ x ≥ 0 ∀ i} = { x: Ax ≥ 0 }

# Arguments
`A::Matrix{Float64}`: Constraint matrix of polyhedral cone
"""
struct PolyhedralCone{T<:Real} <: Cone
    A::Matrix{T}
    n::Int   # Dimension
    m::Int   # Number of constraints
    function PolyhedralCone{T}(A::Matrix{T}) where T<: Real
        return new{T}(A, size(A,2), size(A,1))
    end
end

function PolyhedralCone(A::Matrix{T}, check=false) where T<:Integer
    if check
        A = remove_redundant_constraints(A)
    end
    PolyhedralCone{T}(A)
end

function project(cone::PolyhedralCone, p::Vector{Float64})
    res = quadprog(-2.0*p, 2.0*Array{Float64}(I, cone.n, cone.n), cone.A,  '>',
                   zeros(cone.m), fill(-Inf,cone.n), fill(Inf, cone.n),
                   GurobiSolver(OutputFlag=0, Threads=6))
    return res.sol
end

# Dimension of ambient space
length(cone::PolyhedralCone) = size(cone.A, 2)


text:
platforms_and_sources = Dict(
    Platform("x86_64", "linux") => [
        ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.2.2/10.2_07062021/cudnn-10.2-linux-x64-v8.2.2.26.tgz",
                      "b4a2067774f509e65a1d8ba3bd86162b9e09de5946bb636887b4cc605dddeb6e")],
    Platform("x86_64", "windows") => [
        ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.2.2/10.2_07062021/cudnn-10.2-windows10-x64-v8.2.2.26.zip",
                      "27dee1094c49f4993754ddecd1e3bf4a14b1d4d43951cc5956d2375eb3bcb9f8")],
)

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
using Base.Test
using MIToS
using MIToS.Utils
using MIToS.MSA
using MIToS.Information
using MIToS.PDB
using MIToS.SIFTS
using MIToS.Pfam
using PairwiseListMatrices

# Utils
include("utils.jl")
include("eachlinestring.jl")
# MSA
include("residues.jl")
include("indexedarrays.jl")
include("annotations.jl")
include("rawalnandgaps.jl")
include("multiplesequencealignment.jl")
include("msaannotations.jl")
include("shuffle.jl")
# Clustering
include("clustering.jl")
# PDB
include("pdb.jl")
include("contacts.jl")
include("kabsch.jl")
#  SIFTS
include("sifts.jl")
# Information
include("probabilities.jl")
include("informationmeasures.jl")
include("iterations.jl")
include("buslje09.jl")
# Pfam
include("pfam.jl")
# Scripts
include("scripts.jl")

print("""

----- =D -----

""")


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
using Roots

function withcounter(f, n = 0)
    c = Ref(n)
    F(x...; y...) = (c[] += 1; f(x...; y...))
    getcounter() = c[]
    setcounter!(n) = c[] = n
    F, getcounter, setcounter!
end

f(x::Float64)::Float64 = @ccall sin(x::Float64)::Float64
@show f(1.0);

# %%
F, getcounter, setcounter! = withcounter(x -> f(x) - 0.8414709848078965)
@show find_zero(F, 0.9)
@show getcounter();

# %%
setcounter!(0)
@show find_zeros(F, 0, 2π)
@show getcounter();

# %%
dump(F)

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
using Documenter, TaylorSeries

makedocs(
    modules = [TaylorSeries],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "TaylorSeries.jl",
    authors  = "Luis Benet and David P. Sanders",
    pages = [
        "Home" => "index.md",
        "Background" => "background.md",
        "User guide" => "userguide.md",
        "Examples" => "examples.md",
        "API" => "api.md"
    ]
)

deploydocs(
    repo   = "github.com/JuliaDiff/TaylorSeries.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing
)


text:
using MCMCChains, Test

@testset "Array constructor tests" begin
    @testset "Smoke tests" begin
        val = convert(Array{Union{Missing, Float64},3}, rand(1000, 5, 5))
        chns = Chains(val, ["a", "b", "c", "d", "e"], [:internals => ["d", "e"]])

        # Config.
        main_params = 3
        internal_params = 2
        total_params = main_params + internal_params

        d, p, c = size(chns)

        @test size(Array(chns)) == (d*c, main_params)
        @test size(Array(chns, [:parameters])) == (d*c, main_params)
        @test size(Array(chns, [:parameters])) == size(Array(chns))
        @test size(Array(chns, [:parameters, :internals])) == (d*c, total_params)
        @test size(Array(chns, [:parameters, :internals])) ==
            size(Array(chns, showall=true))
        @test size(Array(chns, [:internals])) == (d*c, internal_params)
        @test size(Array(chns, append_chains=true)) == (d*c, main_params)
        @test size(Array(chns, append_chains=false)) == (5,)
        @test size(Array(chns, append_chains=false)[1]) == (d, main_params)
        @test typeof(Array(chns, append_chains=true)) == Array{Float64, 2}
        @test size(Array(chns, remove_missing_union=false)) == (d*c, main_params)
        @test size(Array(chns, append_chains=false,
            remove_missing_union=false)) == (5,)
        @test size(Array(chns, append_chains=false,
            remove_missing_union=false)[1]) == (d, main_params)
        @test typeof(Array(chns, append_chains=true, remove_missing_union=false)) ==
            Array{Union{Missing, Float64}, 2}
        @test size(Array(chns[:b])) == (d*c,)

        Array(chns)
        Array(chns[:a])
        Array(chns, [:parameters])
        Array(chns, [:parameters, :internals])
    end
    @testset "Accuracy" begin
        nchains = 5
        val = rand(500, 5, nchains)
        long_val = vcat([val[:,:,i] for i in 1:nchains]...)
        squished_val = [val[:,:,i] for i in 1:nchains]
        chn = Chains(val)

        # Test no_append
        arr1 = Array(chn, append_chains=true)
        @test all(arr1 .== long_val)

        # Test append
        arr2 = Array(chn, append_chains=false)
        @test all(vcat([arr2[i] .== squished_val[i] for i in 1:nchains]...))
    end
end


text:
a = 1:10
b = 1:10
@test PearsonCorrelation(a, b) == 1

x = [12,15,17,21,18]
y = [6,9,-1,3,1]

@test PearsonCorrelation(x,y) ==  -0.5201361507569889

max_stars_repo_name: arubhardwaj/Julia

max_issues_repo_name: arubhardwaj/Julia

max_forks_repo_name: arubhardwaj/Julia


text:
@testset "97.interleaving-string.jl" begin
    @test is_interleave("aabcc", "dbbca", "aadbbbaccc") == false
    @test is_interleave("aabcc", "dbbca", "aadbbcbcac") == true
    @test is_interleave("", "", "") == true
end


text:
# jl1 aggregate_results.jl /scratch/midway2/mkolar/diffTest/sim8/res_1 our_res_1

using DiffPrecTest
using JLD

resPath = ARGS[1]
outPrefix = ARGS[2]

NUM_REP = 1000

results = Vector{BootstrapEstimates}(undef, NUM_REP)
for rep=1:NUM_REP
    global results
    fname = "$(resPath)_$(rep).jld"

    try
      file = jldopen(fname, "r")
      res = read(file, "boot_res")
      close(file)

      results[rep] = deepcopy(res)
    catch
        @show fname
    end
end

@save "$(outPrefix).jld" results


text:
using FlightSims
const FS = FlightSims
using DifferentialEquations
using LinearAlgebra
using Plots
using Test
using Transducers

function main()
    # linear system
    A = [0 1;
         0 0]  # 2 x 2
    B = [0 1]'  # 2 x 1
    n, m = 2, 1
    env = LinearSystemEnv(A, B)  # exported from FlightSims
    x0 = State(env)([1.0, 2.0])
    p0 = zero.(x0)  # auxiliary parameter
    # optimal control
    Q = Matrix(I, n, n)
    R = Matrix(I, m, m)
    lqr = LQR(A, B, Q, R)  # exported from FlightSims
    u_lqr = Command(lqr)  # (x, p, t) -> -K*x; minimise J = ∫ (x' Q x + u' R u) from 0 to ∞

    # simulation
    tf = 10.0
    Δt = 0.01
    affect!(integrator) = integrator.p = copy(integrator.u)  # auxiliary callback funciton
    cb = PeriodicCallback(affect!, Δt; initial_affect=true)  # auxiliary callback
    @Loggable function dynamics!(dx, x, p, t)
        @onlylog p  # activate this line only when logging data
        u = u_lqr(x)
        @log x, u
        @nested_log Dynamics!(env)(dx, x, p, t; u=u)  # exported `state` and `input` from `Dynamics!(env)`
    end
    prob, df = sim(
                   x0,  # initial condition
                   dynamics!,  # dynamics with input of LQR
                   p0;
                   tf=tf,  # final time
                   callback=cb,
                   savestep=Δt,
                  )
    ts = df.time
    xs = df.sol |> Map(datum -> datum.x) |> collect
    us = df.sol |> Map(datum -> datum.u) |> collect
    ps = df.sol |> Map(datum -> datum.p) |> collect
    states = df.sol |> Map(datum -> datum.state) |> collect
    inputs = df.sol |> Map(datum -> datum.input) |> collect
    @test xs == states
    @test us == inputs
    p_x = plot(ts, hcat(states...)';
               title="state variable", label=["x1" "x2"], color=[:black :black], lw=1.5,
              )  # Plots
    plot!(p_x, ts, hcat(ps...)';
          ls=:dash, label="param", color=[:red :orange], lw=1.5
         )
    savefig("figures/x_lqr.png")
    plot(ts, hcat(inputs...)'; title="control input", label="u")  # Plots
    savefig("figures/u_lqr.png")
    df
end


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

using FEMBase, MortarContact2DAD, Test, LinearAlgebra, SparseArrays, Statistics,
      DelimitedFiles

include(joinpath("..", "docs", "make.jl"))

@testset "MortarContact2DAD.jl" begin
    @testset "test_01.jl" begin include("test_01.jl") end
    @testset "test_02.jl" begin include("test_02.jl") end
    @testset "test_contact_1.jl" begin include("test_contact_1.jl") end
    @testset "test_contact_2.jl" begin include("test_contact_2.jl") end
    @testset "test_contact_projection.jl" begin include("test_contact_projection.jl") end
end

include(joinpath("..", "docs", "deploy.jl"))


text:
## Goals:
# For the flow-clustering paper, a referee asked about how these ideas could
# be used on images.
# Here, we provide an example of how they can be used on images.
# This is based on the mqi-animation figures I made for the Argonne seminar.
# I copied "code" from that directory
# cp -r  ~/Dropbox/publications/flow-based-argonne/code .

# Originally based on dgleich/Dropbox/research/2021/09-22-mqi-figs-for-image

# Use Nate's flowseed code. (Push Relabel with customizable seeds... )
include("FlowSeed.jl")

## Setup code to work on a grid graph...
using SparseArrays

function grid_graph(m::Int, n::Int; distance::Int=1)
  N = m*n
  imap = reshape(1:N, m, n)
  ei = zeros(Int,0)
  ej = zeros(Int,0)
  for i=1:m
    for j=1:n
      # left-right neighbors
      for di=-distance:distance
        dj = 0
        if i+di >= 1 && i+di <= m && j+dj >= 1 && j+dj <= n
          src = imap[i,j]
          dst = imap[i+di,j+dj]
          if (src != dst)
            push!(ei,src)
            push!(ej,dst)
          end
        end
      end
      # up-down neighbors
      for dj=-distance:distance
        di = 0
        if i+di >= 1 && i+di <= m && j+dj >= 1 && j+dj <= n
          src = imap[i,j]
          dst = imap[i+di,j+dj]
          if (src != dst)
            push!(ei,src)
            push!(ej,dst)
          end
        end
      end
    end
  end
  return sparse(ei,ej,1.0,N,N)
end

A = grid_graph(100,100)
## Now generate a blurry image on a grid graph.
using Plots
function make_T(m::Integer,n::Integer)
  IType = typeof(m)

  # start of the left-right branch
  lstart = floor(IType,m/8)
  rend = m-lstart
  # end of the left-right branch
  tstart = floor(IType, n/8)
  bend = n-tstart

  halflstart = div(lstart,2)

  mstart = floor(IType,m/2)-lstart
  mend = floor(IType,m/2)+lstart

  X = zeros(m,n)
  X[lstart:rend,tstart:(tstart+lstart)] .= 1
  X[mstart:mend,tstart:bend] .= 1
  return X
end
X = make_T(100,100)
heatmap((make_T(100,100)))
## blur this picture.
using Images, Random
Random.seed!(0)
T = imfilter(X, Kernel.moffat(1.5, 1.2, 5))
Tblur = copy(T)
T .+= 0.1*(2*rand(size(T)...).-1)
T = max.(min.(T/0.3,1.0), 0.0)
heatmap(T)
## Solve one step of MQI-like problem on T
# This weights the problem with a volume dependent term
# That biases the results towards the T solution.
using PerceptualColourMaps
function MQITest(A,x,delta)
  volA = sum(A)
  d = vec(sum(A,dims=1))

  R = findall(x .> 0)
  #source = zeros(size(A,1))
  source = delta*d.*x # build edge weights from source to R

  sink = volA*ones(size(A,1)) # build edge weights from Rbar to sink
  sink[R] .= 0 # no edges from R
  #sstart = set_stats(A,R,volA) # get conductance
  #println("starting cond=", sstart[4], " = ", sstart[1], "/", sstart[2])
  S = NonLocalPushRelabel(A,R,source,sink)
  if length(S) > 0
    send = set_stats(A,S,volA)
    println("  ending cond=", send[4], " = ", send[1], "/", send[2])
  else
    send = sstart
    S = R
  end
  return S, send
end
S = MQITest(A,reshape(T,100*100),0.11)[1]
SMat = begin M =zeros(100,100); M[S] .= 1; M end
SMqi = begin S = MQITest(A,reshape(T,100*100),0.04)[1]; M =zeros(100,100); M[S] .= 1; M end
plot(heatmap(1 .-X, title="original",color=cmap("L3")),
  heatmap(1 .-T,title="noisy/blury",color=cmap("L3")),
  heatmap(1 .-(0.5*SMat+0.5*X),color=cmap("L3"),title="MQI like\\n subprob (\\delta=0.11)"),
  heatmap(1 .-(0.5*SMqi+0.5*X),color=cmap("L3"),title="MQI like\\n opt (\\delta=0.04)"), xticks=[],yticks=[],
    colorbar=false, aspect_ratio=:equal, framestyle=:none, size=(600,200),
    layout=(1,4),titlefontsize=12,color=cmap("L3"))

## Save final figures... for manuscript
using PerceptualColourMaps
colorview(RGB,permutedims(applycolourmap(1 .- X, cmap("L3")), (3,1,2)))
##
save("image-1-original.png",
  map(clamp01nan, colorview(RGB,permutedims(applycolourmap(1 .- X, cmap("L3")), (3,1,2))) ))
##
save("image-1-blur.png",
  map(clamp01nan, colorview(RGB,permutedims(applycolourmap(1 .- Tblur, cmap("L3")), (3,1,2))) ))

##
save("image-1-noisy.png",
  map(clamp01nan, colorview(RGB,permutedims(applycolourmap(1 .- T, cmap("L3")), (3,1,2))) ))
##
save("image-1-colormap.png",repeat(reverse(cmap("L3")), 1, 30))
save("image-1-colormap-horiz.png",repeat(reverse(cmap("L3")), 1, 30)')

##
save("image-1-MQI-like-1.png",
  map(clamp01nan, colorview(RGB,permutedims(applycolourmap(1 .- SMat, cmap("L3")), (3,1,2))) ))

##
save("image-1-MQI-opt-1.png",
  map(clamp01nan, colorview(RGB,permutedims(applycolourmap(1 .- SMqi, cmap("L3")), (3,1,2))) ))

##
save("image-1-MQI-like-1-errors.png",
  map(clamp01nan, colorview(RGB,permutedims(applycolourmap(
    1 .- 0.5*xor.(Int.(SMat),Int.(X)), cmap("L3"), [0,1]), (3,1,2))) ))
##
save("image-1-MQI-like-1-errors-and-soln.png",
  map(clamp01nan, colorview(RGB,permutedims(applycolourmap(
    1 .- 0.5.*SMat .- 0.5.*X , cmap("L3"), [0,1]), (3,1,2))) ))

##
save("image-1-MQI-opt-1-errors.png",
  map(clamp01nan, colorview(RGB,permutedims(applycolourmap(
    1 .- 0.5*xor.(Int.(SMqi),Int.(X)), cmap("L3"), [0,1]), (3,1,2))) ))

max_stars_repo_name: dgleich/flowpaper-code

max_issues_repo_name: dgleich/flowpaper-code

max_forks_repo_name: dgleich/flowpaper-code


text:
using TextTableRows
using Base.Test

import TextTableRows: perc, count_and_perc

@testset "perc" begin
    @test perc(pi) == "314.2%"
    @test perc(1/3) == "33.3%"
    @test count_and_perc(10, 20) == "10  (50.0%)"
end

@testset "Extrema" begin
    s = Extrema(Int)
    @test isempty(s)
    @test length(s) == 0
    push!(s, 1)
    @test !isempty(s)
    @test minimum(s) == maximum(s) == length(s) == 1
    push!(s, -7)
    @test extrema(s) == (-7,1)
    @test length(s) == 2
    @test repr(s) == "Extrema{$(Int),$(Int)}: 2 values, -7..1"
end

@testset "UniqueStrings" begin
    a = "foo"
    b = "bar"
    s = UniqueStrings()
    @test length(s) == 0
    @test isempty(s)
    push!(s, a)
    @test length(s) == 1
    @test !isempty(s)
    for _ in 1:9 push!(s, a) end
    for _ in 1:30 push!(s, b) end
    @test length(s) == 40
    @test repr(s) == """UniqueStrings{Int64}: 40 values
    bar => 30  (75.0%)
    foo => 10  (25.0%)
"""
end

@testset "MultipleSummaries" begin
    s = MultipleSummaries((UniqueStrings(), Extrema(Int)))
    s1 = UniqueStrings()
    s2 = Extrema(Int)
    @test length(s) == 0
    @test isempty(s)
    strings1 = ["foo", "bar", "baz"]
    for _ in 1:rand(100:200)
        x = rand(strings1)
        push!(s, 1, x)
        push!(s1, x)
    end
    # for _ in 1:rand(100:200)
    #     x = rand(1:10000)
    #     push!(s, 2, x)
    #     push!(s2, x)
    # end
    # @test length(s) == length(s1) + length(s2)
end


text:
function isPrime(n::Int)
	for i in 2:ceil(sqrt(n))
		if n%i==0 && n!=i return false end
	end
	return true
end
primes = [1,2,3]

ubound = 2000000
for i in 5:2:ubound
	if isPrime(i)
		push!(primes, i)
	end
end

println(sum(primes))

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
reg_function = Union{Function,Dict{Vector{Int8},Int8},Dict{BitArray{1},Bool}}
state_type = Union{Bool,Int8}
sim_type = Union{BitArray{2},Array{Int8,2},Array{Bool,2}}
mutable struct BNode
    id::AbstractString               # Name
    state::Vector{Int8}              # Current state
    def_state::Int8              # Default state
    state_range::Vector{Int8}    # Possible states
    sigma_in::Vector{Int16}
    sigma_in_state::Vector{Int8}
    sigma_in_state_mu::Vector{Float16}
    sigma_out::Vector{Int}
    F::Dict{Tuple{Vararg{Int8}},Int8}
    tau_in::Vector{UInt16}
    tau_out::UInt16
    thresholds::Vector{Float16}
    function BNode()
        this = new()
        this.state = [false]
        this.def_state = false
        this.state_range = [false,true]
        this.sigma_in = Vector{Int16}()
        this.sigma_in_state = Vector{Int8}()
        this.sigma_in_state_mu = Vector{Float16}()
        this.sigma_out = Vector{Int}()
        this.F = Dict{Tuple{Vararg{Int8}},Int8}()
        this.tau_in = Vector{UInt16}()
        this.tau_out = 1
        this.thresholds = Vector{Float16}()
        return this
    end
    function BNode(id::AbstractString)
        this = new()
        this.id = id
        this.state = [false]
        this.def_state = false
        this.state_range = [false,true]
        this.sigma_in = Vector{Int16}()
        this.sigma_in_state = Vector{Int8}()
        this.sigma_in_state_mu = Vector{Float16}()
        this.sigma_out = Vector{Int}()
        this.F = Dict{Tuple{Vararg{Int8}},Int8}()
        this.tau_in = Vector{UInt16}()
        this.tau_out = 1
        this.thresholds = Vector{Float16}()
        return this
    end
    #Node(id::AbstractString,state::Integer,def_state::Integer,state_range::Vector{Integer},sigma_in::Vector{Node})=new(id,state,def_state,state_range,sigma_in,Dict{Tuple{Vararg{Int8}},Integer}())
    #Node(id::AbstractString,state::Integer,def_state::Integer,state_range::Vector{Integer},sigma_in::Vector{Node},F::Dict{Tuple{Vararg{Int8}},Integer})=new(id,state,def_state,state_range,sigma_in,F)
end

# type BNode{T<:Union{Int8,AbstractString}}
#     id::T                           # Name
#     state::Int8                  # Current state
#     def_state::Int8              # Default state
#     state_range::Vector{Int8}    # Possible states
#     sigma_in::Vector{BNode{T}}
#     F::Dict{Tuple{Vararg{Int8}},Int8}
#     BNode() = new(false,false,[false,true],Vector{BNode{T}}(),Dict{Tuple{Vararg{Integer}},Integer}())
# end

mutable struct Net
    state::Vector{BNode}
    graph::DiGraph
    Net(state::Vector{BNode},graph::DiGraph) = new(state,graph)
end

mutable struct BNode2{T<:state_type}
    id::AbstractString               # Name
    def_state::T             # Default state
    state_range::Vector{T}    # Possible states
    rule::reg_function
    thresholds::Vector{Float16}
    sync_order::Int
    function BNode2{T}(id::AbstractString,def_state::Union{Int,Bool},state_range::Union{Vector{Int},Vector{Bool}}) where {T}
        vectorType = T==Bool ? BitArray{1} : Vector{T}
        this = new()
        this.id = id
        this.def_state = def_state
        this.state_range = state_range
        this.rule = Dict{vectorType,T}()::reg_function
        this.thresholds = Vector{Float16}()
        this.sync_order = 1
        return this
    end
    function BNode2{T}(id::AbstractString,def_state::Union{Int,Bool},state_range::Union{Vector{Int},Vector{Bool}},sync_order::Int) where {T}
        vectorType = T==Bool ? BitArray{1} : Vector{T}
        this = new()
        this.id = id
        this.def_state = def_state
        this.state_range = state_range
        this.rule = Dict{vectorType,T}()::reg_function
        this.thresholds = Vector{Float16}()
        this.sync_order = sync_order
        return this
    end
end

mutable struct Net2{T<:state_type}
    nodes::Vector{BNode2{T}}
    graph::DiGraph
    in_taus::SparseMatrixCSC{UInt16,Int}
    in_states::SparseMatrixCSC{T,Int}
    mu_in_states::SparseMatrixCSC{Float16,Int}
    sim_buffer::Union{BitArray{2},Array{T,2}}
    Net2{T}(nodes::Vector{BNode2{T}},graph::DiGraph) where {T<:state_type} = new(nodes,graph,spzeros(nv(graph),nv(graph)),adjacency_matrix(graph),adjacency_matrix(graph),T==Bool ? falses(nv(graph),1) : zeros(T,nv(graph),1))
    Net2{T}(nodes::Vector{BNode2{T}},graph::DiGraph,in_taus::SparseMatrixCSC{Int,Int}) where {T<:state_type}= (nv(graph),nv(graph)) == size(in_taus) ? new(nodes,graph,in_taus,adjacency_matrix(graph),adjacency_matrix(graph),T==Bool ? falses(nv(graph),maximum(in_taus)+1) : zeros(T,nv(graph),maximum(in_taus)+1)) : error("tau_matrix must have the same dimensions as adjacency matrix of network")
end

    #Node(id::AbstractString,state::Integer,def_state::Integer,state_range::Vector{Integer},sigma_in::Vector{Node})=new(id,state,def_state,state_range,sigma_in,Dict{Tuple{Vararg{Int8}},Integer}())
    #Node(id::AbstractString,state::Integer,def_state::Integer,state_range::Vector{Integer},sigma_in::Vector{Node},F::Dict{Tuple{Vararg{Int8}},Integer})=new(id,state,def_state,state_range,sigma_in,F)


text:
# Initialization of some simple polygons

function rectangle(x0::T, y0::T, x1::T, y1::T) where T
    bottomleft_corner = Point{T}(x0, y0)
    bottomright_corner = Point{T}(x1, y0)
    topright_corner = Point{T}(x1, y1)
    topleft_corner = Point{T}(x0, y1)
    ConvexPolygon{T}([bottomleft_corner, topleft_corner, topright_corner, bottomright_corner])
end
rectangle(x0, y0, x1, y1) = rectangle(promote(x0, y0, x1, y1)...)
rectangle(bottom_left, top_right) = rectangle(bottom_left[1], bottom_left[2], top_right[1], top_right[2])
rectangle(; bottom_left, top_right) = rectangle(bottom_left, top_right)

square(bottom_left, side::Number) = rectangle(bottom_left[1], bottom_left[2], bottom_left[1] + side, bottom_left[2] + side)
square(; bottom_left, side) = square(bottom_left, side)

function circle(x0::T, y0::T, radius::T, angle_range::AbstractVector) where T
    vertices = [Point{T}(x0 + radius*cos(θ), y0 + radius*sin(θ)) for θ in reverse(angle_range)]
    push!(vertices, Point{T}(x0, y0))
    return ConvexPolygon(vertices)
end

function circle(x0::T, y0::T, radius::T, nb_vertices::Int) where T
    vertices = [Point{T}(x0 + radius*cos(θ), y0 + radius*sin(θ)) for θ in reverse(LinRange(0.0, 2π, nb_vertices+1)[1:nb_vertices])]
    # Reversed, such that the inner of the domain is on the right of the interface
    return ConvexPolygon(vertices)
end

circle(x0, y0, r, θ_range) = circle(promote(x0, y0, r)..., θ_range)
circle(center::Union{Point, Tuple}, r, θ_range) = circle(center[1], center[2], r, θ_range)

function circle(center::Point, point::Point, θ_range)
    radius = sqrt((point1[1] - center[1])^2 + (point1[2] - center[2])^2)
    return circle(center, radius, θ_range)
end

function _circle_center(point1::Point{T}, point2::Point{T}, curvature::T) where T
    radius = one(T)/curvature
    distance_between_points = sqrt((point2[1] - point1[1])^2 + (point2[2] - point1[2])^2)
    t = (point2 - point1) / distance_between_points
    on_median = sqrt(radius^2 - (distance_between_points/2)^2)
    n = SVector(t[2], -t[1])
    return radius, point1 + distance_between_points/2 * t + on_median * n
end

function circle(point1::Point{T}, point2::Point{T}, curvature::T, θ_range) where T <: Real
    radius, center = _circle_center(point1, point2, curvature)
    return circle(center, radius, θ_range)
end


text: include("/tmp/exp_m6_mono_taylor_2_7.jl");

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
###
# Reading an archive
#
# Basic outline for reading an archive in C:
#     1. Ask `archive_read_new` for an archive reader object.
#     2. Update any global properties as appropriate.
#        In particular, you'll certainly want to call appropriate
#        `archive_read_support_XXX` functions.
#     3. Call `archive_read_open_XXX` to open the archive
#     4. Repeatedly call `archive_read_next_header` to get information about
#        successive archive entries.  Call `archive_read_data` to extract
#        data for entries of interest.
#     5. Call `archive_read_free` to end processing.
#
# Basic outline for reading an archive in Julia:
#     1. Create a LibArchive.Reader
#     2. Update any global properties as appropriate.
#        In particular, you'll certainly want to call appropriate
#        LibArchive.support_XXX functions.
#     3. Read the entries (API TBD). LibArchive.jl will call libarchive
#        open functions if needed.
#     4. Call LibArchive.free to end processing.

abstract type ReaderData end

mutable struct Reader{T<:ReaderData} <: Archive
    data::T
    ptr::Ptr{Cvoid}
    opened::Bool
    passphrase_mgr
    function Reader{T}(data::T) where T
        ptr = ccall((:archive_read_new, libarchive), Ptr{Cvoid}, ())
        ptr == C_NULL && throw(OutOfMemoryError())
        return Reader{T}(data, ptr, false)
    end
    # For pre-openned archive. (e.g. from LibALPM)
    function Reader{T}(data::T, ptr::Ptr{Cvoid}, opened::Bool) where T
        obj = new(data, ptr, opened)
        finalizer(obj, free)
        obj
    end
end

Reader(data::T) where {T<:ReaderData} = Reader{T}(data)

function free(archive::Reader)
    ptr = archive.ptr
    ptr == C_NULL && return
    ccall((:archive_read_free, libarchive), Cint, (Ptr{Cvoid},), ptr)
    archive.ptr = C_NULL
    nothing
end

Base.close(archive::Reader) =
    @_la_call(archive_read_close, (Ptr{Cvoid},), archive)

function Base.cconvert(::Type{Ptr{Cvoid}}, archive::Reader)
    archive.ptr == C_NULL && error("archive already freed")
    archive
end
Base.unsafe_convert(::Type{Ptr{Cvoid}}, archive::Reader) = archive.ptr

function ensure_open(archive::Reader)
    archive.opened && return
    archive.opened = true
    try
        do_open(archive)
    catch
        archive.opened = false
        rethrow()
    end
    nothing
end

###
# Open filename

struct ReadFileName{T} <: ReaderData
    name::T
    block_size::Int
end

Reader(fname::AbstractString; block_size=10240) =
    Reader(ReadFileName(fname, Int(block_size)::Int))
Reader(; block_size=10240) =
    Reader(ReadFileName(nothing, Int(block_size)::Int))

function do_open(archive::Reader{ReadFileName{T}}) where T
    data = archive.data
    @_la_call(archive_read_open_filename, (Ptr{Cvoid}, Cstring, Csize_t),
              archive, data.name, data.block_size)
end

function do_open(archive::Reader{ReadFileName{Cvoid}})
    data = archive.data
    @_la_call(archive_read_open_filename, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t),
              archive, C_NULL, data.block_size)
end

struct ReadFD <: ReaderData
    fd::Cint
    block_size::Int
end

Reader(fd::Integer; block_size=10240) =
    Reader(ReadFD(Cint(fd), Int(block_size)::Int))

function do_open(archive::Reader{ReadFD})
    data = archive.data
    @_la_call(archive_read_open_fd, (Ptr{Cvoid}, Cint, Csize_t),
              archive, data.fd, data.block_size)
end

struct ReadFILE <: ReaderData
    file::Libc.FILE
end

Reader(file::Libc.FILE) = Reader(ReadFILE(file))

function do_open(archive::Reader{ReadFILE})
    data = archive.data
    @_la_call(archive_read_open_FILE, (Ptr{Cvoid}, Ptr{Cvoid}), archive, data.file)
end

struct ReadFileNames{T} <: ReaderData
    names::Vector{T}
    block_size::Int
end

"""
Use this for reading multivolume files by filenames.
NOTE: Sorting is NOT done.
"""
Reader(fnames::AbstractVector{T} where T<:AbstractString; block_size=10240) =
    Reader(ReadFileNames{T}(Vector(fnames), Int(block_size)::Int))

function do_open(archive::Reader{ReadFileNames{T}}) where T
    data = archive.data
    @_la_call(archive_read_open_filenames, (Ptr{Cvoid}, Ptr{Cstring}, Csize_t),
              archive, data.names, data.block_size)
end

###
# Open memory

struct ReadMemory{T,TO} <: ReaderData
    ptr::T
    obj::TO # Reference for GC
    size::Int
end

Reader(ptr::T, size, obj::TO=nothing) where {T<:Ptr,TO} =
    Reader(ReadMemory{T,TO}(ptr, obj, size))
Reader(ary::Vector, size=sizeof(ary)) = Reader(pointer(ary), size, ary)

Base.unsafe_convert(::Type{Ptr{Cvoid}}, mem_data::ReadMemory) = mem_data.ptr

function do_open(archive::Reader{T}) where {T<:ReadMemory}
    data = archive.data
    @_la_call(archive_read_open_memory, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t),
              archive, data, data.size)
end

###
# Generic reader

struct GenericReadData{T} <: ReaderData
    data::T
    buff::Vector{UInt8}
end

Reader(data::T, buff_size=10240) where {T} =
    Reader(GenericReadData{T}(data, Vector{UInt8}(undef, buff_size)))

reader_open(archive::Reader, data) = nothing
function reader_readbytes end
function reader_skip end
function reader_seek end
reader_close(archive::Reader, data) = nothing

function reader_open_callback(c_archive, archive)
    try
        clear_error(archive)
        reader_open(archive, archive.data.data)
        return Libc.errno(archive) == 0 ? Cint(0) : Status.WARN
    catch ex
        return set_exception(archive, ex)
    end
end

function reader_read_callback(c_archive, archive, buff::Ptr{Ptr{Cvoid}})
    try
        clear_error(archive)
        bytes_read = reader_readbytes(archive, archive.data.data,
                                      archive.data.buff)
        unsafe_store!(buff, pointer(archive.data.buff))
        return Cssize_t(bytes_read)
    catch ex
        set_exception(archive, ex)
        return Cssize_t(0)
    end
end

function reader_skip_callback(c_archive, archive, request)
    try
        clear_error(archive)
        return Int64(reader_skip(archive, archive.data.data, request))
    catch ex
        set_exception(archive, ex)
        return Int64(0)
    end
end

function reader_seek_callback(c_archive, archive, request, whence)
    try
        clear_error(archive)
        return Int64(reader_seek(archive, archive.data.data, request, whence))
    catch ex
        set_exception(archive, ex)
        return Int64(0)
    end
end

function reader_close_callback(c_archive, archive)
    try
        clear_error(archive)
        reader_close(archive, archive.data.data)
        return Libc.errno(archive) == 0 ? Cint(0) : Status.WARN
    catch ex
        return set_exception(archive, ex)
    end
end

reader_readbytes(archive, io::IO, buff) = readbytes!(io, buff)
reader_skip(archive, io::IO, sz) =
    (skip(io, sz); sz)
function reader_seek(archive, io::IO, request, whence)
    if whence == Filesystem.SEEK_CUR
        skip(io, request)
        return position(io)
    elseif whence == Filesystem.SEEK_SET
        seek(io, request)
        return request
    elseif whence == Filesystem.SEEK_END
        seekend(io)
        skip(io, request)
        return position(io)
    else
        throw(ArgumentError("Unknown seek type: $whence"))
    end
end

function do_open(archive::Reader{T}) where {T<:GenericReadData}
    # Set various callbacks
    @_la_call(archive_read_set_callback_data,
              (Ptr{Cvoid}, Any), archive, archive)
    @_la_call(archive_read_set_open_callback,
              (Ptr{Cvoid}, Ptr{Cvoid}), archive,
              @to_open_callback(reader_open_callback, Reader{T}))
    @_la_call(archive_read_set_read_callback,
              (Ptr{Cvoid}, Ptr{Cvoid}), archive,
              @to_read_callback(reader_read_callback, Reader{T}))
    @_la_call(archive_read_set_seek_callback,
              (Ptr{Cvoid}, Ptr{Cvoid}), archive,
              @to_seek_callback(reader_seek_callback, Reader{T}))
    @_la_call(archive_read_set_skip_callback,
              (Ptr{Cvoid}, Ptr{Cvoid}), archive,
              @to_skip_callback(reader_skip_callback, Reader{T}))
    @_la_call(archive_read_set_close_callback,
              (Ptr{Cvoid}, Ptr{Cvoid}), archive,
              @to_close_callback(reader_close_callback, Reader{T}))
    @_la_call(archive_read_open1, (Ptr{Cvoid},), archive)
end

function next_header(archive::Reader)
    ensure_open(archive)
    entry = Entry(archive)
    @_la_call(archive_read_next_header2, (Ptr{Cvoid}, Ptr{Cvoid}),
              archive, entry)
    entry
end

"""
Retrieve the byte offset in UNCOMPRESSED data where last-read
header started.
"""
header_position(archive::Reader) =
    ccall((:archive_read_header_position, libarchive),
          Int64, (Ptr{Cvoid},), archive)

"""
Returns `Cint(1)` if the archive contains at least one encrypted entry.
If the archive format not support encryption at all
`ReadFormatEncryption.UNSUPPORTED` is returned.
If for any other reason (e.g. not enough data read so far)
we cannot say whether there are encrypted entries, then
`ReadFormatEncryption.DONT_KNOW` is returned.
In general, this function will return values below zero when the
reader is uncertain or totally incapable of encryption support.
When this function returns `0` you can be sure that the reader
supports encryption detection but no encrypted entries have
been found yet.

NOTE: If the metadata/header of an archive is also encrypted, you
cannot rely on the number of encrypted entries. That is why this
function does not return the number of encrypted entries but
just shows that there are some.
"""
has_encrypted_entries(archive::Reader) =
    ccall((:archive_read_has_encrypted_entries, libarchive), Cint, (Ptr{Cvoid},), archive)

"""
Returns a bitmask of capabilities that are supported by the archive format reader.
If the reader has no special capabilities, `ReadFormatCaps.NONE` is returned.
"""
format_capabilities(archive::Reader) =
    ccall((:archive_read_format_capabilities, libarchive), Cint, (Ptr{Cvoid},), archive)

"Read data from the body of an entry.  Similar to read(2)."
@inline function unsafe_archive_read(archive::Reader, ptr::Ptr{UInt8}, sz::UInt)
    nb = ccall((:archive_read_data, libarchive), Cssize_t,
               (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), archive, ptr, sz)
    nb < 0 && _la_error(Cint(nb), archive)
    nb % Csize_t
end

function Base.eof(archive::Reader)
    b = Ref{UInt8}()
    nb = ccall((:archive_read_data, libarchive), Cssize_t,
               (Ptr{Cvoid}, Ptr{UInt8}, Csize_t), archive, b, 1)
    if nb > 0
        skip(archive, -1)
        return false
    end
    return true
end

"Read data from the body of an entry.  Similar to read(2)."
function unsafe_read(archive::Reader, ptr::Ptr{UInt8}, sz::UInt)
    unsafe_archive_read(archive, ptr, sz) != sz && throw(EOFError())
    nothing
end

"Read data from the body of an entry.  Similar to read(2)."
function Base.read(archive::Reader, ::Type{UInt8})
    b = Ref{UInt8}()
    GC.@preserve b unsafe_read(archive, Base.unsafe_convert(Ptr{UInt8}, b), UInt(1))
    b[]
end

"Read data from the body of an entry.  Similar to read(2)."
function Base.readbytes!(archive::Reader, b::Array{UInt8}, nb=length(b))
    nbread = GC.@preserve b unsafe_archive_read(archive, Base.unsafe_convert(Ptr{UInt8}, b),
                                                UInt(nb))
    nbread < nb && resize!(b, nbread)
    nbread
end

"Read data from the body of an entry.  Similar to read(2)."
function Base.read(archive::Reader)
    # Maybe we can optimize this with block read?
    block_size = UInt(1024)
    res = Vector{UInt8}(undef, block_size)
    pos = 1
    GC.@preserve res while true
        ptr = Base.unsafe_convert(Ptr{UInt8}, Ref(res, pos))
        nbread = unsafe_archive_read(archive, ptr, block_size)
        if nbread < block_size
            resize!(res, nbread + pos - 1)
            return res
        end
        pos = Int(pos + block_size)
        resize!(res, pos + block_size - 1)
    end
    return res
end

"Read data from the body of an entry.  Similar to read(2)."
@inline Base.readavailable(archive::Reader) = read(archive)

"Seek within the body of an entry.  Similar to lseek(2)."
function Base.seek(archive::Reader, offset, what=Filesystem.SEEK_SET)
    res = ccall((:archive_seek_data, libarchive), Int64,
                (Ptr{Cvoid}, Int64, Cint), archive, offset, what)
    if res < 0
        _la_error(Cint(res), archive)
    end
    return archive
end
Base.skip(archive::Reader, request) = seek(archive, request, Filesystem.SEEK_CUR)
Base.seekend(archive::Reader, request=0) = seek(archive, request, Filesystem.SEEK_END)

"Skips entire entry"
Base.skip(archive::Reader) =
    @_la_call(archive_read_data_skip, (Ptr{Cvoid},), archive)

"Writes data to specified filedes"
read_into_fd(archive::Reader, fd::Integer) =
    @_la_call(archive_read_data_into_fd, (Ptr{Cvoid}, Cint), archive, fd)

Reader(f::Function) = archive_guard(f, Reader())
Reader(f::Function, args...; kws...) =
    archive_guard(f, Reader(args...; kws...))

add_passphrase(reader::Reader, passphrase::AbstractString) =
    @_la_call(archive_read_add_passphrase, (Ptr{Cvoid}, Cstring), reader, passphrase)
function set_passphrase(reader::Reader, cb::Cb) where Cb
    mgr = PassphraseMgr(cb)
    fptr = @cfunction(_passphrase_cb, Ptr{UInt8}, (Ptr{Cvoid},Ref{PassphraseMgr{Cb}}))
    GC.@preserve mgr begin
        @_la_call(archive_read_set_passphrase_callback,
                  (Ptr{Cvoid}, Ref{PassphraseMgr{Cb}}, Ptr{Cvoid}),
                  reader, mgr, fptr)
    end
    reader.passphrase_mgr = mgr
    return
end

set_format_option(reader::Reader, m::AbstractString, o::AbstractString, v::AbstractString) =
    @_la_call(archive_read_set_format_option,
              (Ptr{Cvoid}, Cstring, Cstring, Cstring), reader, m, o, v)
set_filter_option(reader::Reader, m::AbstractString, o::AbstractString, v::AbstractString) =
    @_la_call(archive_read_set_filter_option,
              (Ptr{Cvoid}, Cstring, Cstring, Cstring), reader, m, o, v)
set_option(reader::Reader, m::AbstractString, o::AbstractString, v::AbstractString) =
    @_la_call(archive_read_set_option,
              (Ptr{Cvoid}, Cstring, Cstring, Cstring), reader, m, o, v)
set_options(reader::Reader, opts::AbstractString) =
    @_la_call(archive_read_set_options, (Ptr{Cvoid}, Cstring), reader, opts)

# TODO
# /*
#  * A zero-copy version of archive_read_data that also exposes the file offset
#  * of each returned block.  Note that the client has no way to specify
#  * the desired size of the block.  The API does guarantee that offsets will
#  * be strictly increasing and that returned blocks will not overlap.
#  */
# int archive_read_data_block(struct archive *a,
#                    const void **buff, size_t *size, int64_t *offset);

# /*-
#  * Convenience function to recreate the current entry (whose header
#  * has just been read) on disk.
#  *
#  * This does quite a bit more than just copy data to disk. It also:
#  *  - Creates intermediate directories as required.
#  *  - Manages directory permissions:  non-writable directories will
#  *    be initially created with write permission enabled; when the
#  *    archive is closed, dir permissions are edited to the values specified
#  *    in the archive.
#  *  - Checks hardlinks:  hardlinks will not be extracted unless the
#  *    linked-to file was also extracted within the same session. (TODO)
#  */

# /* The "flags" argument selects optional behavior, 'OR' the flags you want. */

# int archive_read_extract(struct archive *, struct archive_entry *,
#                      int flags);
# int archive_read_extract2(struct archive *, struct archive_entry *,
#                      struct archive * /* dest */);
# void archive_read_extract_set_progress_callback(struct archive *,
#      void (*_progress_func)(void *), void *_user_data);

# /* Record the dev/ino of a file that will not be written.  This is
#  * generally set to the dev/ino of the archive being read. */
# void archive_read_extract_set_skip_file(struct archive *,
#      int64_t, int64_t);


text:
## generic `Tables.rows` and `Tables.columns` fallbacks
## if a table provides Tables.rows or Tables.columns,
## we'll provide a default implementation of the dual

# generic row iteration of columns
function rowcount(cols)
    props = propertynames(cols)
    isempty(props) && return 0
    return length(getproperty(cols, props[1]))
end

struct ColumnsRow{T}
    columns::T # a `Columns` object
    row::Int
end

Base.getproperty(c::ColumnsRow, ::Type{T}, col::Int, nm::Symbol) where {T} = getproperty(getfield(c, 1), T, col, nm)[getfield(c, 2)]
Base.getproperty(c::ColumnsRow, nm::Int) = getproperty(getfield(c, 1), nm)[getfield(c, 2)]
Base.getproperty(c::ColumnsRow, nm::Symbol) = getproperty(getfield(c, 1), nm)[getfield(c, 2)]
Base.propertynames(c::ColumnsRow) = propertynames(getfield(c, 1))

@generated function Base.isless(c::ColumnsRow{T}, d::ColumnsRow{T}) where {T <: NamedTuple{names}} where names
    exprs = Expr[]
    for n in names
        var1 = Expr(:., :c, QuoteNode(n))
        var2 = Expr(:., :d, QuoteNode(n))
        bl = quote
            a, b = $var1, $var2
            isless(a, b) && return true
            isequal(a, b) || return false
        end
        push!(exprs, bl)
    end
    push!(exprs, :(return false))
    Expr(:block, exprs...)
end

@generated function Base.isequal(c::ColumnsRow{T}, d::ColumnsRow{T}) where {T <: NamedTuple{names}} where names
    exprs = Expr[]
    for n in names
        var1 = Expr(:., :c, QuoteNode(n))
        var2 = Expr(:., :d, QuoteNode(n))
        push!(exprs, :(isequal($var1, $var2) || return false))
    end
    push!(exprs, :(return true))
    Expr(:block, exprs...)
end

struct RowIterator{T}
    columns::T
    len::Int
end
Base.eltype(x::RowIterator{T}) where {T} = ColumnsRow{T}
Base.length(x::RowIterator) = x.len
schema(x::RowIterator) = schema(x.columns)

function Base.iterate(rows::RowIterator, st=1)
    st > length(rows) && return nothing
    return ColumnsRow(rows.columns, st), st + 1
end

function rows(x::T) where {T}
    if columnaccess(T)
        cols = columns(x)
        return RowIterator(cols, Int(rowcount(cols)))
    elseif IteratorInterfaceExtensions.isiterable(x)
        return nondatavaluerows(x)
    end
    throw(ArgumentError("no default `Tables.rows` implementation for type: $T"))
end

# build columns from rows
"""
    Tables.allocatecolumn(::Type{T}, len) => returns a column type (usually AbstractVector) w/ size to hold `len` elements

    Custom column types can override with an appropriate "scalar" element type that should dispatch to their column allocator.
"""
allocatecolumn(T, len) = DataAPI.defaultarray(T, 1)(undef, len)

@inline function allocatecolumns(::Schema{names, types}, len) where {names, types}
    if @generated
        vals = Tuple(:(allocatecolumn($(fieldtype(types, i)), len)) for i = 1:fieldcount(types))
        return :(NamedTuple{Base.map(Symbol, names)}(($(vals...),)))
    else
        return NamedTuple{Base.map(Symbol, names)}(Tuple(allocatecolumn(fieldtype(types, i), len) for i = 1:fieldcount(types)))
    end
end

# add! will push! or setindex! a value depending on if the row-iterator HasLength or not
@inline add!(val, col::Int, nm, dest::AbstractArray, ::Union{Base.HasLength, Base.HasShape}, row) = setindex!(dest, val, row)
@inline add!(val, col::Int, nm, dest::AbstractArray, L, row) = push!(dest, val)

@inline function buildcolumns(schema, rowitr::T) where {T}
    L = Base.IteratorSize(T)
    len = Base.haslength(T) ? length(rowitr) : 0
    nt = allocatecolumns(schema, len)
    for (i, row) in enumerate(rowitr)
        eachcolumns(add!, schema, row, nt, L, i)
    end
    return nt
end

@inline add!(dest::AbstractArray, val, ::Union{Base.HasLength, Base.HasShape}, row) = setindex!(dest, val, row)
@inline add!(dest::AbstractArray, val, T, row) = push!(dest, val)

replacex(t, col::Int, x) = ntuple(i->i == col ? x : t[i], length(t))

@inline function add_or_widen!(val, col::Int, nm, dest::AbstractArray{T}, row, updated, L) where {T}
    if val isa T
        add!(dest, val, L, row)
        return
    else
        new = allocatecolumn(promote_type(T, typeof(val)), length(dest))
        row > 1 && copyto!(new, 1, dest, 1, row - 1)
        add!(new, val, L, row)
        updated[] = replacex(updated[], col, new)
        return
    end
end

function __buildcolumns(rowitr, st, sch, columns, rownbr, updated)
    while true
        state = iterate(rowitr, st)
        state === nothing && break
        row, st = state
        rownbr += 1
        eachcolumns(add_or_widen!, sch, row, columns, rownbr, updated, Base.IteratorSize(rowitr))
        columns !== updated[] && return __buildcolumns(rowitr, st, sch, updated[], rownbr, updated)
    end
    return updated
end

struct EmptyVector <: AbstractVector{Union{}}
    len::Int
end
Base.IndexStyle(::Type{EmptyVector}) = Base.IndexLinear()
Base.size(x::EmptyVector) = (x.len,)
Base.getindex(x::EmptyVector, i::Int) = throw(UndefRefError())

function _buildcolumns(rowitr, row, st, sch, columns, updated)
    eachcolumns(add_or_widen!, sch, row, columns, 1, updated, Base.IteratorSize(rowitr))
    return __buildcolumns(rowitr, st, sch, updated[], 1, updated)
end

# when Tables.schema(x) === nothing
@inline function buildcolumns(::Nothing, rowitr::T) where {T}
    state = iterate(rowitr)
    state === nothing && return NamedTuple()
    row, st = state
    names = Tuple(propertynames(row))
    len = Base.haslength(T) ? length(rowitr) : 0
    sch = Schema(names, nothing)
    columns = Tuple(EmptyVector(len) for _ = 1:length(names))
    return NamedTuple{Base.map(Symbol, names)}(_buildcolumns(rowitr, row, st, sch, columns, Ref{Any}(columns))[])
end

struct CopiedColumns{T}
    x::T
end

source(x::CopiedColumns) = getfield(x, :x)
istable(::Type{<:CopiedColumns}) = true
columnaccess(::Type{<:CopiedColumns}) = true
columns(x::CopiedColumns) = x
schema(x::CopiedColumns) = schema(source(x))
materializer(x::CopiedColumns) = materializer(source(x))
Base.propertynames(x::CopiedColumns) = propertynames(source(x))
Base.getproperty(x::CopiedColumns, nm::Symbol) = getproperty(source(x), nm)

@inline function columns(x::T) where {T}
    if rowaccess(T)
        r = rows(x)
        return CopiedColumns(buildcolumns(schema(r), r))
    elseif TableTraits.supports_get_columns_copy_using_missing(x)
        return CopiedColumns(TableTraits.get_columns_copy_using_missing(x))
    elseif IteratorInterfaceExtensions.isiterable(x)
        iw = nondatavaluerows(x)
        return CopiedColumns(buildcolumns(schema(iw), iw))
    end
    throw(ArgumentError("no default `Tables.columns` implementation for type: $T"))
end


text:
using PyCall

"""
Converts an array to a Python array.  Given `revdims=true`, convert to row major ordering.
e.g. dims of (1, 2, 3, 4) on Julia side -> (4, 3, 2, 1) on Python side

.  Does not modify underlying memory.

Temporary utility function until next release of PyCall (as of 1.4.0). Adds `revdims` flag.
"""
function to_python_array{T<:PyCall.NPY_TYPES}(a::StridedArray{T}, revdims::Bool=false)
    @PyCall.npyinitialize
    size_a = revdims ? reverse(size(a)) : size(a)
    strides_a = revdims ? reverse(strides(a)) : strides(a)
    p = ccall(PyCall.npy_api[:PyArray_New], PyPtr,
              (PyPtr,Cint,Ptr{Int},Cint, Ptr{Int},Ptr{T}, Cint,Cint,PyPtr),
              PyCall.npy_api[:PyArray_Type],
              ndims(a), Int[size_a...], PyCall.npy_type(T),
              Int[strides_a...] * sizeof(eltype(a)), a, sizeof(eltype(a)),
              PyCall.NPY_ARRAY_ALIGNED | PyCall.NPY_ARRAY_WRITEABLE,
              C_NULL)
    return PyObject(p, a)
end

"""
Convert array to row major, C style arrays.
Does not modify memory during conversion
"""
function to_python_array(arr::AbstractArray)
    to_python_array(arr, size(arr)...)
end

"""
Convert array to row major, C style arrays.
Does not modify buffer during conversion

Specifically for bit arrays -> UInt8

Specify dimensions IN JULIA COORDINATES
Given Julia array with dims: (1, 2, 3)
Return Python array with dims: (3, 2, 1)
"""
function to_python_array(arr::AbstractArray, dims...)
    reshaped = reshape(arr, dims)
    rounded = round(UInt8, reshaped)
    to_python(rounded, true)
end