text:
"""
	fourier_constant(src::SourceParameters)

Define the constant source term for the Fourier Amplitude Spectrum.
Constant set to permit distances to be passed in km, densities in t/m^3, and velocities in km/s
The reference distance is set to 1.0 km (and interpreted to be a rupture distance).
"""
function fourier_constant(src::SourceParameters)
	# note that the 1e-20 factor allows one to pass in distances in km, density in t/m^3, velocity in km/s
	Rref0 = 1.0
	return src.RΘϕ * src.V * src.F * 1e-20 / ( 4π * src.ρ * src.β^3 * Rref0 )
end

fourier_constant(fas::FourierParameters) = fourier_constant(fas.source)

"""
	fourier_source_shape(f::T, m::S, src::SourceParameters) where {S<:Real,T<:Float64}

Source shape of the Fourier Amplitude Spectrum of _displacement_, without the constant term or seismic moment.
This simply includes the source spectral shape.
"""
function fourier_source_shape(f::T, m::S, src::SourceParameters) where {S<:Real,T<:Float64}
	fa, fb, ε = corner_frequency(m, src)
	if src.model == :Brune
		# use single corner, with fa=fc
		return 1.0 / (1.0 + (f/fa)^2)
	elseif src.model == :Atkinson_Silva_2000
		# use double corner model
		return (1.0 - ε)/(1.0 + (f/fa)^2 ) + ε/(1.0 + (f/fb)^2)
	else
		# use single corner, with fa=fc
		return 1.0 / (1.0 + (f/fa)^2)
	end
end

"""
	fourier_source_shape(f, m, fas::FourierParameters)

Source shape of the Fourier amplitude spectrum of _displacement_, without the constant term or seismic moment.
Defined using a `FourierParameters` instance for the source model.

See also: [`fourier_source_shape`](@ref)
"""
fourier_source_shape(f, m, fas::FourierParameters) = fourier_source_shape(f, m, fas.source)

"""
	fourier_source_shape(f::Float64, fa::T, fb::T, ε::T, model::Symbol) where T<:Real

Fourier amplitude spectral shape for _displacement_ defined by corner frequencies.

See also: [`fourier_source_shape`](@ref)
"""
function fourier_source_shape(f::Float64, fa::T, fb::T, ε::T, model::Symbol) where T<:Real
	if model == :Brune
		# use single corner, with fa=fc
		return 1.0 / (1.0 + (f/fa)^2)
	elseif model == :Atkinson_Silva_2000
		# use double corner model
		return (1.0 - ε)/(1.0 + (f/fa)^2 ) + ε/(1.0 + (f/fb)^2)
	else
		# use single corner, with fa=fc
		return 1.0 / (1.0 + (f/fa)^2)
	end
end

"""
	fourier_source(f::T, m::S, src::SourceParameters) where {S<:Real,T<:Float64}

Source Fourier Amplitude Spectrum of displacement, without the constant term.
This simply includes the seismic moment and the source spectral shape.

See also: [`fourier_source_shape`](@ref)
"""
function fourier_source(f::T, m::S, src::SourceParameters) where {S<:Real,T<:Float64}
	Mo = magnitude_to_moment(m)
	return Mo * fourier_source_shape(f, m, src)
end

"""
	fourier_source(f, m, fas::FourierParameters)

Source Fourier Amplitude Spectrum of displacement, without the constant term.
This simply includes the seismic moment and the source spectral shape.
Defined using a `FourierParameters` instance for the source model.

See also: [`fourier_source_shape`](@ref)
"""
fourier_source(f, m, fas::FourierParameters) = fourier_source(f, m, fas.source)

"""
	fourier_path(f::U, r_ps::T, m::S, geo::GeometricSpreadingParameters, sat::NearSourceSaturationParameters, ane::AnelasticAttenuationParameters) where {S<:Real,T<:Real,U<:Float64}

Path scaling of Fourier spectral model -- combination of geometric spreading and anelastic attenuation.

See also: [`geometric_spreading`](@ref), [`anelastic_attenuation`](@ref)
"""
function fourier_path(f::U, r_ps::T, m::S, geo::GeometricSpreadingParameters, sat::NearSourceSaturationParameters, ane::AnelasticAttenuationParameters) where {S<:Real,T<:Real,U<:Float64}
	Zr = geometric_spreading(r_ps, m, geo, sat)
	if ane.rmetric == :Rrup
		r_rup = rupture_distance_from_equivalent_point_source_distance(r_ps, m, sat)
		Qr = anelastic_attenuation(f, r_rup, ane)
	else
		Qr = anelastic_attenuation(f, r_ps, ane)
	end
	return Zr * Qr
end

fourier_path(f, r_ps, m, path::PathParameters) = fourier_path(f, r_ps, m, path.geometric, path.saturation, path.anelastic)
fourier_path(f, r_ps, m, fas::FourierParameters) = fourier_path(f, r_ps, m, fas.path)
# simplified cases where m is not required
fourier_path(f, r_ps, path::PathParameters) = fourier_path(f, r_ps, NaN, path.geometric, path.saturation, path.anelastic)
fourier_path(f, r_ps, fas::FourierParameters) = fourier_path(f, r_ps, fas.path)

"""
	fourier_attenuation(f::S, r::T, ane::AnelasticAttenuationParameters{U,V}, site::SiteParameters{W}) where {S<:Float64,T<:Real,U<:Real,V<:Real,W<:Real}

Combined full-path attenuation, including `Q(f)` effects and `κ0` filter for frequency `f`
Distance defined in terms of an equivalent point source distance `r_ps` or rupture distance `r_rup` depending upon what metric is defined in `ane.rmetric`
"""
function fourier_attenuation(f::S, r::T, ane::AnelasticAttenuationParameters{U,V}, site::SiteParameters{W}) where {S<:Float64,T<:Real,U<:Real,V<:Real,W<:Real}
	if f < eps()
		return oneunit(T)*oneunit(U)*oneunit(V)*oneunit(W)
	else
		Qr = anelastic_attenuation(f, r, ane)
		Kf = kappa_filter(f, site)
		return Qr * Kf
	end
end

fourier_attenuation(f::S, r::T, path::PathParameters, site::SiteParameters) where {S<:Float64,T<:Real} = fourier_attenuation(f, r, path.anelastic, site)
fourier_attenuation(f::S, r::T, fas::FourierParameters) where {S<:Float64,T<:Real} = fourier_attenuation(f, r, fas.path, fas.site)

"""
	fourier_site(f::Float64, site::SiteParameters)

Combined site amplification and kappa filter for frequency `f`
"""
function fourier_site(f::Float64, site::SiteParameters)
	# kappa filter
    Kf = kappa_filter(f, site)
    # site impedance function
    Sf = site_amplification(f, site)
	return Sf * Kf
end

fourier_site(f::Float64, fas::FourierParameters) = fourier_site(f, fas.site)

"""
	fourier_spectral_ordinate(f::S, m::S, r_ps::T, src::SourceParameters, geo::GeometricSpreadingParameters, ane::AnelasticAttenuationParameters, site::SiteParameters) where {S<:Float64,T<:Real}

Fourier acceleration spectral ordinate (m/s) based upon an equivalent point source distance `r_ps`
- `f` is frequency (Hz)
- `m` is magnitude
- `r_ps` is the equivalent point source distance including saturation effects (km)
- `src` are the source parameters `SourceParameters`
- `geo` are the geometric spreading parameters `GeometricSpreadingParameters`
- `sat` are the near source saturation parameters `NearSourceSaturationParameters`
- `ane` are the anelastic attenuation parameters `AnelasticAttenuationParameters`
- `site` are the site parameters `SiteParameters`

See also: [`fourier_spectrum`](@ref), [`fourier_spectrum!`](@ref)
"""
function fourier_spectral_ordinate(f::U, m::S, r_ps::T, src::SourceParameters, geo::GeometricSpreadingParameters, sat::NearSourceSaturationParameters, ane::AnelasticAttenuationParameters, site::SiteParameters) where {S<:Real,T<:Real,U<:Float64}
	# define all constant terms here
  	C = fourier_constant(src)
	# source term
	Ef = fourier_source(f, m, src)
	# geometric spreading
	Gr = geometric_spreading(r_ps, m, geo, sat)
	# combined attenuation of both path (κr) and site (κ0)
	if ane.rmetric == :Rrup
		r_rup = rupture_distance_from_equivalent_point_source_distance(r_ps, m, sat)
		Kf = fourier_attenuation(f, r_rup, ane, site)
	else
		Kf = fourier_attenuation(f, r_ps, ane, site)
	end
	# site impedance
    Sf = site_amplification(f, site)

	# overall displacement spectrum (in cms) is Ef*Pf*Sf
	# the division by 10^2 below is to convert from units of cm/s to m/s
	return (2π*f)^2 * C * Ef * Gr * Kf * Sf / 100.0
