text:
# Functions in this file have derived from some in the PROJ4.jl pacakge https://github.com/JuliaGeo/Proj4.jl
# so part of the credit goes also for the authors of that package but the main function, "geod" was highly
# modified and does a lot more than the proj4 'geod' model.

abstract type _geodesic end
mutable struct null_geodesic <: _geodesic end

struct PJ_INFO
	major::Cint
	minor::Cint
	patch::Cint
	release::Cstring
	version::Cstring
	searchpath::Cstring
	paths::Ptr{Cstring}
	path_count::Csize_t
end

struct PJ_ELLPS
	id::Cstring
	major::Cstring
	ell::Cstring
	name::Cstring
end

mutable struct geod_geodesic <: _geodesic
	a::Cdouble
	f::Cdouble
	f1::Cdouble
	e2::Cdouble
	ep2::Cdouble
	n::Cdouble
	b::Cdouble
	c2::Cdouble
	etol2::Cdouble
	A3x::NTuple{6, Cdouble}
	C3x::NTuple{15,Cdouble}
	C4x::NTuple{21,Cdouble}
	geod_geodesic() = new()
end

# -------------------------------------------------------------------------------------------------
"""
    geod(lonlat, azim, distance; proj::String="", s_srs::String="", epsg::Integer=0, dataset=false, unit=:m)

Solve the direct geodesic problem.

Args:

- `lonlat`:   - longitude, latitude (degrees). This can be a vector or a matrix with one row only.
- `azimuth`:  - azimuth (degrees) ∈ [-540, 540)
- `distance`: - distance to move from (lat,lon); can be negative, Default is meters but see `unit`
- `proj` or `s_srs`:  - the given projection whose ellipsoid we move along. Can be a proj4 string or an WKT
- `epsg`:     - Alternative way of specifying the projection [Default is WGS84]
- `dataset`:  - If true returns a GMTdataset instead of matrix
- `unit`:     - If `distance` is not in meters use one of `unit=:km`, or `unit=:Nautical` or `unit=:Miles`

The `distance` argument can be a scalar, a Vector, a Vector{Vector} or an AbstractRange. The `azimuth` can be
a scalar or a Vector.

When `azimuth` is a Vector we always return a GMTdataset with the multiple lines. Use this together with a
non-scalar `distance` to get lines with multiple points along the line. The number of points along line does
not even need to be the same. For data, give the `distance` as a Vector{Vector} where each element of `distance`
is a vector with the distances of the points along a line. In this case the number of `distance` elements
must be equal to the number of `azimuth`.

### Returns
- dest - destination after moving for [distance] metres in [azimuth] direction.
- azi  - forward azimuth (degrees) at destination [dest].

## Example: Compute two lines starting at (0,0) with lengths 111100 & 50000, heading at 15 and 45 degrees.

    geod([0., 0], [15., 45], [[0, 10000, 50000, 111100.], [0., 50000]])[1]
"""
geod(lonlat::Vector{<:Real}, azim, distance; proj::String="", s_srs::String="", epsg::Integer=0, dataset=false, unit=:m) =
	geod([lonlat[1] lonlat[2]], azim, distance; proj=proj, s_srs=s_srs, epsg=epsg, dataset=dataset, unit=unit)
function geod(lonlat::Matrix{<:Real}, azim, distance; proj::String="", s_srs::String="", epsg::Integer=0, dataset=false, unit=:m)
	f = unit_factor(unit)
	_dist = isa(distance, AbstractRange) ? collect(Float64, distance) .* f : distance .* f

	proj_string, projPJ_ptr, isgeog = helper_geod(proj, s_srs, epsg)
	dest, azi = helper_gdirect(projPJ_ptr, lonlat, azim, _dist, proj_string, isgeog, dataset)
	proj_destroy(projPJ_ptr)
	return dest, azi
end

function unit_factor(unit=:m)
	# Returns the factor that converts a length to meters, or 1 if unit is unknown
	f = 1.0
	if (unit != :m)
		_u = lowercase(string(unit))		# Parse the units arg
		f = (_u[1] == 'k') ? 1000. : ((_u[1] == 'n') ? 1852.0 : (startswith(_u, "mi") ? 1609.344 : 1.0))
	end
	(unit != :m && f == 1.0) && @warn("Unknown unit ($_u). Ignoring it")
	f
end

# -------------------------------------------------------------------------------------------------
"""
    d, az1, az2 = invgeod(lonlat1::Vector{<:Real}, lonlat2::Vector{<:Real}; proj::String="", s_srs::String="", epsg::Integer=0)

Solve the inverse geodesic problem.

Args:

- `lonlat1`:  - coordinates of point 1 in the given projection
- `lonlat2`:  - coordinates of point 2 in the given projection
- `proj` or `s_srs`:  - the given projection whose ellipsoid we move along. Can be a proj4 string or an WKT
- `epsg`:     - Alternative way of specifying the projection [Default is WGS84]

### Returns
dist - distance between point 1 and point 2 (meters).
azi1 - azimuth at point 1 (degrees) ∈ [-180, 180)
azi2 - (forward) azimuth at point 2 (degrees) ∈ [-180, 180)

Remarks:

If either point is at a pole, the azimuth is defined by keeping the longitude fixed,
writing lat = 90 +/- eps, and taking the limit as eps -> 0+.
"""
invgeod(lonlat1::Vector{<:Real}, lonlat2::Vector{<:Real}; proj::String="", s_srs::String="", epsg::Integer=0) =
	invgeod([lonlat1[1] lonlat1[2]], [lonlat2[1] lonlat2[2]]; proj=proj, s_srs=s_srs, epsg=epsg)
