text:
#
# # Credits
# Author: Pawan Bharadwaj
#         (bharadwaj.pawan@gmail.com)
#
# * original code in FORTRAN90: March 2013
# * modified: 11 Sept 2013
# * major update: 25 July 2014
# * code optimization with help from Jan Thorbecke: Dec 2015
# * rewritten in Julia: June 2017
# * added parrallelization over supersources in Julia: July 2017
# * efficient parrallelization using distributed arrays: Sept 2017
# * optimized memory allocation: Oct 2017
#

"""
Forward modeling using a finite-difference simulation of the acoustic wave-equation.
"""
struct Fdtd end

"""
Linearized forward modeling using a finite-difference simulation of the acoustic wave-equation.
"""
struct FdtdBorn end

global const npml = 50
global const nlayer_rand = 0

#
#As forward modeling method, the
#finite-difference method is employed.
#It uses a discrete version of the two-dimensional isotropic acoustic wave equation.
#
#```math
#\pp[\tzero] - \pp[\tmo] = \dt \mB \left({\partial_x\vx}[\tmh]
# + \partial_z \vz[\tmh]  + \dt\sum_{0}^{\tmo}\sfo\right)
# ```
# ```math
#\pp[\tpo] - \pp[\tzero] = \dt \mB \left(\partial_x \vx[\tph]
# + {\partial_z \vz}[\tph]  + \dt\sum_{0}^{\tzero}\sfo\right)
# ```

"""
Modelling parameters common for all supersources
# Keyword Arguments that are modified by the method (some of them are returned as well)

* `gradient::Vector{Float64}=Medium(model.mgrid)` : gradient model modified only if `gmodel_flag`
* `TDout::Vector{Data.TD}=[Data.TD_zeros(rfields,tgridmod,ageom[ip]) for ip in 1:length(findn(rflags))]`
* `illum::Array{Float64,2}=zeros(length(model.mgrid[1]), length(model.mgrid[2]))` : source energy if `illum_flag`
* `boundary::Array{Array{Float64,4},1}` : stored boundary values for first propagating wavefield
* `snaps::Array{Float64,4}=zeros(length(model.mgrid[1]),length(model.mgrid[2]),length(tsnaps),ageom[1].nss)` :snapshots saved at `tsnaps`

# Return (in order)

* modelled data for each propagating wavefield as `Vector{TD}`
* stored boundary values of the first propagating wavefield as `Array{Array{Float64,4},1}` (use for backpropagation)
* final conditions of the first propagating wavefield as `Array{Float64,4}` (use for back propagation)
* gradient model as `Seismic`
* stored snaps shots at tsnaps as Array{Float64,4}
"""
mutable struct PFdtdc{T}
	jobname::Symbol
	attrib_mod::T
	npw::Int64
	activepw::Vector{Int64}
	exmodel::Medium
	model::Medium
	ageom::Vector{AGeom}
	srcwav::Vector{SrcWav}
	abs_trbl::Vector{Symbol}
	sfields::Vector{Vector{Symbol}}
	isfields::Vector{Vector{Int64}}
	sflags::Vector{Int64}
	irfields::Vector{Int64}
	rflags::Vector{Int64}
	δt::Float64
	δxI::Float64
	δzI::Float64
	nx::Int64
	nz::Int64
	nt::Int64
	δtI::Float64
	δx24I::Float64
	δz24I::Float64
	a_x::Vector{Float64}
	b_x::Vector{Float64}
	k_xI::Vector{Float64}
	a_x_half::Vector{Float64}
	b_x_half::Vector{Float64}
	k_x_halfI::Vector{Float64}
	a_z::Vector{Float64}
	b_z::Vector{Float64}
	k_zI::Vector{Float64}
	a_z_half::Vector{Float64}
	b_z_half::Vector{Float64}
	k_z_halfI::Vector{Float64}
	modtt::Matrix{Float64}
	modttI::Matrix{Float64} # just storing inv(modtt) for speed
	modrr::Matrix{Float64}
	modrrvx::Matrix{Float64}
	modrrvz::Matrix{Float64}
	δmodtt::Matrix{Float64}
	δmodrr::Matrix{Float64}
	δmodrrvx::Matrix{Float64}
	δmodrrvz::Matrix{Float64}
	δmod::Vector{Float64} # perturbation vector (KI, rhoI)
	gradient::Vector{Float64}  # output gradient vector w.r.t (KI, rhoI)
	grad_modtt_stack::SharedArrays.SharedArray{Float64,2} # contains gmodtt
	grad_modrrvx_stack::SharedArrays.SharedArray{Float64,2}
	grad_modrrvz_stack::SharedArrays.SharedArray{Float64,2}
	grad_modrr_stack::Array{Float64,2}
	illum_flag::Bool
	illum_stack::SharedArrays.SharedArray{Float64,2}
	backprop_flag::Int64
	snaps_flag::Bool
	itsnaps::Vector{Int64}
	gmodel_flag::Bool
	ibx0::Int64
	ibz0::Int64
	ibx1::Int64
	ibz1::Int64
	isx0::Int64
	isz0::Int64
	datamat::SharedArrays.SharedArray{Float64,3}
	data::Vector{Data}
	verbose::Bool
end

"""
Modelling parameters per every supersource for each worker
"""
mutable struct PFdtdss
	iss::Int64
	wavelets::Matrix{Matrix{Float64}}
	ssprayw::Vector{Matrix{Float64}}
	records::Vector{Array{Float64,3}}
	rinterpolatew::Vector{Matrix{Float64}}
	isx1::Vector{Vector{Int64}}
	isx2::Vector{Vector{Int64}}
	isz1::Vector{Vector{Int64}}
	isz2::Vector{Vector{Int64}}
	irx1::Vector{Vector{Int64}}
	irx2::Vector{Vector{Int64}}
	irz1::Vector{Vector{Int64}}
	irz2::Vector{Vector{Int64}}
        boundary::Vector{Array{Float64,3}}
	snaps::Array{Float64,3}
	illum::Matrix{Float64}
	grad_modtt::Matrix{Float64}
	grad_modrrvx::Matrix{Float64}
	grad_modrrvz::Matrix{Float64}
end

"""
PFdtdeters per every worker, not necessarily for every supersource.
Note that a single worker can take care of multiple supersources.
"""
mutable struct PFdtdp
	ss::Vector{PFdtdss}
	p::Vector{Array{Float64,3}}
	pp::Vector{Array{Float64,3}}
	ppp::Vector{Array{Float64,3}}
	dpdx::Vector{Array{Float64,3}}
	dpdz::Vector{Array{Float64,3}}
	memory_dp_dx::Vector{Array{Float64,2}}
	memory_dp_dz::Vector{Array{Float64,2}}
	memory_dvx_dx::Vector{Array{Float64,2}}
	memory_dvz_dz::Vector{Array{Float64,2}}
	born_svalue_stack::Matrix{Float64}
end

mutable struct PFdtd{T}
	sschunks::Vector{UnitRange{Int64}} # how supersources are distributed among workers
	p::DistributedArrays.DArray{PFdtdp,1,PFdtdp} # distributed parameters among workers
	c::PFdtdc{T} # common parameters
end

"""
```julia
pa = SeisForwExpt(attrib_mod; ageom, srcwav, medium, tgrid);
```
Method to create an instance of `SeisForwExpt`.
The output of this method can be used as an input to the in-place method `update!`, to actually perform a
finite-difference modeling.

# Keyword Arguments

* `attrib_mod` : attribute to choose the type of modeling. Choose from
  * `=Fdtd()` for full wavefield modeling  (finite-difference simulation of the acoustic wave-equation)
  * `=FdtdBorn()` for Born modeling
* `model::Medium` : medium parameters
* `tgrid` : modeling time grid, maximum time in `tgrid`should be greater than or equal to maximum source time, same sampling interval as in `srcwav`
* `ageom` :  acquisition
* `srcwav` : source wavelets

# Optional Keyword Arguments

* `sflags=2` : source related flags
  * `=0` inactive sources
  * `=1` sources with injection rate
  * `=2` volume injection sources
  * `=3` sources input after time reversal (use only during backpropagation)
* `rflags=1` : receiver related flags
  * `=0` receivers do not record (or) inactive receivers
  * `=1` receivers are active only for the second propagating wavefield
* `rfields=[:P]` : multi-component receiver flag. Choose `Vector{Symbol}`
  * `=[:P]` record pressure
  * `=[:Vx]` record horizontal component of particle velocity
  * `=[:Vz]` record vertical component of particle velocity
  * `=[:P, :Vx]` record both pressure and velocity
* `tsnaps` : store snaps at these modeling times (defaults to half time)
  * `=[0.1,0.2,0.3]` record at these instances of tgrid
* `snaps_flag::Bool=false` : return snapshots or not
* `verbose::Bool=false` : verbose flag
"""
SeisForwExpt(attrib_mod::Union{Fdtd,FdtdBorn},args1...;args2...)=PFdtd(attrib_mod,args1...;args2...)

function initialize!(pap::PFdtdp)
	reset_per_ss!(pap)
	for issp in 1:length(pap.ss)
		pass=pap.ss[issp]
		for i in 1:length(pass.records)
			fill!(pass.records[i],0.0)
		end
		fill!(pass.snaps,0.0)
		fill!(pass.illum,0.0)
		fill!(pass.grad_modtt,0.0)
		fill!(pass.grad_modrrvx,0.0)
		fill!(pass.grad_modrrvz,0.0)
	end
end

function iszero_boundary(pa)
	result=false
	@sync begin
		for (ip, p) in enumerate(procs(pa.p))
			@async remotecall_wait(p) do
				pap=localpart(pa.p)
				for issp in 1:length(pap.ss)
					pass=pap.ss[issp]
					for i in 1:length(pass.boundary)
						result = result | iszero(pass.boundary[i])
					end
				end
			end
		end
	end
	return result
end

function initialize_boundary!(pa)
	# zero out results stored per worker
	@sync begin
		for (ip, p) in enumerate(procs(pa.p))
			@async remotecall_wait(p) do
				pap=localpart(pa.p)
				for issp in 1:length(pap.ss)
					pass=pap.ss[issp]
					for i in 1:length(pass.boundary)
						fill!(pass.boundary[i],0.0)
					end
				end
			end
		end
	end
end

function initialize!(pac::PFdtdc)
	fill!(pac.gradient,0.0)
	fill!(pac.grad_modtt_stack,0.0)
	fill!(pac.grad_modrrvx_stack,0.0)
	fill!(pac.grad_modrrvz_stack,0.0)
	fill!(pac.grad_modrr_stack,0.0)
	fill!(pac.illum_stack,0.0)
	for dat in pac.data
		fill!(dat, 0.0)
	end
	fill!(pac.datamat,0.0)
end

function reset_per_ss!(pap::PFdtdp)
	for ipw in 1:length(pap.p)
		fill!(pap.born_svalue_stack,0.0)
		fill!(pap.p[ipw],0.0)
		fill!(pap.pp[ipw],0.0)
		fill!(pap.ppp[ipw],0.0)
		fill!(pap.dpdx[ipw],0.0)
		fill!(pap.dpdz[ipw],0.0)
		fill!(pap.memory_dp_dz[ipw],0.0)
		fill!(pap.memory_dp_dx[ipw],0.0)
		fill!(pap.memory_dvx_dx[ipw],0.0)
		fill!(pap.memory_dvz_dz[ipw],0.0)
	end
end

"""
Primary method to generate Expt variable when Fdtd() and FdtdBorn() are used.

# Some internal arguments
* `jobname::String` : name
* `npw::Int64=1` : number of independently propagating wavefields in `model`
* `backprop_flag::Bool=Int64` : save final state variables and the boundary conditions for later use
  * `=1` save boundary and final values in `boundary`
  * `=-1` use stored values in `boundary` for back propagation
* `gmodel_flag=false` : flag that is used to output gradient; there should be atleast two propagating wavefields in order to do so: 1) forward wavefield and 2) adjoint wavefield
* `illum_flag::Bool=false` : flag to output wavefield energy or source illumination; it can be used as preconditioner during inversion
* `abs_trbl::Vector{Symbol}=[:top, :bottom, :right, :left]` : use absorbing PML boundary conditions or not
  * `=[:top, :bottom]` apply PML conditions only at the top and bottom of the model
  * `=[:bottom, :right, :left]` top is reflecting
"""
function PFdtd(attrib_mod;
	jobname::Symbol=:forward_propagation,
	npw::Int64=1,
	medium::Medium=nothing,
	abs_trbl::Vector{Symbol}=[:top, :bottom, :right, :left],
	tgrid::StepRangeLen=nothing,
	ageom::Union{AGeom,Vector{AGeom}}=nothing,
	srcwav::Union{SrcWav,Vector{SrcWav}}=nothing,
	sflags::Union{Int,Vector{Int}}=fill(2,npw),
	rflags::Union{Int,Vector{Int}}=fill(1,npw),
	rfields::Vector{Symbol}=[:P],
	backprop_flag::Int64=0,
	gmodel_flag::Bool=false,
	illum_flag::Bool=false,
	tsnaps::Vector{Float64}=fill(0.5*(tgrid[end]+tgrid[1]),1),
	snaps_flag::Bool=false,
	verbose::Bool=false,
	nworker=nothing)

	# convert to vectors
	if(typeof(ageom)==AGeom); ageom=[ageom]; end
	if(typeof(srcwav)==SrcWav); srcwav=[srcwav]; end
	if(typeof(sflags)==Int); sflags=[sflags]; end
	if(typeof(rflags)==Int); rflags=[rflags]; end

	# alias
	model=medium
	tgridmod=tgrid

	#println("********PML Removed*************")
	#abs_trbl=[:null]

	# check sizes and errors based on input
	#(length(TDout) ≠ length(findn(rflags))) && error("TDout dimension")
	(length(ageom) ≠ npw) && error("ageom dimension")
	(length(srcwav) ≠ npw) && error("srcwav dimension")
	(length(sflags) ≠ npw) && error("sflags dimension")
	(length(rflags) ≠ npw) && error("rflags dimension")
	(maximum(tgridmod) < maximum(srcwav[1][1].grid)) && error("modeling time is less than source time")
	#(any([getfield(TDout[ip],:tgrid).δx < tgridmod.δx for ip=1:length(TDout)])) && error("output time grid sampling finer than modeling")
	#any([maximum(getfield(TDout[ip],:tgrid).x) > maximum(tgridmod) for ip=1:length(TDout)]) && error("output time > modeling time")

	#! no modeling if source wavelet is zero
	#if(maxval(abs(src_wavelets)) .lt. tiny(rzero_de))
	#        return
	#endif

	# all the propagating wavefields should have same supersources? check that?

	# check dimension of model

	# check if all sources are receivers are inside model
	any(.![(ageom[ip] ∈ model.mgrid) for ip=1:npw]) ? error("sources or receivers not inside model") : nothing

	length(ageom) != npw ? error("ageom size") : nothing
	length(srcwav) != npw ? error("srcwav size") : nothing
	all([issimilar(ageom[ip],srcwav[ip]) for ip=1:npw]) ? nothing : error("ageom and srcwav mismatch")

	# necessary that nss and fields should be same for all nprop
	nss = length(ageom[1]);
	sfields=[names(srcwav[ipw][1].d)[1] for ipw in 1:npw]
	isfields = [Array{Int}(undef,length(sfields[ipw])) for ipw in 1:npw]
	for ipw in 1:npw
		for (i,field) in enumerate(sfields[ipw])
			isfields[ipw][i]=findfirst(isequal(field), [:P,:Vx,:Vz])
		end
	end
	fill(nss, npw) != [length(ageom[ip]) for ip=1:npw] ? error("different supersources") : nothing

	# create acquisition ageometry with each source shooting
	# at every unique receiver position
	irfields = Array{Int}(undef,length(rfields))
	for (i,field) in enumerate(rfields)
		irfields[i]=findfirst(isequal(field), [:P,:Vx,:Vz])
	end

	#(verbose) &&	println(string("\t> number of super sources:\t",nss))

	# find maximum and minimum frequencies in the source wavelets
	freqmin = Utils.findfreq(srcwav[1][1].d[1][:,1],srcwav[1][1].grid,attrib=:min)
	freqmax = Utils.findfreq(srcwav[1][1].d[1][:,1],srcwav[1][1].grid,attrib=:max)

	# minimum and maximum velocities
	vpmin = minimum(broadcast(minimum,[model.bounds[:vp]]))
	vpmax = maximum(broadcast(maximum,[model.bounds[:vp]]))
	#verbose && println("\t> minimum and maximum velocities:\t",vpmin,"\t",vpmax)

	check_fd_stability(vpmin, vpmax, step(model.mgrid[2]), step(model.mgrid[1]), freqmin, freqmax, step(tgridmod), verbose)

	# where to store the boundary values (careful, born scaterrers cannot be inside these!?)
	ibx0=npml-2; ibx1=length(model.mgrid[2])+npml+3
	ibz0=npml-2; ibz1=length(model.mgrid[1])+npml+3