end

"""
	fourier_spectral_ordinate(f::U, m::S, r_ps::T, src::SourceParameters, path::PathParameters, site::SiteParameters) where {S<:Real,T<:Real,U<:Float64}

Fourier acceleration spectral ordinate (m/s) based upon an equivalent point source distance `r_ps`
- `f` is frequency (Hz)
- `m` is magnitude
- `r_ps` is the equivalent point source distance including saturation effects (km)
- `src` are the source parameters `SourceParameters`
- `path` are the path parameters `PathParameters`
- `site` are the site parameters `SiteParameters`
"""
fourier_spectral_ordinate(f::U, m::S, r_ps::T, src::SourceParameters, path::PathParameters, site::SiteParameters) where {S<:Real,T<:Real,U<:Float64} = fourier_spectral_ordinate(f, m, r_ps, src, path.geometric, path.saturation, path.anelastic, site)

"""
	fourier_spectral_ordinate(f::U, m::S, r_ps::T, fas::FourierParameters) where {S<:Real,T<:Real,U<:Float64}

Fourier acceleration spectral ordinate (m/s) based upon an equivalent point source distance `r_ps`
- `f` is frequency (Hz)
- `m` is magnitude
- `r_ps` is the equivalent point source distance including saturation effects (km)
- `fas` are the Fourier spectral parameters `FourierParameters`
"""
fourier_spectral_ordinate(f::U, m::S, r_ps::T, fas::FourierParameters) where {S<:Real,T<:Real,U<:Float64} = fourier_spectral_ordinate(f, m, r_ps, fas.source, fas.path, fas.site)

"""
	squared_fourier_spectral_ordinate(f::U, m::S, r_ps::T, src::SourceParameters, geo::GeometricSpreadingParameters, ane::AnelasticAttenuationParameters, site::SiteParameters) where {S<:Real,T<:Real,U<:Float64}

Squared Fourier acceleration spectral ordinate (m^2/s^2) based upon an equivalent point source distance `r_ps`
- `f` is frequency (Hz)
- `m` is magnitude
- `r_ps` is the equivalent point source distance including saturation effects (km)
- `src` are the source parameters `SourceParameters`
- `geo` are the geometric spreading parameters `GeometricSpreadingParameters`
- `sat` are the near source saturation parameters `NearSourceSaturationParameters`
- `ane` are the anelastic attenuation parameters `AnelasticAttenuationParameters`
- `site` are the site parameters `SiteParameters`

See also: [`fourier_spectral_ordinate`](@ref), [`squared_fourier_spectrum`](@ref)
"""
function squared_fourier_spectral_ordinate(f::S, m::S, r_ps::T, src::SourceParameters, geo::GeometricSpreadingParameters, sat::NearSourceSaturationParameters, ane::AnelasticAttenuationParameters, site::SiteParameters) where {S<:Real,T<:Real,U<:Float64}
	return fourier_spectral_ordinate(f, m, r_ps, src, geo, sat, ane, site)^2
end

squared_fourier_spectral_ordinate(f::U, m::S, r_ps::T, src::SourceParameters, path::PathParameters, site::SiteParameters) where {S<:Real,T<:Real,U<:Float64} = squared_fourier_spectral_ordinate(f, m, r_ps, src, path.geometric, path.saturation, path.anelastic, site)
squared_fourier_spectral_ordinate(f::U, m::S, r_ps::T, fas::FourierParameters) where {S<:Real,T<:Real,U<:Float64} = squared_fourier_spectral_ordinate(f, m, r_ps, fas.source, fas.path, fas.site)

"""
	fourier_spectrum(f::Vector{U}, m::S, r_ps::T, fas::FourierParameters) where {S<:Real,T<:Real,U<:Float64}

Fourier acceleration spectrum (m/s) based upon an equivalent point source distance `r_ps`
- `f` is `Vector` of frequencies (Hz)
- `m` is magnitude
- `r_ps` is the equivalent point source distance including saturation effects (km)
- `fas` are the Fourier spectral parameters `FourierParameters`

See also: [`fourier_spectral_ordinate`](@ref)
"""
function fourier_spectrum(f::Vector{U}, m::S, r_ps::T, fas::FourierParameters) where {S<:Real,T<:Real,U<:Float64}
	numf = length(f)
	V = get_parametric_type(fas)
	if V <: Dual
		W = V
	elseif S <: Dual
		W = S
	else
		W = U
	end
	if numf == 0
		return Vector{W}()
	else
		# define all frequency independent terms here
	  	C = fourier_constant(fas)
		# source amplitude and corner frequencies
		Mo = magnitude_to_moment(m)
		fa, fb, ε = corner_frequency(m, fas)
		# geometric spreading
		Gr = geometric_spreading(r_ps, m, fas)
		factor = 4π^2 * C * Mo * Gr / 100.0
		if fas.path.anelastic.rmetric == :Rrup
			r_rup = rupture_distance_from_equivalent_point_source_distance(r_ps, m, fas)
		end

		Af = Vector{W}(undef, numf)
	  	for i in 1:numf
			@inbounds fi = f[i]
		  	# source term
			Ef = fourier_source_shape(fi, fa, fb, ε, fas.source.model)
	  		# combined attenuation
			if fas.path.anelastic.rmetric == :Rrup
				Kf = fourier_attenuation(fi, r_rup, fas)
			else
				Kf = fourier_attenuation(fi, r_ps, fas)
			end
	  		# site impedance
	      	Sf = site_amplification(fi, fas)
			# apply factor and convert to acceleration in appropriate units
	  		@inbounds Af[i] = Ef * Kf * Sf * factor * fi^2
	    end
	    return Af
	end
end

"""
	squared_fourier_spectrum(f::Vector{U}, m::S, r_ps::T, fas::FourierParameters) where {S<:Real,T<:Real,U<:Float64}

Squared Fourier amplitude spectrum. Useful within spectral moment calculations.

See also: [`fourier_spectrum`](@ref), [`squared_fourier_spectral_ordinate`](@ref)
"""
function squared_fourier_spectrum(f::Vector{U}, m::S, r_ps::T, fas::FourierParameters) where {S<:Real,T<:Real,U<:Float64}
	Af = fourier_spectrum(f, m, r_ps, fas)
	return Af.^2
end

"""
	fourier_spectrum!(Af::Vector{U}, f::Vector{V}, m::S, r_ps::T, fas::FourierParameters) where {S<:Real,T<:Real,U<:Real,V<:Float64}

Fourier acceleration spectrum (m/s) based upon an equivalent point source distance `r_ps`
- `Af` is the vector of fas amplitudes to be filled (m/s)
- `f` is `Vector` of frequencies (Hz)
- `m` is magnitude
- `r_ps` is the equivalent point source distance including saturation effects (km)
- `fas` are the Fourier spectral parameters `FourierParameters`

See also: [`fourier_spectrum`](@ref)
"""
function fourier_spectrum!(Af::Vector{U}, f::Vector{V}, m::S, r_ps::T, fas::FourierParameters) where {S<:Real,T<:Real,U<:Real,V<:Float64}
	numf = length(f)
	numAf = length(Af)
	numf == numAf || error("length of `f` and `Af` must be equal")
	if numf == 0
		Af = Vector{U}()
	else
		# define all frequency independent terms here
	  	C = fourier_constant(fas)
		# source amplitude and corner frequencies
		Mo = magnitude_to_moment(m)
		fa, fb, ε = corner_frequency(m, fas)
		# geometric spreading
		Gr = geometric_spreading(r_ps, m, fas)
		factor = 4π^2 * C * Mo * Gr / 100.0
		if fas.path.anelastic.rmetric == :Rrup
			r_rup = rupture_distance_from_equivalent_point_source_distance(r_ps, m, fas)
		end

		for i in 1:numf
			@inbounds fi = f[i]
			# source term
			Ef = fourier_source_shape(fi, fa, fb, ε, fas.source.model)
	  		# combined attenuation
			if fas.path.anelastic.rmetric == :Rrup
				Kf = fourier_attenuation(fi, r_rup, fas)
			else
				Kf = fourier_attenuation(fi, r_ps, fas)
			end
	  		# site impedance
	      	Sf = site_amplification(fi, fas)
			# apply factor and convert to acceleration in appropriate units
	  		@inbounds Af[i] = Ef * Kf * Sf * factor * fi^2
	  	end
	end
	return nothing
