text:
# TODO: add waveguide sites for matched conditions

"""
	module Simulations

Package data used to define simulation and generate the relevant differential operators.
Exports are `Simulation`, `smooth!`, `update_dielectric!`, `update_pump`
"""
module Simulations

const N_NEAREST_SITES_FOR_TESSELLATION = 15

export Simulation,
smooth!,
update_dielectric!,
update_pump!,
update_sim!

using ...Defaults
using ..Boundaries
using ..Domains
using ..Lattices
using ..Shapes
using ..Tessellations
using Combinatorics
using LinearAlgebra
using ProgressMeter
using RecipesBase
using SparseArrays
using Statistics
using VoronoiDelaunay

"""
	struct Simulation

	Simulation(domains...; verbose=true) -> sim

	(::Simulation)(;kwargs...) -> new_sim

Takes all the given `domains` and builds a simulation object which can be passed
to the various solvers to find resonance frequencies, etc. The where domains overlap,
the one that appears first in the arguments of `Simulation` takes priority.

Contains all the information that defines a simulation (except for frequency).
It contains several kinds of fields:
 * vectors of data about each site, such as `x`, `y`, `ε`, `F`
 * boolean vectors that give geometric properites of each site, such as whether it is an `interior` point
 * vectors of weights and indices giving the coupling between sites
 * vectors of weights and indices giving the self-couplings
 * the `tessellation` used to guide interpolation

 Constructing a `Simulation` is sped up if multiple threads are running,
 i.e. by declaring the environment variable "JULIA\\_NUM\\_THREADS" *before*
 starting Julia.
"""
struct Simulation
	domains::Array{Domain,1}
	nsites::Int
	n::Int

	# site boolean properties, one for each site
	interior::BitArray{1}
	exterior::BitArray{1}
	corner::BitArray{1}
	bulk::BitArray{1}
	surface::BitArray{1}
	lattice_wall::BitArray{1}
	domain_wall::BitArray{1}

	# site properties, indices
	domain_index::Array{Int,1}
	nnxm::Array{Int,1}
	nnxp::Array{Int,1}
	nnym::Array{Int,1}
	nnyp::Array{Int,1}
	ij::Array{CartesianIndex{2},1}
	exterior_from::Array{Int,1}
	exterior_side::Array{Int,1}
	not_this_domains_exterior::Array{Array{Int,1},1}
	exterior_inds::Array{Int,1}

	# site properties, other
	x::Array{Float64,1}
	y::Array{Float64,1}
	r::Array{Float64,1}
	ε::Array{ComplexF64,1}
    F::Array{Float64,1}
	name::Array{Symbol,1}
	distance::Array{Array{Float64,1},1}
	normal::Array{Array{Array{Float64,1},1},1}
	tangent::Array{Array{Array{Float64,1},1},1}
	side::Array{Array{Int,1},1}

	# tessellation::Tessellation
	smoothed::Base.RefValue{Bool}

	# links between sites
	link_start::Array{Int,1}
	link_stop::Array{Int,1}
	link_weight::Array{Float64,1}
	link_half_x::Array{Float64,1}
	link_half_y::Array{Float64,1}
	link_half_r::Array{Float64,1}
	link_x_bool::BitArray{1}
	link_y_bool::BitArray{1}
	link_r_bool::BitArray{1}

	# self-energy terms
	self_index::Array{Int,1}
	self_weight::Array{Float64,1}
	self_half_x::Array{Float64,1}
	self_half_y::Array{Float64,1}
	self_half_r::Array{Float64,1}
	self_x_bool::BitArray{1}
	self_y_bool::BitArray{1}
	self_r_bool::BitArray{1}

	Simulation(sim::Simulation) = Simulation(sim.domains;kwargs...)
    Simulation(domains...; kwargs...) = Simulation(domains; kwargs...)
    function Simulation(domains::Tuple; verbose::Bool=false)

		domain_index = Int[]
		nnxm, nnxp, nnym, nnyp = Int[], Int[], Int[], Int[]
		name = Symbol[]
		x, y, r, θ = Float64[], Float64[], Float64[], Float64[]
		ij = CartesianIndex{2}[]
		ε, F = ComplexF64[], Float64[]
		interior, surface, bulk = BitArray([]), BitArray([]), BitArray([])
		corner = BitArray([])

		verbose ? println("\nBuilding Simulation...") : nothing
		verbose ? println("\tCombining $(length(domains)) domains") : nothing
			append_domains!(domain_index,nnxm,nnxp,nnym,nnyp,name,x,y,r,θ,ij,ε,F,interior,bulk,surface,corner,domains)

			removed = get_sites_to_remove(domains,domain_index,x,y)

			trim_nn!(nnxm,nnxp,nnym,nnyp,removed,bulk)

			remove_sites!((domain_index,name,x,y,r,θ,ij,ε,F,interior,bulk,surface,corner,nnxm,nnxp,nnym,nnyp),removed)

		verbose ? println("\tComputing $(sum(surface)) normals and extending domain") : nothing
			lattice_wall, domain_wall = generate_lattice_and_domain_wall!(interior,bulk,surface,corner,domains,domain_index,x,y,ij)

			exterior_x, exterior_y, exterior_from, normal, tangent, distance, side, exterior_side = generate_exterior!(corner,domains,domain_index,x,y,surface,verbose)

			exterior = falses(length(x))

			append_exterior!(domain_index,nnxm,nnxp,nnym,nnyp,name,x,y,r,θ,ij,ε,F,interior,bulk,surface,corner,
				exterior,domain_wall,lattice_wall,normal,tangent,distance,side,exterior_from,exterior_side,exterior_x,exterior_y)

			not_this_domains_exterior = Array{Array{Int,1}}(undef,length(domains))#falses(length(x),length(domains))
			for d ∈ eachindex(domains)
				not_this_domains_exterior[d] = sum(interior) .+ findall(!isequal(d),domain_index[exterior_from[exterior]])
			end
			exterior_inds = findall(exterior)

		# verbose ? println("\tTessellating $(length(x)) points") : nothing
			# tessellation = Tessellation(x,y)

		bulk_link_size = 4sum(bulk)
		lattice_wall_link_size = 4NNN_BULK*sum(lattice_wall)
		surface_link_size = (NNN_BULK+2NNN_SURFACE)*sum(surface)
		surface_self_size = sum(surface)
		corner_link_size = 2NNN_CORNER*sum(corner)
		corner_self_size = 2sum(corner)

		link_start = Array{Int}(undef,bulk_link_size+lattice_wall_link_size+surface_link_size+corner_link_size)
		link_stop = Array{Int}(undef,length(link_start))
		link_weight = Array{Float64}(undef,length(link_start))
		link_half_x = Array{Float64}(undef,length(link_start))
		link_half_y = Array{Float64}(undef,length(link_start))
		link_half_r = Array{Float64}(undef,length(link_start))
		link_x_bool = falses(length(link_start))
		link_y_bool = falses(length(link_start))
		link_r_bool = falses(length(link_start))

		self_index = Array{Int}(undef,surface_self_size+corner_self_size)
		self_weight = Array{Float64}(undef,length(self_index))
		self_half_x = Array{Float64}(undef,length(self_index))
		self_half_y = Array{Float64}(undef,length(self_index))
		self_half_r = Array{Float64}(undef,length(self_index))
		self_x_bool = falses(length(self_index))
		self_y_bool = falses(length(self_index))
		self_r_bool = falses(length(self_index))

		link_fiduciary = falses(length(link_start))
		self_fiduciary = falses(length(self_index))

		sim = new([domains...],length(x),sum(interior),interior,exterior,corner,bulk,
			surface,lattice_wall,domain_wall,
			domain_index,nnxm,nnxp,nnym,nnyp,ij,
			exterior_from,exterior_side,not_this_domains_exterior,exterior_inds,
			x,y,r,ε,F,name,distance,normal,tangent,side,
			# tessellation,
			Ref(false),
			link_start,link_stop,link_weight,
			link_half_x,link_half_y,link_half_r,link_x_bool,link_y_bool,link_r_bool,
			self_index,self_weight,
			self_half_x,self_half_y,self_half_r,self_x_bool,self_y_bool,self_r_bool)

		verbose ? println("\tLinking $(sum(sim.bulk)) bulk points") : nothing
			link_bulk!(link_fiduciary,sim,verbose)

		verbose ? println("\tLinking $(sum(sim.lattice_wall)) lattice mismatch points") : nothing
			link_lattice_wall!(link_fiduciary,sim,verbose)

		verbose ? println("\tLinking $(sum(sim.surface)) surface points") : nothing
			link_surface!(link_fiduciary,self_fiduciary,sim,verbose)

		verbose ? println("\tLinking $(sum(sim.corner)) corner points") : nothing
			link_corner!(link_fiduciary,self_fiduciary,sim,verbose)

		verbose ? println("\tCleaning up $(sum(link_fiduciary)+sum(self_fiduciary)) links") : nothing
			trim_link!(sim,link_fiduciary)

			trim_self!(sim,self_fiduciary)

		verbose ? println("\tApplying boundary conditions") : nothing
			apply_local_bc!(sim)

			apply_nonlocal_bc!(sim)

		verbose ? println("\tDone!") : nothing
		return sim
    end

	function Base.show(io::IO,sim::Simulation)
		if sim.smoothed[]
	        print(io,"Smoothed Simulation with ", length(sim.domains) ," domain")
		else
			print(io,"Unsmoothed Simulation with ", length(sim.domains)," domain")
		end
		length(sim.domains)>1 ? println(io,"s") : println(io)
        println(io,"\tsites: ",sim.nsites)
        println(io,"\tinterior: ", sim.n)
		println(io,"\texterior: ",sum(sim.exterior))
		println(io,"\tbulk: ",sum(sim.bulk))
        println(io,"\tsurface: ",sum(sim.surface))
		println(io,"\tcorner: ",sum(sim.corner))
		println(io,"\tlattice wall: ",sum(sim.lattice_wall))
		println(io,"\tdomain wall: ",sum(sim.domain_wall))
		println(io,"\t\t====================")
		for d ∈ eachindex(sim.domains)
			println(io,"\t\tDomain ", d,": ",sim.domains[d].type)
		end
	end
end

"""
	smooth!(sim,[num_sub_pixel])

Use sub-pixel sampling along domain walls to smooth the dieletric and pump profiles.

`num_sub_pixel` defaults to `NUM_SUB_PIXEL` in module `Defaults`, which is currently set to 5.
"""
function smooth!(sim::Simulation,num_sub_pixel::Int=NUM_SUB_PIXEL)
    domains = sim.domains
    domain_wall = sim.domain_wall
	lattices = map(z->z.lattice,sim.domains)
    X = Array{Float64}(undef,num_sub_pixel,num_sub_pixel)
    Y = Array{Float64}(undef,num_sub_pixel,num_sub_pixel)
    E = Array{ComplexF64}(undef,num_sub_pixel,num_sub_pixel)
    F = Array{Float64}(undef,num_sub_pixel,num_sub_pixel)
	dw_inds = findall(domain_wall)
    for dw ∈ eachindex(dw_inds)
		d = dw_inds[dw]
		ds = sim.domain_index[d]
	    I,J = Tuple(sim.ij[d])
        is = LinRange(I-.5,I+.5,num_sub_pixel)
        js = LinRange(J-.5,J+.5,num_sub_pixel)
        for i ∈ eachindex(is), j ∈ eachindex(js)
            X[i,j],Y[i,j] = lattices[ds](is[i],js[j])
        end
        simulation_dielectric!(E,domains,X,Y)
        simulation_pump!(F,domains,X,Y)
        sim.ε[d] = mean(E)
        sim.F[d] = mean(F)
    end
	sim.smoothed[] = true
    return nothing
end

# take individually defined domains and concatenate all relevant data
function append_domains!(domain_index,nnxm,nnxp,nnym,nnyp,name,x,y,r,θ,ij,ε,F,interior,bulk,surface,corner,domains)
	count = 1
	for d ∈ domains
		append!(domain_index,fill(count,length(d.x)))
		append!(nnxm,d.nnxm.+length(nnxm))
		append!(nnxp,d.nnxp.+length(nnxp))
		append!(nnym,d.nnym.+length(nnym))
		append!(nnyp,d.nnyp.+length(nnyp))
		append!(name,fill(d.name,length(d.x)))
		append!(x,d.x)
		append!(y,d.y)
		append!(r,d.r)
		append!(θ,d.θ)
		append!(ij,d.ij)
		append!(ε,d.ε)
		append!(F,d.F)
		append!(interior,d.interior)
		append!(bulk,d.bulk)
		append!(surface,d.surface)
		append!(corner,d.corner)
		count += 1
	end
	return nothing
end

# identify which sites are to be removed by determining which sites belong to the topmost domain.
get_sites_to_remove(domains,dom,x,y) = which_domains(domains,x,y).!==dom

# identify which domain is associated with a given spatial position (x,y)
function which_domains(domains,x,y)
    doms = Array{Int}(undef,length(x))
	shapes = map(z->z.boundary.shape,domains)
    which_domains!(doms,domains,shapes,x,y)
    return doms
end
function which_domains!(doms::Array{Int},domains,shapes,x,y)
    for i ∈ eachindex(x)
		doms[i] = which_domain(domains,shapes,x[i],y[i])
    end
    return nothing
end
function which_domain(domains,x,y)
	bnd = map(z->z.boundary,domains)
	shapes = map(z->z.boundary.shape,domains)
    return which_domain(domains,shapes,x,y)
end
@inline function which_domain(domains,shapes,x,y)
	dom = findfirst(map(z->z(x,y),shapes))
    isnothing(dom) ? dom=0 : nothing
    iszero(dom) ? nothing : (isvoid(domains[dom]) ? dom=0 : nothing)
    return dom
end
isvoid(domain::Domain) = domain.type==:Void