#	ibx0=npml+1; ibx1=length(model.mgrid[2])+npml
#	ibz0=npml+1; ibz1=length(model.mgrid[1])+npml
#	println("**** Boundary Storage Changed **** ")

	# for snaps
	isx0, isz0=npml, npml

	# extend models in the PML layers
	exmodel = Medium_pml_pad_trun(model, nlayer_rand, npml);

	"density values on vx and vz stagerred grids"
	modrr=copy(exmodel[:rhoI])
	modrrvx = get_rhovxI(modrr)
	modrrvz = get_rhovzI(modrr)

	modttI = copy(exmodel[:K])
	modtt = copy(exmodel[:KI])

	if(typeof(attrib_mod)==FdtdBorn)
		(npw≠2) && error("born modeling needs npw=2")
	end

	#create some aliases
	nz, nx = length.(exmodel.mgrid)
	nzd, nxd = length.(model.mgrid)
	δx, δz = step(exmodel.mgrid[2]), step(exmodel.mgrid[1])
	δt = step(tgridmod)
	nt=length(tgridmod)

	δx24I, δz24I = inv(δx * 24.0), inv(δz * 24.0)
	δxI, δzI = inv(δx), inv(δz)
	δt = step(tgridmod)
	δtI = inv(δt)
	nt=length(tgridmod)
	mesh_x, mesh_z = exmodel.mgrid[2], exmodel.mgrid[1]

	# perturbation vectors are required on full space
	δmodtt = zeros(nz, nx)
	δmod = zeros(2*nzd*nxd)
	δmodrr = zeros(nz, nx)
	δmodrrvx = zeros(nz, nx)
	δmodrrvz = zeros(nz, nx)

	# pml_variables
	a_x, b_x, k_xI, a_x_half, b_x_half, k_x_halfI = pml_variables(nx, δt, δx, npml-3, vpmax, vpmin, freqmin, freqmax,
							       [any(abs_trbl .== :left), any(abs_trbl .== :right)])
	a_z, b_z, k_zI, a_z_half, b_z_half, k_z_halfI = pml_variables(nz, δt, δz, npml-3, vpmax, vpmin, freqmin, freqmax,
							       [any(abs_trbl .== :top), any(abs_trbl .== :bottom)])

	gradient=zeros(2*nzd*nxd)
	grad_modtt_stack=SharedMatrix{Float64}(zeros(nz,nx))
	# these are full size, as rhoI -> rhovx and rhovz is also full size
	grad_modrrvx_stack=SharedMatrix{Float64}(zeros(nz,nx))
	grad_modrrvz_stack=SharedMatrix{Float64}(zeros(nz,nx))
	grad_modrr_stack=zeros(nz,nx)
	illum_stack=SharedMatrix{Float64}(zeros(nzd, nxd))

	itsnaps = [argmin(abs.(tgridmod .- tsnaps[i])) for i in 1:length(tsnaps)]

	nrmat=[ageom[ipw][iss].nr for ipw in 1:npw, iss in 1:nss]
	datamat=SharedArray{Float64}(nt,maximum(nrmat),nss)
	data=[Data(tgridmod,ageom[ip],rfields) for ip in 1:length(findall(!iszero, rflags))]

	# default is all prpagating wavefields are active
	activepw=[ipw for ipw in 1:npw]
	pac=PFdtdc(jobname,attrib_mod,npw,activepw,
	    exmodel,model,
	    ageom,srcwav,abs_trbl,sfields,isfields,sflags,
	    irfields,rflags,δt,δxI,δzI,
            nx,nz,nt,δtI,δx24I,δz24I,a_x,b_x,k_xI,a_x_half,b_x_half,k_x_halfI,a_z,b_z,k_zI,a_z_half,b_z_half,k_z_halfI,
	    modtt, modttI,modrr,modrrvx,modrrvz,
	    δmodtt,δmodrr,δmodrrvx,δmodrrvz,
	    δmod,
	    gradient,
	    grad_modtt_stack,
	    grad_modrrvx_stack,
	    grad_modrrvz_stack,grad_modrr_stack,
	    illum_flag,illum_stack,
	    backprop_flag,
	    snaps_flag,
	    itsnaps,
	    gmodel_flag,
	    ibx0,ibz0,ibx1,ibz1,
	    isx0,isz0,
	    datamat,
	    data,
	    verbose)

	# dividing the supersources to workers
	if(nworker===nothing)
		nworker = min(nss, Distributed.nworkers())
	end
	work = Distributed.workers()[1:nworker]
	ssi=[round(Int, s) for s in range(0,stop=nss,length=nworker+1)]
	sschunks=Array{UnitRange{Int64}}(undef, nworker)
	for ib in 1:nworker
		sschunks[ib]=ssi[ib]+1:ssi[ib+1]
	end

	# a distributed array of PFdtdp --- note that the parameters for each super source are efficiently distributed here
	papa=ddata(T=PFdtdp, init=I->PFdtdp(sschunks[I...][1],pac), pids=work);

	return PFdtd(sschunks, papa, pac)
end

"""
Create modeling parameters for each worker.
Each worker performs the modeling of supersources in `sschunks`.
The parameters common to all workers are stored in `pac`.
"""
function PFdtdp(sschunks::UnitRange{Int64},pac::PFdtdc)
	nx=pac.nx; nz=pac.nz; npw=pac.npw

	born_svalue_stack = zeros(nz, nx)

	p=[zeros(nz,nx,3) for ipw in 1:npw];
	dpdx=[zeros(nz,nx,3) for ipw in 1:npw];
	dpdz=[zeros(nz,nx,3) for ipw in 1:npw];
	pp=[zeros(nz,nx,3) for ipw in 1:npw];
	ppp=[zeros(nz,nx,3) for ipw in 1:npw];

	memory_dvx_dx=[zeros(nz,nx) for ipw in 1:npw]
	memory_dvx_dz=[zeros(nz,nx) for ipw in 1:npw]
	memory_dvz_dx=[zeros(nz,nx) for ipw in 1:npw]
	memory_dvz_dz=[zeros(nz,nx) for ipw in 1:npw]
	memory_dp_dx=[zeros(nz,nx) for ipw in 1:npw]
	memory_dp_dz=[zeros(nz,nx) for ipw in 1:npw]

	ss=[PFdtdss(iss, pac) for (issp,iss) in enumerate(sschunks)]

	pap=PFdtdp(ss,p,pp,ppp,dpdx,dpdz,memory_dp_dx,memory_dp_dz,memory_dvx_dx,memory_dvz_dz,
	    born_svalue_stack)

	return pap
end

"""
update the perturbation vector using the perturbed model
in this case, model will be treated as the background model
* `δmod` is [δKI, δrhoI]
"""
function update_δmods!(pac::PFdtdc, δmod::Vector{Float64})
	nx=pac.nx; nz=pac.nz
	nznxd=prod(length.(pac.model.mgrid))
	copyto!(pac.δmod, δmod)
	fill!(pac.δmodtt,0.0)
	δmodtt=view(pac.δmodtt,npml+1:nz-npml,npml+1:nx-npml)
	for i in 1:nznxd
		# put perturbation due to KI as it is
		δmodtt[i] = δmod[i]
	end
	fill!(pac.δmodrr,0.0)
	δmodrr=view(pac.δmodrr,npml+1:nz-npml,npml+1:nx-npml)
	for i in 1:nznxd
		# put perturbation due to rhoI here
		δmodrr[i] = δmod[nznxd+i]
	end
	# project δmodrr onto the vz and vx grids
	get_rhovxI!(pac.δmodrrvx, pac.δmodrr)
	get_rhovzI!(pac.δmodrrvz, pac.δmodrr)
	return nothing
end

"""
This method should be executed only after the updating the main model.
Update the `δmods` when a perturbed `model_pert` is input.
The model through which the waves are propagating
is assumed to be the background model.
"""
function update_δmods!(pac::PFdtdc, model_pert::Medium)
	nznxd=prod(length.(pac.model.mgrid))
	fill!(pac.δmod,0.0)
	copyto!(pac.δmod, model_pert, [:KI, :rhoI])

	for i in 1:nznxd
		pac.δmod[i] -= pac.modtt[i] # subtracting the background model
		pac.δmod[nznxd+i] -= pac.modrr[i] # subtracting the background model
	end
	update_δmods!(pac, pac.δmod)
	return nothing
end

"""
```julia
update!(pa,medium_new)
```
Update `pa` with a new bundle of medium parameters `medium_new`, without additional memory allocation.
This routine is used during inversion, where medium parameters are iteratively updated.
The ability to iteratively run the forward modeling task (with no additional memory allocation) on
various subsurface models is necessary while implementing inversion
algorithms.
"""
function update!(pa::PFdtd, model::Medium)
	return update!(pa.c, model)
end
function update!(pac::PFdtdc, model::Medium)

	copyto!(pac.model, model)

	Medium_pml_pad_trun!(pac.exmodel, pac.model, nlayer_rand)

	copyto!(pac.modttI, pac.exmodel, [:K])
	copyto!(pac.modtt, pac.exmodel, [:KI])
	copyto!(pac.modrr, pac.exmodel, [:rhoI])
	get_rhovxI!(pac.modrrvx, pac.modrr)
	get_rhovzI!(pac.modrrvz, pac.modrr)
	return nothing
end

"""
```julia
update!(pa,srcwav_new,sflags)
```
Update `pa` with a new bundle of source wavelets `srcwav_new`, without additional memory allocation.
Optionally, `sflags` can be changed.
"""
function update!(pa::PFdtd, srcwav::SrcWav, sflags::Any=nothing)
	update_srcwav!(pa,[srcwav],sflags)
end
function update_srcwav!(pa::PFdtd, srcwav::Vector{SrcWav}, sflags=nothing)
	# update srcwav in pa.c
	(length(srcwav) ≠ pa.c.npw) && error("cannot update")
	for i in 1:length(srcwav)
		copyto!(pa.c.srcwav[i], srcwav[i])
	end
	if(!(sflags===nothing))
		copyto!(pa.c.sflags, sflags)
	end

	# fill_wavelets for each supersource
	@sync begin
		for (ip, p) in enumerate(procs(pa.p))
			@async remotecall_wait(p) do
				pap=localpart(pa.p)
				for is in 1:length(pap.ss)
					iss=pap.ss[is].iss
					wavelets=pap.ss[is].wavelets
					fill!.(wavelets, 0.0)
					fill_wavelets!(iss, wavelets, pa.c.srcwav, pa.c.sflags)
				end
			end
		end
	end
end

"""
Create modeling parameters for each supersource.
Every worker models one or more supersources.
"""
function PFdtdss(iss::Int64, pac::PFdtdc)

	irfields=pac.irfields
	isfields=pac.isfields
	npw=pac.npw
	nt=pac.nt
	nx=pac.nx; nz=pac.nz
	nzd,nxd=length.(pac.model.mgrid)
	ageom=pac.ageom
	srcwav=pac.srcwav
	sflags=pac.sflags
	mesh_x, mesh_z = pac.exmodel.mgrid[2], pac.exmodel.mgrid[1]

	# records_output, distributed array among different procs
	records = [zeros(nt,pac.ageom[ipw][iss].nr,length(irfields)) for ipw in 1:npw]

	# gradient outputs
	grad_modtt = zeros(nz, nx)
	grad_modrrvx = zeros(nz, nx)
	grad_modrrvz = zeros(nz, nx)

	# saving illum
	illum =  (pac.illum_flag) ? zeros(nz, nx) : zeros(1,1)

	snaps = (pac.snaps_flag) ? zeros(nzd,nxd,length(pac.itsnaps)) : zeros(1,1,1)

	# source wavelets
	wavelets = [zeros(pac.ageom[ipw][iss].ns,length(isfields[ipw])) for ipw in 1:npw, it in 1:nt]
	fill_wavelets!(iss, wavelets, srcwav, sflags)

	# storing boundary values for back propagation
	nz1, nx1=length.(pac.model.mgrid)
	if(pac.backprop_flag ≠ 0)
		boundary=[zeros(3,nx1+6,nt),
		  zeros(nz1+6,3,nt),
		  zeros(3,nx1+6,nt),
		  zeros(nz1+6,3,nt),
		  zeros(nz1+2*npml,nx1+2*npml,3)]
	else
		boundary=[zeros(1,1,1) for ii in 1:5]
	end
	# source_spray_weights per supersource
	ssprayw = [zeros(4,ageom[ipw][iss].ns) for ipw in 1:npw]
	denomsI = [zeros(ageom[ipw][iss].ns) for ipw in 1:npw]
	isx1=[zeros(Int64,ageom[ipw][iss].ns) for ipw in 1:npw]
	isx2=[zeros(Int64,ageom[ipw][iss].ns) for ipw in 1:npw]
	isz1=[zeros(Int64,ageom[ipw][iss].ns) for ipw in 1:npw]
	isz2=[zeros(Int64,ageom[ipw][iss].ns) for ipw in 1:npw]
	for ipw=1:npw
		for is=1:ageom[ipw][iss].ns
			weights=ssprayw[ipw]
			denom=denomsI[ipw]
			Interpolation.get_spray_weights!(view(weights, :,is), view(denom,is),
				    view(isx1[ipw],is), view(isx2[ipw],is),
				    view(isz1[ipw],is), view(isz2[ipw],is),
				    mesh_x, mesh_z, ageom[ipw][iss].sx[is], ageom[ipw][iss].sz[is])
		end
	end

	# receiver interpolation weights per sequential source
	rinterpolatew = [zeros(4,ageom[ipw][iss].nr) for ipw in 1:npw]
	denomrI = [zeros(ageom[ipw][iss].nr) for ipw in 1:npw]
	irx1=[zeros(Int64,ageom[ipw][iss].nr) for ipw in 1:npw]
	irx2=[zeros(Int64,ageom[ipw][iss].nr) for ipw in 1:npw]
	irz1=[zeros(Int64,ageom[ipw][iss].nr) for ipw in 1:npw]
	irz2=[zeros(Int64,ageom[ipw][iss].nr) for ipw in 1:npw]
	for ipw=1:npw
		for ir=1:ageom[ipw][iss].nr
			weights=rinterpolatew[ipw]
			denom=denomrI[ipw]
			Interpolation.get_interpolate_weights!(view(weights, :,ir), view(denom, ir),
				  view(irx1[ipw],ir), view(irx2[ipw],ir),
				  view(irz1[ipw],ir), view(irz2[ipw],ir),
				  mesh_x, mesh_z, ageom[ipw][iss].rx[ir], ageom[ipw][iss].rz[ir])
		end
	end

	pass=PFdtdss(iss,wavelets,ssprayw,records,rinterpolatew,
	      isx1,isx2,isz1,isz2,irx1,irx2,irz1,irz2,boundary,snaps,illum,
	      grad_modtt,grad_modrrvx,grad_modrrvz)

	return pass