end

"""
	squared_fourier_spectrum!(Afsq::Vector{U}, f::Vector{V}, m::S, r_ps::T, fas::FourierParameters) where {S<:Real,T<:Real,U<:Real,V<:Float64}

Fourier acceleration spectrum (m/s) based upon an equivalent point source distance `r_ps`
- `Afsq` is the vector of squared fas amplitudes to be filled (m^2/s^2)
- `f` is `Vector` of frequencies (Hz)
- `m` is magnitude
- `r_ps` is the equivalent point source distance including saturation effects (km)
- `fas` are the Fourier spectral parameters `FourierParameters`

See also: [`squared_fourier_spectrum`](@ref), ['squared_fourier_spectral_ordinate'](@ref)
"""
function squared_fourier_spectrum!(Afsq::Vector{U}, f::Vector{V}, m::S, r_ps::T, fas::FourierParameters) where {S<:Real,T<:Real,U<:Real,V<:Float64}
	numf = length(f)
	numAfsq = length(Afsq)
	numf == numAfsq || error("length of `f` and `Afsq` must be equal")
	if numf == 0
		Af = Vector{U}()
	else
		# define all frequency independent terms here
	  	C = fourier_constant(fas)
		# source amplitude and corner frequencies
		Mo = magnitude_to_moment(m)
		fa, fb, ε = corner_frequency(m, fas)
		# geometric spreading
		Gr = geometric_spreading(r_ps, m, fas)
		factor = (4π^2 * C * Mo * Gr / 100.0)^2
		if fas.path.anelastic.rmetric == :Rrup
			r_rup = rupture_distance_from_equivalent_point_source_distance(r_ps, m, fas)
		end

		for i in 1:numf
			@inbounds fi = f[i]
			# source term
			Ef = fourier_source_shape(fi, fa, fb, ε, fas.source.model)
	  		# combined attenuation
			if fas.path.anelastic.rmetric == :Rrup
				Kf = fourier_attenuation(fi, r_rup, fas)
			else
				Kf = fourier_attenuation(fi, r_ps, fas)
			end
	  		# site impedance
	      	Sf = site_amplification(fi, fas)
			# apply factor and convert to acceleration in appropriate units
	  		@inbounds Afsq[i] = (Ef * Kf * Sf)^2 * factor * fi^4
	  	end
	end
	return nothing
end

"""
	combined_kappa_frequency(r::T, Af2target::Float64, ane::AnelasticAttenuationParameters, site::SiteParameters) where T<:Real

Frequency at which the combined κ_r and κ_0 filters (squared versions) give a value of `Af2target`.
`r` can be either `r_ps` or `r_rup` depending upon what matches `ane.rmetric`
"""
function combined_kappa_frequency(r::T, Af2target::Float64, ane::AnelasticAttenuationParameters, site::SiteParameters) where T<:Real
	Q0_eff, η_eff, cQ_eff = effective_quality_parameters(r, ane)
  	if η_eff < 0.1
    	# a closed form solution exists (for effectively η=0)
    	return log(1.0/Af2target) / ( 2π * ( r / (Q0_eff * cQ_eff) + site.κ0 ) )
  	else
		g(f) = Af2target - (fourier_attenuation(f, r, ane, site)^2)
		f_0 = find_zero(g, (0.01,100.0), Bisection(); xatol=1e-2)
		# fk = min(max(f_0, 0.2), 1.0)
		fk = max(f_0, 0.1)
		U = get_parametric_type(ane)
		V = get_parametric_type(site)
		if T <: Float64
			if U <: Float64
				unit = oneunit(V)
			else
				unit = oneunit(U)
			end
		else
			unit = oneunit(T)
		end
	    return fk * unit
  	end
end

combined_kappa_frequency(r::T, Af2target::Float64, path::PathParameters, site::SiteParameters) where T<:Real = combined_kappa_frequency(r, Af2target, path.anelastic, site)
combined_kappa_frequency(r::T, Af2target::Float64, fas::FourierParameters) where T<:Real = combined_kappa_frequency(r, Af2target, fas.path, fas.site)


text:
using Documenter
using BulkSMS

makedocs(
    format = :html,
    sitename = "BulkSMS.jl",
    pages = [
        "index.md",
    ]
)

deploydocs(
    repo = "github.com/scls19fr/BulkSMS.jl.git"
)

ext: jl


text:
#cd("/home/user/github/TOF-Tracer2/")

include("$(pwd())/startupAPiTOF.jl")

# fp = "./ExampleFiles/APiTOFDATA/" # All files in this path will be processed
#fp = "/media/wiebke/Extreme SSD/PSM_vs_PTR3/Data/apiTOFdata/bCP/"

# fp = "/media/wiebke/Extreme SSD/PSM_vs_PTR3/Data/apiTOFdata/bCP_14-07/1stFile/"
fp = "/media/wiebke/Extreme SSD/PSM_vs_PTR3/Data/apiTOFdata/bCP_13-07/"
# fp = "/media/wiebke/Extreme SSD/PSM_vs_PTR3/Data/apiTOFdata/NO3-/bCP/"

filefilterRegexp = r"\.h5$"
#rf = "$(fp)Data_2019_07_14-15_19_40.h5"  # The mass scale from this file defines the mass scale of all
# rf = "$(fp)Data_2019_07_14-11_19_39.h5"		# 14.7. first file
rf = "$(fp)Data_2019_07_13-17_34_29.h5"	# 13.7.
# rf = "$(fp)Data_2019_07_29-13_25_54.h5" # NO3-/bCP/

#masslist = MasslistFunctions.loadMasslist("./ExampleFiles/MASSLISTS/exampleMasslist.csv")
#masslist = MasslistFunctions.loadMasslist("/media/wiebke/Extreme SSD/PSM_vs_PTR3/Data/apiTOFdata/masslists/bCP_BAG_APItof.csv")
masslist = MasslistFunctions.loadMasslist("/media/wiebke/Extreme SSD/PSM_vs_PTR3/Data/apiTOFdata/bCP_13-07/masslist-JAN2021_2.csv")
# masslist = MasslistFunctions.loadMasslist("$(fp)results/Masslist.csv") # no3-CIMS
# cr = [37.028406 487.3424 693.546 739.515 959.697] # pos bCP, 14.-16.7.
cr = [37.028406 252.17256 634.344]	# pos bCP, 13.7
#cr = [37.028406 252.17256]
#cr = [62 125] # neg
# cr = [61.9878 187.979] # NO3--CIMS

# alternatively: use an auto generated masslist
# masslistMasses, masslistElements, masslistElementsMasses, masslistCompositions = MasslistFunctions.createMassList(;C=13:120, O=0:20, N=0:1, allowRadicals=true) #
(masslistMasses, masslistElements, masslistElementsMasses, masslistCompositions) = masslist
s = (masslistMasses.>17) .& ( masslistMasses.<1200)
masslistMasses = masslistMasses[s]
masslistCompositions = masslistCompositions[s,:]

####################### END OF SETTINGS ###############################################################

####################### Processing sequence ###########################################################

correctMassScaleAndExtractSumSpec(
    fp,
    masslistMasses,
    masslistElements,
    masslistElementsMasses,
    masslistCompositions,
    rf,
    cr,
    filefilterRegexp=filefilterRegexp,
    onlyUseAverages = false,
    plotControlMass = false,
    testRangeStart = 393.5, # the mass shift of this region will be shown if plot control mass is set true. Should not be part of calibRegions
    testRangeEnd = 396.5,
    recalibInterval = 900,
    resolution = 4500,
    firstNFiles= 0,
    lastNFiles = 0,
    UMRpeaks = true, # create Unit Mass Resolution peaks
    UMRmasses = length(s),
    UMRmassLow = 0.1,
    UMRmassHigh = 0.7,
    umrBinWidth =100,
    #binWidth = 10
    )