# remove the pre-computed nearest neighbors, which makes the bulk nn calculation very fast
function trim_nn!(nnxm,nnxp,nnym,nnyp,removed,bulk)
	cr = cumsum(removed)
	for i ∈ eachindex(nnxm)
		if bulk[i]
			nnxm[i] -= cr[nnxm[i]]
			nnxp[i] -= cr[nnxp[i]]
			nnym[i] -= cr[nnym[i]]
			nnyp[i] -= cr[nnyp[i]]
		end
	end
end

# remove the sites determined by get_sites_to_remove
remove_sites!(sites,removed) = map(z->deleteat!(z,findall(removed)),sites)

# compute which sites form the walls along a lattice mismatch and domain walls
function generate_lattice_and_domain_wall!(interior,bulk,surface,corner,domains,domain_index,x,y,ij)
	lattices = map(z->z.lattice,domains)
	shapes = map(z->z.boundary.shape,domains)
	# DX, DY = map(z->z.dx,lattices), map(z->z.dy,lattices)
	domain_wall = falses(size(x)...)
	lattice_wall = falses(size(x)...)
	for i ∈ eachindex(x)
		if interior[i]
			ds = domain_index[i]
			xt,yt = lattices[ds](ij[i]-CartesianIndex(1,0))
			dom1 = which_domain(domains,shapes,xt,yt)
			xt,yt = lattices[ds](ij[i]+CartesianIndex(1,0))
			dom2 = which_domain(domains,shapes,xt,yt)
			xt,yt = lattices[ds](ij[i]-CartesianIndex(0,1))
			dom3 = which_domain(domains,shapes,xt,yt)
			xt,yt = lattices[ds](ij[i]+CartesianIndex(0,1))
			dom4 = which_domain(domains,shapes,xt,yt)

			domain_wall[i] = !(domain_index[i]==dom1==dom2==dom3==dom4)
			domain_wall[i] = domain_wall[i] && (dom1!==0 && dom2!==0 && dom3!==0 && dom4!==0)
			if domain_wall[i]
				lat1 = lattices[dom1]
				lat2 = lattices[dom2]
				lat3 = lattices[dom3]
				lat4 = lattices[dom4]
				lattice_wall[i] = !all(isequal(lattices[domain_index[i]]),(lat1,lat2,lat3,lat4))
				corner[i] = false
				surface[i] = false
				bulk[i] = true
			end
			if lattice_wall[i]
				corner[i] = false
				surface[i] = false
				bulk[i] = false
			end
		end
	end
	return lattice_wall, domain_wall
end

# do a mirror reflection of the surface about the boundary, used in stabilizing the computation of nearest neighbors via a tessellation
function generate_exterior!(corner,domains,domain_index,x,y,surface,verbose::Bool)
	bnds = map(z->z.boundary,domains) ; bcs = map(z->z.bcs,bnds)
	shapes = map(z->z.shape,bnds)
	dxs, dys = map(z->z.lattice.dx,domains), map(z->z.lattice.dy,domains)

	ex_x = NaN*Array{Float64}(undef,sum(surface)+3sum(corner))
	ex_y = NaN*Array{Float64}(undef,sum(surface)+3sum(corner))
	ex_from = zeros(Int,sum(surface)+3sum(corner))
	ex_side = zeros(Int,sum(surface)+3sum(corner))
	normal = Array{Array{Array{Float64,1},1}}(undef,length(x))
	tangent = Array{Array{Array{Float64,1},1}}(undef,length(x))
	distance = Array{Array{Float64,1}}(undef,length(x))
	side = Array{Array{Int,1}}(undef,length(x))

	pg = Progress(sum(surface),.1,"\t\t")
	pg_ref = Threads.Atomic{Int}(0)
	surface_inds = findall(surface .& .!corner)
	for si ∈ eachindex(surface_inds)
		i = surface_inds[si]
		ds = domain_index[i]
        normal[i],tangent[i],distance[i],side[i] = normal_distance(shapes[ds],x[i],y[i])
		ex_x[si] = x[i] - normal[i][1][1]*2distance[i][1]
		ex_y[si] = y[i] - normal[i][1][2]*2distance[i][1]
		ex_from[si] = i
		ex_side[si] = side[i][1]
		verbose ? Threads.atomic_add!(pg_ref,1) : nothing
		if verbose
			update!(pg,pg_ref[])
		end
	end
	corner_inds = findall(corner)
	count = length(surface_inds)+1
	for ci ∈ eachindex(corner_inds)
		i = corner_inds[ci]
		ds = domain_index[i]
        normal[i],tangent[i],distance[i],side[i] = normal_distance(shapes[ds],x[i],y[i])
		dx = min(dxs[ds],dys[ds])
        ex_x[count] = x[i] - normal[i][1][1]*2distance[i][1]
		ex_y[count] = y[i] - normal[i][1][2]*2distance[i][1]
		ex_from[count] = i
		ex_side[count] = side[i][1]
		count += 1
		if corner[i]
			ex_x[count] = x[i] - normal[i][2][1]*2distance[i][2]
			ex_y[count] = y[i] - normal[i][2][2]*2distance[i][2]
			ex_from[count] = i
			ex_side[count] = side[i][2]
			count += 1

			crn_distance = Inf
			crn_index = 0
			for j ∈ eachindex(shapes[ds].corners)
				crn = shapes[ds].corners[j]
				crn_distance_temp = hypot(x[i]-crn[1],y[i]-crn[2])
				if crn_distance_temp < crn_distance
					crn_index=j
					crn_distance = crn_distance_temp
				end
			end
			crn_normal = [x[i]-shapes[ds].corners[crn_index][1],y[i]-shapes[ds].corners[crn_index][2]]
			if norm(crn_normal) > eps()
				crn_normal = crn_normal/norm(crn_normal)
				ex_x[count] = x[i] - crn_normal[1]*2crn_distance
				ex_y[count] = y[i] - crn_normal[2]*2crn_distance
				ex_from[count] = i
				ex_side[count] = length(shapes[ds].corners)+crn_index
				count += 1
			end
		end
	end
	verbose ? finish!(pg) : nothing
	return ex_x[1:count-1], ex_y[1:count-1], ex_from[1:count-1], normal, tangent, distance, side, ex_side[1:count-1]
end

# add the exterior points to the data previously concatenated from the given domains
function append_exterior!(domain_index,nnxm,nnxp,nnym,nnyp,name,x,y,r,θ,ij,ε,F,interior,bulk,surface,corner,
		exterior,domain_wall,lattice_wall,normal,tangent,distance,side,exterior_from,
		exterior_side,exterior_x,exterior_y)

	append!(domain_index,zeros(Int,length(exterior_x)))
	append!(nnxm,zeros(Int,length(exterior_x)))
	append!(nnxp,zeros(Int,length(exterior_x)))
	append!(nnym,zeros(Int,length(exterior_x)))
	append!(nnyp,zeros(Int,length(exterior_x)))
	append!(name,Array{Symbol}(undef,length(exterior_x)))
	append!(ij,Array{CartesianIndex{2}}(undef,length(exterior_x)))
	append!(ε,NaN*ones(ComplexF64,length(exterior_x)))
	append!(F,zeros(Float64,length(exterior_x)))

	append!(interior,falses(length(exterior_x)))
	append!(bulk,falses(length(exterior_x)))
	append!(surface,falses(length(exterior_x)))
	append!(corner,falses(length(exterior_x)))
	append!(exterior,trues(length(exterior_x)))
	append!(domain_wall,falses(length(exterior_x)))
	append!(lattice_wall,falses(length(exterior_x)))

	append!(normal,Array{Array{Array{Float64,1},1}}(undef,length(exterior_x)))
	append!(tangent,Array{Array{Array{Float64,1},1}}(undef,length(exterior_x)))
	append!(distance,Array{Array{Float64,1}}(undef,length(exterior_x)))
	append!(side,Array{Array{Int,1}}(undef,length(exterior_x)))

	prepend!(exterior_from,zeros(Int,length(x)))
	prepend!(exterior_side,zeros(Int,length(x)))
	append!(x,exterior_x)
	append!(y,exterior_y)
	append!(r,NaN*ones(Float64,length(exterior_x)))
	return nothing
end

# KEY BLOCKS: LINKS BULK POINTS via precomputed nn and lookup tables
function link_bulk!(link_fiduciary::BitArray,sim::Simulation,verbose::Bool)

	domains,x,y,ij = sim.domains, sim.x, sim.y, sim.ij
	lattices = map(z->z.lattice,domains)

	DX, DY = map(z->z.dx,lattices), map(z->z.dy,lattices)
	V1, V2 = map(z->z.v1,lattices), map(z->z.v2,lattices)
	DR, DΘ = map(z->z.dr,lattices), map(z->z.dθ,lattices)

	ijlat = map((x,y)->(x,lattices[y]),ij[sim.interior],sim.domain_index[sim.interior])

	pg = Progress(sum(sim.bulk),.1,"\t\t")
	pg_ref = Threads.Atomic{Int}(0)
	c0 = 0
	bulk_inds = findall(sim.bulk)
	# Threads.@threads
	for bi ∈ eachindex(bulk_inds)
		c = c0 + 4(bi-1)+1
		i = bulk_inds[bi]
        begin
			ds = sim.domain_index[i]
			lds = lattices[ds]
			dx, dy = DX[ds], DY[ds]
			dr, dθ = DR[ds], DΘ[ds]

			sim.link_start[c] = i
			sim.link_start[c+1] = i
			sim.link_start[c+2] = i
			sim.link_start[c+3] = i

			if !sim.domain_wall[i]
				sim.link_stop[c  ] = sim.nnxm[i]
	            sim.link_stop[c+1] = sim.nnxp[i]
	            sim.link_stop[c+2] = sim.nnym[i]
	            sim.link_stop[c+3] = sim.nnyp[i]
			else
				sim.link_stop[c  ] = findfirst(isequal((ij[i]+CartesianIndex(-1, 0),lds)),ijlat)
				sim.link_stop[c+1] = findfirst(isequal((ij[i]+CartesianIndex( 1, 0),lds)),ijlat)
				sim.link_stop[c+2] = findfirst(isequal((ij[i]+CartesianIndex( 0,-1),lds)),ijlat)
				sim.link_stop[c+3] = findfirst(isequal((ij[i]+CartesianIndex( 0, 1),lds)),ijlat)
			end

			if lds.type==:Cartesian
				sim.link_weight[c  ] = 1/dx^2
				sim.link_weight[c+1] = 1/dx^2
				sim.link_weight[c+2] = 1/dy^2
				sim.link_weight[c+3] = 1/dy^2
			elseif lds.type==:Polar
				sim.link_weight[c  ] = 1/dr^2
				sim.link_weight[c+1] = 1/dr^2
				sim.link_weight[c+2] = 1/dθ^2
				sim.link_weight[c+3] = 1/dθ^2

				sim.link_half_r[c] = lds.r0 + (sim.ij[i][1]-1/2)*lds.dr
				sim.link_half_r[c+1] = lds.r0 + (sim.ij[i][1]+1/2)*lds.dr
				sim.link_r_bool[c] = true
			else
				throw(LatticeError(lds.type))
			end

			sim.link_half_x[c],sim.link_half_y[c] = lds(Tuple(sim.ij[i]).-(1/2,0))
			sim.link_x_bool[c] = true

			sim.link_half_x[c+1],sim.link_half_y[c+1] = lds(Tuple(sim.ij[i]).+(1/2,0))

			sim.link_x_bool[c+1] = true

			sim.link_half_x[c+2],sim.link_half_y[c+2] = lds(Tuple(sim.ij[i]).-(0,1/2))
			sim.link_y_bool[c+2] = true

			sim.link_half_x[c+3],sim.link_half_y[c+3] = lds(Tuple(sim.ij[i]).+(0,1/2))
			sim.link_y_bool[c+3] = true

			link_fiduciary[c] = true
			link_fiduciary[c+1] = true
			link_fiduciary[c+2] = true
			link_fiduciary[c+3] = true

		end
		verbose ? Threads.atomic_add!(pg_ref,1) : nothing
		if verbose && Threads.threadid()==Threads.nthreads()
			update!(pg,pg_ref[])
		end
    end
	verbose ? finish!(pg) : nothing
    return nothing
end