function invgeod(lonlat1::Matrix{<:Real}, lonlat2::Matrix{<:Real}; proj::String="", s_srs::String="", epsg::Integer=0)
	@assert (size(lonlat1) == size(lonlat2) || size(lonlat2,1) == 1) "Both matrices must have same size or second have one row"
	proj_string, projPJ_ptr, isgeog = helper_geod(proj, s_srs, epsg)
	(!isgeog) && (lonlat1 = xy2lonlat(lonlat1, proj_string);	lonlat2 = xy2lonlat(lonlat2, proj_string))	# Convert to geogd first
	d   = Vector{Float64}(undef, size(lonlat1,1))
	az1 = Vector{Float64}(undef, size(lonlat1,1))
	az2 = Vector{Float64}(undef, size(lonlat1,1))
	dist, azi1, azi2 = Ref{Cdouble}(), Ref{Cdouble}(), Ref{Cdouble}()
	for k = 1:size(lonlat1,1)
		kk = (size(lonlat2,1) == 1) ? 1 : k			# To allow the "all against one" case
		ccall((:geod_inverse, libproj), Cvoid, (Ptr{Cvoid},Cdouble,Cdouble,Cdouble, Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
			  pointer_from_objref(get_ellipsoid(projPJ_ptr)), lonlat1[k,2], lonlat1[k,1], lonlat2[kk,2], lonlat2[kk,1], dist, azi1, azi2)
		d[k], az1[k], az2[k] = dist[], azi1[], azi2[]
	end
	proj_destroy(projPJ_ptr)
	return d, az1, az2
end

# -------------------------------------------------------------------------------------------------
"""
    circgeo(lon, lat; radius=X, proj="", s_srs="", epsg=0, dataset=false, unit=:m, np=120, shape="")
or

    circgeo(lonlat; radius=X, proj="", s_srs="", epsg=0, dataset=false, unit=:m, np=120, shape="")

Args:

- `lonlat`:   - longitude, latitude (degrees). If a Mx2 matrix, returns as many segments as number of rows.
                Use this to compute multiple shapes at different positions. In this case output type is
				always a vector of GMTdatasets.
- `radius`:   - The circle radius in meters (but see `unit`) or circumscribing circle for the other shapes
- `proj` or `s_srs`:  - the given projection whose ellipsoid we move along. Can be a proj4 string or an WKT
- `epsg`:     - Alternative way of specifying the projection [Default is WGS84]
- `dataset`:  - If true returns a GMTdataset instead of matrix (with single shapes)
- `unit`:     - If `radius` is not in meters use one of `unit=:km`, or `unit=:Nautical` or `unit=:Miles`
- `np`:       - Number of points into which the circle is descretized (Default = 120)
- `shape`:    - Optional string/symbol with "triangle", "square", "pentagon" or "hexagon" (or just the first char)
                to compute one of those geometries instead of a circle. `np` is ignored in these cases.

### Returns
- circ - A Mx2 matrix or GMTdataset with the circle coordinates

## Example: Compute circle about the (0,0) point with a radius of 50 km

    c = circgeo([0.,0], radius=50, unit=:k)
"""
circgeo(lon::Real, lat::Real; radius::Real=0., proj::String="", s_srs::String="", epsg::Integer=0, dataset::Bool=false, unit=:m, np::Int=120, shape="") =
	circgeo([lon lat]; radius=radius, proj=proj, s_srs=s_srs, epsg=epsg, dataset=dataset, unit=unit, np=np, shape=shape)
circgeo(lonlat::Vector{<:Real}; radius::Real=0., proj::String="", s_srs::String="", epsg::Integer=0, dataset::Bool=false, unit=:m, np::Int=120, shape="") =
	circgeo([lonlat[1] lonlat[2]]; radius=radius, proj=proj, s_srs=s_srs, epsg=epsg, dataset=dataset, unit=unit, np=np, shape=shape)

function circgeo(lonlat::Matrix{<:Real}; radius::Real=0., proj::String="", s_srs::String="", epsg::Integer=0, dataset::Bool=false, unit=:m, np::Int=120, shape="")
	(radius == 0) && error("Must provide circle Radius. Obvious, not?")
	_shape = lowercase(string(shape))
	if     (_shape == "")     azim = collect(Float64, linspace(0, 360, np))	# The default (circle)
	elseif (_shape[1] == 't') azim = [0.,  120, 240, 0]				# Triangle
	elseif (_shape[1] == 's') azim = [45., 135, 225, 315, 45]		# Square
	elseif (_shape[1] == 'p') azim = [0.,  72, 144, 216, 288, 360]	# Pentagon
	elseif (_shape[1] == 'h') azim = [0.,  60, 120, 180, 240, 300, 360]	# Hexagon
	else   error("Bad shape name ($(shape))")
	end
	n_shapes = size(lonlat, 1)
	if (n_shapes > 1)	# Multiple shapes require a GMTdataset output
		D = Array{GMTdataset}(undef, n_shapes)
		for k = 1:n_shapes
			D[k] = geod(lonlat[k,:], azim, radius; proj=proj, s_srs=s_srs, epsg=epsg, dataset=true, unit=unit)[1]
		end
	else		# Here the output may be a GMTdadaset or a simple Mx2 matrix depending on the 'dataset' value.
		D = geod(lonlat[1,:], azim, radius; proj=proj, s_srs=s_srs, epsg=epsg, dataset=dataset, unit=unit)[1]
	end
	D
end

# -------------------------------------------------------------------------------------------------
"""
    buffergeo(D::GMTdataset; width=0, unit=:m, np=120, flatstart=false, flatend=false, epsg::Integer=0, tol=0.01)
or

    buffergeo(line::Matrix; width=0, unit=:m, np=120, flatstart=false, flatend=false, proj::String="", epsg::Integer=0, tol=0.01)
or

    buffergeo(fname::String; width=0, unit=:m, np=120, flatstart=false, flatend=false, proj::String="", epsg::Integer=0, tol=0.01)

Computes a buffer arround a poly-line. This calculation is performed on a ellipsoidal Earth (or other planet)
using the GeographicLib (via PROJ4) so it should be very accurate.

### Parameters
- `D` | `line` | fname: - the geometry. This can either be a GMTdataset (or vector of it), a Mx2 matrix, the name
                          of file that can be read as a GMTdataset by `gmtread()` or a GDAL AbstractDataset object
- `width`:  - the buffer width to be applied. Expressed meters (the default), km or Miles (see `unit`)
- `unit`:   - If `width` is not in meters use one of `unit=:km`, or `unit=:Nautical` or `unit=:Miles`
- `np`:     - Number of points into which circles are descretized (Default = 120)
- `flatstart` - When computing buffers arround poly-lines make the start *flat* (no half-circle)
- `flatend`   - Same as `flatstart` but for the buffer end
- `proj`  - If line data is in Cartesians but with a known projection pass in a PROJ4 string to allow computing the buffer
- `epsg`  - Same as `proj` but using an EPSG code
- `tol`   - At the end simplify the buffer line with a Douglas-Peucker procedure. Use TOL=0 to NOT do the line
            simplification, or use any other value in degrees. Default computes it as 0.5% of buffer width.

### Returns
A GMT dataset or a vector of it (when input is Vector{GMTdataset})

## Example: Compute a buffer with 50000 m width

    D = buffergeo([0 0; 10 10; 15 20], width=50000);
"""
buffergeo(fname::String; width=0, unit=:m, np=120, flatstart=false, flatend=false, epsg::Integer=0, tol=-1.0) =
	buffergeo(gmtread(fname); width=width, unit=unit, np=np, flatstart=flatstart, flatend=flatend, epsg=epsg, tol=tol)

buffergeo(ds::Gdal.AbstractDataset; width=0, unit=:m, np=120, flatstart=false, flatend=false, epsg::Integer=0, tol=-1.0) =
	buffergeo(gmt2gd(ds); width=width, unit=unit, np=np, flatstart=flatstart, flatend=flatend, epsg=epsg, tol=tol)

buffergeo(D::GMTdataset; width=0, unit=:m, np=120, flatstart=false, flatend=false, epsg::Integer=0, tol=-1.0) =
	buffergeo(D.data; width=width, unit=unit, np=np, flatstart=flatstart, flatend=flatend, proj=D.proj4, epsg=epsg, tol=tol)[1]

function buffergeo(D::Vector{<:GMTdataset}; width=0, unit=:m, np=120, flatstart=false, flatend=false, epsg::Integer=0, tol=-1.0)
	_D = Vector{GMTdataset}(undef, length(D))
	for k = 1:length(D)
		_D[k], = buffergeo(D[k].data; width=width, unit=unit, np=np, flatstart=flatstart, flatend=flatend,
		                   proj=D[1].proj4, epsg=epsg, tol=tol)
	end
	return (length(_D) == 1) ? _D[1] : _D		# Drop the damn Vector singletons
end

function buffergeo(line::Matrix{<:Real}; width=0, unit=:m, np=120, flatstart=false, flatend=false, proj::String="", epsg::Integer=0, tol=-1.0)
	# This function can be a bit optimized with a clever choice of which points to compute on the circle.
	# Points close to the segment axis are ofc nover needed.
	(width == 0) && error("Must provide the buffer width. Obvious, not?")
	dist, azim, = invgeod(line[1:end-1, :], line[2:end, :], proj=proj, epsg=epsg)	# distances and azimuths between the polyline vertices
	#(line[1, 1:2] == line[end, 1:2]) && (flatstart = flatend = false)		# Polygons can't have start/end flat edges
	width *= unit_factor(unit)
	n_seg = length(azim)
	for n = 1:n_seg
		seg = [geod(line[n,:], azim[n], 0:width/4:dist[n], proj=proj, epsg=epsg)[1]; line[n+1:n+1,:]]
		n_vert = size(seg,1)
		for k = 2:n_vert
			if (k == 2)
				c0 = circgeo(seg[1,:], radius=width, np=np, proj=proj, epsg=epsg);	trim_dateline!(c0, seg[1,1])
				c  = circgeo(seg[2,:], radius=width, np=np, proj=proj, epsg=epsg);	trim_dateline!(c , seg[2,1])
				global _D = polyunion(c0, c)
				continue
			end
			c = circgeo(seg[k,:], radius=width, np=np, proj=proj, epsg=epsg)
			trim_dateline!(c, seg[k,1])
			_D = polyunion(_D, c)
		end
		if (n == 1)  global D = _D
		else         D = polyunion(_D, D)
		end
	end

	# At this point we still have a "wavy" buffer resulting from the union of the circles. Ideally we should
	# be able to use GMT's mapproject to find only the points that are at the exact 'width' distance from the
	# line, but at this moment there seems to be an issue in GMT and distances are shorter. However with can
	# do some cheap statistics to get read of the more inner points a get an almost perfec buffer.
	Ddist2line = mapproject(D, L=(line=line, unit=:e))		# Find the distance from buffer points to input line
	d_mean = mean(view(Ddist2line[1].data, :, 3))
	ind = view(Ddist2line[1].data, :, 3) .>= d_mean
	ind[1], ind[end] = true, true			# Ensure first and last are not removed
	D[1].data = D[1].data[ind, :]			# Remove all points that are less 1% above the mean

	(proj == "" && epsg == 0) && (D[1].proj4 == "+proj=longlat +datum=WGS84 +no_defs")
	(proj != "") && (D[1].proj4 = proj)
	(epsg != 0) && (D[1].proj4 = toPROJ4(importEPSG(epsg)))
	(tol < 0) && (tol = width * 0.005 / 111111)		# Tolerance as 0.5% of buffer with converted to degrees
	return (tol > 0) ? simplify(D, tol) : D			# If TOL > 0 do a DP simplify on final buffer
end

function trim_dateline!(mat, lon0)
	# When the the mat polygon (a circle normally) crosses the dateline, trim it from the standpoint of its centroid
	mi, ma = extrema(view(mat,:,1))
	if ((ma - mi) > 180)		# Polygon crosses the dateline
		if (lon0 > 0)  mat[view(mat,:,1) .< 0, 1] .= 180.  # and centroid is at the 11:xxx side
		else           mat[view(mat,:,1) .> 0, 1] .= -180.
		end
	end
end

# -------------------------------------------------------------------------------------------------
"""
    function orthodrome(D; step=0, unit=:m, np=0, proj::String="", epsg::Integer=0)

Generate a orthodrome line(s) (shortest distace) on an ellipsoid. Input data can be two or more points.
In later case each line segment is descretized at `step` increments,

### Parameters

- `D`: - the input points. This can either be a GMTdataset (or vector of it), a Mx2 matrix, the name
         of file that can be read as a GMTdataset by `gmtread()` or a GDAL AbstractDataset object
- `step`:  - Incremental distance at which the segment line is descretized in meters(the default), but see `unit`
- `unit`:  - If `step` is not in meters use one of `unit=:km`, or `unit=:Nautical` or `unit=:Miles`
- `np`:    - Number of intermediate points between poly-line vertices (alternative to `step`)
- `proj`  - If line data is in Cartesians but with a known projection pass in a PROJ4 string
- `epsg`  - Same as `proj` but using an EPSG code

### Returns
A Mx2 matrix with the lon lat of the points along the orthodrome when input is a matrix. A GMT dataset
or a vector of it (when input is Vector{GMTdataset}).

## Example: Compute an orthodrome between points (0,0) and (30,50) discretized at 100 km steps.

    mat = orthodrome([0 0; 30 50], step=100, unit=:k);
"""
orthodrome(fname::String; step=0, unit=:m, np=0, proj::String="", epsg::Integer=0) =
	orthodrome(gmtread(fname); step=step, unit=unit, np=np, proj=proj, epsg=epsg)

orthodrome(lon1::Real, lat1::Real, lon2::Real, lat2::Real; step=0, unit=:m, np=0, proj::String="", epsg::Integer=0) =
	orthodrome([lon1 lat1; lon2 lat2]; step=step, unit=unit, np=np, proj=proj, epsg=epsg)

orthodrome(ds::Gdal.AbstractDataset; step=0, unit=:m, np=0, proj::String="", epsg::Integer=0) =
	orthodrome(gmt2gd(ds); step=step, unit=unit, np=np, proj=proj, epsg=epsg)

function orthodrome(D::Vector{<:GMTdataset}; step=0, unit=:m, np=0, proj::String="", epsg::Integer=0)
	_D = Vector{GMTdataset}(undef, length(D))
	for k = 1:length(D)
		_D[k] = mat2ds(orthodrome(D[k].data; step=step, unit=unit, np=np, proj = (proj == "") ? D[k].proj4 : proj, epsg=epsg))[1]
	end
	return (length(_D) == 1) ? _D[1] : _D		# Drop the damn Vector singletons
end
function orthodrome(D::GMTdataset; step=0, unit=:m, np=0, proj::String="", epsg::Integer=0)
	mat2ds(orthodrome(D.data; step=step, unit=unit, np=np, proj = (proj == "") ? D.proj4 : proj, epsg=epsg))[1]
end

function orthodrome(line::Matrix{<:Real}; step=0, unit=:m, np=0, proj::String="", epsg::Integer=0)
	(step == 0 && np == 0) && error("Must provide either a 'step' or a 'np' (number of points).")
	(size(line, 1) == 1) && error("Lines cannot have a single point.")
	(np > 0 && size(line, 1) != 2) && error("Number of intermediate points cannot be used with polylines.")
	step *= unit_factor(unit)
	dist, azim, = invgeod(line[1:end-1, :], line[2:end, :], proj=proj, epsg=epsg)	# dist and azims between the polyline vertices
	(np > 0) && (step = dist[1] / (np + 1))
	seg = geod(line[1,:], azim[1], 0:step:dist[1], proj=proj, epsg=epsg)[1]
	if (size(line, 1) == 2)
		(hypot((seg[end,:] .- line[2,:])...) > 1) && (seg = [seg; line[2:2, :]])	# If increment does not land exactly in end point
	else
		for n = 2:size(line,1)-1
			_seg = geod(line[n,:], azim[n], 0:step:dist[n], proj=proj, epsg=epsg)[1]
			(hypot((_seg[end,:] .- line[n,:])...) > 1) && (_seg = [_seg; line[n+1:n+1, :]])
			seg = [seg; _seg]			# Horrible but we can't know in advance the number of final points
		end
	end
	seg
end

# -------------------------------------------------------------------------------------------------
"""
    function loxodrome(lon1,lat1,lon2,lat2; step=0, unit=:m, np=0, proj::String="", epsg::Integer=0)
or

    function loxodrome(D; step=0, unit=:m, np=0, proj::String="", epsg::Integer=0)

Generate a loxodrome (rhumb line) on an ellipsoid. Input data can be two or more points.
In later case each line segment is descretized at `step` increments,

### Parameters

- `D`: - the input points. This can either be a 2x2 matrix or a GMTdataset. Note that only the first 2 points are used.
- `step`:  - Incremental distance at which the segment line is descretized in meters(the default), but see `unit`
- `unit`:  - If `step` is not in meters use one of `unit=:km`, or `unit=:Nautical` or `unit=:Miles`
- `np`:    - Number of intermediate points to be generated between end points (alternative to `step`)
- `proj`  - If line data is in Cartesians but with a known projection pass in a PROJ4 string
- `epsg`  - Same as `proj` but using an EPSG code

### Returns
A Mx2 matrix with the on lat of the points along the loxodrome when input is a matrix or the 2 pairs of points.
A GMTdataset when the input is GMTdataset.

## Example: Compute an loxodrome between points (0,0) and (30,50) discretized at 100 km steps.

    loxo = loxodrome([0 0; 30 50], step=100, unit=:k);
"""
loxodrome(lon1::Real, lat1::Real, lon2::Real, lat2::Real; step=0, unit=:m, np=0, proj::String="", epsg::Integer=0) =
	loxodrome([lon1 lat1; lon2 lat2]; step=step, unit=unit, np=np, proj=proj, epsg=epsg)

function loxodrome(D::GMTdataset; step=0, unit=:m, np=0, proj::String="", epsg::Integer=0)
	mat2ds(loxodrome(D.data; step=step, unit=unit, np=np, proj = (proj == "") ? D.proj4 : proj, epsg=epsg))[1]
end

function loxodrome(line::Matrix{<:Real}; step=0, np=0, unit=:m, proj::String="", epsg::Integer=0)
	(step == 0 && np == 0) && error("Must provide either a 'step' or a 'np' (number of points).")
	(size(line, 1) == 1) && error("Lines cannot have a single point.")
	step *= unit_factor(unit)
	dist, azim = loxodrome_inverse(line[1,1], line[1,2], line[2,1], line[2,2])
	(np > 0) && (step = dist / (np + 1))
	distances = collect(1:step:dist)
	((dist - distances[end]) > 1) && append!(distances, dist)	# Make sure we don't loose destiny point
	loxo = Array{Float64,2}(undef, length(distances), 2)
	for k = 1:length(distances)
		loxo[k, :] = loxodrome_direct(line[1,1], line[1,2], azim, distances[k])
	end
	loxo
end

# -------------------------------------------------------------------------------------------------
function helper_gdirect(projPJ_ptr, lonlat, azim, dist, proj_string, isgeog, dataset)
	lonlat = Float64.(lonlat)
	(!isgeog) && (lonlat = xy2lonlat(lonlat, proj_string))		# Need to convert to geogd first
	ggd = get_ellipsoid(projPJ_ptr)

	(isvector(azim) && isa(dist, Real)) && (dist = [dist,])		# multi-points. Just make 'dist' a vector to reuse a case below

	if (isa(azim, Real) && isa(dist, Real))						# One line only with one end-point
		dest, azi = _geod_direct!(ggd, copy(lonlat), azim, dist)
		(!isgeog) && (dest = lonlat2xy(dest, proj_string))
		(dataset) && (dest = helper_gdirect_SRS(dest, proj_string, wkbPoint))
	elseif (isa(azim, Real) && isvector(dist))					# One line only with several points along it
		dest, azi = Array{Float64}(undef, length(dist), 2), Vector{Float64}(undef, length(dist))
		for k = 1:length(dist)
			d, azi[k] = _geod_direct!(ggd, copy(lonlat), azim, dist[k])
			dest[k, :] = (isgeog) ? d : lonlat2xy(d, proj_string)
		end
		(dataset) && (dest = helper_gdirect_SRS([dest azi], proj_string, wkbLineString))	# Return a GMTdataset
	elseif (isvector(azim) && length(dist) == 1)				# A circle (dist has became a vector in 4rth line)
		dest = Array{Float64}(undef, length(azim), 2)
		for np = 1:length(azim)
			d = _geod_direct!(ggd, copy(lonlat), azim[np], dist[1])[1]
			dest[np, 1:2] = (isgeog) ? d : lonlat2xy(d, proj_string)
		end
		(dataset) &&		# Return a GMTdataset
			(dest = GMTdataset(dest, Vector{String}(), "", Vector{String}(), startswith(proj_string, "+proj") ? proj_string : "", "", wkbPolygon))
		azi = nothing
	elseif (isvector(azim) && isvector(dist))					# Multi-lines with variable length and/or number of points
		n_lines = length(azim)							# Number of lines
		(!isa(dist, Vector)) && (dist = vec(dist))
		(!isa(dist, Vector) && !isa(dist, Vector{<:Vector})) && error("The 'distances' input MUST be a Vector or a Vector{Vector}")
		if (!isa(dist, Vector{<:Vector}))				# If not, make it a Vector{Vector} to use the same algo below
			isa(dist, Matrix) && (dist = vec(dist))		# Because we accepted also 1-row or 1-col matrices
			Vdist = Vector{Vector{Float64}}(undef, n_lines)
			[Vdist[k] = dist for k = 1:n_lines]
		else
			(length(dist) != n_lines) && (proj_destroy(projPJ_ptr);error("Number of distance vectors MUST be equal to number of azimuths"))
			Vdist = dist
		end
		D = Vector{GMTdataset}(undef, n_lines)
		for nl = 1:n_lines
			n_pts = length(Vdist[nl])					# Number of points in this line
			dest = Array{Float64}(undef, n_pts, 3)		# Azimuth goes into the D too
			for np = 1:n_pts
				d, azi = _geod_direct!(ggd, copy(lonlat), azim[nl], Vdist[nl][np])
				dest[np, 1:2] = (isgeog) ? d : lonlat2xy(d, proj_string)
				dest[np, 3] = azi		# Fck language that makes it a pain to try anything vectorized
			end
			D[nl] = GMTdataset(dest, Vector{String}(), "", Vector{String}(), "", "", wkbLineString)
			helper_gdirect_SRS(dest, proj_string, wkbLineString, D[nl])	# Just assign the SRS
		end
		return D, nothing		# Here both the point coordinates and the azim are in the GMTdataset
	else
		error("'azimuth' MUST be either a scalar or a 1-dim array, and 'distance' may also be a Vector{Vector}")
	end
	return dest, azi
end

# -------------------------------------------------------------------------------------------------
function helper_gdirect_SRS(mat, proj_string::String, geom, D=GMTdataset())
	# Convert the output of geod_direct into a GMTdataset and, if possible, assign it a SRS
	# If a 'D' is sent in, we only (eventually) assign it an SRS
	isempty(D) && (D = GMTdataset([mat[1] mat[2]], Vector{String}(), "", Vector{String}(), "", "", geom))
	(startswith(proj_string, "+proj")) && (D.proj4 = proj_string)
	D
end

# -------------------------------------------------------------------------------------------------
function _geod_direct!(geod::geod_geodesic, lonlat, azim, distance)
	p = pointer(lonlat)
	azi = Ref{Cdouble}()		# the (forward) azimuth at the destination
	ccall((:geod_direct, libproj), Cvoid, (Ptr{Cvoid},Cdouble,Cdouble,Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble}, Ptr{Cdouble}),
		  pointer_from_objref(geod), lonlat[2], lonlat[1], azim, distance, p+sizeof(Cdouble), p, azi)
	lonlat, azi[]