#
baselineAndPeakshape(
    fp,
    peakshapeRegions=6,
    peakshapeRegionStretch=1.0,
    peakshapeQuantileValue = 0.7,
    peakfindingNoiseThresholdValue = 5,
    peakWindowWidth = 100,
    peakfindingSignalLimit=0.5
    )
    #
#
mtrx = deconvolute(
    fp,
    calcTransposed = false
    )
#


text:
### Edge and Gradient Related Image Operations ###

# Edge/gradient filters

function sobel()
    f = [ -1.0  0.0  1.0
          -2.0  0.0  2.0
          -1.0  0.0  1.0 ]
    return f', f
end

function prewitt()
    f = [ -1.0  0.0  1.0
          -1.0  0.0  1.0
          -1.0  0.0  1.0 ]
    return f', f
end

# Consistent Gradient Operators
# Ando Shigeru
# IEEE Trans. Pat. Anal. Mach. Int., vol. 22 no 3, March 2000
#
# TODO: These coefficients were taken from the paper It would be nice
#       to resolve the optimization problem and use higher precision
#       versions, which might allow better separable approximations of
#       ando4 and ando5.

function ando3()
    f = [ -0.112737  0.0  0.112737
          -0.274526  0.0  0.274526
          -0.112737  0.0  0.112737 ]
    return f', f
end

# Below, the ando4() and ando5() functions return filters with
# the published filter values.  The ando4_sep() and ando5_sep()
# functions return separable approximations to the corresponding
# filters, estimated using the projection of the actual values on the
# eigenvector corresponding to the largest eigenvalue of the SVD of
# the original filter.

function ando4()
    f = [ -0.022116 -0.025526  0.025526  0.022116
          -0.098381 -0.112984  0.112984  0.098381
          -0.098381 -0.112984  0.112984  0.098381
          -0.022116 -0.025526  0.025526  0.022116 ]
    return f', f
end

function ando4_sep()
    f = [-0.022175974729759376 -0.025473821998749126 0.025473821998749126 0.022175974729759376
         -0.09836750569692418  -0.11299599504060115  0.11299599504060115  0.09836750569692418
         -0.09836750569692418  -0.11299599504060115  0.11299599504060115  0.09836750569692418
         -0.022175974729759376 -0.025473821998749126 0.025473821998749126 0.022175974729759376]
    return f', f
end

function ando5()
    f = [ -0.003776 -0.010199  0.0  0.010199  0.003776
          -0.026786 -0.070844  0.0  0.070844  0.026786
          -0.046548 -0.122572  0.0  0.122572  0.046548
          -0.026786 -0.070844  0.0  0.070844  0.026786
          -0.003776 -0.010199  0.0  0.010199  0.003776 ]
    return f', f
end

function ando5_sep()
    f = [-0.0038543900766123762 -0.0101692999709622   0.0  0.0101692999709622   0.0038543900766123762
         -0.026843218687756566  -0.07082229291692607  0.0  0.07082229291692607  0.026843218687756566
         -0.046468878396946627  -0.12260200818803602  0.0  0.12260200818803602  0.046468878396946627
         -0.026843218687756566  -0.07082229291692607  0.0  0.07082229291692607  0.026843218687756566
         -0.0038543900766123762 -0.0101692999709622   0.0  0.0101692999709622   0.0038543900766123762]
    return f', f
end

# Image gradients in the X and Y direction
function imgradients(img::AbstractArray, method::String="ando3", border::String="replicate")
    sx,sy = spatialorder(img)[1] == "x" ? (1,2) : (2,1)
    s = (method == "sobel"     ? sobel() :
         method == "prewitt"   ? prewitt() :
         method == "ando3"     ? ando3() :
         method == "ando4"     ? ando4() :
         method == "ando5"     ? ando5() :
         method == "ando4_sep" ? ando4_sep() :
         method == "ando5_sep" ? ando5_sep() :
         error("Unknown gradient method: $method"))

    grad_x = imfilter(img, s[sx], border)
    grad_y = imfilter(img, s[sy], border)

    return grad_x, grad_y
end

# Magnitude of gradient, calculated from X and Y image gradients
magnitude(grad_x::AbstractArray, grad_y::AbstractArray) = hypot(grad_x, grad_y)

# Phase (angle of steepest gradient ascent), calculated from X and Y gradient images
function phase{T}(grad_x::AbstractArray{T}, grad_y::AbstractArray{T})
    EPS = sqrt(eps(eltype(T)))
    # Set phase to zero when both gradients are close to zero
    reshape([atan2(-grad_y[i], grad_x[i]) * ((abs(grad_x[i]) > EPS) | (abs(grad_y[i]) > EPS))
             for i=1:length(grad_x)], size(grad_x))
end

function phase(grad_x::AbstractImageDirect, grad_y::AbstractImageDirect)
    img = copy(grad_x, phase(data(grad_x), data(grad_y)))
    img["limits"] = (-float(pi),float(pi))
    img
end

# Orientation of the strongest edge at a point, calculated from X and Y gradient images
# Note that this is perpendicular to the phase at that point, except where
# both gradients are close to zero.

function orientation{T}(grad_x::AbstractArray{T}, grad_y::AbstractArray{T})
    EPS = sqrt(eps(eltype(T)))
    # Set orientation to zero when both gradients are close to zero
    # (grad_y[i] should probably be negated here, but isn't for consistency with earlier releases)
    reshape([atan2(grad_x[i], grad_y[i]) * ((abs(grad_x[i]) > EPS) | (abs(grad_y[i]) > EPS))
             for i=1:length(grad_x)], size(grad_x))
end

function orientation(grad_x::AbstractImageDirect, grad_y::AbstractImageDirect)
    img = copy(grad_x, orientation(data(grad_x), data(grad_y)))
    img["limits"] = (-float(pi),float(pi))
    img
end

# Return both the magnitude and phase in one call
magnitude_phase(grad_x::AbstractArray, grad_y::AbstractArray) = (magnitude(grad_x,grad_y), phase(grad_x,grad_y))

# Return the magnituded and phase of the gradients in an image
function magnitude_phase(img::AbstractArray, method::String="ando3", border::String="replicate")
    grad_x, grad_y = imgradients(img, method, border)
    return magnitude_phase(grad_x, grad_y)
end

# Return the x-y gradients and magnitude and phase of gradients in an image
function imedge(img::AbstractArray, method::String="ando3", border::String="replicate")
    grad_x, grad_y = imgradients(img, method, border)
    mag = magnitude(grad_x, grad_y)
    orient = orientation(grad_x, grad_y)
    return (grad_x, grad_y, mag, orient)
end

# Thin edges
thin_edges{T}(img::AbstractArray{T,2}, gradientangles::AbstractArray, border::String="replicate") =
    thin_edges_nonmaxsup(img, gradientangles, border)
thin_edges_subpix{T}(img::AbstractArray{T,2}, gradientangles::AbstractArray, border::String="replicate") =
    thin_edges_nonmaxsup_subpix(img, gradientangles, border)

# Code below is related to non-maximal suppression, and was ported to Julia from
# http://www.csse.uwa.edu.au/~pk/research/matlabfns/Spatial/nonmaxsup.m
# (Please conserve the original copyright below.)

# NONMAXSUP - Non-maxima suppression
#
# Usage:
#          (im,location) = nonmaxsup(img, gradientangles, radius);
#
# Function for performing non-maxima suppression on an image using
# gradient angles.  Gradient angles are assumed to be in radians.
#
# Input:
#   img - image to be non-maxima suppressed.
#
#   gradientangles - image containing gradient angles around each pixel in radians
#                    (-pi,pi)
#
#   radius  - Distance in pixel units to be looked at on each side of each
#             pixel when determining whether it is a local maxima or not.
#             This value cannot be less than 1.
#             (Suggested value about 1.2 - 1.5)
#
# Returns:
#   im        - Non maximally suppressed image.
#   location  - `Graphics.Point` image holding subpixel locations of edge
#               points.
#
# Notes:
#
# This function uses bilinear interpolation to estimate
# intensity values at ideal, real-valued pixel locations on each side of
# pixels to determine if they are local maxima.

# Copyright (c) 1996-2013 Peter Kovesi
# Centre for Exploration Targeting
# The University of Western Australia
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# The Software is provided "as is", without warranty of any kind.