# KEY BLOCKS: LINKS LATTICE WALL POINTS via interpolation
function link_lattice_wall!(link_fiduciary::BitArray,sim::Simulation,verbose::Bool)

	NNN = NNN_BULK

	domains,x,y,ij = sim.domains,sim.x,sim.y,sim.ij

	lattices = map(z->z.lattice,domains)
	DX, DY = map(z->z.dx,lattices), map(z->z.dy,lattices)
	DR, DΘ = map(z->z.dr,lattices), map(z->z.dθ,lattices)

	pg = Progress(sum(sim.lattice_wall),.1,"\t\t")
	pg_ref = Threads.Atomic{Int}(0)
	c0 = findlast(link_fiduciary)
	isnothing(c0) ? c0=0 : nothing

	lattice_wall_inds = findall(sim.lattice_wall)
	# Threads.@threads
	for lwi ∈ eachindex(lattice_wall_inds)
		c = c0 + 4NNN*(lwi-1)+1
		i = lattice_wall_inds[lwi]
		begin
			ds = sim.domain_index[i]
			lds = lattices[ds]

			if lds.type==:Cartesian
				dx, dy = DX[ds], DY[ds]
			elseif lds.type==:Polar
				dr, dθ = DR[ds], DΘ[ds]
				dx = dr
				dy = sim.r[i]*dθ
			else
				throw(LatticeError(lds.type))
			end

			x1, y1 = lds(Tuple(sim.ij[i]).-(1,0))
			x2, y2 = lds(Tuple(sim.ij[i]).+(1,0))
			x3, y3 = lds(Tuple(sim.ij[i]).-(0,1))
			x4, y4 = lds(Tuple(sim.ij[i]).+(0,1))

			hx1, hy1 = lds(Tuple(sim.ij[i]).-(1/2,0))
			hx2, hy2 = lds(Tuple(sim.ij[i]).+(1/2,0))
			hx3, hy3 = lds(Tuple(sim.ij[i]).-(0,1/2))
			hx4, hy4 = lds(Tuple(sim.ij[i]).+(0,1/2))

			it,wt = generate_inds_weights(x1,y1,sim,NNN,i,ds)
			N = length(it)
			for j ∈ 1:N
				sim.link_start[c+j-1] = i
				sim.link_stop[c+j-1] = it[j]
				sim.link_weight[c+j-1] = wt[j]/dx^2
				sim.link_half_x[c+j-1] = hx1
				sim.link_half_y[c+j-1] = hy1
				sim.link_x_bool[c+j-1] = true
				link_fiduciary[c+j-1] = true
			end
			c += NNN

			it,wt = generate_inds_weights(x2,y2,sim,NNN,i,ds)
			N = length(it)
			for j ∈ 1:N
				sim.link_start[c+j-1] = i
				sim.link_stop[c+j-1] = it[j]
				sim.link_weight[c+j-1] = wt[j]/dx^2
				sim.link_half_x[c+j-1] = hx2
				sim.link_half_y[c+j-1] = hy2
				sim.link_x_bool[c+j-1] = true
				link_fiduciary[c+j-1] = true
			end
			c += NNN

			it,wt = generate_inds_weights(x3,y3,sim,NNN,i,ds)
			N = length(it)
			for j ∈ 1:N
				sim.link_start[c+j-1] = i
				sim.link_stop[c+j-1] = it[j]
				sim.link_weight[c+j-1] = wt[j]/dy^2
				sim.link_half_x[c+j-1] = hx3
				sim.link_half_y[c+j-1] = hy3
				sim.link_y_bool[c+j-1] = true
				link_fiduciary[c+j-1] = true
			end
			c += NNN

			it,wt = generate_inds_weights(x4,y4,sim,NNN,i,ds)
			N = length(it)
			for j ∈ 1:N
				sim.link_start[c+j-1] = i
				sim.link_stop[c+j-1] = it[j]
				sim.link_weight[c+j-1] = wt[j]/dy^2
				sim.link_half_x[c+j-1] = hx4
				sim.link_half_y[c+j-1] = hy4
				sim.link_y_bool[c+j-1] = true
				link_fiduciary[c+j-1] = true
			end
		end
		verbose ? Threads.atomic_add!(pg_ref,1) : nothing
		if verbose && Threads.threadid()==Threads.nthreads()
			update!(pg,pg_ref[])
		end
	end
	verbose ? finish!(pg) : nothing
	return nothing
end

# KEY BLOCKS: LINKS SURFACE POINTS via interpolation
function link_surface!(link_fiduciary::BitArray,self_fiduciary::BitArray,sim::Simulation,verbose::Bool)

	domains,x,y = sim.domains,sim.x,sim.y
	lattices = map(z->z.lattice,domains)
	DX, DY = map(z->z.dx,lattices), map(z->z.dy,lattices)
	DR, DΘ = map(z->z.dr,lattices), map(z->z.dθ,lattices)

	surface_inds = findall(sim.surface .& .!sim.corner)
	pg = Progress(length(surface_inds),.1,"\t\t")
	pg_ref = Threads.Atomic{Int}(0)
	cl0 = findlast(link_fiduciary)
	isnothing(cl0) ? cl0=0 : nothing
	cs0 = 0
	# Threads.@threads
	for si ∈ eachindex(surface_inds)

		i = surface_inds[si]
		cs = cs0 + si
		cl = cl0 + (NNN_BULK+2NNN_SURFACE)*(si-1) + 1

		begin
			ds = sim.domain_index[i]
			n,t,d,sides = sim.normal[i],sim.tangent[i],sim.distance[i],sim.side[i]
			lds = lattices[ds]

			dxn = 2d[1]
			NNN = NNN_BULK

			hxn1, hyn1 = x[i] - n[1][1]*dxn/2, y[i] - n[1][2]*dxn/2
			sim.self_index[cs] = i
			sim.self_weight[cs] = -1/dxn^2
			sim.self_half_x[cs] = hxn1
			sim.self_half_y[cs] = hyn1
			if sides[1] ∈ (1,2)
				sim.self_x_bool[cs] = true
			else
				sim.self_y_bool[cs] = true
			end
			self_fiduciary[cs] = true

			if lds.type==:Cartesian
				dxn = max(dxn,min(DX[ds],DY[ds]))
			elseif lds.type==:Polar
				dxn = max(dxn,min(DR[ds],sim.r[i]*DΘ[ds]))
			else
				throw(LatticeError(lds.type))
			end

			xn2, yn2   = x[i] + n[1][1]*dxn  , y[i] + n[1][2]*dxn
			hxn2, hyn2 = x[i] + n[1][1]*dxn/2, y[i] + n[1][2]*dxn/2
			it,wt = generate_inds_weights(xn2,yn2,sim,NNN,i,ds)
			N = length(it)
			for j ∈ 1:N
				sim.link_start[cl+j-1] = i
				sim.link_stop[cl+j-1] = it[j]
				sim.link_weight[cl+j-1] = wt[j]/dxn^2
				sim.link_half_x[cl+j-1] = hxn2
				sim.link_half_y[cl+j-1] = hyn2
				if sides[1] ∈ (1,2)
					sim.link_x_bool[cl+j-1] = true
				else
					sim.link_y_bool[cl+j-1] = true
				end
				link_fiduciary[cl+j-1] = true
			end
			cl += NNN

			if lds.type==:Cartesian
				dxt = min(DX[ds],DY[ds])
			elseif lds.type==:Polar
				dxt = min(DR[ds],sim.r[i]*DΘ[ds])
			else
				throw(LatticeError(lds.type))
			end

			NNN =  NNN_SURFACE

			xt1, yt1   = x[i] - t[1][1]*dxt  , y[i] - t[1][2]*dxt
			hxt1, hyt1 = x[i] - t[1][1]*dxt/2, y[i] - t[1][2]*dxt/2

			it,wt = generate_inds_weights(xt1,yt1,sim,NNN,i,ds)
			N = length(it)
			for j ∈ 1:N
				sim.link_start[cl+j-1] = i
				sim.link_stop[cl+j-1] = it[j]
				sim.link_weight[cl+j-1] = wt[j]/dxt^2
				sim.link_half_x[cl+j-1] = hxt1
				sim.link_half_y[cl+j-1] = hyt1
				if sides[1] ∈ (1,2)
					sim.link_y_bool[cl+j-1] = true
				else
					sim.link_x_bool[cl+j-1] = true
				end
				link_fiduciary[cl+j-1] = true
			end
	        cl += NNN

			xt2, yt2   = x[i] + t[1][1]*dxt  , y[i] + t[1][2]*dxt
			hxt2, hyt2 = x[i] + t[1][1]*dxt/2, y[i] + t[1][2]*dxt/2
			it,wt = generate_inds_weights(xt2,yt2,sim,NNN,i,ds)
			N = length(it)
			for j ∈ 1:N
				sim.link_start[cl+j-1] = i
				sim.link_stop[cl+j-1] = it[j]
				sim.link_weight[cl+j-1] = wt[j]/dxt^2
				sim.link_half_x[cl+j-1] = hxt2
				sim.link_half_y[cl+j-1] = hyt2
				if sides[1] ∈ (1,2)
					sim.link_y_bool[cl+j-1] = true
				else
					sim.link_x_bool[cl+j-1] = true
				end
				link_fiduciary[cl+j-1] = true
			end
		end
		verbose ? Threads.atomic_add!(pg_ref,1) : nothing
		if verbose && Threads.threadid()==Threads.nthreads()
			update!(pg,pg_ref[])
		end
	end
	verbose ? finish!(pg) : nothing
    return nothing
end

# KEY BLOCKS: LINKS CORNER POINTS via interpolation
function link_corner!(link_fiduciary::BitArray,self_fiduciary::BitArray,sim::Simulation,verbose::Bool)
	domains,x,y = sim.domains,sim.x,sim.y

	NNN = NNN_CORNER

	cl0 = findlast(link_fiduciary)
	cs0 = findlast(self_fiduciary)

	lattices = map(z->z.lattice,domains)
	DX, DY = map(z->z.dx,lattices), map(z->z.dy,lattices)
	DR, DΘ = map(z->z.dr,lattices), map(z->z.dθ,lattices)

	corner_inds = findall(sim.corner)
	for ci ∈ eachindex(corner_inds)
		i = corner_inds[ci]
		cl = cl0 + 2NNN*(ci-1) + 1
		cs = cs0 + 2(ci-1) + 1
		begin
			ds = sim.domain_index[i]
			n,t,d,sides = sim.normal[i],sim.tangent[i],sim.distance[i],sim.side[i]
			lds = lattices[ds]

			dxn = 2d[1]

			hxn1, hyn1 = x[i] - n[1][1]*dxn/2, y[i] - n[1][2]*dxn/2
			sim.self_index[cs] = i
			sim.self_weight[cs] = -1/dxn^2
			sim.self_half_x[cs] = hxn1
			sim.self_half_y[cs] = hyn1
			self_fiduciary[cs] = true
			if sides[1] ∈ (1,2)
				sim.self_x_bool[cs] = true
			else
				sim.self_y_bool[cs] = true
			end
			cs += 1

			if lds.type==:Cartesian
				dxn = max(dxn,min(DX[ds],DY[ds]))
			elseif lds.type==:Polar
				dxn = max(dxn,min(DR[ds],DΘ[ds]))
			else
				throw(LatticeError(lds.type))
			end

			xn2, yn2   = x[i] + n[1][1]*dxn  , y[i] + n[1][2]*dxn
			hxn2, hyn2 = x[i] + n[1][1]*dxn/2, y[i] + n[1][2]*dxn/2
			it,wt = generate_inds_weights(xn2,yn2,sim,NNN,i,ds)
			N = length(it)
			for j ∈ 1:N
				sim.link_start[cl+j-1] = i
				sim.link_stop[cl+j-1] = it[j]
				sim.link_weight[cl+j-1] = wt[j]/dxn^2
				sim.link_half_x[cl+j-1] = hxn2
				sim.link_half_y[cl+j-1] = hyn2
				link_fiduciary[cl+j-1] = true
				if sides[1] ∈ (1,2)
					sim.link_x_bool[cl+j-1] = true
				else
					sim.link_y_bool[cl+j-1] = true
				end
			end
			cl += NNN

			dxn = 2d[2]

			hxn1, hyn1 = x[i] - n[2][1]*dxn/2, y[i] - n[2][2]*dxn/2
			sim.self_index[cs] = i
			sim.self_weight[cs] = -1/dxn^2
			sim.self_half_x[cs] = hxn1
			sim.self_half_y[cs] = hyn1
			self_fiduciary[cs] = true
			if sides[1] ∈ (1,2)
				sim.self_y_bool[cs] = true
			else
				sim.self_x_bool[cs] = true
			end

			if lds.type==:Cartesian
				dxn = max(dxn,min(DX[ds],DY[ds]))
			elseif lds.type==:Polar
				dxn = max(dxn,min(DR[ds],DΘ[ds]))
			else
				throw(LatticeError(lds.type))
			end

			xn2, yn2   = x[i] + n[2][1]*dxn  , y[i] + n[2][2]*dxn
			hxn2, hyn2 = x[i] + n[2][1]*dxn/2, y[i] + n[2][2]*dxn/2
			it,wt = generate_inds_weights(xn2,yn2,sim,NNN,i,ds)
			N = length(it)
			for j ∈ 1:N
				sim.link_start[cl+j-1] = i
				sim.link_stop[cl+j-1] = it[j]
				sim.link_weight[cl+j-1] = wt[j]/dxn^2
				sim.link_half_x[cl+j-1] = hxn2
				sim.link_half_y[cl+j-1] = hyn2
				link_fiduciary[cl+j-1] = true
				if sides[1] ∈ (1,2)
					sim.link_y_bool[cl+j-1] = true
				else
					sim.link_x_bool[cl+j-1] = true
				end
			end
		end
	end
    return nothing
end

# trim the preallocated link arrays to the size actually used
function trim_link!(sim::Simulation,link_fiduciary::BitArray)
	deleteat!(sim.link_start,.!link_fiduciary)
	deleteat!(sim.link_stop,.!link_fiduciary)
	deleteat!(sim.link_weight,.!link_fiduciary)
	deleteat!(sim.link_half_x,.!link_fiduciary)
	deleteat!(sim.link_half_y,.!link_fiduciary)
	deleteat!(sim.link_half_r,.!link_fiduciary)

	link_x_bool = Array(sim.link_x_bool)
	deleteat!(link_x_bool,.!link_fiduciary)
	sim.link_x_bool[1:length(link_x_bool)] = link_x_bool
	resize!(sim.link_x_bool,length(link_x_bool))

	link_y_bool = Array(sim.link_y_bool)
	deleteat!(link_y_bool,.!link_fiduciary)
	sim.link_y_bool[1:length(link_y_bool)] = link_y_bool
	resize!(sim.link_y_bool,length(link_y_bool))

	link_r_bool = Array(sim.link_r_bool)
	deleteat!(link_r_bool,.!link_fiduciary)
	sim.link_r_bool[1:length(link_r_bool)] = link_r_bool
	resize!(sim.link_r_bool,length(link_r_bool))

	return nothing
end

# trim the preallocated self-linking arrays to the size actually used
function trim_self!(sim::Simulation,self_fiduciary::BitArray)
	deleteat!(sim.self_index,.!self_fiduciary)
	deleteat!(sim.self_weight,.!self_fiduciary)
	deleteat!(sim.self_half_x,.!self_fiduciary)
	deleteat!(sim.self_half_y,.!self_fiduciary)
	deleteat!(sim.self_half_r,.!self_fiduciary)

	self_x_bool = Array(sim.self_x_bool)
	deleteat!(self_x_bool,.!self_fiduciary)
	sim.self_x_bool[1:length(self_x_bool)] = self_x_bool
	resize!(sim.self_x_bool,length(self_x_bool))

	self_y_bool = Array(sim.self_y_bool)
	deleteat!(self_y_bool,.!self_fiduciary)
	sim.self_y_bool[1:length(self_y_bool)] = self_y_bool
	resize!(sim.self_y_bool,length(self_y_bool))

	self_r_bool = Array(sim.self_r_bool)
	deleteat!(self_r_bool,.!self_fiduciary)
	sim.self_r_bool[1:length(self_r_bool)] = self_r_bool
	resize!(sim.self_r_bool,length(self_r_bool))

	return nothing