end

# -------------------------------------------------------------------------------------------------
function helper_geod(proj::String, s_srs::String, epsg::Integer)::Tuple{String, Ptr{Nothing}, Bool}
	# 'proj' and 's_srs' are synonyms.
	# Return the projection string ans also if the projection is geogs.
	if     (proj  != "")  prj_string = proj
	elseif (s_srs != "")  prj_string = s_srs
	elseif (epsg > 0)     prj_string = toPROJ4(importEPSG(epsg))
	else                  prj_string = "+proj=longlat +datum=WGS84 +no_defs"
	end
	(startswith(prj_string, "GEOGC")) && (prj_string = toPROJ4(importWKT(prj_string)))
	prj_string, proj_create(prj_string), (startswith(prj_string, "+proj=longl") || startswith(prj_string, "+proj=latl") || epsg == 4326)
end

# -------------------------------------------------------------------------------------------------
function geod_geodesic(a::Cdouble, f::Cdouble)::geod_geodesic
	geod = geod_geodesic()
	ccall((:geod_init, libproj), Cvoid, (Ptr{Cvoid}, Cdouble, Cdouble), pointer_from_objref(geod), a, f)
	geod
end

# -------------------------------------------------------------------------------------------------
function get_ellipsoid(projPJ_ptr::Ptr{Cvoid})::geod_geodesic
	#a, ecc2 = pj_get_spheroid_defn(projPJ_ptr)
	#geod_geodesic(a, 1-sqrt(1-ecc2))
	a, inv_f, = proj_ellipsoid_get_parameters(C_NULL, projPJ_ptr)
	geod_geodesic(a, 1/inv_f)
end

# -------------------------------------------------------------------------------------------------
function proj_create(proj_str::String, ctx=C_NULL)
	# Returns an Object that must be unreferenced with proj_destroy()
	# THIS GUY IS VERY EXPENSIVE. TRY TO MINIMIZE ITS USAGE.
	#projPJ_ptr = ccall((:proj_create, libproj), Ptr{Cvoid}, (Ptr{Cvoid}, Cstring), ctx, proj_str)
	projPJ_ptr = ccall((:proj_create, libproj), Ptr{Cvoid}, (Ptr{Cvoid}, Cstring), ctx, proj2wkt(proj_str))
	(projPJ_ptr == C_NULL) && error("Could not parse projection: \"$proj_str\"")
	projPJ_ptr
end

proj_create_crs_to_crs(s_crs, t_crs, area, ctx=C_NULL) =
	ccall((:proj_create_crs_to_crs, libproj), Ptr{Cvoid}, (Ptr{Cvoid}, Cstring, Cstring, Ptr{Cvoid}), ctx, s_crs, t_crs, area)

# -------------------------------------------------------------------------------------------------
function proj_destroy(projPJ_ptr::Ptr{Cvoid})	# Free C datastructure associated with a projection.
	@assert projPJ_ptr != C_NULL
	ccall((:proj_destroy, libproj), Ptr{Cvoid}, (Ptr{Cvoid},), projPJ_ptr)
end

#= -------------------------------------------------------------------------------------------------
function pj_get_spheroid_defn(proj_ptr::Ptr{Cvoid})
	a = Ref{Cdouble}()		# major_axis
	ecc2 = Ref{Cdouble}()	# eccentricity squared
	ccall((:pj_get_spheroid_defn, libproj), Cvoid, (Ptr{Cvoid}, Ptr{Cdouble}, Ptr{Cdouble}), proj_ptr, a, ecc2)
	a[], ecc2[]
end
=#

# -------------------------------------------------------------------------------------------------
function proj_info()
	pji = ccall((:proj_info, libproj), PJ_INFO, ())
	println(unsafe_string(pji.release), "\n", "Located at: ", unsafe_string(pji.searchpath))