# December  1996 - Original version
# September 2004 - Subpixel localization added
# August    2005 - Made Octave compatible
# October   2013 - Final thinning applied to binary image for Octave
#                  compatbility (Thanks to Chris Pudney)
# June      2014 - Ported (and modified significantly) to Julia (Kevin Squire)

import Base.Graphics.Point

if !applicable(zero, Point)
    import Base.zero
    zero(Point) = Point(0.0,0.0)
end

# Used to encode the sign, integral, and fractional components of
# an offset from a coordinate
immutable CoordOffset
    s::Int      # sign
    i::Int      # integer part
    f::Float64  # fractional part
end

CoordOffset(x::Float64) = ((frac,i) = modf(x); CoordOffset(sign(frac), int(i), abs(frac)))
(-)(off::CoordOffset) = CoordOffset(-off.s,-off.i, off.f)
(*)(x::Number, off::CoordOffset) = x*(off.i + off.s*off.f)
(*)(off::CoordOffset, x::Number) = x*(off.i + off.s*off.f)
(+)(x::Number, off::CoordOffset) = x + off.i + off.s*off.f
(+)(off::CoordOffset, x::Number) = x + off.i + off.s*off.f

# Precalculate x and y offsets relative to centre pixel for each orientation angle
function _calc_discrete_offsets(θ, radius, transposed)

    θ_count = iround(2π/θ)
    θ = 2π/θ_count
    angles = (0:θ_count)*θ

    # x and y offset of points at specified radius and angles
    # from each reference position.

    if transposed
        # θ′ = -π/2 - θ
        xoffs = [CoordOffset(-x) for x in  radius*sin(angles)]
        yoffs = [CoordOffset( y) for y in  radius*cos(angles)]
    else
        xoffs = [CoordOffset( x) for x in  radius*cos(angles)]
        yoffs = [CoordOffset(-y) for y in  radius*sin(angles)]
    end

    return θ, xoffs, yoffs
end

_discretize_angle(angle::FloatingPoint, invθ) =
    angle < 0 ? iround((angle + 2π)*invθ)+1 : iround(angle*invθ)+1

# Interpolate the value of an offset from a particular pixel
#
# Returns (interpolated value, min_value of adjacent pixels in direction of offset)
#
# The second value is made available to eliminate double edges; if the value at
# (x,y) is less than the value or values adjacent to it in the direction of the
# gradient (xoff,yoff), then it is not a local maximum

function _interp_offset(img::AbstractArray, x::Integer, y::Integer, xoff::CoordOffset, yoff::CoordOffset, Ix, Iy, pad)
    fx = Ix[x + xoff.i + pad]
    fy = Iy[y + yoff.i + pad]
    cx = Ix[x + xoff.i + xoff.s + pad]
    cy = Iy[y + yoff.i + yoff.s + pad]

    tl = img[fy,fx]    # Value at bottom left integer pixel location.
    tr = img[fy,cx]    # bottom right
    bl = img[cy,fx]    # top left
    br = img[cy,cx]    # top right

    upperavg = tl + xoff.f * (tr - tl)  # Now use bilinear interpolation to
    loweravg = bl + xoff.f * (br - bl)  # estimate value at x,y

    min_adjacent = (fx == x) & (fy == y) ? min(tr,bl) : tl

    return (upperavg + yoff.f * (loweravg - upperavg), min_adjacent)
end

# Core edge thinning algorithm using nonmaximal suppression
function thin_edges_nonmaxsup_core!{T}(out::AbstractArray{T,2}, location::AbstractArray{Point,2},
                                       img::AbstractArray{T,2}, gradientangles, radius, border, theta)
    calc_subpixel = !isempty(location)

    # Error checking
    size(img) == size(gradientangles) == size(out) || error("image, gradient angle, and output image must all be the same size")
    calc_subpixel && size(location) != size(img) && error("subpixel location has a different size than the input image")
    radius < 1.0 && error("radius must be >= 1")

    # Precalculate x and y offsets relative to centre pixel for each orientation angle
    transposed = spatialorder(img)[1] == "x"
    θ, xoffs, yoffs = _calc_discrete_offsets(theta, radius, transposed)
    iθ = 1/θ

    # Indexes to use for border handling
    pad = iceil(radius)
    Ix = Images.padindexes(img, 2, pad, pad, border)
    Iy = Images.padindexes(img, 1, pad, pad, border)

    # Now run through the image interpolating grey values on each side
    # of the centre pixel to be used for the non-maximal suppression.

    (height,width) = size(img)

    for x = 1:width, y = 1:height
        (c = img[y,x]) == 0 && continue  # For thresholded images

        or = _discretize_angle(gradientangles[y,x],iθ)   # Disretized orientation
        v1, n1 = _interp_offset(img, x, y, xoffs[or], yoffs[or], Ix, Iy, pad)

        if (c > v1) & (c >= n1) # We need to check the value on the other side...
            v2, n2 = _interp_offset(img, x, y, -xoffs[or], -yoffs[or], Ix, Iy, pad)

            if (c > v2) & (c >= n2)  # This is a local maximum.
                                     # Record value in the output image.
                if calc_subpixel
                    # Solve for coefficients of parabola that passes through
                    # [-1, v2]  [0, img] and [1, v1].
                    # v = a*r^2 + b*r + c

                    # c = img[y,x]
                    a = (v1 + v2)/2 - c
                    b = a + c - v2

                    # location where maxima of fitted parabola occurs
                    r = -b/2a
                    location[y,x] = transposed ? Point(y + r*yoffs[or], x + r*xoffs[or]) :
                                                 Point(x + r*xoffs[or], y + r*yoffs[or])

                    if T<:FloatingPoint
                        # Store the interpolated value
                        out[y,x] = a*r^2 + b*r + c
                    else
                        out[y,x] = c
                    end
                else
                    out[y,x] = c
                end
            end
        end
    end

    out
end

# Main function call when subpixel location of edges is not needed
function thin_edges_nonmaxsup!{A<:AbstractArray,B<:AbstractArray}(out::A, img::A, gradientangles::B, border::String="replicate";
                                                                  radius::Float64=1.35, theta=pi/180)
    properties(out) != properties(img) && error("Input and output arrays must have the same properties")
    thin_edges_nonmaxsup_core!(data(out), Array(Point,(0,0)), img, gradientangles, radius, border, theta)
    out
end

function thin_edges_nonmaxsup{T}(img::AbstractArray{T,2}, gradientangles::AbstractArray, border::String="replicate";
                                 radius::Float64=1.35, theta=pi/180)
    (height,width) = size(img)
    out = zeros(T, height, width)
    thin_edges_nonmaxsup_core!(out, Array(Point,(0,0)), img, gradientangles, radius, border, theta)
    copy(img, out)
end

# Main function call when subpixel location of edges is desired
function thin_edges_nonmaxsup_subpix!{A<:AbstractArray, B<:AbstractArray, C<:AbstractArray}(out::A, location::B, img::A, gradientangles::C,
                                     border::String="replicate"; radius::Float64=1.35, theta=pi/180)
    properties(out) != properties(img) && error("Input and output arrays must have the same properties")
    eltype(location) != Point && error("Preallocated subpixel location array/image must have element type Graphics.Point")

    thin_edges_nonmaxsup_core!(data(out), data(location), img, gradientangles, radius, border, theta)
    img, location
end

function thin_edges_nonmaxsup_subpix{T}(img::AbstractArray{T}, gradientangles::AbstractArray,
                                        border::String="replicate"; radius::Float64=1.35, theta=pi/180)
    (height,width) = size(img)
    out = zeros(T, height, width)
    location = zeros(Point, height, width)
    thin_edges_nonmaxsup_core!(out, location, img, gradientangles, radius, border, theta)

    copy(img, out), copy(img, location)
end