end

# apply the local boundary conditions (Dirichlet, Neumann, etc)
function apply_local_bc!(sim::Simulation)
	apply_local_bc_link!(sim)
	apply_local_bc_self!(sim)
end

# apply local boundary conditions to the links
function apply_local_bc_link!(sim::Simulation)
	for i ∈ eachindex(sim.link_start)
		if sim.exterior[sim.link_stop[i]]
			new_stop = sim.exterior_from[sim.link_stop[i]]
			side = sim.exterior_side[sim.link_stop[i]]
			bcs = sim.domains[sim.domain_index[new_stop]].boundary.bcs
			shp = sim.domains[sim.domain_index[new_stop]].boundary.shape
			if typeof(shp)<:AbstractQuadrilateral
				if side≤length(shp.corners)
					if typeof(bcs[side])<:DirichletBC
						sim.link_weight[i] *= -1
		            elseif typeof(bcs[side])<:NeumannBC
		                nothing
					elseif typeof(bcs[side])<:noBC
						sim.link_weight[i] *= 0
		            end
				else
					if side-length(shp.corners)==1
						sides = [1,3]
					elseif side-length(shp.corners)==2
						sides = [2,3]
					elseif side-length(shp.corners)==3
						sides = [1,4]
					else
						sides = [2,4]
					end
					for j ∈ sides
						if typeof(bcs[j])<:DirichletBC
							sim.link_weight[i] *= -1
						elseif typeof(bcs[j])<:NeumannBC
							nothing
						elseif typeof(bcs[j])<:noBC
							sim.link_weight[i] *= 0
						end
					end
				end
			else
				if typeof(bcs[side])<:DirichletBC
					sim.link_weight[i] *= -1
	            elseif typeof(bcs[side])<:NeumannBC
	                nothing
				elseif typeof(bcs[side])<:noBC
					sim.link_weight[i] *= 0
	            end
			end
			sim.link_stop[i] = new_stop
		end
	end
end

# apply local boundary conditions to the self-links
function apply_local_bc_self!(sim::Simulation)
	for i ∈ eachindex(sim.self_index)
		si = sim.self_index[i]
		if !sim.corner[si]
			side = sim.side[si][1]
			bcs = sim.domains[sim.domain_index[si]].boundary.bcs
			if typeof(bcs[side])<:DirichletBC
				sim.self_weight[i] *= 2
	        elseif typeof(bcs[side])<:NeumannBC
	            sim.self_weight[i] *= 0
			elseif typeof(bcs[side])<:noBC
				nothing
	        end
		else
			if sim.self_x_bool[i]
				side = sim.side[si][1]
			else
				side = sim.side[si][2]
			end
			bcs = sim.domains[sim.domain_index[si]].boundary.bcs
			if typeof(bcs[side])<:DirichletBC
				sim.self_weight[i] *= 2
			elseif typeof(bcs[side])<:NeumannBC
				sim.self_weight[i] *= 0
			elseif typeof(bcs[side])<:noBC
				nothing
			end
		end
	end
end

# apply the nonlocal boundary conditions (not implemented yet)
function apply_nonlocal_bc!(sim::Simulation)
	return nothing
end

# get the interpolation indices and weights for (X,Y) on sites (x,y), with a maximal interpolation order of `NNN`
# do this by first finding the surrounding sites 2 layers deep (from Tessellations) and then finding the subset
# of N points which give the best-conditioned matrix and computing the weights for those points.
@inline function generate_inds_weights(X::Real,Y::Real,sim::Simulation,NNN::Int,site_index::Int,domain_index::Int)
	# t = sim.tessellation
	NNT = ceil(Int,sqrt(N_NEAREST_SITES_FOR_TESSELLATION/π))
	dx = sim.domains[domain_index].lattice.dx
	dy = sim.domains[domain_index].lattice.dy
	inds = hypot.(sim.x.-X,sim.y.-Y) .< max(dx,dy)*NNT
	deep = false
	tess_inds = get_surrounding_sites(Tessellation(sim.x[inds],sim.y[inds]),X,Y,deep)
	inds = findall(inds)[tess_inds]

	# inds = inds[hypot.(X.-x[inds],Y.-y[inds]).≤2NNN*max(dx,dy)]
	# inds = inds[sim.interior[inds]]
	# surface && (sum(map(z->z.boundary.shape(x[site_index],y[site_index]),sim.domains)) ≤ 1) ? setdiff!(inds,sim.not_this_domains_exterior[domain_index]) : nothing
	# surface=true
	# surface ?
	# setdiff!(inds,sim.not_this_domains_exterior[domain_index])# : setdiff!(inds,sim.exterior_inds)

	NN = min(length(inds),NNN)
	if N_CUBIC ≤ NN
		N = N_CUBIC
	elseif N_QUADRATIC ≤ NN
		N = N_QUADRATIC
	elseif N_BILINEAR ≤ NN
		N = N_BILINEAR
	elseif N_LINEAR ≤ NN
		N = N_LINEAR
	end
	# perm = sortperm(hypot.(X.-sim.x[inds],Y.-sim.y[inds]))
	# inds = inds[perm[1:N]]

	# find the combination of indices with the smallest condition number
	x = sim.x; y=sim.y
	C = collect(combinations(inds,N))
	conds = Ref(Inf); ii = Ref(1)
	for i ∈ eachindex(C)
		temp = cond(generate_interpolating_matrix(X,Y,x[C[i]],y[C[i]]))
		if temp < conds[]
			conds[] = temp
			ii[] = i
		end
	end
	inds = C[ii[]]

    weights = generate_weights(X,Y,x[inds],y[inds])
	return inds, weights
end

const POLYVEC3 = [1,0,0]
const POLYVEC4 = [1,0,0,0]
const POLYVEC6 = [1,0,0,0,0,0]
const POLYVEC10 = [1,0,0,0,0,0,0,0,0,0]
# "invert" the interpolating matrix or do a pseudo-inverse if singular
@inline function generate_weights(X::Real,Y::Real,x::Array,y::Array)
	N = length(x)
	interp_matrix = generate_interpolating_matrix(X,Y,x,y)
	try
		if N==3
			weights = interp_matrix\POLYVEC3
		elseif N==4
			weights = interp_matrix\POLYVEC4
		elseif N==6
			weights = interp_matrix\POLYVEC6
		elseif N==10
			weights = interp_matrix\POLYVEC10
		end
		return weights
	catch
		return weights = pinv(interp_matrix)[:,1]
	end
end

# matrix of powers of relative position, whose inverse contains the interpolation coefficients
@inline function generate_interpolating_matrix(X::Real,Y::Real,x::Array,y::Array)
	N = length(x)
	@assert N≤10 "only coded for N≤10"
	interp_matrix = Matrix{Float64}(undef,N,N)
	@inbounds for i ∈ 1:N
		N≥1  ? interp_matrix[1 ,i] = 1 : nothing
		N≥2  ? interp_matrix[2 ,i] = (y[i]-Y) : nothing
		N≥3  ? interp_matrix[3 ,i] = (x[i]-X) : nothing
		N≥4  ? interp_matrix[4 ,i] = (x[i]-X)*(y[i]-Y) : nothing
		N≥5  ? interp_matrix[5 ,i] = (y[i]-Y)^2 : nothing
		N≥6  ? interp_matrix[6 ,i] = (x[i]-X)^2 : nothing
		N≥7  ? interp_matrix[7 ,i] = (x[i]-X)^1*(y[i]-Y)^2 : nothing
		N≥8  ? interp_matrix[8 ,i] = (x[i]-X)^2*(y[i]-Y)^1 : nothing
		N≥9  ? interp_matrix[9 ,i] = (y[i]-Y)^3 : nothing
		N≥10 ? interp_matrix[10,i] = (x[i]-X)^3 : nothing
	end
	return interp_matrix
end

"""
	update_dielectric!(sim,[ε])

update simulation `sim` by recomputing the dielectric function or by giving it an explicit array `ε`
"""
function update_dielectric!(sim::Simulation,ε::Array)
	sim.ε[:] = ε
	sim.smoothed[] = false
	return nothing
end
update_dielectric!(sim::Simulation) = update_dielectric!(sim,simulation_dielectric.(Ref(sim),sim.x,sim.y))
"""
	update_pump!(sim,[F])

update simulation `sim` by recomputing the pump function or by giving it an explicit array `F`
"""
function update_pump!(sim::Simulation,F::Array)
	sim.F[:] = F
	sim.smoothed[] = false
	return nothing
end
update_pump!(sim::Simulation) = update_pump!(sim,simulation_pump.(Ref(sim),sim.x,sim.y))

"""
	update_sim!(sim)

update simulation `sim` by recomputing the dielectric and pump functions
"""
function update_sim!(sim::Simulation)
	update_dielectric!(sim)
	update_pump!(sim)
	return nothing
end

# used in the smoothing routing smooth!, it computes the dielectric as a function of position (x,y)
function simulation_dielectric!(ε,domains,x,y)
    for i ∈ eachindex(x)
        ε[i] = simulation_dielectric(domains,x[i],y[i])
    end
    return nothing
end
function simulation_dielectric(domains,x,y)
    dom = which_domain(domains,x,y)
    ε = dom==0 ? complex(NaN,NaN) : domains[dom].dielectric(x,y)
    return ε
end
simulation_dielectric(sim::Simulation,x,y) = simulation_dielectric(sim.domains,x,y)

# used in the smoothing routing smooth!, it computes the pump as a function of position (x,y)
function simulation_pump!(F,domains,x,y)
    for i ∈ eachindex(x)
        F[i] = simulation_pump(domains,x[i],y[i])
    end
    return nothing
end
function simulation_pump(domains,x,y)
    dom = which_domain(domains,x,y)
    F = dom==0 ? 0.0 : domains[dom].pump(x,y)
    return F
end
simulation_pump(sim::Simulation,x,y) = simulation_pump(sim.domains,x,y)

####################################################################################
# PLOTTING
####################################################################################

# plotting sim, with by argument, for external use
@recipe function f(sim::Simulation; by=nothing)
	if isnothing(by)
		layout --> (1,3)
	else
		layout --> 1
	end
	(sim, by)
end
@recipe function f(sim::Simulation,by::Union{Symbol,Function,Nothing})
	if isnothing(by)
		bys = [:real,:imag,:F]
		for i ∈ eachindex(bys)
			@series begin
				subplot --> i
				markersize --> MARKERSIZE_SCALE/sqrt(length(sim.x))/3
				(sim,bys[i],1)
			end
		end
	else
		markersize --> MARKERSIZE_SCALE/sqrt(length(sim.x))
		(sim,by,1)
	end
end

# elementary Simulation plot, for internal use
@recipe function f(sim::Simulation,by::Union{Symbol,Function},_1::Int)
	aspect_ratio --> 1
	legend --> false
	if by ∈ [:real,real,:Real,:Re,:re]
		z = real(sqrt.(sim.ε[sim.interior]))
		@series begin
			seriestype --> :scatter
			markershape --> MARKERSHAPE
			markerstrokealpha --> 0
	        clims --> (minimum(z), maximum(z))
			color --> :sequential
			marker_z --> z
			(sim.x[sim.interior], sim.y[sim.interior])
		end
		for i ∈ reverse(eachindex(sim.domains))
			@series begin
				title --> "Re n(x)"
				sim.domains[i].boundary
			end
		end
	elseif by ∈ [:imag,imag,:Imag,:Im,:im]
		z = imag(sqrt.(sim.ε[sim.interior]))
		@series begin
			seriestype --> :scatter
			markershape --> MARKERSHAPE
			markerstrokealpha --> 0
			clims --> (-maximum(abs.(z)), maximum(abs.(z)))
			color --> :diverging
			marker_z --> z
			(sim.x[sim.interior], sim.y[sim.interior])
		end
		for i ∈ reverse(eachindex(sim.domains))
			@series begin
				title --> "Im n(x)"
				sim.domains[i].boundary
			end
		end
	elseif by ∈ [:F,:pump,:f]
		z = sim.F[sim.interior]
		@series begin
			seriestype --> :scatter
			markershape --> MARKERSHAPE
			markerstrokealpha --> 0
			clims --> (-maximum(abs.(z)), maximum(abs.(z)))
			color --> :diverging
	        marker_z --> z
			(sim.x[sim.interior], sim.y[sim.interior])
		end
		for i ∈ reverse(eachindex(sim.domains))
			@series begin
				title --> "F(x)"
				sim.domains[i].boundary
			end
		end
	elseif by ∈ [:links, :link, :Links, :Link, :Hopping, :hopping, :adjacency, :Adjacency]
		@series (sim,true)
	else
		throw("unrecognized `by` keyword $by, must be one of :real, :imag, :F")
	end
end

@recipe function f(sim::Simulation,islink::Bool)
	legend --> false
	aspect_ratio --> 1
	S = sparse(sim.link_start,sim.link_stop,sim.link_weight,length(sim.x),length(sim.x),+)
	I,J,W = findnz(S)
	perm = sortperm(W,by=abs)
	I = I[perm]; J = J[perm]; W = abs.(W[perm])
	for i ∈ eachindex(I)
		inds = findall(isequal(i),I)
		W[inds] = W[inds]/mean(W[inds])
	end
	X = [sim.x[I] sim.x[J] NaN*ones(length(I))]
	Y = [sim.y[I] sim.y[J] NaN*ones(length(I))]
	Z = abs.([W W W])
	for i ∈ reverse(eachindex(sim.domains))
		@series begin
			ms --> LINK_SCALE_REDUCTION*MARKERSIZE_SCALE/sqrt(length(sim.x))
			sim.domains[i]
		end
	end
	@series begin
		clim --> (0,2)
		colorbar --> false
		seriestype --> :line
		aspect_ratio --> 1
		line_z --> Z'
		color--> :sequential_r
		(X',Y')
	end
end

end # module


