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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 571613a1-6b4b-496d-9a68-aac3f6a83a4b
using PlutoUI, HypertextLiteral

# ╔═╡ 97914842-76d2-11eb-0c48-a7eedca870fb
md"""
# Using _JavaScript_ inside Pluto

You have already seen that Pluto is designed to be _interactive_. You can make fantastic explorable documents using just the basic inputs provided by PlutoUI, together with the wide range of visualization libraries that Julia offers.

_However_, if you want to take your interactive document one step further, then Pluto offers a great framework for **combining Julia with HTML, CSS and _JavaScript_**.
"""

# ╔═╡ 168e13f7-2ff2-4207-be56-e57755041d36
md"""
## Prerequisites

This document assumes that you have used HTML, CSS and JavaScript before in another context. If you know Julia, and you want to add these web languages to your skill set, we encourage you to do so! It will be useful knowledge, also outside of Pluto.

"""

# ╔═╡ 28ae1424-67dc-4b76-a172-1185cc76cb59
@htl("""

		Learning HTML and CSS

		It is easy to learn HTML and CSS because they are not 'programming languages' like Julia and JavaScript, they are markup languages: there are no loops, functions or arrays, you only declare how your document is structured (HTML) and what that structure looks like on a 2D color display (CSS).

		As an example, this is what this cell looks like, written in HTML and CSS:

""")

# ╔═╡ ea39c63f-7466-4015-a66c-08bd9c716343
md"""
> My personal favourite resource for learning HTML and CSS is the [Mozilla Developer Network (MDN)](https://developer.mozilla.org/en-US/docs/Web/CSS).
>
> _-fons_
"""

# ╔═╡ 8b082f9a-073e-4112-9422-4087850fc89e
md"""
#### Learning JavaScript
After learning HTML and CSS, you can already spice up your Pluto notebooks by creating custom layouts, generated dynamically from Julia data. To take things to the next level, you can learn JavaScript. We recommend using an online resource for this.

> My personal favourite is [javascript.info](https://javascript.info/), a high-quality, open source tutorial. I use it too!
>
> _-fons_

It is hard to say whether it is easy to _learn JavaScript using Pluto_. On one hand, we highly recommend the high-quality public learning material that already exists for JavaScript, which is generally written in the context of writing traditional web apps. On the other hand, if you have a specific Pluto-related project in mind, then this could be a great motivator to continue learning!

A third option is to learn JavaScript using [observablehq.com](https://observablehq.com), an online reactive notebook for JavaScript (it's awesome!). Pluto's JavaScript runtime is designed to be very close to the way you write code in observable, so the skills you learn there will be transferable!

If you chose to learn JavaScript using Pluto, let me know how it went, and how we can improve! [fons@plutojl.org](mailto:fons@plutojl.org)
"""

# ╔═╡ d70a3a02-ef3a-450f-bf5a-4a0d7f6262e2
TableOfContents()

# ╔═╡ 10cf6ed1-8276-4a4a-ad06-097d10335512
md"""
# Essentials

## Using HTML, CSS and JavaScript

To use web languages inside Pluto, we recommend the small package [`HypertextLiteral.jl`](https://github.com/MechanicalRabbit/HypertextLiteral.jl), which provides an `@htl` macro.

You wrap `@htl` around a string expression to mark it as an *HTML literal*, as we did in the example cell from earlier. When a cell outputs an HTML-showable object, it is rendered directly in your browser.
"""

# ╔═╡ d967cdf9-3df9-40bb-9b08-09cae95a5ca7
@htl("  Hello!  ")

# ╔═╡ 858745a9-cd59-43a6-a296-803515518e57
md"""
### CSS and JavaScript

You can use CSS and JavaScript by including it inside HTML, just like you do when writing a web page.

For example, here we use `<script>` to include some JavaScript, and `<style>` to include CSS.
"""

# ╔═╡ 21a9e3e6-92f4-475d-9c8e-21e15c09336b
@htl("""

Hello!

""")

# ╔═╡ 4a3398be-ee86-45f3-ac8b-f627a38c00b8
md"""
## Interpolation

Julia has a nice feature: _string interpolation_:
"""

# ╔═╡ 2d5fd611-284b-4428-b6a5-8909203990b9
who = "🌍"

# ╔═╡ 82de4674-9ecc-46c4-8a57-0b4453c579c3
"Hello $(who)!"

# ╔═╡ 70a415be-881a-4c01-9f8c-635b8b89e1ad
md"""
With some (frustrating) exceptions, you can also interpolate into Markdown literals:
"""

# ╔═╡ 730a692f-2bf2-4d5b-86da-6ab861e8b8ac
md"""
Hello $(who)!
"""

# ╔═╡ a45fdec4-2d4b-429b-b809-4c256b57fffe
md"""
**However**, you cannot interpolate into an `html"` string:
"""

# ╔═╡ c68ebd7b-5fb6-4527-ac34-33f9730e4587
html"""

Hello $(who)!

"""

# ╔═╡ 8c03139f-a94b-40cc-859f-0d86f1c72143
md"""

😢 Luckily we can perform these kinds of interpolations (and much more) with the `@htl` macro, as we will see next.

### Interpolating into HTML -- HypertextLiteral.jl
"""

# ╔═╡ d8dcb044-0ac8-46d1-a043-1073bb6d1ff1
@htl("""

 Hello $(who)!

	""")

# ╔═╡ e7d3db79-8253-4cbd-9832-5afb7dff0abf
cool_features = [
	md"Interpolate any **HTML-showable object**, such as plots and images, or another `@htl` literal."
	md"Interpolated lists are expanded _(like in this cell!)_."
	"Easy syntax for CSS"
	]

# ╔═╡ bf592202-a9a4-4e9b-8433-fed55e3aa3bc
@htl("""

It has a bunch of very cool features! Including:

$([
		@htl(
			"
$(item)
"
		)
		for item in cool_features
	])

	""")

# ╔═╡ 5ac5b984-8c02-4b8d-a342-d0f05f7909ec
md"""
#### Why not just `HTML(...)`?

You might be thinking, why don't we just use the `HTML` function, together with string interpolation? The main problem is correctly handling HTML _escaping rules_. For example:
"""

# ╔═╡ ef28eb8d-ec98-43e5-9012-3338c3b84f1b
cool_elements = "
 and <marquee>"

# ╔═╡ 1ba370cc-3631-47ea-9db5-75587e8e4ff3
HTML("""

 My favourite HTML elements are $(cool_elements)!

""")

# ╔═╡ 7fcf2f3f-d902-4338-adf0-8ef181e79420
@htl("""

 My favourite HTML elements are $(cool_elements)!

""")

# ╔═╡ 7afbf8ef-e91c-45b9-bf22-24201cbb4828
md"""
### Interpolating into JS -- _HypertextLiteral.jl_

As we see above, using HypertextLiteral.jl, we can interpolate objects (numbers, string, images) into HTML output, great! Next, we want to **interpolate _data_ into _scripts_**. Although you could use `JSON.jl`, HypertextLiteral.jl actually has this ability built-in!

> When you **interpolate Julia objects into a `<script>` tag** using the `@htl` macro, it will be converted to a JS object _automatically_.
"""

# ╔═╡ b226da72-9512-4d14-8582-2f7787c25028
simple_data = (msg="Hello! ", times=3)

# ╔═╡ a6fd1f7b-a8fc-420d-a8bb-9f549842ad3e
@htl("""

""")

# ╔═╡ 965f3660-6ec4-4a86-a2a2-c167dbe9315f
md"""
**Let's look at a more exciting example:**
"""

# ╔═╡ 00d97588-d591-4dad-9f7d-223c237deefd
@bind fantastic_x Slider(0:400)

# ╔═╡ 01ce31a9-6856-4ee7-8bce-7ce635167457
my_data = [
	(name="Cool", coordinate=[100, 100]),
	(name="Awesome", coordinate=[200, 100]),
	(name="Fantastic!", coordinate=[fantastic_x, 150]),
]

# ╔═╡ 21f57310-9ceb-423c-a9ce-5beb1060a5a3
@htl("""

""")

# ╔═╡ 7d9d6c28-131a-4b2a-84f8-5c085f387e85
md"""
#### Future: directly embedding data

In the future, you will be able to embed data directly into JavaScript, using Pluto's built-in, optimized data transfer. See [the Pull Request](https://github.com/fonsp/Pluto.jl/pull/1124) for more info.
"""

# ╔═╡ 0866afc2-fd42-42b7-a572-9d824cf8b83b
md"""
## Custom `@bind` output
"""

# ╔═╡ 75e1a973-7ef0-4ac5-b3e2-5edb63577927
md"""
**You can use JavaScript to write input widgets.** The `input` event can be triggered on any object using

```javascript
obj.value = ...
obj.dispatchEvent(new CustomEvent("input"))
```

For example, here is a button widget that will send the number of times it has been clicked as the value:

"""

# ╔═╡ e8d8a60e-489b-467a-b49c-1fa844807751
ClickCounter(text="Click") = @htl("""

<button>$(text)</button>

""")

# ╔═╡ 9346d8e2-9ba0-4475-a21f-11bdd018bc60
@bind num_clicks ClickCounter()

# ╔═╡ 7822fdb7-bee6-40cc-a089-56bb32d77fe6
num_clicks

# ╔═╡ 701de4b8-42d3-46a3-a399-d7761dccd83d
md"""
As an exercise to get familiar with these techniques, you can try the following:
- 👉 Add a "reset to zero" button to the widget above.
- 👉 Make the bound value an array that increases size when you click, instead of a single number.
- 👉 Create a "two sliders" widget: combine two sliders (`<input type=range>`) into a single widget, where the bound value is the two-element array with both values.
- 👉 Create a "click to send" widget: combine a text input and a button, and only send the contents of the text field when the button is clicked, not on every keystroke.

Questions? Ask them on our [GitHub Discussions](https://github.com/fonsp/Pluto.jl/discussions)!
"""

# ╔═╡ 88120468-a43d-4d58-ac04-9cc7c86ca179
md"""
## Debugging

The HTML, CSS and JavaScript that you write run in the browser, so you should use the [browser's built-in developer tools](https://developer.mozilla.org/en-US/docs/Learn/Common_questions/What_are_browser_developer_tools) to debug your code.
"""

# ╔═╡ ea4b2da1-4c83-4a1f-8fc3-c71a120e58e1
@htl("""

""")

# ╔═╡ 08bdeaff-5bfb-49ab-b4cc-3a3446c63edc
@htl("""

Can you find out which CSS class this is?

	""")

# ╔═╡ 9b6b5da9-8372-4ebf-9c66-ae9fcfc45d47
md"""
## Selecting elements

When writing the javascript code for a widget, it is common to **select elements inside the widgets** to manipulate them. In the number-of-clicks example above, we selected the `<span>` and `<button>` elements in our code, to trigger the input event, and attach event listeners, respectively.

There are a numbers of ways to do this, and the recommended strategy is to **create a wrapper `<span>`, and use `currentScript.parentElement` to select it**.

### `currentScript`

When Pluto runs the code inside `<script>` tags, it assigns a reference to that script element to a variable called `currentScript`. You can then use properties like `previousElementSibling` or `parentElement` to "navigate to other elements".

Let's look at the "wrapper span strategy" again.

```htmlmixed
@htl("\""

<!-- the wrapper span -->
<span>

	<button id="first">Hello</button>
	<button id="second">Julians!</button>

	<script>
		const wrapper_span = currentScript.parentElement
		// we can now use querySelector to select anything we want
		const first_button = wrapper_span.querySelector("button#first")

		console.log(first_button)
	</script>
</span>
"\"")
```
"""

# ╔═╡ f18b98f7-1e0f-4273-896f-8a667d15605b
md"""
#### Why not just select on `document.body`?

In the example above, it would have been easier to just select the button directly, using:
```javascript
// ⛔ do no use:
const first_button = document.body.querySelector("button#first")
```

However, this becomes a problem when **combining using the widget multiple times in the same notebook**, since all selectors will point to the first instance.

Similarly, try not to search relative to the `<pluto-cell>` or `<pluto-output>` element, because users might want to combine multiple instances of the widget in a single cell.
"""

# ╔═╡ d83d57e2-4787-4b8d-8669-64ed73d79e73
md"""
## Script loading

To use external javascript dependencies, you can load them from a CDN, such as:
- [jsdelivr.com](https://www.jsdelivr.com/)
- [esm.sh](https://esm.sh)

Just like when writing a browser app, there are two ways to import JS dependencies: a `<script>` tag, and the more modern ES6 import.

### Loading method 1: ES6 imports

We recommend that you use an [**ES6 import**](https://developer.mozilla.org/en-US/docs/Web/JavaScript/Guide/Modules) if the library supports it. (If it does not, you might be able to still get it using [esm.sh](https://esm.sh) or [skypack.dev](https://www.skypack.dev)!)

##### Awkward note about syntax

Normally, you can import libraries inside JS using the import syntax:
```javascript
// ⛔ do no use:
import confetti from "https://esm.sh/canvas-confetti@1.4.0"
import { html, render, useEffect } from "https://cdn.jsdelivr.net/npm/htm@3.0.4/preact/standalone.mjs"
```

In Pluto, this is [currently not yet supported](https://github.com/fonsp/Pluto.jl/issues/992), and you need to use a different syntax as workaround:
```javascript
// ✔ use:
const { default: confetti } = await import("https://esm.sh/canvas-confetti@1.4.0")
const { html, render, useEffect } = await import("https://cdn.jsdelivr.net/npm/htm@3.0.4/preact/standalone.mjs")
```
"""