end

include("source.jl")
include("receiver.jl")
include("core.jl")
include("rho_projection.jl")
include("gradient.jl")
include("born.jl")
include("boundary.jl")
include("gallery.jl")

"""
```julia
update!(pa)
```
In-place method to perform the experiment and update `pa` after wave propagation. After update, see
`pa[:data]` and `pa[:snaps]`.
"""
@fastmath function update!(pa::PFdtd)

	global to

	reset_timer!(to)

	@timeit to "initialize!" begin
		# zero out all the results stored in pa.c
		initialize!(pa.c)

		# zero out results stored per worker
		@sync begin
			for (ip, p) in enumerate(procs(pa.p))
				@async remotecall_wait(p) do
					initialize!(localpart(pa.p))
				end
			end
		end
	end

	@timeit to "mod_per_proc!" begin
		# all localparts of DArray are input to this method
		# parallelization over shots
		@sync begin
			for (ip, p) in enumerate(procs(pa.p))
				@async remotecall_wait(p) do
					mod_per_proc!(pa.c, localpart(pa.p))
				end
			end
		end
	end

	@timeit to "stack_grads!" begin
		# stack gradients and illum over sources
		@sync begin
			for (ip, p) in enumerate(procs(pa.p))
				@sync remotecall_wait(p) do
					(pa.c.gmodel_flag) && stack_grads!(pa.c, localpart(pa.p))
					(pa.c.illum_flag) && stack_illums!(pa.c, localpart(pa.p))
				end
			end
		end
	end

	@timeit to "update gradient" begin
		# update gradient model using grad_modtt_stack, grad_modrr_stack
		(pa.c.gmodel_flag) && update_gradient!(pa.c)
	end

	@timeit to "record data" begin
	for ipw in pa.c.activepw
		if(pa.c.rflags[ipw] ≠ 0) # record only if rflags is non-zero
			for ifield in 1:length(pa.c.irfields)

				fill!(pa.c.datamat, 0.0)
				@sync begin
					for (ip, p) in enumerate(procs(pa.p))
						@sync remotecall_wait(p) do
							update_datamat!(ifield, ipw, pa.c, localpart(pa.p))
						end
					end
				end
				update_data!(ifield, ipw, pa.c)
			end
		end
	end
	end
	if(pa.c.verbose)
		show(to)
		println("  ")
	end
	return nothing
end

function update_datamat!(ifield, ipw, pac::PFdtdc, pap::PFdtdp)
	datamat=pac.datamat
        pass=pap.ss
	for issp in 1:length(pass)
		iss=pass[issp].iss
		records=pass[issp].records[ipw]
		for ir in 1:pac.ageom[ipw][iss].nr
			for it in 1:pac.nt
				datamat[it,ir,iss]=records[it,ir,ifield]
			end
		end
        end
end

function update_data!(ifield, ipw, pac::PFdtdc)
	datamat=pac.datamat
	for iss in 1:length(pac.ageom[1])
		data=pac.data[ipw][iss].d[ifield]
		for ir in 1:pac.ageom[ipw][iss].nr
			for it in 1:pac.nt
				data[it,ir]=datamat[it,ir,iss]
			end
		end
	end
end

# update TDout after forming a vector and resampling
#	ipropout=0;
#	for iprop in 1:pac.npw
#		if(pac.rflags[iprop] ≠ 0)
#			ipropout += 1
##			Data.TD_resamp!(pac.data[ipropout], Data.TD_urpos((Array(records[:,:,iprop,:,:])), rfields, tgridmod, ageom[iprop],
##				ageom_urpos[1].nr[1],
##				(ageom_urpos[1].rz[1], ageom_urpos[1].rx[1])
##				))
#		end
#	end
	# return without resampling for testing
	#return [Data.TD(reshape(records[1+(iprop-1)*nd : iprop*nd],tgridmod.nx,recv_n,nss),
	#		       tgridmod, ageom[1]) for iprop in 1:npw]

function stack_illums!(pac::PFdtdc, pap::PFdtdp)
	nx, nz=pac.nx, pac.nz
	illums=pac.illum_stack
	pass=pap.ss
	for issp in 1:length(pass)
		gs=pass[issp].illum
		gss=view(gs,npml+1:nz-npml,npml+1:nx-npml)
		@. illums += gss
	end
end

# modelling for each processor
function mod_per_proc!(pac::PFdtdc, pap::PFdtdp)
	# source_loop
	for issp in 1:length(pap.ss)
		reset_per_ss!(pap)

		iss=pap.ss[issp].iss

		if(pac.backprop_flag==-1)
			"initial conditions from boundary for first propagating field only"
			boundary_force_snap_p!(issp,pac,pap.ss,pap)
			boundary_force_snap_vxvz!(issp,pac,pap.ss,pap)
		end

		prog = Progress(pac.nt, dt=1, desc="\tmodeling supershot $iss/$(length(pac.ageom[1])) ",
		  		color=:white)
		# time_loop
		"""
		* don't use shared arrays inside this time loop, for speed when using multiple procs
		"""
		for it=1:pac.nt

			pac.verbose && next!(prog, :white)

			advance!(pac,pap)

			# force p[1] on boundaries
			(pac.backprop_flag==-1) && boundary_force!(it,issp,pac,pap.ss,pap)

			add_source!(it, issp, iss, pac, pap.ss, pap, Source_B1())

			# no born flag for adjoint modelling
			if(!pac.gmodel_flag)
				(typeof(pac.attrib_mod)==FdtdBorn) && add_born_sources!(issp, pac, pap.ss, pap)
			end

			# record boundaries after time reversal already
			(pac.backprop_flag==1) && boundary_save!(pac.nt-it+1,issp,pac,pap.ss,pap)

			record!(it, issp, iss, pac, pap.ss, pap, Receiver_B1())

			(pac.gmodel_flag) && compute_gradient!(issp, pac, pap.ss, pap)

			(pac.illum_flag) && compute_illum!(issp, pap.ss, pap)

			if(pac.snaps_flag)
				iitsnaps=findall(x->==(x,it),pac.itsnaps)
				for itsnap in iitsnaps
					snaps_save!(itsnap,issp,pac,pap.ss,pap)
				end
			end

		end # time_loop
		"now pressure is at [nt], velocities are at [nt-1/2]"

		"one more propagating step to save pressure at [nt+1] -- for time revarsal"
		advance!(pac,pap)

		"save last snap of pressure field"
		(pac.backprop_flag==1) && boundary_save_snap_p!(issp,pac,pap.ss,pap)

		"one more propagating step to save velocities at [nt+3/2] -- for time reversal"
		advance!(pac,pap)

		"save last snap of velocity fields with opposite sign for adjoint propagation"
		(pac.backprop_flag==1) && boundary_save_snap_vxvz!(issp,pac,pap.ss,pap)

		"scale gradients for each issp"
		(pac.gmodel_flag) && scale_gradient!(issp, pap.ss, step(pac.model.mgrid[2])*step(pac.model.mgrid[1]))

	end # source_loop
end # mod_per_shot

# Need illumination to estimate the approximate diagonal of Hessian
@inbounds @fastmath function compute_illum!(issp::Int64, pass::Vector{PFdtdss}, pap::PFdtdp)
	# saving illumination to be used as preconditioner
	p=pap.p
	illum=pass[issp].illum
	pp=view(p[1],:,:,1)
	for i in eachindex(illum)
		illum[i] += abs2(pp[i])
	end
end

@fastmath @inbounds function snaps_save!(itsnap::Int64,issp::Int64,pac::PFdtdc,pass::Vector{PFdtdss},pap::PFdtdp)
	isx0=pac.isx0
	isz0=pac.isz0
	p=pap.p
	snaps=pass[issp].snaps
	for ix=1:size(snaps,2)
		@simd for iz=1:size(snaps,1)
			snaps[iz,ix,itsnap]=p[1][isz0+iz,isx0+ix,1]
		end
	end
end

function check_fd_stability(vpmin::Float64, vpmax::Float64, δx::Float64, δz::Float64, freqmin::Float64, freqmax::Float64, δt::Float64, verbose::Bool)

	# check spatial sampling
	δs_temp=round(vpmin/5.0/freqmax,digits=2);
	δs_max = maximum([δx, δz])
	str1=@sprintf("%0.2e",δs_max)
	str2=@sprintf("%0.2e",δs_temp)
	if(str1 ≠ str2)
		if(all(δs_max .> δs_temp))
			@warn "decrease spatial sampling ($str1) below $str2"
		else
			@info "spatial sampling ($str1) can be as high as $str2"
		end
	end

	# check time sampling
	δs_min = minimum([δx, δz])
	δt_temp=0.5*δs_min/vpmax
	str1=@sprintf("%0.2e", δt)
	str2=@sprintf("%0.2e", δt_temp)
	if(str1 ≠ str2)
		if(all(δt .> δt_temp))
			@warn "decrease time sampling ($str1) below $str2"
		else
			@info "time sampling ($str1) can be as high as $str2"
		end
	end

end

"""
Output vectors related to PML boundaries.
"""
function pml_variables(
		nx::Int64,
		δt::Float64,
		δx::Float64, mesh_na_pml::Int64,
		velmax::Float64, velmin::Float64,
		freqmin::Float64, freqmax::Float64,
	        flags::Vector{Bool}
		)

	NPOWER = 2.e0
	K_MAX_PML = 1.e0 # from Gedney page 8.11
	ALPHA_MAX_PML = 2.e0 * pi * ((freqmin + freqmax))/4.e0 # from Festa and Vilotte

	# thickness of the PML layer in meters
	thickness_PML_x = mesh_na_pml * δx

	"reflection coefficient (INRIA report section 6.1) http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf"
	Rcoef = 0.001e0

	#! check that NPOWER is okay
	#if(NPOWER < 1) stop "NPOWER must be greater than 1"

	# compute d0 from INRIA report section 6.1 http://hal.inria.fr/docs/00/07/32/19/PDF/RR-3471.pdf
	d0_x = - (NPOWER + 1) * (velmax+velmin)/2.0 * log(Rcoef) / (2.e0 * thickness_PML_x)

	d_x=zeros(nx);	d_x_half=zeros(nx)
	k_x=ones(nx);	k_xI=ones(nx)
	k_x_half=ones(nx);	k_x_halfI=ones(nx)
	alpha_x=zeros(nx);	alpha_x_half=zeros(nx)
	a_x=zeros(nx);	a_x_half=zeros(nx)
	b_x=zeros(nx);	b_x_half=zeros(nx)

	"damping in the X direction"
	"origin of the PML layer (position of right edge minus thickness, in meters)"
	xoriginleft = thickness_PML_x
	xoriginright = (nx-1)*δx - thickness_PML_x

	for ix=1:nx
		# abscissa of current grid point along the damping profile
		xval = δx * real(ix-1)

		#---------- left edge
		if(flags[1])

			# define damping profile at the grid points
			abscissa_in_PML = xoriginleft - xval
			if(abscissa_in_PML >= 0.0)
				abscissa_normalized = abscissa_in_PML / thickness_PML_x
				d_x[ix] = d0_x * abscissa_normalized.^NPOWER
				# this taken from Gedney page 8.2
				k_x[ix] = 1.e0 + (K_MAX_PML - 1.e0) * abscissa_normalized.^NPOWER
				alpha_x[ix] = ALPHA_MAX_PML * (1.e0 - abscissa_normalized) + 0.1e0 * ALPHA_MAX_PML
			end

			# define damping profile at half the grid points
			abscissa_in_PML = xoriginleft - (xval + δx/2.e0)
			if(abscissa_in_PML >= 0.0)
				abscissa_normalized = abscissa_in_PML / thickness_PML_x
				d_x_half[ix] = d0_x * abscissa_normalized.^NPOWER
				# this taken from Gedney page 8.2
				k_x_half[ix] = 1.e0 + (K_MAX_PML - 1.e0) * abscissa_normalized.^NPOWER
				alpha_x_half[ix] = ALPHA_MAX_PML * (1.e0 - abscissa_normalized) + 0.1e0 * ALPHA_MAX_PML
			end
		end

		#---------- right edge
		if(flags[2])

			# define damping profile at the grid points
			abscissa_in_PML = xval - xoriginright
			if(abscissa_in_PML >= 0.0)
				abscissa_normalized = abscissa_in_PML / thickness_PML_x
				d_x[ix] = d0_x * abscissa_normalized^NPOWER
				# this taken from Gedney page 8.2
				k_x[ix] = 1.e0 + (K_MAX_PML - 1.e0) * abscissa_normalized^NPOWER
				alpha_x[ix] = ALPHA_MAX_PML * (1.e0 - abscissa_normalized) + 0.1e0 * ALPHA_MAX_PML
			end

			# define damping profile at half the grid points
			abscissa_in_PML = xval + δx/2.e0 - xoriginright
			if(abscissa_in_PML >= 0.0)
				abscissa_normalized = abscissa_in_PML / thickness_PML_x
				d_x_half[ix] = d0_x * abscissa_normalized^NPOWER
				# this taken from Gedney page 8.2
				k_x_half[ix] = 1.e0 + (K_MAX_PML - 1.e0) * abscissa_normalized^NPOWER
				alpha_x_half[ix] = ALPHA_MAX_PML * (1.e0 - abscissa_normalized) + 0.1e0 * ALPHA_MAX_PML
			end

		end