text:
module Visualize_Assembly

using Modia3D

# Basic definitions
vmat1 = Modia3D.Material(color="LightBlue" , transparency=0.5)
vmat2 = Modia3D.Material(color="Red")

@assembly Bar(;Lx = 0.1, Ly=Lx/5, Lz=Ly) begin
   frame0 = Modia3D.Object3D(Modia3D.Solid(Modia3D.SolidBeam(Lx,Ly,Lz), nothing, vmat1))
   frame1 = Modia3D.Object3D(frame0; r=[-Lx/2, 0.0, 0.0])
   frame2 = Modia3D.Object3D(frame0; r=[ Lx/2, 0.0, 0.0])

   cyl    = Modia3D.Cylinder(Ly/2,1.2*Ly; material=vmat2)
   cyl1   = Modia3D.Object3D(frame1, cyl)
   cyl2   = Modia3D.Object3D(frame2, cyl)
end

Modia3D.visualizeAssembly!( Bar(Lx = 1.0, sceneOptions=Modia3D.SceneOptions(visualizeFrames=true,
                                                                            defaultFrameLength=0.3,
                                                                            enableContactDetection=false) ) )

println("... success of Visualize_Assembly.jl!")
end

ext: jl


text:
using PGFPlots
using GridInterpolations
using Colors
using ColorBrewer
using POMDPs
using AutomotivePOMDPs
using AutomotivePOMDPs: xv_to_state, yv_to_state
using POMDPModelTools
using LinearAlgebra
using ProgressMeter

struct SARSOPPolicy{P<:Policy}
    p::P
end

struct QMDPPolicy{P<:Policy}
    p::P
end

function POMDPs.value(policy::QMDPPolicy, b::SparseCat)
    val = zeros(length(actions(policy.pomdp)))
    for (i,s) in enumerate(b.vals)
        val += value(policy, s)*b.p[i]
    end
    return val
end

function POMDPs.value(p::SARSOPPolicy, b::SparseCat)
    policy = p.p
    pomdp = policy.pomdp
    vectors = policy.alphas
    n = length(vectors)
    utilities = zeros(length(actions(pomdp)), n)
    for i=1:n
        act = actionindex(pomdp, policy.action_map[i])
        for (j, s) in enumerate(b.vals)
            s_itps, w_itps = interpolate_state(s, pomdp)
            for (k, s_itp) in enumerate(s_itps)
                si = stateindex(pomdp, s_itp)
                utilities[act, i] += vectors[i][si]*b.probs[j]*w_itps[k]
            end
        end
    end
    return maximum(utilities, dims=2)
end

function POMDPs.value(p::QMDPPolicy, s::SingleOCState)
    policy = p.p
    pomdp = policy.pomdp
    val = zeros(length(actions(pomdp)))
    alphas = hcat(policy.alphas...)

    ego_grid = RectangleGrid(get_X_grid(pomdp), get_V_grid(pomdp)) #XXX must not allocate the grid at each function call, find better implementation
    ego_indices, ego_weights = interpolants(ego_grid, [s.ego.posG.x, s.ego.v])

    if !isinf(s.ped.v)
        ped_grid = RectangleGrid(get_Y_grid(pomdp), get_V_ped_grid(pomdp))
        ped_indices, ped_weights = interpolants(ped_grid, [s.ped.posG.y, s.ped.v])
    else
        ped_indices = Int64[0]
        ped_weights = Float64[1.0]
    end

    sum_weight = 0.
    for i=1:length(ego_indices)
        for j=1:length(ped_indices)
            xg, vg = ind2x(ego_grid, ego_indices[i])
            ego = xv_to_state(pomdp, xg, vg)
            if ped_indices[1] == 0
                ped = get_off_the_grid(pomdp)
            else
                yg, v_pedg = ind2x(ped_grid, ped_indices[j])
                ped = yv_to_state(pomdp, yg, v_pedg)
            end
            state = SingleOCState(collision_checker(ego, ped, pomdp.ego_type, pomdp.ped_type),ego,ped)
            si = stateindex(pomdp, state)
            sum_weight += ego_weights[i]*ped_weights[j]
            val .+= view(alphas, si, :)*ego_weights[i]*ped_weights[j]
        end
    end
    val ./= sum_weight
    return val
end

function get_XY_grid(pomdp::SingleOCPOMDP)
    env = pomdp.env
    rl = env.params.roadway_length
    cw = env.params.crosswalk_width
    lw = env.params.lane_width
    cl = env.params.crosswalk_length
    x_ped = rl/2
    Y = LinRange(-cl/4, cl/4, Int(floor(cl/2/pomdp.pos_res)) + 1)
    X = LinRange(5., rl/2 + 2*cw , Int(floor((rl/2 + 2*cw -5)/pomdp.pos_res)) + 1)
    return X, Y
end

function get_grid_data(p::QMDPPolicy, grid::RectangleGrid, pomdp::SingleOCPOMDP, v_ego::Float64, v_ped::Float64)
    policy = p.p
    values = zeros(length(actions(pomdp)),length(grid))
    alphas = hcat(policy.alphas...)
    for i=1:length(grid)
        x, y = ind2x(grid, i)
        ego = xv_to_state(pomdp, x, v_ego)
        ped = yv_to_state(pomdp, y, v_ped)
        s = SingleOCState(collision_checker(ego, ped, pomdp.ego_type, pomdp.ped_type), ego, ped)
        si = stateindex(pomdp, s)
        values[:,i] =  alphas[si, :]
    end
    return values
end

function get_belief(y::Float64, sig::Float64, n_bins::Int64)
    y_noise = LinRange(y-2*sig, y+2*sig, n_bins)
    weights = zeros(n_bins)
    for i=1:n_bins
        weights[i] = 1/(sig*sqrt(2*pi))*exp(-(y_noise[i] - y)^2/(2*sig^2))
    end
    normalize!(weights, 1)
    return y_noise, weights
end

function get_belief(x::Float64, y::Float64, sig::Float64, n_bins::Int64, v_ego::Float64, v_ped::Float64)
    y_noise = LinRange(y-2*sig, y+2*sig, n_bins)
    belief = SparseCat(Vector{SingleOCState}(undef, n_bins),Vector{Float64}(undef, n_bins))
    for i=1:n_bins
        ego = xv_to_state(pomdp, x, v_ego)
        ped = yv_to_state(pomdp, y_noise[i], v_ped)
        s = SingleOCState(collision_checker(ego, ped, pomdp.ego_type, pomdp.ped_type), ego, ped)
        weight = 1/(sig*sqrt(2*pi))*exp(-(y_noise[i] - y)^2/(2*sig^2))
        belief.vals[i] = s
        belief.probs[i] = weight
    end
    normalize!(belief.probs, 1)
    return belief
end

function compute_acts(pomdp::SingleOCPOMDP, p::QMDPPolicy, sig::Float64 = 0.01, v_ego::Float64 = 5., v_ped::Float64 = 1., n_bins::Int64 = 20, n_pts::Int64 = 100)
    X, Y = get_XY_grid(pomdp)
    X = X[1:21]
    grid = RectangleGrid(X, Y)
    values = get_grid_data(p, grid, pomdp, v_ego, v_ped)

    policy = p.p
    xrange = LinRange(X[1], X[end], n_pts)
    yrange = LinRange(Y[1], Y[end], n_pts)
    action_map = zeros(n_pts, n_pts)
    n_acts = length(actions(pomdp))
    pomdp_amap = ordered_actions(pomdp)

    for (i,x) in enumerate(xrange)
        for (j,y) in enumerate(yrange)
            y_noise, weights = get_belief(y, sig, n_bins)
            acts_xy = zeros(n_acts)
            for k=1:n_acts
                for l=1:n_bins
                    val = interpolate(grid, values[k,:], [x,y_noise[l]])
                    acts_xy[k] += weights[l]*val
                end
            end
            max, indmax = findmax(acts_xy)
            action_map[i,j] = pomdp_amap[indmax].acc
        end
    end
    acts = action_map
    acts = acts'
    acts = acts[end:-1:1,1:end]
    return acts
end

function compute_acts(pomdp::SingleOCPOMDP, policy::SARSOPPolicy, sig::Float64 = 0.01, v_ego::Float64 = 5., v_ped::Float64 = 1., n_bins::Int64 = 20, n_pts::Int64 = 100)
    X, Y = get_XY_grid(pomdp)
    X = X[1:21]
    grid = RectangleGrid(X, Y)

    xrange = LinRange(X[1], X[end], n_pts)
    yrange = LinRange(Y[1], Y[end], n_pts)
    action_map = zeros(n_pts, n_pts)
    n_acts = length(actions(pomdp))
    pomdp_amap = ordered_actions(pomdp)

    @showprogress for (i,x) in enumerate(xrange)
        for (j,y) in enumerate(yrange)
            belief = get_belief(x, y, sig, n_bins, v_ego, v_ped)
            max_val, ind_max = findmax(value(policy, belief))
            action_map[i,j] = pomdp_amap[ind_max].acc
        end
    end
    acts = action_map
    acts = acts'
    acts = acts[end:-1:1,1:end]
    return acts
end

function policy_plot(pomdp::SingleOCPOMDP, policy::Union{QMDPPolicy, SARSOPPolicy}, ;
                     sig::Float64 = 0.01, v_ego::Float64 = 5., v_ped::Float64 = 1., n_bins::Int64 = 20, n_pts::Int64 = 100)
    acts = compute_acts(pomdp, policy, sig, v_ego, v_ped, n_bins, n_pts)
    cmap = ColorMaps.RGBArrayMap(colormap("RdBu"))
    X, Y = get_XY_grid(pomdp)
    X = X[1:21]
    ax =Axis(Plots.Image(acts, (X[1], X[end]), (Y[1], Y[end]), colormap = cmap),
    xlabel = "Ego car position on the road",
    ylabel = "Pedestrian position along the crosswalk",
    title = "Policy Plot for a fixed ego velocity at $v_ego m/s and fixed pedestrian speed at $v_ped m/s")
    return ax
end

function value_plot(pomdp::SingleOCPOMDP, policy::SARSOPPolicy;
                    sig::Float64 = 0.01, v_ego::Float64 = 5., v_ped::Float64 = 1., n_bins::Int64 = 20, n_pts::Int64 = 100)
    X, Y = get_XY_grid(pomdp)
    X = X[1:21]
    grid = RectangleGrid(X, Y)

    xrange = LinRange(X[1], X[end], n_pts)
    yrange = LinRange(Y[1], Y[end], n_pts)
    vals = zeros(n_pts, n_pts)
    @showprogress for (i,x) in enumerate(xrange)
        for (j,y) in enumerate(yrange)
            belief = get_belief(x, y, sig, n_bins, v_ego, v_ped)
            vals[i,j] = maximum(value(policy, belief))
        end
    end
    vals = vals'
    vals = vals[end:-1:1,1:end]

    cmap = ColorMaps.RGBArrayMap(colormap("RdBu"))
    ax = Axis(Plots.Image(vals, (X[1], X[end]), (Y[1], Y[end]), colormap = cmap),
    xlabel = "Ego car position on the road",
    ylabel = "Pedestrian position along the crosswalk",
    title = "Probability of success for a fixed ego velocity at $v_ego m/s and fixed pedestrian speed at $v_ped m/s")
    return ax
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
# single HH neuron
using NeuralNets

include("../apma2821v/stimuli.jl")

net = SpikingNet()
add_layer!(net,HHLayer("single",1))
T = 300
dt = 0.01

measurements = ["single"=>["v","didSpike","m","h","n"]]
settings = RunSettings(T,dt,step_stim,measurements)

out = runSpikingNet(net,settings)

# analysis
using PyPlot

time_trace = 0:dt:T
v = out.data["single"]["v"]'
didSpike = out.data["single"]["didSpike"]'
println("Number of spikes:", length(find(didSpike)))

plot(time_trace[find(didSpike)],v[find(didSpike)],"r.")
plot(time_trace,v)
ylim([-80, 0])
show()