# ╔═╡ 077c95cf-2a1b-459f-830e-c29c11a2c5cc
md"""

### Loading method 2: script tag

`<script src="...">` tags with a `src` attribute set, like this tag to import the d3.js library:

```html
<script src="https://cdn.jsdelivr.net/npm/d3@6.2.0/dist/d3.min.js"></script>
```

will work as expected. The execution of other script tags within the same cell is delayed until a `src` script finished loading, and Pluto will make sure that every source file is only loaded once.
"""

# ╔═╡ 80511436-e41f-4913-8a30-d9e113cfaf71
md"""
### Pinning versions

When using a CDN almost **never** want to use an unpinned import. Always version your CDN imports!
```js
// ⛔ do no use:
"https://esm.sh/canvas-confetti"
"https://cdn.jsdelivr.net/npm/htm/preact/standalone.mjs"

// ✔ use:
"https://esm.sh/canvas-confetti@1.4.0"
"https://cdn.jsdelivr.net/npm/htm@3.0.4/preact/standalone.mjs"
```
"""

# ╔═╡ 8388a833-d535-4cbd-a27b-de323cea60e8
md"""
# Advanced
"""

# ╔═╡ 4cf27df3-6a69-402e-a71c-26538b2a52e7
md"""
## Script output & `observablehq/stdlib`

Pluto's original inspiration was [observablehq.com](https://observablehq.com/), and online reactive notebook for JavaScript. _(It's REALLY good, try it out!)_ We design Pluto's JavaScript runtime to be close to the way you write code in observable.

Read more about the observable runtime in their (interactive) [documentation](https://observablehq.com/@observablehq/observables-not-javascript). The following is also true for JavaScript-inside-scripts in Pluto:
- ⭐ If you return an HTML node, it will be displayed.
- ⭐ The [`observablehq/stdlib`](https://observablehq.com/@observablehq/stdlib) library is pre-imported, you can use `DOM`, `html`, `Promises`, etc.
- ⭐ When a cell re-runs reactively, `this` will be set to the previous output (with caveat, see the later section)
- The variable `invalidation` is a Promise that will get resolved when the cell output is changed or removed. You can use this to remove event listeners, for example.
- You can use top-level `await`, and a returned HTML node will be displayed when ready.
- Code is run in "strict mode", use `let x = 1` instead of `x = 1`.

The following is different in Pluto:
- JavaScript code is not reactive, there are no global variables.
- Cells can contain multiple script tags, and they will run consecutively (also when using `await`)
- We do not (yet) support async generators, i.e. `yield`.
- We do not support the observable keywords `viewof` and `mutable`.
"""

# ╔═╡ 5721ad33-a51a-4a91-adb2-0915ea0efa13
md"""
### Example:
(Though using `HypertextLiteral.jl` would make more sense for this purpose.)
"""

# ╔═╡ fc8984c8-4668-418a-b258-a1718809470c

# ╔═╡ 846354c8-ba3b-4be7-926c-d3c9cc9add5f
films = [
	(title="Frances Ha", director="Noah Baumbach", year=2012),
	(title="Portrait de la jeune fille en feu", director="Céline Sciamma", year=2019),
	(title="De noorderlingen", director="Alex van Warmerdam", year=1992),
];

# ╔═╡ c857bb4b-4cf4-426e-b340-592cf7700434
@htl("""

	""")

# ╔═╡ a33c7d7a-8071-448e-abd6-4e38b5444a3a
md"""
## Stateful output with `this`

Just like in observablehq, if a cell _re-runs reactively_, then the javascript variable `this` will take the value of the last thing that was returned by the script. If the cell runs for the first time, then `this == undefined`. In particular, if you return an HTML node, and the cell runs a second time, then you can access the HTML node using `this`. Two reasons for using this feature are:
- Stateful output: you can persist some state in-between re-renders.
- Performance: you can 'recycle' the previous DOM and update it partially (using d3, for example). _When doing so, Pluto guarantees that the DOM node will always be visible, without flicker._

##### 'Re-runs reactively'?
With this, we mean that the Julia cell re-runs not because of user input (Ctrl+S, Shift+Enter or clicking the play button), but because it was triggered by a variable reference.

##### ☝ Caveat
This feature is **only enabled** for `<script>` tags with the `id` attribute set, e.g. `<script id="first">`. Think of setting the `id` attribute as saying: "I am a Pluto script". There are two reasons for this:
- One Pluto cell can output multiple scripts, Pluto needs to know which output to assign to which script.
- Some existing scripts assume that `this` is set to `window` in toplevel code (like in the browser). By hiding the `this`-feature behind this caveat, we still support libraries that output such scripts.

"""

# ╔═╡ 91f3dab8-5521-44a0-9890-8d988a994076
trigger = "edit me!"

# ╔═╡ dcaae662-4a4f-4dd3-8763-89ea9eab7d43
let
	trigger

	html"""

	"""
end

# ╔═╡ e77cfefc-429d-49db-8135-f4604f6a9f0b
md"""
### Example: d3.js transitions

Type the coordinates of the circles here!
"""

# ╔═╡ 2d5689f5-1d63-4b8b-a103-da35933ad26e
@bind positions TextField(default="100, 300")

# ╔═╡ 6dd221d1-7fd8-446e-aced-950512ea34bc
dot_positions = try
	parse.([Int], split(replace(positions, ',' => ' ')))
catch e
	[100, 300]
end

# ╔═╡ 0a9d6e2d-3a41-4cd5-9a4e-a9b76ed89fa9
# dot_positions = [100, 300] # edit me!

# ╔═╡ 0962d456-1a76-4b0d-85ff-c9e7dc66621d
md"""
Notice that, even though the cell below re-runs, we **smoothly transition** between states. We use `this` to maintain the d3 transition states in-between reactive runs.
"""

# ╔═╡ bf9b36e8-14c5-477b-a54b-35ba8e415c77
@htl("""

""")

# ╔═╡ 781adedc-2da7-4394-b323-e508d614afae
md"""
### Example: Preact with persistent state
"""

# ╔═╡ de789ad1-8197-48ae-81b2-a21ec2340ae0
md"""
Modify `x`, add and remove elements, and notice that preact maintains its state.
"""

# ╔═╡ 85483b28-341e-4ed6-bb1e-66c33613725e
x = ["hello pluton!", 232000,2,2,12 ,12,2,21,1,2, 120000]

# ╔═╡ 3266f9e6-42ad-4103-8db3-b87d2c315290
state = Dict(
	:x => x
	)

# ╔═╡ 9e37c18c-3ebb-443a-9663-bb4064391d6e
@htl("""

""")

# ╔═╡ da7091f5-8ba2-498b-aa8d-bbf3b4505b81
md"""
# Appendix
"""

# ╔═╡ 64cbf19c-a4e3-4cdb-b4ec-1fbe24be55ad
details(x, summary="Show more") = @htl("""

		$(summary)
		$(x)

	""")

# ╔═╡ 93abe0dc-f041-475f-9ef7-d8ee4408414b
details(md"""
	```htmlmixed

	<article class="learning">
		<h4>
			Learning HTML and CSS
		</h4>
		<p>
			It is easy to learn HTML and CSS because they are not 'programming languages' like Julia and JavaScript, they are <em>markup languages</em>: there are no loops, functions or arrays, you only <em>declare</em> how your document is structured (HTML) and what that structure looks like on a 2D color display (CSS).
		</p>
		<p>
			As an example, this is what this cell looks like, written in HTML and CSS:
		</p>
	</article>

	<style>

		article.learning {
			background: #fde6ea4c;
			padding: 1em;
			border-radius: 5px;
		}

		article.learning h4::before {
			content: "☝";
		}

		article.learning p::first-letter {
			font-size: 1.5em;
			font-family: cursive;
		}

	</style>
	```
	""", "Show with syntax highlighting")

# ╔═╡ d12b98df-8c3f-4620-ba3c-2f3dadac521b
details(md"""
	```htmlmixed
	<script>

	// interpolate the data 🐸
	const data = $(simple_data)

	const span = document.createElement("span")
	span.innerText = data.msg.repeat(data.times)

	return span
	</script>
	```
	""", "Show with syntax highlighting")

# ╔═╡ 94561cb1-2325-49b6-8b22-943923fdd91b
details(md"""
	```htmlmixed
	<script src="https://cdn.jsdelivr.net/npm/d3@6.2.0/dist/d3.min.js"></script>

	<script>

	// interpolate the data 🐸
	const data = $(my_data)

	const svg = DOM.svg(600,200)
	const s = d3.select(svg)

	s.selectAll("text")
		.data(data)
		.join("text")
		.attr("x", d => d.coordinate[0])
		.attr("y", d => d.coordinate[1])
		.style("fill", "red")
		.text(d => d.name)

	return svg
	</script>
	```
	""", "Show with syntax highlighting")

# ╔═╡ b0c246ed-b871-461b-9541-280e49b49136
details(md"""
```htmlmixed
<div>
<button>$(text)</button>

<script>

	// Select elements relative to `currentScript`
	const div = currentScript.parentElement
	const button = div.querySelector("button")

	// we wrapped the button in a `div` to hide its default behaviour from Pluto

	let count = 0

	button.addEventListener("click", (e) => {
		count += 1

		// we dispatch the input event on the div, not the button, because
		// Pluto's `@bind` mechanism listens for events on the **first element** in the
		// HTML output. In our case, that's the div.

		div.value = count
		div.dispatchEvent(new CustomEvent("input"))
		e.preventDefault()
	})

	// Set the initial value
	div.value = count

</script>
</div>
```
""", "Show with syntax highlighting")

# ╔═╡ d121e085-c69b-490f-b315-c11a9abd57a6
details(md"""
	```htmlmixed
	<script>

	let data = $(films)

	// html`...` is from https://github.com/observablehq/stdlib
	// note the escaped dollar signs:
	let Film = ({title, director, year}) => html`
		<li class="film">
			<b>\${title}</b> by <em>\${director}</em> (\${year})
		</li>
	`

	// the returned HTML node is rendered
	return html`
		<ul>
			\${data.map(Film)}
		</ul>
	`

	</script>
	```
	""", "Show with syntax highlighting")

# ╔═╡ d4bdc4fe-2af8-402f-950f-2afaf77c62de
details(md"""
	```htmlmixed
	<script id="something">

	console.log("'this' is currently:", this)

	if(this == null) {
		return html`<blockquote>I am running for the first time!</blockqoute>`
	} else {
		return html`<blockquote><b>I was triggered by reactivity!</b></blockqoute>`
	}

	</script>
	```
	""", "Show with syntax highlighting")

# ╔═╡ e910982c-8508-4729-a75d-8b5b847918b6
details(md"""
```htmlmixed
<script src="https://cdn.jsdelivr.net/npm/d3@6.2.0/dist/d3.min.js"></script>

<script id="hello">

const positions = $(dot_positions)

const svg = this == null ? DOM.svg(600,200) : this
const s = this == null ? d3.select(svg) : this.s

s.selectAll("circle")
	.data(positions)
	.join("circle")
	.transition()
	.duration(300)
	.attr("cx", d => d)
	.attr("cy", 100)
	.attr("r", 10)
	.attr("fill", "gray")

const output = svg
output.s = s
return output
</script>
```
""", "Show with syntax highlighting")

# ╔═╡ 05d28aa2-9622-4e62-ab39-ca4c7dde6eb4
details(md"""
	```htmlmixed
	<script type="module" id="asdf">
		//await new Promise(r => setTimeout(r, 1000))

		const { html, render, Component, useEffect, useLayoutEffect, useState, useRef, useMemo, createContext, useContext, } = await import( "https://cdn.jsdelivr.net/npm/htm@3.0.4/preact/standalone.mjs")

		const node = this ?? document.createElement("div")

		const new_state = $(state)

		if(this == null){

			// PREACT APP STARTS HERE

			const Item = ({value}) => {
				const [loading, set_loading] = useState(true)

				useEffect(() => {
					set_loading(true)

					const handle = setTimeout(() => {
						set_loading(false)
					}, 1000)

					return () => clearTimeout(handle)
				}, [value])

				return html`<li>\${loading ?
					html`<em>Loading...</em>` :
					value
				}</li>`
			}

			const App = () => {

				const [state, set_state] = useState(new_state)
				node.set_app_state = set_state

				return html`<h5>Hello world!</h5>
					<ul>\${
					state.x.map((x,i) => html`<\${Item} value=\${x} key=\${i}/>`)
				}</ul>`;
			}

			// PREACT APP ENDS HERE

			render(html`<\${App}/>`, node);

		} else {

			node.set_app_state(new_state)
		}
		return node
	</script>
	```
	""", "Show with syntax highlighting")

# ╔═╡ cc318a19-316f-4fd9-8436-fb1d42f888a3
demo_img = let
	url = "https://user-images.githubusercontent.com/6933510/116753174-fa40ab80-aa06-11eb-94d7-88f4171970b2.jpeg"
	data = read(download(url))
	PlutoUI.Show(MIME"image/jpg"(), data)
end

# ╔═╡ 7aacdd8c-1571-4461-ba6e-0fd65dd8d788
demo_html = @htl("Hello!")