#
		# just in case, for -5 at the end
		(alpha_x[ix] < 0.0) ? alpha_x[ix] = 0.0 : nothing
		(alpha_x_half[ix] < 0.0) ? alpha_x_half[ix] = 0.0 : nothing

		b_x[ix] = exp(- (d_x[ix] / k_x[ix] + alpha_x[ix]) * δt)
		b_x_half[ix] = exp(- (d_x_half[ix] / k_x_half[ix] + alpha_x_half[ix]) * δt)

		# this to avoid division by zero outside the PML
		(abs(d_x[ix]) > 1.e-6) ? a_x[ix] = d_x[ix] * (b_x[ix] - 1.e0) / (k_x[ix] * (d_x[ix] + k_x[ix] * alpha_x[ix])) : nothing
		(abs(d_x_half[ix]) > 1.e-6) ? a_x_half[ix] = d_x_half[ix] * (b_x_half[ix] - 1.e0) / (k_x_half[ix] * (d_x_half[ix] + k_x_half[ix] * alpha_x_half[ix])) : nothing

	end
	k_xI = k_x.^(-1.e0)
	k_x_halfI = k_x_half.^(-1.e0)

	return a_x, b_x, k_xI, a_x_half, b_x_half, k_x_halfI
end

include("getprop.jl")


text:
b(y) = a(y) +3
a(x) = x+1

@show b(2)

ext: jl

max_stars_repo_name: GabrielTrettel/Portran

max_issues_repo_name: GabrielTrettel/Portran

max_forks_repo_name: GabrielTrettel/Portran


text:
# # Vertical Collision Avoidance (VCAS)
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/models/VertCAS.ipynb)
#
# This example [^JK] [^ABKL] considers commercial aircraft, which are required to
# operate with a collision avoidance system that gives vertical climbrate advisories
# to pilots to prevent near midair collisions (NMACs). An NMAC occurs when the
# aircraft are separated by less than 100 ft vertically and 500 ft horizontally.

module VertCAS  #jl

using NeuralNetworkAnalysis, LinearAlgebra

# ## Model
#
# This benchmark is a closed-loop variant of aircraft collision avoidance
# system ACAS X. The scenario involves two aircraft, the ownship and the
# intruder, where the ownship is equipped with a collision avoidance system
# referred to as VerticalCAS [3]. VerticalCAS once every second issues vertical
# climbrate advisories to the ownship pilot to avoid a near mid-air collision
# (NMAC), a region where the ownship and intruder are separated by less than
# 100ft vertically and 500ft horizontally. The ownship (black) is assumed
# to have a constant horizontal speed, and the intruder (red) is assumed to
# follow a constant horizontal trajectory towards ownship, see the figure below.
#
# ![](https://raw.githubusercontent.com/JuliaReach/NeuralNetworkAnalysis.jl/master/models/VertCAS/nmac.png)
#
# The current geometry of the system is described by:
#
# 1) Intruder altitude relative to ownship $h(ft)$
# 2) Ownship vertical climbrate $\dot h_0 (ft/min)$
# 3) The seconds until the ownship (black) and intruder (red) are no longer horizontally separated $τ(s)$
#
# We can, therefore, assume that the intruder is static and the horizontal
# separation $\tau$ decreases by one each second.
# There are 9 advisories and each of them instructs the pilot to accelerate
# until the vertical climbrate of the ownship complies with the advisory:
#
# 1) COC: Clear of Conflict
# 2) DNC: Do Not Climb
# 3) DND: Do Not Descend
# 4) DES1500: Descend at least 1500 ft/min
# 5) CL1500: Climb at least 1500 ft/min
# 6) SDES1500: Strengthen Descent to at least 1500 ft/min
# 7) SCL1500: Strengthen Climb to at least 1500 ft/min
# 8) SDES2500: Strengthen Descent to at least 2500 ft/min
# 9) SCL2500: Strengthen Climb to at least 2500 ft/min
#
# In addition to the parameters describing the geometry of the encounter, the
# current state of the system stores the advisory $adv$ issued to the ownship
# at the previous time step. VerticalCAS is implemented as nine ReLU networks
# $N_i$, one for each (previous) advisory, with three inputs
# $(h,\dot{h}_0,\tau)$, five fully-connected hidden layers of 20 units each,
# and nine outputs representing the score of each possible advisory. Therefore,
# given a current state $(h,\dot{h}_0,\tau,\text{adv})$, the new advisory
# $adv$ is obtained by computing the argmax of the output of $N_{\text{adv}}$
# on $(h,\dot{h}_0,\tau)$.
# Given the new advisory, if the current climbrate does not comply with it, the
# pilot can choose acceleration $\ddot{h}_0$ from the given set:
#
# 1) COC: $\{-\frac{g}{8}, 0, \frac{g}{8}\}$
# 2) DNC: $\{-\frac{g}{3}, -\frac{7g}{24}, -\frac{g}{4}\}$
# 3) DND: $\{\frac{g}{4}, \frac{7g}{24}, \frac{g}{3}\}$
# 4) DES1500: $\{-\frac{g}{3}, -\frac{7g}{24}, -\frac{g}{4}\}$
# 5) CL1500: $\{\frac{g}{4}, \frac{7g}{24}, \frac{g}{3}\}$
# 6) SDES1500: $\{-\frac{g}{3}\}$
# 7) SCL1500: $\{\frac{g}{3}\}$
# 8) SDES2500: $\{-\frac{g}{3}\}$
# 9) SCL2500: $\{\frac{g}{3}\}$
#
# where $g$ represents the gravitational constant $32.2 \ \text{ft/s}^2$.
# If the new advisory is COC(1), then it can be any acceleration from the set ${−g/8, 0, g/8}$.
# For all remaining advisories, if the previous advisory coincides with the new one and the
# current climb rate complies with the new advisory (e.g., $\dot{h}_0$ is non-positive for DNC and
# $\dot{h}_0 ≥ 1500$ for CL1500) the acceleration `\ddot{h}_0$ is $0$.
#
# Given the current system state $(h,\dot{h}_0,\tau,\text{adv})$, the new
# advisory $\text{adv}'$ and the acceleration $\ddot{h}_0$, the new state
# of the system $(h(t+1),\dot{h}_0(t+1),\tau(t+1),\text{adv}(t+1))$ can
# be computed as follows:
#
# ```math
# \begin{aligned}
# h(t+1) &=& h - \dot{h}_0 \Delta\tau - 0.5 \ddot{h}_0 \Delta\tau^2 \\
# \dot{h}_0(t+1) &=& \dot{h}_0 + \ddot{h}_0 \Delta\tau \\
# \tau(t+1) &=& \tau - \Delta\tau \\
# \text{adv}(t+1) &=& \text{adv}'
# \end{aligned}
# ```
# where $\Delta\tau=1$.

const g = 32.2  # gravitational constant

## accelerations (middle)
const ACC_MIDDLE = Dict(:COC => 0.0, :DNC => -7g/24, :DND => 7g/24,
                        :DES1500 => -7g/24, :CL1500 => 7g/24, :SDES1500 => -g/3,
                        :SCL1500 => g/3, :SDES2500 => -g/3, :SCL2500 => g/3)

## continuous dynamics matrix (h, hdot0)
const Δτ = 1.0
const A = [1  -Δτ; 0  1.];

# We load the controllers in a dictionary with the keys being the advisories.

CONTROLLERS = Dict{Symbol, Any}()

CTRL_IDX = [:COC, :DNC, :DND, :DES1500,
            :CL1500, :SDES1500, :SCL1500,
            :SDES2500, :SCL2500]

path = @modelpath("VertCAS", "")
for i = 1:9
    file = joinpath(path, "VertCAS_noResp_pra0$(i)_v9_20HU_200.nnet")
    adv = CTRL_IDX[i]
    CONTROLLERS[adv] = read_nnet(file)
end

struct State{T, N}
    state::T # state for (h, hdot0) variables
    τ::N
    adv::Symbol
end

ADVISORIES = Dict{Symbol, Any}()

## for every advisory, return a function to check
## whether the current climbrate complies with the advisory or not
## if the advisory *does not* comply => it is changed according to the ACC dictionary
## otherwise, the new acceleration is zero
using Symbolics
@variables x

ADVISORIES[:COC] = EmptySet(1)
ADVISORIES[:DNC] = HalfSpace(x <= 0)
ADVISORIES[:DND] = HalfSpace(x >= 0)
ADVISORIES[:DES1500] = HalfSpace(x <= -1500)
ADVISORIES[:CL1500] = HalfSpace(x >= 1500)
ADVISORIES[:SDES1500] = HalfSpace(x <= -1500)
ADVISORIES[:SCL1500] = HalfSpace(x >= 1500)
ADVISORIES[:SDES2500] = HalfSpace(x <= -2500)
ADVISORIES[:SCL2500] = HalfSpace(x >= 2500)

## this function receives X = [h, hdot0, τ, adv′] and the
## *previous* advisory adv
function get_acceleration(X::State, adv; ACC=ACC_MIDDLE)

    ## obtain projection on hdot
    hdot = _interval(X.state, 2)

    ## transform units from ft/s to ft/min
    hdot = 60 * hdot

    ## new advisory
    adv′ = X.adv

    ## check whether the current state complies with the advisory
    comply = hdot ⊆ ADVISORIES[adv′]

    if adv == adv′ && comply
        return 0.0
    else
        return ACC[adv′]
    end
end

const normalization_additive = -[0.0, 0, 20]
const normalization_multiplicative = 1.0 ./ [16000.0, 5000, 40]

function normalize(x::AbstractVector)
    y = x .+ normalization_additive
    z = y .* normalization_multiplicative
    return z
end

function normalize(X::AbstractZonotope)
    Y = minkowski_sum(X, Singleton(normalization_additive))
    Z = linear_map(diagm(normalization_multiplicative), Y)
    return Z
end

## scalar case; alg is ignored
function forward_adv(X::Singleton, τ, adv; alg=nothing)
    v = vcat(element(X), τ)
    v = normalize(v)
    u = forward(CONTROLLERS[adv], v)
    imax = argmax(u)
    return CTRL_IDX[imax]
end

## set-based case
function forward_adv(X::AbstractZonotope, τ, adv; alg=Ai2())
    Y = cartesian_product(X, Singleton([τ]))

    Y = normalize(Y)
    out = forward_network(alg, CONTROLLERS[adv], Y)

    imax = argmax(high(out))
    return CTRL_IDX[imax]
end

function VCAS!(out::Vector{State{T, N}}, KMAX; ACC=ACC_MIDDLE, alg_nn=Ai2()) where {T, N}

    ## unpack initial state
    X0 = first(out)
    S = X0.state
    τ = X0.τ
    adv = X0.adv

    ## get initial acceleration
    hddot = ACC[adv]

    for i in 1:KMAX
        ## compute next state
        b = [-hddot*Δτ^2 / 2, hddot * Δτ]
        S′ = affine_map(A, S, b)
        τ′ = τ - 1
        adv′ = forward_adv(S′, τ′, adv, alg=alg_nn)

        ## store new state
        X′ = State(S′, τ′, adv′)
        push!(out, X′)

        ## get acceleration from network
        ## this logic only works for ACC_MIDDLE
        hddot = get_acceleration(X′, adv; ACC=ACC)

        ## update current state
        S = S′
        τ = τ′
        adv = adv′
    end
    return out
end;

# ## Specifications
#
# For this benchmark the aim is to verify that the ownship avoids entering the
# NMAC zone after $k \in \{1, \dots, 10\}$ time steps, i.e., $h(k) > 100$ or
# $h(k) < -100$, for all possible choices of acceleration by the pilot. The
# set of initial states considered is as follows: $h(0) \in [-133, -129]$,
# $\dot{h}_0(0) \in \{-19.5, -22.5, -25.5, -28.5\}$, $\tau(0) = 25$ and
# $\text{adv}(0) = \text{COC}$.

bad_states = HalfSpace([1.0, 0.0], 100.) ∩ HalfSpace([-1.0, 0.0], 100.)

## property for guaranteed violation
predicate = X -> X ⊆ bad_states
predicate_sol = sol -> any(predicate(R) for F in sol for R in F);

# ## Results

# ### Simulation

const h0 = Interval(-133, -129)
const hdot0 = [-19.5,-22.5, -25.5, -28.5]
const τ0 = 25.0
const adv0 = :COC

function _random_states(k=1, include_vertices::Bool=false, rand_h0::Bool=true)
    N = Float64
    T = Singleton{N, Vector{N}}
    states = Vector{State{T, N}}()
    xs = sample(h0, k, include_vertices=include_vertices)
    for x in xs
        if rand_h0
            ## use a random value for y
            y = hdot0[rand(1:4)]
            S0 = State(Singleton([x[1], y]), τ0, adv0)
            push!(states, S0)
            continue
        end
        ## use all possible values for y
        for i in 1:4
            y = hdot0[i]
            S0 = State(Singleton([x[1], y]), τ0, adv0)
            push!(states, S0)
        end
    end
    return states
end

function _all_states()
    S0 = [convert(Zonotope, concretize(h0 × Singleton([hdot0[i]]))) for i in 1:4]
    return [State(S0i, τ0, adv0) for S0i in S0]
end

function simulate_VCAS(X0::State; KMAX=10)
    out = [X0]
    sizehint!(out, KMAX+1)

    VCAS!(out, KMAX, ACC=ACC_MIDDLE)
    return out
end

## project onto the h variable
function _project(X::Vector{State{T, N}}) where {T<:Singleton, N}
    return [Singleton([Xi.state.element[1], Xi.τ]) for Xi in X]
end

_interval(X::LazySet, i) = overapproximate(Projection(X, (i,)), Interval)

function _project(X::Vector{State{T, N}}) where {T<:Zonotope, N}
    Xint = [_interval(Xi.state, 1) × Singleton([Xi.τ]) for Xi in X]
end

function run(X0)
    ensemble = [simulate_VCAS(X0i) for X0i in X0]
    res = _project.(ensemble)
    return res
end

function check(sol)
    println("property checking")
    res_pred = @timed predicate_sol(sol)
    print_timed(res_pred)
    if res_pred.value
        println("The property is violated.")
    else
        println("The property may be satisfied.")
    end
end;

# Simulation result for random choice of velocity:

X0 = _random_states(10, true, false)  # randomly sampled points (incl. vertices)
println("$(length(X0)) simulations with central advisories")
res = @timed begin
    res1 = run(X0)
    check(res1)
end
println("total analysis time")
print_timed(res);

# Simulation result for all choices of velocity:

println("flowpipe construction (unsound) with central advisories")
res = @timed begin
    res2 = run(_all_states())
    check(res1)
end
println("total analysis time")
print_timed(res);

# Finally we plot the results:

using Plots
import DisplayAs

function plot_helper()
    fig = plot(xlims=(-200, 200), ylims=(14, 26), xlab="h (vertical distance)",
               ylab="τ (time to reach horizontally)")
    plot!(fig, bad_states, alpha=0.2, c=:red, lab="unsafe states")
    return fig
end

fig = plot_helper()
for o in res1
    plot!(fig, o, alpha=1.0)
end
fig = DisplayAs.Text(DisplayAs.PNG(fig))
## savefig("VertCAS-rand.png")
fig

#-

fig = plot_helper()
for (i, c) in [(1, :blue), (2, :orange), (3, :green), (4, :cyan)]
    [plot!(fig, o, lw=2.0, alpha=1., markershape=:none, seriestype=:shape, c=c) for o in res2[i]]
end
fig = DisplayAs.Text(DisplayAs.PNG(fig))
## savefig("VertCAS-sets.png")
fig

# ## References

# [^JK]: Julian, K. D., & Kochenderfer, M. J. (2019). *A reachability method for
#        verifying dynamical systems with deep neural network controllers.*
#        [arXiv preprint arXiv:1903.00520.](https://arxiv.org/pdf/1903.00520.pdf).
#
# [^ABKL]: Akintunde, M. E., Botoeva, E., Kouvaros, P., & Lomuscio, A. (2020).
#          *Formal Verification of Neural Agents in Non-deterministic Environments.*
#          In [Proceedings of the 19th International Conference on Autonomous Agents and
#          Multiagent Systems (pp. 25-33)](http://ifaamas.org/Proceedings/aamas2020/pdfs/p25.pdf).

end  #jl
nothing  #jl


text:
module CartesianGeometry

# Dependencies of this module

using StaticArrays: SVector, MVector, @SVector
using Gridap
using Gridap.Helpers

# Functionality provided

export CartesianGrid
export CartesianDiscreteModel
import Base: size, getindex, IndexStyle
import Gridap: points, cells, celltypes, cellorders
import Gridap: Grid, FaceLabels
import Gridap: UnstructuredGrid
import Gridap: FlexibleUnstructuredGrid
import Gridap: GridGraph
import Gridap: FullGridGraph

struct CartesianGrid{D} <: Grid{D,D}
  dim_to_limits::NTuple{D,NTuple{2,Float64}}
  dim_to_ncells::NTuple{D,Int}
  extrusion::NTuple{D,Int}
  order:: Int
  map
end

const identity = (x -> x)

function CartesianGrid(;partition::NTuple{D,Int},domain=nothing,
                        order::Int=1,map=identity) where D
  _cartesiangrid(partition,domain,order,map)
end

function points(self::CartesianGrid)
  dim_to_npoint = tuple([ i+1 for i in self.dim_to_ncells ]...)
  CartesianGridPoints(self.map,self.dim_to_limits,dim_to_npoint)
end

cells(self::CartesianGrid) = CartesianGridCells(self.dim_to_ncells)

celltypes(self::CartesianGrid) = ConstantCellValue(self.extrusion,prod(self.dim_to_ncells))

cellorders(self::CartesianGrid) = ConstantCellValue(self.order,prod(self.dim_to_ncells))

"""
Construct an `UnstructuredGrid` from a `CartesianGrid`
"""
function UnstructuredGrid(grid::CartesianGrid{D}) where D
  ps = _compute_points(grid)
  ts = celltypes(grid)
  os = cellorders(grid)
  data, ptrs = _compute_cells(grid,UnstructuredGrid{D,D})
  UnstructuredGrid(ps,data,ptrs,ts,os)
end

"""
Construct an `FlexibleUnstructuredGrid` from a `CartesianGrid`
"""
function FlexibleUnstructuredGrid(grid::CartesianGrid{D}) where D
  ps = _compute_points(grid)
  ts = celltypes(grid)
  os = cellorders(grid)
  cs = _compute_cells(grid,FlexibleUnstructuredGrid{D,D})
  FlexibleUnstructuredGrid(ps,cs,ts,os)
end

"""
DiscreteModel associated with a CartesianGrid
"""
struct CartesianDiscreteModel{D} <: DiscreteModel{D}
  cgrid::CartesianGrid{D}
  ugrid::UnstructuredGrid
  gridgraph::FullGridGraph
end

function CartesianDiscreteModel(; args...)
  cgrid = CartesianGrid(; args...)
  ugrid = UnstructuredGrid(cgrid)
  gridgraph = FullGridGraph(ugrid)
  CartesianDiscreteModel(cgrid, ugrid, gridgraph)
end

Grid(model::CartesianDiscreteModel{D},::Val{D}) where D = model.cgrid

function Grid(model::CartesianDiscreteModel{D},::Val{Z}) where {D,Z}

  ugrid = model.ugrid
  gridgraph = model.gridgraph
  face_to_vertices = connections(gridgraph,Z,0)
  nfaces = length(face_to_vertices)
  fcode = tuple([HEX_AXIS for i in 1:Z]...)

  order = 1
  @notimplementedif order != (cellorders(model.cgrid)).value

  _points = points(ugrid)
  _cells_data = face_to_vertices.data
  _cells_ptrs = face_to_vertices.ptrs
  _ctypes = ConstantCellValue(fcode,nfaces)
  _corders = ConstantCellValue(order,nfaces)

  UnstructuredGrid(
    _points, _cells_data, _cells_ptrs, _ctypes, _corders)
end

GridGraph(model::CartesianDiscreteModel) = model.gridgraph

FullGridGraph(model::CartesianDiscreteModel) = model.gridgraph

#@fverdugo precompute this result
function FaceLabels(model::CartesianDiscreteModel{D}) where D
  dim_to_face_to_geolabel = Vector{Vector{Int}}(undef,D+1)
  dim_to_offset = _generate_dim_to_offset(D)
  interior_id = dim_to_offset[end]+1
  boundary_id = -1
  for d in 0:(D-1)
    face_to_cells = connections(model.gridgraph,d,D)
    cell_to_faces = connections(model.gridgraph,D,d)
    offset = dim_to_offset[d+1]
    face_to_geolabel = _generate_pre_geolabel(
      face_to_cells,
      cell_to_faces,
      offset,
      interior_id,boundary_id,d,D)
    dim_to_face_to_geolabel[d+1] = face_to_geolabel
  end
  for d in 0:(D-2)
    for j in (d+1):(D-1)
      dface_to_jfaces = connections(model.gridgraph,d,j)
      dface_to_geolabel = dim_to_face_to_geolabel[d+1]
      jface_to_geolabel = dim_to_face_to_geolabel[j+1]
      _fix_dface_geolabels!(
        dface_to_geolabel,
        jface_to_geolabel,
        dface_to_jfaces,
        interior_id,boundary_id)
    end
  end
  _ncells = ncells(model.cgrid)
  dim_to_face_to_geolabel[end] = ConstantCellValue(interior_id,_ncells)
  phys_labels = [ [i] for i in 1:interior_id ]
  push!(phys_labels,[i for i in 1:(interior_id-1)])
  tag_to_name = ["physical_tag_$i" for i in 1:interior_id-1]
  push!(tag_to_name,"interior")
  push!(tag_to_name,"boundary")
  FaceLabels(dim_to_face_to_geolabel, phys_labels, tag_to_name)
end

# Helpers

function _generate_dim_to_offset(D)
  code = tuple(fill(HEX_AXIS,D)...)
  polytope = Polytope(code)
  dim_to_offset = zeros(Int,D+1)
  for d in 1:D
    dim_to_offset[d+1] = dim_to_offset[d] + num_nfaces(polytope,d-1)
  end
  dim_to_offset
end

function _generate_pre_geolabel(
  face_to_cells,
  cell_to_faces,
  offset,
  interior_id,
  boundary_id,d,D)

  nfaces = length(face_to_cells)
  face_to_geolabel = fill(interior_id,nfaces)

  max_ncells_around = 2^(D-d)

  _generate_pre_geolabel_kernel!(
    face_to_geolabel,
    face_to_cells.data,
    face_to_cells.ptrs,
    cell_to_faces.data,
    cell_to_faces.ptrs,
    offset,boundary_id,max_ncells_around)

  face_to_geolabel
end

function _generate_pre_geolabel_kernel!(
  face_to_geolabel,
  face_to_cells_data,
  face_to_cells_ptrs,
  cell_to_faces_data,
  cell_to_faces_ptrs,
  offset,boundary_id,max_ncells_around)

  nfaces = length(face_to_geolabel)
  for face in 1:nfaces
    a = face_to_cells_ptrs[face]-1
    ncells_around = face_to_cells_ptrs[face+1] - (a+1)
    if ncells_around == 1
      icell_around = 1
      cell = face_to_cells_data[a+icell_around]
      b = cell_to_faces_ptrs[cell]-1
      nlfaces = cell_to_faces_ptrs[cell+1] - (b+1)
      for lface in 1:nlfaces
        face2 = cell_to_faces_data[b+lface]
        if face == face2
          face_to_geolabel[face] = lface + offset
          break
        end
      end
    elseif ncells_around != max_ncells_around
      face_to_geolabel[face] = boundary_id
    end
  end

end

function _fix_dface_geolabels!(
  dface_to_geolabel,
  jface_to_geolabel,
  dface_to_jfaces,
  interior_id,boundary_id)

  for (dface, jfaces) in enumerate(dface_to_jfaces)
    if dface_to_geolabel[dface] != boundary_id
      continue
    end
    for jface in jfaces
      geolabel = jface_to_geolabel[jface]
      if geolabel != interior_id && geolabel != boundary_id
        dface_to_geolabel[dface] = geolabel
        break
      end
    end
  end

end

function _cartesiangrid(partition::NTuple{D,Int},domain,order,map) where D
  if domain === nothing
    _domain = [ i*j for i in ones(D) for j in 0:1 ]
    # _domain = [ i*(-1)^j for i in ones(D) for j in 1:2 ]
  else
    _domain = domain
  end
  dim_to_limits = tuple([(_domain[2*i-1],_domain[2*i]) for i in 1:D ]...)
  extrusion = tuple(fill(HEX_AXIS,D)...)
  dim_to_ncells = partition
  @notimplementedif order != 1
  CartesianGrid{D}(dim_to_limits,dim_to_ncells,extrusion,order,map)
end

struct CartesianGridPoints{D} <: IndexCellValue{Point{D,Float64},D}
  map
  dim_to_limits::NTuple{D,NTuple{2,Float64}}
  dim_to_npoint::NTuple{D,Int}
end

size(self::CartesianGridPoints) = self.dim_to_npoint

IndexStyle(::Type{CartesianGridPoints{D}} where D) = IndexCartesian()

function getindex(self::CartesianGridPoints{D}, I::Vararg{Int, D}) where D
  p = zero(MVector{D,Float64})
  @inbounds for d in 1:D
    xa = self.dim_to_limits[d][1]
    xb = self.dim_to_limits[d][2]
    p[d] =  xa + (I[d]-1)*(xb-xa)/(self.dim_to_npoint[d]-1)
  end
  # map(Point(p)
  # self.map != identity ? self.map(Point(p)) : Point(p)
  self.map(Point(p))
end

struct CartesianGridCells{D,L} <: IndexCellArray{Int,1,SVector{L,Int},D}
  dim_to_ncell::SVector{D,Int}
end

function CartesianGridCells(dim_to_ncell::NTuple{D,Int}) where D
  CartesianGridCells{D,2^D}(dim_to_ncell)
end

size(self::CartesianGridCells) = self.dim_to_ncell.data

IndexStyle(::Type{CartesianGridCells{D,L}} where {D,L}) = IndexCartesian()

function getindex(self::CartesianGridCells{D,L}, I::Vararg{Int, D}) where {D,L}
  dim_to_ngpoint = 1 .+ self.dim_to_ncell
  dim_to_nlpoint = @SVector fill(2,D)
  offset = @SVector fill(1,D)
  pointgids = LinearIndices(dim_to_ngpoint.data)
  cellpointlids = CartesianIndices(dim_to_nlpoint.data)
  cellgid = CartesianIndex(I...) - CartesianIndex(offset.data)
  cellpointgids = cellpointlids .+ cellgid
  ids = zero(MVector{L,Int})
  @inbounds for (l,pgid) in enumerate(cellpointgids)
    ids[l] = pointgids[pgid]
  end
  SVector{L,Int}(ids)
end

function _compute_points(grid::CartesianGrid{D}) where D
  ps = Array{Point{D,Float64},1}(undef,(length(points(grid)),))
  for (i,xi) in enumerate(points(grid))
    ps[i] = xi
  end
  ps
end

function _compute_cells(grid::CartesianGrid{D},::Type{FlexibleUnstructuredGrid{D,D}}) where D
  cs = [ Array{Int,1}(undef,(2^D,)) for i in 1:length(cells(grid)) ]
  for (i,ci) in enumerate(cells(grid))
    cs[i] .= ci
  end
  cs
end

function _compute_cells(grid::CartesianGrid{D},::Type{UnstructuredGrid{D,D}}) where D
  ptrs = fill(2^D,(length(cells(grid))+1,))
  length_to_ptrs!(ptrs)
  data = zeros(Int,ptrs[end]-1)
  _fill_cell_data!(data,cells(grid))
  (data, ptrs)
end

function _fill_cell_data!(data,cells)
  k = 1
  for v in cells
    for vi in v
      @inbounds data[k] = vi
      k +=1
    end
  end
end

end # module


text:
# See https://github.com/jeff-regier/Celeste.jl/wiki/About-SDSS-and-Stripe-82
struct SkyIntensity
    sky_small::Matrix{Float32} # background flux per pixel, in DNs
    sky_x::Vector{Float32} # interpolation coordinates
    sky_y::Vector{Float32}
    calibration::Vector{Float64} # nMgy per DN for each row
end

