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

using Markdown
using InteractiveUtils

# ╔═╡ 882c47a6-13c5-11ec-2eee-5f38073a3824
begin
  using TikzPictures
  using LaTeXStrings
end

# ╔═╡ 0391208b-b2ab-4240-a681-37bc7ded21ab
let
	a = "A"
	"ab$(a)"
end

# ╔═╡ 52dfc51b-bcc7-44ed-93f3-5a616cb7ac5d
let
  TikzPicture(
    L"""
    \def\fillcolor{gray}
    \def\radius{3ex}
    \def\height{4.5ex};
    \draw[fill=\fillcolor,\fillcolor] (2*\radius, 0) circle (\radius);
    \node[\fillcolor] at (6*\radius, 0) {\Huge $\cdots$};
    \draw[fill=\fillcolor,\fillcolor] (10*\radius, 0) circle (\radius);
    """,
    #options="scale=1.5",
    preamble="")
end

# ╔═╡ 9734ee4b-1a80-4de0-8293-cfd7817663b6
let
  TikzPicture(
    L"""
    \def\fillcolor{gray}
    \def\radius{2ex}
    \def\height{4.5ex}
    \node (left parent) at (5*\radius, 0) [circle, draw=gray, fill=gray, inner sep=\radius] {};
    \node [white, above] at (left parent.north) {\large $2^{n}$};
    \node (right parent) at (10*\radius, 0) [circle, draw=gray, fill=gray, inner sep=\radius] {};
    \node[\fillcolor] at (6*\radius, 0) {\Huge $\cdots$};
    """,
    #options="scale=1",
    preamble="")
end

# ╔═╡ 06e834ff-8904-46a6-b554-014ca15541c9
let
  TikzPicture(
    L"""
    \node {parent} [circle, fill=green, sibling distance = 2.5cm, green]
        child [red] {node {child 1} edge from parent [dashed, green]}
        child {node [red] {child 2}
        child {node [xshift = -1cm] {grandchild 1}}
        child {node [yshift = .5cm] {grandchild 2}}
        edge from parent [green] node [right, brown] {x}};
    """,
    #options="scale=1",
    preamble="")
end

# ╔═╡ 35950c6b-d4d0-4be2-8b41-786146187c6f
let
  TikzPicture(
    L"""
    \node {parent} [sibling distance = 2.5cm, green]
      child [red] {node {child 1} edge from parent [dashed, green]}
      child {node [red] {child 2}
      child {node [xshift = -1cm] {grandchild 1}}
      child {node [yshift = .5cm] {grandchild 2}}
      edge from parent [green] node [right, brown] {x}};
    """,
    #options="scale=1",
    preamble="")
end

# ╔═╡ acab0352-d24d-4300-92f5-65beb1e62e00
let
  TikzPicture(
    L"""
    \node [circle, draw=green, inner sep=3ex] {parent} [sibling distance = 1.5cm, green]
      child [red] {node {child 1} edge from parent [dashed, green]}
      child {node [red] {child 2}
      child {node [xshift = -1cm] {grandchild 1}}
      child {node [yshift = .5cm] {grandchild 2}}
      edge from parent [green] node [right, brown] {x}};
    """,
    #options="scale=1",
    preamble="")
end

# ╔═╡ 09f4b3f1-e7f0-4ab1-9597-0fdbaaf8e615
let
  TikzPicture(
    L"""
    \def\radius{1em}
    \def\sep{10em}
    \def\color{gray}
    \def\sib_dist{7em}
    \node (left) [circle, fill=\color, inner sep=\radius] {} [sibling distance=\sib_dist]
      child [draw=white, line width=0.5mm] { node (left left) [circle, fill=\color, inner sep=\radius] {} edge from parent }
      child [draw=white, line width=0.5mm] { node (left right) [circle, fill=\color, inner sep=\radius] {} edge from parent };

    \node [white, above] at (left.north) {\large $2^{n}$};
    \node [white, above] at (left left.north) {\large $2^{n+1}$};
    \node [white, above] at (left right.north) {\large $2^{n+1}$};

    \node (ellipsis) [\color] at (\sep, 0) {\Huge $\cdots$};
    \node (right) at (2*\sep, 0) [circle, fill=\color, inner sep=\radius] {} [sibling distance=\sib_dist]
      child [draw=white, line width=0.5mm] { node (right left) [circle, fill=\color, inner sep=\radius] {} edge from parent }
      child [draw=white, line width=0.5mm] { node (right right) [circle, fill=\color, inner sep=\radius] {} edge from parent };

    \node [white, above] at (right.north) {\large $2^{n}+k$};
    \node [rectangle, draw=white, white, above] at (right left.north) {\large\textbf ?};
    \node [rectangle, draw=white, white, above] at (right right.north) {\large\textbf ?};
    """,
    options="scale=1.0",
    preamble="")
end

# ╔═╡ 4b4afc4c-ceb3-423f-b4f2-d7688350d812
  # morceaux = TikzPicture(
  #   L"""
  #   \def\d{0.8}
  #   \def\L{10}
  #   \draw[orange, ultra thick] (0,0) -- (\L,0);
  #   \draw[purple, ultra thick] (0,-0.2) -- (0,0.2) node [anchor=south] {$\frac{0}{M-1}$};
  #   \draw[purple, ultra thick] (\d,-0.2) -- (\d,0.2) node [anchor=south] {$\frac{1}{M-1}$};
  #   \draw[purple, ultra thick] (2*\d,-0.2) -- (2*\d,0.2) node [anchor=south] {$\frac{2}{M-1}$};
  #   \draw[purple, ultra thick] (\L/2,-0.2) -- (\L/2,0.2) node [anchor=south] {$\cdots$};
  #   %node (\L/2,0.2) [anchor=south] {$\cdots$};
  #   \draw[purple, ultra thick] (\L,-0.2) -- (\L,0.2) node [anchor=south] {$\frac{M-1}{M-1}$};
  #   \draw[purple, ultra thick] (\L-\d,-0.2) -- (\L-\d,0.2) node [anchor=south] {$\frac{M-2}{M-1}$};
  #   """,
  #   options="scale=1.5",
  #   preamble="",
  # )

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
TikzPictures = "37f6aa50-8035-52d0-81c2-5a1d08754b2d"

[compat]
LaTeXStrings = "~1.2.1"
TikzPictures = "~3.4.1"
"""

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

[[Artifacts]]
deps = ["Pkg"]
git-tree-sha1 = "c30985d8821e0cd73870b17b0ed0ce6dc44cb744"
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
version = "1.3.0"

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

[[Bzip2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "c3598e525718abcc440f69cc6d5f60dda0a1b61e"
uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
version = "1.0.6+5"

[[Cairo_jll]]
deps = ["Artifacts", "Bzip2_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "e2f47f6d8337369411569fd45ae5753ca10394c6"
uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
version = "1.16.0+6"

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

[[Expat_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "1402e52fcda25064f51c77a9655ce8680b76acf0"
uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
version = "2.2.7+6"

[[Fontconfig_jll]]
deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "35895cf184ceaab11fd778b4590144034a167a2f"
uuid = "a3f928ae-7b40-5064-980b-68af3947d34b"
version = "2.13.1+14"

[[FreeType2_jll]]
deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "cbd58c9deb1d304f5a245a0b7eb841a2560cfec6"
uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7"
version = "2.10.1+5"

[[Gettext_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "8c14294a079216000a0bdca5ec5a447f073ddc9d"
uuid = "78b55507-aeef-58d4-861c-77aaff3498b1"
version = "0.20.1+7"

[[Glib_jll]]
deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "04690cc5008b38ecbdfede949220bc7d9ba26397"
uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
version = "2.59.0+4"

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

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JpegTurbo_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "9aff0587d9603ea0de2c6f6300d9f9492bbefbd3"
uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
version = "2.0.1+3"

[[LZO_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "f128cd6cd05ffd6d3df0523ed99b90ff6f9b349a"
uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac"
version = "2.10.0+3"

[[LaTeXStrings]]
git-tree-sha1 = "c7f1c695e06c01b95a67f0cd1d34994f3e7db104"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.2.1"

[[LibGit2]]
deps = ["Printf"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[Libffi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "a2cd088a88c0d37eef7d209fd3d8712febce0d90"
uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490"
version = "3.2.1+4"

[[Libgcrypt_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll", "Pkg"]
git-tree-sha1 = "b391a18ab1170a2e568f9fb8d83bc7c780cb9999"
uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4"
version = "1.8.5+4"

[[Libgpg_error_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "ec7f2e8ad5c9fa99fc773376cdbc86d9a5a23cb7"
uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8"
version = "1.36.0+3"

[[Libiconv_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "cba7b560fcc00f8cd770fa85a498cbc1d63ff618"
uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
version = "1.16.0+8"

[[Libmount_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "51ad0c01c94c1ce48d5cad629425035ad030bfd5"
uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9"
version = "2.34.0+3"

[[Libtiff_jll]]
deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Pkg", "Zlib_jll", "Zstd_jll"]
git-tree-sha1 = "291dd857901f94d683973cdf679984cdf73b56d0"
uuid = "89763e89-9b03-5906-acba-b20f662cd828"
version = "4.1.0+2"

[[Libuuid_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "f879ae9edbaa2c74c922e8b85bb83cc84ea1450b"
uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700"
version = "2.34.0+7"

[[LittleCMS_jll]]
deps = ["JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pkg"]
git-tree-sha1 = "e6ea89d915cdad8d264f7f9158c6664f879edcde"
uuid = "d3a379c0-f9a3-5b72-a4c0-6bf4d2e8af0f"
version = "2.9.0+0"

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

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

[[OpenJpeg_jll]]
deps = ["Libdl", "Libtiff_jll", "LittleCMS_jll", "Pkg", "libpng_jll"]
git-tree-sha1 = "e330ffff1c6a593fa44cc40c29900bee82026406"
uuid = "643b3616-a352-519d-856d-80112ee9badc"
version = "2.3.1+0"

[[PCRE_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "1b556ad51dceefdbf30e86ffa8f528b73c7df2bb"
uuid = "2f80f16e-611a-54ab-bc61-aa92de5b98fc"
version = "8.42.0+4"

[[Pixman_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "6a20a83c1ae86416f0a5de605eaea08a552844a3"
uuid = "30392449-352a-5448-841d-b1acce4e97dc"
version = "0.40.0+0"

[[Pkg]]
deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[Poppler_jll]]
deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "OpenJpeg_jll", "Pkg", "libpng_jll"]
git-tree-sha1 = "e11443687ac151ac6ef6699eb75f964bed8e1faa"
uuid = "9c32591e-4766-534b-9725-b71a8799265b"
version = "0.87.0+2"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

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

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

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

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

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

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[TOML]]
deps = ["Dates"]
git-tree-sha1 = "44aaac2d2aec4a850302f9aa69127c74f0c3787e"
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
version = "1.0.3"

[[Tectonic]]
deps = ["Pkg"]
git-tree-sha1 = "acf12eccb390a78653ee805cd527898f01f78a85"
uuid = "9ac5f52a-99c6-489f-af81-462ef484790f"
version = "0.6.1"

[[TikzPictures]]
deps = ["LaTeXStrings", "Poppler_jll", "Requires", "Tectonic"]
git-tree-sha1 = "a08671c0979063a437378f6410bb75a465f3cd1c"
uuid = "37f6aa50-8035-52d0-81c2-5a1d08754b2d"
version = "3.4.1"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

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

[[XML2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "be0db24f70aae7e2b89f2f3092e93b8606d659a6"
uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
version = "2.9.10+3"

[[XSLT_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "2b3eac39df218762d2d005702d601cd44c997497"
uuid = "aed1982a-8fda-507f-9586-7b0439959a61"
version = "1.1.33+4"

[[Xorg_libX11_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll", "Xorg_xtrans_jll"]
git-tree-sha1 = "5be649d550f3f4b95308bf0183b82e2582876527"
uuid = "4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"
version = "1.6.9+4"

[[Xorg_libXau_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4e490d5c960c314f33885790ed410ff3a94ce67e"
uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec"
version = "1.0.9+4"

[[Xorg_libXdmcp_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4fe47bd2247248125c428978740e18a681372dd4"
uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05"
version = "1.1.3+4"

[[Xorg_libXext_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "b7c0aa8c376b31e4852b360222848637f481f8c3"
uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3"
version = "1.3.4+4"

[[Xorg_libXrender_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "19560f30fd49f4d4efbe7002a1037f8c43d43b96"
uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa"
version = "0.9.10+4"

[[Xorg_libpthread_stubs_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "6783737e45d3c59a4a4c4091f5f88cdcf0908cbb"
uuid = "14d82f49-176c-5ed1-bb49-ad3f5cbd8c74"
version = "0.1.0+3"

[[Xorg_libxcb_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"]
git-tree-sha1 = "daf17f441228e7a3833846cd048892861cff16d6"
uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b"
version = "1.13.0+3"

[[Xorg_xtrans_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "79c31e7844f6ecf779705fbc12146eb190b7d845"
uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10"
version = "1.4.0+3"

[[Zlib_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "320228915c8debb12cb434c59057290f0834dbf6"
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
version = "1.2.11+18"

[[Zstd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "2c1332c54931e83f8f94d310fa447fd743e8d600"
uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
version = "1.4.8+0"

[[libpng_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "6abbc424248097d69c0c87ba50fcb0753f93e0ee"
uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
version = "1.6.37+6"
"""