figure()
for var in ["m","h","n"]
    plot(time_trace,out.data["single"][var]')
end
legend(["m","h","n"])
show()

println("Press Enter to Continue...")
readline()


text:
```
oilrate(args...)  -->Plot

Recipe to Plot the Oil, Liquid Production.
Optionally you can make decline curve analysis, estimate a Forecast to either specific time or economic limit

Exponential Decline curve equation, Np estimation, Time to Economic Limit
`` { Q }_{ t }={ Q }_{ i }{ e }^{ -{ D }_{ i }t }``
``Np=\int _{ { t }_{ i } }^{ { t }_{ f } }{ { Q }_{ i }{ e }^{ -{ D }_{ i }t }dt }``
``Np=\frac { { -Q }_{ i } }{ { { D }_{ i } } } \left[ { e }^{ { -D }_{ i }{ t }_{ f } }-{ e }^{ { -D }_{ i }{ t }_{ i } } \right] ``
``{ t }_{ L }=\frac { -Ln\left( \frac { { Q }_{ L } }{ { Q }_{ i } }  \right)  }{ D } ``

The next table show the list of variables allowed:

|PropertyName|Args|Default|Input|Description
|---|---|---|---|---|
|Date|Mandatory|--|Array{Date,1}| Array of dates|
|OilRate|Mandatory|--|Array{Number,1}| Array of Rates|
|LiquidRate|Optional|--|Array{Number,1}| Array of Rates|
|Dlim|Optional|Dlim=false|Dlim=[Date(y,m,d)]|Initial Date to plot|
|mrange|Optional|mrange=false|mrange=number| Time interval in months to set xticks|
|dtf|Optional|dft="u yy"|dtf=DateFormat|Date format of xticks
|Liquid|Optional|Liquid=false|Liquid=Bool|If liquid array provided, choose if plot Liquid Rate|
|Dec|Optional|Dec=false|Dec=Bool|Set decline curve analysis to desired interval |
|Drange|Optional|Drange=false|Drange=[Date(y,m,d) Date(y,m,d)]|Decline curve range to fit|
|Forecast|Optional|Forecast=false|Forecast=Bool|Set to plot forecast estimation|
|Mforecast|Optional|Mforecast=12|Mforcast=Number|Set number of months to forecast <br> from last date of Declination range|
|EconLimit|Optional|EconLimit=false|EconLimit=Number|Set the Economic Limit to estimate time to forecast
|Ad|Optional|Ad=true|Ad=Bool|Set if apply Anormaly detection to decline curve|
```
@userplot OilRate

@recipe function f(h::OilRate; Dlim=false,Monthly=true, mrange=6, dft="u yy",
                                Liquid=false,
                                Dec=false, Drange=false,
                                Forecast=false, Mforecast=12,
                                Ad=true, EconLimit=false)
                  legend --> false
                  ylabel := "Rate [bbl/d]"

    if length(h.args)==2
        x, y = h.args
    elseif length(h.args)==3
        x, y, z = h.args
    end

    if Dec==true
        if Drange==false
            DaysDeclination=x
            RateDeclination=y
        else
            DaysDeclination=x[(x.>=Drange[1]) .& (x.<=Drange[2])]
            RateDeclination=y[(x.>=Drange[1]) .& (x.<=Drange[2])]
        end

        (S, DecTime, DecLine)=Declination(DaysDeclination,RateDeclination, AD=Ad)

        if Forecast==true
            if EconLimit!=false
                MforecastEL=round(Int,((log(EconLimit/S.Qi)/S.Di)/30.42))
                ForecastTime=DaysDeclination[end]:Dates.Month(1):DaysDeclination[1]+Dates.Month(MforecastEL)

            else
                ForecastTime=DaysDeclination[end]:Dates.Month(1):DaysDeclination[end]+Dates.Month(Mforecast)
            end
            ForecastLine=DeclinationForecast(S,ForecastTime)
            Np=(ForecastLine[1]-ForecastLine[end])/(-S.Di*1000)
            f=ForecastTime[end]+Dates.Month(1)
        else
            f=x[end]+Dates.Month(1)
        end
    else
        f=x[end]+Dates.Month(1)
    end

####
if Dlim==false
    if Monthly==true
      r = x[1]:Dates.Month(mrange):f
  else
     r = x[1]:Dates.Day(mrange):f
   end
    d = Dates.format.(r,dft)
else
    if Monthly==true
        r = Dlim[1]:Dates.Month(mrange):f
    else
        r= Dlim[1]:Dates.Day(mrange):f

    end
d = Dates.format.(r,dft)
end
xticks := (r,d)
if Dlim!=false
xlim := (Dates.value(Dlim[1]),Dates.value(Dlim[2]))
end
###
#=
    if Dlim==false
        r = x[1]:Dates.Month(mrange):f
        d = Dates.format.(r,dft)
    else
        r = Dlim[1]:Dates.Month(mrange):f
        d = Dates.format.(r,dft)
        xlim := (Dates.value(Dlim[1]),Dates.value(f))
    end
        xticks := (r,d)
=#
    #Oil rate Series
    @series begin
        seriestype := :path
        linewidth --> 3
        seriescolor --> :green

        label := "Oil Rate [bbl/d]"
        x, y
        end

    if Liquid==true
    #Liquid rate Series
    @series begin
        seriestype := :path
        linewidth --> 1
        seriescolor --> RGB(0.5, 0.5, 0.5)
        subplot := 1
        label --> "Liquid Rate [bbl/d]"
        x, z
        end
    end

    if Dec==true
        @series begin
            seriestype := :path
            linewidth --> 3
            seriescolor --> :red
            lab := "Di = $(round(S.Di*365/12,digits=2)) Monthly // $(round(S.Di*365,digits=2)) Annually"
            DecTime, DecLine
        end
    end

        if Forecast==true

        @series begin
            seriestype := :path
            linestyle --> :dash
            linewidth --> 1
            seriescolor --> :red
                lab := "$(Dates.format(ForecastTime[1],dft))--> $(Dates.format(ForecastTime[end],dft))// Np= $(round(Np,digits=1)) Mbbl))"
                ForecastTime, ForecastLine
        end

        end
end


text:
## materialize
mtab_materialize(x, ::Integer) = x
mtab_materialize(bc::Broadcasted, TN::Integer) = mtab_copy(instantiate(bc), TN)

mtb_materialize!(dest::Tuple{Vararg{AbstractArray}}, bc::Broadcasted, TN::Integer) =
    (mtb_materialize!(TupleDummy(dest), bc, TN); dest)

## copyto
mtab_copy(bc::Broadcasted{Style{Tuple}}, ::Integer) = tab_copy(bc)
mtab_copy(bc::FilledBC, ::Integer) = copy(bc)
mtab_copy(bc::Broadcasted, TN::Integer) = mtb_copyto!(toa_similar(bc), bc, TN).arrays


text:
module BanyanArrays

using Banyan

export Array, Vector, Matrix
export read_hdf5, write_hdf5
export ndims, size, length, eltype
export fill, zeros, ones, trues, falses
export map, mapslices, reduce, sort, sortlices

include("array.jl")

end # module

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


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

macro test_nothrow(ex)
    quote
        @test begin
            $(esc(ex))
            true
        end
    end
end

using IPython
using IPython: @compatattr, _setproperty!
using Compat


text:
using RNGPool

include("demo.jl")
include("lgModel.jl")

import Statistics: mean, var

setRNGs(12345)

println("Using a LG Model with n = 100")
model, lM, _ = setupLGModel(100, 1.0)
vs = CoupledConditionalSMC.couplingTimes(model, lM, 32, model.maxn, 100, :BS, true)
println(mean(vs), ", ", var(vs))

println("Using a LG Model with n = 200")
model, lM, _ = setupLGModel(200, 1.0)
vs = CoupledConditionalSMC.couplingTimes(model, lM, 32, model.maxn, 100, :BS, true)
println(mean(vs), ", ", var(vs))

println("Using a LG Model with n = 400")
model, lM, _ = setupLGModel(400, 1.0)
vs = CoupledConditionalSMC.couplingTimes(model, lM, 32, model.maxn, 100, :BS, true)
println(mean(vs), ", ", var(vs))


text:
map_str = ""
passports = []
open("day04/day04.in") do io
    global map_str, passports
    for line in readlines(io)
        line = line * " "
        if cmp(line, " ") == 0
            push!(passports, map_str)
            map_str = ""
        else
            map_str = map_str * line
        end
    end
end

pass_map = Dict(
    "byr" => 2, "iyr" => 3, "eyr" => 5, "hgt" => 7, "hcl" => 11, "ecl" => 13,
    "pid" => 17, "cid" => 19
    )
pass_valid = Int64(0)
for passport in passports
    global pass_valid
    checkmult = Int64(1)
    for item in split(passport, ' ')
        if length(item) > 2
            checkmult *= pass_map[item[1:3]]
        end
    end
    if checkmult == 9699690 || checkmult == 510510
        pass_valid += 1
    end
end
println(pass_valid)

max_stars_repo_name: detrin/advent-of-code

max_issues_repo_name: detrin/advent-of-code

max_forks_repo_name: detrin/advent-of-code


text:
export PTBLM
module PTBLM

const defdir = joinpath(Pkg.dir("MLDatasets"), "datasets/ptblm")

function traindata(dir=defdir)
    xs = readdata(dir, "ptb.train.txt")
    ys = makeys(xs)
    xs, ys
end

function testdata(dir=defdir)
    xs = readdata(dir, "ptb.test.txt")
    ys = makeys(xs)
    xs, ys
end

function readdata(dir, filename)
    mkpath(dir)
    url = "https://raw.githubusercontent.com/tomsercu/lstm/master/data"
    path = joinpath(dir, filename)
    isfile(path) || download("$(url)/$(filename)", path)
    lines = open(readlines, path)
    map(l -> Vector{String}(split(chomp(l))), lines)
end

function makeys(xs::Vector{Vector{String}})
    map(xs) do x
        y = copy(x)
        shift!(y)
        push!(y, "<eos>")
    end
end

end


text:
# using Revise
using Test
using Logging
using LarSurf

@testset "Triangulate quad" begin
    FV = [[1,2,3,4], [4,5,6,7]]
    FVtri = LarSurf.triangulate_quads(FV)
    @test size(FVtri[1],1) == 3
end


text:
const DEFAULT_VOIGT_ORDER = ([1], [1 3; 4 2], [1 6 5; 9 2 4; 8 7 3])
"""
    tovoigt(A::Union{SecondOrderTensor, FourthOrderTensor}; kwargs...)

Converts a tensor to "Voigt"-format.

Keyword arguments:
 - `offdiagscale`: determines the scaling factor for the offdiagonal elements.
   This argument is only applicable for `SymmetricTensor`s. `tomandel` can also
   be used for the "Mandel"-format which sets `offdiagscale = √2` for `SymmetricTensor`s,
   and is equivalent to `tovoigt` for `Tensor`s.
 - `order`: matrix of the linear indices determining the Voigt order. The default
   index order is `[11, 22, 33, 23, 13, 12, 32, 31, 21]`, corresponding to
   `order = [1 6 5; 9 2 4; 8 7 3]`.

See also [`tovoigt!`](@ref) and [`fromvoigt`](@ref).

```jldoctest
julia> tovoigt(Tensor{2,3}(1:9))
9-element Vector{Int64}:
 1
 5
 9
 8
 7
 4
 6
 3
 2

julia> tovoigt(SymmetricTensor{2,3}(1:6); offdiagscale = 2)
6-element Vector{Int64}:
  1
  4
  6
 10
  6
  4

julia> tovoigt(Tensor{4,2}(1:16))
4×4 Matrix{Int64}:
 1  13   9  5
 4  16  12  8
 3  15  11  7
 2  14  10  6
```
"""
function tovoigt end
@inline function tovoigt(A::Tensor{2, dim, T, M}; order=DEFAULT_VOIGT_ORDER[dim]) where {dim, T, M}
    @inbounds tovoigt!(Vector{T}(undef, M), A; order=order)
end
@inline function tovoigt(A::Tensor{4, dim, T, M}; order=DEFAULT_VOIGT_ORDER[dim]) where {dim, T, M}
    @inbounds tovoigt!(Matrix{T}(undef, Int(√M), Int(√M)), A; order=order)
end
@inline function tovoigt(A::SymmetricTensor{2, dim, T, M}; offdiagscale::T=one(T), order=DEFAULT_VOIGT_ORDER[dim]) where {dim, T, M}
    @inbounds tovoigt!(Vector{T}(undef, M), A; offdiagscale=offdiagscale, order=order)
end
@inline function tovoigt(A::SymmetricTensor{4, dim, T, M}; offdiagscale::T=one(T), order=DEFAULT_VOIGT_ORDER[dim]) where {dim, T, M}
    @inbounds tovoigt!(Matrix{T}(undef, Int(√M), Int(√M)), A; offdiagscale=offdiagscale, order=order)
end

"""
    tovoigt!(v::Array, A::Union{SecondOrderTensor, FourthOrderTensor}; kwargs...)

Converts a tensor to "Voigt"-format using the following index order:
`[11, 22, 33, 23, 13, 12, 32, 31, 21]`.

Keyword arguments:
 - `offset`: offset index for where in the array `A` the tensor should be stored.
   For 4th order tensors the keyword arguments are `offset_i` and `offset_j`,
   respectively. Defaults to `0`.
 - `offdiagscale`: determines the scaling factor for the offdiagonal elements.
   This argument is only applicable for `SymmetricTensor`s. `frommandel!` can also
   be used for the "Mandel"-format which sets `offdiagscale = √2` for `SymmetricTensor`s,
   and is equivalent to `fromvoigt!` for `Tensor`s.
 - `order`: matrix of the linear indices determining the Voigt order. The default
   index order is `[11, 22, 33, 23, 13, 12, 32, 31, 21]`.

See also [`tovoigt`](@ref) and [`fromvoigt`](@ref).

```jldoctest
julia> T = rand(Tensor{2,2})
2×2 Tensor{2, 2, Float64, 4}:
 0.590845  0.566237
 0.766797  0.460085

julia> x = zeros(4);

julia> tovoigt!(x, T)
4-element Vector{Float64}:
 0.5908446386657102
 0.4600853424625171
 0.5662374165061859
 0.7667970365022592

julia> x = zeros(5);

julia> tovoigt!(x, T; offset=1)
5-element Vector{Float64}:
 0.0
 0.5908446386657102
 0.4600853424625171
 0.5662374165061859
 0.7667970365022592
```
"""
function tovoigt! end
Base.@propagate_inbounds function tovoigt!(v::AbstractVector, A::Tensor{2, dim}; offset::Int=0, order=DEFAULT_VOIGT_ORDER[dim]) where {dim}
    for j in 1:dim, i in 1:dim
        v[offset + order[i, j]] = A[i, j]
    end
    return v
end
Base.@propagate_inbounds function tovoigt!(v::AbstractMatrix, A::Tensor{4, dim}; offset_i::Int=0, offset_j::Int=0, order=DEFAULT_VOIGT_ORDER[dim]) where {dim}
    for l in 1:dim, k in 1:dim, j in 1:dim, i in 1:dim
        v[offset_i + order[i, j], offset_j + order[k, l]] = A[i, j, k, l]
    end
    return v
end
Base.@propagate_inbounds function tovoigt!(v::AbstractVector{T}, A::SymmetricTensor{2, dim}; offdiagscale::T=one(T), offset::Int=0, order=DEFAULT_VOIGT_ORDER[dim]) where {T, dim}
    for j in 1:dim, i in 1:j
        v[offset + order[i, j]] = i == j ? A[i, j] : A[i, j] * offdiagscale
    end
    return v
end
Base.@propagate_inbounds function tovoigt!(v::AbstractMatrix{T}, A::SymmetricTensor{4, dim}; offdiagscale::T=one(T), offset_i::Int=0, offset_j::Int=0, order=DEFAULT_VOIGT_ORDER[dim]) where {T, dim}
    for l in 1:dim, k in 1:l, j in 1:dim, i in 1:j
        v[offset_i + order[i, j], offset_j + order[k, l]] =
            (i == j && k == l) ? A[i, j, k, l] :
            (i == j || k == l) ? A[i, j, k, l] * offdiagscale :
                                 A[i, j, k, l] * (offdiagscale * offdiagscale)
    end
    return v
end

@inline tomandel(A::SymmetricTensor{o, dim, T}; kwargs...) where{o,dim,T} = tovoigt(A; offdiagscale=√(T(2)), kwargs...)
@inline tomandel(A::Tensor; kwargs...) = tovoigt(A; kwargs...)

Base.@propagate_inbounds tomandel!(v::AbstractVecOrMat{T}, A::SymmetricTensor; kwargs...) where{T} = tovoigt!(v, A; offdiagscale=√(T(2)), kwargs...)
Base.@propagate_inbounds tomandel!(v::AbstractVecOrMat, A::Tensor; kwargs...) = tovoigt!(v, A; kwargs...)

"""
    fromvoigt(S::Type{<:AbstractTensor}, A::Array{T}; kwargs...)

Converts an array `A` stored in Voigt format to a Tensor of type `S`.

Keyword arguments:
 - `offset`: offset index for where in the array `A` the tensor starts. For 4th order
   tensors the keyword arguments are `offset_i` and `offset_j`, respectively.
   Defaults to `0`.
 - `offdiagscale`: determines the scaling factor for the offdiagonal elements.
   This argument is only applicable for `SymmetricTensor`s. `frommandel` can also
   be used for the "Mandel"-format which sets `offdiagscale = √2` for `SymmetricTensor`s,
   and is equivalent to `fromvoigt` for `Tensor`s.
 - `order`: matrix of the linear indices determining the Voigt order. The default
   index order is `[11, 22, 33, 23, 13, 12, 32, 31, 21]`, corresponding to
   `order = [1 6 5; 9 2 4; 8 7 3]`.

See also [`tovoigt`](@ref).

```jldoctest
julia> fromvoigt(Tensor{2,3}, 1.0:1.0:9.0)
3×3 Tensor{2, 3, Float64, 9}:
 1.0  6.0  5.0
 9.0  2.0  4.0
 8.0  7.0  3.0
```
"""
Base.@propagate_inbounds function fromvoigt(TT::Type{<: Tensor{2, dim}}, v::AbstractVector; offset::Int=0, order=DEFAULT_VOIGT_ORDER[dim]) where {dim}
    return TT(function (i, j); return v[offset + order[i, j]]; end)
end
Base.@propagate_inbounds function fromvoigt(TT::Type{<: Tensor{4, dim}}, v::AbstractMatrix; offset_i::Int=0, offset_j::Int=0, order=DEFAULT_VOIGT_ORDER[dim]) where {dim}
    return TT(function (i, j, k, l); return v[offset_i + order[i, j], offset_j + order[k, l]]; end)
end
Base.@propagate_inbounds function fromvoigt(TT::Type{<: SymmetricTensor{2, dim}}, v::AbstractVector{T}; offdiagscale::T = T(1), offset::Int=0, order=DEFAULT_VOIGT_ORDER[dim]) where {dim, T}
    return TT(function (i, j)
            i > j && ((i, j) = (j, i))
            i == j ? (return v[offset + order[i, j]]) :
                     (return v[offset + order[i, j]] / offdiagscale)
        end)
end
Base.@propagate_inbounds function fromvoigt(TT::Type{<: SymmetricTensor{4, dim}}, v::AbstractMatrix{T}; offdiagscale::T = T(1), offset_i::Int=0, offset_j::Int=0, order=DEFAULT_VOIGT_ORDER[dim]) where {dim, T}
    return TT(function (i, j, k, l)
            i > j && ((i, j) = (j, i))
            k > l && ((k, l) = (l, k))
            i == j && k == l ? (return v[offset_i + order[i, j], offset_j + order[k, l]]) :
            i == j || k == l ? (return v[offset_i + order[i, j], offset_j + order[k, l]] / offdiagscale) :
                               (return v[offset_i + order[i, j], offset_j + order[k, l]] / (offdiagscale * offdiagscale))
        end)
end

Base.@propagate_inbounds frommandel(TT::Type{<: SymmetricTensor}, v::AbstractVecOrMat{T}; kwargs...) where{T} = fromvoigt(TT, v; offdiagscale=√(T(2)), kwargs...)
Base.@propagate_inbounds frommandel(TT::Type{<: Tensor}, v::AbstractVecOrMat; kwargs...) = fromvoigt(TT, v; kwargs...)


text:
using Test
@testset "AutoGrad" begin
using Random; Random.seed!(1)
#TODO include("addto.jl")
@time include("base.jl")
@time include("cat.jl")
@time include("core.jl")
@time include("getindex.jl")
@time include("highorder.jl")
@time include("iterate.jl")
@time include("linearalgebra.jl")
#TODO include("macros.jl")
@time include("math.jl")
@time include("neuralnet.jl")
#TODO include("params.jl")
@time include("rosenbrock.jl")
#TODO include("show.jl")
@time include("sparse.jl")
@time include("specialfunctions.jl")
@time include("statistics.jl")
#TODO include("unbroadcast.jl")
end


text:
module TestDoctest

using Documenter
using Test
using Transducers

@testset "doctest" begin
    include("../docs/utils.jl")
    transducers_rm_examples()
    doctest(Transducers; manual=true)
end

end  # module


text:
edges = Dict()
isabove = Set()
f = open("input")
for line in map(chomp, eachline(f))
    m = match(r"([a-z]+) \(\d+\)(?: -> (.*))?", line)
    if m == nothing
        continue
    end
    bottom = m.captures[1]
    top = m.captures[2]
    if top == nothing
        edges[bottom] = []
    else
        above = split(top, ", ")
        map(p -> push!(isabove, p), above)
        edges[bottom] = above
    end
end
close(f)

for vertex in keys(edges)
    if !in(vertex, isabove)
        println("$vertex")
    end
end

max_stars_repo_name: nikklassen/advent-of-code


text:
function lagrangian_derivatives(mass_matrix, dynamics_bias, q, v)
    D1L = -1.0 * dynamics_bias(q, v)
    D2L = mass_matrix(q) * v
    return D1L, D2L
end

function dynamics(model::Model, mass_matrix, dynamics_bias, h, q0, q1, u1, w1, λ1, q2)
    # evalutate at midpoint
    qm1 = 0.5 * (q0 + q1)
    vm1 = (q1 - q0) / h[1]
    qm2 = 0.5 * (q1 + q2)
    vm2 = (q2 - q1) / h[1]

    D1L1, D2L1 = lagrangian_derivatives(mass_matrix, dynamics_bias, qm1, vm1)
    D1L2, D2L2 = lagrangian_derivatives(mass_matrix, dynamics_bias, qm2, vm2)

    d = 0.5 * h[1] * D1L1 + D2L1 + 0.5 * h[1] * D1L2 - D2L2 # variational integrator (midpoint)
    d .+= transpose(input_jacobian(model, qm2)) * u1             # control inputs
    d .+= λ1                                                # contact impulses
    d[(floating_base_dim(model) + 1):end] .-= model.friction_joint .* vm2[(floating_base_dim(model) + 1):end] # joint friction

    return d
end


text:
module Example

import Base.sho
import Base.+

export MyType, foo, +

struct MyType
    x
end

bar(x) = 2x
foo(a :: MyType) = bar(a.x) + 1
+(a :: MyType, b :: MyType) = a.x + b.x

show(io :: IO, a :: MyType)  = println(io, "MyType $(a.x)")

end

max_stars_repo_name: YinpuLi/SparsityBoosting

max_issues_repo_name: YinpuLi/SparsityBoosting

max_forks_repo_name: YinpuLi/SparsityBoosting


text:
# Various ways to query subgraphs.

using JLD, PyPlot, LinearAlgebra, SparseArrays

# Load data
file_dir = joinpath(Base.source_dir(),"./Entity.jld")
E = load(file_dir)["E"]
#E = loadassoc(joinpath(Base.source_dir(),"./Entity.jld2"))
E = logical(E)

# Compute entity (all facet pairs).
A = sqIn(E)
d = diag(getadj(A))
A = putAdj(A,getadj(A)-Diagonal(d))

# Compute normalized correlation.
i,j,v = findnz(getadj(A))
An = putAdj(A, sparse(i,j,v ./ min.(d[i],d[j])))

# Multi-facet queries.
x = "LOCATION/new york,"
p = StartsWith("PERSON/,")
printFull( (A[p,x] > 4) & (An[p,x] > 0.3))

# Triangles.
p0 = "PERSON/john kennedy,"

p1 = getrow(A[p,p0] + A[p0,p])
figure()
spy(A[p1,p1])

p2 = getrow( A[p1,p1] - (A[p,p0]+ A[p0,p]))
A[p2,p2] > 1


text:
using .BaseAgent
using Params
using Traffic.Simple
using Channel.Gilbert
using Util
using Params.QParams

export IndividualQ

type IndividualQ <: Agent
    s :: AgentState
    Q :: Array{Float64, 3}
    visit :: Array{Float64, 3}
    a :: Int
    state :: (Int, Int)
    P_tx :: Float64
    bitrate :: Float64
    status :: Status
    buf_interval :: Int
    beta_idle :: Float64
end

function IndividualQ(i, P, pos)
    Q = rand(int(Params.n_channel), P.buf_levels + 1, idle_action)
    Q *= Params.P_tx * Params.t_slot # a good initial randomization
    visit = zeros(Params.n_channel, P.buf_levels + 1, idle_action)
    IndividualQ(AgentState(i, P, pos), Q, visit, 0, (0, 0), 0, 0, Initialized,  div(Params.B + 1, P.buf_levels), P.beta_idle)
end

function policy!(a :: IndividualQ)
    # check whether we have packets first
    if rand() < epsilon
        a.a = rand(1:idle_action)
    else
        a.a = indmax(a.Q[a.s.chan, div(a.s.B_empty, a.buf_interval) + 1, :])
    end
    nothing
end

function BaseAgent.act(a :: IndividualQ, env, t)
    if a.status == Initialized
        if a.s.B_max - a.s.B_empty == 0
            a.a = -1
            return idle(a)
        end

        policy!(a)
        a.state = (a.s.chan, div(a.s.B_empty, a.buf_interval) + 1)
        a.visit[a.s.chan, a.state[2], a.a] += 1

        if a.a == idle_action
            return idle(a)
        end

        chan = int8(fld(a.a - 1, n_p_levels) + 1)
        return switch(a, chan)
    elseif a.status == Switched
        return sense(a, env, detect_traffic)
    elseif a.status == Sensed
        if detect_traffic(a.s, env.traffics[a.s.chan], a.s.t_remaining)
            return idle(a)
        else
            a.P_tx = Params.P_levels[(a.a - 1) % n_p_levels + 1]
            a.bitrate = capacity(env.channels[a.s.chan], a.P_tx, a.s.pos)
            pkgs_to_send = min(a.s.B_max - a.s.B_empty, floor(a.s.t_remaining * a.bitrate / Params.pkt_size))

            if pkgs_to_send == 0
                return idle(a)
            end
            return transmit!(a, a.P_tx, env, pkgs_to_send)
        end
    elseif a.status == Transmitted
        return idle(a)
    end
end

function alpha(a :: IndividualQ)
    return 0.2 + 0.8 / (1 + a.visit[a.state[1], a.state[2], a.a])
end

function BaseAgent.feedback(a :: IndividualQ, res :: Result, idle :: Bool = false, n_pkt :: Int = 0)
    feedback(a.s, res, idle, n_pkt)
    # if we didn't take any action that's worth learning, return immediately
    if a.a == -1
        return nothing
    end
    r :: Float64 = 0
    if idle
        r = - a.beta_idle * a.bitrate * Params.t_slot / a.s.E_slot
    elseif res == Success
        K = 1 # a.P_tx ^ 2 * Params.t_slot / a.bitrate
        r = K * Params.pkt_size * n_pkt / a.s.E_slot
    elseif res == Collision
        r = - beta_md * a.bitrate * Params.t_slot / a.s.E_slot
    elseif res == LostInChannel
        r = - beta_loss * a.bitrate * Params.t_slot / a.s.E_slot
    else
        return nothing
    end
    # UPDATE Q
    max_Q_next = maximum(a.Q[a.s.chan, floor(a.s.B_empty / a.buf_interval) + 1])
    Q_now = a.Q[a.state[1], a.state[2], a.a]
    a.Q[a.state[1], a.state[2], a.a] += alpha(a) * (r + discount * max_Q_next - Q_now)
    nothing
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

@testset "basic properties" begin

    @test typemin(Char) == Char(0)
    @test ndims(Char) == 0
    @test getindex('a', 1) == 'a'
    @test_throws BoundsError getindex('a', 2)
    # This is current behavior, but it seems questionable
    @test getindex('a', 1, 1, 1) == 'a'
    @test_throws BoundsError getindex('a', 1, 1, 2)

    @test 'b' + 1 == 'c'
    @test typeof('b' + 1) == Char
    @test 1 + 'b' == 'c'
    @test typeof(1 + 'b') == Char
    @test 'b' - 1 == 'a'
    @test typeof('b' - 1) == Char

    @test widen('a') === 'a'
    # just check this works
    @test_throws Base.CodePointError Base.code_point_err(UInt32(1))
end

@testset "ASCII conversion to/from Integer" begin
    numberchars = ['0', '1', '2', '3', '4', '5', '6', '7', '8', '9']
    lowerchars = ['a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o', 'p', 'q', 'r', 's', 't', 'u', 'v', 'w', 'x', 'y', 'z']
    upperchars = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X', 'Y', 'Z']
    plane1_playingcards = ['🂠', '🂡', '🂢', '🂣', '🂤', '🂥', '🂦', '🂧', '🂨', '🂩', '🂪', '🂫', '🂬', '🂭', '🂮']
    plane2_cjkpart1 = ['𠀀', '𠀁', '𠀂', '𠀃', '𠀄', '𠀅', '𠀆', '𠀇', '𠀈', '𠀉', '𠀊', '𠀋', '𠀌', '𠀍', '𠀎', '𠀏']

    testarrays = [numberchars; lowerchars; upperchars; plane1_playingcards; plane2_cjkpart1]

    #Integer(x::Char) = Int(x)
    #tests ASCII 48 - 57
    counter = 48
    for x in numberchars
        @test Integer(x) == counter
        counter += 1
    end

    #tests ASCII 65 - 90
    counter = 65
    for x in upperchars
        @test Integer(x) == counter
        counter += 1
    end

    #tests ASCII 97 - 122
    counter = 97
    for x in lowerchars
        @test Integer(x) == counter
        counter += 1
    end

    #tests Unicode plane 1: 127136 - 127150
    counter = 127136
    for x in plane1_playingcards
        @test Integer(x) == counter
        counter += 1
    end

    #tests Unicode plane 2: 131072 - 131087
    counter = 131072
    for x in plane2_cjkpart1
        @test Integer(x) == counter
        counter += 1
    end

    #convert(::Type{Char}, x::Float16) = char(convert(UInt32, x))
    #convert(::Type{Char}, x::Float32) = char(convert(UInt32, x))
    #convert(::Type{Char}, x::Float64) = char(convert(UInt32, x))
    for x = 1:9
        @test convert(Char, Float16(x)) == convert(Char, Float32(x)) == convert(Char, Float64(x)) == Char(x)
    end

    #size(c::Char) = ()
    for x in testarrays
        @test size(x) == ()
        @test_throws BoundsError size(x,0)
        @test size(x,1) == 1
    end

    #ndims(c::Char) = 0
    for x in testarrays
        @test ndims(x) == 0
    end

    #length(c::Char) = 1
    for x in testarrays
        @test length(x) == 1
    end

    #lastindex(c::Char) = 1
    for x in testarrays
        @test lastindex(x) == 1
    end

    #getindex(c::Char) = c
    for x in testarrays
        @test getindex(x) == x
    end

    #first(c::Char) = c
    for x in testarrays
        @test first(x) == x
    end

    #last(c::Char) = c
    for x in testarrays
        @test last(x) == x
    end

    #eltype(c::Char) = Char
    for x in testarrays
        @test eltype(x) == Char
    end

    #iterate(c::Char)
    for x in testarrays
        @test iterate(x)[1] == x
        @test iterate(x, iterate(x)[2]) == nothing
    end

    #isless(x::Char, y::Integer) = isless(UInt32(x), y)
    for x in upperchars
        @test isless(x, Char(91)) == true
    end

    for x in lowerchars
        @test isless(x, Char(123)) == true
    end

    for x in numberchars
        @test isless(x, Char(66)) == true
    end

    for x in plane1_playingcards
        @test isless(x, Char(127151)) == true
    end

    for x in plane2_cjkpart1
        @test isless(x, Char(131088)) == true
    end

    #isless(x::Integer, y::Char) = isless(x, UInt32(y))
    for x in upperchars
        @test isless(Char(64), x) == true
    end

    for x in lowerchars
        @test isless(Char(96), x) == true
    end

    for x in numberchars
        @test isless(Char(47), x) == true
    end

    for x in plane1_playingcards
        @test isless(Char(127135), x) == true
    end

    for x in plane2_cjkpart1
        @test isless(Char(131071), x) == true
    end

    @test !isequal('x', 120)
    @test convert(Signed, 'A') === Int32(65)
    @test convert(Unsigned, 'A') === UInt32(65)
end #end of let block

@testset "issue #14573" begin
    array = ['a', 'b', 'c'] + [1, 2, 3]
    @test array == ['b', 'd', 'f']
    @test eltype(array) == Char

    array = [1, 2, 3] + ['a', 'b', 'c']
    @test array == ['b', 'd', 'f']
    @test eltype(array) == Char

    array = ['a', 'b', 'c'] - [0, 1, 2]
    @test array == ['a', 'a', 'a']
    @test eltype(array) == Char
end

@testset "sprint, repr" begin
    @test sprint(show, "text/plain", '$') == "'\$': ASCII/Unicode U+0024 (category Sc: Symbol, currency)"
    @test repr('$') == "'\$'"
end

@testset "read incomplete character at end of stream or file" begin
    local file = tempname()
    local iob = IOBuffer([0xf0])
    local bytes(c::Char) = codeunits(string(c))
    @test bytes(read(iob, Char)) == [0xf0]
    @test eof(iob)
    try
        write(file, 0xf0)
        open(file) do io
            @test bytes(read(io, Char)) == [0xf0]
            @test eof(io)
        end
        let io = Base.Filesystem.open(file, Base.Filesystem.JL_O_RDONLY)
            @test bytes(read(io, Char)) == [0xf0]
            @test eof(io)
            close(io)
        end
    finally
        rm(file, force=true)
    end
end

@testset "overlong codes" begin
    function test_overlong(c::Char, n::Integer, rep::String)
        if isvalid(c)
            @test Int(c) == n
        else
            @test_throws Base.InvalidCharError UInt32(c)
        end
        @test sprint(show, c) == rep
        if Base.isoverlong(c)
            @test occursin(rep*": [overlong]", sprint(show, "text/plain", c))
        end
    end

    # TODO: use char syntax once #25072 is fixed
    test_overlong('\0', 0, "'\\0'")
    test_overlong("\xc0\x80"[1], 0, "'\\xc0\\x80'")
    test_overlong("\xe0\x80\x80"[1], 0, "'\\xe0\\x80\\x80'")
    test_overlong("\xf0\x80\x80\x80"[1], 0, "'\\xf0\\x80\\x80\\x80'")

    test_overlong('\x30', 0x30, "'0'")
    test_overlong("\xc0\xb0"[1], 0x30, "'\\xc0\\xb0'")
    test_overlong("\xe0\x80\xb0"[1], 0x30, "'\\xe0\\x80\\xb0'")
    test_overlong("\xf0\x80\x80\xb0"[1], 0x30, "'\\xf0\\x80\\x80\\xb0'")

    test_overlong('\u8430', 0x8430, "'萰'")
    test_overlong("\xf0\x88\x90\xb0"[1], 0x8430, "'\\xf0\\x88\\x90\\xb0'")
end

# create a new AbstractChar type to test the fallbacks
primitive type ASCIIChar <: AbstractChar 8 end
ASCIIChar(c::UInt8) = reinterpret(ASCIIChar, c)
ASCIIChar(c::UInt32) = ASCIIChar(UInt8(c))
Base.codepoint(c::ASCIIChar) = reinterpret(UInt8, c)

@testset "abstractchar" begin
    @test AbstractChar('x') === AbstractChar(UInt32('x')) === 'x'

    @test isascii(ASCIIChar('x'))
    @test ASCIIChar('x') < 'y'
    @test ASCIIChar('x') == 'x' === Char(ASCIIChar('x')) === convert(Char, ASCIIChar('x'))
    @test ASCIIChar('x')^3 == "xxx"
    @test repr(ASCIIChar('x')) == "'x'"
    @test string(ASCIIChar('x')) == "x"
    @test_throws MethodError write(IOBuffer(), ASCIIChar('x'))
    @test_throws MethodError read(IOBuffer('x'), ASCIIChar)
end

@testset "ncodeunits(::Char)" begin
    # valid encodings
    @test ncodeunits('\0')       == 1
    @test ncodeunits('\x1')      == 1
    @test ncodeunits('\x7f')     == 1
    @test ncodeunits('\u80')     == 2
    @test ncodeunits('\uff')     == 2
    @test ncodeunits('\u7ff')    == 2
    @test ncodeunits('\u800')    == 3
    @test ncodeunits('\uffff')   == 3
    @test ncodeunits('\U10000')  == 4
    @test ncodeunits('\U10ffff') == 4
    # invalid encodings
    @test ncodeunits(reinterpret(Char, 0x80_00_00_00)) == 1
    @test ncodeunits(reinterpret(Char, 0x01_00_00_00)) == 1
    @test ncodeunits(reinterpret(Char, 0x00_80_00_00)) == 2
    @test ncodeunits(reinterpret(Char, 0x00_01_00_00)) == 2
    @test ncodeunits(reinterpret(Char, 0x00_00_80_00)) == 3
    @test ncodeunits(reinterpret(Char, 0x00_00_01_00)) == 3
    @test ncodeunits(reinterpret(Char, 0x00_00_00_80)) == 4
    @test ncodeunits(reinterpret(Char, 0x00_00_00_01)) == 4
end

@testset "reinterpret(Char, ::UInt32)" begin
    for s = 0:31
        u = one(UInt32) << s
        @test reinterpret(UInt32, reinterpret(Char, u)) === u
    end
end

@testset "broadcasting of Char" begin
    @test identity.('a') == 'a'
end

max_stars_repo_name: milesfrain/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
contains(a::UnitRange{Int}, b::UnitRange{Int}) = intersect(a, b) == b

function complimentary(a::UnitRange{Int}, b::UnitRange{Int})
    lower = min(first(a), first(b)):(max(first(a), first(b)) - 1)
    upper = (min(last(a), last(b)) + 1):max(last(a), last(b))
    return lower, upper
end

function pick_available_subregion(region::UnitRange{Int}, size::Int)
    rgnstart = pick_start_point(region, size)
    rgnstop = rgnstart + size - 1
    return rgnstart:rgnstop
end

function pick_start_point(region::UnitRange{Int}, size::Int)
    if length(region) < size
        throw(ArgumentError("region must be longer or equal to size."))
    end
    laststart = last(region) - size + 1
    return rand(first(region):laststart)
end


text:
module Geo
import SymPy

export
    Point2D,
    Point4D,
    Point4To2D,
    Distance2D,
    Geometry,
    Length2XY,
    Length2XYθ,
    ArcCenter,
    GeosIntersect
include("Common.jl")
include("LineArc.jl")
include("LineLine.jl")
include("ArcArc.jl")

end # module


text:
export Multiplication, *

"""
Description:

    For continuous random variables, the multiplication node acts
    as a (matrix) multiplication constraint, with node function

    f(out, in1, a) = δ(out - a*in1)

Interfaces:

    1. out
    2. in1
    3. a

Construction:

    Multiplication(out, in1, a, id=:some_id)
"""
mutable struct Multiplication <: DeltaFactor
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol,Interface}

    function Multiplication(out, in1, a; id=generateId(Multiplication))
        @ensureVariables(out, in1, a)
        self = new(id, Array{Interface}(undef, 3), Dict{Int,Interface}())
        addNode!(currentGraph(), self)
        self.i[:out] = self.interfaces[1] = associate!(Interface(self), out)
        self.i[:in1] = self.interfaces[2] = associate!(Interface(self), in1)
        self.i[:a] = self.interfaces[3] = associate!(Interface(self), a)

        return self
    end
