text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    localanisotropies(Geometric, searcher; simplify=true)

Extract `LocalAnisotropy` from a searcher object. Based on the spatial data and
the parameters from the searcher object, a group of neighbor points is collected
at each place. PCA is applied to their spatial coordinates, returning eigenvectors
that represent a local ellipse/ellipsoid aligned with the directions with more
clustered points. The magnitude is the ratio of the PCA eigenvalues. It's useful
to extract local anisotropies from 3-D surface points. If `simplify = true` and
data is 3-D, only the minimum direction is extracted from PCA; the other axes
are calculated (main direction will be always along max dip direction and the
intermediate axis will be orthogonal to the others). This is default for 3-D.

## Example

```julia
searcher = KNearestSearch(pts3dsurface, 10)
prelpars = localanisotropies(Geometric, searcher) # extract local pars at surface
lpars = idwpars(prelpars, searcher, grid) # interpolate local pars to a grid
```
"""
function localanisotropies(::Type{Geometric}, searcher::NeighborSearchMethod;
	simplify::Bool=true)
	D = searcher.domain
	X = coords(D)
	N, len = size(X)

    quat = Array{Quaternion}(undef,len)
    m = Array{Vector}(undef,len)

    Threads.@threads for i in 1:len
		neighids = search(centroid(D,i), searcher)
		λ, v = pca(view(X,:,neighids), simplify)

		det(v) < 0 && (v = Diagonal(SVector{3}([-1,1,1])) * v)

        q = dcm_to_quat(v)
        quat[i] = q
        m[i] = λ/λ[1]
    end

	# deal with -1 eigvals
	m = reduce(hcat,vec(m))
	posm = m .> 0
	if sum(posm) < N*len
		for d in 2:N
			posd = view(posm,d,:)
			sum(posd) == len && continue
			minx = minimum(view(m,d,findall(posd)))
			m[d,findall(.!posd)] .= minx
		end
	end

    LocalAnisotropy(quat, m)
end

function pca(X, simplify)
	N = size(X,1)
	M = fit(PCA, X, maxoutdim=N, pratio=1)
	λ = principalvars(M) ./ principalvars(M)[1]
	nv = length(λ)
	v = N == 3 ? projection(M) : vcat(projection(M),[0 0 1][1:nv])

	if nv == 1
		vx = [-v[2,1]; v[1,1]; v[3,1]]
		v = hcat(v,vx,cross(v[:,1],vx))
		append!(λ,[λ[1],λ[1]])
	elseif nv == 2 && N==3
		v = hcat(v,cross(v[:,1],v[:,2]))
		push!(λ,minimum(λ))
	end

	if simplify && N==3
		# use maxdip vector and cross(maxdip,normal) as the main vectors
		abs(v[3,3]) > 1 && (v[3,3] = round(v[3,3], digits=1))
		az, dp = [atan(v[1,3], v[2,3]), -asin(v[3,3])]
		if dp < 0
			dp += pi/2
		else
			dp -= pi/2
		end
		v[:,1] .= [sin(az)*cos(dp), cos(az)*cos(dp), -sin(dp)]
		v[:,2] .= cross(v[:,1], v[:,3])
		λ[1:2] .= [1.0, 1.0]
	end

	λ[1:N], SMatrix{3,3}(v')
end


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ fca3485c-552a-4bf2-8cfd-bf7f548a6d1a
begin
	using PlutoUI
	PlutoUI.TableOfContents()
end

# ╔═╡ 9f0ba450-ca8a-11eb-2d40-594936524dda
md"""
# First steps to Julia
"""

# ╔═╡ 8972e6cf-5abd-4e34-9016-d316d3b94fd7
md"""
# Comments
"""

# ╔═╡ 183e7365-5aa3-4ebe-8757-ff4ca66da186
# You can leave comments on a single line using the pound/hash key

#=

For multi-line comments,
use the '#= =#' sequence.

=#

# ╔═╡ ff7496bc-5dd8-4b31-a8e2-910f4df257fb
md"""
# Variables and expressions

The assignment operator `=` binds a variable name to the data.

To see the content in the data:

- `@show exp` to show the expression and its result
- [`println(stuff)`](https://docs.julialang.org/en/v1/base/io-network/#Base.println) for regular print
- [`@printf`](https://docs.julialang.org/en/v1/stdlib/Printf/) for C-like formatted output.
- Inline display for the last expression
- A semicolon `;` at the end will suppress the output

There are some differences and limitations in `Pluto` notebooks:

- The results are on the top of the code cells instead at the bottom (in Jupyter notebooks)
- It does not allow *reassignment* nor *multiple expression per cell*.
- [Standard output](https://en.wikipedia.org/wiki/Standard_streams#Standard_output_(stdout)) is not displayed in the Pluto notebook.
- `Distributed.jl` (multiprocessing) is not avialable since Pluto notebooks run on worker processes.

If you really want `println()` output in Pluto notebooks, use `with_terminal()` in the `PlutoUI.jl` package.
"""

# ╔═╡ fb3654b2-a35e-4db8-8b8f-d8c8e337d14d
# Integer, 64 bit
x = 1

# ╔═╡ 6bde72a8-e159-4415-85de-c08bc186cff6
# Floating-point, 64-bit
y = 1.0

# ╔═╡ dedfaa4d-859a-41c0-b92f-de63e02bcf7a
# Complex number
z = 3 + 5im

# ╔═╡ 22a0a72f-d5e3-440a-b1ea-a8d667f12863
α = 3.74 # Unicode variable names — type it by "\alpha<tab>"

# ╔═╡ d32b9ce8-8aed-46f6-bdf8-22868710cfd7
# Strings (Text) are surrounded by double quotes
s = "Julia"

# ╔═╡ 764b4d0a-f643-4f59-9060-bda0be7a6afe
# Rational type
station = 9 + 3//4

# ╔═╡ 3b832f33-a264-4b10-ba4a-14b9f8a3e212
# Characters are surrounded by single quotes
c = ' '

# ╔═╡ bb231574-ce25-4396-a32b-e5ace12937ec
# Constants will emit a warning / error if you try to change it
const theUltimateAnswer = 42

# ╔═╡ a2d2a8c2-6afd-43fb-bb41-7a605f4c6d74
# Use with_terminal() to capture println() outputs
with_terminal() do
	println("Hello World")
	println("Hello ", s)
	println(1, 2, 3)
end

# ╔═╡ 2d1eeab7-93ae-412f-a7f1-96f3c023094a
with_terminal() do
	# @show will print x = val
	@show x	1+1 2-2 3*3
end

# ╔═╡ 7a35ab55-173c-4468-8ca0-0e510d96e0d4
md"""
# Data types
"""

# ╔═╡ 73c1c8ce-c227-42be-b5a3-d88745a3e60e
with_terminal() do
	# `typeof(val)` shows the data type.
	@show typeof(x)
	@show typeof(y)
	@show typeof(c)
	@show typeof(s)
	@show typeof(z)
	@show typeof(1//2)
end

# ╔═╡ 48793eaf-cc99-4fac-99d8-d59e2c07b335
# convert(T,x) to concert x to type T
typeof(convert(Float64, x))

# ╔═╡ 24a0d4ac-4f72-4127-8afa-87bcbb16c876
# Often you can use Type(x)
typeof(Float64(x))

# ╔═╡ 0bd93bdc-901d-438b-965e-416032f8dd88
# Or even shorthand function
typeof(float(x))

# ╔═╡ 6539f01d-0b1a-42b9-ac03-882ce9f419f6
md"""
# Compound expressions

- A [begin block](https://docs.julialang.org/en/v1/base/base/#begin) `begin` ... `end` squashes a multiline expressions into one.
- A [let block](https://docs.julialang.org/en/v1/manual/variables-and-scoping/#Let-Blocks) `let` ... `end` is similar to a begin block but Variables in the let block will be discarded outside the block.
"""

# ╔═╡ aaa4115e-4dd9-4f48-9a79-27f209970dce
# Compound expression: begin ... end
# Variables in the begin block will be kept outside the block
begin
	a1 = 2
	a2 = 35
	a1 = a1 * a2
end

# ╔═╡ 966de6c6-c941-4493-a3eb-51ee7064e538
with_terminal() do
	@show a1 a2
end

# ╔═╡ 56b16f1b-6809-48f5-85dd-72cf86e7d416
# Compound expression: let ... end.
# Variables in the let block will be discarded outside the block
let
	x = 1
	y = 2
	z = 3

	x + y * z
end

# ╔═╡ b7af6203-833c-4526-bf40-5635a1aac795
md"""
# Math expressions

Julia supports basic arithmatic operations and essentail math functions by default.

## Basic arithmetic
"""

# ╔═╡ 1b2e41b0-3e5e-4e1b-ae4b-820f27c1cc70
a = 2

# ╔═╡ 3636e2ba-f034-4b73-ba78-48fd65198e7a
b = 3

# ╔═╡ 3929c584-f12b-47f6-913b-4b0ab2c06994
a + b

# ╔═╡ 5478fc76-c2b6-41ca-8648-b8437caacce4
a - b

# ╔═╡ 8551ddbe-096d-4c14-bf81-9a2a4e648116
a * b

# ╔═╡ bfa38d3f-e655-4b13-be22-88fab5fbdfc7
b / a

# ╔═╡ de20682a-1410-40c5-b115-5b5470de870e
typeof(b / a)

# ╔═╡ 5a80edaf-bbd0-45b0-af9d-c135749318cc
b // a

# ╔═╡ b04a59ff-9def-4799-b47e-1e1e948b2e3f
typeof(b // a)

# ╔═╡ e5e519d3-1d44-4490-a97b-6fb397741a21
# integer division: \div<tab>`, the same as div(a, b)`
a÷b

# ╔═╡ 3d4846ab-b537-416b-bf7e-364215e7bd7e
typeof(a÷b)

# ╔═╡ 826c913c-a88b-4e4b-8e0a-83888c8df767
b % a

# ╔═╡ 63dda274-30f3-489d-b7fd-eae5266f097c
a^a

# ╔═╡ 05f68beb-d349-4d01-aa0d-6e879da8b14e
md"""
## Comparison

Returns a Boolean value (true / false)
"""

# ╔═╡ 968bad8a-6df8-4b92-9eff-f808e0220bef
a < 1

# ╔═╡ be6283eb-1265-42f6-9eae-df35edaf9306
b > 4

# ╔═╡ c5db3a91-0fba-4abd-9fb7-3d7b4b7f53e8
a <= b

# ╔═╡ 6b354bd0-1e53-45f6-a783-6b560019e8d7
a != b

# ╔═╡ 116c21df-d466-4723-9513-ee83a73fc4bf
a == b + 1

# ╔═╡ a63d784c-71e5-4f74-8186-c6887f829900
# Chained comparisons are supported
1 < a <= b

# ╔═╡ 5f949ee8-3694-459f-be29-63e1b4371d30
# Approximation operator `\approx` for floating point number equivalence
1e10 + 1.0 ≈ 1e10

# ╔═╡ 9644b2e4-b270-488b-a6d0-a4ab4c8aa73a
isapprox(1e10 + 1.0, 1e10)

# ╔═╡ f53e2490-f4b0-44dd-9ac9-907e60458faf
md"""
## Math functions

"""

# ╔═╡ 51bc09d3-a05c-46be-ab94-319ede805291
sin(0.5*π)

# ╔═╡ 7e9cecbe-8669-44c8-80f0-2b70a3739cc3
cos(0.5*π)

# ╔═╡ c026471c-3b69-4089-bdd1-44d2b4c19d6c
cospi(1//2)  # More accurate version to cos(0.5*π)

# ╔═╡ 5844f1f0-bcb9-4b2d-97cd-3e221e4edd17
sqrt(π) ≈ √π

# ╔═╡ fb567d8c-e27f-4ca9-87c4-782f8ab0a8e2
log(10)   # Natual log

# ╔═╡ e6d3cc69-39e4-4da7-86ec-000a8ad24cb9
log10(10) # Common log

# ╔═╡ 0d1863ac-ada5-456c-acce-9ffc5028670b
exp(-5)

# ╔═╡ e9008f32-904f-4b6b-bdf6-b04b15e9bc73
# expm1(x) is more accurate that exp(x) - 1 when x is very small
exp(1e-12)-1, expm1(1e-12)

# ╔═╡ 52f5238a-e018-4612-baa0-3c63165f348e
md"""
# Strings

- [Think Julia ch8](https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#chap08)
- [Introduction_to_Julia_tutorials](https://nbviewer.jupyter.org/github/xorJane/Introduction_to_Julia_tutorials/blob/master/02.Strings.ipynb)

A `string` is a sequence of characters.

- `" ... "` for one line strings.
- Three double quotes to surround multiline string.
- `str1*str2*...` to concatenate strings
- `string(str1, str2, ...)` to convert the data (if neede) and make a string.
- `^` to repeat a string: `str^3` to repeat `str` three times.
- `[idx]` to access an individual character.
- `$` to insert (or interpolate) a value into a string.

Although `string(x, y)` looks less streamlined, it is generally faster than interpolation `$` or concatenation `*` and is most recommended.
"""

# ╔═╡ db06884a-7a79-4972-bd80-7f5bd7410461
"I am a string."

# ╔═╡ 9237e760-8aa5-4b1d-bcca-9f911f6542a8
"""I am also a string."""

# ╔═╡ 8b888ee5-b5ca-48e5-932b-95faf87912b7
'a' == "a"

# ╔═╡ 675dd38b-8ca3-4742-a635-7b74534333d8
str1 = "BEBI"

# ╔═╡ 513552fb-8e07-460c-a2c8-93bbd522beba
str2 = "5009"

# ╔═╡ c0687e61-9c09-42af-a081-4daa7cd5f9a5
string("The class is ", str1, '-', str2, ".")

# ╔═╡ 11ef8f09-a560-4599-b274-2263ea69c854
"The class is $str1-$str2."

# ╔═╡ 0350508e-3c8b-4e21-81c8-61a2d38f8331
str2^5

# ╔═╡ ef124768-6486-4a4b-b06c-74dd966d0569
str1*"-"*str2

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.9"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "81690084b6198a2e1da36fcfda16eeca9f9f24e4"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.1"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "c8abc88faa3f7a3950832ac5d6e690881590d6dc"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "1.1.0"

[[PlutoUI]]
deps = ["Base64", "Dates", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "Suppressor"]
git-tree-sha1 = "44e225d5837e2a2345e69a1d1e01ac2443ff9fcb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.9"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.1.0"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
"""

# ╔═╡ Cell order:
# ╠═fca3485c-552a-4bf2-8cfd-bf7f548a6d1a
# ╠═9f0ba450-ca8a-11eb-2d40-594936524dda
# ╠═8972e6cf-5abd-4e34-9016-d316d3b94fd7
# ╠═183e7365-5aa3-4ebe-8757-ff4ca66da186
# ╠═ff7496bc-5dd8-4b31-a8e2-910f4df257fb
# ╠═fb3654b2-a35e-4db8-8b8f-d8c8e337d14d
# ╠═6bde72a8-e159-4415-85de-c08bc186cff6
# ╠═dedfaa4d-859a-41c0-b92f-de63e02bcf7a
# ╠═22a0a72f-d5e3-440a-b1ea-a8d667f12863
# ╠═d32b9ce8-8aed-46f6-bdf8-22868710cfd7
# ╠═764b4d0a-f643-4f59-9060-bda0be7a6afe
# ╠═3b832f33-a264-4b10-ba4a-14b9f8a3e212
# ╠═bb231574-ce25-4396-a32b-e5ace12937ec
# ╠═a2d2a8c2-6afd-43fb-bb41-7a605f4c6d74
# ╠═2d1eeab7-93ae-412f-a7f1-96f3c023094a
# ╠═7a35ab55-173c-4468-8ca0-0e510d96e0d4
# ╠═73c1c8ce-c227-42be-b5a3-d88745a3e60e
# ╠═48793eaf-cc99-4fac-99d8-d59e2c07b335
# ╠═24a0d4ac-4f72-4127-8afa-87bcbb16c876
# ╠═0bd93bdc-901d-438b-965e-416032f8dd88
# ╠═6539f01d-0b1a-42b9-ac03-882ce9f419f6
# ╠═aaa4115e-4dd9-4f48-9a79-27f209970dce
# ╠═966de6c6-c941-4493-a3eb-51ee7064e538
# ╠═56b16f1b-6809-48f5-85dd-72cf86e7d416
# ╠═b7af6203-833c-4526-bf40-5635a1aac795
# ╠═1b2e41b0-3e5e-4e1b-ae4b-820f27c1cc70
# ╠═3636e2ba-f034-4b73-ba78-48fd65198e7a
# ╠═3929c584-f12b-47f6-913b-4b0ab2c06994
# ╠═5478fc76-c2b6-41ca-8648-b8437caacce4
# ╠═8551ddbe-096d-4c14-bf81-9a2a4e648116
# ╠═bfa38d3f-e655-4b13-be22-88fab5fbdfc7
# ╠═de20682a-1410-40c5-b115-5b5470de870e
# ╠═5a80edaf-bbd0-45b0-af9d-c135749318cc
# ╠═b04a59ff-9def-4799-b47e-1e1e948b2e3f
# ╠═e5e519d3-1d44-4490-a97b-6fb397741a21
# ╠═3d4846ab-b537-416b-bf7e-364215e7bd7e
# ╠═826c913c-a88b-4e4b-8e0a-83888c8df767
# ╠═63dda274-30f3-489d-b7fd-eae5266f097c
# ╠═05f68beb-d349-4d01-aa0d-6e879da8b14e
# ╠═968bad8a-6df8-4b92-9eff-f808e0220bef
# ╠═be6283eb-1265-42f6-9eae-df35edaf9306
# ╠═c5db3a91-0fba-4abd-9fb7-3d7b4b7f53e8
# ╠═6b354bd0-1e53-45f6-a783-6b560019e8d7
# ╠═116c21df-d466-4723-9513-ee83a73fc4bf
# ╠═a63d784c-71e5-4f74-8186-c6887f829900
# ╠═5f949ee8-3694-459f-be29-63e1b4371d30
# ╠═9644b2e4-b270-488b-a6d0-a4ab4c8aa73a
# ╠═f53e2490-f4b0-44dd-9ac9-907e60458faf
# ╠═51bc09d3-a05c-46be-ab94-319ede805291
# ╠═7e9cecbe-8669-44c8-80f0-2b70a3739cc3
# ╠═c026471c-3b69-4089-bdd1-44d2b4c19d6c
# ╠═5844f1f0-bcb9-4b2d-97cd-3e221e4edd17
# ╠═fb567d8c-e27f-4ca9-87c4-782f8ab0a8e2
# ╠═e6d3cc69-39e4-4da7-86ec-000a8ad24cb9
# ╠═0d1863ac-ada5-456c-acce-9ffc5028670b
# ╠═e9008f32-904f-4b6b-bdf6-b04b15e9bc73
# ╠═52f5238a-e018-4612-baa0-3c63165f348e
# ╠═db06884a-7a79-4972-bd80-7f5bd7410461
# ╠═9237e760-8aa5-4b1d-bcca-9f911f6542a8
# ╠═8b888ee5-b5ca-48e5-932b-95faf87912b7
# ╠═675dd38b-8ca3-4742-a635-7b74534333d8
# ╠═513552fb-8e07-460c-a2c8-93bbd522beba
# ╠═c0687e61-9c09-42af-a081-4daa7cd5f9a5
# ╠═11ef8f09-a560-4599-b274-2263ea69c854
# ╠═0350508e-3c8b-4e21-81c8-61a2d38f8331
# ╠═ef124768-6486-4a4b-b06c-74dd966d0569
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

ext: jl


text:
#########################################################################
# Extract shape and other information about arrays that support Python's
# Buffer Interface/Protocol (PEP 3118)
#########################################################################
struct PyArray_Info{T,N}
    native::Bool # native byte order?
    sz::NTuple{N,Int}
    st::NTuple{N,Int} # strides, in multiples of bytes!
    data::Ptr{T}
    readonly::Bool
    pybuf::PyBuffer
end

function PyArray_Info(o::PyObject)
    # n.b. the pydecref(::PyBuffer) finalizer handles releasing the PyBuffer
    pybuf = PyBuffer(o, PyBUF_ND_STRIDED)
    T, native_byteorder = array_format(pybuf)
    sz = size(pybuf)
    strd = strides(pybuf)
    length(strd) == 0 && (sz = ())
    N = length(sz)
    isreadonly = pybuf.buf.readonly==1
    return PyArray_Info{T,N}(native_byteorder, sz, strd, pybuf.buf.buf, isreadonly, pybuf)
end

aligned(i::PyArray_Info{T,N}) where {T,N} = #  FIXME: also check pointer alignment?
  all(m -> m == 0, mod.(i.st, sizeof(T))) # strides divisible by elsize

eltype(i::PyArray_Info{T,N}) where {T,N} = T
ndims(i::PyArray_Info{T,N})  where {T,N} = N

function default_stride(sz::NTuple{N, Int}, ::Type{T}) where {T,N}
    stv = Vector{Int}(N)
    stv[end] = sizeof(T)
    for i = N-1:-1:1
        stv[i] = stv[i+1]*sz[i+1]
    end
    ntuple(i->stv[i], N)
end

"""
`f_contiguous(T::Type, sz::NTuple{N,Int}, st::NTuple{N,Int})::Bool`

Whether an array is in column-major (Fortran, Julia) order, and
has elements stored contiguously. Any array that is f_contiguous will have
identical memory layout to a Julia `Array` of the same size.

`sz` should be the dimensions of the array in number of elements (i.e. what
     `size(a)` would return)
`st` should be the stride(s) *in bytes* between elements in each dimension
"""
function f_contiguous(::Type{T}, sz::NTuple{N,Int}, st::NTuple{N,Int}) where {T,N}
    N == 0 && return true # 0-dimensional arrays have 1 element, always contiguous
    if st[1] != sizeof(T)
        # not contiguous
        return false
    end
    if prod(sz) == 1 || length(sz) == 1
        # 0 or 1-dim arrays (with correct stride) should default to f-contiguous
        # in julia
        return true
    end
    for j = 2:N
        # check stride[cur_dim] == stride[prev_dim]*sz[prev_dim] for all dims>1,
        # implying contiguous column-major storage (n.b. st[1] == sizeof(T) here)
        if st[j] != st[j-1] * sz[j-1]
            return false
        end
    end
    return true
end

f_contiguous(T::Type, sz::NTuple{N1,Int}, st::NTuple{N2,Int}) where {N1,N2} =
    error("stride and size are different lengths, size: $sz, strides: $sz")

f_contiguous(i::PyArray_Info{T,N}) where {T,N} = f_contiguous(T, i.sz, i.st)
c_contiguous(i::PyArray_Info{T,N}) where {T,N} =
    f_contiguous(T, reverse(i.sz), reverse(i.st))

#########################################################################
# PyArray: no-copy wrapper around NumPy ndarray
#
# Hopefully, in the future this can be a subclass of StridedArray (see
# Julia issue #2345), which will allow it to be used with most Julia
# functions, but that is not possible at the moment.  So, to use this
# with Julia linalg functions etcetera a copy is still required.

"""
    PyArray(o::PyObject)

This converts an `ndarray` object `o` to a PyArray.

This implements a nocopy wrapper to a NumPy array (currently of only numeric types only).

If you are using `pycall` and the function returns an `ndarray`, you can use `PyArray` as the return type to directly receive a `PyArray`.
"""
mutable struct PyArray{T,N} <: AbstractArray{T,N}
    o::PyObject
    info::PyArray_Info
    dims::NTuple{N,Int}
    st::NTuple{N,Int}
    f_contig::Bool
    c_contig::Bool
    data::Ptr{T}

    function PyArray{T,N}(o::PyObject, info::PyArray_Info) where {T,N}
        if !aligned(info)
            throw(ArgumentError("only NPY_ARRAY_ALIGNED arrays are supported"))
        elseif !info.native
            throw(ArgumentError("only native byte-order arrays are supported"))
        elseif eltype(info) != T
            throw(ArgumentError("inconsistent type in PyArray constructor"))
        elseif length(info.sz) != N || length(info.st) != N
            throw(ArgumentError("inconsistent ndims in PyArray constructor"))
        end
        return new{T,N}(o, info, tuple(info.sz...), div.(info.st, sizeof(T)),
                        f_contiguous(info), c_contiguous(info),
                        convert(Ptr{T}, info.data))
    end
end

function PyArray(o::PyObject)
    info = PyArray_Info(o)
    return PyArray{eltype(info), length(info.sz)}(o, info)
end

size(a::PyArray) = a.dims
ndims(a::PyArray{T,N}) where {T,N} = N

similar(a::PyArray, ::Type{T}, dims::Dims) where {T} = Array{T}(undef, dims)

"""
Update the data ptr of the `a` to point to the buffer exposed by `o` through
the Python buffer interface
"""
function setdata!(a::PyArray{T,N}, o::PyObject) where {T,N}
    pybufinfo = a.info.pybuf
    PyBuffer!(pybufinfo, o, PyBUF_ND_STRIDED)
    dataptr = pybufinfo.buf.buf
    a.data = reinterpret(Ptr{T}, dataptr)
    a
end

function copy(a::PyArray{T,N}) where {T,N}
    # memcpy is ok iff `a` is f_contig (implies same memory layout as the equiv
    # `Array`) otherwise we do a regular `copyto!`, such that A[I...] == a[I...]
    A = Array{T}(undef, a.dims)
    if a.f_contig
        ccall(:memcpy, Cvoid, (Ptr{T}, Ptr{T}, Int), A, a, sizeof(T)*length(a))
    else
        copyto!(A, a)
    end
    return A
end

unsafe_data_load(a::PyArray, i::Integer) = GC.@preserve a unsafe_load(a.data, i)

@inline data_index(a::PyArray{<:Any,N}, i::CartesianIndex{N}) where {N} =
    1 + sum(ntuple(dim -> (i[dim]-1) * a.st[dim], Val{N}())) # Val lets julia unroll/inline
data_index(a::PyArray{<:Any,0}, i::CartesianIndex{0}) = 1

# handle passing fewer/more indices than dimensions by canonicalizing to M==N
@inline function fixindex(a::PyArray{<:Any,N}, i::CartesianIndex{M}) where {M,N}
    if M == N
        return i
    elseif M < N
        @boundscheck(all(ntuple(k -> size(a,k+M)==1, Val{N-M}())) ||
                     throw(BoundsError(a, i))) # trailing sizes must == 1
        return CartesianIndex(Tuple(i)..., ntuple(k -> 1, Val{N-M}())...)
    else # M > N
        @boundscheck(all(ntuple(k -> i[k+N]==1, Val{M-N}())) ||
                     throw(BoundsError(a, i))) # trailing indices must == 1
        return CartesianIndex(ntuple(k -> i[k], Val{N}()))
    end
end

@inline function getindex(a::PyArray, i::CartesianIndex)
    j = fixindex(a, i)
    @boundscheck checkbounds(a, j)
    unsafe_data_load(a, data_index(a, j))
end
@inline getindex(a::PyArray, i::Integer...) = a[CartesianIndex(i)]
@inline getindex(a::PyArray{<:Any,1}, i::Integer) = a[CartesianIndex(i)]

# linear indexing
function getindex(a::PyArray, i::Integer)
    @boundscheck checkbounds(a, i)
    if a.f_contig
        return unsafe_data_load(a, i)
    else
        @inbounds return a[CartesianIndices(a)[i]]
    end
end

function writeok_assign(a::PyArray, v, i::Integer)
    if a.info.readonly
        throw(ArgumentError("read-only PyArray"))
    else
        GC.@preserve a unsafe_store!(a.data, v, i)
    end
    return v
end

@inline function setindex!(a::PyArray, v, i::CartesianIndex)
    j = fixindex(a, i)
    @boundscheck checkbounds(a, j)
    writeok_assign(a, v, data_index(a, j))
end
@inline setindex!(a::PyArray, v, i::Integer...) = setindex!(a, v, CartesianIndex(i))
@inline setindex!(a::PyArray{<:Any,1}, v, i::Integer) = setindex!(a, v, CartesianIndex(i))

# linear indexing
function setindex!(a::PyArray, v, i::Integer)
    @boundscheck checkbounds(a, i)
    if a.f_contig
        return writeok_assign(a, v, i)
    else
        @inbounds return setindex!(a, v, CartesianIndices(a)[i])
    end
end

stride(a::PyArray, i::Integer) = a.st[i]

Base.unsafe_convert(::Type{Ptr{T}}, a::PyArray{T}) where {T} = a.data

pointer(a::PyArray, i::Int) = pointer(a, ind2sub(a.dims, i))

function pointer(a::PyArray{T}, is::Tuple{Vararg{Int}}) where T
    offset = 0
    for i = 1:length(is)
        offset += (is[i]-1)*a.st[i]
    end
    return a.data + offset*sizeof(T)
end

summary(a::PyArray{T}) where {T} = string(Base.dims2string(size(a)), " ",
                                          string(T), " PyArray")

#########################################################################
# PyArray <-> PyObject conversions

const PYARR_TYPES = Union{Bool,Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64,Float16,Float32,Float64,ComplexF32,ComplexF64,PyPtr,PyObject}

PyObject(a::PyArray) = a.o

convert(::Type{PyArray}, o::PyObject) = PyArray(o)

# PyObject arrays are created by taking a NumPy array of PyPtr and converting
pyo2ptr(T::Type) = T
pyo2ptr(::Type{PyObject}) = PyPtr
pyocopy(a) = copy(a)
pyocopy(a::AbstractArray{PyPtr}) = GC.@preserve a map(pyincref, a)

function convert(::Type{Array{T, 1}}, o::PyObject) where T<:PYARR_TYPES
    try
        return pyocopy(PyArray{pyo2ptr(T), 1}(o, PyArray_Info(o))) # will check T and N vs. info
    catch
        return py2vector(T, o)
    end
end

function convert(::Type{Array{T}}, o::PyObject) where T<:PYARR_TYPES
    try
        info = PyArray_Info(o)
        try
            return pyocopy(PyArray{pyo2ptr(T), length(info.sz)}(o, info)) # will check T == eltype(info)
        catch
            return py2array(T, Array{T}(undef, info.sz...), o, 1, 1)
        end
    catch
        return py2array(T, o)
    end
end

function convert(::Type{Array{T,N}}, o::PyObject) where {T<:PYARR_TYPES,N}
    try
        info = PyArray_Info(o)
        try
            pyocopy(PyArray{pyo2ptr(T),N}(o, info)) # will check T,N == eltype(info),ndims(info)
        catch
            nd = length(info.sz)
            nd == N || throw(ArgumentError("cannot convert $(nd)d array to $(N)d"))
            return py2array(T, Array{T}(undef, info.sz...), o, 1, 1)
        end
    catch
        A = py2array(T, o)
        ndims(A) == N || throw(ArgumentError("cannot convert $(ndims(A))d array to $(N)d"))
        return A
    end
end

array_format(o::PyObject) = array_format(PyBuffer(o, PyBUF_ND_STRIDED))

"""
```
NoCopyArray(o::PyObject)
```
Convert a Python array-like object, to a Julia `Array` or `PermutedDimsArray`
without making a copy of the data.

If the data is stored in row-major format
(the default in Python/NumPy), then the returned array `nca` will be a
`PermutedDimsArray` such that the arrays are indexed the same way in Julia and
Python. i.e. `nca[idxs...] == o[idxs...]`

If the data is stored in column-major format then a regular Julia `Array` will
be returned.

Warning: This function is only lightly tested, and should be considered
experimental - it may cause segmentation faults on conversion or subsequent
array access, or be subtly broken in other ways. Only dense/contiguous, native
endian arrays that support the Python Buffer protocol are likely be converted
correctly.
"""
function NoCopyArray(o::PyObject)
    # n.b. the pydecref(::PyBuffer) finalizer handles releasing the PyBuffer
    pybuf = PyBuffer(o, PyBUF_ND_STRIDED)
    T, native_byteorder = array_format(pybuf)
    !native_byteorder && throw(ArgumentError(
      "Only native endian format supported, format string: '$(get_format_str(pybuf))'"))
    T == Nothing && throw(ArgumentError(
      "Array datatype '$(get_format_str(pybuf))' not supported"))
    # TODO more checks on strides etc
    sz = size(pybuf)
    arr = unsafe_wrap(Array, convert(Ptr{T}, pybuf.buf.buf), sz, own=false)

    !f_contiguous(T, sz, strides(pybuf)) &&
        (arr = PermutedDimsArray(reshape(arr, reverse(sz)), (pybuf.buf.ndim:-1:1)))
    return arr
end


text:
include("line_scatter.jl")

function subplots1()
    p1 = linescatter1()
    p2 = linescatter2()
    p = [p1 p2]
    p
end

function subplots2()
    p1 = linescatter1()
    p2 = linescatter2()
    p = [p1; p2]
    p
end

function subplots3()
    p1 = linescatter6()
    p2 = linescatter2()
    p3 = linescatter3()
    p4 = batman()
    p = [p1 p2; p3 p4]
    p.plot.layout["showlegend"] = false
    p.plot.layout["width"] = 1000
    p.plot.layout["height"] = 600
    p
end

function subplots_withcomprehension()
    hcat([plot(scatter(x=1:5, y=rand(5))) for i in 1:3]...)
end

function subplots_withsharedaxes()
    data =  [
        scatter(x=1:3, y=2:4),
        scatter(x=20:10:40, y=fill(5, 3), xaxis="x2", yaxis="y"),
        scatter(x=2:4, y=600:100:800, xaxis="x", yaxis="y3"),
        scatter(x=4000:1000:6000, y=7000:1000:9000, xaxis="x4", yaxis="y4")
    ]
    layout = Layout(
        xaxis_domain=[0, 0.45],
        yaxis_domain=[0, 0.45],
        xaxis4=attr(domain=[0.55, 1.0], anchor="y4"),
        xaxis2_domain=[0.55, 1],
        yaxis3_domain=[0.55, 1],
        yaxis4=attr(domain=[0.55, 1], anchor="x4")
    )
    plot(data, layout)
end

function with_make_subplots1()

    # The `shared_xaxes` argument to `make_subplots` can be used to link the x
    # axes of subplots in the resulting figure. The `vertical_spacing` argument
    # is used to control the vertical spacing between rows in the subplot grid.

    # Here is an example that creates a figure with 3 vertically stacked
    # subplots with linked x axes. A small vertical spacing value is used to
    # reduce the spacing between subplot rows.

    p = make_subplots(rows=3, cols=1, shared_xaxes=true, vertical_spacing=0.02)
    add_trace!(p, scatter(x=0:2, y=10:12), row=3, col=1)
    add_trace!(p, scatter(x=2:4, y=100:10:120), row=2, col=1)
    add_trace!(p, scatter(x=3:5, y=1000:100:1200), row=1, col=1)
    relayout!(p, title_text="Stacked Subplots with Shared X-Axes")
    p
end

function with_make_subplots2()
    # The `shared_yaxes` argument to `make_subplots` can be used to link the y
    # axes of subplots in the resulting figure.

    # Here is an example that creates a figure with a 2 x 2 subplot grid, where
    # the y axes of each row are linked.

    p = make_subplots(rows=3, cols=2, shared_yaxes=true)
    add_trace!(p, scatter(x=0:2, y=10:12), row=1, col=1)
    add_trace!(p, scatter(x=20:10:40, y=1:3), row=1, col=2)
    add_trace!(p, scatter(x=3:5, y=600:100:800), row=2, col=1)
    add_trace!(p, scatter(x=3:5, y=1000:100:1200), row=2, col=2)
    relayout!(p, title_text="Multiple Subplots with Shared Y-Axes")
    p
end

function with_make_subplots3()
    # The `specs` argument to `make_subplots` is used to configure per-subplot
    # options.  `specs` must be a `Matrix` with dimensions that match those
    # provided as the `rows` and `cols` arguments. The elements of `specs` may
    # either be `missing`, indicating no subplot should be initialized starting
    # with this grid cell, or an instance of `Spec` containing subplot options.
    # The `colspan` subplot option specifies the number of grid columns that the
    # subplot starting in the given cell should occupy.  If unspecified,
    # `colspan` defaults to 1.

    # Here is an example that creates a 2 by 2 subplot grid containing 3
    # subplots. The subplot `specs` element for position (2, 1) has a `colspan`
    # value of 2, causing it to span the full figure width. The subplot `specs`
    # element f or position (2, 2) is `None` because no subplot begins at this
    # location in the grid.
    p = make_subplots(
        rows=2, cols=2,
        specs=[Spec() Spec(); Spec(colspan=2) missing],
        subplot_titles=["First Subplot" "Second Subplot"; "Third Subplot" missing]
    )

    add_trace!(p, scatter(x=[1, 2], y=[1, 2]), row=1, col=1)
    add_trace!(p, scatter(x=[1, 2], y=[1, 2]), row=1, col=2)
    add_trace!(p, scatter(x=[1, 2, 3], y=[2, 1, 2]), row=2, col=1)

    relayout!(p, showlegend=false, title_text="Specs with Subplot Title")
    p
end

function with_make_subplots4()
    # Here is an example that uses the `rowspan` and `colspan` subplot options
    # to create a custom subplot layout with subplots of mixed sizes.
    p = make_subplots(
        rows=5, cols=2,
        specs=[Spec() Spec(rowspan=2)
               Spec() missing
               Spec(rowspan=2, colspan=2) missing
               missing missing
               Spec() Spec()]
    )

    add_trace!(p, scatter(x=[1, 2], y=[1, 2], name="(1,1)"), row=1, col=1)
    add_trace!(p, scatter(x=[1, 2], y=[1, 2], name="(1,2)"), row=1, col=2)
    add_trace!(p, scatter(x=[1, 2], y=[1, 2], name="(2,1)"), row=2, col=1)
    add_trace!(p, scatter(x=[1, 2], y=[1, 2], name="(3,1)"), row=3, col=1)
    add_trace!(p, scatter(x=[1, 2], y=[1, 2], name="(5,1)"), row=5, col=1)
    add_trace!(p, scatter(x=[1, 2], y=[1, 2], name="(5,2)"), row=5, col=2)

    relayout!(p, height=600, width=600, title_text="specs examples")
    p
end


text:
## tests of find_zero interface
using Roots
using Test
using ForwardDiff; Base.adjoint(f) = x -> ForwardDiff.derivative(f, float(x))

# for a user-defined method
import Roots.Setfield
import Roots.Setfield: @set!
struct Order3_Test <: Roots.AbstractSecant end

## Test the interface
@testset "find_zero interface tests" begin

    meths = [Order0(),
             Order1(), Roots.Order1B(), Roots.King(),
             Order2(), Roots.Steffensen(), Roots.Order2B(), Roots.Esser(),
             Order5(), Roots.KumarSinghAkanksha(),
             Order8(), Roots.Thukral8(),
             Order16(), Roots.Thukral16(),
             Roots.LithBoonkkampIJzerman(3,0), Roots.LithBoonkkampIJzerman(4,0)
             ]

    ## different types of initial values
    for m in meths
        @test find_zero(sin, 3, m) ≈ pi
        @test find_zero(sin, 3.0, m) ≈ pi
        @test find_zero(sin, big(3), m) ≈ pi
        @test find_zero(sin, big(3.0), m) ≈ pi
    end

    ## defaults for method argument
    @test find_zero(sin, 3.0) ≈ pi    # order0()
    @test @inferred(find_zero(sin, (3,4))) ≈ π   # Bisection()
    @test @inferred(find_zero(sin, [3,4])) ≈ π   # Bisection()

    ## test tolerance arguments
    ## xatol, xrtol, atol, rtol, maxevals, maxfneval, strict
    fn, xstar = x -> sin(x) - x + 1, 1.9345632107520243
    x0, M = 20.0, Order2()
    @test find_zero(fn, x0, M) ≈ xstar   # needs 16 iterations, 33 fn evaluations, difference is exact

    # test of maxevals
    @test_throws Roots.ConvergenceFailed find_zero(fn, x0, M, maxevals=2)
    # test of maxfneval REMOVED
    # @test_throws Roots.ConvergenceFailed find_zero(fn, x0, M, maxfnevals=2)

    # tolerance on f, atol, rtol: f(x) ~ 0
    M = Order2()
    h = 1e-2
    rt = find_zero(fn, x0, M, atol=h, rtol=0.0)
    @test abs(fn(rt)) > h^2 / 100
    rt = find_zero(fn, x0, M, atol=0.0, rtol=h)
    @test abs(fn(rt)) > h^2 / 100

    ## test of tolerances xatol, xrtol with bisection
    a,b = 1.5, 2.0
    h = 1e-6
    M = Roots.Bisection()
    tracks = Roots.Tracks(Float64, Float64);
    @inferred(find_zero(fn, (a,b), M, tracks=tracks, xatol=h, xrtol=0.0))
    u,v = tracks.xs[end-1:end]
    @test h >= abs(u-v) >= h/2

    ## test of strict
    fn, x0 = x -> cos(x) - 1, pi/4
    @test fn(find_zero(fn, x0, Order5())) <= 1e-8
    @test_throws Roots.ConvergenceFailed find_zero(fn, x0, Order5(), strict=true)

    # xn increment needs atol setting for zeros near 0.0 if strict=true
    M = Order1()
    fn = x -> x * exp(x) + nextfloat(0.0)
    @test_throws Roots.ConvergenceFailed find_zero(fn, 1.0, M, atol=0.0, rtol=0.0, strict=true, xatol=0.0)
    @test abs(find_zero(fn, 1.0, M, atol=0.0, rtol=0.0, strict=true)) <= eps()

    ## test of extreme values for fn, bisection
    c = pi
    fn = x -> Inf * sign(x - c)
    @inferred(find_zero(fn, (-Inf, Inf))) ≈ c

    fn = x -> Inf * x/abs(x) # stop at NaN values
    @inferred(find_zero(fn, (-Inf, Inf))) ≈ 0

    bracketing_meths = (Roots.Bisection(),
                        Roots.A42(), Roots.AlefeldPotraShi(),
                        Roots.Brent(),
                        Roots.FalsePosition(), Roots.FalsePosition(2)
                        )

    # test flexbility in interval specification
    for M ∈ bracketing_meths
        @test @inferred(find_zero(sin, (3,4))) ≈ pi
        @test @inferred(find_zero(sin, [3,4])) ≈ pi
        @test @inferred(find_zero(sin, 3:4)) ≈ pi
        @test @inferred(find_zero(sin, SomeInterval(3,4))) ≈ pi
        @test @inferred(find_zero(sin, range(3,stop=4, length=20))) ≈ pi
    end
end

@testset "non simple zeros" begin

    Ms = (Roots.Order1B(), Roots.Order2B(),Roots.Schroder(),
          Roots.Thukral2B(),Roots.Thukral3B(),Roots.Thukral4B(),Roots.Thukral5B())

    g(x) = exp(x) + x - 2
    f(x) = g(x)^2
    x₀ = 1/4

    α = find_zero(g, x₀)
    fs = (f,f',f'',f''',f'''',f''''',f'''''')
    for M ∈ Ms
        @test find_zero(fs, x₀, M) ≈ α atol=1e-6
    end
end

@testset "find_zero internals" begin

    ##  init_state method
    g1 = x -> x^5 - x - 1
    x0_, xstar_ = (1.0, 2.0), 1.1673039782614187
    M = Roots.A42()
    G1 = Roots.Callable_Function(M, g1)
    state = @inferred(Roots.init_state(M, G1, x0_))
    options = @inferred(Roots.init_options(M, state))
    for M in (Roots.A42(), Roots.Bisection(), Roots.FalsePosition())
        Gₘ = Roots.Callable_Function(M, G1)
        stateₘ = @inferred(Roots.init_state(M, state, Gₘ))
        @test @inferred(solve(M, Gₘ, stateₘ)) ≈ xstar_
    end

    # iterator interface (ZeroProblem, solve; init, solve!)
    meths = [Order1(), Roots.Order1B(), Roots.King(),
             Order2(), Roots.Steffensen(), Roots.Order2B(), Roots.Esser(),
             Order5(), Roots.KumarSinghAkanksha(),
             Order8(), Roots.Thukral8(),
             Order16(), Roots.Thukral16()]
    g1(x) = x^5 - x - 1
    x0_, xstar_ = 1.16, 1.1673039782614187
    fx = ZeroProblem(g1, x0_)
    for M in meths
        @test solve(fx, M) ≈ xstar_
        P = init(fx, M)
        @test solve!(P) ≈ xstar_
    end

    # solve and parameters
    # should be positional, but named supported for now
    g2 = (x,p) -> cos(x) - x/p
    fx = ZeroProblem(g2, (0,pi/2))
    @test solve(fx, 2) ≈ @inferred(find_zero(x -> cos(x) - x/2, (0, pi/2)))
    @test solve(fx, p=2) ≈ @inferred(find_zero(x -> cos(x) - x/2, (0, pi/2)))
    @test @inferred(solve(fx, p=3)) ≈ @inferred(find_zero(x -> cos(x) - x/3, (0, pi/2)))
    g3 = (x, p) -> cos(x) + p[1]*x - p[2]
    fx = ZeroProblem(g3, (0,pi/2))
    @test @inferred(solve(fx, p=[-1/10, 1/10])) ≈ @inferred(find_zero(x -> cos(x) - x/10 - 1/10, (0, pi/2)))

    ## test with early evaluation of bracket
    f = x -> sin(x)
    xs = (3.0, 4.0)
    fxs = f.(xs)
    M = Bisection()
    state = @inferred(Roots.init_state(M, f, xs..., fxs..., m=3.5, fm=f(3.5)))
    @test @inferred(find_zero(M, f, state)) ≈ π

    #     ## hybrid
    g1 = x -> exp(x) - x^4
    x0_, xstar_ = (5.0, 20.0), 8.613169456441398
    M = Roots.Bisection()
    G1 = Roots.Callable_Function(M, g1)
    state = @inferred(Roots.init_state(M, G1, x0_))
    options = @inferred(Roots.init_options(M, state, xatol=1/2))
    ZPI = @inferred(init(M,G1,state,options))
    ϕ = iterate(ZPI)
    while ϕ !== nothing
        val, st = ϕ
        state, ctr = st
        ϕ = iterate(ZPI, st)
    end

    N = Roots.Order1() # switch to N
    G2 = Roots.Callable_Function(N, G1)
    stateₙ = Roots.init_state(N, state, G2)
    options = Roots.init_options(N, stateₙ)
    x = solve(N, G2, stateₙ, options)
    @test x ≈ xstar_

    ## test creation of new methods
    ## xn - f/f' - f'' * f^2 / 2(f')^3 = xn - r1 - r1^2/r2 is third order,
    # had to previousely define:
    function Roots.update_state(M::Order3_Test, f, o::Roots.AbstractUnivariateZeroState{T,S}, options, l=Roots.NullTracks()) where {T, S}
        # xn - f/f' - f'' * f^2 / 2(f')^3 = xn - r1 - r1^2/r2 is third order
        xn_1, xn = o.xn0, o.xn1
        fxn_1, fxn = o.fxn0, o.fxn1

        f_10 = (fxn - fxn_1) / (xn - xn_1)
        xn1::T = xn - fxn / f_10
        fxn1::S = f(xn1)

        f01 = (fxn1 - fxn) /  (xn1 - xn)

        if isnan(f_10) || iszero(f_10) || isnan(f01) || iszero(f01)
            return (o, true)
        end

        r1 = fxn1 / f01
        r2 = f01 / ((f01 - f_10) / (xn1 - xn_1))

        wn = xn1 - r1 - r1^2/r2
        fwn::S = f(wn)

        @set! o.xn0 = xn
        @set! o.xn1 = wn
        @set! o.fxn0 = fxn
        @set! o.fxn1 = fwn

        return (o, false)

    end

    g1 = x -> exp(x) - x^4
    @test find_zero(g1, 8.3, Order3_Test()) ≈ find_zero(g1, 8.3, Order1())

    # test many different calling styles
    f(x) = (sin(x), sin(x)/cos(x)) # x -> (f(x), f(x)/f′(x))
    fs(x) = (sin, cos) # (f, f′)
    x0 = (3,4)
    g(x,p) = begin
        fx = cos(x) - x/p
        (fx, fx/(-sin(x) - 1/p))
    end
    x0a = (0.0, pi/2)
    α₂, α₃ = 1.0298665293222589, 1.1701209500026262
    @test find_zero(f, x0) ≈ π
    @test find_zero(f, first(x0)) ≈ π
    @test find_zero(g, x0a, p=2) ≈ α₂
    @test find_zero(g, first(x0a), p=2) ≈ α₂
    Z = ZeroProblem(f, x0)
    Za = ZeroProblem(g, x0a)
    @test solve(Z) ≈ π
    @test solve(Za, 3) ≈ α₃
    @test solve(Za, p=2) ≈ α₂
    @test solve!(init(Z)) ≈ π
    @test solve!(init(Za, 3)) ≈ α₃
    @test solve!(init(Za, p=3)) ≈ α₃
    Ms = (Roots.Secant(), Roots.Bisection(), Roots.Newton())
    for M ∈ Ms
        @test find_zero(f, x0, M) ≈ π
        @test solve(Z, M) ≈ π
        @test solve!(init(Z, M)) ≈ π
        @test find_zero(g, x0a, M, p=2) ≈ α₂
        @test solve(Za, M, 2) ≈ α₂
        @test solve(Za, M, p=2) ≈ α₂
        @test solve!(init(Za, M, 2)) ≈ α₂
    end

end

@testset "find_zero issue tests" begin

    ## Misc tests
    Ms = [Order0(), Order1(), Roots.Order1B(), Order2(), Roots.Order2B(),
          Order5(), Order8(), Order16()]

    ## issues with starting near a maxima. Some bounce out of it, but
    ## one would expect all to have issues
    fn, xstar = x -> x^3 + 4x^2 -10,  1.365230013414097
    for M in [Order1(), Roots.Order1B(), Order2(), Roots.Order2B(), Order5()]
        @test_throws Roots.ConvergenceFailed find_zero(fn, -1.0, M)
    end
    for M in [Order0(),  Roots.Thukral8(), Roots.Thukral16()]
        @test find_zero(fn, -1.0, M) ≈ xstar
    end

    ## non-simple root
    ## convergence can depend on relaxed convergence checked after an issue
    fn, xstar, x0 = x -> cos(x) - 1, 0.0, 0.1
    for M in Ms
        xn = find_zero(fn, x0, M)
        @test abs(fn(xn)) <= 1e-10
    end
    for M in [Order2(), Order5(), Order8(), Order16()]
        @test_throws Roots.ConvergenceFailed find_zero(fn, x0, M, strict=true)
    end

    ## issue with large steps
    fn, x0 = x -> x^20 - 1, 0.5
    for M in Ms[2:end] # not 0, as it uses bracket
        @test_throws Roots.ConvergenceFailed find_zero(fn, x0, M)
    end

    ## issue with large f''
    fn, x0 = cbrt, 1.0
    for M in [Order1(), Order2(), Order5()]
        @test_throws Roots.ConvergenceFailed find_zero(fn, x0, M)
    end
    ### these stop but only because rtol is used for checking f(xn) ~ 0
    for M in [Roots.Thukral8(), Roots.Thukral16()]
        @test abs(find_zero(fn, x0, M) -  0.0) >= 100
    end

    ## similar (http://people.sc.fsu.edu/~jburkardt/cpp_src/test_zero/test_zero.html)
    function newton_baffler(x)
        a = 1/10
        m, b = 1/4, 1/8

        if x < -a
            m*x - b
        elseif x > a
            m*x + b
        else
            (m*a+b)/a * (x+a) + (-m*a -b)
        end
    end
    for M in (Order0(), Order1(), Roots.Order1B(), Order2(), Roots.Order2B(), Order5(), Order8())
        @test abs(find_zero(newton_baffler, 1.0, M)) <= 1e-15
    end
    for M in (Roots.KumarSinghAkanksha(), Roots.Thukral8(), Roots.Thukral16())
        @test_throws Roots.ConvergenceFailed find_zero(newton_baffler, 1.0, M)
    end

    ## Closed issues ###
    ## issue tests: put in tests to ensure closed issues don't reappear.

    ## issue #94; tolerances not matching documentation
    function test_94(;kwargs...)
        g, T = 1.62850, 14.60000
        α, t1, tf = 0.00347, 40.91375, 131.86573
        y, ya, yf = 0.0, 9000.0, 10000.0
        vy = sqrt(2g*(ya-y))
        θ0, θ1 = atan(α*tf), atan(α*(tf-t1))
        I_sintan(x) = tan(x)/2cos(x) - atanh(tan(x/2))
        I_sintan(x, y) = I_sintan(y) - I_sintan(x)
        function lhs(θ)
            tRem = (vy - T/α*(sec(θ1) - sec(θ))) / g
            val = -yf + y + vy*tRem - 0.5g*tRem^2 - T/α^2*I_sintan(θ, θ1)
            val
        end

        M = Roots.FalsePosition()
        x0 = [atan(α*tf), atan(α*(tf-t1))]
        F = Roots.Callable_Function(M, lhs, nothing) #Roots.DerivativeFree(lhs)
        state = Roots.init_state(M, F, x0)
        options = Roots.init_options(M, state)
        l = Roots.Tracks(state);
        find_zero(M, F, state, options, l)

        @test l.steps <= 45 # 15
    end
    test_94()

    ## Issue with quad_step after truncated M-step PR #140
    @test find_zero(x -> tanh(x) - tan(x), 7.36842 , Order0()) ≈ 7.068582745628732

    ## Use tolerance on f, not x with bisectoin
    atol = 0.01
    u = @inferred(find_zero(sin, (3, 4), atol=atol))
    @test atol >= abs(sin(u)) >= atol^2

    ## issue #159 bracket with zeros should be found
    @test @inferred(find_zero(x->x+1,(-1,1))) == -1

    ## issue #178 passinig through method
    @test fzero(sin, 3, 4, Roots.Brent()) ≈ π

    ## issue #188 with A42
    f = let a=0.18
        x -> x*(1-x^2)/((x^2+a^2)*(1+a^2*x^2))
    end
    r = 0.05
    xs = (r + 1e-12, 1.0)
    @test @inferred(find_zero(x -> f(r)-f(x), xs, Roots.A42())) ≈ 0.4715797678171889
end

struct _SampleCallableObject end
(::_SampleCallableObject)(x) = x^5 - x - 1

mutable struct Cnt
    cnt::Int
    f
    Cnt(f) = new(0, f)
end
(f::Cnt)(x) = (f.cnt += 1; f.f(x))

@testset "find_zero with other callable types" begin

    Ms = [Order0(),
          Order1(), Roots.Order1B(),
          Order2(), Roots.Order2B(),
          Order5(),
          Order8(),
          Order16()]

    for M in Ms
        @test find_zero(_SampleCallableObject(), 1.1, M) ≈ 1.1673039782614187
    end

    for M in Ms
        g = Cnt(x -> x^5 - x - 1)
        @test find_zero(g, 1.1, M) ≈ 1.1673039782614187
        @test g.cnt <= 30
    end

end

@testset "find_bracket test" begin
    fs_xs = ((x -> x^5 - x - 1, (0,2)),
             (sin, (3,4)),
             (x -> exp(x) - x^4, (5, 20))
             )

    for (f, x0) in fs_xs
        out = Roots.find_bracket(f, x0)
        @test prod(f.(out.bracket)) <= 0
    end

    ## test size of  bracket
    for M in  (Roots.BisectionExact(),  Roots.A42(), Roots.AlefeldPotraShi())
        for (fn, b) in ((x  -> x == 0. ? 0. : x/exp(1/(x*x)), (-1.0, 4.0)),
                        (x -> exp(-15*x)*(x-1) + x^15, (0., 1.)),
                        (x -> (-40.0)*x*exp(-x),  (-9, 31))
                        )
            l =  Roots.find_bracket(fn, b, M)
            @test  l.exact ||  abs(l.bracket[2] -  l.bracket[1])  <=  eps(l.xstar)
        end
    end

    ## subnormal
    x0 = nextfloat(nextfloat(0.0))
    fn = x  ->  (x-x0)
    b =  (0.0, 1.0)
    m =  Roots.__middle(b...)  # 1e-154
    M = Roots.BisectionExact()
    l = Roots.find_bracket(fn,  b,  M)
    ## Here l.exact=true, but  this checks the bracket size
    ab = abs(l.bracket[2] - l.bracket[1])
    #@test  !(ab) <=  eps(l.xstar))
    @test  abs(-(l.bracket...)) <=  m

    M =  Roots.A42()
    l = Roots.find_bracket(fn,  b,  M)
    ab = abs(l.bracket[2] - l.bracket[1])
    #@test  !(ab <=  eps(l.xstar))
    @test  ab <= m

    M =  Roots.AlefeldPotraShi()
    l = Roots.find_bracket(fn,  b,  M)
    ab = abs(l.bracket[1] - l.bracket[1])
    #@test  !(ab <=  eps(l.xstar))
    @test  ab ≤ m

    x0 = nextfloat(0.0)
    fn  = x ->  x <=  0.0 ?  x-x0  : x  +  x0
    b = (-1.0, 1.0)
    M = Roots.BisectionExact()
    l = Roots.find_bracket(fn,  b,  M)
    @test !l.exact &&   abs(-(l.bracket...))  <= maximum(eps.(l.bracket))

end

@testset "function evalutions" begin

    function wrapper(f)
        cnt = 0
        x -> begin
            cnt += 1
            f(x)
       end
    end

    # as of v"1.3.0", no more maxfnevals for stopping, just maxevals
    # this is an alternative
    function fz(f, x0::Number, M; maxfnevals=10, kwargs...)
        F = wrapper(f)
        ZPI = init(ZeroProblem(F, x0), M; kwargs...)
        x = NaN * float(x0)
        ϕ = iterate(ZPI)
        while ϕ !== nothing
            x, st = ϕ
            F.cnt.contents >= maxfnevals && return NaN*float(x0)
            ϕ = iterate(ZPI, st)
        end
        x
    end
    f(x) = x^20 - 1
    x0 = 0.9
    M = Order1()
    @test isnan(fz(f, x0, M))  # takes 19 fn evals, not 10

    # test that for update state, fnevals are correctly counted for simpler
    # methods
    fn = (x) -> sin(x); x0 = (3,4); M = Order1()
    state = Roots.init_state(M, Roots.Callable_Function(M,fn), x0)
    options = Roots.init_options(M, state)

    for M ∈ (Order1(), Order2(), Order5(),Order8(), Order16(),
             Roots.Order1B(), Roots.Order2B(),
             Roots.BisectionExact(), Roots.Brent(),
             Roots.A42(), Roots.AlefeldPotraShi())

        # test initial count
        g = wrapper(fn)
        G = Roots.Callable_Function(M,g)
        Roots.init_state(M, G, x0)
        @test  g.cnt.contents ≤ Roots.initial_fncalls(M)

        # test update state
        g = wrapper(fn)
        stateₘ = Roots.init_state(M, state, Roots.Callable_Function(M,f))
        G = Roots.Callable_Function(M,g)
        l = Roots.Tracks(Float64, Float64)
        Roots.update_state(M, G, stateₘ, options,l)
        @test g.cnt.contents == l.fncalls

    end

end

@testset "_extrema" begin
    @test @inferred(Roots._extrema((π, 0))) === (0.0, Float64(π))
    @test @inferred(Roots._extrema([π, 0])) === (0.0, Float64(π))
    @test_throws ArgumentError Roots._extrema(π)
    @test_throws ArgumentError Roots._extrema((π, π))
    @test_throws ArgumentError Roots._extrema([π, π])
end


text:
function symset(v::String, j::Int64)
    return Sym[Sym("$v$i") for i in 1:j]
end

unique_var_count = 0

function uniquevar(v="v")
    global unique_var_count += 1
    return sympify(Symbol("$v$unique_var_count"))
end

function replace(a::Array{T}, d::Dict{T,T}) where T
    b = copy(a)
    replace!(b, d...)
    return b
end

function symfunctions(expr::Sym)
    return Sym.(collect(expr.atoms(AppliedUndef)))
end

# override show for SymPy.SymFunction
# function Base.show(io::IO, f::SymPy.SymFunction)
#     return show(io, Sym(f.x))
# end

function coeff_rem(expr::Sym, t::Sym)
    c = coeff(expr, t)
    return c, expr - c*t
end

function summands(expr::Sym)
    expr = expand(expr)
    if funcname(expr) == "Add"
        return expr.args
    end
    expr
end

function clear_denom(expr::Sym)
    ls = summands(expr)
    ds = denom.(ls)
    val = lcm2(ds...)
    expr *= val
    simplify(expr)
end

function lcm2(n::Sym, rest::Sym...)
    SymPy.lcm(n, lcm2(rest...))
end

function lcm2()
    1
end


text:
function roll_non_degenerate_tet()
    v1 = randn(SVector{3,Float64})
    v2 = randn(SVector{3,Float64})
    v3 = randn(SVector{3,Float64})
    v4 = randn(SVector{3,Float64})
    if volume(v1, v2, v3, v4) < 0.25
        return roll_non_degenerate_tet()
    else
        return v1, v2, v3, v4
    end
end

function make_4_sided()
    v1, v2, v3, v4 = roll_non_degenerate_tet()
    (0.0 < volume(v1, v2, v3, v4)) || error("bad tet")
    A = asMatOnePad(SVector{4,SVector{3,Float64}}(v1, v2, v3, v4))
    inv_A = inv(A)
    n̂ = normalize(randn(SVector{3,Float64}))
    plane = SMatrix{1,4,Float64,4}(n̂[1], n̂[2], n̂[3], randn())
    c = clip_plane_tet(plane, A)
    if length(c) == 4
        return c
    else
        return make_4_sided()
    end
end

dist_from_plane(plane::SMatrix{1,4,Float64,4}, v::SVector{3,Float64}) = dot(v, unPad(plane)) + plane[4]
project_into_plane(plane::SMatrix{1,4,Float64,4}, v::SVector{3,Float64}) = v - dist_from_plane(plane, v) * unPad(plane)


text:
const six_three_one_plus_g_d_split = Dict(
"H" =>
" S   3
 1        18.7311370              0.03349460
 2         2.8253937              0.23472695
 3         0.6401217              0.81375733
 S   1
 1         0.1612778              1.00000000
",

"O" =>
" S   6
 1      5484.6717000              0.0018311
 2       825.2349500              0.0139501
 3       188.0469600              0.0684451
 4        52.9645000              0.2327143
 5        16.8975700              0.4701930
 6         5.7996353              0.3585209
 S   3
 1        15.5396160             -0.1107775
 2         3.5999336             -0.1480263
 3         1.0137618              1.1307670
 S   1
 1         0.2700058              1.0000000
 S   1
 1         0.0845000              1.0000000
 P   3
 1        15.5396160              0.0708743
 2         3.5999336              0.3397528
 3         1.0137618              0.7271586
 P   1
 1         0.2700058              1.0000000
 P   1
 1         0.0845000              1.0000000
 D   1
 1         0.8000000              1.0000000
"
)


text:
using DeepConvert
using Test

include("deepconvert_test.jl")


text:
const COUNTS = Hwloc.histmap(TOPOLOGY);
# TODO: Makes topological assumptions that aren't right for
# multiple nodes or with >3 or <3 levels of Cache.

const CACHE_COUNT = (
    COUNTS[:L1Cache],
    COUNTS[:L2Cache],
    COUNTS[:L3Cache],
    COUNTS[:L4Cache]
)
const NUM_CORES = COUNTS[:Core]

const CACHE_LEVELS = something(findfirst(isequal(0), CACHE_COUNT), length(CACHE_COUNT) + 1) - 1

function define_cache(N)
    if N > CACHE_LEVELS
        return (
            size = nothing,
            depth = nothing,
            linesize = nothing,
            associativity = nothing,
            type = nothing
        )
    end
    cache_name = (:L1Cache, :L2Cache, :L3Cache, :L4Cache)[N]
    c = first(t for t in TOPOLOGY if t.type_ == cache_name && t.attr.depth == N).attr
    (
        size = c.size,
        depth = c.depth,
        linesize = c.linesize,
        associativity = c.associativity,
        type = c.type_
    )
end

const L₁CACHE = define_cache(1)
const L₂CACHE = define_cache(2)
const L₃CACHE = define_cache(3)
const L₄CACHE = define_cache(4)
"""
L₁, L₂, L₃, L₄ cache size
"""
const CACHE_SIZE = (
    L₁CACHE.size,
    L₂CACHE.size,
    L₃CACHE.size,
    L₄CACHE.size
)


text:
#Written By Johnathan Bizzano

include("../../src/HyperSphere.jl")

using .HyperSphere
import .HyperSphere.Functions.Optimizer
import .HyperSphere.Functions.NeuralNet.Layers
import .HyperSphere.Functions.NeuralNet.Layers.Nodes
import .HyperSphere.Functions.NeuralNet.Layers.Functional
import .HyperSphere.Functions.NeuralNet.Activation
import .HyperSphere.Functions.NeuralNet.Initializer
import .HyperSphere.Functions.Error


text:
"""
Produces a publication-quality regression table, similar to Stata's `esttab` and R's `stargazer`.

### Arguments
* `rr::FixedEffectModel...` are the `FixedEffectModel`s from `FixedEffectModels.jl` that should be printed. Only required argument.
* `regressors` is a `Vector` of regressor names (`String`s) that should be shown, in that order. Defaults to an empty vector, in which case all regressors will be shown.
* `fixedeffects` is a `Vector` of FE names (`String`s) that should be shown, in that order. Defaults to an empty vector, in which case all FE's will be shown. Note that the string needs to match the display label exactly, otherwise it will not be shown.
* `labels` is a `Dict` that contains displayed labels for variables (`String`s) and other text in the table. If no label for a variable is found, it default to variable names. See documentation for special values.
* `estimformat` is a `String` that describes the format of the estimate. Defaults to "%0.3f".
* `estim_decoration` is a `Function` that takes the formatted string and the p-value, and applies decorations (such as the beloved stars). Defaults to (* p<0.05, ** p<0.01, *** p<0.001).
* `statisticformat` is a `String` that describes the format of the number below the estimate (se/t). Defaults to "%0.3f".
* `below_statistic` is a `Symbol` that describes a statistic that should be shown below each point estimate. Recognized values are `:blank`, `:se`, and `:tstat`. Defaults to `:se`.
* `below_decoration` is a `Function` that takes the formatted statistic string, and applies a decorations. Defaults to round parentheses.
* `regression_statistics` is a `Vector` of `Symbol`s that describe statistics to be shown at the bottom of the table. Recognized symbols are `:nobs`, `:r2`, `:adjr2`, `:r2_within`, `:f`, `:p`, `:f_kp`, `:p_kp`, and `:dof`. Defaults to `[:nobs, :r2]`.
* `number_regressions` is a `Bool` that governs whether regressions should be numbered. Defaults to `true`.
* `number_regressions_decoration` is a `Function` that governs the decorations to the regression numbers. Defaults to `s -> "(\$s)"`.
* `print_fe_section` is a `Bool` that governs whether a section on fixed effects should be shown. Defaults to `true`.
* `print_estimator_section`  is a `Bool` that governs whether to print a section on which estimator (OLS/IV/NL) is used. Defaults to `true`.
* `standardize_coef` is a `Bool` that governs whether the table should show standardized coefficients. Note that this only works with `TableRegressionModel`s, and that only coefficient estimates and the `below_statistic` are being standardized (i.e. the R^2 etc still pertain to the non-standardized regression).
* `out_buffer` is an `IOBuffer` that the output gets sent to (unless an output file is specified, in which case the output is only sent to the file).
* `renderSettings::RenderSettings` is a `RenderSettings` composite type that governs how the table should be rendered. Standard supported types are ASCII (via `asciiOutput(outfile::String)`) and LaTeX (via `latexOutput(outfile::String)`). If no argument to these two functions are given, the output is sent to STDOUT. Defaults to ASCII with STDOUT.
* `transform_labels` is a `Function` that is used to transform labels. It takes the `String` to be transformed as an argument. See `README.md` for an example.

### Details
A typical use is to pass a number of `FixedEffectModel`s to the function, along with a `RenderSettings` object.
```
regtable(regressionResult1, regressionResult2; renderSettings = asciiOutput())
```
Pass a string to the functions that create a `RenderSettings` to divert output to a file. For example, using LaTeX output,
```
regtable(regressionResult1, regressionResult2; renderSettings = latexOutput("myoutfile.tex"))
```
See the full argument list for details.

### Examples
```julia
using RegressionTables, DataFrames, RDatasets, FixedEffectModels
df = dataset("datasets", "iris")
df[:SpeciesDummy] = pool(df[:Species])
df[:isSmall] = pool(df[:SepalWidth] .< 2.9)
rr1 = reg(df, @model(SepalLength ~ SepalWidth))
rr2 = reg(df, @model(SepalLength ~ SepalWidth + PetalLength   , fe = SpeciesDummy))
rr3 = reg(df, @model(SepalLength ~ SepalWidth + PetalLength + PetalWidth  , fe = SpeciesDummy  + isSmall))
rr4 = reg(df, @model(SepalWidth ~ SepalLength + PetalLength + PetalWidth  , fe = SpeciesDummy))
rr5 = reg(df, @model(SepalWidth ~ SepalLength + (PetalLength ~ PetalWidth)  , fe = SpeciesDummy))
# default
regtable(rr1,rr2,rr3,rr4; renderSettings = asciiOutput())
# display of statistics below estimates
regtable(rr1,rr2,rr3,rr4; renderSettings = asciiOutput(), below_statistic = :blank)
regtable(rr1,rr2,rr3,rr4; renderSettings = asciiOutput(), below_decoration = s -> "[\$(s)]")
# ordering of regressors, leaving out regressors
regtable(rr1,rr2,rr3,rr4; renderSettings = asciiOutput(), regressors = ["SepalLength";"PetalWidth";"SepalWidth"])
# format of the estimates
regtable(rr1,rr2,rr3,rr4; renderSettings = asciiOutput(), estimformat = "%02.5f")
# replace some variable names by other strings
regtable(rr1,rr2,rr3; renderSettings = asciiOutput(), labels = Dict(:SepalLength => "My dependent variable: SepalLength", :PetalLength => "Length of Petal", :PetalWidth => "Width of Petal", Symbol("(Intercept)") => "Const." , :isSmall => "isSmall Dummies", :SpeciesDummy => "Species Dummies"))
# do not print the FE block
regtable(rr1,rr2,rr3,rr4; renderSettings = asciiOutput(), print_fe_section = false)
# re-order fixed effects
regtable(rr1,rr2,rr3,rr4; renderSettings = asciiOutput(), fixedeffects = ["isSmall", "SpeciesDummy"])
# change the yes/no labels in the fixed effect section, and statistics labels
regtable(rr1,rr2,rr3,rr4; renderSettings = asciiOutput(), labels = Dict("__LABEL_FE_YES__" => "Mhm.", "__LABEL_FE_NO__" => "Nope.", "__LABEL_STATISTIC_N__" => "Number of observations", "__LABEL_STATISTIC_R2__" => "R Squared"))
# full set of available statistics
regtable(rr1,rr2,rr3,rr5; renderSettings = asciiOutput(), regression_statistics = [:nobs, :r2, :adjr2, :r2_within, :f, :p, :f_kp, :p_kp, :dof])
# LaTeX output
regtable(rr1,rr2,rr3,rr4; renderSettings = latexOutput())
# LaTeX output to file
regtable(rr1,rr2,rr3,rr4; renderSettings = latexOutput("myoutfile.tex"))
```
"""
function regtable(rr::Union{FixedEffectModel,TableRegressionModel}...;
    regressors::Vector{String} = Vector{String}(),
    fixedeffects::Vector{String} = Vector{String}(),
    labels::Dict{String,String} = Dict{String,String}(),
    estimformat::String = "%0.3f",
    estim_decoration::Function = default_ascii_estim_decoration,
    statisticformat::String = "%0.3f",
    below_statistic::Symbol = :se,
    below_decoration::Function = s::String -> "($s)",
    regression_statistics::Vector{Symbol} = [:nobs, :r2],
    number_regressions::Bool = true,
    number_regressions_decoration::Function = i::Int64 -> "($i)",
    print_fe_section = true,
    print_estimator_section = true,
    standardize_coef = false,
    out_buffer = IOBuffer(),
    transform_labels::Function = identity,
    renderSettings::RenderSettings = asciiOutput()
    )

    # define some functions that makes use of StatsModels' RegressionModels
    coefnames(r::TableRegressionModel) = StatsModels.coefnames(r.mf)
    coefnames(r::FixedEffectModel) = String.(r.coefnames) # this will need to be updated when we move
                                                                  # to FixedEffectModels 0.8.2
    # if standardize_coef == true
    #     function coef(r::TableRegressionModel)
    #         cc = StatsModels.coef(r)
    #         return [ cc[i]*std(r.model.pp.X[:,i])/std(r.model.rr.y) for i in 1:length(cc) ]
    #     end
    #     function vcov(r::TableRegressionModel)
    #         vc = StatsModels.vcov(r)
    #         mul = [ std(r.model.pp.X[:,i])*std(r.model.pp.X[:,j])/(std(r.model.rr.y)*std(r.model.rr.y)) for i in 1:size(cc,1), j in 1:size(cc,1)  ]
    #         return mul .* vc
    #     end
    # else # do not standardize
    #     function coef(r::TableRegressionModel)
    #         return StatsModels.coef(r)
    #     end
    #     function vcov(r::TableRegressionModel)
    #         return StatsModels.vcov(r)
    #     end
    # end
    coef(r::FixedEffectModel) = r.coef
    vcov(r::FixedEffectModel) = r.vcov
    coef(r::TableRegressionModel) = StatsModels.coef(r)
    vcov(r::TableRegressionModel) = StatsModels.vcov(r)
    df_residual(r::FixedEffectModel) = dof_residual(r)
    df_residual(r::TableRegressionModel) = dof_residual(r)
    yname(r::FixedEffectModel) = r.yname
    yname(r::TableRegressionModel) = r.mf.f.lhs.sym # returns a Symbol
    ther2(r::FixedEffectModel) = r.r2
    ther2(r::TableRegressionModel) = isa(r.model, LinearModel) ? r2(r) : NaN

    # print a warning message if standardize_coef == true but one
    # of the regression results is not a TableRegressionModel
    if standardize_coef && any(.!isa.(rr,StatsModels.TableRegressionModel))
        @warn("Standardized coefficients are only shown for TableRegressionModel regression results.")
    end

    numberOfResults = size(rr,1)

    # Create an RegressionTable from the regression results

    # ordering of regressors:
    if length(regressors) == 0
        # construct default ordering: from ordering in regressions (like in Stata)
        regressorList = Vector{String}()
        for r in rr # FixedEffectModel
            names = coefnames(r)
            for regressorIndex = 1:length(names)
                if !(any(regressorList .== names[regressorIndex]))
                    # add to list
                    push!(regressorList, names[regressorIndex])
                end
            end
        end
    else
        # take the list of regressors from the argument
        regressorList = regressors
    end

    # for each regressor, check each regression result, calculate statistic, and construct block
    estimateBlock = Array{String}(undef,0,numberOfResults+1)
    for regressor in regressorList
        estimateLine = fill("", 2, numberOfResults+1)
        for resultIndex = 1:numberOfResults
            thiscnames = coefnames(rr[resultIndex])
            thiscoef = coef(rr[resultIndex])
            thisvcov = vcov(rr[resultIndex])
            if standardize_coef && isa(rr[resultIndex],StatsModels.TableRegressionModel)
                thiscoef = [ thiscoef[i]*std(rr[resultIndex].model.pp.X[:,i])/std(rr[resultIndex].model.rr.y) for i in 1:length(thiscoef) ]
                mul = [ std(rr[resultIndex].model.pp.X[:,i])*std(rr[resultIndex].model.pp.X[:,j])/(std(rr[resultIndex].model.rr.y)*std(rr[resultIndex].model.rr.y)) for i in 1:length(thiscoef), j in 1:length(thiscoef)  ]
                thisvcov  = mul .* thisvcov
            end
            thisdf_residual = df_residual(rr[resultIndex])
            index = findall(regressor .== thiscnames)
            if !isempty(index)
                pval = ccdf(FDist(1, thisdf_residual ), abs2(thiscoef[index[1]]/sqrt(thisvcov[index[1],index[1]])))
                estimateLine[1,resultIndex+1] = estim_decoration(sprintf1(estimformat,thiscoef[index[1]]),pval)
                if below_statistic == :tstat
                    s = sprintf1(statisticformat, thiscoef[index[1]]/sqrt(thisvcov[index[1],index[1]]))
                    estimateLine[2,resultIndex+1] = below_decoration(s)
                elseif below_statistic == :se
                    s = sprintf1(statisticformat, sqrt(thisvcov[index[1],index[1]]))
                    estimateLine[2,resultIndex+1] = below_decoration(s)
                elseif below_statistic == :blank
                    estimateLine[2,resultIndex+1] = "" # for the sake of completeness
                end
            end
        end
        # check if the regressor was not found
        if estimateLine == fill("", 2, numberOfResults+1)
           @warn("Regressor $(String(regressor)) not found among regression results.")
        else
            # add label on the left:
            estimateLine[1,1] = haskey(labels,regressor) ? labels[regressor] : transform_labels(regressor)
            # add to estimateBlock
            estimateBlock = [estimateBlock; estimateLine]
        end
    end

    # Regressand block
    #   needs to be separately rendered
    regressandBlock = fill("", 1, numberOfResults+1)
    for rIndex = 1:numberOfResults
        # keep in mind that yname is a Symbol
        regressandBlock[1,rIndex+1] = haskey(labels,string(yname(rr[rIndex]))) ? labels[string(yname(rr[rIndex]))] : transform_labels(string(yname(rr[rIndex])))
    end

    # Regression numbering block (if we do it)
    if number_regressions
        regressionNumberBlock = fill("", 1, numberOfResults + 1)
        for rIndex = 1:numberOfResults
            regressionNumberBlock[1,rIndex+1] = number_regressions_decoration(rIndex)
        end
    end

    # Fixed effects block
    print_fe_block = print_fe_section && any(isFERegressionResult.(rr))
    if print_fe_block

        # construct list of fixed effects for display
        feList = Vector{AbstractTerm}()
        for r in rr if isFERegressionResult(r)

            for term in eachterm(r.formula.rhs)
                if has_fe(term) && !(any(name.(feList) .== name(term)))
                    push!(feList, term)
                end
            end

            # if isa(r.feformula, Symbol)
            #     if !(any(feList .== string(r.feformula)))
            #         # add to list
            #         push!(feList, string(r.feformula))
            #     end
            # elseif r.feformula.args[1] == :+
            #     x = r.feformula.args
            #     for i in 2:length(x)
            #         if isa(x[i], Symbol) | isa(x[i], Expr) # if expression, push the whole expression
            #             if !(any(feList .== string(x[i])))
            #                 # add to list
            #                 push!(feList, string(x[i]))
            #             end
            #         end
            #     end
            # elseif r.feformula.args[1] == :*
            #     # push the whole interaction
            #     if !(any(feList .== string(r.feformula)))
            #         push!(feList, string(r.feformula))
            #     end
            # end

        end end
        # in case the user supplies a list of FE's, cut down the list
        if !isempty(fixedeffects)
            # compare the user-supplied list to `feList`, and keep only those that match
            for i = length(feList):-1:1
                if !any(name(feList[i]) .== fixedeffects)
                    deleteat!(feList, [i])
                end
            end
        end

        # construct a list of fixed effects (Term's) for each RegressionResult
        febyrr = Vector{Vector{AbstractTerm}}()
        for r in rr
            fe = Vector{AbstractTerm}()
            if isFERegressionResult(r)
                for term in eachterm(r.formula.rhs)
                    if has_fe(term)
                        push!(fe, term)
                    end
                end
                # if isa(r.feformula, Symbol)
                #     # add to list
                #     push!(fe, string(r.feformula))
                # elseif r.feformula.args[1] == :+
                #     x = r.feformula.args
                #     for i in 2:length(x)
                #         if isa(x[i], Symbol) | isa(x[i], Expr) # if expression, push the whole expression
                #             # add to list
                #             push!(fe, string(x[i]))
                #         end
                #     end
                # elseif r.feformula.args[1] == :*
                #     push!(fe, string(r.feformula))
                # end
            end
            push!(febyrr, fe)
        end

        # construct FE block
        feBlock = Array{String}(undef,0,numberOfResults+1)
        for fe in feList
            feLine = fill("", 1, numberOfResults+1)
            for resultIndex = 1:numberOfResults if isFERegressionResult(rr[resultIndex])
                index = findall(name(fe) .== name.(febyrr[resultIndex]))
                if !isempty(index)
                    feLine[1,resultIndex+1] = haskey(labels, "__LABEL_FE_YES__") ? labels["__LABEL_FE_YES__"] : renderSettings.label_fe_yes
                else
                    feLine[1,resultIndex+1] = haskey(labels, "__LABEL_FE_NO__") ? labels["__LABEL_FE_NO__"] : renderSettings.label_fe_no
                end
            end end
            # check if the regressor was not found
            if feLine == fill("", 1, numberOfResults+1)
               @warn("Fixed effect $fe not found in any regression results.")
            else
                # add label on the left:
                feLine[1,1] = haskey(labels,name(fe)) ? labels[name(fe)] : transform_labels(name(fe))
                # add to estimateBlock
                feBlock = [feBlock; feLine]
            end
        end
    end

    if print_estimator_section
        estimatorBlock = fill("", 1, numberOfResults+1)
        estimatorBlock[1,1] = haskey(labels, "__LABEL_ESTIMATOR__") ? labels["__LABEL_ESTIMATOR__"] : renderSettings.label_estimator
        for i = 1:numberOfResults
            if isOLSRegressionResult(rr[i])
                estimatorBlock[1,i+1] =  haskey(labels, "__LABEL_ESTIMATOR_OLS__") ? labels["__LABEL_ESTIMATOR_OLS__"] : renderSettings.label_estimator_ols
            elseif isIVRegressionResult(rr[i])
                estimatorBlock[1,i+1] =  haskey(labels, "__LABEL_ESTIMATOR_IV__") ? labels["__LABEL_ESTIMATOR_IV__"] : renderSettings.label_estimator_iv
            else
                estimatorBlock[1,i+1] =  haskey(labels, "__LABEL_ESTIMATOR_NL__") ? labels["__LABEL_ESTIMATOR_NL__"] : renderSettings.label_estimator_nl
            end
        end
    end

    if length(regression_statistics)>0
        # we have a statistics block (N, R^2, etc)
        print_statistics_block = true

        # one line for each statistic
        statisticBlock = fill("", length(regression_statistics), numberOfResults+1)
        for i = 1:length(regression_statistics)
            if regression_statistics[i] == :r2_a
                @warn "Use :adjr2 instead of :r2_a"
                regression_statistics[i] = :adjr2
            end
            if regression_statistics[i] == :nobs
                statisticBlock[i,1] = haskey(labels, "__LABEL_STATISTIC_N__") ? labels["__LABEL_STATISTIC_N__"] : renderSettings.label_statistic_n
                for resultIndex = 1:numberOfResults
                    statisticBlock[i,resultIndex+1] = sprintf1("%'i",nobs(rr[resultIndex]))
                end
            elseif regression_statistics[i] == :r2
                statisticBlock[i,1] = haskey(labels, "__LABEL_STATISTIC_R2__") ? labels["__LABEL_STATISTIC_R2__"] : renderSettings.label_statistic_r2
                for resultIndex = 1:numberOfResults
                    statisticBlock[i,resultIndex+1] = isnan(ther2(rr[resultIndex])) ? "" : sprintf1(statisticformat, ther2(rr[resultIndex]))
                end
            elseif regression_statistics[i] == :adjr2
                statisticBlock[i,1] = haskey(labels, "__LABEL_STATISTIC_adjr2__") ? labels["__LABEL_STATISTIC_adjr2__"] : renderSettings.label_statistic_adjr2
                for resultIndex = 1:numberOfResults
                    statisticBlock[i,resultIndex+1] = isdefined(rr[resultIndex], :adjr2) && !isnothing(rr[resultIndex].adjr2) ? sprintf1(statisticformat, rr[resultIndex].adjr2) : ""
                end
            elseif regression_statistics[i] == :r2_within
                statisticBlock[i,1] = haskey(labels, "__LABEL_STATISTIC_R2_WITHIN__") ? labels["__LABEL_STATISTIC_R2_WITHIN__"] : renderSettings.label_statistic_r2_within
                for resultIndex = 1:numberOfResults
                    statisticBlock[i,resultIndex+1] = isdefined(rr[resultIndex], :r2_within) && !isnothing(rr[resultIndex].r2_within) ? sprintf1(statisticformat, rr[resultIndex].r2_within) : ""
                end
            elseif regression_statistics[i] == :f
                statisticBlock[i,1] = haskey(labels, "__LABEL_STATISTIC_F__") ? labels["__LABEL_STATISTIC_F__"] : renderSettings.label_statistic_f
                for resultIndex = 1:numberOfResults
                    statisticBlock[i,resultIndex+1] = isdefined(rr[resultIndex], :F) ? sprintf1(statisticformat, rr[resultIndex].F) : ""
                end
            elseif regression_statistics[i] == :p
                statisticBlock[i,1] = haskey(labels, "__LABEL_STATISTIC_P__") ? labels["__LABEL_STATISTIC_P__"] : renderSettings.label_statistic_p
                for resultIndex = 1:numberOfResults
                    statisticBlock[i,resultIndex+1] = isdefined(rr[resultIndex], :p) ? sprintf1(statisticformat, rr[resultIndex].p) : ""
                end
            elseif regression_statistics[i] == :f_kp
                statisticBlock[i,1] = haskey(labels, "__LABEL_STATISTIC_F_KP__") ? labels["__LABEL_STATISTIC_F_KP__"] : renderSettings.label_statistic_f_kp
                for resultIndex = 1:numberOfResults
                    statisticBlock[i,resultIndex+1] = isdefined(rr[resultIndex], :F_kp) && !isnothing(rr[resultIndex].F_kp) ? sprintf1(statisticformat, rr[resultIndex].F_kp) : ""
                end
            elseif regression_statistics[i] == :p_kp
                statisticBlock[i,1] = haskey(labels, "__LABEL_STATISTIC_P_KP__") ? labels["__LABEL_STATISTIC_P_KP__"] : renderSettings.label_statistic_p_kp
                for resultIndex = 1:numberOfResults
                    statisticBlock[i,resultIndex+1] = isdefined(rr[resultIndex], :p_kp) && !isnothing(rr[resultIndex].p_kp) ? sprintf1(statisticformat, rr[resultIndex].p_kp) : ""
                end
            elseif regression_statistics[i] == :dof
                statisticBlock[i,1] = haskey(labels, "__LABEL_STATISTIC_DOF__") ? labels["__LABEL_STATISTIC_DOF__"] : renderSettings.label_statistic_dof
                for resultIndex = 1:numberOfResults
                    statisticBlock[i,resultIndex+1] = isdefined(rr[resultIndex], :dof_residual) ? sprintf1("%i",rr[resultIndex].dof_residual) : ""
                end
            end

        end
    else
        print_statistics_block = false
    end

    # construct alignment string:
    align = "l" * ("r" ^ numberOfResults)

    bodyBlocks = [estimateBlock]

    if print_fe_block
        push!(bodyBlocks, feBlock)
    end

    if print_estimator_section
        push!(bodyBlocks, estimatorBlock)
    end

    if print_statistics_block
        push!(bodyBlocks, statisticBlock)
    end

    # if we're numbering the regression columns, add a block before the other stuff

    if number_regressions
        insert!(bodyBlocks,1,regressionNumberBlock)
    end

    # create RegressionTable
    tab = RegressionTable(numberOfResults+1, "", regressandBlock, bodyBlocks , "")

    # create output stream
    if renderSettings.outfile == ""
        outstream = out_buffer
    else
        try
            outstream = open(renderSettings.outfile, "w")
        catch ex
            error("Error opening file $(renderSettings.outfile): $(ex)")
        end
    end

    render(outstream, tab , align, renderSettings )

    # if we're writing to a file, close it
    if renderSettings.outfile != ""
        close(outstream)
    else # else print the table
        println(Compat.String(take!(copy(outstream))))
    end

end


text:
export isadmtup

# Is sequence an admissible prime k-tuple.
function isadmtup(x)
    mx = sort( @compat rem.(x,length(x)))
    @inbounds for i in 0:length(x)-1
        mx[i+1] == i || return true
    end
    false
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Constants to compute the Equation of Origins according to the IAU-2006
#   theory.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Obtained from the files:
#
#   ftp://tai.bipm.org/iers/conv2010/chapter5/tab5.2e.txt
#
# accessed on 18 Mar 2021.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

const COEFS_IAU2006_EO_0 = Float64[
    1        2640.96         -0.39    0    0    0    0    1    0    0    0    0    0    0    0    0    0
    2          63.52         -0.02    0    0    0    0    2    0    0    0    0    0    0    0    0    0
    3          11.75          0.01    0    0    2   -2    3    0    0    0    0    0    0    0    0    0
    4          11.21          0.01    0    0    2   -2    1    0    0    0    0    0    0    0    0    0
    5          -4.55          0.00    0    0    2   -2    2    0    0    0    0    0    0    0    0    0
    6           2.02          0.00    0    0    2    0    3    0    0    0    0    0    0    0    0    0
    7           1.98          0.00    0    0    2    0    1    0    0    0    0    0    0    0    0    0
    8          -1.72          0.00    0    0    0    0    3    0    0    0    0    0    0    0    0    0
    9          -1.41         -0.01    0    1    0    0    1    0    0    0    0    0    0    0    0    0
   10          -1.26         -0.01    0    1    0    0   -1    0    0    0    0    0    0    0    0    0
   11          -0.63          0.00    1    0    0    0   -1    0    0    0    0    0    0    0    0    0
   12          -0.63          0.00    1    0    0    0    1    0    0    0    0    0    0    0    0    0
   13           0.46          0.00    0    1    2   -2    3    0    0    0    0    0    0    0    0    0
   14           0.45          0.00    0    1    2   -2    1    0    0    0    0    0    0    0    0    0
   15           0.36          0.00    0    0    4   -4    4    0    0    0    0    0    0    0	  0    0
   16          -0.24         -0.12    0    0    1   -1    1    0   -8   12    0    0    0    0    0    0
   17           0.32          0.00    0    0    2    0    0    0    0    0    0    0    0    0    0    0
   18           0.28          0.00    0    0    2    0    2    0    0    0    0    0    0    0    0    0
   19           0.27          0.00    1    0    2    0    3    0    0    0    0    0    0    0    0    0
   20           0.26          0.00    1    0    2    0    1    0    0    0    0    0    0    0    0    0
   21          -0.21          0.00    0    0    2   -2    0    0    0    0    0    0    0    0    0    0
   22           0.19          0.00    0    1   -2    2   -3    0    0    0    0    0    0    0    0    0
   23           0.18          0.00    0    1   -2    2   -1    0    0    0    0    0    0    0    0    0
   24          -0.10          0.05    0    0    0    0    0    0    8  -13    0    0    0    0    0   -1
   25           0.15          0.00    0    0    0    2    0    0    0    0    0    0    0    0    0    0
   26          -0.14          0.00    2    0   -2    0   -1    0    0    0    0    0    0    0    0    0
   27           0.14          0.00    1    0    0   -2    1    0    0    0    0    0    0    0    0    0
   28          -0.14          0.00    0    1    2   -2    2    0    0    0    0    0    0    0    0    0
   29           0.14          0.00    1    0    0   -2   -1    0    0    0    0    0    0    0    0    0
   30           0.13          0.00    0    0    4   -2    4    0    0    0    0    0    0    0    0    0
   31          -0.11          0.00    0    0    2   -2    4    0    0    0    0    0    0    0    0    0
   32           0.11          0.00    1    0   -2    0   -3    0    0    0    0    0    0    0    0    0
   33           0.11          0.00    1    0   -2    0   -1    0    0    0    0    0    0    0    0    0
]'

const COEFS_IAU2006_EO_1 = Float64[
  34          -0.87          0.00    0    0    0    0    1    0    0    0    0    0    0    0    0    0
]'


text:
# Note: implement as generic, in order to also support views which are of type SubArray
@inline function swap!(a::AbstractVector, i, j)
    tmp = a[i]
    a[i] = a[j]
    a[j] = tmp
    a
end


text:
using Test
using LinearAlgebra

@testset "tt_up_rks" begin
	x = randn(4,4,4,4)
	x_tt = ttv_decomp(x,tol=0.1)
	x = ttv_to_tensor(x_tt)
	y_tt = tt_up_rks(x_tt,20)
	@test isapprox(x,ttv_to_tensor(y_tt),atol=1e-10)
	z_tt = tt_up_rks(x_tt,20,ϵ_wn=1e-10)
	@test isapprox(x,ttv_to_tensor(z_tt),atol=1e-6)
end

@testset "TT operations" begin
    n = 3
    d = 3
    L = randn(n,n,n,n,n,n)
    S = randn(n,n,n,n,n,n)
    L_tto = tto_decomp(L)
    S_tto = tto_decomp(S)
    x = randn(n,n,n)
    y = randn(n,n,n)
    x_tt = ttv_decomp(x)
    y_tt = ttv_decomp(y)
    a = randn()
    #mult_a_tt
    @test isapprox(a.*x, ttv_to_tensor(a*x_tt))
    #tt_add
    z_tt = x_tt+y_tt
    @test(isapprox(ttv_to_tensor(z_tt),x+y))
    #tt_dot
    @test isapprox(x[:]'*(x+y)[:], dot(x_tt,z_tt))
    #tto_add
    A_tto = L_tto+S_tto
    @test(isapprox(tto_to_tensor(A_tto),L+S)) #tto_add test
    #tto*ttv multiplication
    z = reshape(L,n^d,:)*x[:]
    z_tt = L_tto*x_tt
    @test(isapprox(ttv_to_tensor(z_tt)[:],z))
    #tto*tto multiplication
    A = reshape(L,n^d,:)*reshape(S,n^d,:)
    A_tto = L_tto*S_tto
    @test(isapprox(reshape(tto_to_tensor(A_tto),n^d,:),A))
end

@testset "TT rounding" begin
    n=5
    d=3
    L = randn(n,n,n,n,n,n)
    x = randn(n,n,n)
    y = reshape(L,n^d,:)*x[:]
    L_tt = tto_decomp(L)
    x_tt = ttv_decomp(x)
    y_tt = L_tt*x_tt
    @test(isapprox(ttv_to_tensor(tt_rounding(y_tt))[:],y))
    @test(isapprox(dot(x_tt,y_tt),dot(x[:],y)))
#    @test(isapprox(ttv_to_tensor(tt_compression_par(y_tt))[:],y))
end

@testset "TT orthogonalize" begin
    dims = 2*ones(Int,8)
    rks = [1,2,4,8,8,8,4,2,1]
    x_tt = rand_tt(dims,rks)
    y_tt = orthogonalize(x_tt)
    i = rand(2:8)
    @test(isapprox(sum(y_tt.ttv_vec[i][k,:,:]*y_tt.ttv_vec[i][k,:,:]' for k in 1:dims[i]),Matrix(I,rks[i],rks[i])))
    @test(size(y_tt.ttv_vec[i])==(dims[i],rks[i],rks[i+1]))
    @test(isapprox(ttv_to_tensor(x_tt),ttv_to_tensor(y_tt)))
end

max_stars_repo_name: msdupuy/Tensor-Train-Julia

max_issues_repo_name: msdupuy/Tensor-Train-Julia

max_forks_repo_name: msdupuy/Tensor-Train-Julia


text:
function SeisProcess(in::AbstractString,out::AbstractString,operators,parameters;key=[])
    # Run processing flows that read and write from disk
    #
    # f is a function that has the following syntax: d2,h2 = f(d1,h1,param),
    # where
    # param is list of keyword arguments for the function.
    # note that f can be a vector of functions. They will be executed
    # sequentially on the same group of traces.
    #
    # get list of gather lengths
    ext = ReadTextHeader(in)
    nx = ext.n2*ext.n3*ext.n4*ext.n5
    filename_headers = ParseHeaderName(in)
    stream_h = open(filename_headers)
    curr = zeros(length(key),1)
    itrace = 1
    h1 = GrabHeader(stream_h,1)
    for ikey = 1 : length(key)
	curr[ikey] = getfield(h1,Symbol(key[ikey]))
    end
    prev = 1*curr
    L = []
    for j = 1 : nx
	h1 = GrabHeader(stream_h,j)
	for ikey = 1 : length(key)
	    curr[ikey] = getfield(h1,Symbol(key[ikey]))
	end
	if curr != prev
	    push!(L,j - itrace)
	    itrace = j
	end
	prev = 1*curr
    end
    push!(L,nx - itrace + 1)
    close(stream_h)

    itrace_in = 1
    itrace_out = 1
    for igather = 1 : length(L)
	d1,h1,e1 = SeisRead(in, group="some",itrace=itrace_in,ntrace=L[igather])
	num_traces_in = size(d1,2)
	for j = 1 : length(operators)
	    op = operators[j]
	    d2,h2 = op(d1,h1,e1;parameters[j]...)
	    d1 = copy(d2)
	    h1 = copy(h2)
	end
	num_traces_out = size(d1,2)
	SeisWrite(out,d1,h1,e1,itrace=itrace_out)
	itrace_in += num_traces_in
	itrace_out += num_traces_out
    end

end

function SeisProcess(in::Array{AbstractString,1}, out::Array{AbstractString,1}, operators,
                     parameters;key=[])

    for j = 1 : length(in)
	SeisProcess(in[j],out[j],operators,parameters;key=key)
    end

end

function SeisProcess(in1::AbstractString, in2::AbstractString, out::AbstractString, operators,
                     parameters; key=[])
    # Run processing flows that read 2 inputs and write 1 output
    #
    # f is a function that has the following syntax:
    # d3,h3 = f(d1,d2,h1,h2,param), where
    # param is list of keyword arguments for the function.
    # note that f can be a vector of functions.
    # They will be executed sequentially on the same group of traces.

    # get list of gather lengths
    ext = ReadTextHeader(in1)
    nx = ext.n2*ext.n3*ext.n4*ext.n5
    filename_headers = ParseHeaderName(in1)
    stream_h = open(filename_headers)
    curr = zeros(length(key),1)
    itrace = 1
    h1 = GrabHeader(stream_h,1)
    for ikey = 1 : length(key)
	curr[ikey] = getfield(h1,Symbol(key[ikey]))
    end
    prev = 1*curr
    L = []
    for j = 1 : nx
	h1 = GrabHeader(stream_h,j)
	for ikey = 1 : length(key)
	    curr[ikey] = getfield(h1,Symbol(key[ikey]))
	end
	if curr != prev
	    push!(L,j - itrace)
	    itrace = j
	end
	prev = 1*curr
    end
    push!(L,nx - itrace + 1)
    close(stream_h)

    itrace_in = 1
    itrace_out = 1
    for igather = 1 : length(L)
	d1,h1,e1 = SeisRead(in1,group="some",itrace=itrace_in,ntrace=L[igather])
	d2,h2,e2 = SeisRead(in2,group="some",itrace=itrace_in,ntrace=L[igather])
	num_traces_in = size(d1,2)
	for j = 1 : length(operators)
	    op = operators[j]
	    d3,h3 = op(d1,d2,h1,h2;parameters[j]...)
	    d1 = copy(d3)
	    h1 = copy(h3)
	end
	num_traces_out = size(d1,2)
	SeisWrite(out,d1,h1,e1,itrace=itrace_out)
	itrace_in += num_traces_in
	itrace_out += num_traces_out
    end

end


text:
# Autogenerated wrapper script for Check_jll for x86_64-apple-darwin
export libcheck

JLLWrappers.@generate_wrapper_header("Check")
JLLWrappers.@declare_library_product(libcheck, "@rpath/libcheck.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libcheck,
        "lib/libcheck.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
__precompile__(true)

module NumericalTricks

using StaticArrays
using LinearAlgebra
using Rotations
using ForwardDiff: Dual, value

include("geometry_kernel.jl")
include("poly_approx.jl")
include("div_by_zero.jl")
include("rotations.jl")
include("vector_projections.jl")
include("utility.jl")
include("basic_dh.jl")
# include("matrix_factor_derivatives.jl")

export
  # geometry_kernel.jl
  area,
  centroid,
  triangleNormal,
  triangle_area,
  vector_area,
  volume,

  # # poly_approx.jl
  # fastSoftPlus,
  # fastSigmoid,
  # soft_clamp,
  # # smooth_c1_ramp,
  # # smooth_c2_ramp,

  # # div_by_zero.jl
  # safe_normalize,
  # unsafe_normalize,
  # safe_inv_norm²,
  # safe_inv_norm,
  # safe_norm,
  # norm²,
  # safe_scalar_divide,
  # unsafe_norm,
  # unsafe_inv_norm,

  # rotations.jl
  quatErr,
  cheapRV,
  components,

  # vector_projections.jl
  vec_proj,
  vec_sub_vec_proj,
  a_dot_one_pad_b,

  # utility.jl
  first_3_of_6,
  last_3_of_6,
  getTop,
  unPad,
  onePad,
  zeroPad,
  mul_then_un_pad,
  one_pad_then_mul,
  weightPoly,
  make_pd_gains,
  rand_pd,

  # basic_dh.jl
  basic_dh,
  dh_R_t,
  povray_12,
  dh_vector_mul

  # # matrix_factor_derivatives.jl
  # cholesky_U_deravitive

end


text:
module RocketFunctionObservableTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "FunctionObservable" begin

    println("Testing: make")

    @testset begin
        source = make(Int) do actor
            next!(actor, 0)
            complete!(actor)
        end

        io = IOBuffer()

        show(io, source)

        printed = String(take!(io))

        @test occursin("FunctionObservable", printed)
        @test occursin(string(eltype(source)), printed)

        subscription = subscribe!(source, void())

        show(io, subscription)

        printed = String(take!(io))

        @test occursin("FunctionObservableSubscription", printed)

        unsubscribe!(subscription)
    end

    source1 = make(Int) do actor
        complete!(actor)
    end

    source2 = make(Float64) do actor
        setTimeout(100) do
            complete!(actor)
        end
    end

    source3 = make(Int) do actor
        next!(actor, 1)
        setTimeout(100) do
            complete!(actor)
        end
    end

    source4 = make(Int) do actor
        next!(actor, 1)
        setTimeout(100) do
            next!(actor, 2)
            complete!(actor)
        end
    end

    source5 = make(Int) do actor
        next!(actor, 1)
        setTimeout(100) do
            next!(actor, 2)
            error!(actor, "Error")
        end
    end

    run_testset([
        (
            source = source1,
            values = @ts(c),
            source_type = Int
        ),
        (
            source = source2,
            values = @ts(100 ~ c),
            source_type = Float64
        ),
        (
            source = source3,
            values = @ts([ 1 ] ~ 100 ~ c),
            source_type = Int
        ),
        (
            source = source4,
            values = @ts([ 1 ] ~ 100 ~ [ 2, c ]),
            source_type = Int
        ),
        (
            source = source5,
            values = @ts([ 1 ] ~ 100 ~ [ 2, e("Error") ]),
            source_type = Int
        )
    ])

end

end


text:
#
#   Gaussian Beam
#
#   V: 2019-04-05
#
#   Stefan Knirck
#

export gauss_profile

@doc raw"""
    gauss_profile(coords::CoordinateSystem; z = 0.000000001, omega0 = 0.1, f = 10e9)

Compute normalized 2D gaussian beam profile at given z position (waist at ```z = 0```).

# Arguments
- `coords`: contains X,Y and Kx,Ky coordinates
- `z`: Distance from beam waist
- `omega0` ```> 0```: Beam waist radius
- `f` ```> 0```: Frequency of monochromatic EM radiation
"""
function gauss_profile(coords::CoordinateSystem; z = 0.000000001, omega0 = 0.1, f = 10e9)
    radius = [sqrt(x^2 + y^2) for x in coords.X, y in coords.Y]

    lambda = wavelength(f)

    #Gaussian Beam profile:
    k = 2*pi/lambda
    zR = pi*omega0^2/lambda
    omega(z) = omega0*sqrt(1+(z/zR)^2)
    R(z)     = z*(1+(zR/z)^2)
    Psi(z)   = atan(z/zR)

   return omega0/omega(z) .* exp.(-radius.^2/omega(z)^2) .* exp.(-1im .* (k*z .+ k*radius.^2. /(2*R.(z)) .- Psi.(z)))
end


text:
import Base: muladd, fma

function fma(xhi::T, xlo::T, yhi::T, ylo::T, zhi::T, zlo::T) where {T<:AbstractFloat}
   chi, c1 = two_prod(xhi, yhi)
   t0 = xlo * ylo
   t1 = fma(xhi, ylo, t0)
   c2 = fma(xlo, yhi, t1)
   c3 = c1 + c2
   dhi, dlo = two_hilo_sum(chi, c3)

   shi, slo = two_sum(zhi, dhi)
   thi, tlo = two_sum(zlo, dlo)
   c = slo + thi
   vhi, vlo = two_hilo_sum(shi, c)
   w = tlo + vlo
   hi, lo = two_hilo_sum(vhi, w)
   return DoubleFloat{T}(hi, lo)
end

@inline function fma(x::DoubleFloat{T}, y::DoubleFloat{T}, z::DoubleFloat{T}) where {T<:AbstractFloat}
   return fma(x.hi, x.lo, y.hi, y.lo, z.hi, z.lo)
end

function fma(x::T, y::T, zhi::T, zlo::T) where {T<:AbstractFloat}
   chi, c1 = two_prod(x, y)
   shi, slo = two_sum(zhi, chi)
   thi, tlo = two_sum(zlo, c1)
   c = slo + thi
   vhi, vlo = two_hilo_sum(shi, c)
   w = tlo + vlo
   hi, lo = two_hilo_sum(vhi, w)
   return DoubleFloat{T}(hi, lo)
end

@inline function fma(x::T, y::T, z::DoubleFloat{T}) where {T<:AbstractFloat}
   return fma(x, y, z.hi, z.lo)
end

function fma(xhi::T, xlo::T, yhi::T, ylo::T, z::T) where {T<:AbstractFloat}
   chi, c1 = two_prod(xhi, yhi)
   t0 = xlo * ylo
   t1 = fma(xhi, ylo, t0)
   c2 = fma(xlo, yhi, t1)
   c3 = c1 + c2
   dhi, dlo = two_hilo_sum(chi, c3)
   shi, slo = two_sum(z, dhi)
   c = slo + dlo
   hi, lo = two_hilo_sum(shi, c)
   return DoubleFloat{T}(hi, lo)
end

@inline function fma(x::DoubleFloat{T}, y::DoubleFloat{T}, z::T) where {T<:AbstractFloat}
   return fma(x.hi, x.lo, y.hi, y.lo, z)
end

@inline function fma(x::DoubleFloat{T}, y::DoubleFloat{T}, z::Tuple{T,T}) where {T<:AbstractFloat}
   return fma(x.hi, x.lo, y.hi, y.lo, z[1], z[2])
end
@inline function fma(x::DoubleFloat{T}, y::Tuple{T,T}, z::DoubleFloat{T}) where {T<:AbstractFloat}
   return fma(x.hi, x.lo, y[1], y[2], z.hi, z.lo)
end
@inline function fma(x::Tuple{T,T}, y::DoubleFloat{T}, z::DoubleFloat{T}) where {T<:AbstractFloat}
   return fma(x[1], x[2], y.hi, y.lo, z.hi, z.lo)
end
@inline function fma(x::DoubleFloat{T}, y::Tuple{T,T}, z::Tuple{T,T}) where {T<:AbstractFloat}
   return fma(x.hi, x.lo, y[1], y[2], z[1], z[2])
end
@inline function fma(x::Tuple{T,T}, y::Tuple{T,T}, z::DoubleFloat{T}) where {T<:AbstractFloat}
   return fma(x[1], x[2], y[1], y[2], z.hi, z.lo)
end
@inline function fma(x::Tuple{T,T}, y::DoubleFloat{T}, z::Tuple{T,T}) where {T<:AbstractFloat}
   return fma(x[1], x[2], y.hi, y.lo, z[1], z[2])
end
@inline function fma(x::Tuple{T,T}, y::Tuple{T,T}, z::Tuple{T,T}) where {T<:AbstractFloat}
   return fma(x[1], x[2], y[1], y[2], z[1], z[2])
end

function muladd(xhi::T, xlo::T, yhi::T, ylo::T, zhi::T, zlo::T) where {T<:AbstractFloat}
   chi, c1 = two_prod(xhi, yhi)
   t0 = xlo * ylo
   c2 = fma(xlo, yhi, t0)
   c3 = c1 + c2
   dhi, dlo = two_hilo_sum(chi, c3)

   shi, slo = two_sum(zhi, dhi)
   thi, tlo = two_sum(zlo, dlo)
   c = slo + thi
   vhi, vlo = two_hilo_sum(shi, c)
   w = tlo + vlo
   hi, lo = two_hilo_sum(vhi, w)
   return DoubleFloat{T}(hi, lo)
end

@inline function muladd(x::DoubleFloat{T}, y::DoubleFloat{T}, z::DoubleFloat{T}) where {T<:AbstractFloat}
   return muladd(x.hi, x.lo, y.hi, y.lo, z.hi, z.lo)
end

function muladd(x::T, y::T, zhi::T, zlo::T) where {T<:AbstractFloat}
   chi, clo = two_prod(x, y)
   shi, slo = two_sum(zhi, chi)
   thi, tlo = two_sum(zlo, clo)
   c = slo + thi
   vhi, vlo = two_hilo_sum(shi, c)
   w = tlo + vlo
   hi, lo = two_hilo_sum(vhi, w)
   return DoubleFloat{T}(hi, lo)
end

@inline function muladd(x::T, y::T, z::DoubleFloat{T}) where {T<:AbstractFloat}
   return muladd(x, y, z.hi, z.lo)
end

function muladd(xhi::T, xlo::T, yhi::T, ylo::T, z::T) where {T<:AbstractFloat}
   chi, c1 = two_prod(xhi, yhi)
   t0 = xlo * ylo
   c2 = fma(xlo, yhi, t0)
   c3 = c1 + c2
   dhi, dlo = two_hilo_sum(chi, c3)

   shi, slo = two_sum(z, dhi)
   c = slo + dlo
   hi, lo = two_hilo_sum(shi, c)
   return DoubleFloat{T}(hi, lo)
end

@inline function muladd(x::DoubleFloat{T}, y::DoubleFloat{T}, z::T) where {T<:AbstractFloat}
   return muladd(x.hi, x.lo, y.hi, y.lo, z)
end

@inline function muladd(x::DoubleFloat{T}, y::DoubleFloat{T}, z::Tuple{T,T}) where {T<:AbstractFloat}
   return muladd(x.hi, x.lo, y.hi, y.lo, z[1], z[2])
end
@inline function muladd(x::DoubleFloat{T}, y::Tuple{T,T}, z::DoubleFloat{T}) where {T<:AbstractFloat}
   return muladd(x.hi, x.lo, y[1], y[2], z.hi, z.lo)
end
@inline function muladd(x::Tuple{T,T}, y::DoubleFloat{T}, z::DoubleFloat{T}) where {T<:AbstractFloat}
   return muladd(x[1], x[2], y.hi, y.lo, z.hi, z.lo)
end
@inline function muladd(x::DoubleFloat{T}, y::Tuple{T,T}, z::Tuple{T,T}) where {T<:AbstractFloat}
   return muladd(x.hi, x.lo, y[1], y[2], z[1], z[2])
end
@inline function muladd(x::Tuple{T,T}, y::Tuple{T,T}, z::DoubleFloat{T}) where {T<:AbstractFloat}
   return muladd(x[1], x[2], y[1], y[2], z.hi, z.lo)
end
@inline function muladd(x::Tuple{T,T}, y::DoubleFloat{T}, z::Tuple{T,T}) where {T<:AbstractFloat}
   return muladd(x[1], x[2], y.hi, y.lo, z[1], z[2])
end
@inline function muladd(x::Tuple{T,T}, y::Tuple{T,T}, z::Tuple{T,T}) where {T<:AbstractFloat}
   return muladd(x[1], x[2], y[1], y[2], z[1], z[2])
end