"""
Interpolate the 2-d array `sky_small` at the grid of array coordinates spanned
by the vectors `sky_x` and `sky_y` using bilinear interpolation.
The output array will have size `(length(sky_x), length(sky_y))`.
For example, if `x[1] = 3.3` and `y[2] = 4.7`, the element `out[1, 2]`
will be a result of linear interpolation between the values
`sky_small[3:4, 4:5]`.

For coordinates that are out-of-bounds (e.g., `sky_x[i] < 1.0` or
`sky_x[i] > size(sky_small,1)` where the interpolation would index sky_small values
outside the array, the nearest values in the sky_small array are used. (This is
constant extrapolation.)

Return value has units of nMgy.
"""
function interp_sky_kernel(sky::SkyIntensity, i::Int, j::Int)
    nx, ny = size(sky.sky_small)

    y0 = floor(Int, sky.sky_y[j])
    y1 = y0 + 1
    yw0 = sky.sky_y[j] - y0
    yw1 = 1.0f0 - yw0

    # modify out-of-bounds indicies to 1 or ny
    y0 = min(max(y0, 1), ny)
    y1 = min(max(y1, 1), ny)

    x0 = floor(Int, sky.sky_x[i])
    x1 = x0 + 1
    xw0 = sky.sky_x[i] - x0
    xw1 = 1.0f0 - xw0

    # modify out-of-bounds indicies to 1 or nx
    x0 = min(max(x0, 1), nx)
    x1 = min(max(x1, 1), nx)

    # bi-linear interpolation
    sky_dns = (xw0 * yw0 * sky.sky_small[x0, y0]
             + xw1 * yw0 * sky.sky_small[x1, y0]
             + xw0 * yw1 * sky.sky_small[x0, y1]
             + xw1 * yw1 * sky.sky_small[x1, y1])

    # return sky intensity in nMgy
    sky_dns * sky.calibration[i]
end

import Base.getindex

function getindex(sky::SkyIntensity, i::Int, j::Int)
    interp_sky_kernel(sky, i, j)
end

"""An image, taken though a particular filter band"""
mutable struct Image
    # The image height.
    H::Int

    # The image width.
    W::Int

    # An HxW matrix of pixel intensities, in raw electron counts (nelec).
    pixels::Matrix{Float32}

    # The band id (takes on values from 1 to 5).
    b::Int

    # World coordinates
    wcs::WCSTransform

    # The components of the point spread function.
    psf::Vector{PsfComponent}

    # SDSS-specific identifiers. A field is a particular region of the sky.
    # A Camcol is the output of one camera column as part of a Run.
    run_num::Int
    camcol_num::Int
    field_num::Int

    # The background noise in nanomaggies. (varies by position)
    sky::SkyIntensity

    # The expected number of photons contributed to this image
    # by a source 1 nanomaggie in brightness. (varies by row)
    nelec_per_nmgy::Array{Float32, 1}

    # storing a RawPSF here isn't ideal, because it's an SDSS type
    # not a Celeste type
    raw_psf_comp::RawPSF
end


text:
using Lynx
using Documenter

DocMeta.setdocmeta!(Lynx, :DocTestSetup, :(using Lynx); recursive=true)

makedocs(;
    modules=[Lynx],
    authors="Jorge Brito <jorge.brito.json@gmail.com>",
    repo="https://github.com/jorge-brito/Lynx.jl/blob/{commit}{path}#{line}",
    sitename="Lynx.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        assets=String[],
    ),
    pages=[
        "Introduction" => "index.md",
        "Tutorials" => [
            "Quick Start" => "tutorials/quickstart.md"
        ],
        "Reference" => [
            "Alphabetical function list"   =>  "reference/functionindex.md"
            "Function reference"           =>  "reference/api.md"
        ],
    ],
)

deploydocs(
    repo = "github.com/jorge-brito/Lynx.jl.git",
    target = "build"
)


text:
import Base: length, append!, sizehint!

@params struct TopOptTrace{T, TI <: Integer}
	c_hist::AbstractVector{T}
    v_hist::AbstractVector{T}
	x_hist::AbstractVector{<:AbstractVector{T}}
	add_hist::AbstractVector{TI}
	rem_hist::AbstractVector{TI}
end
TopOptTrace{T, TI}() where {T, TI<:Integer} = TopOptTrace(Vector{T}(), Vector{T}(), Vector{Vector{T}}(), Vector{TI}(), Vector{TI}())
length(t::TopOptTrace) = length(t.v_hist)

topopt_trace_fields = fieldnames(TopOptTrace)

macro append_fields_t1_t2()
    return esc(Expr(:block, [Expr(:call, :append!, :(t1.$f), :(t2.$f)) for f in topopt_trace_fields]...))
end
function append!(t1::TopOptTrace, t2::TopOptTrace)
    @append_fields_t1_t2()
end

macro sizehint!_fields_t()
    return esc(Expr(:block, [Expr(:call, :sizehint!, :(t.$f), :n) for f in topopt_trace_fields]...))
end
function sizehint!(t::TopOptTrace, n)
    @sizehint!_fields_t()
    return
end

function append!(ts::Vector{<:TopOptTrace})
    sizehint!(ts[1], sum(length.(ts)))
    for i in 2:length(ts)
        append!(ts[1], ts[i])
    end
    return ts[1]
end


text:
module FilterTestHelpers
using DSP, Test

export tffilter_eq, zpkfilter_eq, tffilter_accuracy, zpkfilter_accuracy,
       matrix_to_sosfilter, sosfilter_to_matrix

function lt(a, b)
    if abs(real(a) - real(b)) > 1e-10
        isless(real(a), real(b))
    else
        isless(imag(a), imag(b))
    end
end

function tffilter_eq(f1, f2)
    b1, a1 = (coefb(f1), coefa(f1))
    b2, a2 = (coefb(f2), coefa(f2))
    @test map(Float64, b1) ≈ map(Float64, b2)
    @test map(Float64, a1) ≈ map(Float64, a2)
end

function zpkfilter_eq(f1, f2)
    if !isempty(f1.z) || !isempty(f2.z)
        @test map(ComplexF64, sort(f1.z, lt=lt)) ≈ map(ComplexF64, sort(f2.z, lt=lt))
    end
    @test map(ComplexF64, sort(f1.p, lt=lt)) ≈ map(ComplexF64, sort(f2.p, lt=lt))
    @test map(Float64, f1.k) ≈ map(Float64, f2.k)
end

function zpkfilter_eq(f1, f2, eps)
    if !isempty(f1.z) || !isempty(f2.z)
        @test ≈(map(ComplexF64, sort(f1.z, lt=lt)), map(ComplexF64, sort(f2.z, lt=lt)), atol=eps)
    end
    @test ≈(map(ComplexF64, sort(f1.p, lt=lt)), map(ComplexF64, sort(f2.p, lt=lt)), atol=eps)
    @test ≈(map(Float64, f1.k), map(Float64, f2.k), atol=eps)
end

loss(x::Real, y::Real) = abs(float(x) - float(y))/eps(float(x))
loss(x::Union{Real,Complex}, y::Union{Real,Complex}) = loss(real(x), real(y)) + loss(imag(x), imag(y))
loss(x::AbstractVector, y::AbstractVector) = sum(map(loss, x, y))

function accuracy_check(err1, err2, part, relerr=1)
    try
        @test err1 <= relerr*err2
    catch e
        println("Filter 1 $part error (ULP): ", map(Float64, err1))
        println("Filter 2 $part error (ULP): ", map(Float64, err2))
        println("Ratio: ", map(Float64, err1/err2))
        rethrow(e)
    end
end

function tffilter_accuracy(f1, f2, accurate_f)
    b1, a1 = (coefb(f1), coefa(f1))
    b2, a2 = (coefb(f2), coefa(f2))
    accurate_b, accurate_a = (coefb(accurate_f), coefa(accurate_f))
    accuracy_check(loss(b1, accurate_b), loss(b2, accurate_b), "b")
    accuracy_check(loss(a1, accurate_a), loss(a2, accurate_a), "a")
end

function zpkfilter_accuracy(f1, f2, accurate_f; relerr=1, compare_gain_at=nothing, eps=nothing)
    z1, p1 = sort(f1.z, lt=lt), sort(f1.p, lt=lt)
    z2, p2 = sort(f2.z, lt=lt), sort(f2.p, lt=lt)
    accurate_z, accurate_p = sort(accurate_f.z, lt=lt), sort(accurate_f.p, lt=lt)
    if !isempty(z1) || !isempty(z2) || !isempty(accurate_z)
        if eps != nothing
            @test ≈(z1, accurate_z, atol=eps)
            @test ≈(z2, accurate_z, atol=eps)
        else
            @test z1 ≈ accurate_z
            @test z2 ≈ accurate_z
        end
        accuracy_check(loss(z1, accurate_z), loss(z2, accurate_z), "z", relerr)
    end
    if eps != nothing
        @test ≈(p1, accurate_p, atol=eps)
        @test ≈(p2, accurate_p, atol=eps)
        @test ≈(f1.k, accurate_f.k, atol=eps)
        @test ≈(f2.k, accurate_f.k, atol=eps)
    else
        @test p1 ≈ accurate_p
        @test p2 ≈ accurate_p
        @test f1.k ≈ accurate_f.k
        @test f2.k ≈ accurate_f.k
    end
    accuracy_check(loss(p1, accurate_p), loss(p2, accurate_p), "p", relerr)
    if compare_gain_at === nothing
        accuracy_check(loss(f1.k, accurate_f.k), loss(f2.k, accurate_f.k), "k", relerr)
    else
        jω = compare_gain_at*im
        gain1 = Float64(abs(f1.k * prod(Complex{BigFloat}.(z1) .- jω) / prod(Complex{BigFloat}.(p1) .- jω)))
        gain2 = Float64(abs(f2.k * prod(Complex{BigFloat}.(z2) .- jω) / prod(Complex{BigFloat}.(p2) .- jω)))
        accurate_gain = abs(accurate_f.k * prod(Complex{BigFloat}.(accurate_z) .- jω) / prod(Complex{BigFloat}.(accurate_p) .- jω))
        accuracy_check(loss(gain1, accurate_gain), loss(gain2, accurate_gain), "k", relerr)
    end
end

# Convert an SOS matrix, as returned by MATLAB, to an SecondOrderSections
matrix_to_sosfilter(sos::Matrix, g::Number) =
    SecondOrderSections([Biquad(sos[i, 1], sos[i, 2], sos[i, 3],
                            sos[i, 4], sos[i, 5], sos[i, 6])
               for i = 1:size(sos, 1)], g)

# Convert an SecondOrderSections to an SOS matrix as returned by MATLAB
function sosfilter_to_matrix(sos::SecondOrderSections)
    A = ones(length(sos.biquads), 6)
    for (i, biquad) in enumerate(sos.biquads)
        A[i, 1] = biquad.b0
        A[i, 2] = biquad.b1
        A[i, 3] = biquad.b2
        A[i, 5] = biquad.a1
        A[i, 6] = biquad.a2
    end
    A
end

# Show the poles and zeros in each biquad
# This is not currently used for testing, but is useful for debugging
function sosfilter_poles_zeros(sos::SecondOrderSections)
    z = fill(complex(nan(T)), 2, length(sos.biquads))
    p = fill(complex(nan(T)), 2, length(sos.biquads))
    for (i, biquad) in enumerate(sos.biquads)
        zpk = convert(ZeroPoleGain, biquad)
        z[1:length(zpk.z), i] = zpk.z
        p[1:length(zpk.p), i] = zpk.p
    end
    (z, p)
end
end


text:
### derivation free

@testset ":LN_NELDERMEAD, no scale" begin
    res0 = get_endpoint(
        [2.,2.,2.],
        (x) -> x[1]^2, # scan func
        (x) -> (x[1]-5.)^2 + (x[2] - 6.)^2, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 64.,
        local_alg = :LN_NELDERMEAD,
        silent = true
    )
    @test res0.status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0.counter < 1000
    @test isapprox(res0.value, 169.; atol=1e-2)
end

@testset ":LN_NELDERMEAD, no scale, !unpredictable" begin
    res0 = get_endpoint(
        [2.,2.,2.],
        (x) -> x[3]^2, # scan func
        (x) -> (x[1]-5.)^2 + (x[2] - 6.)^2, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 64.,
        scan_bound = 200.,
        local_alg = :LN_NELDERMEAD,
        silent = true
    )
    @test res0.status == :SCAN_BOUND_REACHED
    @test res0.counter < 1000
end

@testset ":LN_NELDERMEAD, log scale" begin
    res0 = get_endpoint(
        [2.,2.,2.],
        (x) -> x[1]^2, # scan func
        (x) -> (x[1]-5.)^2 + (x[2] - 6.)^2, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 64.,
        scale = [:log, :log, :log],
        local_alg = :LN_NELDERMEAD,
        silent = true
    )
    @test res0.status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0.counter < 1000
    @test isapprox(res0.value, 169.; atol=1e-2)
end

@testset ":LN_NELDERMEAD, log scale, !unpredictable" begin
    res0 = get_endpoint(
        [2.,2.,2.],
        (x) -> x[3]^2, # scan func
        (x) -> (x[1]-5.)^2 + (x[2] - 6.)^2, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 64.,
        scale = [:log, :log, :log],
        scan_bound = 200.,
        local_alg = :LN_NELDERMEAD,
        silent = true
    )
    @test res0.status == :SCAN_BOUND_REACHED
    @test res0.counter < 1000
end

### gradient-based

@testset ":LD_SLSQP, no scale, no scan_grad => error" begin
    @test_throws ArgumentError get_endpoint(
        [2.,2.,2.],
        (x) -> x[1]^2, # scan func
        (x) -> (x[1]-5.)^2 + (x[2] - 6.)^2, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 64.,
        local_alg = :LD_SLSQP,
        silent = true
    )
end

@testset ":LD_SLSQP, log scale, no scan_grad => error" begin
    @test_throws ArgumentError get_endpoint(
        [2.,2.,2.],
        (x) -> x[1]^2, # scan func
        (x) -> (x[1]-5.)^2 + (x[2] - 6.)^2, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 64.,
        scale = [:log, :log, :log],
        local_alg = :LD_SLSQP,
        silent = true
    )
end

## autodiff

@testset ":LD_SLSQP, no scale, :AUTODIFF, :AUTODIFF" begin
    res0 = get_endpoint(
        [2.,2.,2.],
        (x) -> x[1]^2, # scan func
        (x) -> (x[1]-5.)^2 + (x[2] - 6.)^2, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 64.,
        local_alg = :LD_SLSQP,
        scan_grad = :AUTODIFF,
        loss_grad = :AUTODIFF,
        silent = true
    )
    @test res0.status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0.counter < 1000
    @test isapprox(res0.value, 169.; atol=1e-2)
end

@testset ":LD_SLSQP, no scale, :AUTODIFF, :AUTODIFF, !unpredictable" begin
    res0 = get_endpoint(
        [2.,2.,2.],
        (x) -> x[3]^2, # scan func
        (x) -> (x[1]-5.)^2 + (x[2] - 6.)^2, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 64.,
        scan_bound = 200.,
        local_alg = :LD_SLSQP,
        scan_grad = :AUTODIFF,
        loss_grad = :AUTODIFF,
        silent = true
    )
    @test res0.status == :SCAN_BOUND_REACHED
    @test res0.counter < 1000
end