text:
###Measurement dataset for the FitzHugh-Nagumo model
### Column 1: Membrane Potential
### Column 2: Recovery Variable
function fitzhughnagumodataset1()
    #create the data
    timepoints = linspace(0.0,20.0,100)
    measurements =
        [-0.921082915224310	1.05807254284326;
         -0.500535384521994	0.998164068376517;
         -0.768325064924386	1.10680611674722;
         0.255419084728969	1.07302916025702;
         1.12280059296848	1.09162325946439;
         1.66551088744314	0.926274890184664;
         1.97284465158118	0.869719057122973;
         2.07855316641910	0.752860189269569;
         2.51598964168118	0.689590370696254;
         2.35479975777723	0.431865384763390;
         1.70613356517142	0.176083501497936;
         2.30450527630536	0.149013019696289;
         1.91907950350955	0.188445901705940;
         1.75555328728763	-0.0867899360211768;
         1.82036449971724	-0.0491487726541759;
         1.63422125572360	-0.205098633238317;
         1.59284756728848	-0.207705813006382;
         1.77399661344685	-0.530915138581783;
         1.70346553016991	-0.492490677428846;
         1.64231131404376	-0.640553246555571;
         1.46595478617355	-0.429028904022557;
         1.11738627815334	-0.640211608888171;
         1.31928638894147	-0.663184131319016;
         1.36187817721493	-0.886284832848348;
         1.08712200870823	-0.892889498373029;
         1.03531619730304	-0.918113366624786;
         0.816419611670972	-0.851546780148877;
         0.413322482483550	-0.960495284497250;
         0.0893474521057900	-0.884920203434804;
         -0.781753401696825	-1.03052396959950;
         -1.39939458536829	-0.813280569722894;
         -2.07881606398160	-0.703639111094413;
         -2.10792570250433	-0.603675493940061;
         -2.03308033114145	-0.310170074385900;
         -2.29462712665774	-0.181446208278314;
         -1.59498462780282	-0.0598454144616705;
         -1.69895988929073	-0.0212779721991443;
         -1.79491682075546	0.275980915725031;
         -1.41668329284817	0.343442623375234;
         -1.79816273939288	0.413551182167112;
         -1.48536900468108	0.544310465847005;
         -1.42171923369266	0.626321374333452;
         -1.24523884390136	0.617872545135805;
         -1.13765992644468	0.805456961068046;
         -1.18037478939099	0.845038980651403;
         -0.892424006196079	0.901298709234302;
         -0.683180624613448	0.941669548300683;
         -0.209026714126607	1.01752211945379;
         0.490557914203931	0.916343766410746;
         1.46133449146722	1.06690514510804;
         1.78903121067960	0.924702975861434;
         2.03577383530041	0.755632602347172;
         1.82707695269093	0.622861434073005;
         1.80250538832746	0.442271009186925;
         1.92175719502324	0.444481551519144;
         2.10304938083406	0.246133961520414;
         1.71922477276955	0.0918338337045172;
         1.83949186415677	0.125105477330890;
         1.70319184432776	-0.0886301066507938;
         1.84994155095814	-0.213905983891063;
         1.47376581265411	-0.288833963613749;
         1.58376241695376	-0.417610633824004;
         1.60269934897514	-0.521947994160419;
         1.62242302825940	-0.350425779810062;
         1.62266660256747	-0.485746419666361;
         1.33849205873554	-0.648547299607231;
         1.02991634748470	-0.820478569127451;
         1.05347385926553	-0.850919872222649;
         0.906784172804639	-0.853847425452138;
         1.28746854395457	-0.829510375043008;
         0.698841106938930	-1.00917940847391;
         0.686326820818500	-1.09859162505725;
         0.217983199788489	-1.04021467327920;
         -0.193867396744838	-0.890219994291912;
         -1.30880509976363	-0.810624062099305;
         -2.02327338556172	-0.677672145225874;
         -2.16047832471496	-0.567329653518438;
         -1.86124938568770	-0.394550555341886;
         -1.91096822989846	-0.293576567693192;
         -1.85887732018248	-0.0599748773654953;
         -1.56359949530767	-0.0781842605789201;
         -1.66791797726924	0.177009282023116;
         -1.61702322794403	0.210628369060525;
         -1.34419493206595	0.363731689674640;
         -1.62128634877271	0.536324589690918;
         -1.32026468435208	0.674527970437903;
         -1.21038513194256	0.623010943616307;
         -1.26686385215658	0.877194445843865;
         -1.07952532445313	0.916954470065889;
         -1.13397562899673	0.937921646353368;
         -0.933989145160547	0.988039291007753;
         -0.458145891632811	1.02862284526891;
         0.125712026070199	1.03886478304373;
         1.25725486092870	1.03269416125531;
         1.64694338282048	0.943683614618281;
         2.03647161336170	0.869785318820964;
         2.00596683034952	0.785529247325093;
         1.69948310060718	0.582659095081282;
         1.87663874931099	0.410658466868349;
         1.63352834363561	0.347679080252495]

    #return the data
    timepoints,measurements
end


text:
# beamform 2d util

function lastindmax(x)
   k = 1
   m = x[1]
   @inbounds for i in eachindex(x)
       if x[i]>=m
           k = i
           m = x[i]
       end
   end
   return k
end

function findMaxEnvelope(timeseries::Vector{Float64};
                         windowlen::Int=2^16,
                         fs::Int=192000,
                         tukeyalpha::Float64=0.7  )
  #
  window = tukey(windowlen, tukeyalpha)
  convres = conv(abs.(timeseries), window)
  stw = lastindmax(convres) - 2^16
  spw = stw + 2^16-1
  return stw, spw
end

function loadPs3WavEnvelopes(wavfilepath::String; variable="wavdata", cutmax::Bool=true)
  wavdata = load(wavfilepath, variable)
  if cutmax
    stw, spw = findMaxEnvelope(wavdata[0])
    for (k,va) in wavdata
      wavdata[k] = wavdata[k][stw:spw]
    end
  end
  return wavdata
end