end

#=
function _transform!(src_ptr::Ptr{Cvoid}, dest_ptr::Ptr{Cvoid}, point_count::Integer, point_stride::Integer,
                     x::Ptr{Cdouble}, y::Ptr{Cdouble}, z::Ptr{Cdouble})
	@assert src_ptr != C_NULL && dest_ptr != C_NULL
	err = ccall((:pj_transform, libproj), Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Cint, Ptr{Cdouble}, Ptr{Cdouble},
                Ptr{Cdouble}), src_ptr, dest_ptr, point_count, point_stride, x, y, z)
	err != 0 && error("transform error: $(_strerrno(err))")
end

_transform!(s_ptr::Ptr{Cvoid}, t_ptr::Ptr{Cvoid}, pos::Vector{Cdouble}) = _transform!(s_ptr, t_ptr, reshape(pos, 1, length(pos)))
function _transform!(s_ptr::Ptr{Cvoid}, t_ptr::Ptr{Cvoid}, position::Array{Cdouble,2})
	@assert s_ptr != C_NULL && t_ptr != C_NULL
	npoints, ndim = size(position)
	@assert ndim >= 2

	x = pointer(position)
	y = x + sizeof(Cdouble)*npoints
	z = (ndim == 2) ? Ptr{Cdouble}(C_NULL) : x + 2*sizeof(Cdouble)*npoints

	_transform!(s_ptr, t_ptr, npoints, 1, x, y, z)
	position
end

mutable struct Projection
    #ctx::Context   # Projection context object
    rep::Ptr{Cvoid} # Pointer to internal projPJ struct
    geod::_geodesic
end

function Projection(proj_ptr::Ptr{Cvoid})
    proj = Projection(proj_ptr, null_geodesic())
    finalizer(freeProjection, proj)
    proj
end

Projection(proj_string::String) = Projection(proj_create(proj_string))
function freeProjection(proj::Projection)
    proj_destroy(proj.rep)
    proj.rep = C_NULL
end

function is_latlong(proj_ptr::Ptr{Cvoid})
	@assert proj_ptr != C_NULL
	ccall((:pj_is_latlong, libproj), Cint, (Ptr{Cvoid},), proj_ptr) != 0
end
=#

function proj_get_ellipsoid(proj_ptr::Ptr{Cvoid}, ctx=C_NULL)
	# Returned object must be unreferenced with proj_destroy()
	ccall((:proj_get_ellipsoid, libproj), Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}), ctx, proj_ptr)
end

function proj_ellipsoid_get_parameters(ellipsoid::Ptr{Cvoid}=NULL, proj_ptr::Ptr{Cvoid}=NULL, ctx=NULL)
	# proj_ellipsoid_get_parameters(PJ_CONTEXT *ctx, const PJ *ellipsoid, double *out_semi_major_metre, double *out_semi_minor_metre, int *out_is_semi_minor_computed, double *out_inv_flattening)
	semi_major = Ref{Cdouble}()		# semi-major axis in meters
	semi_minor = Ref{Cdouble}()		# semi-minor axis in meters
	inv_flattening = Ref{Cdouble}()		# semi-minor axis in meters
	is_semi_minor_computed = Ref{Cint}()		#
	need_destroy = (ellipsoid == NULL)
	(ellipsoid == NULL) && (ellipsoid = proj_get_ellipsoid(proj_ptr))
	ccall((:proj_ellipsoid_get_parameters, libproj), Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cdouble}), ctx, ellipsoid, semi_major, semi_minor, is_semi_minor_computed, inv_flattening)
	need_destroy && proj_destroy(ellipsoid)		# If it was created above
	semi_major[], inv_flattening[], semi_minor[]
end

"""
    proj2wkt(proj4_str::String, pretty::Bool=false)

Convert a PROJ4 string into the WKT form. Use `pretty=true` to return a more human readable text.
"""
function proj2wkt(proj4_str::String; pretty::Bool=false)
	!startswith(proj4_str, "+proj=") && error("$(proj4_str) is not a valid proj4 string")
	toWKT(importPROJ4(proj4_str), pretty)
end

"""
    wkt2proj(wkt_str::String)

Convert a WKT SRS string into the PROJ4 form.
"""
wkt2proj(wkt_str::String) = toPROJ4(importWKT(wkt_str))

"""
    epsg2proj(code::Integer)

Convert a EPSG code into the PROJ4 form.
"""
epsg2proj(code::Integer)  = toPROJ4(importEPSG(code))

"""
    epsg2wkt(code::Integer, pretty::Bool=false)

Convert a EPSG code into the WKT form. Use `pretty=true` to return a more human readable text.
"""
epsg2wkt(code::Integer; pretty::Bool=false) = toWKT(importEPSG(code), pretty)


text:
module Athesis

using TOML
using Logging
using BasicModelInterface

const BMI = BasicModelInterface

include("athesis/utils.jl")
include("athesis/initialize.jl")
include("athesis/run.jl")
include("athesis/equations.jl")
include("athesis/external_forcing.jl")
include("athesis/postprocessing.jl")
include("athesis/model.jl")
include("athesis/model_input.jl")
include("athesis/grids.jl")
include("athesis/boundary_conditions.jl")
include("athesis/gridloop.jl")
include("athesis/kernels_pressure_equation.jl")
include("athesis/kernels_darcy_equation.jl")
include("athesis/kernels.jl")
include("athesis/operators.jl")
include("athesis/fields.jl")
include("athesis/bmi.jl")

export getDefaultInput,
    initSimulation,
    runSimulation!,
    doExplicitIter!,
    plotSimulation,

    # Data structures
    Grid,
    Model,
    ModelInput,
    Parameters,
    State,
    TimeData,
    SolverData,
    Simulation,

    # Macros
    @synctimeit,
    @withCUDA

end

ext: jl


text:
unorm(U) = sum(map((x->x.^2),U))
# unorm(U::NTuple{1,T}) where {T} = first(U).^2 # specialize for singleton tuples
# unorm(U::AbstractArray) = U.^2 # scalar case
# unorm(U::AbstractFloat) = U^2 # scalar case

"function primitive_to_conservative_nd(rho,u,v,p)

    convert primitive variables (ρ,U,p) to conservative vars (ρ,ρU,E).
    n-dimensional version where U = tuple(u1,...,u_d)"
function primitive_to_conservative_nd(rho,U,p)
    rhoU = (x->rho.*x).(U)
    Unorm = unorm(U)
    E = @. p/(γ-1) + .5*rho*Unorm
    return (rho,rhoU,E)
end

#####
##### functions of conservative variables
#####

"function pfun_nd(rho, rhoU, E)
    pressure as a function of conservative variables (n-dimensional version).
    n-dimensional version where U = tuple(u1,...,u_d)"
function pfun_nd(rho, rhoU, E)
    rhoUnorm2 = unorm(rhoU)./rho
    return @. (γ-1)*(E - .5*rhoUnorm2)
end

"function betafun_nd(rho,rhoU,E)
    inverse temperature (used in entropy conservative fluxes)"
function betafun_nd(rho,rhoU,E)
    p = pfun_nd(rho,rhoU,E)
    return (@. rho/(2*p))
end

# "function rhoe_ufun_nd(rho, rhoU, E)
#     specific energy as a function of conservative variables"
# function rhoe_ufun_nd(rho, rhoU, E)
#     return pfun_nd(rho, rhoU, E) / (γ-1)
# end

"function sfun(rho, rhoU, E)
    Specific entropy as a function of conservative variables"
function sfun_nd(rho, rhoU, E)
    p = pfun_nd(rho, rhoU, E)
    return @. log(p/(rho^γ))
end

"function Sfun(rho,rhoU,E)
    Mathematical entropy as a function of conservative variables"
function Sfun_nd(rho, rhoU, E)
    return -rho.*sfun_nd(rho, rhoU, E)*entropy_scaling
end

"function v_ufun(rho, rhoU, E)
    Entropy variables as functions of conservative vars"
function v_ufun_nd(rho, rhoU, E)
    s = sfun_nd(rho,rhoU,E)
    p = pfun_nd(rho,rhoU,E)

    v1 = (@. (γ + 1 - s) - (γ-1)*E/p)
    vU = (x->@. x*(γ-1)/p).(rhoU)
    vE = (x->@. x*(γ-1)/p)(-rho)

    # v1,vU,vE = scale_entropy_output(v1, vU, vE)
    return v1,vU,vE
end

#####
##### functions of entropy variables
#####
"function s_vfun(v1,vU,vE)
    entropy as function of entropy variables"
function s_vfun_nd(v1,vU,vE)
    vUnorm = unorm(vU)
    return @. γ - v1 + vUnorm/(2*vE)
end

"function rhoe_vfun(v1,vU,vE)
    specific energy as function of entropy variables"
function rhoe_vfun_nd(v1,vU,vE)
    s = s_vfun_nd(v1,vU,vE)
    return (@. ((γ-1)/((-vE)^γ))^(1/(γ-1)) * exp(-s/(γ-1)))
end

"function u_vfun(v1,vU,vE)
    Conservative vars as functions of entropy variables"
function u_vfun_nd(v1,vU,vE)
    # v1,vU,vE = scale_entropy_input(v1,vU,vE)
    rhoeV     = rhoe_vfun_nd(v1,vU,vE)
    vUnorm    = unorm(vU)
    rho       = (@. rhoeV*(-vE))
    rhoU      = (x->rhoeV.*x).(vU)
    E         = (@. rhoeV*(1-vUnorm/(2*vE)))
    return (rho,rhoU,E)
end