# ╔═╡ ebec177c-4c33-45a4-bdbd-f16944631aff
md"""
## Embeddable output

Pluto has a multimedia object viewer, which is used to display the result of a cell's output. Depending on the _type_ of the resulting object, the richest possible viewer is used. This includes:
- an interactive structure viewer for arrays, tables, dictionaries and more: $(embed_display([1,2,(a=3, b=4)]))
- an `<img>` tag with optimized data transfer for images: $(embed_display(demo_img))
- raw HTML for HTML-showable objects: $(embed_display(demo_html))

Normally, you get this object viewer for the _output_ of a cell. However, as demonstrated in the list above, you can also **embed Pluto's object viewer in your own HTML**. To do so, Pluto provides a function:
```julia
embed_display(x)
```

#### Example

As an example, here is how you display two arrays side-by-side:

```julia
@htl("\""

<div style="display: flex;">
$(embed_display(rand(4)))
$(embed_display(rand(4)))
</div>

"\"")

```

You can [learn more](https://github.com/fonsp/Pluto.jl/pull/1126) about how this feature works, or how to use it inside packages.
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
HypertextLiteral = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
HypertextLiteral = "~0.9.3"
PlutoUI = "~0.7.34"
"""

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

[[AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

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

[[ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

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

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

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

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

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

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

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

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

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

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

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

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

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

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

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

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

[[PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "670e559e5c8e191ded66fa9ea89c97f10376bb4c"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.38"

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

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

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

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

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

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

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

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

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

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

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─97914842-76d2-11eb-0c48-a7eedca870fb
# ╠═571613a1-6b4b-496d-9a68-aac3f6a83a4b
# ╟─168e13f7-2ff2-4207-be56-e57755041d36
# ╠═28ae1424-67dc-4b76-a172-1185cc76cb59
# ╟─93abe0dc-f041-475f-9ef7-d8ee4408414b
# ╟─ea39c63f-7466-4015-a66c-08bd9c716343
# ╟─8b082f9a-073e-4112-9422-4087850fc89e
# ╟─d70a3a02-ef3a-450f-bf5a-4a0d7f6262e2
# ╟─10cf6ed1-8276-4a4a-ad06-097d10335512
# ╠═d967cdf9-3df9-40bb-9b08-09cae95a5ca7
# ╟─858745a9-cd59-43a6-a296-803515518e57
# ╠═21a9e3e6-92f4-475d-9c8e-21e15c09336b
# ╟─4a3398be-ee86-45f3-ac8b-f627a38c00b8
# ╠═2d5fd611-284b-4428-b6a5-8909203990b9
# ╠═82de4674-9ecc-46c4-8a57-0b4453c579c3
# ╟─70a415be-881a-4c01-9f8c-635b8b89e1ad
# ╠═730a692f-2bf2-4d5b-86da-6ab861e8b8ac
# ╟─a45fdec4-2d4b-429b-b809-4c256b57fffe
# ╠═c68ebd7b-5fb6-4527-ac34-33f9730e4587
# ╟─8c03139f-a94b-40cc-859f-0d86f1c72143
# ╠═d8dcb044-0ac8-46d1-a043-1073bb6d1ff1
# ╠═bf592202-a9a4-4e9b-8433-fed55e3aa3bc
# ╟─e7d3db79-8253-4cbd-9832-5afb7dff0abf
# ╟─5ac5b984-8c02-4b8d-a342-d0f05f7909ec
# ╠═ef28eb8d-ec98-43e5-9012-3338c3b84f1b
# ╠═1ba370cc-3631-47ea-9db5-75587e8e4ff3
# ╠═7fcf2f3f-d902-4338-adf0-8ef181e79420
# ╟─7afbf8ef-e91c-45b9-bf22-24201cbb4828
# ╠═b226da72-9512-4d14-8582-2f7787c25028
# ╠═a6fd1f7b-a8fc-420d-a8bb-9f549842ad3e
# ╟─d12b98df-8c3f-4620-ba3c-2f3dadac521b
# ╟─965f3660-6ec4-4a86-a2a2-c167dbe9315f
# ╠═01ce31a9-6856-4ee7-8bce-7ce635167457
# ╠═00d97588-d591-4dad-9f7d-223c237deefd
# ╠═21f57310-9ceb-423c-a9ce-5beb1060a5a3
# ╟─94561cb1-2325-49b6-8b22-943923fdd91b
# ╟─7d9d6c28-131a-4b2a-84f8-5c085f387e85
# ╟─0866afc2-fd42-42b7-a572-9d824cf8b83b
# ╟─75e1a973-7ef0-4ac5-b3e2-5edb63577927
# ╠═e8d8a60e-489b-467a-b49c-1fa844807751
# ╟─b0c246ed-b871-461b-9541-280e49b49136
# ╠═9346d8e2-9ba0-4475-a21f-11bdd018bc60
# ╠═7822fdb7-bee6-40cc-a089-56bb32d77fe6
# ╟─701de4b8-42d3-46a3-a399-d7761dccd83d
# ╟─88120468-a43d-4d58-ac04-9cc7c86ca179
# ╠═ea4b2da1-4c83-4a1f-8fc3-c71a120e58e1
# ╟─08bdeaff-5bfb-49ab-b4cc-3a3446c63edc
# ╟─9b6b5da9-8372-4ebf-9c66-ae9fcfc45d47
# ╟─f18b98f7-1e0f-4273-896f-8a667d15605b
# ╟─d83d57e2-4787-4b8d-8669-64ed73d79e73
# ╟─077c95cf-2a1b-459f-830e-c29c11a2c5cc
# ╟─80511436-e41f-4913-8a30-d9e113cfaf71
# ╟─8388a833-d535-4cbd-a27b-de323cea60e8
# ╟─4cf27df3-6a69-402e-a71c-26538b2a52e7
# ╟─5721ad33-a51a-4a91-adb2-0915ea0efa13
# ╠═c857bb4b-4cf4-426e-b340-592cf7700434
# ╟─d121e085-c69b-490f-b315-c11a9abd57a6
# ╟─fc8984c8-4668-418a-b258-a1718809470c
# ╠═846354c8-ba3b-4be7-926c-d3c9cc9add5f
# ╟─a33c7d7a-8071-448e-abd6-4e38b5444a3a
# ╠═91f3dab8-5521-44a0-9890-8d988a994076
# ╠═dcaae662-4a4f-4dd3-8763-89ea9eab7d43
# ╟─d4bdc4fe-2af8-402f-950f-2afaf77c62de
# ╟─e77cfefc-429d-49db-8135-f4604f6a9f0b
# ╠═2d5689f5-1d63-4b8b-a103-da35933ad26e
# ╠═6dd221d1-7fd8-446e-aced-950512ea34bc
# ╠═0a9d6e2d-3a41-4cd5-9a4e-a9b76ed89fa9
# ╟─0962d456-1a76-4b0d-85ff-c9e7dc66621d
# ╠═bf9b36e8-14c5-477b-a54b-35ba8e415c77
# ╟─e910982c-8508-4729-a75d-8b5b847918b6
# ╟─781adedc-2da7-4394-b323-e508d614afae
# ╟─de789ad1-8197-48ae-81b2-a21ec2340ae0
# ╠═85483b28-341e-4ed6-bb1e-66c33613725e
# ╠═9e37c18c-3ebb-443a-9663-bb4064391d6e
# ╟─05d28aa2-9622-4e62-ab39-ca4c7dde6eb4
# ╠═3266f9e6-42ad-4103-8db3-b87d2c315290
# ╟─ebec177c-4c33-45a4-bdbd-f16944631aff
# ╟─da7091f5-8ba2-498b-aa8d-bbf3b4505b81
# ╠═64cbf19c-a4e3-4cdb-b4ec-1fbe24be55ad
# ╟─cc318a19-316f-4fd9-8436-fb1d42f888a3
# ╟─7aacdd8c-1571-4461-ba6e-0fd65dd8d788
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
module Aras

import GreyDecision.GreyNumbers: GreyNumber
import GreyDecision.Utility: unitize, apply_columns

function aras(decisionMat::Array{GreyNumber{T}, 2}, weights::Array{Float64,1},
    fns::Array{Function,1}) where {T <: Real}

    nrows, ncols = size(decisionMat)
    w = unitize(weights)
    referenceRow = apply_columns(fns, decisionMat)
    extendMat = [decisionMat; referenceRow']

    for col in 1:ncols
        if fns[col] == minimum
            extendMat[:, col] = map(x -> GreyNumber(1.0) / x, extendMat[:, col])
        end
    end

    normalizedMat = similar(extendMat)

    for col in 1:ncols
        for row in 1:nrows + 1
            normalizedMat[row, col] = extendMat[row, col] / sum(extendMat[:, col])
        end
    end

    optimality_degrees = zeros(GreyNumber{Float64}, nrows + 1)
    for i in 1:(nrows + 1)
        optimality_degrees[i] = (w .* normalizedMat[i, :]) |> sum
    end

    utility_degrees = zeros(GreyNumber{Float64}, nrows)
    for i in 1:nrows
        utility_degrees[i] = optimality_degrees[i] /  optimality_degrees[end]
    end

    orderings = sortperm(utility_degrees)
    bestIndex = orderings |> last

    result = Dict(
        :referenceRow       => referenceRow,
        :extendMat          => extendMat,
        :normalizedMat      => normalizedMat,
        :optimality_degrees => optimality_degrees,
        :scores    => utility_degrees,
        :orderings          => orderings,
        :bestIndex          => bestIndex
    )

    return result

end

end # end of module

ext: jl


text:
using StaticArrays

typealias Vec2d SVector{2,Float64};

type CartesianMesh
    nx
    ny
    Δx
    Δy
end

CartesianMesh(nx, ny) = CartesianMesh(nx, ny, 1.0/nx, 1.0/ny);

immutable Cell
    id  :: Int
    x   :: Vec2d
    vol :: Float64
end

immutable Face
    x     :: Vec2d
    s     :: Vec2d
    owner :: Int
    neigh :: Int
end;

immutable Patch
    name
    faces
end

function cell(msh::CartesianMesh, id)
    i,j = rem(id-1, msh.nx)+1, div(id-1, msh.nx)+1
    return Cell(id, Vec2d( (i-1/2)*msh.Δx, (j-1/2)*msh.Δy ), msh.Δx*msh.Δy )
end;

cells(msh::CartesianMesh) = @task begin
    id = 1
    vol = msh.Δx * msh.Δy
    for y in range(msh.Δy/2, msh.Δy, msh.ny)
        for x in range(msh.Δx/2, msh.Δx, msh.nx)
            produce( Cell(id, Vec2d(x,y), vol) )
            id += 1
        end
    end
end;

internalfaces(msh::CartesianMesh) = @task begin
    id(i,j) = i + (j-1)*msh.nx
    Δx, Δy = msh.Δx, msh.Δy

    s = Vec2d(Δy, 0.0)
    for j=1:msh.ny
        for i=1:msh.nx-1
            owner, neigh = id(i,j), id(i+1,j)
            x = Vec2d( i*Δx, (j-1/2)*Δy )
            produce( Face(x, s, owner, neigh) )
        end
    end

    s = Vec2d(0.0, Δx)
    for j=1:msh.ny-1
        for i=1:msh.nx
            owner, neigh = id(i,j), id(i,j+1)
            x = Vec2d( (i-1/2)*Δx, j*Δy )
            produce( Face(x, s, owner, neigh) )
        end
    end

end;

function patches(msh::CartesianMesh)
    id(i,j) = i + (j-1)*msh.nx
    Δx, Δy = msh.Δx, msh.Δy

    leftfaces() = @task begin
        s = Vec2d(-Δy, 0.0)
        for j=1:msh.ny
            produce( Face(Vec2d(0.0, (j-1/2)*Δy), s, id(1,j), 0) )
        end
    end

    rightfaces() = @task begin
        s = Vec2d(Δy, 0.0)
        for j=1:msh.ny
            produce( Face(Vec2d(msh.nx*Δx, (j-1/2)*Δy), s, id(msh.nx,j), 0) )
        end
    end

    bottomfaces() = @task begin
        s = Vec2d(0.0,-Δx)
        for i=1:msh.nx
            produce( Face(Vec2d((i-1/2)*Δx, 0.0), s, id(i,1), 0) )
        end
    end

    topfaces() = @task begin
        s = Vec2d(0.0,Δx)
        for i=1:msh.nx
            produce( Face(Vec2d((i-1/2)*Δx, msh.ny*Δy), s, id(i,msh.ny), 0) )
        end
    end

    return [ Patch("left", leftfaces), Patch("right", rightfaces),
        Patch("bottom", bottomfaces), Patch("top", topfaces) ]
end;

name(p::Patch) = p.name;
boundaryfaces(p::Patch) = p.faces();

max_stars_repo_name: furstj/NSP

max_issues_repo_name: furstj/NSP

max_forks_repo_name: furstj/NSP


text:
using Documenter, Literate, LikelihoodfreeInference, CMAEvolutionStrategy, BlackBoxOptim, SimultaneousPerturbationStochasticApproximation, BayesianOptimization, StatsPlots

OUTDIR = joinpath(@__DIR__, "src", "generated")
for example in ("toyexample.jl", "blowfly.jl")
    infile = joinpath(@__DIR__, "literate", example)
    Literate.markdown(infile, OUTDIR, documenter = true)
    Literate.notebook(infile, OUTDIR, documenter = true, execute = true)
end

makedocs(
    modules = [LikelihoodfreeInference],
    format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
    authors = "Johanni Brea",
    sitename = "LikelihoodfreeInference.jl",
    pages = Any[
                "index.md",
                "generated/toyexample.md",
                "generated/blowfly.md",
                "reference.md"
               ]
    # strict = true,
    # clean = true,
    # checkdocs = :exports,
)

deploydocs(
    repo = "github.com/jbrea/LikelihoodfreeInference.jl.git",
    push_preview = true
)


text:
# This file is part of the TaylorIntegration.jl package; MIT licensed

using TaylorIntegration, DiffEqBase
using Test
using LinearAlgebra: norm
using Elliptic

@testset "Testing `taylorize.jl`" begin

    # Constants for the integrations
    local _order = 20
    local _abstol = 1.0E-20
    local t0 = 0.0
    local tf = 1000.0

    # Scalar integration
    @testset "Scalar case: xdot(x, p, t) = b-x^2" begin
        b1 = 3.0
        @taylorize xdot1(x, p, t) = b1-x^2
        @test length(methods(TaylorIntegration.jetcoeffs!)) == 3
        @test (@isdefined xdot1)

        x0 = 1.0
        tv1, xv1 = taylorinteg( xdot1, x0, t0, tf, _order, _abstol, maxsteps=1000,
            parse_eqs=false)
        tv1p, xv1p = taylorinteg( xdot1, x0, t0, tf, _order, _abstol, maxsteps=1000)

        @test length(tv1) == length(tv1p)
        @test iszero( norm(tv1-tv1p, Inf) )
        @test iszero( norm(xv1-xv1p, Inf) )

        # Now using `local` constants
        @taylorize xdot2(x, p, t) = (local b2 = 3; b2-x^2)
        @test length(methods(TaylorIntegration.jetcoeffs!)) == 4
        @test (@isdefined xdot2)

        tv2, xv2 = taylorinteg( xdot2, x0, t0, tf, _order, _abstol, maxsteps=1000,
            parse_eqs=false)
        tv2p, xv2p = taylorinteg( xdot2, x0, t0, tf, _order, _abstol, maxsteps=1000)

        @test length(tv2) == length(tv2p)
        @test iszero( norm(tv2-tv2p, Inf) )
        @test iszero( norm(xv2-xv2p, Inf) )

        # Passing a parameter
        @taylorize xdot3(x, p, t) = p-x^2
        @test length(methods(TaylorIntegration.jetcoeffs!)) == 5
        @test (@isdefined xdot3)

        tv3, xv3 = taylorinteg( xdot3, x0, t0, tf, _order, _abstol, b1, maxsteps=1000,
            parse_eqs=false)
        tv3p, xv3p = taylorinteg( xdot3, x0, t0, tf, _order, _abstol, b1, maxsteps=1000)

        @test length(tv3) == length(tv3p)
        @test iszero( norm(tv3-tv3p, Inf) )
        @test iszero( norm(xv3-xv3p, Inf) )

        # Comparing integrations
        @test length(tv1p) == length(tv2p) == length(tv3p)
        @test iszero( norm(tv1p-tv2p, Inf) )
        @test iszero( norm(tv1p-tv3p, Inf) )
        @test iszero( norm(xv1p-xv2p, Inf) )
        @test iszero( norm(xv1p-xv3p, Inf) )

        xv4 = taylorinteg( xdot2, x0, t0:0.5:tf, _order, _abstol, maxsteps=1000,
            parse_eqs=false)
        xv4p = taylorinteg( xdot2, x0, t0:0.5:tf, _order, _abstol, maxsteps=1000)
        @test iszero( norm(xv4-xv4p, Inf) )

        # Compare to exact solution
        exact_sol(t, b, x0) = sqrt(b)*((sqrt(b)+x0)-(sqrt(b)-x0)*exp(-2sqrt(b)*t)) /
            ((sqrt(b)+x0)+(sqrt(b)-x0)*exp(-2sqrt(b)*t))
        @test norm(xv1p[end] - exact_sol(tv1p[end], b1, x0), Inf) < 1.0e-15
    end

    @testset "Scalar case: xdot(x, p, t) = -10" begin
        xdot1(x, p, t) = -10 + zero(t) # `zero(t)` is needed; cf #20
        @taylorize xdot1_parsed(x, p, t) = -10 + zero(t) # `zero(t)` can be avoided here

        @test (@isdefined xdot1_parsed)
        tv1, xv1   = taylorinteg( xdot1, 10, 1, 20.0, _order, _abstol)
        tv1p, xv1p = taylorinteg( xdot1_parsed, 10.0, 1.0, 20.0, _order, _abstol)

        @test length(tv1) == length(tv1p)
        @test iszero( norm(tv1-tv1p, Inf) )
        @test iszero( norm(xv1-xv1p, Inf) )

        # Now using `local`
        @taylorize function xdot2(x, p, t)
            local ggrav = -10 + zero(t)  # zero(t) is needed for the parse_eqs=false case
            tmp = ggrav  # needed to avoid an error when parsing
        end

        @test (@isdefined xdot2)

        tv2, xv2   = taylorinteg( xdot2, 10.0, 1.0, 20.0, _order, _abstol, parse_eqs=false)
        tv2p, xv2p = taylorinteg( xdot2, 10.0, 1, 20.0, _order, _abstol)

        @test length(tv2) == length(tv2p)
        @test iszero( norm(tv2-tv2p, Inf) )
        @test iszero( norm(xv2-xv2p, Inf) )

        # Passing a parameter
        @taylorize xdot3(x, p, t) = p + zero(t) # `zero(t)` can be avoided here

        @test (@isdefined xdot3)

        tv3, xv3   = taylorinteg( xdot3, 10.0, 1.0, 20.0, _order, _abstol, -10,
            parse_eqs=false)
        tv3p, xv3p = taylorinteg( xdot3, 10, 1.0, 20.0, _order, _abstol, -10)

        @test length(tv3) == length(tv3p)
        @test iszero( norm(tv3-tv3p, Inf) )
        @test iszero( norm(xv3-xv3p, Inf) )

        # Comparing integrations
        @test length(tv1p) == length(tv2p) == length(tv3p)
        @test iszero( norm(tv1p-tv2p, Inf) )
        @test iszero( norm(tv1p-tv3p, Inf) )
        @test iszero( norm(xv1p-xv2p, Inf) )
        @test iszero( norm(xv1p-xv3p, Inf) )

        # Compare to exact solution
        exact_sol(t, g, x0) = x0 + g*(t-1.0)
        @test norm(xv1p[end] - exact_sol(20, -10, 10), Inf) < 10eps(exact_sol(20, -10, 10))
    end

    # Pendulum integrtf = 100.0
    @testset "Integration of the pendulum and DiffEqs interface" begin
        @taylorize function pendulum!(dx, x, p, t)
            dx[1] = x[2]
            dx[2] = -sin( x[1] )
            nothing
        end

        @test (@isdefined pendulum!)
        q0 = [pi-0.001, 0.0]
        tv2, xv2 = taylorinteg(pendulum!, q0, t0, tf, _order, _abstol, parse_eqs=false,
            maxsteps=5000)
        tv2p, xv2p = taylorinteg(pendulum!, q0, t0, tf, _order, _abstol, maxsteps=5000)

        @test length(tv2) == length(tv2p)
        @test iszero( norm(tv2-tv2p, Inf) )
        @test iszero( norm(xv2-xv2p, Inf) )

        prob = ODEProblem(pendulum!, q0, (t0, tf), nothing) # no parameters
        sol1 = solve(prob, TaylorMethod(_order), abstol=_abstol, parse_eqs=true)
        sol2 = solve(prob, TaylorMethod(_order), abstol=_abstol, parse_eqs=false)
        sol3 = solve(prob, TaylorMethod(_order), abstol=_abstol)

        @test sol1.t == sol2.t == sol3.t == tv2p
        @test sol1.u[end] == sol2.u[end] == sol3.u[end] == xv2p[end,1:2]
    end

    # Complex dependent variables
    @testset "Complex dependent variable" begin
        cc = complex(0.0,1.0)
        @taylorize eqscmplx(x, p, t) = cc*x

        @test (@isdefined eqscmplx)
        cx0 = complex(1.0, 0.0)
        tv1, xv1 = taylorinteg(eqscmplx, cx0, t0, tf, _order, _abstol, maxsteps=1500,
            parse_eqs=false)
        tv1p, xv1p = taylorinteg(eqscmplx, cx0, t0, tf, _order, _abstol, maxsteps=1500)

        @test length(tv1) == length(tv1p)
        @test iszero( norm(tv1-tv1p, Inf) )
        @test iszero( norm(xv1-xv1p, Inf) )

        # Using `local` for the constant value
        @taylorize eqscmplx2(x, p, t) = (local cc1 = Complex(0.0,1.0); cc1*x)

        @test (@isdefined eqscmplx2)
        tv2, xv2 = taylorinteg(eqscmplx2, cx0, t0, tf, _order, _abstol, maxsteps=1500,
            parse_eqs=false)
        tv2p, xv2p = taylorinteg(eqscmplx2, cx0, t0, tf, _order, _abstol, maxsteps=1500)

        @test length(tv2) == length(tv2p)
        @test iszero( norm(tv2-tv2p, Inf) )
        @test iszero( norm(xv2-xv2p, Inf) )

        # Passing a parameter
        @taylorize eqscmplx3(x, p, t) = p*x
        @test (@isdefined eqscmplx3)
        tv3, xv3 = taylorinteg(eqscmplx3, cx0, t0, tf, _order, _abstol, cc, maxsteps=1500,
            parse_eqs=false)
        tv3p, xv3p = taylorinteg(eqscmplx3, cx0, t0, tf, _order, _abstol, cc, maxsteps=1500)

        @test length(tv3) == length(tv3p)
        @test iszero( norm(tv3-tv3p, Inf) )
        @test iszero( norm(xv3-xv3p, Inf) )

        # Comparing integrations
        @test length(tv1p) == length(tv2p) == length(tv3p)
        @test iszero( norm(tv1p-tv2p, Inf) )
        @test iszero( norm(tv1p-tv3p, Inf) )
        @test iszero( norm(xv1p-xv2p, Inf) )
        @test iszero( norm(xv1p-xv3p, Inf) )

        # Compare to exact solution
        exact_sol(t, p, x0) = x0*exp(p*t)
        nn = norm.(abs2.(xv1p)-abs2.(exact_sol.(tv1p, cc, cx0)), Inf)
        @test maximum(nn) < 1.0e-14
    end

    @testset "Time-dependent integration (with and without `local` vars)" begin
        @taylorize function integ_cos1(x, p, t)
            y = cos(t)
            return y
        end
        @taylorize function integ_cos2(x, p, t)
            local y = cos(t)  # allows to calculate directly `cos(t)` *once*
            yy = y   # needed to avoid an error
            return yy
        end

        @test (@isdefined integ_cos1)
        @test (@isdefined integ_cos2)

        tv11, xv11 = taylorinteg(integ_cos1, 0.0, 0.0, pi, _order, _abstol, parse_eqs=false)
        tv12, xv12 = taylorinteg(integ_cos1, 0.0, 0.0, pi, _order, _abstol)
        @test length(tv11) == length(tv12)
        @test iszero( norm(tv11-tv12, Inf) )
        @test iszero( norm(xv11-xv12, Inf) )

        tv21, xv21 = taylorinteg(integ_cos2, 0.0, 0.0, pi, _order, _abstol, parse_eqs=false)
        tv22, xv22 = taylorinteg(integ_cos2, 0.0, 0.0, pi, _order, _abstol)
        @test length(tv21) == length(tv22)
        @test iszero( norm(tv21-tv22, Inf) )
        @test iszero( norm(xv21-xv22, Inf) )

        @test iszero( norm(tv12-tv22, Inf) )
        @test iszero( norm(xv12-xv22, Inf) )

        # Compare to exact solution
        @test norm(xv11[end] - sin(tv11[end]), Inf) < 1.0e-15
    end

    @testset "Time-dependent integration (vectorial case)" begin
        @taylorize function integ_vec(dx, x, p, t)
            dx[1] = cos(t)
            dx[2] = -sin(t)
            return dx
        end

        @test (@isdefined integ_vec)
        x0 = [0.0, 1.0]
        tv11, xv11 = taylorinteg(integ_vec, x0, 0.0, pi, _order, _abstol, parse_eqs=false)
        tv12, xv12 = taylorinteg(integ_vec, x0, 0.0, pi, _order, _abstol)
        @test length(tv11) == length(tv12)
        @test iszero( norm(tv11-tv12, Inf) )
        @test iszero( norm(xv11-xv12, Inf) )

        # Compare to exact solution
        @test norm(xv12[end, 1] - sin(tv12[end]), Inf) < 1.0e-15
        @test norm(xv12[end, 2] - cos(tv12[end]), Inf) < 1.0e-15
    end

    # Simple harmonic oscillator
    @testset "Simple harmonic oscillator" begin
        @taylorize function harm_osc!(dx, x, p, t)
            local ω = p[1]
            dx[1] = x[2]
            dx[2] = -ω^2 * x[1]
            return nothing
        end

        @test (@isdefined harm_osc!)
        q0 = [1.0, 0.0]
        p = [2.0]
        tv1, xv1 = taylorinteg(harm_osc!, q0, t0, tf, _order, _abstol,
            p, maxsteps=1000, parse_eqs=false)
        tv1p, xv1p = taylorinteg(harm_osc!, q0, t0, tf, _order, _abstol,
            p, maxsteps=1000)

        @test length(tv1) == length(tv1p)
        @test iszero( norm(tv1-tv1p, Inf) )
        @test iszero( norm(xv1-xv1p, Inf) )

        # Comparing to exact solution
        @test norm(xv1p[end, 1] - cos(p[1]*tv1p[end]), Inf) < 1.0e-12
        @test norm(xv1p[end, 2] + p[1]*sin(p[1]*tv1p[end]), Inf) < 3.0e-12
    end

    @testset "Multiple pendula" begin
        NN = 3
        nnrange = 1:3
        @taylorize function multpendula1!(dx, x, p, t)
            for i in p[2]
                dx[i] = x[p[1]+i]
                dx[i+p[1]] = -sin( x[i] )
            end
            return nothing
        end

        @test (@isdefined multpendula1!)
        q0 = [pi-0.001, 0.0, pi-0.001, 0.0,  pi-0.001, 0.0]
        tv1, xv1 = taylorinteg(multpendula1!, q0, t0, tf, _order, _abstol,
            [NN, nnrange], maxsteps=1000, parse_eqs=false)
        tv1p, xv1p = taylorinteg(multpendula1!, q0, t0, tf, _order, _abstol,
            [NN, nnrange], maxsteps=1000)

        @test length(tv1) == length(tv1p)
        @test iszero( norm(tv1-tv1p, Inf) )
        @test iszero( norm(xv1-xv1p, Inf) )

        @taylorize function multpendula2!(dx, x, p, t)
            local NN = 3
            local nnrange = 1:NN
            for i in nnrange
                dx[i] = x[NN+i]
                dx[i+NN] = -sin( x[i] )
            end
            return nothing
        end

        @test (@isdefined multpendula2!)
        tv2, xv2 = taylorinteg(multpendula2!, q0, t0, tf, _order, _abstol,
            maxsteps=1000, parse_eqs=false)
        tv2p, xv2p = taylorinteg(multpendula2!, q0, t0, tf, _order, _abstol,
            maxsteps=1000)

        @test length(tv2) == length(tv2p)
        @test iszero( norm(tv2-tv2p, Inf) )
        @test iszero( norm(xv2-xv2p, Inf) )

        @taylorize function multpendula3!(dx, x, p, t)
            nn, nnrng = p    # `local` is not needed to reassign `p`; internally is treated as local)
            for i in nnrng
                dx[i] = x[nn+i]
                dx[i+nn] = -sin( x[i] )
            end
            return nothing
        end

        @test (@isdefined multpendula3!)
        tv3, xv3 = taylorinteg(multpendula3!, q0, t0, tf, _order, _abstol,
            [NN, nnrange], maxsteps=1000, parse_eqs=false)
        tv3p, xv3p = taylorinteg(multpendula3!, q0, t0, tf, _order, _abstol,
            [NN, nnrange], maxsteps=1000)

        # Comparing integrations
        @test length(tv1) == length(tv2) == length(tv3)
        @test iszero( norm(tv1-tv2, Inf) )
        @test iszero( norm(tv1-tv3, Inf) )
        @test iszero( norm(xv1-xv2, Inf) )
        @test iszero( norm(xv1-xv3, Inf) )
    end

    # Kepler problem
    # Redefining some constants
    local _order = 28
    local tf = 2π*100.0
    @testset "Kepler problem (using `^`)" begin
        @taylorize function kepler1!(dq, q, p, t)
            μ = p
            r_p3d2 = (q[1]^2+q[2]^2)^1.5

            dq[1] = q[3]
            dq[2] = q[4]
            dq[3] = μ * q[1]/r_p3d2
            dq[4] = μ * q[2]/r_p3d2

            return nothing
        end

        @taylorize function kepler2!(dq, q, p, t)
            nn, μ = p
            r2 = zero(q[1])
            for i = 1:nn
                r2_aux = r2 + q[i]^2
                r2 = r2_aux
            end
            r_p3d2 = r2^(3/2)
            for j = 1:nn
                dq[j] = q[nn+j]
                dq[nn+j] = μ*q[j]/r_p3d2
            end

            nothing
        end

        @test (@isdefined kepler1!)
        @test (@isdefined kepler2!)
        pars = (2, -1.0)
        q0 = [0.2, 0.0, 0.0, 3.0]
        tv1, xv1 = taylorinteg(kepler1!, q0, t0, tf, _order, _abstol, -1.0,
            maxsteps=500000, parse_eqs=false)
        tv1p, xv1p = taylorinteg(kepler1!, q0, t0, tf, _order, _abstol, -1.0,
            maxsteps=500000)

        tv6p, xv6p = taylorinteg(kepler2!, q0, t0, tf, _order, _abstol, pars,
            maxsteps=500000, parse_eqs=false)
        tv7p, xv7p = taylorinteg(kepler2!, q0, t0, tf, _order, _abstol, pars,
            maxsteps=500000)

        @test length(tv1) == length(tv1p) == length(tv6p)
        @test iszero( norm(tv1-tv1p, Inf) )
        @test iszero( norm(xv1-xv1p, Inf) )
        @test iszero( norm(tv1-tv6p, Inf) )
        @test iszero( norm(xv1-xv6p, Inf) )
        @test iszero( norm(tv7p-tv6p, Inf) )
        @test iszero( norm(xv7p-xv6p, Inf) )

        @taylorize function kepler3!(dq, q, p, t)
            local mμ = -1.0
            r_p3d2 = (q[1]^2+q[2]^2)^1.5

            dq[1] = q[3]
            dq[2] = q[4]
            dq[3] = mμ*q[1]/r_p3d2
            dq[4] = mμ*q[2]/r_p3d2

            return nothing
        end

        @taylorize function kepler4!(dq, q, p, t)
            local mμ = -1.0
            local NN = 2
            r2 = zero(q[1])
            for i = 1:NN
                r2_aux = r2 + q[i]^2
                r2 = r2_aux
            end
            r_p3d2 = r2^(3/2)
            for j = 1:NN
                dq[j] = q[NN+j]
                dq[NN+j] = mμ*q[j]/r_p3d2
            end

            nothing
        end

        @test (@isdefined kepler3!)
        @test (@isdefined kepler4!)
        tv3, xv3 = taylorinteg(kepler3!, q0, t0, tf, _order, _abstol,
            maxsteps=500000, parse_eqs=false)
        tv3p, xv3p = taylorinteg(kepler3!, q0, t0, tf, _order, _abstol,
            maxsteps=500000)

        tv4, xv4 = taylorinteg(kepler4!, q0, t0, tf, _order, _abstol,
            maxsteps=500000, parse_eqs=false)
        tv4p, xv4p = taylorinteg(kepler4!, q0, t0, tf, _order, _abstol,
            maxsteps=500000)

        @test length(tv3) == length(tv3p) == length(tv4)
        @test iszero( norm(tv3-tv3p, Inf) )
        @test iszero( norm(xv3-xv3p, Inf) )
        @test iszero( norm(tv3-tv4, Inf) )
        @test iszero( norm(xv3-xv4, Inf) )
        @test iszero( norm(tv4p-tv4, Inf) )
        @test iszero( norm(xv4p-xv4, Inf) )

        # Comparing both integrations
        @test iszero( norm(tv1p-tv3p, Inf) )
        @test iszero( norm(xv1p-xv3p, Inf) )
    end

    @testset "Kepler problem (using `sqrt`)" begin
        @taylorize function kepler1!(dq, q, p, t)
            local μ = -1.0
            r = sqrt(q[1]^2+q[2]^2)
            r_p3d2 = r^3

            dq[1] = q[3]
            dq[2] = q[4]
            dq[3] = μ * q[1] / r_p3d2
            dq[4] = μ * q[2] / r_p3d2

            return nothing
        end

        @taylorize function kepler2!(dq, q, p, t)
            local NN = 2
            μ = p
            r2 = zero(q[1])
            for i = 1:NN
                r2_aux = r2 + q[i]^2
                r2 = r2_aux
            end
            r = sqrt(r2)
            r_p3d2 = r^3
            for j = 1:NN
                dq[j] = q[NN+j]
                dq[NN+j] =  μ * q[j] / r_p3d2
            end

            nothing
        end

        @test (@isdefined kepler1!)
        @test (@isdefined kepler2!)
        q0 = [0.2, 0.0, 0.0, 3.0]
        tv1, xv1 = taylorinteg(kepler1!, q0, t0, tf, _order, _abstol,
            maxsteps=500000, parse_eqs=false)
        tv1p, xv1p = taylorinteg(kepler1!, q0, t0, tf, _order, _abstol,
            maxsteps=500000)

        tv2, xv2 = taylorinteg(kepler2!, q0, t0, tf, _order, _abstol, -1.0,
            maxsteps=500000, parse_eqs=false)
        tv2p, xv2p = taylorinteg(kepler2!, q0, t0, tf, _order, _abstol, -1.0,
            maxsteps=500000)

        @test length(tv1) == length(tv1p) == length(tv2)
        @test iszero( norm(tv1-tv1p, Inf) )
        @test iszero( norm(xv1-xv1p, Inf) )
        @test iszero( norm(tv1-tv2, Inf) )
        @test iszero( norm(xv1-xv2, Inf) )
        @test iszero( norm(tv2p-tv2, Inf) )
        @test iszero( norm(xv2p-xv2, Inf) )

        @taylorize function kepler3!(dq, q, p, t)
            μ = p
            x, y, px, py = q
            r = sqrt(x^2+y^2)
            r_p3d2 = r^3

            dq[1] = px
            dq[2] = py
            dq[3] = μ * x / r_p3d2
            dq[4] = μ * y / r_p3d2

            return nothing
        end

        @taylorize function kepler4!(dq, q, p, t)
            local NN = 2
            local μ = -1.0
            r2 = zero(q[1])
            for i = 1:NN
                r2_aux = r2 + q[i]^2
                r2 = r2_aux
            end
            r = sqrt(r2)
            r_p3d2 = r^3
            for j = 1:NN
                dq[j] = q[NN+j]
                dq[NN+j] = μ * q[j] / r_p3d2
            end

            nothing
        end

        @test (@isdefined kepler3!)
        @test (@isdefined kepler4!)
        tv3, xv3 = taylorinteg(kepler3!, q0, t0, tf, _order, _abstol, -1.0,
            maxsteps=500000, parse_eqs=false)
        tv3p, xv3p = taylorinteg(kepler3!, q0, t0, tf, _order, _abstol, -1.0,
            maxsteps=500000)

        tv4, xv4 = taylorinteg(kepler4!, q0, t0, tf, _order, _abstol,
            maxsteps=500000, parse_eqs=false)
        tv4p, xv4p = taylorinteg(kepler4!, q0, t0, tf, _order, _abstol,
            maxsteps=500000)

        @test length(tv3) == length(tv3p) == length(tv4)
        @test iszero( norm(tv3-tv3p, Inf) )
        @test iszero( norm(xv3-xv3p, Inf) )
        @test iszero( norm(tv3-tv4, Inf) )
        @test iszero( norm(xv3-xv4, Inf) )
        @test iszero( norm(tv4p-tv4, Inf) )
        @test iszero( norm(xv4p-xv4, Inf) )

        # Comparing both integrations
        @test iszero( norm(tv1p-tv3p, Inf) )
        @test iszero( norm(xv1p-xv3p, Inf) )
    end

    local tf = 20.0
    @testset "Lyapunov spectrum and `@taylorize`" begin
        #Lorenz system parameters
        params = [16.0, 45.92, 4.0]

        #Lorenz system ODE:
        @taylorize function lorenz1!(dq, q, p, t)
            x, y, z = q
            σ, ρ, β = p
            dq[1] = σ*(y-x)
            dq[2] = x*(ρ-z) - y
            dq[3] = x*y - β*z
            nothing
        end

        #Lorenz system Jacobian (in-place):
        function lorenz1_jac!(jac, q, p, t)
            x, y, z = q
            σ, ρ, β = p
            jac[1,1] = -σ+zero(x)
            jac[2,1] = ρ-z
            jac[3,1] = y
            jac[1,2] = σ+zero(x)
            jac[2,2] = -one(y)
            jac[3,2] = x
            jac[1,3] = zero(z)
            jac[2,3] = -x
            jac[3,3] = -β+zero(x)
            nothing
        end

        q0 = [19.0, 20.0, 50.0] #the initial condition
        xi = set_variables("δ", order=1, numvars=length(q0))

        tv1, lv1, xv1 = lyap_taylorinteg(lorenz1!, q0, t0, tf, _order, _abstol,
            params, maxsteps=2000, parse_eqs=false);

        tv1p, lv1p, xv1p = lyap_taylorinteg(lorenz1!, q0, t0, tf, _order, _abstol,
            params, maxsteps=2000);

        @test tv1 == tv1p
        @test lv1 == lv1p
        @test xv1 == xv1p

        tv2, lv2, xv2 = lyap_taylorinteg(lorenz1!, q0, t0, tf, _order, _abstol,
            params, lorenz1_jac!, maxsteps=2000, parse_eqs=false);

        tv2p, lv2p, xv2p = lyap_taylorinteg(lorenz1!, q0, t0, tf, _order, _abstol,
            params, lorenz1_jac!, maxsteps=2000,  parse_eqs=true);

        @test tv2 == tv2p
        @test lv2 == lv2p
        @test xv2 == xv2p

        # Comparing both integrations (lorenz1)
        @test tv1 == tv2
        @test lv1 == lv2
        @test xv1 == xv2

        #Lorenz system ODE:
        @taylorize function lorenz2!(dx, x, p, t)
            #Lorenz system parameters
            local σ = 16.0
            local β = 4.0
            local ρ = 45.92

            dx[1] = σ*(x[2]-x[1])
            dx[2] = x[1]*(ρ-x[3])-x[2]
            dx[3] = x[1]*x[2]-β*x[3]
            nothing
        end

        #Lorenz system Jacobian (in-place):
        function lorenz2_jac!(jac, x, p, t)
            #Lorenz system parameters
            local σ = 16.0
            local β = 4.0
            local ρ = 45.92

            jac[1,1] = -σ+zero(x[1])
            jac[2,1] = ρ-x[3]
            jac[3,1] = x[2]
            jac[1,2] = σ+zero(x[1])
            jac[2,2] = -1.0+zero(x[1])
            jac[3,2] = x[1]
            jac[1,3] = zero(x[1])
            jac[2,3] = -x[1]
            jac[3,3] = -β+zero(x[1])
            nothing
        end

        tv3, lv3, xv3 = lyap_taylorinteg(lorenz2!, q0, t0, tf, _order, _abstol,
            maxsteps=2000, parse_eqs=false);

        tv3p, lv3p, xv3p = lyap_taylorinteg(lorenz2!, q0, t0, tf, _order, _abstol,
            maxsteps=2000);

        @test tv3 == tv3p
        @test lv3 == lv3p
        @test xv3 == xv3p

        tv4, lv4, xv4 = lyap_taylorinteg(lorenz2!, q0, t0, tf, _order, _abstol,
            lorenz2_jac!, maxsteps=2000, parse_eqs=false);

        tv4p, lv4p, xv4p = lyap_taylorinteg(lorenz2!, q0, t0, tf, _order, _abstol,
            lorenz2_jac!, maxsteps=2000,  parse_eqs=true);

        @test tv4 == tv4p
        @test lv4 == lv4p
        @test xv4 == xv4p

        # Comparing both integrations (lorenz2)
        @test tv3 == tv4
        @test lv3 == lv4
        @test xv3 == xv4

        # Comparing both integrations (lorenz1! vs lorenz2!)
        @test tv1 == tv3
        @test lv1 == lv3
        @test xv1 == xv3

        # Using ranges
        lv5, xv5 = lyap_taylorinteg(lorenz2!, q0, t0:0.125:tf, _order, _abstol,
            maxsteps=2000, parse_eqs=false);

        lv5p, xv5p = lyap_taylorinteg(lorenz2!, q0, t0:0.125:tf, _order, _abstol,
            maxsteps=2000, parse_eqs=true);

        @test lv5 == lv5p
        @test xv5 == xv5p

        lv6, xv6 = lyap_taylorinteg(lorenz2!, q0, t0:0.125:tf, _order, _abstol,
            lorenz2_jac!, maxsteps=2000, parse_eqs=false);

        lv6p, xv6p = lyap_taylorinteg(lorenz2!, q0, t0:0.125:tf, _order, _abstol,
            lorenz2_jac!, maxsteps=2000,  parse_eqs=true);

        @test lv6 == lv6p
        @test xv6 == xv6p
    end

    @testset "Tests for throwing errors" begin
        # Wrong number of arguments
        ex = :(function f_p!(dx, x, p, t, y)
            dx[1] = x[2]
            dx[2] = -sin( x[1] )
        end)
        @test_throws ArgumentError TaylorIntegration._make_parsed_jetcoeffs(ex)

        # `&&` is not yet implemented
        ex = :(function f_p!(x, p, t)
            true && x
        end)
        @test_throws ArgumentError TaylorIntegration._make_parsed_jetcoeffs(ex)

        # a is not an Expr; String
        ex = :(function f_p!(x, p, t)
            "a"
        end)
        @test_throws ArgumentError TaylorIntegration._make_parsed_jetcoeffs(ex)

        # KeyError: key :fname not found
        ex = :(begin
            x=1
            x+x
        end)
        @test_throws KeyError TaylorIntegration._make_parsed_jetcoeffs(ex)

        # BoundsError; no return variable defined
        ex = :(function f_p!(x, p, t)
            local cos(t)
        end)
        @test_throws BoundsError TaylorIntegration._make_parsed_jetcoeffs(ex)
    end

    @testset "Jet transport with @taylorize macro" begin
        @taylorize function pendulum!(dx, x, p, t)
            dx[1] = x[2]
            dx[2] = -sin( x[1] )
            nothing
        end

        varorder = 2 #the order of the variational expansion
        p = set_variables("ξ", numvars=2, order=varorder) #TaylorN steup
        q0 = [1.3, 0.0] #the initial conditions
        q0TN = q0 + p #parametrization of a small neighbourhood around the initial conditions
        # T is the librational period == 4Elliptic.K(sin(q0[1]/2)^2)
        T = 4Elliptic.K(sin(q0[1]/2)^2) # equals 7.019250311844546
        integstep = 0.25*T #the time interval between successive evaluations of the solution vector

        #the time range
        tr = t0:integstep:T;
        #note that as called below, taylorinteg uses the parsed jetcoeffs! method by default
        xvp = taylorinteg(pendulum!, q0, tr, _order, _abstol, maxsteps=100)

        # "warmup" for jet transport integration
        xvTN = taylorinteg(pendulum!, q0TN, tr, _order, _abstol, maxsteps=1, parse_eqs=false)
        xvTN = taylorinteg(pendulum!, q0TN, tr, _order, _abstol, maxsteps=1)
        @test size(xvTN) == (5,2)
        #jet transport integration with parsed jetcoeffs!
        xvTNp = taylorinteg(pendulum!, q0TN, tr, _order, _abstol, maxsteps=100)
        #jet transport integration with non-parsed jetcoeffs!
        xvTN = taylorinteg(pendulum!, q0TN, tr, _order, _abstol, maxsteps=100, parse_eqs=false)
        @test xvTN == xvTNp
        @test norm(xvTNp[:,:]() - xvp, Inf) < 1E-15

        dq = 0.0001rand(2)
        q1 = q0 + dq
        yv = taylorinteg(pendulum!, q1, tr, _order, _abstol, maxsteps=100)
        yv_jt = xvTNp[:,:](dq)
        @test norm(yv-yv_jt, Inf) < 1E-11

        dq = 0.001
        t = Taylor1([0.0, 1.0], 10)
        x0T1 = q0+[0t,t]
        q1 = q0+[0.0,dq]
        tv, xv = taylorinteg(pendulum!, q1, t0, 2T, _order, _abstol)
        tvT1, xvT1 = taylorinteg(pendulum!, x0T1, t0, 2T, _order, _abstol, parse_eqs=false)
        tvT1p, xvT1p = taylorinteg(pendulum!, x0T1, t0, 2T, _order, _abstol)
        @test tvT1 == tvT1p
        @test xvT1 == xvT1p
        xv_jt = xvT1p[:,:](dq)
        @test norm(xv_jt[end,:]-xv[end,:]) < 20eps(norm(xv[end,:]))
    end
end


text:
module Analyzer
using  ..Helper
export analyze

include("SlurpOptimizations.jl")
include("Function.jl")

import .Function: analyze

end


text:
y = 10
e = :(x + $(y / 2))     # y/2 evaluated here
@assert e == :(x + 5)

max_stars_repo_name: idrougge/ple

max_issues_repo_name: idrougge/ple

max_forks_repo_name: idrougge/ple


text:
"
https://arxiv.org/pdf/cond-mat/0203500.pdf
"
function dynamicaldmrg(A::Union{MPSComoving,FiniteMPS},z,ham::MPOHamiltonian;init=copy(A),solvtol=Defaults.tol,tol=solvtol*length(A)*2,maxiter=Defaults.maxiter,verbose=Defaults.verbose)
    w=real(z);eta=imag(z)

    envs1 = environments(init,ham) #environments for h
    (ham2,envs2) = squaredenvs(init,ham,envs1) #environments for h^2
    mixedenvs = environments(init,A); #environments for <init | A>

    delta=2*tol

    numit = 0
    while delta>tol && numit<maxiter
        numit+=1
        delta=0

        for i in [1:(length(A)-1);length(A):-1:2]

            @plansor tos[-1 -2;-3] := leftenv(mixedenvs,i,init)[-1;1]*A.AC[i][1 -2;2]*rightenv(mixedenvs,i,init)[2;-3]

            (res,convhist) = @closure linsolve(-eta*tos,init.AC[i],GMRES(tol=solvtol)) do x
                y=(eta*eta+w*w)*x
                y-=2*w*ac_prime(x,i,init,envs1)
                y+=ac_prime(x,i,init,envs2)
            end

            delta = max(delta,norm(res-init.AC[i]))
            init.AC[i] = res

            convhist.converged == 0 && @info "($(i)) failed to converge $(convhist.normres)"
        end

        verbose && @info "ddmrg sweep delta : $(delta)"
    end

    a = @plansor leftenv(mixedenvs,1,init)[-1;1]*A.AC[1][1 -2;2]*rightenv(mixedenvs,1,init)[2;-3]*conj(init.AC[1][-1 -2;-3])
    a = a';

    cb = leftenv(envs1,1,A);
    for i in 1:length(A)
        cb = transfer_left(cb,ham[i],init.AL[i],A.AL[i]);
    end

    b = 0*a
    for i in 1:length(cb)
        b+= @plansor cb[i][1 2;3]*A.CR[end][3;4]*rightenv(envs1,length(A),A)[i][4 2;5]*conj(init.CR[end][1;5]);
    end

    v = b/eta-w/eta*a+1im*a
    return v,init
end

function squaredenvs(state::Union{MPSComoving,FiniteMPS},ham::MPOHamiltonian,envs=environments(state,ham))
    nham = conj(ham)*ham

    #to construct the squared caches we will first initialize environments
    #then make all data invalid so it will be recalculated
    #then initialize the right caches at the edge
    ncocache = environments(state,nham)

    #make sure the dependencies are incorrect, so data will be recalculated
    for i in 1:length(state)
        poison!(ncocache,i)
    end

    #impose the correct boundary conditions (important for comoving mps, should do nothing for finite mps)
    indmap = LinearIndices((ham.odim,ham.odim));

    nleft = leftenv(ncocache,1,state)
    nright = rightenv(ncocache,length(state),state)

    for i in 1:ham.odim,j in 1:ham.odim
        f1 = isomorphism(space(nleft[indmap[i,j]],2),space(leftenv(envs,1,state)[i],2)'*space(leftenv(envs,1,state)[j],2))
        @plansor nleft[indmap[i,j]][-1 -2;-3]:=leftenv(envs,1,state)[j][1 3;-3]*conj(leftenv(envs,1,state)[i][1 2;-1])*f1[-2;2 3]

        f2 = isomorphism(space(nright[indmap[i,j]],2),space(rightenv(envs,length(state),state)[j],2)*space(rightenv(envs,length(state),state)[i],2)');
        @plansor nright[indmap[i,j]][-1 -2;-3]:=rightenv(envs,length(state),state)[j][-1 2;1]*conj(rightenv(envs,length(state),state)[i][-3 3;1])*f2[-2;2 3]
    end

    return nham,ncocache
end


text:
#=
UI_main.jl
    Written by Matthew White
    11/14/2019
=#
function runUI()
    println("\33[2J") # clear command window
    commandHist = [] # User command history log
    pdef = ConstDef(0) # pre-define kinematics
    udef = ConstDef(1)
    alphadef = ConstDef(0)
    alphadef = ConstDef(0)
    curfield = UnsteadyFlowSolvers.TwoDFlowField()
    delvort = UnsteadyFlowSolvers.delNone()
    nsteps = 0 # pre-define global vars
    trange = collect(1:10)
    lockTime = false
    mat = []
    frames = []

    ## Dictonary of Cmdmands
    # Stage 1 commands
    list1 = ["airfoil","ndiv","lespcrit","alphadef","oper","pvt","runtime","pdef","veldef","dt","rho","cambercalc","motionfile"]
    help1 = ["Define airfoil designation: coordinate file or FlatPlate","Define airfoil panel division count","Define critical LESP value for LEV shedding","Define pitch motion parameters","Change to operation phase","Define airfoil pivot location","Define simulation run time","Define plunge motion parameters","Define velocity profile parameters","Define time step","Define operating density","Define camber calculation method","Use a file to describe the motion"]
    help1 = hcat(list1,help1)
    help1[1,1] = "load" # overwrite to show command
    help1 = sortslices(help1,dims = 1)

    # Stage 2 commands
    list2 = ["wflag","sflag","wcount","run","save","anim","animstep","plot","vortplot","infoplot"]
    help2 = ["If enabled, writes out iteration files to default save location","If enabled, start simulation from previous output files","Define number of output files","Run a simulation","Save the polar data to a file","If enabled, shows an animation of polars and motion during LVE simulation","Define step count between animation frames","Plot simulation output data","Plot vortex location data for each write step (Must have wflag enabled)","Plot step parameters for each write step (Must have wflag enabled)"]
    help2 = sortslices(hcat(list2,help2),dims = 1)

    commandDict = commands(list1,list2)
    commandDict = setDefaults(commandDict)

    stage = 1 # set preliminary stage

    # Run 'HELP' on startup
    helpMe(stage,help1,help2)

    while true # Loop forever or unitl user exit
        if stage == 1 # Set stage prompts
            prompt = "UNSflow   c> "
            cmds = commandDict.s1
        elseif stage == 2
            prompt = "OPER      c> "
            cmds = commandDict.s2
        end
        ## User Input
        print(prompt)
        usrCmd = ask()
        commandHist = vcat(commandHist,usrCmd)

        ## check for multiple commands in one line
        if occursin(" ",usrCmd) && length(usrCmd) > 2 # sub command
            temp = split(usrCmd," ")
            if length(temp) <= 2 # extra command added
                usrCmd = lowercase(temp[1])
                subCmd = temp[2]
            else
                errorHandler(usrCmd,"many")
                continue
            end
        else # no sub command
            usrCmd = lowercase(usrCmd)
            subCmd = ""
        end

        ## Top level Commands
        if stage == 1 && usrCmd == "exit"
            close("all")
            dirvec = readdir() # Clear step files if prompted
            if "Step Files" in dirvec
                if ynAsk("Clear step file cache?")
                    cleanWrite()
                end
            end
            break # exit program
        end
        if stage == 2 && (usrCmd == "" || usrCmd == "back")
            stage = 1
            close("all")
            continue # return to stage 1
        end
        if usrCmd == "?"
            cmdValues(commandDict,stage)
            continue # show variables
        end
        if usrCmd == "clc"
            println("\33[2J") # clear command window
            continue # clear screen
        end
        if usrCmd == "help"
            helpMe(stage,help1,help2)
            continue
        end

        if stage == 1 && usrCmd == "load" # change from load command
            usrCmd = "airfoil"
        end
        if stage == 2 && (usrCmd == "run") # Calculate surf and curfield
            global surf = UnsteadyFlowSolvers.TwoDSurf(commandDict.s1["airfoil"],commandDict.s1["pvt"],kinem,[commandDict.s1["lespcrit"];]; ndiv = commandDict.s1["ndiv"], camberType = commandDict.s1["cambercalc"], rho = commandDict.s1["rho"])
            curfield = UnsteadyFlowSolvers.TwoDFlowField()
        end

        ## General Commands
        if haskey(cmds,usrCmd) # Test if command is available
            ## Evaluate Cmdmands
            if stage == 1 # Initilization Stage
                if usrCmd == "airfoil" # load airfoil
                    o = loadGeo(subCmd)
                    if o != nothing
                        commandDict.s1[usrCmd] = o
                        makeInitPlot(commandDict.s1,trange,alphadef,pdef,udef)
                    end

                elseif usrCmd == "ndiv" # Set panel divsions
                    o = svPrompt("i", "Panel Count: ", subCmd)
                    if o != nothing
                        commandDict.s1[usrCmd] = o
                        makeInitPlot(commandDict.s1,trange,alphadef,pdef,udef)
                    end

                elseif usrCmd == "lespcrit" # Set panel divsions
                    o = svPrompt("f", "LESP Critical: ", subCmd)
                    if o != nothing ; commandDict.s1[usrCmd] = o end

                elseif usrCmd == "pvt" # Set pivot location
                    o = svPrompt("f", "Pivot Point: ", subCmd)
                    if o != nothing ; commandDict.s1[usrCmd] = o end

                elseif usrCmd == "dt" # Set time step
                    if !lockTime
                        o = svPrompt("f", "Time Step: ", subCmd)
                        if o != nothing
                            commandDict.s1[usrCmd] = o
                            if commandDict.s1["runtime"] != "---"
                                nsteps, trange = iterCalc(commandDict.s1["runtime"],o)
                                println(" The simulation will take $nsteps iterations.")
                            end
                        end
                    else
                        errorHandler(1,"The file '$(commandDict.s1["motionfile"])' prevents this from being changed.")
                    end

                elseif usrCmd == "runtime" # Set total time
                    if !lockTime
                         o = svPrompt("f", "Total Time: ", subCmd)
                         if o != nothing
                             commandDict.s1[usrCmd] = o
                             nsteps, trange = iterCalc(o,commandDict.s1["dt"])
                             println(" The simulation will take $nsteps iterations.")
                         end
                     else
                         errorHandler(1,"The file '$(commandDict.s1["motionfile"])' prevents this from being changed.")
                    end

                elseif usrCmd == "rho" # Set density
                    o = svPrompt("f", "Density [kg/m^3]: ", subCmd)
                    if o != nothing ; commandDict.s1[usrCmd] = o end

                elseif usrCmd == "cambercalc" # Set camber calculation method
                    o = setCamberCalc(subCmd)
                    if o != nothing
                        commandDict.s1[usrCmd] = o
                        makeInitPlot(commandDict.s1,trange,alphadef,pdef,udef)
                    end

                elseif usrCmd == "motionfile" # load time, alphadef, h, and vel from file
                    t,a,h,u,fn = fileInput(subCmd)
                    if t != nothing
                        commandDict.s1["motionfile"] = fn

                        alphadef = FileDef(t,a)
                        commandDict.s1["alphadef"] = "file"

                        pdef = FileDef(t,h)
                        commandDict.s1["pdef"] = "file"

                        udef = FileDef(t,u)
                        commandDict.s1["veldef"] = "file"

                        commandDict.s1["dt"] = t[2]-t[1]
                        commandDict.s1["runtime"] = t[end]
                        lockTime = true # might need a way to toggle off if file is completely overrided

                        trange = t
                        makeInitPlot(commandDict.s1,trange,alphadef,pdef,udef)
                    end

                elseif usrCmd == "alphadef" # Set alpha definition
                    alphavec = modef(subCmd)
                    if alphavec != nothing
                        alphadef = alphaCorr(alphavec)
                        param = join(alphavec[2:end], ", ")
                        commandDict.s1[usrCmd] = "$(alphavec[1])($param)"
                        rt = endTime(alphadef)
                        if rt != "---"
                            if  commandDict.s1["runtime"] == "---" || rt > commandDict.s1["runtime"] && !lockTime
                                commandDict.s1["runtime"] = rt
                                nsteps, trange = iterCalc(commandDict.s1["runtime"],commandDict.s1["dt"])
                            end
                        end
                        makeInitPlot(commandDict.s1,trange,alphadef,pdef,udef)
                    end

                elseif usrCmd == "pdef" || usrCmd == "veldef"# Set height and velocity definition
                    o = modef(subCmd)
                    if o != nothing
                        if o[1] == "const"
                            odef = ConstDef(o[2])
                        elseif o[1] == "lin"
                            odef = LinearDef(o[2],o[3],o[4],o[5])
                        elseif o[1] == "eld"
                            odef = EldRampReturntstartDef(o[2],o[3],o[4],o[5])
                        elseif o[1] == "sin"
                            odef = SinDef(o[2],o[3],o[4],o[5])
                        elseif o[1] == "cos"
                            odef = CosDef(o[2],o[3],o[4],o[5])
                        elseif o[1] == "gust"
                            odef = StepGustDef(o[2],o[3],o[4])
                        end
                        param = join(o[2:end], ", ")
                        commandDict.s1[usrCmd] = "$(o[1])($param)"

                        rt = endTime(odef)

                        if rt != "---"
                            if commandDict.s1["runtime"] == "---" || rt > commandDict.s1["runtime"] && !lockTime
                                commandDict.s1["runtime"] = rt
                                nsteps, trange = iterCalc(commandDict.s1["runtime"],commandDict.s1["dt"])
                            end
                        end

                        if usrCmd == "pdef"
                            pdef = odef
                        else
                            udef = odef
                        end
                        makeInitPlot(commandDict.s1,trange,alphadef,pdef,udef)
                    end

                elseif usrCmd == "oper" # Change to stage 2
                    o = gotoOper(commandDict.s1)
                    if o
                        stage = 2
                        global kinem = KinemDef(alphadef, pdef, udef)
                        nsteps = Int(ceil(commandDict.s1["runtime"]/commandDict.s1["dt"]))
                    end
                    close()
                end

            elseif stage == 2 # Operation Stage
                if usrCmd == "wcount" # Set step output file count
                    o = svPrompt("i","Step output file count: ", subCmd)
                    if o != nothing ; commandDict.s2[usrCmd] = o end

                elseif usrCmd == "wflag" # Write step output during sim
                    commandDict.s2[usrCmd] = ynAsk("Write step output?",subCmd)

                elseif usrCmd == "sflag" # Start from previous step output
                    commandDict.s2[usrCmd] = ynAsk("Start from previous step output?",subCmd)

                elseif usrCmd == "anim" # Set LVE animation state
                    commandDict.s2[usrCmd] = ynAsk("Display LVE Animation?", subCmd)

                elseif usrCmd == "animstep" # set animation step
                    o = svPrompt("i", "Animation Step: ", subCmd)
                    if o != nothing ; commandDict.s2[usrCmd] = o end

                elseif usrCmd == "run" # Run simulation
                    if commandDict.s2["wflag"]
                        cleanWrite()
                    end
                    o = runSim(commandDict,surf,curfield,nsteps,delvort,subCmd)
                    if o != nothing ; mat = o end

                elseif usrCmd == "plot"
                    plotResults(mat,subCmd)

                elseif usrCmd == "save" # Set step output file count
                    if mat != [] # mat has been defined
                        saveMat(subCmd,mat)
                    else
                        errorHandler(mat,"A simulation must be ran before saving the results.")
                    end

                elseif usrCmd == "infoplot"
                    if mat != []
                        tt = plotTime(subCmd)
                        makeKinemClVortPlots2D(mat,tt)
                    else
                        errorHandler(1,"A simulation must be ran before infoplots can be generated.")
                    end

                elseif usrCmd == "vortplot"
                    tt = plotTime(subCmd)
                    makeVortPlots2D(tt)

                end
            end
        elseif usrCmd != "" # return error message
            errorHandler(usrCmd)
            continue
        end

    end
    #return commandDict
end


text:
function solve(tree::NDTree, x::Vector{Float64})
    y = x[tree.p]
    maxLevel = length(tree.t)
    # Forward
    for l = 1:maxLevel
        for s = 1:2^(maxLevel-l)
            Lss = tree.t[l][s].APiv[1,1]
            p   = get_dofs(tree.t[l][s])
            # Pivot
            BLAS.trsv!('L', 'N', 'U', Lss, view(y, p)) # y[p] <- Lss^-1 y[p]
            # Panel
            for (Lnp, n) in zip(tree.t[l][s].ALow[:,1], tree.t[l][s].NbrRanges)
                BLAS.gemv!('N', -1.0, Lnp, view(y, p), 1.0, view(y, n)) # y[n] <- y[n] - Lnp * y[p]
            end
        end
    end
    # Backward
    for l = maxLevel:-1:1
        for s = 1:2^(maxLevel-l)
            Uss = tree.t[l][s].APiv[1,1]
            p   = get_dofs(tree.t[l][s])
            # Panel
            for (Upn, n) in zip(tree.t[l][s].AUpp[1,:], tree.t[l][s].NbrRanges)
                BLAS.gemv!('N', -1.0, Upn, view(y, n), 1.0, view(y, p)) # y[p] <- y[p] - Upn * y[n]
            end
            # Pivot
            BLAS.trsv!('U', 'N', 'N', Uss, view(y, p)) # y[p] <- Uss^-1 y[p]
        end
    end
    return y[invperm(tree.p)]
end


text:
#=

  Euler #29 in Julia.

  Problem 29:
  """
  Consider all integer combinations of a^b for 2 <= a <= 5 and 2 <= b <= 5:

      2^2=4, 2^3=8, 2^4=16, 2^5=32
      3^2=9, 3^3=27, 3^4=81, 3^5=243
      4^2=16, 4^3=64, 4^4=256, 4^5=1024
      5^2=25, 5^3=125, 5^4=625, 5^5=3125

  If they are then placed in numerical order, with any repeats removed, we get the
  following sequence of 15 distinct terms:

  4, 8, 9, 16, 25, 27, 32, 64, 81, 125, 243, 256, 625, 1024, 3125

  How many distinct terms are in the sequence generated by a^b for
  2 <= a <= 100 and 2 <= b <= 100?
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

# 0.00826962s
function euler29a()
    min = 2
    max = 100
    h = Dict()
    for a = min:max
        for b = min:max
            h[BigInt(a)^BigInt(b)] = 1
        end
    end

    return keys(h)|>length
end

# using Set
# 0.20847729s
function euler29b()
    s = (2:100).|>a->(2:100).|>b->BigInt(a)^BigInt(b)
    return Set(vcat(s...))|>length
end

# Using Dict instead
# 0.23657872s
function euler29c()
    s = (2:100).|>a->(2:100).|>b->BigInt(a)^BigInt(b)=>1
    return Dict(vcat(s...))|>length
end

# run_euler(euler29a)
run_euler(euler29b)
# run_euler(euler29c)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
@reexport module AccessControl

export AccessControlEnv

using ReinforcementLearningBase

using Distributions

const N_SERVERS = 10
const PRIORITIES = [1.0, 2.0, 4.0, 8.0]
const FREE_PROB = 0.06
const ACTIONS = (:accept, :reject)
const CUSTOMERS = 1:length(PRIORITIES)

const TRANSFORMER = LinearIndices((0:N_SERVERS, CUSTOMERS))

Base.@kwdef mutable struct AccessControlEnv <: AbstractEnv
    n_servers::Int = 10
    n_free_servers::Int = 0
    customer::Int = rand(CUSTOMERS)
    reward::Float64 = 0.0
end

RLBase.get_observation_space(env::AccessControlEnv) = DiscreteSpace(length(TRANSFORMER))
RLBase.get_action_space(env::AccessControlEnv) = DiscreteSpace(2)

function (env::AccessControlEnv)(a)
    action, reward = ACTIONS[a], 0.0
    if env.n_free_servers > 0 && action == :accept
        env.n_free_servers -= 1
        reward = PRIORITIES[env.customer]
    end

    env.n_free_servers += rand(Binomial(env.n_servers - env.n_free_servers, FREE_PROB))
    env.customer = rand(CUSTOMERS)
    env.reward = reward

    nothing
end

RLBase.observe(env::AccessControlEnv) =
    (
        reward = env.reward,
        terminal = false,
        state = TRANSFORMER[CartesianIndex(env.n_free_servers + 1, env.customer)],
    )

function RLBase.reset!(env::AccessControlEnv)
    env.n_free_servers = env.n_servers
    env.customer = rand(CUSTOMERS)
    env.reward = 0.0
    nothing
end

end


text:
include("tst0.jl")
include("tst1.jl")
include("tst2.jl")
include("tst3.jl")


text:
#
# Copyright (c) 2014-2016, imec
# All rights reserved.
#

using MAT, Distributions, MatrixMarket

macro default_value(sym, value)
  quote
    symsym = $(QuoteNode(sym))
    v = $value
    if ! isdefined(symsym)
      $(esc(sym)) = v
      println("Setting default value for $symsym = $v")
    end
  end
end

@default_value(num_latent, 10)
@default_value(alpha, 2) # observation noise (precision)
@default_value(skip_load, false)
@default_value(nsims, 10)
@default_value(burnin, 5)
@default_value(clamp_lo, 1)
@default_value(clamp_hi, 5)

if length(ARGS) < 2
  println("Usage bpmf.jl <matrix_train.mtx> <matrix_test.mtx>")
  exit()
end

if ! skip_load
  println("Loading data")
  println(ARGS[1])
  Am = MatrixMarket.mmread(ARGS[1])
  P = MatrixMarket.mmread(ARGS[2])
  (I, J, V) = findnz(P)
  P = sparse(I, J, V)
  num_p, num_m = size(Am)
  Au = Am'
  mean_rating = sum(Am)/nnz(Am)
end

sample_u = zeros(num_p, num_latent)
sample_m = zeros(num_m, num_latent)

# Initialize hierarchical priors
mu_u = zeros(num_latent,1)
mu_m = zeros(num_latent,1)
Lambda_u = eye(num_latent)
Lambda_m = eye(num_latent)

# parameters of Inv-Whishart distribution (see paper for details)
WI_u = eye(num_latent)
b0_u = 2
df_u = num_latent
mu0_u = zeros(num_latent,1)

WI_m = eye(num_latent)
b0_m = 2
df_m = num_latent
mu0_m = zeros(num_latent,1)

function pred(P, sample_m, sample_u, mean_rating)
  (I, J, Vin) = findnz(P)
  pr = sum(sample_m[J,:].*sample_u[I,:],2) + mean_rating;
  Vout = clamp(pr[:,1], clamp_lo, clamp_hi)
  sparse(I, J, Vout)
end

function ConditionalNormalWishart(U::Matrix{Float64}, mu::Vector{Float64}, kappa::Real, T::Matrix{Float64}, nu::Real)
  N = size(U, 1)
  Ū = mean(U,1)
  S = cov(U, mean=Ū)
  Ū = Ū'

  mu_c = (kappa*mu + N*Ū) / (kappa + N)
  kappa_c = kappa + N
  T_c = inv( inv(T) + N * S + (kappa * N)/(kappa + N) * (mu - Ū) * (mu - Ū)' )
  nu_c = nu + N

  NormalWishart(vec(mu_c), kappa_c, T_c, nu_c)
end

function grab_col{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, col::Integer)
  r = A.colptr[col]:A.colptr[col+1]-1
  A.rowval[r], A.nzval[r]
end

function sample_movie(mm, Am, mean_rating, sample_u, alpha, mu_m, Lambda_m)
  ff, v = grab_col(Am, mm)
  rr = v - mean_rating
  MM = sample_u[ff,:]

  # Sample from normal distribution
  L = chol(Lambda_m + alpha * MM'*MM, :L)
  L' \ (randn(num_latent) + (L \ (alpha * MM'*rr + Lambda_m * mu_m)))
end

function sample_user(uu, Au, mean_rating, sample_m, alpha, mu_u, Lambda_u)
  ff, v = grab_col(Au, uu)
  rr = v - mean_rating
  MM = sample_m[ff,:]

  # Sample from normal distribution
  L = chol(Lambda_u + alpha * MM'*MM, :L)
  L' \ (randn(num_latent) + (L \ (alpha * MM'*rr + Lambda_u * mu_u)))
end

tic()
println("Sampling")
for i in 1:nsims
  # Sample from movie hyperparams
  mu_m, Lambda_m = rand( ConditionalNormalWishart(sample_m, vec(mu0_m), b0_m, WI_m, df_m) )

  # Sample from user hyperparams
  mu_u, Lambda_u = rand( ConditionalNormalWishart(sample_u, vec(mu0_u), b0_u, WI_u, df_u) )

  for mm = 1:num_m
    sample_m[mm, :] = sample_movie(mm, Am, mean_rating, sample_u, alpha, mu_m, Lambda_m)
  end

  for uu = 1:num_p
    sample_u[uu, :] = sample_user(uu, Au, mean_rating, sample_m, alpha, mu_u, Lambda_u)
  end

  probe_rat = pred(P, sample_m, sample_u, mean_rating)

  if i > burnin
    probe_rat_all = (counter_prob*probe_rat_all + probe_rat)/(counter_prob+1)
    counter_prob = counter_prob + 1
  else
    probe_rat_all = probe_rat
    counter_prob = 1
  end

  err_avg = sqrt( sum((P .- probe_rat_all).^2) / nnz(P) );
  err = sqrt( sum((P .- probe_rat).^2) / nnz(P) );

  @printf("Iteration %d:\t avg RMSE %6.4f RMSE %6.4f FU(%6.4f) FM(%6.4f)\n", i, err_avg, err, vecnorm(sample_u), vecnorm(sample_m))
end
toc()

max_stars_repo_name: IoanAnev/BPMF-ArgoDSM

max_issues_repo_name: IoanAnev/argodsm-bpmf

max_forks_repo_name: IoanAnev/argodsm-bpmf


text:
using Reexport

export @dependencies

macro dependencies()
  :(using Genie, SearchLight, Helpers, ControllerHelper, Renderer, Cache, Router, Util, Logger, Macros, Flax)
end

max_stars_repo_name: essenciary/genie-demo-tweet-stats

max_issues_repo_name: essenciary/genie-demo-tweet-stats

max_forks_repo_name: essenciary/genie-demo-tweet-stats


text:
export PixelPicture, draw, randomize!, set_point!, get_points

"""
`point2pixel(n,xy)` determines which pixel in an `n`-by-`n`
grid corresponds to the point with coordinates given by `xy`.
Retuns a tuple.
"""
function point2pixel{T<:Real}(n::Int, xy::Vector{T})
  x = xy[1]
  y = xy[2]
  i = Int(round((n-1)*(1-y)+1))
  j = Int(round((n-1)*x+1))
  return (i,j)
end

"""
`pixel2point(n,ij)` determines the `xy`-coordinates of the point
in the middle of the `i,j`-pixel in an `n`-by-`n` grid. Returns
a vector.
"""
function pixel2point(n::Int,ij::Tuple{Int,Int})
  i,j = ij
  x = (j-1)/(n-1)
  y = 1 - (i-1)/(n-1)
  return [x,y]
end

"""
`PixelPicture` is a square grid of pixels representing a picture
in the unit square. Create a new one with `PixelPicture(n)`.
"""
type PixelPicture
  data::Matrix{Bool}
  function PixelPicture(n::Int)
    new(zeros(Bool,n,n))
  end
end

"""
`size(P)` returns the size of a `PixelPicture`. Since these
images are square, we return a single number.
"""
function size(P::PixelPicture)
  (r,c) = size(P.data)
  return r
end

function show(io::IO,P::PixelPicture)
  n = size(P)
  print(io,"$n-by-$n PixelPicture")
end

"""
`get_points(P)` returns the set of x,y-vectors in the unit square
in this `PixelPicture`.
"""
function get_points(P::PixelPicture)
  n = size(P)
  pix_coords = ( (i,j) for i=1:n for j=1:n if P.data[i,j] )
  result = ( pixel2point(n,ij) for ij in pix_coords )
  return result
end

"""
`randomize!(P)` fills the `PixelPicture` with random noise.
"""
function randomize!(P::PixelPicture)
  n = size(P)
  P.data = rand(Bool,n,n)
  nothing
end

function set_point!(P::PixelPicture,xy::Array{Float64,1})
  n = size(P)
  i,j = point2pixel(n,xy)
  if 1<=i<=n && 1<=j<=n
    @inbounds P.data[i,j]=true
  end
  nothing
end

"""
`draw(P)` draws the `PixelPicture` on the screen.
"""
function draw(P::PixelPicture)
  clf()
  spy(P.data)
  axis("off")
  nothing
end


text:
"""
    module Lattice

for defining discrete Cartesian and Polar lattice grids used in Iris.
"""
module Lattices

# using ...Defaults
using ..Points
using Formatting
# using LinearAlgebra
# using RecipesBase
using StaticArrays

export Lattice
export latticeindex

struct Lattice{N,M,P}
    constants::NTuple{N,Float64}
    e::NTuple{N,Point{N}}
    type::Symbol
    origin::Point{N}
    angles::SVector{M,Float64}
    r0::Float64 # radius of first lattice point for Polar/Spherical
    R::SMatrix{N,N,Float64,P}

    function Lattice(constants::NTuple{N,Float64},
                type::Symbol,
                origin::Point{N},
                angles::SVector{M,Float64},
                r0::Float64) where {N,M}

        1≤N≤3 || throw(ErrorException("lattice only defined for dimensions ≤ 3, d=$N"))
        M==(N-1)N÷2 || throw(ErrorException("incorrect number of angles for dimension, should be $((N-1)N÷2), but is $M"))
        e,R = lattice_primitives(angles...,constants)
        return new{N,M,N*N}(constants,e,type,origin,angles,r0,R)
    end
end
Lattice(dx::Tuple; kwargs...) = Lattice(dx...; kwargs...)
Lattice(dx...; kwargs...) = Lattice(float.(dx); kwargs...)
function Lattice(constants::NTuple{N,Float64};
            type::Symbol = :Cartesian,
            origin = Point(ntuple(i->0.0,N)),
            angles = SVector{(N-1)N÷2,Float64}(ntuple(i->0.0,(N-1)N÷2)),
            r0::Real = constants[1]/2
            ) where N

    type = _lattice_type(type)
    origin = _lattice_origin(origin,constants)
    n = length(origin)
    (n-1)n÷2 == length(angles) || throw(ErrorException("incorrect number of angles for dimension, should be $((n-1)n÷2), but is $(length(angles))"))
    angles = _lattice_angles(angles,constants)
    return Lattice(constants,type,Point(origin),angles,r0)
end

function _lattice_type(type::Symbol)
    if type ∈ (:r,:R,:rectangular, :Rectangular, :rect, :Rect, :rectilinear, :Rectilinear, :cart, :Cart, :cartesian, :Cartesian)
        return :Cartesian
    elseif type ∈ (:p,:P,:polar, :Polar, :cyl, :Cyl, :cylindrical, :Cylindrical,:pol,:Pol)
        return :Polar
    elseif type ∈ (:s,:S,:spherical, :Spherical, :sph, :Sph)
        return :Spherical
    else
        throw("unrecognized coordinate type $type")
    end
end
_lattice_origin(origin,constants::NTuple{N}) where N = Point(origin)
_lattice_angles(angles,constants::NTuple{N}) where N = SVector{(N-1)N÷2,Float64}(angles)

include("1D/Lattices.jl")
# include("2D/Lattices.jl")
# include("3D/Lattices.jl")

function Base.getproperty(lat::Lattice,sym::Symbol)
    if sym == :dx
        return getfield(lat,:constants)[1]
    elseif sym == :dy
        return getfield(lat,:constants)[2]
    elseif sym == :dz
        return getfield(lat,:constants)[3]
    elseif sym == :dr
        return getfield(lat,:constants)[1]
    elseif sym == :x0
        return getfield(lat,:origin)[1]
    elseif sym == :y0
        return getfield(lat,:origin)[2]
    elseif sym == :z0
        return getfield(lat,:origin)[3]
    elseif sym == :e1
        return getfield(lat,:e)[1]
    elseif sym == :e2
        return getfield(lat,:e)[2]
    elseif sym == :e3
        return getfield(lat,:e)[3]
    elseif Base.sym_in(sym,(:dϕ,:dφ,:dphi))
        return getfield(lat,:constants)[2]
    elseif Base.sym_in(sym,(:dθ,:dϑ,:dtheta))
        return getfield(lat,:constants)[3]
    elseif Base.sym_in(sym,(:ϕ,:φ,:phi))
        return getfield(lat,:angles)[1]
    elseif Base.sym_in(sym,(:θ,:ϑ,:theta))
        return getfield(lat,:angles)[1]
    else
        return getfield(lat,sym)
    end
end

latticeindex(lat::Lattice,x...) = latticeindex(lat,Point(x...))

end # module