@testset ":LD_SLSQP, log scale, :AUTODIFF, :AUTODIFF" begin
    res0 = get_endpoint(
        [2.,2.,2.],
        (x) -> x[1]^2, # scan func
        (x) -> (x[1]-5.)^2 + (x[2] - 6.)^2, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 64.,
        scale = [:log, :log, :log],
        local_alg = :LD_SLSQP,
        scan_grad = :AUTODIFF,
        loss_grad = :AUTODIFF,
        silent = true
    )
    @test res0.status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0.counter < 1000
    @test isapprox(res0.value, 169.; atol=1e-2)
end

@testset ":LD_SLSQP, log scale, :AUTODIFF, :AUTODIFF, !unpredictable" begin
    res0 = get_endpoint(
        [2.,2.,2.],
        (x) -> x[3]^2, # scan func
        (x) -> (x[1]-5.)^2 + (x[2] - 6.)^2, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 64.,
        scan_bound = 200.,
        scale = [:log, :log, :log],
        local_alg = :LD_SLSQP,
        scan_grad = :AUTODIFF,
        loss_grad = :AUTODIFF,
        silent = true
    )
    @test res0.status == :SCAN_BOUND_REACHED
    @test res0.counter < 1000
end

@testset ":LD_SLSQP, no scale, :FINITE, :FINITE" begin
    res0 = get_endpoint(
        [2.,2.,2.],
        (x) -> x[1]^2, # scan func
        (x) -> (x[1]-5.)^2 + (x[2] - 6.)^2, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 64.,
        local_alg = :LD_SLSQP,
        scan_grad = :FINITE,
        loss_grad = :FINITE,
        silent = true
    )
    @test res0.status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0.counter < 1000
    @test isapprox(res0.value, 169.; atol=1e-2)
end

@testset ":LD_SLSQP, no scale, :FINITE, :FINITE, !unpredictable" begin
    res0 = get_endpoint(
        [2.,2.,2.],
        (x) -> x[3]^2, # scan func
        (x) -> (x[1]-5.)^2 + (x[2] - 6.)^2, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 64.,
        scan_bound = 200.,
        local_alg = :LD_SLSQP,
        scan_grad = :FINITE,
        loss_grad = :FINITE,
        silent = true
    )
    @test res0.status == :SCAN_BOUND_REACHED
    @test res0.counter < 1000
end

@testset ":LD_SLSQP, log scale, :FINITE, :FINITE" begin
    res0 = get_endpoint(
        [2.,2.,2.],
        (x) -> x[1]^2, # scan func
        (x) -> (x[1]-5.)^2 + (x[2] - 6.)^2, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 64.,
        scale = [:log, :log, :log],
        local_alg = :LD_SLSQP,
        scan_grad = :FINITE,
        loss_grad = :FINITE,
        silent = true
    )
    @test res0.status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0.counter < 1000
    @test isapprox(res0.value, 169.; atol=1e-2)
end

@testset ":LD_SLSQP, log scale, :FINITE, :FINITE, !unpredictable" begin
    res0 = get_endpoint(
        [2.,2.,2.],
        (x) -> x[3]^2, # scan func
        (x) -> (x[1]-5.)^2 + (x[2] - 6.)^2, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 64.,
        scan_bound = 200.,
        scale = [:log, :log, :log],
        local_alg = :LD_SLSQP,
        scan_grad = :FINITE,
        loss_grad = :FINITE,
        silent = true
    )
    @test res0.status == :SCAN_BOUND_REACHED
    @test res0.counter < 1000
end

@testset ":LD_SLSQP, no scale, Function, Function" begin
    res0 = get_endpoint(
        [2.,2.,2.],
        (x) -> x[1]^2, # scan func
        (x) -> (x[1]-5.)^2 + (x[2] - 6.)^2, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 64.,
        local_alg = :LD_SLSQP,
        scan_grad = (x) -> [2.0 * x[1], 0., 0.],
        loss_grad = (x) -> [2.0 * (x[1] - 5.), 2.0 * (x[2] - 6.), 0.],
        silent = true
    )
    @test res0.status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0.counter < 1000
    @test isapprox(res0.value, 169.; atol=1e-2)
end

@testset ":LD_SLSQP, no scale, Function, Function, !unpredictable" begin
    res0 = get_endpoint(
        [2.,2.,2.],
        (x) -> x[3]^2, # scan func
        (x) -> (x[1]-5.)^2 + (x[2] - 6.)^2, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 64.,
        scan_bound = 200.,
        local_alg = :LD_SLSQP,
        scan_grad = (x) -> [0., 0., 2.0 * x[3]],
        loss_grad = (x) -> [2.0 * (x[1] - 5.), 2.0 * (x[2] - 6.), 0.],
        silent = true
    )
    @test res0.status == :SCAN_BOUND_REACHED
    @test res0.counter < 1000
end

@testset ":LD_SLSQP, log scale, Function, Function" begin
    res0 = get_endpoint(
        [2.,2.,2.],
        (x) -> x[1]^2, # scan func
        (x) -> (x[1]-5.)^2 + (x[2] - 6.)^2, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 64.,
        scale = [:log, :log, :log],
        local_alg = :LD_SLSQP,
        scan_grad = (x) -> [2.0 * x[1], 0., 0.],
        loss_grad = (x) -> [2.0 * (x[1] - 5.), 2.0 * (x[2] - 6.), 0.],
        silent = true
    )
    @test res0.status == :BORDER_FOUND_BY_SCAN_TOL
    @test res0.counter < 1000
    @test isapprox(res0.value, 169.; atol=1e-2)
end

@testset ":LD_SLSQP, log scale, Function, Function, !unpredictable" begin
    res0 = get_endpoint(
        [2.,2.,2.],
        (x) -> x[3]^2, # scan func
        (x) -> (x[1]-5.)^2 + (x[2] - 6.)^2, # loss_func
        :CICO_ONE_PASS,
        loss_crit = 64.,
        scan_bound = 200.,
        scale = [:log, :log, :log],
        local_alg = :LD_SLSQP,
        scan_grad = (x) -> [0., 0., 2.0 * x[3]],
        loss_grad = (x) -> [2.0 * (x[1] - 5.), 2.0 * (x[2] - 6.), 0.],
        silent = true
    )
    @test res0.status == :SCAN_BOUND_REACHED
    @test res0.counter < 1000
end


text:
using Documenter
using Neet

makedocs(
    sitename = "Neet",
    format = Documenter.HTML(),
    modules = [Neet],
    authors="Harrison Smith and Douglas G. Moore",
    pages = Any[
        "Home" => "index.md",
        "core.md"
    ]
)

deploydocs(
    repo = "github.com/dglmoore/Neet.jl.git"
)


text:
using Gnuplot
t = -3.1:0.01:3.1
u = 16*sin.(t).^3
v = 13*cos.(t) .- 5*cos.(2t) .- 2*cos.(3t) .- cos.(4t)
@gp xlab = "u" ylab = "v" :-
@gp :- u v "w filledcu tit 'heart' lc 'black'" "set size square"
@gp :- "set style fill transparent solid 0.8"
save(term = "pngcairo size 600,400", output = "heart.png")

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
function MultivariateFunctionMappingAdapter(arg0::MultivariateFunction, arg1::Vector{jdouble}, arg2::Vector{jdouble})
    return MultivariateFunctionMappingAdapter((MultivariateFunction, Vector{jdouble}, Vector{jdouble}), arg0, arg1, arg2)
end

function bounded_to_unbounded(obj::MultivariateFunctionMappingAdapter, arg0::Vector{jdouble})
    return jcall(obj, "boundedToUnbounded", Vector{jdouble}, (Vector{jdouble},), arg0)
end

function unbounded_to_bounded(obj::MultivariateFunctionMappingAdapter, arg0::Vector{jdouble})
    return jcall(obj, "unboundedToBounded", Vector{jdouble}, (Vector{jdouble},), arg0)
end

function value(obj::MultivariateFunctionMappingAdapter, arg0::Vector{jdouble})
    return jcall(obj, "value", jdouble, (Vector{jdouble},), arg0)
end


text:
function __init__()
    @require TaylorModels = "314ce334-5f6e-57ae-acf6-00b6e903104a" include("init_TaylorModels.jl")
    @require IntervalMatrices = "5c1f47dc-42dd-5697-8aaa-4d102d140ba9" include("init_IntervalMatrices.jl")
    @require IntervalConstraintProgramming = "138f1668-1576-5ad7-91b9-7425abbf3153" include("init_IntervalConstraintProgramming.jl")
end


text:
module KalmanFilter

using Calculus

function predict_int!(statevec, xhat_next::Vector, F::Matrix, Q::Matrix)
    statevec.P = F*statevec.P*F' + Q
    statevec.xhat = xhat_next
end

function predict!(statevec, F::Matrix, Q::Matrix)
    predict_int!(statevec, F*statevec.xhat, F, Q)
end

function predict!(statevec, f::Function, fp::Function, Q::Matrix)
    predict_int!(statevec, f(statevec.xhat), fp(statevec.xhat), Q)
end

function predict!(statevec, f::Function, Q::Matrix)
    predict_int!(statevec, f(statevec.xhat), jacobian(f)(statevec.xhat), Q)
end

function update_int!(statevec, y::Vector, yhat::Vector, H::Matrix, R::Matrix)
    e = y - yhat
    S = H*statevec.P*H' + R
    K = statevec.P*H'/S
    statevec.P = statevec.P - K*H*statevec.P
    statevec.xhat = statevec.xhat + K*e
end

function update!(statevec, y::Vector, H::Matrix, R::Matrix)
    update_int!(statevec, y, H*statevec.xhat, H, R)
end

function update!(statevec, y::Vector, h::Function, hp::Function, R::Matrix)
    update_int!(statevec, y, h(statevec.xhat), hp(statevec.xhat), R)
end

function update!(statevec, y::Vector, h::Function, R::Matrix)
    update_int!(statevec, y, h(statevec.xhat), jacobian(h)(statevec.xhat), R)
end

end


text:
using SafeTestsets

@safetestset "Test Cell Maps" begin
    include("test_cell_map.jl")
end

@safetestset "Test Hooke Stiffness" begin
    include("test_elasticity.jl")
end

@safetestset "Test Mesh Routines" begin
    include("test_mesh_routines.jl")
end

@safetestset "Test Weak Form" begin
    include("test_weak_form.jl")
end

@safetestset "Test Assembly" begin
    include("test_assembly.jl")
end

@safetestset "Test Full Displacement BC Convergence" begin
    include("test_full_displacement_convergence.jl")
end

@safetestset "Test Uniform Tension" begin
    include("test_traction_bc.jl")
end

@safetestset "Test Displacement + Traction BC Convergence" begin
    include("test_displacement_traction_convergence.jl")
end

@safetestset "Test CutMesh Initialization" begin
    include("test_cut_mesh.jl")
end

@safetestset "Test DG Mesh" begin
    include("test_dg_mesh.jl")
end

@safetestset "Test Cut Cell Quadratures" begin
    include("test_cell_quadratures.jl")
end

@safetestset "Test Face Quadrature Rules" begin
    include("test_face_quadratures.jl")
end

@safetestset "Test Interface Quadratures" begin
    include("test_interface_quadratures.jl")
end

@safetestset "Test Cut Mesh Bilinear Forms" begin
    include("test_cut_mesh_bilinear_forms.jl")
end

@safetestset "Test Mean Traction Interface Condition" begin
    include("test_mean_traction_flux.jl")
end

@safetestset "Test Interface Conditions" begin
    include("test_interface_conditions.jl")
    include("test_incoherent_condition.jl")
end

@safetestset "Test Conceptual Merging" begin
    include("test_extended_mean_flux.jl")
end

@safetestset "Test Cell Merging" begin
    include("test_cell_merging.jl")
end

@safetestset "Test Penalty Displacement BC" begin
    include("test_penalty_displacement_bc.jl")
end

@safetestset "Test Penalty Displacement Component BC" begin
    include("test_displacement_component_bc.jl")
end

@safetestset "Test Cut Mesh Convergence" begin
    include("test_cut_convergence.jl")
end

@safetestset "Test Penalty Displacement Component BC Convergence" begin
    include("test_displacement_component_bc_convergence.jl")
end

@safetestset "Test Cut Mesh Stress Convergence" begin
    include("test_cut_stress_convergence.jl")
end

@safetestset "Test Transformation Strain Simple Tension" begin
    include("test_transformation_strain_simple_tension.jl")
end

@safetestset "Test Transformation Strain - Displacement Convergence" begin
    include("test_transformation_strain_displacement_convergence.jl")
    include("test_incoherent_transformation_strain_displacement_convergence.jl")
end

@safetestset "Test Transformation Strain - Stress Convergence" begin
    include("test_transformation_strain_stress_convergence.jl")
    include("test_incoherent_transformation_strain_stress_convergence.jl")
end

@safetestset "Test Transformation Strain - Displacement + Traction BC convergence" begin
    include("test_transformation_strain_displacement_traction_bc_convergence.jl")
    include("test_incoherent_transformation_strain_displacement_traction_bc_convergence.jl")
end

@safetestset "Test Transformation Strain - Mixed BC convergence" begin
    include("test_transformation_strain_mixed_bc_convergence.jl")
    include("test_incoherent_transformation_strain_mixed_bc_convergence.jl")
end

@safetestset "Test Levelset Reinitialization" begin
    include("test_levelset_reinitialize.jl")
end

@safetestset "Test Levelset Propagate" begin
    include("test_levelset_propagate.jl")
end


text:
# Autogenerated wrapper script for HDF5_jll for x86_64-linux-gnu
export libhdf5, libhdf5_hl