end

slug(::Type{Multiplication}) = "×"

function *(a::Variable, in1::Variable)
    out = Variable()
    Multiplication(out, in1, a)
    return out
end

function *(in1::Variable, a)
    @ensureVariables(a)
    *(a, in1)
end

function *(a, in1::Variable)
    @ensureVariables(a)
    *(a, in1)
end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function planck_freq(frequency::T, temperature::T) where {T<:AbstractFloat}
    c1  = 1.47449944028424e-50 # = 2*h/c*c
    c2  = 4.79924466221135e-11 # = h/k
    return c1*frequency^3/expm1(c2*frequency/temperature)
end

"""
    planck_freq(frequency, temperature) -> black_body_flux

### Purpose ###

Calculate the flux of a black body per unit frequency.

### Explanation ###

Return the spectral radiance of a black body per unit frequency using [Planck's
law](https://en.wikipedia.org/wiki/Planck%27s_law)

`` B_\\nu(\\nu, T) = \\frac{2h\\nu ^3}{c^2} \\frac{1}{e^\\frac{h\\nu}{k_\\mathrm{B}T} - 1} ``

### Arguments ###

* `frequency`: frequency at which the flux is to be calculated, in Hertz.
* `temperature`: the equilibrium temperature of the black body, in Kelvin.

### Output ###

The spectral radiance of the black body, in units of W/(sr·m²·Hz).

### Example ###

Plot the spectrum of a black body in \$[10^{12}, 10^{15.4}]\$ Hz at \$8000\$ K.
Use [PyPlot.jl](https://github.com/JuliaPlots/Plots.jl/) for plotting.

```julia
using PyPlot
frequency = exp10.(range(12, stop=15.4, length=1000));
temperature = ones(frequency)*8000;
flux = planck_freq.(frequency, temperature);
plot(frequency, flux)
```

### Notes ###

`planck_wave` calculates the flux of a black body per unit wavelength.
"""
planck_freq(f::Real, t::Real) = planck_freq(promote(float(f), float(t))...)