# ╔═╡ Cell order:
# ╠═882c47a6-13c5-11ec-2eee-5f38073a3824
# ╠═0391208b-b2ab-4240-a681-37bc7ded21ab
# ╟─52dfc51b-bcc7-44ed-93f3-5a616cb7ac5d
# ╠═9734ee4b-1a80-4de0-8293-cfd7817663b6
# ╠═06e834ff-8904-46a6-b554-014ca15541c9
# ╠═35950c6b-d4d0-4be2-8b41-786146187c6f
# ╠═acab0352-d24d-4300-92f5-65beb1e62e00
# ╠═09f4b3f1-e7f0-4ab1-9597-0fdbaaf8e615
# ╟─4b4afc4c-ceb3-423f-b4f2-d7688350d812
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
# This file is to get the Cl and Cd of Blade Element

#the excel interpolation method
function cla(α,ma)
  # NACA0012
  x = α/π*180
  if ma<0.1
    cla = 2e-17*x^6 + 1e-10*x^5 - 8e-13*x^4 - 4e-06*x*3 + 1e-08*x^2 + 0.0351*x + 1e-06
  elseif ma<0.2
    cla = 9e-18*x^6 + 1e-10*x^5 - 5e-13*x^4 - 4e-06*x^3 + 6e-09*x^2 + 0.0351*x - 3e-06
  elseif ma<0.3
    cla = 2e-17*x^6 + 1e-10*x^5 - 8e-13*x^4 - 4e-06*x^3 + 1e-08*x^2 + 0.0351*x + 1e-06
  elseif ma<0.4
    cla = -1e-15*x^6 + 1e-10*x^5 + 6e-11*x^4 - 4e-06*x^3 - 8e-07*x^2 + 0.0349*x + 0.0025
  elseif ma<0.5
    cla = -3e-17*x^6 + 1e-10*x^5 + 1e-12*x^4 - 4e-06*x^3 - 3e-09*x^2 + 0.0349*x - 0.0001
  elseif ma<0.6
    cla = -6e-17*x^6 + 1e-10*x^5 + 2e-12*x^4 - 4e-06*x^3 - 9e-09*x^2 + 0.0347*x - 0.0002
  elseif ma<0.7
    cla = -6e-17*x^6 + 1e-10*x^5 + 2e-12*x^4 - 4e-06*x^3 - 9e-09*x^2 + 0.0345*x - 0.0002
  elseif ma<0.8
    cla = -6e-17*x^6 + 1e-10*x^5 + 2e-12*x^4 - 4e-06*x^3 - 9e-09*x^2 + 0.0345*x - 0.0002
  elseif ma<0.9
    cla = 2e-16*x^6 + 9e-11*x^5 - 7e-12*x^4 - 4e-06*x^3 + 4e-08*x^2 + 0.0313*x + 0.0004
  elseif ma<1.0
    cla = 1e-16*x^6 + 9e-11*x^5 - 5e-12*x^4 - 4e-06*x^3 - 4e-09*x^2 + 0.03*x + 0.0007
  else
    cla = 1e-16*x^6 + 9e-11*x^5 - 5e-12*x^4 - 4e-06*x^3 - 7e-09*x^2 + 0.03*x + 0.0007
  end
  return cla
end

#the simple interpolation method
function liftcoffi(α,Ma)
  if abs(α)>(20.0*π/180.0)
    Cl = 0.0
  else
    Cl = 5.73*α+0.085
  end
  return Cl
end

function dragcoffi(α,Ma)
  if abs(α)<(3.0*π/180.0)
    Cd = 0.0076
  elseif abs(α)<(5.0*π/180.0)
    Cd = 0.008
  elseif abs(α)<(6.0*π/180.0)
    Cd = 0.0086
  elseif abs(α)<(7.0*π/180.0)
    Cd = 0.0093
  elseif abs(α)<(8.0*π/180.0)
    Cd = 0.0103
  elseif abs(α)<(9.0*π/180.0)
    Cd = 0.012
  elseif abs(α)<(10.0*π/180.0)
    Cd = 0.0146
  elseif abs(α)<(11.0*π/180.0)
    Cd = 0.0177
  elseif abs(α)<(12.0*π/180.0)
    Cd = 0.0215
  elseif abs(α)<(13.0*π/180.0)
    Cd = 0.0266
  elseif abs(α)<(14.0*π/180.0)
    Cd = 0.055
  elseif abs(α)<(15.0*π/180.0)
    Cd = 0.125
  elseif abs(α)<(16.0*π/180.0)
    Cd = 0.21
  else
    Cd = 0.5
  end

  Cd = Cd*1.05  #The real blade correction
  return Cd
end

function clcdget(α_aero,vall_r)
  Ma = Array{Float64}(NR,Nbe)
  for k in 1:NR #take Mach Number into consideration
    for i in 1:Nbe
      Ma[k,i] = norm(vall_r[k,i])/v_sound
    end
  end
  clift = Array{Float64}(NR,Nbe)
  cdrag = Array{Float64}(NR,Nbe)
  for k in 1:NR #solve the local clcd
    for i in 1:Nbe
      clift[k,i] = liftcoffi(α_aero[k,i],Ma[k,i])
      cdrag[k,i] = dragcoffi(α_aero[k,i],Ma[k,i])
    end
  end
  return clift,cdrag,Ma
end

ext: jl

max_stars_repo_name: dingzhw/RotorTrim

max_issues_repo_name: helicopterliu/RotorTrim

max_forks_repo_name: helicopterliu/RotorTrim


text:
"""
    gridRectangular(xdim, ydim)

Create coordinates of all neurons on a rectangular SOM.

The return-value is an array of size (Number-of-neurons, 2) with
x- and y- coordinates of the neurons in the first and second
column respectively.
The distance between neighbours is 1.0.
The point of origin is bottom-left.
The first neuron sits at (0,0).

# Arguments
- `xdim`: number of neurons in x-direction
- `ydim`: number of neurons in y-direction
"""
function gridRectangular(xdim, ydim)

    grid = zeros(Float64, (xdim*ydim, 2))
    for ix in 1:xdim
        for iy in 1:ydim

            grid[ix+(iy-1)*xdim, 1] = ix-1
            grid[ix+(iy-1)*xdim, 2] = iy-1
        end
    end
    return grid
end

"""
    gaussianKernel(x, r::Float64)

Return Gaussian(x) for μ=0.0 and σ = r/3.
(a value of σ = r/3 makes the training results comparable between different kernels
for same values of r).

# Arguments

"""
function gaussianKernel(x, r::Float64)

    return Distributions.pdf.(Distributions.Normal(0.0,r/3), x)
end

"""
    distMatrix(grid::Array, toroidal::Bool)::Array{Float64, 2}

Return the distance matrix for a non-toroidal or toroidal SOM.

# Arguments
- `grid`: coordinates of all neurons as generated by one of the `grid-`functions
with x-coordinates in 1st column and y-coordinates in 2nd column.
- `toroidal`: true for a toroidal SOM.
"""
function distMatrix(grid::Array, toroidal::Bool)::Array{Float64, 2}

    X = 1
    Y = 2
    xdim = maximum(grid[:,X]) - minimum(grid[:,X]) + 1.0
    ydim = maximum(grid[:,Y]) - minimum(grid[:,Y]) + 1.0

    numNeurons = size(grid,1)

    dm = zeros(Float64, (numNeurons,numNeurons))
    for i in 1:numNeurons
        for j in 1:numNeurons
            Δx = abs(grid[i,X] - grid[j,X])
            Δy = abs(grid[i,Y] - grid[j,Y])

            if toroidal
                Δx = min(Δx, xdim-Δx)
                Δy = min(Δy, ydim-Δy)
            end

            dm[i,j] = √(Δx^2 + Δy^2)
        end
    end
    # show(STDOUT, "text/plain",  grid)
    # show(STDOUT, "text/plain",  dm)
    return dm
end

"""
    classFrequencies(som::Som, data, classes)

Return a DataFrame with class frequencies for all neurons.

# Arguments:
- `som`: a trained SOM
- `data`: data with row-wise samples and class information in each row
- `classes`: Name of column with class information.

Data must have the same number of dimensions as the training dataset.
The column with class labels is given as `classes` (name or index).
Returned DataFrame has the columns:
* X-, Y-indices and index: of winner neuron for every row in data
* population: number of samples mapped to the neuron
* frequencies: one column for each class label.
"""
function classFrequencies(som::Som, data, classes)

    if size(data,2) != size(som.codes,2) + 1
        println("    data: $(size(data,2)-1), codes: $(size(som.codes,2))")
        error(SOM_ERRORS[:ERR_COL_NUM])
    end

    x = deepcopy(data)
    deletecols!(x, classes)
    classes = data[classes]
    vis = visual(som.codes, x)

    df = makeClassFreqs(som, vis, classes)
    return df
end

"""
    visual(codes::Array{Float64,2}, x::Array{Float64,2})

Return the index of the winner neuron for each training pattern
in x (row-wise).

# Arguments
- `codes`: Codebook
- `x`: training Data
"""
function visual(codes::Array{Float64,2}, x::Array{Float64,2})

    vis = zeros(Int, size(x,1))
    for i in 1:size(x,1)

        vis[i] = findBmu(codes, x[i, : ])
    end

    return vis
end

"""
    makePopulation(numCodes, vis)

Return a vector of neuron populations.

# Arguments
- `numCodes`: total number of neurons
- `vis`: index of the winner neuron for each training pattern in x
"""
function makePopulation(numCodes, vis)

    population = zeros(Int, numCodes)
    for i in 1:size(vis,1)
        population[vis[i]] += 1
    end

    return population
end

"""
    makeClassFreqs(som, vis, classes)

Return a DataFrame with class frequencies for all neurons.

# Arguments
- `som`: a trained SOM
- `vis`: index of the winner neuron for each training pattern in x
- `classes`: name of column with class information
"""
function makeClassFreqs(som, vis, classes)

    # count classes and construct DataFrame:
    #
    classLabels = sort(unique(classes))
    classNum = size(classLabels,1)

    cfs = DataFrame(index = 1:som.numCodes)
    cfs[:X] = som.indices[:X]
    cfs[:Y] = som.indices[:Y]

    cfs[:Population] = zeros(Int, som.numCodes)

    for class in classLabels
        cfs[Symbol(class)] = zeros(Float64, som.numCodes)
    end

    # loop vis and count:
    #
    for i in 1:size(vis,1)

        cfs[vis[i], :Population] += 1
        class = Symbol(classes[i])
        cfs[vis[i], class] += 1
    end

    # make frequencies from counts:
    #
    for i in 1:size(cfs,1)

        counts = [cfs[i, col] for col in 5:size(cfs, 2)]
        total = cfs[i,:Population]
        if total == 0
            freqs = counts * 0.0
        else
            freqs = counts ./ total
        end

        for c in 1:classNum
            class = Symbol(classLabels[c])
            cfs[i,class] = freqs[c]
        end
    end

    return cfs
end