function beamformBasic(csvWaveData::Array{Float64,2}, chirpIn; soundSpeed=330.0)::Vector{Float64}

  fFloor = 3950;
  fCeil = 8050;
  fSampling = 192000;
  nFFT_full = 2^16;
  nFFT_czt = 2^14;
  azimuthDivs = 180;
  # rngSNRfloor = cfgd["range_snr_floor"];
  # aziSNRfloor = cfgd["azi_snr_floor"];

  #Params subject to change/ Heuristics
  totalPhones = 4
  azimuths = LinRange(0,360,azimuthDivs)*pi/180;

  #Data preprocessing MF + CZT
  w = exp(-2im*pi*(fCeil-fFloor)/(nFFT_czt*fSampling))
  a = exp(2im*pi*fFloor/fSampling)

  arrayPos = collect([0.0 0.0 0.0 0.0; 0.03 0.01 -0.01 -0.03]')

  # Matched Filter on Data In
  nFFT_full = nextpow(2,size(csvWaveData,1))  # MF
  mfData = zeros(Complex{Float64}, nFFT_full, totalPhones)
  mf = prepMF(chirpIn, nFFT_full,totalPhones) # MF
  mf(csvWaveData, mfData) # MF
  cztData = zeros(Complex{Float64}, nFFT_czt, totalPhones)
  filterCZT = prepCZTFilter(nFFT_full, totalPhones, w, nFFT_czt,a)
  filterCZT(mfData, cztData)

  # CBF Filter shared by
  FFTfreqs = collect(LinRange(fFloor,fCeil,nFFT_czt))
  cfgTotal = CBFFilterConfig(fFloor,fCeil,nFFT_czt,totalPhones,azimuths,soundSpeed, FFTfreqs)

  sourceXY = zeros(2)
  BFTemp = zeros(Complex{Float64}, totalPhones, nFFT_czt)
  CBFFull = zeros(Complex{Float64}, getCBFFilter2Dsize(cfgTotal))
  BFOutFull = zeros(length(azimuths))

  constructCBFFilter2D!(cfgTotal, arrayPos, CBFFull, sourceXY, BFTemp)

  temp1 = zeros(Complex{Float64}, nFFT_czt)
  temp2 = zeros(Complex{Float64}, nFFT_czt, totalPhones)

  CBF2D_DelaySum!(cfgTotal, cztData, BFOutFull, temp1, temp2, CBFFull)

  intensity = norm.(BFOutFull)
  intensity ./= sum(intensity)

  return intensity
end

function intensityToCircularKDE(intensity::Vector{Float64}; N=100, SNRfloor::Float64=0.8)
  # increase resolution for aliasing sampler
  azi_weights = (repeat(norm.(intensity), 1,8)')[:]  # BFOutFull
  azi_weights ./= sum(azi_weights)

  azi_domain = collect(LinRange(0,2pi-1e-10,length(azi_weights)))

  bss = AliasingScalarSampler(azi_domain, azi_weights, SNRfloor=SNRfloor)

  ptsC = TransformUtils.wrapRad.(rand(bss,N))

  return manikde!(reshape(ptsC,1,:), (:Circular,))
end

#


text:
module AbstractStabilizers

using ..AbstractOps

export AbstractStabilizer

abstract type AbstractStabilizer <: AbstractOp end

end # module AbstractStabilizers


text:
@inline function rand_point()
    return 2.0 .* rand(3) .- 1.0
end

function sample_in_sphere(R::Number, r::Number, n::Integer, max_retry_times::Integer = 100)
    r > R && throw(ArgumentError("particle radius must be smaller than region radius!"))

    points = Vector{Float64}[]
    r′ = r / R
    for _ = 1:n
        retry_times = 0
        new_point = rand_point()
        while norm(new_point) + r′ > 1.0
            new_point = rand_point()
        end
        min_dist = minimum(norm(point - new_point) for point in points; init = Inf)
        while min_dist < 2r′ && retry_times < max_retry_times
            retry_times += 1
            new_point = rand_point()
            while norm(new_point) + r′ > 1.0
                new_point = rand_point()
            end
            min_dist = minimum(norm(point - new_point) for point in points; init = Inf)
        end
        if min_dist >= 2r′
            push!(points, new_point)
        end
    end

    length(points) < n &&
        @warn "Failed to generate enough points. $(length(points))/$n generated."

    return [point ./ r′ for point in points]
end

export sample_in_sphere


text:
__precompile__()

module ExprRules

import TreeView: walk_tree
using StatsBase
using AbstractTrees
using DataStructures  #NodeRecycler

include("interpreter.jl")
using .Interpreter

include("rulenode.jl")
include("grammar_base.jl")
include("expreval.jl")
include("nodelocation.jl")
include("sampling.jl")
include("iterators.jl")
#include("utils.jl")
include("rulenode_operators.jl")

include("cfg.jl")

include("csgrammar/context.jl")
include("csgrammar/constraints.jl")
include("csgrammar/propagators.jl")
include("csgrammar/csg.jl")
include("csgrammar/iterators.jl")
include("csgrammar/sampling.jl")

include("utils.jl")

export
        Grammar,
        ContextFreeGrammar,
        RuleNode,
        NodeLoc,

        ExpressionIterator,
        ContextFreeIterator,

        @cfgrammar,
        max_arity,
        depth,
        node_depth,
        isterminal,
        iseval,
        return_type,
        contains_returntype,
        nchildren,
        child_types,
        get_childtypes,
        nonterminals,
        get_executable,
        sample,
        root_node_loc,
        count_expressions,
        iscomplete,

        SymbolTable,
        interpret,

        change_expr,
	swap_node,
	get_rulesequence,
	rulesoftype,
	rulesonleft,

        NodeRecycler,
        recycle!,

        GrammarContext,
	addparent!,
        copy_and_insert,

        Constraint,
        ValidatorConstraint,
        PropagatorConstraint,
        ComesAfter,
        Ordered,
        Forbidden,

        propagate,
        propagate_index,

        ContextSensitiveGrammar,
        @csgrammar,
        addconstraint!,

        ContextSensitiveIterator,
        propagate_constraints,

        mindepth_map,
        mindepth,
        containedin,
        subsequenceof,
        has_children
        #accumulate_if!

end # module


text:
struct BracketingInterval
    x_lo
    x_mid
    x_hi
    f_lo
    f_mid
    f_hi
    iter
    success::Bool
end

function Base.show(io::IO, bracket::BracketingInterval)
    if bracket.success
        print("Successful seach")
    else
        print("Failed seach")
    end
    print(" with $(bracket.iter)-iterations\n")
    println("\tx\tf(x)")
    println( "lo:\t$(bracket.x_lo)\t$(bracket.f_lo)" )
    println( "mid:\t$(bracket.x_mid)\t$(bracket.f_mid)" )
    println( "hi:\t$(bracket.x_hi)\t$(bracket.f_hi)" )
end

function bracket_min(
    fcn,
    x_lo,
    x_hi;
    step = abs(x_hi-x_lo) * 1e-2,
    step_factor = 1.61,
    max_iter = 1000)

    f_x_lo = fcn(x_lo)
    f_x_hi = fcn(x_hi)

    if f_x_lo < f_x_hi
        # swap the orders
        x_lo, x_hi = x_hi, x_lo
        f_x_lo, f_x_hi = f_x_hi, f_x_lo
        step = -step
    end

    iter = 1
    while iter <= max_iter
        x_next = x_hi + step
        f_next = fcn(x_next)
        if f_next > f_x_hi
            # bracketed interval found
            return BracketingInterval(x_lo, x_hi, x_next, f_x_lo, f_x_hi, f_next, iter, true)
        end

        x_lo, x_hi = x_hi, x_next
        f_x_lo, f_x_hi = f_x_hi, f_next

        step *= step_factor
        iter += 1
    end
    println("Failed to find a bracketing interval in the prescribed number of iterations")
    return BracketingInterval(x_lo, x_hi, x_next, f_x_lo, f_x_hi, f_next, iter, false)

end

struct GoldenSectionSearchStatistics
    x_sol
    y_sol
    iter
end

function golden_section_search(fcn, x_lo, x_hi; max_iter = 100, x_tol = 1e-10)

    rho = 0.5 * (sqrt(5)-1)
    if x_lo > x_hi
        # swap
        x_lo, x_hi = x_hi, x_lo
    end
    a = x_lo
    b = x_hi

    c_fcn(lo, hi) = rho * lo + (1 - rho) * hi
    d_fcn(lo, hi) = (1 - rho) * lo + rho * hi

    y_a = fcn(a)
    y_b = fcn(b)

    c = c_fcn(a, b)
    d = d_fcn(a, b)

    y_c = fcn(c)
    y_d = fcn(d)

    x_sol = NaN
    y_sol = NaN
    iter = 0
    for outer iter = 1:max_iter

        if y_c < y_d
            # choose (a, c, d)
            a, d, b = a, c, d
            y_a, y_d, y_b = y_a, y_c, y_d
            c = c_fcn(a, b)
            y_c = fcn(c)
        else
            # choose (c, d, b)
            a, c, b = c, d, b
            y_a, y_c, y_b = y_c, y_d, y_b
            d = d_fcn(a, b)
            y_d = fcn(d)
        end

        if abs(b - a) <= x_tol
            (x_sol, y_sol) = (y_c < y_d) ? (c, y_c) : (d, y_d)
            break
        end

    end
    if iter == max_iter
        (x_sol, y_sol) = (y_c < y_d) ? (c, y_c) : (d, y_d)
    end

    return GoldenSectionSearchStatistics(x_sol, y_sol, iter)

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
"""
  Model(expr, parts::Vector{FormulaTerm})

Representation of a multipart model composed of a vector of FormulaTerm's.
It is inspired by the [https://github.com/matthieugomez/FixedEffectModels.jl](FixedEffectModels.jl) package
"""
struct Model
  expr
  parts::Vector{FormulaTerm}
end

"""
    @model(formula1, formula2, ...)

Parses a sequence of comma-separated FormulaTerm's representing a multipart model.
"""
macro model(args...)
  Expr(:call, :model_helper, (esc(Base.Meta.quot(a)) for a in args)...)
end

function model_helper(args...)
  expr = args
  parts = FormulaTerm[]
  for part in args
    push!(parts, @eval(@formula($(part.args[2]) ~ $(part.args[3]))))
  end
  Model(expr,parts)
end

function Base.show(io::IO, m::Model)
  print(io, "$(length(m.parts))-part model: [")
  notfirst = false
  for p in m.parts
    if notfirst
      print(io, ", ")
    end
    print(io, p)
    notfirst = true
  end
  print(io, "]")
  nothing
end

function mergerhsterms(a::Term, b::Term)
  if a == b
    a, [1], [1]
  else
    (a, b), [1], [2]
  end
end

mergerhsterms(a, b::Term) = mergerhsterms(a, (b,))
mergerhsterms(a::Term, b) = mergerhsterms((a,), b)

function mergerhsterms(a, b)
  terms = union(a, b)

  ina = findall((in)(a), terms)
  inb = findall((in)(b), terms)

  (terms...,), ina, inb
end

"maps inzero and inpos to ModelMatrix columns (important with factor variables)"
function mapins(inzero, inpos, appliedschema)
  mappedinzero = Int[]
  mappedinpos = Int[]
  from = 1
  to = 0
  for t=1:length(appliedschema)
    to = from + width(appliedschema[t]) - 1
    if t in inzero
      [push!(mappedinzero, i) for i in from:to]
    end
    if t in inpos
      [push!(mappedinpos, i) for i in from:to]
    end
    from = to + 1
  end
  mappedinzero, mappedinpos
end

# special case of a single variable on both model parts
mapins(inzero, inpos, singleterm::AbstractTerm) = [1], [1]

"""
  getrhsterms(m, lhs)

Returns the right-hand-side Terms of the part of model `m` with left-hand-side
symbol `lhs`.
"""
function getrhsterms(m::Model, lhs::Symbol)
  ix = findfirst(p->p.lhs==term(lhs),m.parts)
  @assert ix > 0 "The model is missing a formula with $lhs on its left-hand-side."
  f = m.parts[ix]
  f.rhs
end

StatsModels.missing_omit(data::T, formula::TupleTerm) where T<:ColumnTable =
    missing_omit(NamedTuple{tuple(termvars(formula)...)}(data))

function StatsModels.modelcols(trms, data, counts::AbstractMatrix;
    model::Type{M}=DCR, contrasts=Dict{Symbol,Any}()) where M

  cols = columntable(data)
  cols, nonmissing = missing_omit(cols, trms)
  s = schema(trms, cols, contrasts)
  as = apply_schema(trms, s, M)
  modelcols(as, nonmissing, cols, counts)
end

function StatsModels.modelcols(as, nonmissing, cols, counts::AbstractMatrix)

  covars = modelcols(MatrixTerm(as), cols)

  if !all(nonmissing)
    counts = counts[nonmissing,:]
  end

  covars, counts, as
end

"Similar to StatsModels.TableRegressionModel, but also holds a counts matrix and Model"
struct TableCountsRegressionModel{M,D,C} <: RegressionModel
  model::M    # actual fitted model (result of fit())
  data::D     # data table used
  counts::C   # counts matrix used
  f::Model    # @model specification
  schema      # applied schema
  sprojdir::Union{Nothing,Symbol}   # lhs variable in forward regression
end
TableCountsRegressionModel(model, data, counts, f, schema) = TableCountsRegressionModel(model, data, counts, f, schema, nothing)

function Base.show(io::IO, model::TableCountsRegressionModel)
  println(io, typeof(model))
  println(io)
  println(io, model.f)

  sprojdir = model.sprojdir
  if !isnothing(sprojdir)
    # CIR case show show forward model
    println(io, "\nForward model coefficients for predicting $sprojdir:")

    fwdm = model.model.fwdm

    ct = StatsBase.coeftable(fwdm)

    # all coefnames including prediction target
    cfnames = coefnames(model.schema)
    if isa(cfnames, String)
      cfnames = [cfnames]
    end

    # remove projdir
    projdir = ixprojdir(model.schema, sprojdir)
    deleteat!(cfnames, projdir)

    # prepend intercept (always included in forward regression)
    prepend!(cfnames, ["(Intercept)"])

    # resolve srproj / Z coef names
    znames = srprojnames(model.model.bwdm, length(ct.rownms) - length(cfnames))
    append!(cfnames, znames)

    if length(ct.rownms) == length(cfnames)
        ct.rownms = cfnames
    end

    show(io, ct)
  end
end


text:
ENV["GKSwstype"] = "nul"
using Test
using IntervalSets

import Plots

import ClimaCore
import ClimaCorePlots

OUTPUT_DIR = mkpath(get(ENV, "CI_OUTPUT_DIR", tempname()))

@testset "spectral element 2D cubed-sphere" begin
    R = 6.37122e6

    domain = ClimaCore.Domains.SphereDomain(R)
    mesh = ClimaCore.Meshes.EquiangularCubedSphere(domain, 6)
    grid_topology = ClimaCore.Topologies.Topology2D(mesh)
    quad = ClimaCore.Spaces.Quadratures.GLL{5}()
    space = ClimaCore.Spaces.SpectralElementSpace2D(grid_topology, quad)
    coords = ClimaCore.Fields.coordinate_field(space)

    u = map(coords) do coord
        u0 = 20.0
        α0 = 45.0
        ϕ = coord.lat
        λ = coord.long

        uu = u0 * (cosd(α0) * cosd(ϕ) + sind(α0) * cosd(λ) * sind(ϕ))
        uv = -u0 * sind(α0) * sind(λ)
        ClimaCore.Geometry.UVVector(uu, uv)
    end

    field_fig = Plots.plot(u.components.data.:1)
    @test field_fig !== nothing

    fig_png = joinpath(OUTPUT_DIR, "2D_cubed_sphere_field.png")
    Plots.png(field_fig, fig_png)
    @test isfile(fig_png)
end

@testset "spectral element rectangle 2D" begin
    domain = ClimaCore.Domains.RectangleDomain(
        ClimaCore.Geometry.XPoint(0) .. ClimaCore.Geometry.XPoint(2π),
        ClimaCore.Geometry.YPoint(0) .. ClimaCore.Geometry.YPoint(2π),
        x1periodic = true,
        x2periodic = true,
    )

    n1, n2 = 2, 2
    Nq = 4
    mesh = ClimaCore.Meshes.RectilinearMesh(domain, n1, n2)
    grid_topology = ClimaCore.Topologies.Topology2D(mesh)
    #quad = ClimaCore.Spaces.Quadratures.GLL{Nq}()
    quad = ClimaCore.Spaces.Quadratures.ClosedUniform{Nq + 1}()
    space = ClimaCore.Spaces.SpectralElementSpace2D(grid_topology, quad)
    coords = ClimaCore.Fields.coordinate_field(space)

    space_fig = Plots.plot(space)
    @test space_fig !== nothing

    sinxy = map(coords) do coord
        cos(coord.x + coord.y)
    end

    field_fig = Plots.plot(sinxy)
    @test field_fig !== nothing

    space_png = joinpath(OUTPUT_DIR, "2D_rectangle_space.png")
    field_png = joinpath(OUTPUT_DIR, "2D_rectangle_field.png")
    Plots.png(space_fig, space_png)
    Plots.png(field_fig, field_png)
    @test isfile(space_png)
    @test isfile(field_png)
end

@testset "hybrid finite difference / spectral element 2D" begin
    FT = Float64
    helem = 10
    velem = 40
    npoly = 4

    vertdomain = ClimaCore.Domains.IntervalDomain(
        ClimaCore.Geometry.ZPoint{FT}(0),
        ClimaCore.Geometry.ZPoint{FT}(1000);
        boundary_tags = (:bottom, :top),
    )
    vertmesh = ClimaCore.Meshes.IntervalMesh(vertdomain, nelems = velem)
    vert_center_space = ClimaCore.Spaces.CenterFiniteDifferenceSpace(vertmesh)

    horzdomain = ClimaCore.Domains.IntervalDomain(
        ClimaCore.Geometry.XPoint{FT}(-500) ..
        ClimaCore.Geometry.XPoint{FT}(500),
        periodic = true,
    )
    horzmesh = ClimaCore.Meshes.IntervalMesh(horzdomain; nelems = helem)
    horztopology = ClimaCore.Topologies.IntervalTopology(horzmesh)

    quad = ClimaCore.Spaces.Quadratures.GLL{npoly + 1}()
    horzspace = ClimaCore.Spaces.SpectralElementSpace1D(horztopology, quad)

    hv_center_space = ClimaCore.Spaces.ExtrudedFiniteDifferenceSpace(
        horzspace,
        vert_center_space,
    )
    hv_face_space =
        ClimaCore.Spaces.FaceExtrudedFiniteDifferenceSpace(hv_center_space)

    coords = ClimaCore.Fields.coordinate_field(hv_center_space)

    xcoords_fig = Plots.plot(coords.x)
    @test xcoords_fig !== nothing

    zcoords_fig = Plots.plot(coords.z)
    @test zcoords_fig !== nothing

    xcoords_png = joinpath(OUTPUT_DIR, "hybrid_xcoords_center_field.png")
    zcoords_png = joinpath(OUTPUT_DIR, "hybrid_zcoords_center_field.png")
    Plots.png(xcoords_fig, xcoords_png)
    Plots.png(zcoords_fig, zcoords_png)
    @test isfile(xcoords_png)
    @test isfile(zcoords_png)
end