"function conservative_to_primitive_beta_nd(rho,rhoU,E)
    converts conservative variables to `primitive' variables which make
    evaluating EC fluxes simpler."
function conservative_to_primitive_beta_nd(rho,rhoU,E)
    return rho, (x->x./rho).(rhoU), betafun_nd(rho,rhoU,E)
end

# function dUdV_explicit(v1,vU1,vU2,vE)
#     rho,rhou,rhov,E = u_vfun(v1,vU1,vU2,vE)
#     u,v = (x->x./rho).((rhou,rhov))
#     p = pfun(rho,rhou,rhov,E)
#     a2 = γ*p/rho
#     H = a2/(γ-1) + (u^2+v^2)/2
#
#     dUdV = @SMatrix [rho  rhou        rhov        E;
#                      rhou rhou*u + p  rhou*v      rhou*H;
#                      rhov rhov*u      rhov*v + p  rhov*H;
#                      E    rhou*H      rhov*H      rho*H^2-a2*p/(γ-1)]
#
#     return dUdV*(1/(γ-1))
# end
#
# function dVdU_explicit(rho,rhou,rhov,E)
#     rhoe = rhoe_ufun(rho,rhou,rhov,E)
#     V = v_ufun(rho,rhou,rhov,E)
#     k = .5*(V[2]^2+V[3]^2)/V[4]
#
#     dVdU = @SMatrix [γ+k^2      k*V[2]          k*V[3]         V[4]*(k+1);
#                     k*V[2]      V[2]^2-V[4]     V[2]*V[3]      V[2]*V[4];
#                     k*V[3]      V[2]*V[3]       V[3]^2-V[4]    V[3]*V[4]
#                     V[4]*(k+1)  V[2]*V[4]       V[3]*V[4]      V[4]^2]
#     return -dVdU/(rhoe*V[4])
# end

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
## utility functions
import OffsetArrays
import ImageFiltering
import ShiftedArrays

export cov_avg!
export boxsmooth!
export outest_bounds #

"""
    outest_bounds(cx,sx) -> px0

Helper function to find maximum padding in pixels required to accomodate all query points `cx` outside of the image size 1:`sx`.

# Arguments:
- `cx`: list of integer star centers (in either x or y)
- `sx`: image dimension along the axis indexed by `cx`

# Outputs:
- `px0`: maximum padding in pixels required to accomodate all query points
"""
function outest_bounds(cx,sx)
    px0 = 0
    sortcx = sort(cx)
    if sortcx[1] < 1
        px0 = abs(sortcx[1]-1)
    end
    if sortcx[end] > sx
        if px0 < (sortcx[end]-sx)
            px0 = (sortcx[end]-sx)
        end
    end
    return px0
end

"""
    boxsmooth!(out::AbstractArray, arr::AbstractArray, tot::Array{T,1}, widx::Int, widy::Int)

Boxcar smooths an input image (or paddedview) `arr` with window size `widx` by
`widy`. We pass the original image size `sx` and `sy` to help handle image views.

# Arguments:
- `out::AbstractArray`: preallocated output array for the boxcar smoothed image
- `arr::AbstractArray`: input array for which boxcar smoothing is computed (generally paddedview)
- `tot::Array{T,1}`: preallocated array to hold moving sums along 1 dimension
- `widx::Int`: size of boxcar smoothing window in x
- `widy::Int`: size of boxcar smoothing window in y
"""
function boxsmooth!(out::AbstractArray, arr::AbstractArray, tot::Array{T,1}, widx::Int, widy::Int) where T
    (sx, sy) = size(arr)
    Δx = (widx-1)÷2
    Δy = (widy-1)÷2
    for j=1:(sy-widy+1)
        if (j==1)
            for n = 1:widy
                @simd for m = 1:sx
                    @inbounds tot[m] += arr[m,n]
                end
            end
        else
            @simd for m = 1:sx
                @inbounds tot[m] += arr[m,j+widy-1]-arr[m,j-1]
            end
        end
        tt=zero(eltype(out))
        for i=1:(sx-widx+1)
            if (i==1)
                @simd for n=1:widx
                    @inbounds tt += tot[n]
                end
            else
                @inbounds tt += tot[i+widx-1]-tot[i-1]
            end
            @inbounds out[i,j] = tt
        end
    end
end

"""
    cov_avg!(bimage, ism, bism, in_image; Np::Int=33, widx::Int=129, widy::Int=129, ftype::Int=32)

Key function for constructing the (shifted and multiplied) versions of the input image used to quickly
estimate the local covariance matrix at a large number of locations. The main output is in the preallocated
`bism` which is used as an input to `build_cov!`.

# Arguments:
- `bimage`: preallocated output array for the boxcar smoothed unshifted image
- `ism`: preallocated intermediate array for the input image times itself shifted
- `bism`: preallocated output array to store boxcar-smoothed image products for all shifts
- `in_image`: input image the local covariance of which we want to estimate

# Keywords:
- `Np::Int`: size of local covariance matrix in pixels (default 33)
- `widx::Int`: width of boxcar window in x which determines size of region used for samples for the local covariance estimate (default 129)
- `widy::Int`: width of boxcar window in y which determines size of region used for samples for the local covariance estimate (default 129)
- `ftype::Int`: determine the Float precision, 32 is Float32, otherwise Float64
"""
function cov_avg!(bimage, ism, bism, in_image; Np::Int=33, widx::Int=129, widy::Int=129, ftype::Int=32)
    if ftype == 32
        T = Float32
    else
        T = Float64
    end

    Δx = (widx-1)÷2
    Δy = (widy-1)÷2
    halfNp = (Np-1) ÷ 2

    (sx1, sy1) = size(in_image)
    tot = zeros(T,sx1);
    boxsmooth!(bimage,in_image,tot,widx,widy)
    # loop over shifts
    for dc=0:Np-1       # column shift loop
        for dr=1-Np:Np-1   # row loop, incl negatives
            if (dr < 0) & (dc == 0)
                continue
            end
            # ism = image, shifted and multipled
            @inbounds ism .= in_image .* ShiftedArrays.circshift(in_image,(-dr, -dc))
            fill!(tot,0)
            boxsmooth!(view(bism,:,:,dr+Np,dc+1),ism,tot,widx,widy) # bism = boxcar(ism)
        end
    end
    return
end


text:
#tests the framework of components and connections without actually running the model

using Base.Test
using Mimi

@defcomp A begin
  varA = Variable(index=[time])
  parA = Parameter(index=[time])
end

@defcomp B begin
  varB = Variable(index=[time])
  parB = Parameter(index=[time])
end

m = Model()
setindex(m, :time, [2015:5:2110])

addcomponent(m, A)
addcomponent(m, B, before=:A)

connectparameter(m, :A, :parA, :B, :varB)

@test length(m.components)==2
@test length(m.connections)==1
@test Mimi.get_connections(m, :A, :incoming)[1].source_component_name == :B
@test length(Mimi.get_connections(m, :B, :incoming)) == 0
@test Mimi.get_connections(m, :B, :outgoing)[1].target_component_name == :A
@test length(Mimi.get_connections(m, :A, :all)) == 1


text:
using Rainflow
using Test
using Random

Random.seed!(1234)
signal = randn(100)
ext, t = sort_peaks(signal)
result = count_cycles(ext,t)
bins = sum_cycles(result,13,1)
bins1 = sum_cycles(result, collect(range(0, 100, length=13+1)), collect(range(0, 100, length=1+1)))
bins2 = sum_cycles(result, range(0, 100, length=13+1), range(0, 100, length=1+1))
@test bins==bins1==bins2
@test dropdims(bins,dims=2) == [4.0, 6.0, 7.5, 3.0, 2.5, 5.0, 0.0, 2.5, 1.0, 0.0, 0.5, 0.0, 1.0]


text:
###############################################################################
#
# Function to find zero in the target function
#
###############################################################################
"""
Function to find the first root that gives a target function result of zero. If
    the root does not exist, the function returns the point where the target
    function is most close to zero.

$(METHODLIST)

"""
function find_zero end

"""
This method uses [`BisectionMethod`](@ref) method:

    find_zero(f::Function,
              ms::BisectionMethod{FT},
              tol::Union{ResidualTolerance{FT}, SolutionTolerance{FT}};
              stepping::Bool = false
    ) where {FT<:AbstractFloat}

Returns the solution where target function is zero, given
- `f` Function to solve
- `ms` [`BisectionMethod`](@ref) type method struct
- `tol` [`ResidualTolerance`](@ref) or [`SolutionTolerance`](@ref) type
    tolerance struct
- `stepping` Optional. If true, save the optimization steps to the history
    field in method struct.
"""
find_zero(f::Function,
          ms::BisectionMethod{FT},
          tol::Union{ResidualTolerance{FT}, SolutionTolerance{FT}};
          stepping::Bool = false
) where {FT<:AbstractFloat} =
(
    # create variable to store steps
    # _count for iterations
    _count::Int = 0;

    # calculate the values for x_min and x_max first
    _x_min = ms.x_min;
    _x_max = ms.x_max;
    _x_mid = _x_min   ;
    _y_min = f(_x_min);
    _y_max = f(_x_max);

    # record the history
    if stepping
        push!(ms.history, [_x_min, _y_min]);
        push!(ms.history, [_x_max, _y_max]);
    end;

    # if y_min and y_max are on the same side
    if _y_min * _y_max > 0
        if abs(_y_max) < abs(_y_min)
            _solution = _x_max;
        else
            _solution = _x_min;
        end;

    # if y_min is 0
    elseif _y_min==0
        _solution = _x_min;

    # if y_min and y_max are on different side, or one of them is 0
    else
        while true
            _x_mid = (_x_min + _x_max) / 2;
            _y_mid = f(_x_mid);

            # record the history
            if stepping
                push!(ms.history, [_x_mid, _y_mid]);
            end;

            # if difference is lower than the tolerance
            if if_break(tol, _x_min, _x_max, _y_mid, _count)
                break
            end;

            # if y_mid is on the same side with one of them
            if _y_mid * _y_min > 0
                _x_min = _x_mid;
                _y_min = _y_mid;
            elseif _y_mid * _y_max > 0
                _x_max = _x_mid;
                _y_max = _y_mid;

            # if y_mid = 0 and one of y_min and y_max is 0
            elseif (_y_min==0) && (_y_mid==0)
                _x_min = _x_mid;
                _y_min = _y_mid;
            elseif (_y_max==0) && (_y_mid==0)
                _x_max = _x_mid;
                _y_max = _y_mid;

            # else set x_max to x_mid (takes the fisrt one)
            else
                _x_max = _x_mid;
                _y_max = _y_mid;
            end;

            # _count ++
            _count += 1
        end;

        _solution = _x_mid;
    end;

    # return results
    return _solution
)

"""
This method uses [`NewtonBisectionMethod`](@ref) method:

    find_zero(f::Function,
              ms::NewtonBisectionMethod{FT},
              tol::Union{ResidualTolerance{FT}, SolutionTolerance{FT}};
              stepping::Bool = false
    ) where {FT<:AbstractFloat}

Returns the solution where target function is zero, given
- `f` Function to solve
- `ms` [`NewtonBisectionMethod`](@ref) type method struct
- `tol` [`ResidualTolerance`](@ref) or [`SolutionTolerance`](@ref) type
    tolerance struct
- `stepping` Optional. If true, save the optimization steps to the history
    field in method struct.
"""
find_zero(f::Function,
          ms::NewtonBisectionMethod{FT},
          tol::Union{ResidualTolerance{FT}, SolutionTolerance{FT}};
          stepping::Bool = false
) where {FT<:AbstractFloat} =
(
    # _count for iterations
    _count::Int = 0;

    # calculate the values for x_min and x_max first
    _x_dif = tol.tol ;
    _x_min = ms.x_min;
    _x_max = ms.x_max;
    _y_min = f(_x_min);
    _y_max = f(_x_max);

    # record the history
    if stepping
        push!(ms.history, [_x_min, _y_min]);
        push!(ms.history, [_x_max, _y_max]);
    end;

    # if y_min and y_max are on the same side
    if _y_min * _y_max > 0
        if abs(_y_max) < abs(_y_min)
            _solution = _x_max;
        else
            _solution = _x_min;
        end;

    # if y_min is 0
    elseif _y_min==0
        _solution = _x_min;

    # if y_min and y_max are on different side, or one of them is 0
    else
        _x_ntr = ms.x_ini;
        while true
            _y_ntr = f(_x_ntr);

            # record the history
            if stepping
                push!(ms.history, [_x_ntr, _y_ntr]);
            end;

            # if difference is lower than the tolerance
            if if_break(tol, _x_min, _x_max, _y_ntr, _count)
                break
            end;

            # if y_ntr is on the same side with one of them
            if _y_ntr * _y_min > 0
                _x_min = _x_ntr;
                _y_min = _y_ntr;
            elseif _y_ntr * _y_max > 0
                _x_max = _x_ntr;
                _y_max = _y_ntr;

            # if y_ntr = 0 and one of y_min and y_max is 0
            elseif (_y_min==0) && (_y_ntr==0)
                _x_min = _x_ntr;
                _y_min = _y_ntr;
            elseif (_y_max==0) && (_y_ntr==0)
                _x_max = _x_ntr;
                _y_max = _y_ntr;

            # else set x_max to x_ntr (takes the fisrt one)
            else
                _x_max = _x_ntr;
                _y_max = _y_ntr;
            end;

            # update x_ntr using Newton Raphson
            _x_dx   = _x_ntr + _x_dif * 10;
            _y_dx   = f(_x_dx);
            _slope  = (_y_dx - _y_ntr) / (_x_dx - _x_ntr);
            _x_ntr -= _y_ntr / _slope;

            # set x_ntr in the correct range
            if (_x_ntr>=_x_max) || (_x_ntr<=_x_min) || isnan(_x_ntr)
                _x_ntr = (_x_min + _x_max) / 2;
            end;

            # _count ++
            _count += 1;
        end;

        _solution = _x_ntr;
    end;

    # return results
    return _solution
)

"""
This method uses [`NewtonRaphsonMethod`](@ref) method:

    find_zero(f::Function,
              ms::NewtonRaphsonMethod{FT},
              tol::Union{ResidualTolerance{FT}, SolutionTolerance{FT}};
              stepping::Bool = false
    ) where {FT<:AbstractFloat}

Returns the solution where target function is zero, given
- `f` Function to solve
- `ms` [`NewtonRaphsonMethod`](@ref) type method struct
- `tol` [`ResidualTolerance`](@ref) or [`SolutionTolerance`](@ref) type
    tolerance struct
- `stepping` Optional. If true, save the optimization steps to the history
    field in method struct.
"""
find_zero(f::Function,
          ms::NewtonRaphsonMethod{FT},
          tol::Union{ResidualTolerance{FT}, SolutionTolerance{FT}};
          stepping::Bool = false
) where {FT<:AbstractFloat} =
(
    # _count for iterations
    _count::Int = 0;

    # calculate the values for x_min and x_max first
    _x_dif = tol.tol;

    # find the solution
    _x_lst = FT(Inf);
    _x_ntr = ms.x_ini;
    while true
        _y_ntr = f(_x_ntr);

        # record the history
        if stepping
            push!(ms.history, [_x_ntr, _y_ntr]);
        end;

        # if difference is lower than the tolerance
        if if_break(tol, _x_lst, _x_ntr, _y_ntr, _count)
            break
        end;

        # update x_ntr using Newton Raphson
        _x_dx   = _x_ntr + _x_dif;
        _y_dx   = f(_x_dx);
        _slope  = (_y_dx - _y_ntr) / (_x_dx - _x_ntr);
        _x_lst  = _x_ntr;
        _x_ntr -= _y_ntr / _slope;
        _x_dif  = tol.tol / abs(_slope);

        # _count ++
        _count += 1;
    end;

    # return results
    return _x_ntr
)

"""
This method uses [`ReduceStepMethod`](@ref) method:

    find_zero(f::Function,
              ms::ReduceStepMethod{FT},
              tol::Union{ResidualTolerance{FT}, SolutionTolerance{FT}};
              stepping::Bool = false
    ) where {FT<:AbstractFloat}

Returns the solution where target function is zero, given
- `f` Function to solve
- `ms` [`ReduceStepMethod`](@ref) type method struct
- `tol` [`ResidualTolerance`](@ref) or [`SolutionTolerance`](@ref) type
    tolerance struct
- `stepping` Optional. If true, save the optimization steps to the history
    field in method struct.
"""
find_zero(f::Function,
          ms::ReduceStepMethod{FT},
          tol::Union{ResidualTolerance{FT}, SolutionTolerance{FT}};
          stepping::Bool = false
) where {FT<:AbstractFloat} =
(
    # _count iterations
    _count = 0;

    # define the initial step
    @unpack x_ini, x_max, x_min, Δ_ini = ms;

    # initialize the y
    _tar_x = x_ini;
    _tar_y = abs( f(_tar_x) );

    # record the history
    if stepping
        push!(ms.history, [_tar_x, _tar_y]);
    end;

    _new_x = x_min;
    _new_y = abs( f(x_min) );
    _new_y < _tar_y ? (_tar_x=_new_x; _tar_y=_new_y;) : nothing;

    # record the history
    if stepping
        push!(ms.history, [_new_x, _new_y]);
    end;

    _new_x = x_max;
    _new_y = abs( f(x_max) );
    _new_y < _tar_y ? (_tar_x=_new_x; _tar_y=_new_y;) : nothing;

    # record the history
    if stepping
        push!(ms.history, [_new_x, _new_y]);
    end;

    # find the solution
    _Δx::FT = Δ_ini;
    while true
        # 1. increase the x by Δx till tar_y is bigger
        _count_inc = 0
        while true
            _new_x = _tar_x + _Δx;
            _new_x > x_max ? break : nothing;
            _new_y = abs( f(_new_x) );
            _new_y <  _tar_y ? (_tar_x=_new_x; _tar_y=_new_y;) : break;

            # record the history
            if stepping
                push!(ms.history, [_new_x, _new_y]);
            end;

            _count_inc += 1;
            _count += 1;
        end;

        # 2. decrease the x by Δx till tar_y is bigger
        _count_dec = 0
        while _count_inc == 0
            _new_x = _tar_x - _Δx;
            _new_x < x_min ? break : nothing;
            _new_y = abs( f(_new_x) );
            _new_y <= _tar_y ? (_tar_x=_new_x; _tar_y=_new_y;) : break;

            # record the history
            if stepping
                push!(ms.history, [_new_x, _new_y]);
            end;

            _count_dec += 1;
            _count += 1;
        end;

        # 3. if break
        if if_break(tol, FT(0), _Δx, _tar_y, _count)
            # record the history
            if stepping
                push!(ms.history, [_tar_x, _tar_y]);
            end;

            break
        end;

        # 4. if no update, then 10% the Δx
        if _count_inc + _count_dec == 0
            _Δx /= 10;
        end;
    end;

    return _tar_x
)


text:
module valueFunctions
function logit(x, fixedParam)
    return 1 ./ (1.0 .+ exp.(-fixedParam.lambda .* x))
end
#Functions for updating male value functions
#Male No Child
function getPmNC(x, y,v0,fxp,par,Pi,gem,SmNC,UmNC,PmD_hat,PmPL_hat,PmNC_hat,epsm,alpha,output,tau)
    expectedPmNC = 0.0
    dmNC1 = SmNC[x, :] .- SmNC[x, y]
    #Calculate the value function if no hc increment in the next round
    val1 = fxp.deltam_NC * (Pi[y] + UmNC[x]) +
           fxp.lambdae * par.sigma * sum(v0 .* dmNC1 .* logit(dmNC1, fxp)) * (SmNC[x, y] > 0) +
           fxp.gamma * PmD_hat[x, y] + fxp.chi * PmPL_hat[x, y] +
           (1 - fxp.deltam_NC - fxp.gamma - fxp.chi) * PmNC_hat[x, y]
    #Check if it is possible for HC to grow
    if (x % 7 == 0)
        expectedPmNC = val1
    else
        dmNC2 = SmNC[x+1, :] .- SmNC[x+1, y]
        #Calculate the value function if hc increments in the next round
        val2 = fxp.deltam_NC * (Pi[y] + UmNC[x+1]) +
               fxp.lambdae * par.sigma * sum(v0 .* dmNC2 .* logit(dmNC2, fxp)) * (SmNC[x+1, y] > 0) +
               fxp.gamma * PmD_hat[x+1, y] + fxp.chi * PmPL_hat[x+1, y] +
               (1 - fxp.deltam_NC - fxp.gamma - fxp.chi) * PmNC_hat[x+1, y]
        #weight the expected value by probability of HC growth
        if (y % 7 == 0)
            expectedPmNC = gem[7] * val2 + (1 - gem[7]) * val1
        else
            expectedPmNC = gem[y%7] * val2 + (1 - gem[y%7]) * val1
        end
    end
    #return sum of value at current round and expected value for next perioid
    return (1 - tau[1]) * output[x, y] + epsm[x] * alpha[y] + fxp.beta * expectedPmNC
end

#Male Parental Leave
function getPmPL(x, y,Pi,fxp,UmPL,PmD_hat,PmYC_hat,PmPL_hat,epsm,alpha,output)
    val = fxp.deltam_YC * (Pi[y] + UmPL[x]) +
          fxp.gamma * PmD_hat[x, y] +
          fxp.etam * PmYC_hat[x, y] +
          (1 - fxp.deltam_YC - fxp.gamma - fxp.etam) * PmPL_hat[x, y]

    return fxp.Transfer * output[x, y] + epsm[x] * alpha[y] + fxp.beta * val
end

#Male Young Child
function getPmYC(x, y,v0,Pi,fxp,par,epsm,gem,tau,output,alpha,UmYC,SmYC,PmD_hat,PmPL_hat,PmYC_hat)
    expectedPmYC = 0.0
    dmYC1 = SmYC[x, :] .- SmYC[x, y]
    #Calculate the value function if no hc increment in the next round
    val1 = fxp.deltam_YC * (Pi[y] + UmYC[x]) +
           fxp.lambdae * par.sigma * sum(v0 .* dmYC1 .* logit(dmYC1, fxp)) * (SmYC[x, y] > 0) +
           fxp.gamma * PmD_hat[x, y] + fxp.chi * PmPL_hat[x, y] +
           (1 - fxp.deltam_YC - fxp.gamma - fxp.chi) * PmYC_hat[x, y]
    #Check if it is possible for HC to grow
    if (x % 7 == 0)
        expectedPmYC = val1
    else
        dmYC2 = SmYC[x+1, :] .- SmYC[x+1, y]
        #Calculate the value function if hc increments in the next round
        val2 = fxp.deltam_YC * (Pi[y] + UmYC[x+1]) +
               fxp.lambdae * par.sigma * sum(v0 .* dmYC2 .* logit(dmYC2, fxp)) * (SmYC[x+1, y] > 0) +
               fxp.gamma * PmD_hat[x+1, y] + fxp.chi * PmPL_hat[x+1, y] +
               (1 - fxp.deltam_YC - fxp.gamma - fxp.chi) * PmYC_hat[x+1, y]
        if (y % 7 == 0)
            expectedPmYC = gem[7] * val2 + (1 - gem[7]) * val1
        else
            expectedPmYC = gem[y%7] * val2 + (1 - gem[y%7]) * val1
        end
    end
    return (1 - tau[1]) * output[x, y] + epsm[x] * alpha[y] + fxp.beta * expectedPmYC
end

#Male Done with Children
function getPmD(x,y,v0,Pi,fxp,par,output,gem,tau,epsm,alpha,SmD,PmD_hat,UmD)
    expectedPmD = 0.0
    dmD1 = SmD[x, :] .- SmD[x, y]
    #Calculate the value function if no hc increment in the next round
    val1 = fxp.phi * Pi[y] + fxp.delta * (Pi[y] + UmD[x]) +
           (1 - fxp.phi - fxp.delta) * PmD_hat[x, y] +
           fxp.lambdae * par.sigma * sum(v0 .* dmD1 .* logit(dmD1, fxp)) * (SmD[x, y] > 0)
    #Check if it is possible for HC to grow
    if (x % 7 == 0)
        expectedPmD = val1
    else
        dmD2 = SmD[x+1, :] .- SmD[x+1, y]
        #Calculate the value function if hc increments in the next round
        val2 = fxp.phi * Pi[y] + fxp.delta * (Pi[y] + UmD[x+1]) +
               (1 - fxp.phi - fxp.delta) * PmD_hat[x+1, y] +
               fxp.lambdae * par.sigma * sum(v0 .* dmD2 .* logit(dmD2, fxp)) * (SmD[x+1, y] > 0)
        if (y % 7 == 0)
            expectedPmD = gem[7] * val2 + (1 - gem[7]) * val1
        else
            expectedPmD = gem[y%7] * val2 + (1 - gem[y%7]) * val1
        end
    end
    return (1 - tau[1]) * output[x, y] + epsm[x] * alpha[y] + fxp.beta * expectedPmD
end

#Functions for updating female value functions
#Female No Child
function getPfNC(x, y,v0,fxp,par,Pi,gef,SfNC,UfNC,PfD_hat,PfPL_hat,PfNC_hat,epsf,alpha,output,tau)
    expectedPfNC = 0.0
    dfNC1 = SfNC[x, :] .- SfNC[x, y]
    #Calculate the value function if no hc increment in the next round
    val1 = fxp.deltaf_NC * (Pi[y] + UfNC[x]) +
           fxp.lambdae * par.sigma * sum(v0 .* dfNC1 .* logit(dfNC1, fxp)) * (SfNC[x, y] > 0) +
           fxp.gamma * PfD_hat[x, y] + fxp.chi * PfPL_hat[x, y] +
           (1 - fxp.deltaf_NC - fxp.gamma - fxp.chi) * PfNC_hat[x, y]
    #Check if it is possible for HC to grow
    if (x % 7 == 0)
        expectedPfNC = val1
    else
        dfNC2 = SfNC[x+1, :] .- SfNC[x+1, y]
        #Calculate the value function if hc increments in the next round
        val2 = fxp.deltaf_NC * (Pi[y] + UfNC[x+1]) +
               fxp.lambdae * par.sigma * sum(v0 .* dfNC2 .* logit(dfNC2, fxp)) * (SfNC[x+1, y] > 0) +
               fxp.gamma * PfD_hat[x+1, y] + fxp.chi * PfPL_hat[x+1, y] +
               (1 - fxp.deltaf_NC - fxp.gamma - fxp.chi) * PfNC_hat[x+1, y]
        if (y % 7 == 0)
            expectedPfNC = gef[7] * val2 + (1 - gef[7]) * val1
        else
            expectedPfNC = gef[y%7] * val2 + (1 - gef[y%7]) * val1
        end
    end
    #return value at current round and expected value for next perioid
    return (1 - tau[1]) * output[x, y] + epsf[x] * alpha[y] + fxp.beta * expectedPfNC
end

#Female Parental Leave
function getPfPL(x, y,Pi,par,fxp,UfPL,PfD_hat,PfYC_hat,PfPL_hat,epsf,alpha,output)
    val = fxp.deltaf_YC * (Pi[y] + UfPL[x]) +
          fxp.gamma * PfD_hat[x, y] +
          fxp.etaf * PfYC_hat[x, y] +
          (1 - fxp.deltaf_YC - fxp.gamma - fxp.etaf) * PfPL_hat[x, y]

    return fxp.Transfer * output[x, y] + (par.M + epsf[x]) * alpha[y] + fxp.beta * val
end

#Female Young Child
function getPfYC(x, y,v0,Pi,fxp,par,epsf,gef,tau,output,alpha,UfYC,SfYC,PfD_hat,PfPL_hat,PfYC_hat)
    expectedPfYC = 0.0
    dfYC1 = SfYC[x, :] .- SfYC[x, y]
    #Calculate the value function if no hc increment in the next round
    val1 = fxp.deltaf_YC * (Pi[y] + UfYC[x]) +
           fxp.lambdae * par.sigma * sum(v0 .* dfYC1 .* logit(dfYC1, fxp)) * (SfYC[x, y] > 0) +
           fxp.gamma2 * PfD_hat[x, y] + fxp.chi * PfPL_hat[x, y] +
           (1 - fxp.deltaf_YC - fxp.gamma - fxp.chi) * PfYC_hat[x, y]
    #Check if it is possible for HC to grow
    if (x % 7 == 0)
        expectedPfYC = val1
    else
        dfYC2 = SfYC[x+1, :] .- SfYC[x+1, y]
        #Calculate the value function if hc increments in the next round
        val2 = fxp.deltaf_YC * (Pi[y] + UfYC[x+1]) +
               fxp.lambdae * par.sigma * sum(v0 .* dfYC2 .* logit(dfYC2, fxp)) * (SfYC[x+1, y] > 0) +
               fxp.gamma2 * PfD_hat[x+1, y] + fxp.chi * PfPL_hat[x+1, y] +
               (1 - fxp.deltaf_YC - fxp.gamma - fxp.chi) * PfYC_hat[x+1, y]
        if (y % 7 == 0)
            expectedPfYC = gef[7] * val2 + (1 - gef[7]) * val1
        else
            expectedPfYC = gef[y%7] * val2 + (1 - gef[y%7]) * val1
        end
    end
    return (1 - tau[1]) * output[x, y] + (par.M + epsf[x]) * alpha[y] + fxp.beta * expectedPfYC
end

#Female Done with Children
function getPfD(x,y,v0,Pi,fxp,par,output,gef,tau,epsf,alpha,SfD,PfD_hat,UfD)
    expectedPfD = 0.0
    dfD1 = SfD[x, :] .- SfD[x, y]
    #Calculate the value function if no hc increment in the next round
    val1 = fxp.phi * Pi[y] + fxp.delta * (Pi[y] + UfD[x]) +
           (1 - fxp.phi - fxp.delta) * PfD_hat[x, y] +
           fxp.lambdae * par.sigma * sum(v0 .* dfD1 .* logit(dfD1, fxp)) * (SfD[x, y] > 0)
    #Check if it is possible for HC to grow
    if (x % 7 == 0)
        expectedPfD = val1
    else
        dfD2 = SfD[x+1, :] .- SfD[x+1, y]
        #Calculate the value function if hc increments in the next round
        val2 = fxp.phi * Pi[y] + fxp.delta * (Pi[y] + UfD[x+1]) +
               (1 - fxp.phi - fxp.delta) * PfD_hat[x+1, y] +
               fxp.lambdae * par.sigma * sum(v0 .* dfD2 .* logit(dfD2, fxp)) * (SfD[x+1, y] > 0)
        if (y % 7 == 0)
            expectedPfD = gef[7] * val2 + (1 - gef[7]) * val1
        else
            expectedPfD = gef[y%7] * val2 + (1 - gef[y%7]) * val1
        end
    end
    return (1 - tau[1]) * output[x, y] + epsf[x] * alpha[y] + fxp.beta * expectedPfD
end
export getPmNC,getPmYC,getPmPL,getPmD,getPfNC,getPfYC,getPfPL,getPfD
end


text:
using Test
using CSV
using DataFrames
using RealNeuralNetworks.Neurons.Segments.Synapses

@testset "test Synapses module" begin
    df = CSV.read(joinpath(Base.@__DIR__,
                           "../asset/77625.pre.synapses.csv"))
    synapseList = Vector{Synapse{Float32}}()
    for row in DataFrames.eachrow(df)
        push!(synapseList, Synapse(row))
    end
    @test !isempty( synapseList )

    synapse = synapseList[1]
    @test Synapses.get_psd_size(synapse) > 0
    @test Synapses.get_psd_segmentation_id(synapse) > 0
    @test Synapses.get_psd_coordinate( synapse ) |> length > 0
    Synapses.get_psd_bounding_box( synapse )

    @test Synapses.get_pre_synaptic_segmentation_id( synapse ) > 0
    @test Synapses.get_pre_synaptic_coordinate( synapse ) |> length > 0
    @test Synapses.get_pre_synaptic_weight( synapse ) > 0.0

    @test Synapses.get_post_synaptic_segmentation_id( synapse ) > 0
    @test Synapses.get_post_synaptic_coordinate( synapse ) |> length > 0
    @test Synapses.get_post_synaptic_weight( synapse )  > 0.0
end


text:
module SIIPExamples

export print_struct
export print_tree
export notebook
export JuliaExamples
export PSYExamples
export PSIExamples
export PSDExamples

using Literate
using JSON3
using AbstractTrees
using InteractiveUtils
import IJulia

include("definitions.jl")
include("utils.jl")
include("launch.jl")

end # module


text:
## KL-Distribution Test
## ====================
holesizes = logspace(-1,-10,3)
problems = [Simplex(8,        #ndims
                    [:sample_distribution,:accumulative_KL, :total_time,],
                    5,        #nsamples
                    0.01,
                    mvuniform(-2,2,8))]    #holesize

mh_captures = [:samples]
all_splits = [rand_partial_split]

algorithms = [SigmaSMT(mh_captures, randvartype, sampler, nprocs, split)
  for nprocs = [1],
      randvartype = [DRealRandVar],
      split = all_splits,
      sampler = [Sigma.point_sample_mc]][:]

function kl()
  record(algorithms, problems;
         runname = "kl",
         prefix=benchdir,
         exceptions=false)
end

kl()

## Analysis
## This needs a rethink
# How do I want to use this in the end
# a = [run(analysis) or analysis in all_analyses]
# # layer...
DiscreteDistribution = Dict{Int, Float64}
#
"returns an output analysis of type O"
abstract Analysis{O}

"Data of relational data X vs Y"
immutable XY{X,Y}
  xdata::Vector{X}
  ydata::Vector{Y}
  xlabel::String
  ylabel::String
end

"""A KL Analysis computes the KL Divergence.
Takes in samples over time and produces output of KL vs time"""
immutable KLAnalysis <: Analysis{XY{Float64, Int}}
  f::XY{Float64, Int}
end

function KLAnalysis(allsamples::Vector{Float64}, truth::DiscreteDistribution)
  for i = 1:length(allsamples)
    samples = allsamples[1:i]
    v = vertex_distribution(samples,ndims,0.01)
    truth = groundtruth(ndims)
    push!(kls,KLsafe(truth,v))
  end
  kls
  XY{Float64,Int64}(samples, )
end

run(x::KLAnalysis) = x.XY

## Analyses
## ========

function count(ndims)
  results = DynamicAnalysis.all_records() |> r->DynamicAnalysis.where(r,:runname,j->j=="kl") |> r->DynamicAnalysis.where(r,:problem,j->j.ndims == ndims)
  @show size(results,1)

end

function accumulative_kl(i, ndims)
  records = DynamicAnalysis.all_records() |> r->DynamicAnalysis.where(r,:runname,j->j=="kl") |> r->DynamicAnalysis.where(r,:problem,j->j.ndims == ndims)
  results = records[:results]
  @show length(results)
  allsamples = Lens.get(results[i];lensname=:samples)
  kls = Float64[]
  for i = 1:length(allsamples)
    samples = allsamples[1:i]
    v = vertex_distribution(samples,ndims,0.01)
    truth = groundtruth(ndims)
    push!(kls,KLsafe(truth,v))
  end
  kls
end

function accumulative_mean_kls(dims)
  trials = [accumulative_kl(i, dims) for i = 1:count(dims)]
  mean_kls = Array(Float64, length(trials[1]))
  stds_kls = Float64[]
  for i = 1:length(trials[1])
    total = 0.0
    for j = 1:length(trials)
      total += trials[j][i]
    end
    mean_kls[i] = total/length(trials)
  end
  return mean_kls
end

function plot_kls(dims::Vector{Int})
  dfs = DataFrame[]
  for dim in dims
    akl = accumulative_mean_kls(dim)
    df = DataFrame(x=1:length(akl), y = akl, label = "Dim:$dim")
    push!(dfs,df)
  end
  Gadfly.plot(vcat(dfs...), x="x", y="y", color="label",Gadfly.Geom.line,
     Guide.xlabel("Number of Samples"),
     Guide.ylabel("KL Divergence"))
end

# splitkey = [rand_partial_split => "rand", weighted_mid_split => "mid", weighted_partial_split => "partial"]

# # # Analysis
# # ## KL
# using Lens
# using SQLite
# using Gadfly
# using DataFrames
# q = all_records() |> r->Lens.where(r,"runname",j->j=="kl")
# g = group(q)

# plot_db(y = tree, x = :iterate, groupby=)

# function plot_kl2(groups)
#   dfs = DataFrame[]
#   i = 1
#   for (k,v) in groups
#     results = Result[x[7] for x in v]
#     collated = Lens.collate(results, :samples, :x1)
#     @show length(collated)
#     samples = collated
#     ndims = k[2].ndims
#     holesize = k[2].holesize
#     gt = groundtruth(ndims)
#     akl = accumulative_KL(samples, ndims, gt,holesize)
#     df = DataFrame(x=1:length(akl), y = akl, label = "SMT$i")
#     push!(dfs,df)
#     i += 1
#   end
# #   return  dfs
#   plot(vcat(dfs...), x="x", y="y", color="label",Geom.line,
#        Guide.xlabel("Number of Samples"),
#        Guide.ylabel("KL Divergence"))
# end

# plot_kl2(g)

# draw(PDF("/home/zenna/Dropbox/writing/conferences/icml2015/figures/smt1.pdf",4*1inch,3*1inch), smtplot)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""ResourceQuotaSpec defines the desired hard limits to enforce for Quota.

    IoK8sApiCoreV1ResourceQuotaSpec(;
        hard=nothing,
        scopeSelector=nothing,
        scopes=nothing,
    )

    - hard::Dict{String, IoK8sApimachineryPkgApiResourceQuantity} : hard is the set of desired hard limits for each named resource. More info: https://kubernetes.io/docs/concepts/policy/resource-quotas/
    - scopeSelector::IoK8sApiCoreV1ScopeSelector : scopeSelector is also a collection of filters like scopes that must match each object tracked by a quota but expressed using ScopeSelectorOperator in combination with possible values. For a resource to match, both scopes AND scopeSelector (if specified in spec), must be matched.
    - scopes::Vector{String} : A collection of filters that must match each object tracked by a quota. If not specified, the quota matches all objects.
"""
mutable struct IoK8sApiCoreV1ResourceQuotaSpec <: SwaggerModel
    hard::Any # spec type: Union{ Nothing, Dict{String, IoK8sApimachineryPkgApiResourceQuantity} } # spec name: hard
    scopeSelector::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ScopeSelector } # spec name: scopeSelector
    scopes::Any # spec type: Union{ Nothing, Vector{String} } # spec name: scopes

    function IoK8sApiCoreV1ResourceQuotaSpec(;hard=nothing, scopeSelector=nothing, scopes=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ResourceQuotaSpec, Symbol("hard"), hard)
        setfield!(o, Symbol("hard"), hard)
        validate_property(IoK8sApiCoreV1ResourceQuotaSpec, Symbol("scopeSelector"), scopeSelector)
        setfield!(o, Symbol("scopeSelector"), scopeSelector)
        validate_property(IoK8sApiCoreV1ResourceQuotaSpec, Symbol("scopes"), scopes)
        setfield!(o, Symbol("scopes"), scopes)
        o
    end
end # type IoK8sApiCoreV1ResourceQuotaSpec

const _property_map_IoK8sApiCoreV1ResourceQuotaSpec = Dict{Symbol,Symbol}(Symbol("hard")=>Symbol("hard"), Symbol("scopeSelector")=>Symbol("scopeSelector"), Symbol("scopes")=>Symbol("scopes"))
const _property_types_IoK8sApiCoreV1ResourceQuotaSpec = Dict{Symbol,String}(Symbol("hard")=>"Dict{String, IoK8sApimachineryPkgApiResourceQuantity}", Symbol("scopeSelector")=>"IoK8sApiCoreV1ScopeSelector", Symbol("scopes")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiCoreV1ResourceQuotaSpec }) = collect(keys(_property_map_IoK8sApiCoreV1ResourceQuotaSpec))
Swagger.property_type(::Type{ IoK8sApiCoreV1ResourceQuotaSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1ResourceQuotaSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ResourceQuotaSpec }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ResourceQuotaSpec[property_name]

function check_required(o::IoK8sApiCoreV1ResourceQuotaSpec)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ResourceQuotaSpec }, name::Symbol, val)
end


text:
using JFVM
using PyPlot
include("CRT.jl")

#Calculating error as a function of timesteps

# define the domain and constants
ϕ  = 0.4 # porosity
u  = 1.0/(3600*24) # m/s, eq. 1 m/day
D  = [1e-9, 1e-9, 1e-9] # m^2/s diffusivity
K_eq = 10.0
k_f = 100.0
k_b = k_f/K_eq
nu_A = 1.0 # stoichiometric coefficient of component A
nu_B = 1.0 # stoichiometric coefficient of component B
nu_C = 1.0 # stoichiometric coefficient of component C
a = 1.0 #order of reaction
b = 1.0 #order of reaction
c = 1.0 #order of reaction
N_comp = 3 # number of components

# Defining time-stepping parameters
tol = 1e-10 # convergence tolerance

# boundary conditions
C_bc = [1.0, 0.5, 0.0] # boundary concentration
C_0 = zeros(N_comp) # initial concentration

NX = 100:100:1000
time_array = zeros(5, length(NX))

for i=1:length(NX)
  Lx = 1.0 # [m]
  Nx = NX[i]  # number of grids
  m = createMesh1D(Nx, Lx)  # mesh

  #Time-stepping parameters
  dt = Lx/(u/ϕ)/Nx
  T = Nx/2*dt

  tic()
  C_SNIA = CRT.CRT_SNIA(m, Nx, ϕ, u, D, K_eq, nu_A, nu_B, nu_C, C_bc, C_0, N_comp, dt, T)
  time_array[1, i] = toc()

  tic()
  C_SIA = CRT.CRT_SIA(m, Nx, ϕ, u, D, K_eq, nu_A, nu_B, nu_C, C_bc, C_0, N_comp, dt, T, tol)
  time_array[2, i] = toc()

  tic()
  C_eps = CRT.CRT_GIA_eps(m, Nx, ϕ, u, D, K_eq, nu_A, nu_B, nu_C, C_bc, C_0, N_comp, dt, T, tol)
  time_array[3, i] = toc()

  tic()
  C_GIA_kin_eq = CRT.CRT_GIA_kin_eq(m, Nx, ϕ, u, D, K_eq, nu_A, nu_B, nu_C, a, b, c, C_bc, C_0, N_comp, dt, T, tol)
  time_array[4, i] = toc()

  tic()
  C_GIA_kin_fb = CRT.CRT_GIA_kin_fb(m, Nx, ϕ, u, D, k_f, k_b, nu_A, nu_B, nu_C, a, b, c, C_bc, C_0, N_comp, dt, T, tol)
  time_array[5, i] = toc()

  println(i)
end

figure()
plot(NX, time_array[1,:])
plot(NX, time_array[2,:])
plot(NX, time_array[3,:])
plot(NX, time_array[4,:])
plot(NX, time_array[5,:])
grid("on")
legend(["SNIA", "SIA", "GIA ("L"$\epsilon$"")", "GIA ("L"$K_{eq}$"")", "GIA ("L"$k_fk_b$"")"])
xlabel(L"$N_x$")
ylabel("Time [s]""")
savefig("time_grid_plot.png")

#CRT_SIA title:
#L"$\phi\frac{\partial C_{\alpha}}{\partial t}+\nabla\cdot(uC_{\alpha})-\nabla\cdot(D_{\alpha}\nabla C_{\alpha})+\nu_{\alpha}\epsilon=0,\quad K_{eq}=$""$K_eq"" (SIA)"

#CRT_SNIA title:
#L"$\phi\frac{\partial C_{\alpha}}{\partial t}+\nabla\cdot(uC_{\alpha})-\nabla\cdot(D_{\alpha}\nabla C_{\alpha})+\nu_{\alpha}\epsilon=0,\quad K_{eq}=$""$K_eq"" (SNIA)"

#CRT_GIA_eps title:
#L"$\phi\frac{\partial C_{\alpha}}{\partial t}+\nabla\cdot(uC_{\alpha})-\nabla\cdot(D_{\alpha}\nabla C_{\alpha})+\nu_{\alpha}\epsilon=0,\quad K_{eq}=$""$K_eq"" (FIM)"

#CRT_GIA_kin_eq title:
#L"$\phi\frac{\partial C_{\alpha}}{\partial t}+\nabla\cdot(uC_{\alpha})-\nabla\cdot(D_{\alpha}\nabla C_{\alpha})+\nu_{\alpha}(K_{eq}C_A^aC_B^b-C_C^c)=0,\quad K_{eq}=$""$K_eq"" (FIM)"

#CRT_GIA_kin_fb title:
#L"$\phi\frac{\partial C_{\alpha}}{\partial t}+\nabla\cdot(uC_{\alpha})-\nabla\cdot(D_{\alpha}\nabla C_{\alpha})+\nu_{\alpha}(k_{f}C_A^aC_B^b-k_{b}C_C^c)=0,\quad K_{eq}=$""$K_eq"" (FIM)"

max_stars_repo_name: FrederikLehn/coupled-solvers

max_issues_repo_name: FrederikLehn/coupled-solvers

max_forks_repo_name: FrederikLehn/coupled-solvers