"""
    findBmu(codes::Array{Float64,2}, sample::Array{Float64,1})::Int64

Find the best matching unit for a given vector, row_t, in the SOM

Returns: Index of bmu
Best Matching Unit and bmuIdx is the index of this vector in the SOM

# Arguments
- `codes`: 2D-array of codebook vectors. One vector per row
- `sample`: row in dataset / trainingsset

"""
function findBmu(codes::Array{Float64,2}, sample::Array{Float64,1})::Int64

    x = Distances.colwise(Euclidean(), sample, codes')

    return argmin(x)

end

"""
    normTrainData(x::Array{Float64,2}, normParams)

Normalise every column of training data with the params.

# Arguments
- `x`: DataFrame with training Data
- `normParams`: Shift and scale parameters for each attribute column.
"""
function normTrainData(x::Array{Float64,2}, normParams)

    for i in 1:size(x,2)
        x[:,i] = (x[:,i] .- normParams[1,i]) ./ normParams[2,i]
    end

    return x
end

"""
    normTrainData(train::Array{Float64, 2}, norm::Symbol)

Normalise every column of training data.

# Arguments
- `train`: DataFrame with training Data
- `norm`: type of normalisation; one of `minmax, zscore, none`
"""
function normTrainData(train::Array{Float64, 2}, norm::Symbol)

    normParams = zeros(2, size(train,2))

    if  norm == :minmax
        for i in 1:size(train,2)
            normParams[1,i] = minimum(train[:,i])
            normParams[2,i] = maximum(train[:,i]) - minimum(train[:,i])
        end
    elseif norm == :zscore
        for i in 1:size(train,2)
            normParams[1,i] = mean(train[:,i])
            normParams[2,i] = std(train[:,i])
        end
    else
        for i in 1:size(train,2)
            normParams[1,i] = 0.0  # shift
            normParams[2,i] = 1.0  # scale
        end
    end

    # do the scaling:
    if norm == :none
        x = train
    else
        x = normTrainData(train, normParams)
    end

    return x, normParams
end

"""
    convertTrainingData(data)::Array{Float64,2}

Converts the training data to an Array of type Float64.

# Arguments:
- `data`: Data to be converted

"""
function convertTrainingData(data)::Array{Float64,2}

    if typeof(data) == DataFrame
        train = convert(Matrix{Float64}, data)

    elseif typeof(data) != Matrix{Float64}
        try
            train = convert(Matrix{Float64}, data)
        catch ex
            Base.showerror(STDERR, ex, backtrace())
            error("Unable to convert training data to Array{Float64,2}!")
        end
    else
        train = data
    end

    return train
end

prettyPrintArray(arr) = println("$(show(IOContext(STDOUT, limit=true), "text/plain", arr))")

"""
    checkDir()

Checks if the `pwd()` is the `/test` directory, and if not it changes to it.

"""
function checkDir()

    files = readdir()
    if !in("runtests.jl", files)
        cd(dirname(dirname(pathof(GigaSOM))))
    end
end


text:
using Documenter
using Stopping

makedocs(
  sitename = "Stopping.jl",
  format = Documenter.HTML(
    assets = ["assets/style.css"],
    prettyurls = get(ENV, "CI", nothing) == "true",
  ),
  modules = [Stopping],
  pages = [
    "Home" => "index.md",
    "API" => "api.md",
    "Stopping's ID" => "idcard.md",
    "State's ID" => "idcard-state.md",
    "Meta's ID" => "idcard-stoppingmeta.md",
    "Optimality in Stopping" => "howstopcheckoptimality.md",
    "Stopping in action" => "example-basic-Newton.md",
    "Stop remote control" => "idcard-stopremote.md",
    "Stopping workflow" => "stop-workflow.md",
    "Speak to stopping" => "speak-to-stopping.md",
    "NLPStopping" => "nlpstopping.md",
    "LAStopping" => "lastopping.md",
    "Readme" => "index_tuto.md",
    "How to State" => "howtostate.md",
    "How to State for NLPs" => "howtostate-nlp.md",
    "How to Stop" => "howtostop.md",
    "How to Stop 2" => "howtostop-2.md",
    "How to Stop for NLPs" => "howtostop-nlp.md",
    "Solve linear algebra" => "linear-algebra.md",
    "Use a buffer function" => "buffer.md",
    "A fixed point algorithm" => "fixed-point.md",
    "Backtracking linesearch algorithm" => "backls.md",
    "Unconstrained optimization algorithm" => "uncons.md",
    "Active set algorithm" => "active-set.md",
    "Quadratic penalty algorithm" => "penalty.md",
    "Run optimization algorithms" => "run-optimsolver.md",
    "Benchmark optimization algorithms" => "benchmark.md",
    "Overfitting" => "overfitting.md",
    "Checkpointing" => "checkpointing.md",
    "Mix algorithms" => "gradient-lbfgs.md",
  ],
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(repo = "github.com/vepiteski/Stopping.jl")
#https://juliadocs.github.io/Documenter.jl/stable/man/hosting/ ?


text:
#!/usr/bin/env julia

"""
Fixing VCF files generated from variant calling
to load into IGV
"""

function fixHeader(vcf=ARGS[1])
    file = open(vcf)
    lines = readlines(file)
    info = "##fileformat=VCFv4.1
##fileDate=20151202
##source=Variant Calling
##reference=/wrk/annalam/organisms/homo_sapiens/hg19_nochr.fa
##INFO=Juliette
##INFO=--alt-reads=2 --alt-frac=0.05"
    println(info)
    header = "#" * shift!(lines)
    println(header)
    for i in lines#[1:5]
        println(chomp(i))
    end
end

fixHeader()


text:
using .AxisArrays

function rcopy(::Type{AxisArray}, r::Ptr{S}) where {S<:VectorSxp}
    dnames = getattrib(r, Const.DimNamesSymbol)
    isnull(dnames) && error("r has no dimnames")
    dsym = rcopy(Array{Symbol}, getnames(dnames))
    for i in 1:length(dsym)
        if dsym[i] == Symbol("")
            dsym[i] = i == 1 ? (:row) : i == 2 ? (:col) : i == 3 ? (:page) : Symbol(:dim_, i)
        end
    end
    as = Any[]
    for (i,n) in enumerate(dnames)
        push!(as, Axis{dsym[i]}(isnull(n) ? (1:(size(r)[i])) : rcopy(n)))
    end
    AxisArray(rcopy(AbstractArray, r), as...)
end

function sexp(::Type{RClass{:array}}, aa::AxisArray)
    rv = protect(sexp(aa.data))
    try
        d = OrderedDict(
            k => v.val for (k, v) in zip(axisnames(aa), AxisArrays.axes(aa)))
        setattrib!(rv, Const.DimSymbol, collect(size(aa)))
        setattrib!(rv, Const.DimNamesSymbol, d)
    finally
        unprotect(1)
    end
    rv
end

# AxisArray
sexpclass(v::AxisArray) = RClass{:array}


text:
const CuDense{ElT,VecT} = Dense{ElT,VecT} where {VecT<:CuVector}
const CuDenseTensor{ElT,N,StoreT,IndsT} = Tensor{ElT,N,StoreT,IndsT} where {StoreT<:CuDense}

function Dense{T,SA}(x::Dense{T,SB}) where {T<:Number,SA<:CuArray,SB<:Array}
  return Dense{T,SA}(CuArray(x))
end
function Dense{T,SA}(x::Dense{T,SB}) where {T<:Number,SA<:Array,SB<:CuArray}
  return Dense{T,SA}(collect(x.data))
end
Dense{T,S}(size::Integer) where {T,S<:CuArray{<:T}} = Dense{T,S}(CUDA.zeros(T, size))
function Dense{T,S}(x::T, size::Integer) where {T,S<:CuArray{<:T}}
  arr = CuArray{T}(undef, size)
  fill!(arr, x)
  return Dense{T,S}(arr)
end
cpu(x::CuDense{T}) where {T<:Number} = Dense(collect(x.data))
function Base.complex(::Type{Dense{ElT,VT}}) where {ElT,VT<:CuArray}
  return Dense{complex(ElT),CuVector{complex(ElT)}}
end

similartype(::Type{<:CuArray{<:Any,N}}, eltype::Type) where {N} = CuArray{eltype,N}
function NDTensors.similar(::Type{<:CuArray{T}}, dims) where {T}
  return CuArray{T,length(dims)}(undef, dims)
end

CuArray(x::CuDense{ElT}) where {ElT} = CuVector{ElT}(data(x))
function CuArray{ElT,N}(x::CuDenseTensor{ElT,N}) where {ElT,N}
  return CuArray{ElT,N}(reshape(data(store(x)), dims(inds(x))...))
end
CuArray(x::CuDenseTensor{ElT,N}) where {ElT,N} = CuArray{ElT,N}(x)

*(D::Dense{T,AT}, x::S) where {T,AT<:CuArray,S<:Number} = Dense(x .* data(D))

function Base.:(==)(
  ::Type{<:CuDense{ElT1,CVec1}}, ::Type{<:CuDense{ElT2,CVec2}}
) where {ElT1,ElT2,CVec1,CVec2}
  return (ElT1 == ElT2)
end
Base.getindex(D::CuDense{<:Number}) = collect(data(D))[]
Base.getindex(D::CuDenseTensor{<:Number,0}) = store(D)[]
LinearAlgebra.norm(T::CuDenseTensor) = norm(data(store(T)))

# This is for type promotion for Scalar*Dense
function Base.promote_rule(
  ::Type{<:Dense{ElT1,CuVector{ElT1}}}, ::Type{ElT2}
) where {ElT1,ElT2<:Number}
  ElR = promote_type(ElT1, ElT2)
  VecR = CuVector{ElR}
  return Dense{ElR,VecR}
end

function Base.permutedims(T::CuDenseTensor{<:Number,N}, perm::NTuple{N,Int}) where {N}
  Tp = NDTensors.similar(T, ITensors.NDTensors.permute(inds(T), perm))
  #Tp = permute(T,perm; always_copy=true)
  permute!(Tp, T)
  return Tp
end

function Base.permutedims!(
  R::CuDenseTensor{<:Number,N}, T::CuDenseTensor{<:Number,N}, perm::NTuple{N,Int}
) where {N}
  return permutedims!!(R, T, perm)
end

function permutedims!!(
  B::Tensor{ElT,N,StoreT,IndsB},
  A::Tensor{ElT,N,StoreT,IndsA},
  perm::NTuple{N,Int},
  f::Function=(r, t) -> permute!(r, t),
) where {N,ElT,IndsB,IndsA,StoreT<:CuDense{ElT}}
  Ais = inds(A)
  Bis = ITensors.NDTensors.permute(inds(A), perm)
  B = f(B, A)
  return B
end

function Base.similar(::Type{<:CuDenseTensor{ElT}}, inds::Dims) where {ElT}
  storage_arr = CuVector{ElT}(undef, dim(inds))
  return Tensor(Dense(storage_arr), inds)
end

function Base.similar(::Type{<:CuDenseTensor{ElT}}, inds) where {ElT}
  storage_arr = CuVector{ElT}(undef, dim(inds))
  return Tensor(Dense(storage_arr), inds)
end

function outer!(R::CuDenseTensor, T1::CuDenseTensor, T2::CuDenseTensor)
  R_dat = vec(array(T1)) * transpose(vec(array(T2)))
  copyto!(data(store(R)), vec(R_dat))
  inds_outer = unioninds(inds(T1), inds(T2))
  return R
end

function contract!!(
  R::CuDenseTensor{<:Number,NR},
  labelsR::NTuple{NR},
  T1::CuDenseTensor{<:Number,N1},
  labelsT1::NTuple{N1},
  T2::CuDenseTensor{<:Number,N2},
  labelsT2::NTuple{N2},
) where {NR,N1,N2}
  if N1 == 0
    # TODO: replace with an add! function?
    # What about doing `R .= T1[] .* PermutedDimsArray(T2,perm)`?
    perm = getperm(labelsR, labelsT2)
    newT2 = Tensor(Dense(data(store(T1)) .* data(store(T2))), inds(T2))
    permute!(R, newT2)
  elseif N2 == 0
    perm = getperm(labelsR, labelsT1)
    newT1 = Tensor(Dense(data(store(T2)) .* data(store(T1))), inds(T1))
    permute!(R, newT1)
  elseif N1 + N2 == NR
    # TODO: permute T1 and T2 appropriately first (can be more efficient
    # then permuting the result of T1⊗T2)
    # TODO: implement the in-place version directly
    R = outer!!(R, T1, T2)
    inds_outer = unioninds(inds(T1), inds(T2))
    R = Tensor(store(R), inds_outer)
  else
    R = _contract!!(R, labelsR, T1, labelsT1, T2, labelsT2)
  end
  return R
end

function permutedims!!(
  B::CuDenseTensor{ElT,0},
  A::CuDenseTensor{ElT,0},
  perm::NTuple{0,Int},
  f=(r, t) -> permute!(r, t),
) where {ElT<:Number}
  Cs = f(B, A)
  return Tensor(Dense(vec(Cs)), IndexSet{0}())
end

function permutedims!!(
  B::CuDenseTensor{ElT,N},
  A::CuDenseTensor{ElT,0},
  perm::NTuple{N,Int},
  f=(r, t) -> permute!(r, t),
) where {N,ElT<:Number}
  Cis = ITensors.NDTensors.permute(inds(B), perm)
  Cs = f(B, A)
  return Tensor(Dense(vec(Cs)), Cis)
end

function _contract_scalar!(
  R::CuDenseTensor{ElR},
  labelsR,
  T1::Number,
  labelsT1,
  T2::Number,
  labelsT2,
  α=one(ElR),
  β=zero(ElR),
) where {ElR}
  if iszero(β)
    copyto!(data(R), [α * T1 * T2])
  elseif iszero(α)
    copyto!(data(R), β .* data(R))
  else
    copyto!(data(R), [α * T1 * T2] .+ β .* data(R))
  end
  return R
end

function _contract_scalar!(
  R::CuDenseTensor{ElR,NR},
  labelsR,
  T₁::CuDenseTensor,
  labelsT₁,
  T₂::CuDenseTensor,
  labelsT₂,
  α=one(ElR),
  β=zero(ElR),
) where {ElR,NR}
  if nnz(T₁) == nnz(T₂) == 1
    new_R = Tensor(Dense(data(store(T₁)) .* data(store(T₂))), inds(R))
    copyto!(store(R), store(new_R))
  elseif nnz(T₁) == 1
    props = ContractionProperties(labelsT₁, labelsT₂, labelsR)
    compute_contraction_properties!(props, T₁, T₂, R)
    R = _contract!(R, T₁, T₂, props, α, β)
    #perm = getperm(labelsR,labelsT₂)
    #newT2 = Tensor(Dense(data(store(T₁)).*data(store(T₂))), inds(T₂))
    #permute!(R,newT2)
  elseif nnz(T₂) == 1
    props = ContractionProperties(labelsT₁, labelsT₂, labelsR)
    compute_contraction_properties!(props, T₁, T₂, R)
    R = _contract!(R, T₁, T₂, props, α, β)
    #perm = getperm(labelsR,labelsT₁)
    #newT1 = Tensor(Dense(data(store(T₁)).*data(store(T₂))), inds(T₁))
    #permute!(R,newT1)
  else
    error("In _contract_scalar!, one tensor must be a scalar")
  end
  return R
end

function _gemm_contract!(
  CT::DenseTensor{El,NC},
  AT::DenseTensor{El,NA},
  BT::DenseTensor{El,NB},
  props::ContractionProperties,
  α::Number=one(El),
  β::Number=zero(El),
) where {El,NC,NA,NB}
  # TODO: directly use Tensor instead of Array
  C = array(CT)
  A = array(AT)
  B = array(BT)

  tA = 'N'
  if props.permuteA
    pA = NTuple{NA,Int}(props.PA)
    Ap = permutedims(A, pA)
    AM = reshape(Ap, props.dmid, props.dleft)
    tA = 'T'
  else
    #A doesn't have to be permuted
    if Atrans(props)
      AM = reshape(A, props.dmid, props.dleft)
      tA = 'T'
    else
      AM = reshape(A, props.dleft, props.dmid)
    end
  end

  tB = 'N'
  if props.permuteB
    pB = NTuple{NB,Int}(props.PB)
    Bp = permutedims(B, pB)
    BM = reshape(Bp, props.dmid, props.dright)
  else
    if Btrans(props)
      BM = reshape(B, props.dright, props.dmid)
      tB = 'T'
    else
      BM = reshape(B, props.dmid, props.dright)
    end
  end

  #TODO: this logic may be wrong
  if props.permuteC
    #Need to copy here since we will be permuting
    #into C later
    CM = reshape(copy(C), props.dleft, props.dright)
  else
    if Ctrans(props)
      CM = reshape(C, props.dleft, props.dright)
      (AM, BM) = (BM, AM)
      if tA == tB
        tA = tB = (tA == 'T' ? 'N' : 'T')
      end
    else
      CM = reshape(C, props.dleft, props.dright)
    end
  end

  CM = CUBLAS.gemm!(tA, tB, El(α), AM, BM, El(β), CM)

  if props.permuteC
    pC = NTuple{NC,Int}(props.PC)
    Cr = reshape(CM, props.newCrange...)
    @strided C .= permutedims(Cr, pC)
  end
  return C
end

function _contract!(
  CT::CuDenseTensor{El,NC},
  AT::CuDenseTensor{El,NA},
  BT::CuDenseTensor{El,NB},
  props::ContractionProperties,
  α::Number=one(El),
  β::Number=zero(El),
) where {El,NC,NA,NB}
  if ndims(CT) > 12 || ndims(BT) > 12 || ndims(AT) > 12
    return _gemm_contract!(CT, AT, BT, props, α, β)
  end
  Ainds = inds(AT)
  Adims = dims(Ainds)
  Binds = inds(BT)
  Bdims = dims(Binds)
  Cinds = inds(CT)
  Cdims = dims(Cinds)
  Adata = reshape(data(store(AT)), Adims...)
  Bdata = reshape(data(store(BT)), Bdims...)
  Cdata = reshape(data(store(CT)), Cdims...)
  contracted = commoninds(Ainds, Binds)
  A_only = uniqueinds(Ainds, Binds)
  B_only = uniqueinds(Binds, Ainds)
  ind_dict = Vector{Index}()
  for (idx, i) in enumerate(contracted)
    push!(ind_dict, i)
  end
  if length(A_only) > 0
    for (idx, i) in enumerate(A_only)
      push!(ind_dict, i)
    end
  end
  if length(B_only) > 0
    for (idx, i) in enumerate(B_only)
      push!(ind_dict, i)
    end
  end
  ctainds = zeros(Int, length(Ainds))
  ctbinds = zeros(Int, length(Binds))
  ctcinds = zeros(Int, length(Cinds))
  for (ii, ia) in enumerate(Ainds)
    ctainds[ii] = findfirst(x -> x == ia, ind_dict)
  end
  for (ii, ib) in enumerate(Binds)
    ctbinds[ii] = findfirst(x -> x == ib, ind_dict)
  end
  for (ii, ic) in enumerate(Cinds)
    ctcinds[ii] = findfirst(x -> x == ic, ind_dict)
  end
  id_op = CUDA.CUTENSOR.CUTENSOR.CUTENSOR_OP_IDENTITY
  dict_key = ""
  for cc in zip(ctcinds, Cdims)
    dict_key *= string(cc[1]) * "," * string(cc[2]) * ","
  end
  for aa in zip(ctainds, Adims)
    dict_key *= string(aa[1]) * "," * string(aa[2]) * ","
  end
  for bb in zip(ctbinds, Bdims)
    dict_key *= string(bb[1]) * "," * string(bb[2]) * ","
  end
  if haskey(ENV, "CUTENSOR_AUTOTUNE") && tryparse(Int, ENV["CUTENSOR_AUTOTUNE"]) == 1
    if haskey(ContractionPlans, dict_key)
      dict_val = ContractionPlans[dict_key]
      algo = dict_val
      #plan  = dict_val[2]
      Cdata = CUDA.CUTENSOR.contraction!(
        α,
        Adata,
        Vector{Char}(ctainds),
        id_op,
        Bdata,
        Vector{Char}(ctbinds),
        id_op,
        β,
        Cdata,
        Vector{Char}(ctcinds),
        id_op,
        id_op;
        algo=algo,
      )
    else
      # loop through all algos
      # pick the fastest one
      # store that plan!
      best_time = 1e6
      best_plan = nothing
      best_algo = nothing
      max_algos = Ref{Int32}(C_NULL)
      CUDA.CUTENSOR.cutensorContractionMaxAlgos(max_algos)
      # fix once the other options are documented
      #algos = collect(Cint(CUDA.CUTENSOR.CUTENSOR_ALGO_GETT):Cint(max_algos[] - 1))
      algos = collect(Cint(CUDA.CUTENSOR.CUTENSOR_ALGO_GETT):Cint(-1))
      for algo in reverse(algos)
        try
          Cdata, this_time, bytes, gctime, memallocs = @timed CUDA.CUTENSOR.contraction!(
            α,
            Adata,
            Vector{Char}(ctainds),
            id_op,
            Bdata,
            Vector{Char}(ctbinds),
            id_op,
            β,
            Cdata,
            Vector{Char}(ctcinds),
            id_op,
            id_op;
            algo=CUDA.CUTENSOR.cutensorAlgo_t(algo),
          )
          if this_time < best_time
            best_time = this_time
            #best_plan = this_plan
            best_algo = CUDA.CUTENSOR.cutensorAlgo_t(algo)
          end
        catch err
          @warn "Algorithm $algo not supported"
        end
      end
      ContractionPlans[dict_key] = best_algo
    end
  else
    Cdata = CUDA.CUTENSOR.contraction!(
      α,
      Adata,
      Vector{Char}(ctainds),
      id_op,
      Bdata,
      Vector{Char}(ctbinds),
      id_op,
      β,
      Cdata,
      Vector{Char}(ctcinds),
      id_op,
      id_op,
    )
  end
  return parent(Cdata)
end

function Base.:+(B::CuDenseTensor, A::CuDenseTensor)
  opC = CUTENSOR.CUTENSOR_OP_IDENTITY
  opA = CUTENSOR.CUTENSOR_OP_IDENTITY
  opAC = CUTENSOR.CUTENSOR_OP_ADD
  Ais = inds(A)
  Bis = inds(B)
  ind_dict = Vector{Index}()
  for (idx, i) in enumerate(inds(A))
    push!(ind_dict, i)
  end
  Adata = data(store(A))
  Bdata = data(store(B))
  reshapeBdata = reshape(Bdata, dims(Bis)...)
  reshapeAdata = reshape(Adata, dims(Ais)...)
  ctainds = zeros(Int, length(Ais))
  ctbinds = zeros(Int, length(Bis))
  for (ii, ia) in enumerate(Ais)
    ctainds[ii] = findfirst(x -> x == ia, ind_dict)
  end
  for (ii, ib) in enumerate(Bis)
    ctbinds[ii] = findfirst(x -> x == ib, ind_dict)
  end
  ctcinds = copy(ctbinds)
  C = CUDA.zeros(eltype(Bdata), dims(Bis)...)
  CUTENSOR.elementwiseBinary!(
    one(eltype(Adata)),
    reshapeAdata,
    ctainds,
    opA,
    one(eltype(Bdata)),
    reshapeBdata,
    ctbinds,
    opC,
    C,
    ctcinds,
    opAC,
  )
  copyto!(data(store(B)), vec(C))
  return B
end

function Base.:+(B::CuDense, Bis::IndexSet, A::CuDense, Ais::IndexSet)
  opA = CUTENSOR.CUTENSOR_OP_IDENTITY
  opC = CUTENSOR.CUTENSOR_OP_IDENTITY
  opAC = CUTENSOR.CUTENSOR_OP_ADD
  ind_dict = Vector{Index}()
  for (idx, i) in enumerate(Ais)
    push!(ind_dict, i)
  end
  Adata = data(A)
  Bdata = data(B)
  reshapeBdata = reshape(Bdata, dims(Bis)...)
  reshapeAdata = reshape(Adata, dims(Ais)...)
  ctainds = zeros(Int, length(Ais))
  ctbinds = zeros(Int, length(Bis))
  for (ii, ia) in enumerate(Ais)
    ctainds[ii] = findfirst(x -> x == ia, ind_dict)
  end
  for (ii, ib) in enumerate(Bis)
    ctbinds[ii] = findfirst(x -> x == ib, ind_dict)
  end
  ctcinds = copy(ctbinds)
  C = CUDA.zeros(eltype(Bdata), dims(Bis)...)
  Cis = Bis
  C = CUTENSOR.elementwiseBinary!(
    1, reshapeAdata, ctainds, opA, 1, reshapeBdata, ctbinds, opC, C, ctcinds, opAC
  )
  copyto!(data(B), vec(C))
  return C
end

function Base.:-(B::CuDenseTensor, A::CuDenseTensor)
  opC = CUTENSOR.CUTENSOR_OP_IDENTITY
  opA = CUTENSOR.CUTENSOR_OP_IDENTITY
  opAC = CUTENSOR.CUTENSOR_OP_ADD
  Ais = inds(A)
  Bis = inds(B)
  ind_dict = Vector{Index}()
  for (idx, i) in enumerate(inds(A))
    push!(ind_dict, i)
  end
  Adata = data(store(A))
  Bdata = data(store(B))
  reshapeBdata = reshape(Bdata, dims(Bis)...)
  reshapeAdata = reshape(Adata, dims(Ais)...)
  ctainds = zeros(Int, length(Ais))
  ctbinds = zeros(Int, length(Bis))
  for (ii, ia) in enumerate(Ais)
    ctainds[ii] = findfirst(x -> x == ia, ind_dict)
  end
  for (ii, ib) in enumerate(Bis)
    ctbinds[ii] = findfirst(x -> x == ib, ind_dict)
  end
  ctcinds = copy(ctbinds)
  C = CUDA.zeros(eltype(Bdata), dims(Bis))
  CUTENSOR.elementwiseBinary!(
    one(eltype(Adata)),
    reshapeAdata,
    ctainds,
    opA,
    -one(eltype(Bdata)),
    reshapeBdata,
    ctbinds,
    opC,
    C,
    ctcinds,
    opAC,
  )
  copyto!(data(store(B)), vec(C))
  return B
end

function Base.:-(A::CuDense, Ais::IndexSet, B::CuDense, Bis::IndexSet)
  opA = CUTENSOR.CUTENSOR_OP_IDENTITY
  opC = CUTENSOR.CUTENSOR_OP_IDENTITY
  opAC = CUTENSOR.CUTENSOR_OP_ADD
  ind_dict = Vector{Index}()
  for (idx, i) in enumerate(Ais)
    push!(ind_dict, i)
  end
  Adata = data(A)
  Bdata = data(B)
  reshapeBdata = reshape(Bdata, dims(Bis)...)
  reshapeAdata = reshape(Adata, dims(Ais)...)
  ctainds = zeros(Int, length(Ais))
  ctbinds = zeros(Int, length(Bis))
  for (ii, ia) in enumerate(Ais)
    ctainds[ii] = findfirst(x -> x == ia, ind_dict)
  end
  for (ii, ib) in enumerate(Bis)
    ctbinds[ii] = findfirst(x -> x == ib, ind_dict)
  end
  ctcinds = copy(ctbinds)
  C = CUDA.zeros(eltype(Bdata), dims(Bis)...)
  Cis = Bis
  C = CUTENSOR.elementwiseBinary!(
    one(eltype(Adata)),
    reshapeAdata,
    ctainds,
    opA,
    -one(eltype(Bdata)),
    reshapeBdata,
    ctbinds,
    opC,
    C,
    ctcinds,
    opAC,
  )
  copyto!(data(B), vec(C))
  return C
end

function Base.permute!(B::CuDenseTensor, A::CuDenseTensor)
  Ais = inds(A)
  Bis = inds(B)
  ind_dict = Vector{Index}()
  for (idx, i) in enumerate(Ais)
    push!(ind_dict, i)
  end
  Adata = data(store(A))
  Bdata = data(store(B))
  reshapeBdata = reshape(Bdata, dims(Bis)...)
  reshapeAdata = reshape(Adata, dims(Ais)...)
  if ndims(A) < 12 # use CUTENSOR
    ctainds = zeros(Int, length(Ais))
    ctbinds = zeros(Int, length(Bis))
    for (ii, ia) in enumerate(Ais)
      ctainds[ii] = findfirst(x -> x == ia, ind_dict)
    end
    for (ii, ib) in enumerate(Bis)
      ctbinds[ii] = findfirst(x -> x == ib, ind_dict)
    end
    CUDA.CUTENSOR.permutation!(
      one(eltype(Adata)),
      reshapeAdata,
      Vector{Char}(ctainds),
      reshapeBdata,
      Vector{Char}(ctbinds),
    )
  else # use GPUArrays
    perm = Int[]
    for aix in Ais
      b_pos = findfirst(bix -> bix == aix, Bis)
      push!(perm, b_pos)
    end
    @assert isperm(perm)
    permutedims!(reshapeBdata, reshapeAdata, invperm(perm))
  end
  return vec(reshapeBdata)
end

function Base.permute!(B::CuDense, Bis::IndexSet, A::CuDense, Ais::IndexSet)
  ind_dict = Vector{Index}()
  for (idx, i) in enumerate(Ais)
    push!(ind_dict, i)
  end
  Adata = data(A)
  Bdata = data(B)
  reshapeBdata = reshape(Bdata, dims(Bis)...)
  reshapeAdata = reshape(Adata, dims(Ais)...)
  ctainds = zeros(Int, length(Ais))
  ctbinds = zeros(Int, length(Bis))
  for (ii, ia) in enumerate(Ais)
    ctainds[ii] = findfirst(x -> x == ia, ind_dict)
  end
  for (ii, ib) in enumerate(Bis)
    ctbinds[ii] = findfirst(x -> x == ib, ind_dict)
  end

  CUDA.CUTENSOR.permutation!(
    one(eltype(Adata)),
    reshapeAdata,
    Vector{Char}(ctainds),
    reshapeBdata,
    Vector{Char}(ctbinds),
  )
  return vec(reshapeBdata)
end

Base.:/(A::CuDenseTensor, x::Number) = A * inv(x)


text:
"""
    const Point{D,T} = VectorValue{D,T}

Type representing a point of D dimensions with coordinates of type T.
Fields are evaluated at vectors of `Point` objects.
"""
const Point{D,T} = VectorValue{D,T}

"""
    abstract type Field <: Kernel

Abstract type representing physical fields, bases of fields, and other related objects.
These different cases are distinguished by the return value obtained when evaluating them. E.g.,
a physical field returns a vector of values when evaluated at a vector of points, and a basis of `nf` fields
returns a 2d matrix (`np` x `nf`) when evaluated at a vector of `np` points.

The following functions need to be overloaded:

- [`evaluate_field!(cache,f,x)`](@ref)
- [`field_cache(f,x)`](@ref)

The following functions can be also provided optionally

- [`field_gradient(f)`](@ref)
- [`field_return_type(f,x)`](@ref)

Moreover, if the [`field_gradient(f)`](@ref) is not provided, a default implementation that uses the following
functions will be used.

- [`evaluate_gradient!(cache,f,x)`](@ref)
- [`gradient_cache(f,x)`](@ref)

In order to be able to call `field_gradient` again on the resulting object the following methods have to be
provided

- [`evaluate_hessian!(cache,f,x)`](@ref)
- [`hessian_cache(f,x)`](@ref)

These four methods are only designed to be called by the default implementation of [`field_gradient(f)`](@ref) and thus
cannot be assumed that they are available for an arbitrary field. For this reason, these functions are not
exported. The general way of evaluating a gradient of a field is to
build the gradient with [`field_gradient(f)`](@ref) and evaluating the resulting object. For evaluating
the hessian, use two times `field_gradient`.

The interface can be tested with

- [`test_field`](@ref)

Most of the functionality implemented in terms of this interface relies in duck typing (this is why all functions in the interface
have the word "field").  Thus, it is not strictly needed to work with types
that inherit from `Field`. This is specially useful in order to accommodate
existing types into this framework without the need to implement a wrapper type that inherits from `Field`.
For instance, a default implementation is available for numbers, which behave like "constant" fields, or arrays of numbers, which behave like
"constant" bases of fields.  However, we recommend that new types inherit from `Field`.

"""
abstract type Field <: Kernel end

"""
$(SIGNATURES)

Returns the cache object needed to evaluate field `f` at the vector of points `x`.
"""
function field_cache(f,x)
  @abstractmethod
end

"""
$(SIGNATURES)

Returns an array containing the values of evaluating the field `f` at the vector of points
`x` by (possibly) using the scratch data in the `cache` object.  The returned value is
an array,  for which the length of the first axis is `length(x)`, i.e., the number of points where the field has
been evaluated.E.g.,
a physical field returns a vector of `np` values when evaluated at a vector of `np` points, and a basis of `nf` fields
returns a 2d matrix (`np` x `nf`) when evaluated at a vector of `np` points.

This choice is made

- for performance reasons when integrating fields (i.e., adding contributions at different points) since the added values are closer in memory with this layout.
- In order to simplify operations between field objects. E.g., the result of evaluating a physical field and a basis of `nf`fields at a vector of `np` points (which leads to a vector and a matrix of size `(np,)` and `(np,nf)` respectively)  can be conveniently added with the broadcasted sum `.+` operator.

The `cache` object is computed
with the [`field_cache`](@ref) function.

"""
function evaluate_field!(cache,f,x)
  @abstractmethod
end

# Default return type

"""
$(SIGNATURES)

Computes the type obtained when evaluating field `f` at point `x`.
It returns `typeof(evaluate_field(f,x))` by default.
"""
function field_return_type(f,x)
  typeof(evaluate_field(f,x))
end

# Default gradient and hessian implementation

"""
$(SIGNATURES)

Returns another field that represents the gradient of the given one
"""
function field_gradient(f)
  FieldGrad(f)
end

struct FieldGrad{F} <: Field
  field::F
  FieldGrad(f) = new{typeof(f)}(f)
end

@inline function evaluate_field!(cache,f::FieldGrad,x)
  evaluate_gradient!(cache,f.field,x)
end

@inline function field_cache(f::FieldGrad,x)
  gradient_cache(f.field,x)
end

@inline function field_gradient(f::FieldGrad)
  FieldHessian(f.field)
end

struct FieldHessian{F} <: Field
  field::F
  FieldHessian(f) = new{typeof(f)}(f)
end

@inline function evaluate_field!(cache,f::FieldHessian,x)
  evaluate_hessian!(cache,f.field,x)
end

@inline function field_cache(f::FieldHessian,x)
  hessian_cache(f.field,x)
end

@inline function field_gradient(f::FieldHessian)
  @unreachable "Default implementation of 3rt order derivatives not available"
end

"""
$(SIGNATURES)
"""
function evaluate_gradient!(cache,f,x)
  @abstractmethod
end

"""
$(SIGNATURES)
"""
function gradient_cache(cache,x)
  @abstractmethod
end

"""
$(SIGNATURES)
"""
function evaluate_hessian!(cache,f,x)
  @abstractmethod
end

"""
$(SIGNATURES)
"""
function hessian_cache(cache,x)
  @abstractmethod
end

# Implement kernel interface

function kernel_return_type(f::Field,x)
  field_return_type(f,x)
end

function kernel_cache(f::Field,x)
  field_cache(f,x)
end

@inline function apply_kernel!(cache,f::Field,x)
  evaluate_field!(cache,f,x)
end

# Testers

"""
    test_field(
      f,
      x::AbstractVector{<:Point},
      v::AbstractArray,cmp=(==);
      grad=nothing,
      hessian=nothing)

Function used to test the field interface. `v` is an array containing the expected
result of evaluating the field `f` at the vector of points `x`. The comparison is performed using
the `cmp` function. For fields objects that support the `field_gradient` function, the key-word
argument `grad` can be used. It should contain the result of evaluating `field_gradient(f)` at x.
Idem for `hessian`.
The checks are performed with the `@test` macro.
"""
function test_field(
  f,
  x::AbstractVector{<:Point},
  v::AbstractArray,cmp=(==);
  grad=nothing,
  hessian=nothing)

  w = evaluate_field(f,x)

  np, = size(w)
  @test length(x) == np
  @test cmp(w,v)
  @test typeof(w) == field_return_type(f,x)

  cf = field_cache(f,x)
  r = evaluate_field!(cf,f,x)
  @test cmp(r,v)

  _x = vcat(x,x)
  _v = vcat(v,v)
  _w = evaluate_field!(cf,f,_x)
  @test cmp(_w,_v)

  if isa(f,Field)
    test_kernel(f,(x,),v,cmp)
  end

  if grad != nothing
    g = field_gradient(f)
    test_field(g,x,grad,cmp,grad=hessian)
  end

end

# Some API

"""
    gradient_type(::Type{T},x::Point) where T
"""
function gradient_type(::Type{T},x::Point) where T
  typeof(outer(zero(x),zero(T)))
end

"""
    evaluate(f::Field,x)

Equivalent to

    evaluate_field(f,x)

But only for types that inherit from `Field`. Types that implement
the field interface but not inherit from `Field` (e.g., numbers and arrays of numbers)
cannot use this function. Use `evaluate_field` instead.
"""
function evaluate(f::Field,x)
  evaluate_field(f,x)
end

"""
    evaluate!(cache,f::Field,x)

Equivalent to

    evaluate_field!(cache,f,x)

But only for types that inherit from `Field`. Types that implement
the field interface but not inherit from `Field` (e.g., numbers and arrays of numbers)
cannot use this function. Use `evaluate_field!` instead.
"""
@inline function evaluate!(cache,f::Field,x)
  evaluate_field!(cache,f,x)
end

"""
$(SIGNATURES)

Evaluates the field `f` at the vector of points `x` by creating a temporary cache internally.
Equivalent to

    c = field_cache(f,x)
    evaluate_field!(c,f,x)
"""
function evaluate_field(f,x)
  c = field_cache(f,x)
  evaluate_field!(c,f,x)
end

function gradient end

"""
    const ∇ = gradient

Alias for the `gradient` function.
"""
const ∇ = gradient

"""
    gradient(f::Field)

Equivalent to

    field_gradient(f)

But only for types that inherit from `Field`. Types that implement
the field interface but not inherit from `Field` (e.g., numbers and arrays of numbers)
cannot use this function. Use `field_gradient` instead.
"""
function gradient(f::Field)
  field_gradient(f)
end

"""
    field_return_types(f::Tuple,x) -> Tuple

Computes a tuple with the return types of the fields in the tuple `f` when evaluated at the vector
of points `x`

Equivalent to

    tuple(( field_return_type(fi,x) for fi in f)...)
"""
function field_return_types(f::Tuple,x)
  _field_return_types(x,f...)
end

function _field_return_types(x,a,b...)
  Ta = field_return_type(a,x)
  Tb = field_return_types(b,x)
  (Ta,Tb...)
end

function _field_return_types(x,a)
  Ta = field_return_type(a,x)
  (Ta,)
end

"""
    field_caches(f::Tuple,x) -> Tuple

Equivalent to

    tuple((field_cache(fi,x) for fi in f)...)
"""
function field_caches(f::Tuple,x)
  _field_caches(x,f...)
end

function _field_caches(x,a,b...)
  ca = field_cache(a,x)
  cb = field_caches(b,x)
  (ca,cb...)
end

function _field_caches(x,a)
  ca = field_cache(a,x)
  (ca,)
end

"""
    evaluate_fields(f::Tuple,x) -> Tuple

Equivalent to

    tuple((evaluate_fields(fi,x) for fi in f)...)
"""
function evaluate_fields(f::Tuple,x)
  cf = field_caches(f,x)
  evaluate_fields!(cf,f,x)
end

"""
    evaluate_fields!(cf::Tuple,f::Tuple,x) -> Tuple

Equivalent to

    tuple((evaluate_fields!(ci,fi,x) for (ci,fi) in zip(c,f))...)
"""
@inline function evaluate_fields!(cf::Tuple,f::Tuple,x)
  _evaluate_fields!(cf,x,f...)
end

function _evaluate_fields!(c,x,a,b...)
  ca, cb = _split(c...)
  ax = evaluate_field!(ca,a,x)
  bx = evaluate_fields!(cb,b,x)
  (ax,bx...)
end

function _evaluate_fields!(c,x,a)
  ca, = c
  ax = evaluate_field!(ca,a,x)
  (ax,)
end

@inline function _split(a,b...)
  (a,b)
end

"""
    field_gradients(b...) -> Tuple

Equivalent to

    map(field_gradient,b)
"""
@inline function field_gradients(a,b...)
  ga = field_gradient(a)
  gb = field_gradients(b...)
  (ga,gb...)
end

@inline function field_gradients(a)
  ga = field_gradient(a)
  (ga,)
end

"""
    evaluate_all(f::Tuple,x) -> Tuple

Equivalent to

    tuple((evaluate(fi,x) for fi in f)...)
"""
function evaluate_all(f::Tuple,x)
  _evaluate_all(x,f...)
end

function _evaluate_all(x,a,b...)
  ax = evaluate(a,x)
  bx = evaluate_all(b,x)
  (ax, bx...)
end

function _evaluate_all(x,a)
  ax = evaluate(a,x)
  (ax,)
end

"""
    gradient_all(b...) -> Tuple

Equivalent to

    map(gradient,b)
"""
function gradient_all(a,b...)
  ga = gradient(a)
  gb = gradient_all(b...)
  (ga,gb...)
end

function gradient_all(a)
  ga = gradient(a)
  (ga,)
end


text:
using Distributed
import Nanosoldier, GitHub, Sockets

nodes = Dict(Any => addprocs(1))
@everywhere import Nanosoldier

auth = GitHub.authenticate("GITHUB_AUTH00000000000000000000000000000")
secret = GITHUB_SECRET
user = GITHUB_USER
port = 0xffff

config = Nanosoldier.Config(user, nodes, auth, secret;
                            workdir = joinpath(homedir(), "workdir"),
                            trackrepo = "JuliaLang/julia",
                            reportrepo = "JuliaCI/NanosoldierReports",
                            testmode = false)

server = Nanosoldier.Server(config)
run(server, Sockets.IPv4(0), port)


text:
module TetGen
using DocStringExtensions
using TetGen_jll

using GeometryBasics
using GeometryBasics: Polygon, MultiPolygon, Point, LineFace, Polytope, Line,
    Simplex, connect, Triangle, NSimplex, Tetrahedron,
    TupleView, TriangleFace, SimplexFace, LineString, Mesh, TetrahedronP, TriangleP,
    NgonFace, Ngon, faces, coordinates, metafree, meta, faces

using Printf

using StaticArrays

include("cpptetgenio.jl")
include("jltetgenio.jl")
include("rawtetgenio.jl")
include("meshes.jl")
include("api.jl")

export tetrahedralize
export tetunsuitable
export TetGenError
export RawTetGenIO, facetlist!, RawFacet
export numberofpoints,numberoftetrahedra,numberoftrifaces,numberofedges
export volumemesh,surfacemesh

end # module


text:
import java.awt.BorderLayout;
import java.awt.Cursor;
import java.awt.event.ActionEvent;
import java.awt.event.ActionListener;
import java.awt.event.KeyAdapter;
import java.awt.event.KeyEvent;
import java.awt.event.WindowAdapter;
import java.awt.event.WindowEvent;
import java.beans.PropertyChangeEvent;
import java.beans.PropertyChangeListener;
import java.net.MalformedURLException;
import java.net.URL;
import java.util.ArrayList;
import java.io.IOException;

import javax.swing.JButton;
import javax.swing.JEditorPane;
import javax.swing.JFrame;
import javax.swing.JMenu;
import javax.swing.JMenuBar;
import javax.swing.JMenuItem;
import javax.swing.JOptionPane;
import javax.swing.JPanel;
import javax.swing.JScrollPane;
import javax.swing.JTextField;
import javax.swing.event.HyperlinkEvent;
import javax.swing.event.HyperlinkListener;
import javax.swing.text.html.HTMLDocument;
import javax.swing.text.html.HTMLFrameHyperlinkEvent;

// The Simple Web Browser.
public class Browser extends JFrame implements HyperlinkListener {
    private static final long serialVersionUID = 6482307664780440525L;

    // These are the buttons for iterating through the page list.
    private JButton backButton, forwardButton;

    // Page location text field.
    private JTextField locationTextField;

    // Editor pane for displaying pages.
    private JEditorPane displayEditorPane;

    // Browser's list of pages that have been visited.
    private ArrayList pageList = new ArrayList();

    // Constructor for Mini Web Browser.
    public Browser() {
        // Set application title.
        super("Browser");

        // Set window size.
        setSize(1024, 768);

        // Handle closing events.
        addWindowListener(new WindowAdapter() {
            public void windowClosing(WindowEvent e) {
                actionExit();
            }
        });

        // Set up file menu.
        JMenuBar menuBar = new JMenuBar();
        JMenu fileMenu = new JMenu("File");
        fileMenu.setMnemonic(KeyEvent.VK_F);
        JMenuItem fileExitMenuItem = new JMenuItem("Exit", KeyEvent.VK_X);
        fileExitMenuItem.addActionListener(new ActionListener() {
            public void actionPerformed(ActionEvent e) {
                actionExit();
            }
        });
        fileMenu.add(fileExitMenuItem);
        menuBar.add(fileMenu);
        setJMenuBar(menuBar);

        // Set up button panel.
        JPanel buttonPanel = new JPanel();
        backButton = new JButton("< Back");
        backButton.addActionListener(new ActionListener() {
            public void actionPerformed(ActionEvent e) {
                actionBack();
            }
        });
        backButton.setEnabled(false);
        buttonPanel.add(backButton);
        forwardButton = new JButton("Forward >");
        forwardButton.addActionListener(new ActionListener() {
            public void actionPerformed(ActionEvent e) {
                actionForward();
            }
        });
        forwardButton.setEnabled(false);
        buttonPanel.add(forwardButton);
        locationTextField = new JTextField(35);
        locationTextField.addKeyListener(new KeyAdapter() {
            public void keyReleased(KeyEvent e) {
                if (e.getKeyCode() == KeyEvent.VK_ENTER) {
                    actionGo();
                }
            }
        });
        buttonPanel.add(locationTextField);
        JButton goButton = new JButton("GO");
        goButton.addActionListener(new ActionListener() {
            public void actionPerformed(ActionEvent e) {
                actionGo();
            }
        });
        buttonPanel.add(goButton);

        // Set up page display.
	try {
	    displayEditorPane = new CachedJEditorPane();
	} catch (Exception e) {
	    System.err.println("Initializing CachedJEditorPane failed, falling back on JEditorPane...\n");
	    e.printStackTrace();
	    displayEditorPane = new JEditorPane();
	}
        displayEditorPane.setContentType("text/html");
        displayEditorPane.setEditable(false);
        displayEditorPane.addHyperlinkListener(this);
        displayEditorPane.addPropertyChangeListener("page",
                                                    new PropertyChangeListener() {
                                                        public void propertyChange(
                                                                PropertyChangeEvent evt) {
                                                            actionUpdate();
                                                        }
                                                    });

        getContentPane().setLayout(new BorderLayout());
        getContentPane().add(buttonPanel, BorderLayout.NORTH);
        getContentPane().add(new JScrollPane(displayEditorPane),
                             BorderLayout.CENTER);
    }

    // Exit this program.
    private void actionExit() {
        System.exit(0);
    }

    // Go back to the page viewed before the current page.
    private void actionBack() {
        URL currentUrl = displayEditorPane.getPage();
        int pageIndex = pageList.indexOf(currentUrl);
        try {
            showPage((URL) pageList.get(pageIndex - 1));
        }
        catch (Exception e) {
        }
    }

    // Go forward to the page viewed after the current page.
    private void actionForward() {
        URL currentUrl = displayEditorPane.getPage();
        int pageIndex = pageList.indexOf(currentUrl);
        try {
            showPage((URL) pageList.get(pageIndex + 1));
        }
        catch (Exception e) {
        }
    }

    // Load and show the page specified in the location text field.
    private void actionGo() {
        URL verifiedUrl = verifyUrl(locationTextField.getText());
        if (verifiedUrl != null) {
            showPage(verifiedUrl);
        }
        else {
            showError("Invalid URL");
        }
    }

    // Update the history
    private void actionUpdate() {
        // Get URL of new page being displayed.
        URL newUrl = displayEditorPane.getPage();

        int listSize = pageList.size();
        while (listSize > 20) {
            pageList.remove(0);
        }
        pageList.add(newUrl);

        // Update location text field with URL of current page.
        locationTextField.setText(newUrl.toString());

        // Update buttons based on the page being displayed.
        updateButtons();
    }

    // Show dialog box with error message.
    private void showError(String errorMessage) {
        JOptionPane.showMessageDialog(this,
                                      errorMessage,
                                      "Error",
                                      JOptionPane.ERROR_MESSAGE);
    }

    // Verify URL format.
    private URL verifyUrl(String url) {
        // Only allow HTTP URLs.
        if (!url.toLowerCase().startsWith("http://")) return null;

        // Verify format of URL.
        URL verifiedUrl = null;
        try {
            verifiedUrl = new URL(url);
        }
        catch (Exception e) {
            return null;
        }

        return verifiedUrl;
    }

    /* Show the specified page and add it to
       the page list if specified. */
    private void showPage(URL pageUrl) {
        // Show hour glass cursor while crawling is under way.
        setCursor(Cursor.getPredefinedCursor(Cursor.WAIT_CURSOR));

        try {
            // Get URL of page currently being displayed.
            URL currentUrl = displayEditorPane.getPage();

            // Load and display specified page.
            displayEditorPane.setPage(pageUrl);
        }
        catch (Exception e) {
            // Show error messsage.
            showError("Unable to load page");
        }
        finally {
            // Return to default cursor.
            setCursor(Cursor.getDefaultCursor());
        }
    }

    /* Update back and forward buttons based on
       the page being displayed. */
    private void updateButtons() {
        if (pageList.size() < 2) {
            backButton.setEnabled(false);
            forwardButton.setEnabled(false);
        }
        else {
            URL currentUrl = displayEditorPane.getPage();
            int pageIndex = pageList.indexOf(currentUrl);
            backButton.setEnabled(pageIndex > 0);
            forwardButton.setEnabled(pageIndex < (pageList.size() - 1));
        }
    }

    // Handle hyperlink's being clicked.
    public void hyperlinkUpdate(HyperlinkEvent event) {
        HyperlinkEvent.EventType eventType = event.getEventType();
        if (eventType == HyperlinkEvent.EventType.ACTIVATED) {
            if (event instanceof HTMLFrameHyperlinkEvent) {
                HTMLFrameHyperlinkEvent linkEvent =
                        (HTMLFrameHyperlinkEvent) event;
                HTMLDocument document =
                        (HTMLDocument) displayEditorPane.getDocument();
                document.processHTMLFrameHyperlinkEvent(linkEvent);
            }
            else {
                showPage(event.getURL());
            }
        }
    }

    // Run the Mini Browser.
    public static void main(String[] args) throws MalformedURLException {
        Browser browser = new Browser();
        browser.setVisible(true);

	HyperlinkEvent he = new HyperlinkEvent(browser, HyperlinkEvent.EventType.ACTIVATED, new URL("http://en.wikipedia.org/wiki/Main_Page"));
	browser.hyperlinkUpdate(he);
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
@recipe(RidgePlot, x) do scene
    return Attributes(;
                      conf_level=0.95,
                      vline_at_zero=true,
                      show_intercept=true,
                      )
end

"""
    ridgeplot(x::MixedModelBootstrap;
              conf_level=0.95, vline_at_zero=true, show_intercept=true,
              attributes...)

Create a ridge plot for the bootstrap samples of the fixed effects.

Densities are normalized so that the maximum density is always 1.

The highest density interval correspoding to `conf_level` is marked with a bar at the bottom of each density.
Setting `conf_level=missing` removes the markings for the highest density interval.

!!! note
    This functionality is implemented using [Makie recipes](https://makie.juliaplots.org/v0.15.0/recipes.html)
    and thus there are also additional auto-generated methods for `ridgeplot` and `ridgeplot!` that may be useful
    when constructing more complex figures.
"""
function ridgeplot(
    x::MixedModelBootstrap;
    conf_level=0.95,
    vline_at_zero=true,
    show_intercept=true,
    attributes...,
)
    # need to guarantee a min height of 200
    fig = Figure(; resolution=(640, max(200, 100 * _npreds(x; show_intercept))))
    ax = Axis(fig[1, 1])
    if !ismissing(conf_level)
        pl = coefplot!(ax, x; conf_level, vline_at_zero, show_intercept, attributes...)
    end
    pl = ridgeplot!(ax, x; vline_at_zero, conf_level, show_intercept, attributes...)
    return Makie.FigureAxisPlot(fig, ax, pl)
end

function Makie.plot!(ax::Axis, P::Type{<:RidgePlot}, allattrs::Makie.Attributes, x)
    plot = Makie.plot!(ax.scene, P, allattrs, x)

    if haskey(allattrs, :title)
        ax.title = allattrs.title[]
    end
    if haskey(allattrs, :xlabel)
        ax.xlabel = allattrs.xlabel[]
    else
        lab = if !ismissing(allattrs.conf_level[])
           @sprintf  "Normalized bootstrap density and %g%% confidence interval" (allattrs.conf_level[] * 100)
        else
            "Normalized bootstrap density"
        end
        ax.xlabel = lab
    end
    if haskey(allattrs, :ylabel)
        ax.ylabel = allattrs.ylabel[]
    end
    reset_limits!(ax)
    show_intercept = allattrs.show_intercept[]

    # check conf_level so that we don't double print
    # if coefplot took care of it for us
    if ismissing(allattrs.conf_level[])
        cn = _coefnames(x; show_intercept)
        nticks = _npreds(x; show_intercept)
        ax.yticks = (nticks:-1:1, cn)
        ylims!(ax, 0, nticks + 1)
        allattrs.vline_at_zero[] && vlines!(ax, 0; color=(:black, 0.75), linestyle=:dash)
    end
    return plot
end

function Makie.plot!(plot::RidgePlot{<:Tuple{MixedModelBootstrap}})
    boot, conf_level = plot[1][], plot.conf_level[]
    df = transform!(DataFrame(boot.β), :coefname => ByRow(string) => :coefname)
    plot.show_intercept[] || filter!(:coefname => !=("(Intercept)"), df)
    group = :coefname
    densvar = :β
    gdf = groupby(df, group)

    y = length(gdf):-1:1

    dens = combine(gdf, densvar => kde => :kde)

    for (offset, row) in enumerate(reverse(eachrow(dens)))
        # multiply by 0.95 so that the ridges don't overlap
        dd = 0.95 * row.kde.density ./ maximum(row.kde.density)
        lower = Node(Point2f.(row.kde.x, offset))
        upper = Node(Point2f.(row.kde.x, dd .+ offset))
        band!(plot, lower, upper; color=(:black, 0.3))
        lines!(plot, upper; color=(:black, 1.0))
    end

    return plot
end

# """
#     ridgeplot!(ax::Axis, df::AbstractDataFrame, densvar::Symbol, group::Symbol; normalize=false)
#     ridgeplot!(f::Figure, args...; pos=(1,1) kwargs...)
#     ridgeplot(args...; kwargs...)

# Create a "ridge plot".

# A ridge plot is stacked plot of densities for a given variable (`densvar`) grouped by a different variable (`group`). Because densities can very widely in scale, it is sometimes useful to `normalize` the densities so that each density has a maximum of 1.

# The non-mutating method creates a Figure before calling the method for Figure.
# The method for Figure places the ridge plot in the grid position specified by `pos`, default is (1,1).
# """
# function ridgeplot!(ax::Axis, df::AbstractDataFrame, densvar::Symbol, group::Symbol; sort_by_group=false, vline_at_zero=true, normalize=false)
#     # `normalize` makes it so that the max density is always 1
#     # `normalize` works on the density not the area/mass
#     gdf = groupby(df, group)
#     dens = combine(gdf, densvar => kde => :kde)
#     sort_by_group && sort!(dens, group)
#     spacing = normalize ? 1.0 :  0.9 * maximum(dens[!, :kde]) do val
#         return maximum(val.density)
#     end

#     nticks = length(gdf)

#     for (idx, row) in enumerate(eachrow(dens))
#         dd = normalize ? row.kde.density ./ maximum(row.kde.density) : row.kde.density

#         offset =  idx * spacing

#         lower = Node(Point2f.(row.kde.x, offset))
#         upper = Node(Point2f.(row.kde.x, dd .+ offset))
#         band!(ax, lower, upper; color=(:black, 0.3))
#         lines!(ax, upper; color=(:black, 1.0))
#     end

#     vline_at_zero && vlines!(ax, 0; color=(:black, 0.75), linestyle=:dash)

#     ax.yticks[] = (1:spacing:(nticks*spacing), string.(dens[!, group]))
#     ylims!(ax, 0, (nticks + 2) * spacing)
#     ax.xlabel[] = string(densvar)
#     ax.ylabel[] = string(group)

#     ax
# end

# function ridgeplot!(f::Figure, args...; pos=(1,1), kwargs...)
#     ridgeplot!(Axis(f[pos...]), args...; kwargs...)
#     return f
# end

# """
#     ridgeplot(args...; kwargs...)

# See [ridgeplot!](@ref).
# """
# ridgeplot(args...; kwargs...) = ridgeplot!(Figure(), args...; kwargs...)

# """
#     ridgeplot!(Union{Figure, Axis}, bstrp::MixedModelBootstrap, args...; show_intercept=true, kwargs...)
#     ridgeplot(bstrp::MixedModelBootstrap, args...; show_intercept=true, kwargs...)

# Convenience methods that call `DataFrame(bstrp.β)` and pass that onto the primary `ridgeplot[!]` methods.
# By default, the intercept is shown, but this can be disabled.
# """
# function ridgeplot!(axis::Axis, bstrp::MixedModelBootstrap, args...;
#                     show_intercept=true,  normalize=true, kwargs...)

#     df = DataFrame(bstrp.β)
#     show_intercept || filter!(:coefname => !=(Symbol("(Intercept)")), df)
#     ridgeplot!(axis, df, :β, :coefname, args...; sort_by_group=false, normalize, kwargs...)
# end


text:
"""
    struct Sparsity.Ideal{S <: Sparsity.Pattern, C <: AbstractIdealCertificate} <: SumOfSquares.Certificate.AbstractIdealCertificate
        sparsity::S
        certificate::C
    end

Same certificate as `certificate` except that the Sum-of-Squares polynomial `σ`
is modelled as a sum of Sum-of-Squares polynomials with smaller bases
using the sparsity reduction `sparsity`.
"""
struct Ideal{S <: Pattern, C <: SumOfSquares.Certificate.AbstractIdealCertificate} <: SumOfSquares.Certificate.AbstractIdealCertificate
    sparsity::S
    certificate::C
end

function Ideal(sp::Variable, basis, cone, maxdegree::Nothing, newton_polytope)
    error("`maxdegree` cannot be `nothing` when `sparsity` is `Sparsity.Variable`.")
end
function Ideal(sp::Variable, basis, cone, maxdegree::Integer, newton_polytope)
    return Ideal(sp, SumOfSquares.Certificate.MaxDegree(cone, basis, maxdegree))
end
function Ideal(sp::Union{Monomial, SignSymmetry}, basis, cone, maxdegree, newton_polytope)
    return Ideal(sp, SumOfSquares.Certificate.Newton(cone, basis, newton_polytope))
end

function sparsity(poly::MP.AbstractPolynomial, ::Variable, certificate::SumOfSquares.Certificate.MaxDegree)
    H, cliques = chordal_csp_graph(poly, SemialgebraicSets.FullSpace())
    return map(cliques) do clique
        return SumOfSquares.Certificate.maxdegree_gram_basis(certificate.basis, clique, certificate.maxdegree)
    end
end
function sparsity(monos, sp::Union{SignSymmetry, Monomial}, gram_basis::MB.MonomialBasis)
    return MB.MonomialBasis.(sparsity(monos, sp, gram_basis.monomials))
end
function sparsity(poly::MP.AbstractPolynomial, sp::Union{SignSymmetry, Monomial}, certificate::SumOfSquares.Certificate.AbstractIdealCertificate)
    return sparsity(MP.monomials(poly), sp, SumOfSquares.Certificate.gram_basis(certificate, poly))
end
function sparsity(v::SumOfSquares.Certificate.WithVariables, sp, certificate)
    return sparsity(v.inner, sp, certificate)
end
function SumOfSquares.Certificate.gram_basis(certificate::Ideal, poly)
    return sparsity(poly, certificate.sparsity, certificate.certificate)
end
function SumOfSquares.Certificate.gram_basis_type(::Type{Ideal{S, C}}) where {S, C}
    return Vector{<:SumOfSquares.Certificate.gram_basis_type(C)}
end
function SumOfSquares.Certificate.reduced_polynomial(certificate::Ideal, poly, domain)
    return SumOfSquares.Certificate.reduced_polynomial(certificate.certificate, poly, domain)
end
function SumOfSquares.Certificate.cone(certificate::Ideal)
    return SumOfSquares.Certificate.cone(certificate.certificate)
end
SumOfSquares.matrix_cone_type(::Type{Ideal{S, C}}) where {S, C} = SumOfSquares.matrix_cone_type(C)

SumOfSquares.Certificate.zero_basis(certificate::Ideal) = SumOfSquares.Certificate.zero_basis(certificate.certificate)
SumOfSquares.Certificate.zero_basis_type(::Type{Ideal{S, C}}) where {S, C} = SumOfSquares.Certificate.zero_basis_type(C)

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
"""
This functions gives the state inside the Halo orbit given by its
position inside the orbit (first argument), the orbit in the family
(second argument) and the family (third argument)

```julia
intH(x::Real,orb::Real,family::Int)
```

The first and second argument must be between 0 and 1,
and the third argument a 1 or a 2 (H1 or H2), this function
gives back an array of 6 elements, the position and velocity.
"""
function intH(x::Real,orb::Real,family::Int)

    if 0 > x || 1 < x
        error("first argument 'x' must be 0<=x<=1")
    elseif 0 > orb || 1 < orb
        error("second argument 'orb' must be 0<=orb<=1")
    end

    if family == 1 # H1
        i = 1 + x*(size(ORB1,1)-1)
        j = 1 + orb*(size(ORB1,3)-1)
        return interpH1(i,2:7,j)
    elseif family == 2 # H2
        i = 1 + x*(size(ORB2,1)-1)
        j = 1 + orb*(size(ORB2,3)-1)
        return interpH2(i,2:7,j)
    else
        error("last argument has to be either 1 (H1 family) or 2 (H2 family)")
    end

end


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

using Markdown
using InteractiveUtils

# ╔═╡ c1b0a650-3c72-11eb-1fc1-113073ca4477
MLG=["Machine","Learning","Geek"]

# ╔═╡ 7df7a79e-3c73-11eb-3e1c-275b2aaf2f78
for x ∈ MLG
	print(length(x))
end

# ╔═╡ 8b9db1e0-3c75-11eb-095d-471ff4502539
map(MLG) do x length(x) end

# ╔═╡ 4d4c8500-3c76-11eb-3a7b-d711492b441e
map(x->length(MLG[x]),eachindex(MLG))

# ╔═╡ ea18bf70-3c76-11eb-1b50-6192481023f8
[length(x) for x∈MLG]

# ╔═╡ Cell order:
# ╟─c1b0a650-3c72-11eb-1fc1-113073ca4477
# ╠═7df7a79e-3c73-11eb-3e1c-275b2aaf2f78
# ╠═8b9db1e0-3c75-11eb-095d-471ff4502539
# ╠═4d4c8500-3c76-11eb-3a7b-d711492b441e
# ╠═ea18bf70-3c76-11eb-1b50-6192481023f8

max_stars_repo_name: punkmath/Misc

max_issues_repo_name: punkmath/Misc

max_forks_repo_name: punkmath/Misc


text:
function _isavailable_wlclipboard()
    return success(pipeline(`which wl-copy`))
end

"""
    _wlclipboard() -> AbstractMatrix{<:Colorant}

Paste an image from clipboard using `wl-paste`
"""
function _wlclipboard()
    img_buf = IOBuffer()

    # Pipe clipboard image to buffer
    if success(pipeline(`wl-paste -t image/png`; stdout=img_buf))
        # Load image from buffer
        img = load(img_buf)
        return img
    else
        error("No image in clipboard")
    end
end

"""
    _wlclipboard(img::AbstractMatrix{<:Colorant})

Copy an image to clipboard using `wl-copy`
"""
function _wlclipboard(img::AbstractMatrix{<:Colorant})
    img_buf = IOBuffer()

    # Save given image to buffer
    save(Stream{format"PNG"}(img_buf), img)

    # Copy to clipboard
    open(`wl-copy -t image/png`, "w", stdout) do f
        print(f, String(take!(img_buf)))
    end
end


text:
function makedistances(gisregion; scenarioyear="ssp2_2050", res=0.01)
    regions, offshoreregions, regionlist, lonrange, latrange = loadregions(gisregion)
    numreg = length(regionlist)

    println("Finding area- and population-weighted region centers...")
    geocenters, popcenters = getregioncenters(regions, numreg, lonrange, latrange, res, scenarioyear)
    println("\nCalculating distances between centers...")
    distances = [greatcircledistance(Tuple(popcenters[r1,:]), Tuple(popcenters[r2,:])) for r1 = 1:numreg, r2 = 1:numreg]

    println("\nFinding neighboring regions for transmission connections...")
    println("...onshore...")
    connected = connectedregions(regions, numreg)
    println("...offshore...")
    connectedoffshore = connectedregions(offshoreregions, numreg)
    connectedoffshore[connected] .= false
    # regionpop = [sum(pop.==r) for r = 1:numreg]

    println("\nSaving results...")
    matopen(in_datafolder("output", "distances_$gisregion.mat"), "w") do file
        write(file, "distances", distances)
        write(file, "connected", connected)
        write(file, "connectedoffshore", connectedoffshore)
        write(file, "regionlist", string.(regionlist))
        # write(file, "population", pop)
        # write(file, "demand", demand)
        write(file, "regioncenters_lon", popcenters[:,1])
        write(file, "regioncenters_lat", popcenters[:,2])
    end

    # Also save the region file in Matlab format
    # regions2matlab(gisregion)
    nothing
end

# returns great circle distance in km between points given as (lat,lon) tuples (in degrees).
greatcircledistance(point1::Tuple, point2::Tuple) = haversine(point1, point2, 6371.0)

function getregioncenters(regions, numreg, lonrange, latrange, res, scenarioyear)
    lats = (90-res/2:-res:-90+res/2)[latrange]          # latitude values (pixel center)
    cellarea = rastercellarea.(lats, res)

    mkpath(in_datafolder("output"))

    pop = JLD.load(in_datafolder("population_$scenarioyear.jld"), "population")[lonrange,latrange]
    popdens = pop ./ cellarea'

    geocenters, popcenters = regioncenters(regions, numreg, popdens, lonrange, latrange, res)
    return geocenters, popcenters   # column order (lat,lon)
end

function regioncenters(regions, numreg, popdens, lonrange, latrange, res)
    lons = (-180+res/2:res:180-res/2)[lonrange]         # longitude values (pixel center)
    lats = (90-res/2:-res:-90+res/2)[latrange]          # latitude values (pixel center)
    rows, cols = size(regions)
    geocenters = zeros(numreg,2)
    popcenters = zeros(numreg,2)
    counts = zeros(Int, numreg)
    popdenssum = zeros(numreg)
    for c = 1:cols
        lat = lats[c]
        for r = 1:rows
            reg = regions[r,c]
            (reg == 0 || reg == NOREGION) && continue
            lon = lons[r]
            geocenters[reg,:] += [lat, lon]
            popcenters[reg,:] += popdens[r,c] .* [lat, lon]
            counts[reg] += 1
            popdenssum[reg] += popdens[r,c]
        end
    end
    geocenters ./= counts
    popcenters ./= popdenssum
    return geocenters, popcenters   # column order (lat,lon)
end

# find a suitable spot to put the region name label on the map
function find_landarea_near_popcenter(regions, numreg, popcenters, lonrange, latrange, res)
    lons = (-180+res/2:res:180-res/2)[lonrange]         # longitude values (pixel center)
    lats = (90-res/2:-res:-90+res/2)[latrange]          # latitude values (pixel center)
    rows, cols = size(regions)
    filterdistance = max(rows, cols)*res    # about 50 km for a Europe map
    km_per_degree = π*2*6371/360
    disk = diskfilterkernel(filterdistance/km_per_degree/res)
    mostlyland = imfilter((regions.>0) .& (regions.!=NOREGION), disk)     # 0-1
    landcenters = zeros(numreg, 2)
    landfactors = fill(Inf, numreg)
    weight = 0.9
    for c = 1:cols
        lat = lats[c]
        for r = 1:rows
            reg = regions[r,c]
            (reg == 0 || reg == NOREGION) && continue
            lon = lons[r]
            distpop = greatcircledistance((lat,lon), Tuple(popcenters[reg,:]))
            closelandfactor = (1 - weight)*distpop/10000 + weight*(1 - mostlyland[r,c])
            if closelandfactor < landfactors[reg]
                landfactors[reg] = closelandfactor
                landcenters[reg,:] .= lat, lon
            end
        end
    end
    return landcenters
end

function connectedregions(regions, numreg)
    connected = zeros(Bool, numreg, numreg)
    rows, cols = size(regions)
    for c = 1:cols
        for r = 1:rows
            reg = regions[r,c]
            (reg == 0 || reg == NOREGION) && continue
            for adjc = max(1, c-1):min(cols, c+1), adjr = max(1, r-1):min(rows, r+1)
                reg2 = regions[adjr,adjc]
                (reg2 == 0 || reg2 == NOREGION) && continue
                connected[reg,reg2] = true
            end
        end
    end
    for c = 1:numreg, r = 1:numreg
        connected[r,c] = connected[r,c] || connected[c,r]
    end
    return connected
end

max_stars_repo_name: henrikhodel/GlobalEnergyGIS

max_issues_repo_name: henrikhodel/GlobalEnergyGIS

max_forks_repo_name: henrikhodel/GlobalEnergyGIS


text:
function printstyle(node::Node{T,V,S,P,N}) where {T,V,S,P,N<:NodeType}
    s = printnodetype(N)
    s = issmooth(node) ? BLUE_FG(s) : isdecomposable(node) ? RED_FG(s) : s
    s = isdeterministic(node) ? NEGATIVE(s) : s
    return s
end

function prettyprint(io::IO, node::Node{T,V,S,P,SumNode}, level::Int) where {T,V,S,P}
    println(io, string(repeat('\t', level), printstyle(node), " ("))
    for (w, child) in zip(node.params[1:end-1], node.children[1:end-1])
        print(io, string(repeat('\t', level), "$(round(exp(w), digits=3)) × "))
        prettyprint(io, child, level)
        println(io, ", ")
    end
    w = node.params[end]
    print(io, string(repeat('\t', level), "$(round(exp(w), digits=3)) × "))
    prettyprint(io, node.children[end], level)
    print(io, string(repeat('\t', level), ")"))
end

function prettyprint(io::IO, node::Node{T,V,S,P,N}, level::Int) where {T,V,S,P,N}
    println(io, string(repeat('\t', level), printstyle(node), " ("))
    for child in node.children[1:end-1]
        prettyprint(io, child, level+1)
        println(io, ", ")
    end
    prettyprint(io, node.children[end], level+1)
    print(io, string(repeat('\t', level), ")"))
end

function prettyprint(io::IO, node::T, level::Int) where {T<:AbstractLeaf}
    print(io, string(repeat('\t', level), printnodetype(T), "[",params(node),"] - scope: $(node.scope)"))
end

function Base.show(io::IO, node::AbstractNode)
    prettyprint(io, node, 0)
end


text:
module FFTW

using LinearAlgebra, Reexport, Preferences
@reexport using AbstractFFTs
using Base.Threads

import AbstractFFTs: Plan, ScaledPlan,
                     fft, ifft, bfft, fft!, ifft!, bfft!,
                     plan_fft, plan_ifft, plan_bfft, plan_fft!, plan_ifft!, plan_bfft!,
                     rfft, irfft, brfft, plan_rfft, plan_irfft, plan_brfft,
                     fftshift, ifftshift,
                     rfft_output_size, brfft_output_size,
                     plan_inv, normalization

export dct, idct, dct!, idct!, plan_dct, plan_idct, plan_dct!, plan_idct!

include("providers.jl")

function __init__()
    # If someone is trying to set the provider via the old environment variable, warn them that they
    # should instead use `set_provider!()` instead.
    if haskey(ENV, "JULIA_FFTW_PROVIDER")
        Base.depwarn("JULIA_FFTW_PROVIDER is deprecated; use FFTW.set_provider!() instead", :JULIA_FFTW_PROVIDER)
    end

    # Hook FFTW threads up to our partr runtime
    @static if fftw_provider == "fftw"
        fftw_init_threads()
    end
end

# most FFTW calls other than fftw_execute should be protected by a lock to be thread-safe
const fftwlock = ReentrantLock()

# protect thread-safety of expressions or whole functions by fftwlock:
macro exclusive(ex)
    if Meta.isexpr(ex, :function) || (Meta.isexpr(ex, :(=)) && Meta.isexpr(ex.args[1], :call))
        newbody = quote
            lock(fftwlock)
            try
                $(esc(ex.args[2]))
            finally
                unlock(fftwlock)
                destroy_deferred() # deallocate plans queued to be destroyed while we held the lock
            end
        end
        Expr(:function, esc(ex.args[1]), newbody)
    else
        return quote
            lock(fftwlock)
            try
                $(esc(ex))
            finally
                unlock(fftwlock)
                destroy_deferred() # deallocate plans queued to be destroyed while we held the lock
            end
        end
    end
end

include("fft.jl")
include("dct.jl")

include("precompile.jl")
_precompile_()

end # module


text:
@testset "Stress Row/Strings (Equality)" begin
    left = ( "Lorem ipsum dolor sit amet, consectetur adipiscing elit."
             * " Ut bibendum a tortor quis suscipit."
             * " Nulla iaculis malesuada arcu, at blandit purus pellentesque sed."
             * " Etiam interdum feugiat quam ut rutrum."
             * " Mauris at leo leo."
             * " Sed sodales sem sit amet facilisis imperdiet."
             * " Nam laoreet, arcu sed dignissim cursus, ex nulla viverra dui, eget mollis lacus libero a mi."
             * " Curabitur mi felis, tincidunt in tempus vitae, consequat non nisl."
             * " Nulla suscipit lacus sed tincidunt tristique."
             * " Cras nunc orci, luctus non fermentum eget, lobortis ultrices velit."
             * " Vestibulum finibus mi sed varius porttitor."
             * " Praesent id lorem quis massa consectetur laoreet."
             * " Etiam ornare augue eu ligula feugiat, sed ornare neque tempus."
             * " Cras maximus leo ante, nec pharetra libero fringilla sit amet." )

    right = ( "Lorem ipsum dolor sit amet, consectetur adipiscing elit."
              * " Ut bibendum a tortor quis suscipit."
              * " Nulla iaculis malesuada arcu, at blandit purus pellentesque sed."
              * " Etiam interdum feugiat quam ut rutrum."
              * " Mauris at leo leo."
              * " Sed sodales sem sit amet facilisis imperdiet."
              * " Nam laoreet, arcu sed dignissim cursus, ex nulla viverra dui, eget mollis lacus libero a mi."
              * " Curabitur mi felis, tincidunt in tempus vitae, consequat non nisl."
              * " Nulla suscipit lacus sed tincidunt tristique."
              * " Cras nunc orci, luctus non fermentum eget, lobortis ultrices velit."
              * " Vestibulum finibus mi sed varius porttitor."
              * " Praesent id lorem quis massa consectetur laoreet."
              * " Etiam ornare augue eu ligula feugiat, sed ornare neque tempus."
              * " Cras maximus leo ante, nec pharetra libero fringilla sit amet." )

    estimated = (length(right) + 1) * 2 * sizeof(Int64)
    allocated = @allocated distance_row(left, right)
    @test isapprox(estimated, allocated, atol = 0.05 * estimated)
end

@testset "Stress Row/Strings (Inequality)" begin
    left = ( "Lorem ipsum dolor sit amet, consectetur adipiscing elit."
             * " Ut bibendum a tortor quis suscipit."
             * " Nulla iaculis malesuada arcu, at blandit purus pellentesque sed."
             * " Etiam interdum feugiat quam ut rutrum."
             * " Mauris at leo leo."
             * " Sed sodales sem sit amet facilisis imperdiet."
             * " Nam laoreet, arcu sed dignissim cursus, ex nulla viverra dui, eget mollis lacus libero a mi."
             * " Curabitur mi felis, tincidunt in tempus vitae, consequat non nisl."
             * " Nulla suscipit lacus sed tincidunt tristique."
             * " Cras nunc orci, luctus non fermentum eget, lobortis ultrices velit."
             * " Vestibulum finibus mi sed varius porttitor."
             * " Praesent id lorem quis massa consectetur laoreet."
             * " Etiam ornare augue eu ligula feugiat, sed ornare neque tempus."
             * " Cras maximus leo ante, nec pharetra libero fringilla sit amet." )

    right = ( "Sed venenatis, nibh eu hendrerit fermentum, odio sapien imperdiet massa, in sodales lorem mi vitae metus."
              * " Sed nec mi sem."
              * " Cras eu sapien pulvinar, laoreet augue vel, pharetra lacus."
              * " Nam a quam lorem."
              * " Vivamus ornare sem sem, at consectetur mauris tincidunt vitae."
              * " Phasellus eu eros leo."
              * " In scelerisque et ante non varius."
              * " Nullam ante ante, egestas eget lacinia id, commodo quis leo."
              * " Praesent in nunc a odio interdum luctus."
              * " Vestibulum vitae tristique nunc, et tincidunt dolor."
              * " Proin congue lectus ut metus rhoncus, vitae ultricies dolor placerat."
              * " Morbi quam lectus, tincidunt vitae libero condimentum, rhoncus pharetra risus."
              * " Phasellus felis orci, laoreet eget auctor sed, luctus quis nisl."
              * " Vestibulum non rhoncus arcu."
              * " Sed ut vestibulum eros, eget aliquet nibh."
              * " Etiam blandit, orci sit amet hendrerit fringilla, mi purus sagittis velit, et luctus nibh felis at diam." )

    estimated = (length(right) + 1) * 2 * sizeof(Int64)
    allocated = @allocated distance_row(left, right)
    @test isapprox(estimated, allocated, atol = 0.05 * estimated)
end

max_stars_repo_name: jlegare/LevenshteinToolkit

max_issues_repo_name: jlegare/LevenshteinToolkit

max_forks_repo_name: jlegare/LevenshteinToolkit


text:
struct Knnr
    X::Array{Float64}
    Y::Array{Float64}
    fm
    nlvdis::Int
    metric::String
    h::Real
    k::Int
    tol::Real
end

"""
    knnr(X, Y; nlvdis = 0, metric = "eucl", h = Inf, k = 1, tol = 1e-4)
k-Nearest-Neighbours regression (KNNR).
* `X` : X-data.
* `Y` : Y-data.
* `nlvdis` : Number of latent variables (LVs) to consider in the global PLS used for the dimension reduction before
    calculating the dissimilarities. If `nlvdis = 0`, there is no dimension reduction.
* `metric` : Type of dissimilarity used to select the neighbors. Possible values are "eucl" (default; Euclidean distance)
    and "mahal" (Mahalanobis distance).
* `h` : A scalar defining the shape of the weight function. Lower is h, sharper is the function. See function `wdist`.
* `k` : The number of nearest neighbors to select for each observation to predict.
* `tol` : For stabilization when very close neighbors.

For each new observation to predict, the prediction is the weighted mean
over the selected (training) neighborhood. The weights are defined from
the dissimilarities between the new observation and the selected neighborhood,
and are computed from function 'wdist'.

The dissimilarities can be calculated from the original X-data
or after a dimension reduction, depending on argument `nlvdis`.
In the last case, global PLS scores (LVs) are computed from (`X`, `Y`)
and the dissimilarities are computed over these scores.
In general, for high dimensional X-data, using the Mahalanobis distance requires
psuch a reliminary dimensionality reduction.
"""
function knnr(X, Y; nlvdis = 0, metric = "eucl", h = Inf, k = 1, tol = 1e-4)
    X = ensure_mat(X)
    Y = ensure_mat(Y)
    if nlvdis == 0
        fm = nothing
    else
        fm = plskern(X, Y; nlv = nlvdis)
    end
    return Knnr(X, Y, fm, nlvdis, metric, h, k, tol)
end

"""
    predict(object::Knnr, X)
Compute the Y-predictions from the fitted model.
* `object` : The fitted model.
* `X` : X-data for which predictions are computed.
"""
function predict(object::Knnr, X)
    X = ensure_mat(X)
    m = size(X, 1)
    q = size(object.Y, 2)
    # Getknn
    if isnothing(object.fm)
        res = getknn(object.X, X; k = object.k, metric = object.metric)
    else
        res = getknn(object.fm.T, transform(object.fm, X); k = object.k, metric = object.metric)
    end
    listw = copy(res.d)
    for i = 1:m
        w = wdist(res.d[i]; h = object.h)
        w[w .< object.tol] .= object.tol
        listw[i] = w
    end
    # End
    pred = zeros(m, q)
    @inbounds for i = 1:m
        s = res.ind[i]
        w = mweights(listw[i])
        pred[i, :] .= colmeans(vrow(object.Y, s), w)
    end
    (pred = pred, listnn = res.ind, listd = res.d, listw = listw)
end