text:
@testset "Manopt Cyclic Proximal Point" begin
    using Dates
    n = 100
    N = PowerManifold(Circle(), n)
    f = artificial_S1_signal(n)
    F = x -> costL2TV(N, f, 0.5, x)
    proxes = ((λ, x) -> prox_distance(N, λ, f, x), (λ, x) -> prox_TV(N, 0.5 * λ, x))
    o = cyclic_proximal_point(
        N,
        F,
        proxes,
        f;
        λ=i -> π / (2 * i),
        stopping_criterion=StopAfterIteration(100),
        debug=[
            DebugIterate(), " ", DebugCost(), " ", DebugProximalParameter(), "\n", 10000
        ],
        record=[RecordProximalParameter(), RecordIterate(f), RecordCost()],
        return_options=true,
    )
    fR = get_solver_result(o)
    fR2 = cyclic_proximal_point(
        N, F, proxes, f; λ=i -> π / (2 * i), stopping_criterion=StopAfterIteration(100)
    )
    @test fR == fR2
    rec = get_record(o)
    @test F(f) > F(fR)
    #
    o = CyclicProximalPointOptions(f, StopAfterIteration(1), i -> π / (2 * i))
    p = ProximalProblem(N, F, proxes, [1, 2])
    @test_throws ErrorException get_proximal_map(p, 1.0, f, 3)
    @test_throws ErrorException ProximalProblem(N, F, proxes, [1, 2, 2])
end