using Zlib_jll
using OpenSSL_jll
using LibCURL_jll
JLLWrappers.@generate_wrapper_header("HDF5")
JLLWrappers.@declare_library_product(libhdf5, "libhdf5-00e8fae8.so.200.0.0")
JLLWrappers.@declare_library_product(libhdf5_hl, "libhdf5_hl-383c339f.so.200.0.0")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll, OpenSSL_jll, LibCURL_jll)
    JLLWrappers.@init_library_product(
        libhdf5,
        "lib/libhdf5.so.200.0.0",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libhdf5_hl,
        "lib/libhdf5_hl.so.200.0.0",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
struct LinSolveWs{T<:Number,U<:Integer}
    lu::Vector{T}
    ipiv::Vector{BlasInt}
    function LinSolveWs{T,U}(n::U) where {T<:Number,U<:Integer}
        lu = zeros(T, n * n)
        ipiv = zeros(BlasInt, n)
        new(lu, ipiv)
    end
end

# Float64 is the default
LinSolveWs(n) = LinSolveWs{Float64,Int64}(n)

strides(a::Adjoint) = strides(a.parent)

for (getrf, getrs, elty) in (
    (:dgetrf_, :dgetrs_, :Float64),
    (:sgetrf_, :sgetrs_, :Float32),
    (:zgetrf_, :zgetrs_, :ComplexF64),
    (:cgetrf_, :cgetrs_, :ComplexF32),
)
    @eval begin
        function lu!(a::StridedMatrix{$elty}, ws::LinSolveWs)
            copyto!(ws.lu, a)
            mm, nn = size(a)
            m = Ref{BlasInt}(mm)
            n = Ref{BlasInt}(nn)
            # ws.lu isn't a view and has continuous storage
            lda = Ref{BlasInt}(max(1, mm))
            info = Ref{BlasInt}(0)
            ccall(
                (@blasfunc($getrf), liblapack),
                Cvoid,
                (
                    Ref{BlasInt},
                    Ref{BlasInt},
                    Ptr{$elty},
                    Ref{BlasInt},
                    Ptr{BlasInt},
                    Ref{BlasInt},
                ),
                m,
                n,
                ws.lu,
                lda,
                ws.ipiv,
                info,
            )
            if info[] != 0
                chklapackerror(info[])
            end
        end

        t1 = StridedMatrix{$elty}
        t2 = Transpose{$elty,<:StridedMatrix}
        t3 = Adjoint{$elty,<:StridedMatrix}
    end

    for (elty2, transchar) in ((t1, 'N'), (t2, 'T'), (t3, 'C'))

        @eval begin
            function linsolve_core_no_lu!(
                a::$elty2,
                b::StridedVecOrMat{$elty},
                ws::LinSolveWs,
            )
                mm, nn = size(a)
                m = Ref{BlasInt}(mm)
                n = Ref{BlasInt}(nn)
                nhrs = Ref{BlasInt}(size(b, 2))
                # ws.lu isn't a view and has continuous storage
                lda = Ref{BlasInt}(max(1, mm))
                ldb = Ref{BlasInt}(max(1, stride(b, 2)))
                info = Ref{BlasInt}(0)
                ccall(
                    (@blasfunc($getrs), liblapack),
                    Cvoid,
                    (
                        Ref{UInt8},
                        Ref{BlasInt},
                        Ref{BlasInt},
                        Ptr{$elty},
                        Ref{BlasInt},
                        Ptr{BlasInt},
                        Ptr{$elty},
                        Ref{BlasInt},
                        Ref{BlasInt},
                    ),
                    $transchar,
                    n,
                    nhrs,
                    ws.lu,
                    lda,
                    ws.ipiv,
                    b,
                    ldb,
                    info,
                )
                if info[] != 0
                    chklapackerror(info[])
                end
            end

            function linsolve_core!(a::$elty2, b::StridedVecOrMat{$elty}, ws::LinSolveWs)

                lu!(a, ws)
                linsolve_core_no_lu!(a, b, ws)
            end
        end
    end
end

# if a is an adjoint matrix, we compute lu decomposition
# for its parent as the transposition is done at the
# solution stage
lu!(a::Adjoint, ws::LinSolveWs) = lu!(a.parent, ws)


text:
struct VehiclePath{N}
    x_coefs::AbstractVector
    y_coefs::AbstractVector
    s_max::AbstractFloat
    function VehiclePath(x_coefs::AbstractVector, y_coefs::AbstractVector, s_max::AbstractFloat)
        @assert length(x_coefs)==length(y_coefs) "Path length mismatch"
        N = length(x_coefs)
        new{N}(x_coefs, y_coefs, s_max)
    end
end

# Define the dynamics model of the model.
struct RoundaboutGame{N,M,P,SVu,SVx,SVz} <: AbstractGameModel
    n::Int  # Number of states
    m::Int  # Number of controls
    p::Int  # Number of players
    ni::Vector{Int}  # Number of states for each player
    mi::Vector{Int}  # Number of controls for each player
    pu::SVu # Indices of the each player's controls
    px::SVx # Indices of the each player's x and y positions
    pz::SVz # Indices of the each player's states

    paths::Vector{<:VehiclePath} # List of Path objects for each vehicle
    uc::Int # Number of uncontrolled players (default at the end of the player list)
end

function RoundaboutGame(paths::Vector{VehiclePath{N}}, uc::Int) where {N}
    # p = number of players
    p = length(paths)
    n = 4p
    m = p
    pu = [SVector{1,Int}([i]) for i=1:p]
    px = [SVector{2,Int}([i + (j-1)*p for j=1:2]) for i=1:p]
    pz = [SVector{4,Int}([i + (j-1)*p for j=1:4]) for i=1:p]
    TYPE = typeof.((pu,px,pz))
    ni = 4*ones(Int,p)
    mi = 1*ones(Int,p)
    return RoundaboutGame{n,m,p,TYPE...}(n,m,p,ni,mi,pu,px,pz,paths,uc)
end

#@generated function RobotDynamics.dynamics(model::UnicycleGame{N,M,P}, x, u) where {N,M,P}
#    xd  = [:(cos(x[$i])*x[$i+P]) for i=M+1:M+P]
#    yd  = [:(sin(x[$i])*x[$i+P]) for i=M+1:M+P]
#    qdd = [:(u[$i]) for i=1:M]
#    return :(SVector{$N}($(xd...), $(yd...), $(qdd...)))
#end

function powerseries(coefs::AbstractVector, s)
    maxd = length(coefs)-1
    ss = SVector{maxd+1}([s^i for i=0:maxd])
    return dot(coefs, ss)
end

function powerseriesd(coefs::AbstractVector, s, v)
    maxd = length(coefs)-1
    ss = SVector{maxd+1}([i*s^(i-1) for i=0:maxd])
    return v*dot(coefs, ss)
end

function RobotDynamics.discrete_dynamics(::Type{PassThrough}, model::RoundaboutGame{N,M,P},
        x::StaticVector, u::StaticVector, t, dt) where {N,M,P}

    sidxs = SVector{P}([i for i=(2P+1):3P])
    st = x[sidxs]
    vt = x[sidxs.+P]

    vn = vt + u*dt
    sn = st + dt*(vt+vn)/2

    xn = SVector{P}([powerseries(model.paths[i].x_coefs, sn[i]) for i=1:P])
    yn = SVector{P}([powerseries(model.paths[i].y_coefs, sn[i]) for i=1:P])
    return SVector{N}(xn..., yn..., sn..., vn...)
end

function RobotDynamics.dynamics(model::RoundaboutGame{N,M,P}, x, u) where {N,M,P}
    sidxs = SVector{P}([i for i=(2P+1):3P])
    st = x[sidxs]
    vt = x[sidxs.+P]

    vd = u
    sd = vt

    xd = SVector{P}([powerseriesd(model.paths[i].x_coefs, st[i], vt[i]) for i=1:P])
    yd = SVector{P}([powerseriesd(model.paths[i].y_coefs, st[i], vt[i]) for i=1:P])
    return SVector{N}(xd..., yd..., sd..., vd...)
end

max_stars_repo_name: lassepe/DecNashPlanning

max_issues_repo_name: lassepe/DecNashPlanning

max_forks_repo_name: lassepe/DecNashPlanning


text:
function ksintegrateInplace(u, Lx, dt, Nt, nsave);
    u = (1+0im)*u                       # force u to be complex
    Nx = length(u)                      # number of gridpoints
    kx = vcat(0:Nx/2-1, 0:0, -Nx/2+1:-1)# integer wavenumbers: exp(2*pi*kx*x/L)
    alpha = 2*pi*kx/Lx                  # real wavenumbers:    exp(alpha*x)
    D = 1im*alpha                       # spectral D = d/dx operator
    L = alpha.^2 - alpha.^4             # spectral L = -D^2 - D^4 operator
    G = -0.5*D                          # spectral -1/2 D operator

    Nsave = div(Nt, nsave)+1            # number of saved time steps
    t = (0:Nsave)*(dt*nsave)            # t timesteps
    U = zeros(Nsave, Nx)                # matrix of u(xⱼ, tᵢ) values
    U[1,:] = u                          # assign initial condition to U
    s = 2                               # counter for saved data

    # some convenience variables
    dt2  = dt/2
    dt32 = 3*dt/2
    A_inv = (ones(Nx) - dt2*L).^(-1)
    B     =  ones(Nx) + dt2*L

    # compute in-place FFTW plans
    FFT! = plan_fft!(u, flags=FFTW.ESTIMATE)
    IFFT! = plan_ifft!(u, flags=FFTW.ESTIMATE)

    # compute nonlinear term Nn = -u u_x
    Nn  = G.*fft(u.^2);    # Nn = -1/2 d/dx u^2 = -u u_x
    Nn1 = copy(Nn);        # Nn1 = Nn at first time step
    FFT!*u;                # transform u physical -> spectral

    # timestepping loop
    for n = 1:Nt

        Nn1 .= Nn       # shift nonlinear term in time

        Nn .= u         # put u into Nn
        IFFT!*Nn;       # transform Nn = u to gridpt values, in place
        Nn .= Nn.*Nn;   # collocation calculation, set Nn = u^2
        FFT!*Nn;        # transform Nn = u^2 back to spectral coeffs
        Nn .= G.*Nn;    # apply G = -1/2d/dx  to compute Nn = -1/2 d/dx (u^2)

        # loop fusion! Julia translates this line into a single for-loop on
        #   u[i] = A_inv[i] * (B[i]*u[i] + dt32*Nn[i] - dt2*Nn1[i];
        # no temporary vectors!

        u .= A_inv .* (B .* u .+ dt32.*Nn .- dt2.*Nn1);

        if mod(n, nsave) == 0
            U[s,:] = real(ifft(u))
            s += 1
        end
    end

    t,U
end

max_stars_repo_name: raphbacher/julia-intro

max_issues_repo_name: raphbacher/julia-intro

max_forks_repo_name: raphbacher/julia-intro


text:
include("/Users/guilhermebodin/Documents/PUC/MIPTests_Branch_Guilherme/runmiptests.jl")
include("/Users/guilhermebodin/Documents/PUC/MIPSolver2017.jl/src/branch_and_bound.jl")
using Gurobi

runtests(solveMIP, GurobiSolver(OutputFlag = 0))


text:
# matrix
# SymMatrix has methods to expose.
# XXX has issue with @inferred
# XXX link into pyconvert...

Base.promote_op(::T, ::Type{S}, ::Type{Sym}) where {T, S <: Number} = Sym
Base.promote_op(::T, ::Type{Sym}, ::Type{S}) where {T, S <: Number} = Sym
Base.promote_op(::T, ::Type{Sym}, ::Type{Sym}) where {T} = Sym # This helps out linear alg

function unSym(M::AbstractMatrix{T}) where {T <: Sym}
    sympy.py.Matrix(Tuple(Mᵢ for Mᵢ ∈ eachrow(M)))
end

call_matrix_meth(M::Matrix{T}, meth::Symbol, args...; kwargs...) where {T <: Sym} =
    ↑(getproperty(↓(M), meth)(↓(args)...; unSymkwargs(kwargs)...))

# include those exported from base or LinearAlgebra with corresponding sympy
matrix_meths = (
    (:adjoint, :Base, :adjoint),
    (:exp, :Base, :exp),
    (:inv, :Base, :inv),
    (:transpose, :Base, :transpose),
    (:det, :LinearAlgebra, :det),
    (:diag, :LinearAlgebra, :diag),
    (:eigenvects, :LinearAlgebra, :eigvecs),
    (:eigenvals, :LinearAlgebra, :eigvals),
    (:norm, :LinearAlgebra, :norm),
    (:pinv, :LinearAlgebra, :pinv),
    (:rank, :LinearAlgebra, :rank),
)

for (smeth, jmod, jmeth) ∈ matrix_meths
    @eval begin
        ($(jmod).$(jmeth))(M::Matrix{T}, args...; kwargs...) where {T <: Sym} =
            ↑(↓(M).$(smeth)(↓(args)...; unSymkwargs(kwargs)...))
    end
end

function Base.getproperty(M::Matrix{Sym}, prop::Symbol)
    val = pygetattr(↓(M), string(prop))
    meth = pygetattr(val, "__call__", nothing)
    meth == nothing && return Sym(val)
    return SymbolicCallable(meth)
end

function Base.getproperty(X::Vector{Sym}, prop::Symbol)
    val = pygetattr(sympy.Matrix(Tuple(xᵢ for xᵢ ∈ X)), string(prop))
    meth = pygetattr(val, "__call__", nothing)
    meth == nothing && return Sym(val)
    return SymbolicCallable(meth)
end

# solve Ax=b for x, avoiding generic `lu`, which can be very slow for bigger n values
# fix suggested by @olof3 in issue 355
function LinearAlgebra.:\(A::AbstractArray{Sym,2}, b::AbstractArray{S,1}) where {S}

    m,n  = size(A)
    x =  [Sym("x$i") for  i in 1:n]
    out = solve(A*x-b, x)
    isempty(out) && throw(SingularException(0)) # Could also return out here?
    ret = Vector{Sym}(undef, n)
    for (i,xᵢ)  in enumerate(x)
        ret[i] =  get(out,  xᵢ, xᵢ)
    end

    return ret

end

function LinearAlgebra.:\(A::AbstractArray{T,2}, B::AbstractArray{S,2}) where {T <: Sym, S}
    hcat([A \ bⱼ for bⱼ in eachcol(B)]...)
end

LinearAlgebra.norm(x::AbstractVector{T}, args...; kwargs...) where {T <: SymbolicObject} =
    ↑(getproperty(sympy.Matrix(Tuple(xᵢ for xᵢ ∈ x)), :norm)(↓(args)...; unSymkwargs(kwargs)...))

## Issue #359 so that A  +  λI is of type Sym
Base.:+(A::AbstractMatrix{T}, J::UniformScaling) where {T <: SymbolicObject} = (n=LinearAlgebra.checksquare(A); A .+ J.λ*I(n))
Base.:+(A::AbstractMatrix, J::UniformScaling{T}) where {T <: SymbolicObject} = (n=LinearAlgebra.checksquare(A); A .+ J.λ*I(n))
Base.:+(A::AbstractMatrix{T}, J::UniformScaling{T}) where {T <: SymbolicObject} = (n=LinearAlgebra.checksquare(A); A .+ J.λ*I(n))

Base.:-(J::UniformScaling, A::AbstractMatrix{T}) where {T <: SymbolicObject} = (-A) + J
Base.:-(J::UniformScaling{T}, A::AbstractMatrix) where {T <: SymbolicObject} = (-A) + J
Base.:-(J::UniformScaling{T}, A::AbstractMatrix{T}) where {T <: SymbolicObject} = (-A) + J

max_stars_repo_name: jverzani/SymPyCall

max_issues_repo_name: jverzani/SymPyCall

max_forks_repo_name: jverzani/SymPyCall


text:
module DTMarkovChain

include("coretypes.jl")
export MarkovChain,
  MarkovChainFitParams,
  DTFiniteMarkovChain,
  DTFiniteMarkovChainWithFit,
  DTFiniteMarkovChainMleFitParams

include("finite.jl")
export fit,
  fit!

end # module