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

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

# ╔═╡ 393b2f5e-6556-11eb-2119-cf7309ee7392
begin
    import Pkg
    # activate the shared project environment
    Pkg.activate(Base.current_project())

	using Chain
	using CSV
	using DataFrames
	using HTTP
	using LsqFit
	using Measurements
	using Plots
	using PlutoUI
	using Statistics
	using Unitful
	using UnitfulRecipes
	using URIs
end

# ╔═╡ 412eedac-658a-11eb-2326-93e5bf3d1a2c
md"""
# Fitting of equilibrium binding data

This notebook plots an equilibrium binding dataset (with error bars, for datasets containing replicates) and performs non-linear curve fitting of a model to the data.

The following reference explains very well the theory of equilibrium binding experiments, as well as many important practical considerations:

> Jarmoskaite I, AlSadhan I, Vaidyanathan PP & Herschlag D (2020) How to measure and evaluate binding affinities. *eLife* **9**: e57264 <https://doi.org/10.7554/eLife.57264>
"""

# ╔═╡ 270cc0cc-660f-11eb-241e-b75746a39cc7
md"""
## Load data

The data file must be in CSV format. The first row is assumed to contain column names. The first column is assumed to be the $X$ values, all other columns are assumed to be replicate $Y$ values that will be averaged (fitting will be done against the mean values, weighted by their standard deviations). In addition, there must not be any row with an $X = 0$ value (this would result in an error when attempting to plot with a logarithmic scale). I always perform two measurements at $X = 0$, and I always sort rows by descending $X$ values, so this notebook automatically skips the last two rows of the CSV file; adjust accordingly if you don't measure at $X = 0$ and want to keep all rows (see section [Data processing](#1884912a-6aeb-11eb-2b4a-d14d4a321dc5) below). The data does *not* need to be scaled such that $Y$ takes values between $0$ and $1$: the binding models can account for arbitrary minimum and maximum $Y$ values (see section [Model functions](#663a4cae-658a-11eb-382f-cf256c08c9d1) below). Scaling data will hide differences in signal change between datasets, while these differences may tell you something about the system under study, so scaling should never be done "blindly"; always look at the raw data.

Indicate below which data files to process:

- if the path given is not absolute, it is assumed to be relative to the notebook file (wherever the notebook is located)
- a list of files can be provided with one file path per line and separated by commas
- files can be located by local path or URL, but each type of location should be in the dedicated list
"""

# ╔═╡ a2c02fcf-9402-4938-bc3d-819b45a66afa
dataURLs = [
	"https://raw.githubusercontent.com/Guillawme/julia-curve-fitting/main/datasets/dataset_006.csv",
	"https://raw.githubusercontent.com/Guillawme/julia-curve-fitting/main/datasets/dataset_007.csv"
]

# ╔═╡ 14baf100-660f-11eb-1380-ebf4a860eed8
dataFiles = [
	"datasets/dataset_008.csv",
	"datasets/dataset_009.csv"
]

# ╔═╡ 2ce72e97-0133-4f15-bf1d-7fd04ccf3102
md"""
**Number of rows to ignore at the end of files:** $(@bind footerRows PlutoUI.NumberField(0:30, default = 2))
"""

# ╔═╡ 77316a92-425a-4902-9828-52a7c4a74f27
md"""
**Unit of your concentration values ($X$ axis):** $(
@bind chosenConcUnit PlutoUI.Select(
	[
		"pM" => "pM",
		"nM" => "nM",
		"μM" => "μM",
		"mM" => "mM"
	],
	default = "nM"
))
"""

# ╔═╡ 214acce6-6ae8-11eb-3abf-492e50140317
md"""
Your data should appear below shortly, check that it looks normal. In addition to the columns present in your CSV file, you should see four columns named `mean`, `std`, `measurement` and `conc` (these values will be used for fitting and plotting).
"""

# ╔═╡ d5b0e2a1-865c-489c-9d0d-c4ae043828fb
# By defaults, use file names and URLs to identify datasets.
datasetNames = vcat(dataURLs, dataFiles)
# But one can also use custom names.
#datasetNames = ["Monday", "Tuesday", "Wednesday", "Thursday"]

# ╔═╡ 5c5e0392-658a-11eb-35be-3d940d4504cb
md"""
## Visualizations

Your data and fit should appear below shortly. Take a good look at the [data and fit](#3dd72c58-6b2c-11eb-210f-0b547bf38ebe), make sure you check the [residuals](#4f4000b4-6b2c-11eb-015f-d76a0adda0a0). Once you're happy with it, check the [numerical results](#be17b97e-663a-11eb-2158-a381c19ece3f).
"""

# ╔═╡ 3dd72c58-6b2c-11eb-210f-0b547bf38ebe
md"""
### Data and fit

Select binding model:
"""

# ╔═╡ 3da83f72-6a11-11eb-1a74-49b66eb39c96
@bind chosenModel PlutoUI.Radio(
	[
		"Hill" => :Hill,
		"Hyperbolic" => :Hyperbolic,
		"Quadratic" => :Quadratic
	],
	default = "Hill"
)

# ╔═╡ d15cba72-6aeb-11eb-2c80-65702b48e859
md"""
Show fit line with initial parameters?
$@bind showInitialFit PlutoUI.CheckBox(default = false)
"""

# ╔═╡ 01b59d8a-6637-11eb-0da0-8d3e314e23af
md"""
For the quadratic model, indicate receptor concentration (the receptor is the binding partner kept at constant, low concentration across the titration series).
Parameter $R_0 =$
$@bind R0 PlutoUI.Slider(0.01:0.1:500.0, default = 5.0, show_value = true)
"""

# ╔═╡ 4f4000b4-6b2c-11eb-015f-d76a0adda0a0
md"""
### Residuals
"""

# ╔═╡ c50cf18c-6b11-11eb-07d3-0b8e332ec5bc
md"""
The fit residuals should follow a random normal distribution around $0$. If they show a systematic trend, it means the fit systematically deviates from your data, and therefore the model you chose might not be justified (but be careful when considering alternative models: introducing more free parameters will likely get the fit line closer to the data points and yield a lower [sum of squared residuals](#124c4f94-6b99-11eb-2921-d7c2cd00b893), but this is not helpful if these additional parameters don't contribute to explaining the physical phenomenon being modeled). Another possibility is a problem with your data. The most common problems are:

- the data does not cover the proper concentration range
- the concentration of receptor is too high relative to the $K_D$

In either case, your best option is to design a new experiment and collect new data.
"""

# ╔═╡ 5a36fc3f-ce74-42c8-8284-19321e0d687f
md"""
#### Scatter plot
"""

# ╔═╡ 5392d99b-70f9-48cd-90b4-58cba5fc9681
md"""
#### Histogram
"""

# ╔═╡ be17b97e-663a-11eb-2158-a381c19ece3f
md"""
## Numerical results

### Model parameters

The parameter $K_D$ is in unit of $(chosenConcUnit).
"""

# ╔═╡ 124c4f94-6b99-11eb-2921-d7c2cd00b893
md"""
### Sum of squared residuals
"""

# ╔═╡ 7e7a9dc4-6ae8-11eb-128d-83544f01b78b
md"""
## Code

The code doing the actual work is in this section. Do not edit unless you know what you are doing.
"""

# ╔═╡ 512e3028-6ae9-11eb-31b4-1bc9fc66b322
md"### Necessary packages and notebook setup"

# ╔═╡ abc03f64-6a11-11eb-0319-ed7cea455cb5
PlutoUI.TableOfContents()

# ╔═╡ 1884912a-6aeb-11eb-2b4a-d14d4a321dc5
md"### Data processing"

# ╔═╡ c94ef72d-bd12-4434-935e-01e94d5a4588
md"""
#### Concentration unit selection
"""

# ╔═╡ 54272681-dfeb-4034-b127-7e68c19fd576
md"""
First, we need an alias of M (molar) for mol/L, since this is the notation most widely used in the field:
"""

# ╔═╡ 731492c6-95c7-449f-8c19-53e22ab438b8
begin
	Unitful.register(@__MODULE__)
	@unit M "M" Molar 1u"mol/L" true
end

# ╔═╡ 93ef0431-643b-4f6f-8c09-beabff59e0c6
md"""
Check that our alias works:
"""

# ╔═╡ 1c10231d-3bea-4468-a2d8-886c05c6474c
typeof(M)

# ╔═╡ 730b693d-cca8-46de-8382-c151b5f63352
1u"nM" == 1u"nmol/l"

# ╔═╡ 24152a84-523b-4027-9b5a-7e3524b9c659
dimension(1u"nM") == dimension(1u"mol/l")

# ╔═╡ 97f73ee4-3db7-43ba-93eb-17025b485f4f
md"""
This dictionary maps dropdown menu options (in section [Load data](#270cc0cc-660f-11eb-241e-b75746a39cc7) above) to their corresponding unit:
"""

# ╔═╡ 8f2d959e-5e61-48a9-bbf0-538bc1d478a8
concUnits = Dict(
	"pM" => u"pM",
	"nM" => u"nM",
	"μM" => u"μM",
	"mM" => u"mM"
)

# ╔═╡ c866d213-4680-4b08-8ecc-1faefc8661a4
md"""
The following cells simply check which unit is selected in section [Load data](#008f4f8c-9d21-11eb-0fea-f3b2e58957d1) above.
"""

# ╔═╡ 8e9cd30f-1722-4f2f-a26b-2f558805d4a1
chosenConcUnit

# ╔═╡ 13417555-be95-4662-ad11-eca4dece81b5
concUnits[chosenConcUnit]

# ╔═╡ 9c7e922e-c82f-41a4-9513-4462a0559c3f
md"""
#### Processing
"""

# ╔═╡ 4f4b580d-507c-4ad0-b1d5-5967c8ed829e
md"""
The `commonProcessing()` function computes the mean and standard deviation of replicates, defines measurements as mean ± std, and returns a DataFrame containing all the data. It is used by all methods of the following `processData()` function, which handle various inputs (path to a local file, URL to a remote file, loaded CSV file, loaded DataFrame).
"""

# ╔═╡ 5eb607c7-172b-4a6c-a815-acbc195108f0
function commonProcessing(data::DataFrame)
	df = @chain data begin
		# Rename column 1, so we can always call it by the
		# same name regardless of its name in the input file.
		rename(1 => :concentration)
		@aside cols = ncol(_)
		transform(
			# Calculate mean and stddev of replicates
			# (all columns in input except first one).
			AsTable(2:cols) => ByRow(mean) => :mean,
			AsTable(2:cols) => ByRow(std) => :std
		)
		transform(
			# Mean and stddev together define a measurement
			# (this is only for plotting; fitting uses the two
			# original columns separately).
			[:mean, :std] => ByRow(measurement) => :measurement
		)
		transform(
			# Assign unit to concentration column, store as
			# a new column callec conc (only for plotting).
			:concentration => ByRow(x -> x * concUnits[chosenConcUnit]) => :conc
		)
	end
	return df
end

# ╔═╡ 992c86a2-6b13-11eb-1e00-95bdff2736d0
md"""
The `processData()` function loads one data file, computes the mean and standard deviation of replicates, defines measurements as mean ± std, and returns a DataFrame containing all the data.
"""

# ╔═╡ 1fe4d112-6a11-11eb-37a6-bf95fbe032b1
function processData(dataFile::String)
	df = @chain dataFile begin
		CSV.read(footerskip=footerRows, DataFrame)
		commonProcessing()
	end
	return df
end

# ╔═╡ 5cc812af-4f8e-444b-9652-fb063cc6be06
md"""
This functions should also work if passed an URL to a remote file:
"""

# ╔═╡ 7c2af794-a0ed-4827-bc8e-1f15ad205eca
function processData(dataFile::URI)
	df = @chain dataFile begin
		HTTP.get(_).body
		CSV.read(footerskip=footerRows, DataFrame)
		commonProcessing()
	end
	return df
end

# ╔═╡ a08d1b29-fa5e-4312-977c-39b050de4516
md"""
This functions should also work if passed an already loaded CSV file:
"""

# ╔═╡ 0181a0b1-ec6c-4325-8b3a-851a3fe33846
function processData(dataFile::CSV.File)
	df = @chain dataFile begin
		DataFrame()
		commonProcessing()
	end
	return df
end

# ╔═╡ 247c2416-6c67-11eb-01df-8dac01cdbf8f
md"""
This functions should also work if passed an already loaded data frame (for example, if the user wants to load data and pre-process it in a different way before averaging replicates):
"""

# ╔═╡ 36ebe112-6c66-11eb-11f0-7fdb946865e4
function processData(data::DataFrame)
	return commonProcessing(data)
end

# ╔═╡ d904fd76-6af1-11eb-2352-837e03072137
begin
	allData = vcat(URI.(dataURLs), dataFiles)
	dfs = [ processData(df) for df in allData ]
end

# ╔═╡ 0e8af3be-7ae7-4ec2-8d7a-670878cd52ee
md"""
We will need to keep track of dataset names. You can edit them here if you want to use something more meaningful than the file name or URL; changes will propagate to plot legends. This list **must** contain the same number of elements as you have datasets: **$(length(allData))** in the present case.
"""

# ╔═╡ 8e105fae-6aec-11eb-1471-83aebb776241
md"### Plotting"

# ╔═╡ 97c5020c-6aec-11eb-024b-513b1e603d98
md"The `initMainPlot()` function initializes a plot, the `plotOneDataset!()` function plots one dataset (call it repeatedly to plot more datasets on the same axes)."

# ╔═╡ caf4a4a2-6aec-11eb-2765-49d67afa47dd
function initMainPlot()
	plot(
		xlabel = "Concentration",
		ylabel = "Signal",
		xscale = :log10,
		legend = :topleft
	)
end

# ╔═╡ fc194672-6aed-11eb-0a06-2d967ec094b1
md"The `initResidualPlot()` function initializes a plot, the `plotOneResiduals!()` function plots the fit residuals from one dataset (call it repeatedly to plot more datasets on the same axes)."

# ╔═╡ 14db987c-6aee-11eb-06cf-a11987b98f1e
function initResidualPlot()
	plot(
		xlabel = "Concentration",
		ylabel = "Fit residual",
		xscale = :log10,
		legend = :topleft
	)
	hline!([0], label = nothing, color = :red)
end

# ╔═╡ 7f83b838-6a11-11eb-3652-bdff24f3473e
function plotOneResiduals!(plt, df, fit, filePath)
	title = split(filePath, "/")[end]
	scatter!(
		plt,
		df.conc,
		fit.resid,
		label = "$title: $chosenModel fit residual"
	)
end

# ╔═╡ 9020fa5d-7408-4161-a52f-df37b3c2e6f5
md"The `initResidualHistogram()` function initializes a histogram, the `plotOneResidualsHistogram!()` function plots a histogram of the fit residuals from one dataset (call it repeatedly to plot more datasets on the same axes)."

# ╔═╡ 3184d209-1cc9-40ed-a9ec-9f094c5c94b5
function initResidualHistogram()
	histogram(
		xlabel = "Fit residual",
		ylabel = "Count",
		legend = :topleft
	)
	vline!([0], label = nothing, color = :red)
end

# ╔═╡ b38cd229-64e2-4ca4-a78b-8881ec166b09
function plotOneResidualsHistogram!(plt, df, fit, filePath)
	title = split(filePath, "/")[end]
	histogram!(
		plt,
		fit.resid,
		bins = length(fit.resid),
		label = "$title: $chosenModel fit residual"
	)
end

# ╔═╡ 663a4cae-658a-11eb-382f-cf256c08c9d1
md"### Model functions"

# ╔═╡ 594e7534-6aeb-11eb-1254-3b92b71877ed
md"#### Model selection"

# ╔═╡ 32dce844-6aee-11eb-3cf2-3ba420d311d3
md"""
This dictionary maps radio button options (in section [Visualizations](#5c5e0392-658a-11eb-35be-3d940d4504cb) above) to their corresponding model function:
"""

# ╔═╡ a1f56b0a-6aeb-11eb-0a44-556fad58f368
md"""
The remaining cells in this section are only meant to check that the model selection buttons work. This first cell should return the name of the selected binding model (corresponding to the active radio button in section [Visualizations](#5c5e0392-658a-11eb-35be-3d940d4504cb) above):
"""

# ╔═╡ 5be2e5d2-6a11-11eb-1421-492f5af16f9c
chosenModel

# ╔═╡ 58617378-6aee-11eb-23e8-c13d89b4c57f
md"""
This other cell should return the model function corresponding to the selected binding model (the active radio button in section [Visualizations](#5c5e0392-658a-11eb-35be-3d940d4504cb) above):
"""

# ╔═╡ 88d941e6-658a-11eb-08a2-0f021e5ae3a4
md"""
#### Hill model

This is the [Hill equation](https://en.wikipedia.org/wiki/Hill_equation_(biochemistry)):

$S = S_{min} + (S_{max} - S_{min}) \times \frac{L^h}{{K_D}^h + L^h}$

In which $S$ is the measured signal ($Y$ value) at a given value of ligand concentration $L$ ($X$ value), $S_{min}$ and $S_{max}$ are the minimum and maximum values the observed signal can take, respectively, $K_D$ is the equilibrium dissociation constant and $h$ is the Hill coefficient."""

# ╔═╡ e9fc3d44-6559-11eb-2da7-314e8fc76ee9
@. hill(conc, p) = p[1] + (p[2] - p[1]) * conc^p[4] / (p[3]^p[4] + conc^p[4])

# ╔═╡ 9d1f24cc-6a0f-11eb-3b16-35f89aff5d4a
md"""
#### Hyperbolic model

The hyperbolic equation is a special case of the Hill equation, in which $h = 1$:
"""

# ╔═╡ 78e664d0-6618-11eb-135b-5574bb05ddef
@. hyperbolic(conc, p) = p[1] + (p[2] - p[1]) * conc / (p[3] + conc)

# ╔═╡ b0b17206-6a0f-11eb-2f5e-5fc8fa06cd36
md"""
#### Quadratic model

Unlike the Hill and hyperbolic models, the quadratic model does not make the approximation that the concentration of free ligand at equilibrium is equal to the total ligand concentration:

$S = S_{min} + (S_{max} - S_{min}) \times \frac{(K_{D} + R_{tot} + L_{tot}) - \sqrt{(- K_{D} - R_{tot} - L_{tot})^2 - 4 \times R_{tot} \times L_{tot}}}{2 \times R_{tot}}$

Symbols have the same meaning as in the previous equations, except here $L_{tot}$ is the total concentration of ligand, not the concentration of free ligand at equilibrium. $R_{tot}$ is the total concentration of receptor.

In principle, $R_{tot}$ could be left as a free parameter to be determined by the fitting procedure, but in general it is known accurately enough from the experimental set up, and one should replicate the same experiment with different concentrations of receptor to check its effect on the results. $R_{tot}$ should be set in the experiment to be smaller than $K_D$, ideally, or at least of the same order of magnitude than $K_D$. It might take a couple experiments to obtain an estimate of $K_D$ before one can determine an adequately small concentration of receptor at which to perform a definite experiment.
"""

# ╔═╡ 5694f1da-6636-11eb-0fed-9fee5c48b114
@. quadratic(conc, p) = p[1] + (p[2] - p[1]) * ( (p[3] + R0 + conc) - sqrt((- p[3] - R0 - conc) ^ 2 - 4 * R0 * conc) ) / (2 * R0)

# ╔═╡ 605f06ae-6a11-11eb-0b2f-eb81f6526829
bindingModels = Dict(
	"Hill" => hill,
	"Hyperbolic" => hyperbolic,
	"Quadratic" => quadratic
)

# ╔═╡ 7c03fcbe-6a11-11eb-1b7b-cbad863156a6
function plotOneDataset!(plt, df, fit, filePath, showInitial = false, initialValues = nothing)
	title = split(filePath, "/")[end]
	scatter!(
		plt,
		df.conc,
		df.measurement,
		label = "$title: data"
	)
	if showInitial
		plot!(
			plt,
			df.conc,
			bindingModels[chosenModel](df.concentration, initialValues),
			label = "$title: $chosenModel fit (initial)",
			color = :grey
		)
		plot!(
			plt,
			df.conc,
			bindingModels[chosenModel](df.concentration, fit.param),
			label = "$title: $chosenModel fit (converged)",
			color = :red
		)
	else
		plot!(
			plt,
			df.conc,
			bindingModels[chosenModel](df.concentration, fit.param),
			label = "$title: $chosenModel fit",
			color = :red
		)
	end
end

# ╔═╡ 6668c49a-6a11-11eb-2abf-5feecaee8972
bindingModels[chosenModel]

# ╔═╡ 0f960a8a-6a0f-11eb-04e2-b543192f6354
md"### Parameters and their initial values"

# ╔═╡ 2eb890a4-658c-11eb-1fc4-af645d74109d
md"""
The `findInitialValues()` function takes the measured data and returns an array containing initial values for the model parameters (in this order): $S_{min}$, $S_{max}$, $K_D$ and $h$ (for the Hill model only, so the function needs to know which model was selected).

Initial values for $S_{min}$ and $S_{max}$ are simply taken as the minimal and maximal values found in the data. The initial estimate for $K_D$ is the concentration of the data point that has a signal closest to halfway between $S_{min}$ and $S_{max}$ (if the experiment was properly designed, this is a reasonable estimate and close enough to the true value for the fit to converge). The initial estimate of $h$ is $1.0$, meaning we assume no cooperativity.
"""

# ╔═╡ ca2d2f12-6a1a-11eb-13ca-1f93df2b8e4a
function findInitialValues(df, model)
	halfSignal = minimum(df.mean) + (maximum(df.mean) - minimum(df.mean)) / 2
	if model == "Hill"
		params = [
			minimum(df.mean),
			maximum(df.mean),
			df.concentration[findmin(abs.(halfSignal .- df.mean))[2]],
			1.0
		]
	else
		params = [
			minimum(df.mean),
			maximum(df.mean),
			df.concentration[findmin(abs.(halfSignal .- df.mean))[2]]
		]
	end
end

# ╔═╡ 9be41e32-6af0-11eb-0904-d1cf3c288cab
md"Determine initial values of the selected model's parameters from the currently loaded datasets:"

# ╔═╡ 67b538f6-6a1b-11eb-3004-2d89c2f941e8
initialParams = [ findInitialValues(df, chosenModel) for df in dfs ]

# ╔═╡ 213e8ffa-6a0f-11eb-357e-638146193c5d
md"### Fitting"

# ╔═╡ babcb896-6af0-11eb-194a-15922bc2df83
md"Perform fit of the selected model to the measurements' mean values using initial values for the model parameters determined previously. If the dataset contains replicates, the fit will be weighted by the measurements' standard deviations."

# ╔═╡ 47426056-6af2-11eb-17f8-6d27d35003ca
begin
	fits = Vector{LsqFit.LsqFitResult}(undef, length(allData))
	for (df, initialValues, i) in zip(dfs, initialParams, 1:length(allData))
		if ncol(df) > 5
			# If the dataset has more than 5 columns, it means it has
			# replicate Y values, so we weight the fit by their stddev.
			fits[i] = curve_fit(bindingModels[chosenModel],
						df.concentration,
						df.mean,
						df.std,
						initialValues)
		else
			# If the dataset has only 5 columns, it means it doesn't have
			# replicate Y values, so there are no stddev we can use as weights.
			fits[i] = curve_fit(bindingModels[chosenModel],
						df.concentration,
						df.mean,
						initialValues)
		end
	end
	fits
end

# ╔═╡ 264bf9ec-6af5-11eb-1ffd-79fb3466f596
begin
	dataPlot = initMainPlot()
	for (df, fit, title, initialVals) in zip(dfs, fits, datasetNames, initialParams)
		plotOneDataset!(dataPlot, df, fit, title, showInitialFit, initialVals)
	end
	dataPlot
end

# ╔═╡ a951b5dc-6af7-11eb-2401-5d11a14e3067
begin
	residualPlot = initResidualPlot()
	for (df, fit, title) in zip(dfs, fits, datasetNames)
		plotOneResiduals!(residualPlot, df, fit, title)
	end
	residualPlot
end

# ╔═╡ 7625d41a-dab1-4b10-947c-3667c03f85aa
begin
	residualHistogram = initResidualHistogram()
	for (df, fit, title) in zip(dfs, fits, datasetNames)
		plotOneResidualsHistogram!(residualHistogram, df, fit, title)
	end
	residualHistogram
end

# ╔═╡ 2109f516-6b99-11eb-05a0-99b9ecfd0f9d
PlutoUI.with_terminal() do
	println("Dataset\t\t\t\tSum of squared residuals")
	for (dataset, fit) in zip(datasetNames, fits)
		println(
			split(dataset, "/")[end],
			"\t\t",
			round(sum(fit.resid.^2), digits = 2)
		)
	end
end

# ╔═╡ 799680d0-6af1-11eb-321d-b7758a40f931
md"Degrees of freedom:"

# ╔═╡ 1f0384de-659b-11eb-043e-5b86fcdd36e6
dof.(fits)

# ╔═╡ 54501a10-6b9c-11eb-29de-77afc3772fb7
md"Best fit parameters:"

# ╔═╡ 5ed3ab64-6b9c-11eb-149e-43a1ef12ac7d
coef.(fits)

# ╔═╡ 8643b03c-6af1-11eb-0aa7-67acee28d2c0
md"Standard errors of best-fit parameters:"

# ╔═╡ a74998b4-659c-11eb-354d-09ff62710b87
paramsStdErrors = stderror.(fits)

# ╔═╡ 090347fc-6b8e-11eb-0e17-9d9d45749c0b
PlutoUI.with_terminal() do
	if length(initialParams[1]) == 3
		# No Hill coefficient to report.
		println("Dataset\t\t\t\tKd\t\t\t\tSmin\t\t\tSmax")
		for (dataset, fit, stderr) in zip(datasetNames, fits, paramsStdErrors)
			println(
				split(dataset, "/")[end],
				"\t\t",
				round(fit.param[3], digits = 1),
				" ± ",
				round(stderr[3], digits = 1),
				"\t\t",
				round(fit.param[1], digits = 1),
				" ± ",
				round(stderr[1], digits = 1),
				"\t\t",
				round(fit.param[2], digits = 1),
				" ± ",
				round(stderr[2], digits = 1)
			)
		end
	elseif length(initialParams[1]) == 4
		# There is a Hill coefficient to report.
		println("Dataset\t\t\t\tKd\t\t\t\tSmin\t\t\tSmax\t\t\th")
		for (dataset, fit, stderr) in zip(datasetNames, fits, paramsStdErrors)
			println(
				split(dataset, "/")[end],
				"\t\t",
				round(fit.param[3], digits = 1),
				" ± ",
				round(stderr[3], digits = 1),
				"\t\t",
				round(fit.param[1], digits = 1),
				" ± ",
				round(stderr[1], digits = 1),
				"\t\t",
				round(fit.param[2], digits = 1),
				" ± ",
				round(stderr[2], digits = 1),
				"\t\t",
				round(fit.param[4], digits = 1),
				" ± ",
				round(stderr[4], digits = 1),
			)
		end
	else
		# Other number of values in the parameters array make no sense.
		println("Error.")
	end
end

# ╔═╡ Cell order:
# ╟─412eedac-658a-11eb-2326-93e5bf3d1a2c
# ╟─270cc0cc-660f-11eb-241e-b75746a39cc7
# ╠═a2c02fcf-9402-4938-bc3d-819b45a66afa
# ╠═14baf100-660f-11eb-1380-ebf4a860eed8
# ╟─2ce72e97-0133-4f15-bf1d-7fd04ccf3102
# ╟─77316a92-425a-4902-9828-52a7c4a74f27
# ╟─214acce6-6ae8-11eb-3abf-492e50140317
# ╟─d904fd76-6af1-11eb-2352-837e03072137
# ╟─0e8af3be-7ae7-4ec2-8d7a-670878cd52ee
# ╠═d5b0e2a1-865c-489c-9d0d-c4ae043828fb
# ╟─5c5e0392-658a-11eb-35be-3d940d4504cb
# ╟─3dd72c58-6b2c-11eb-210f-0b547bf38ebe
# ╟─3da83f72-6a11-11eb-1a74-49b66eb39c96
# ╟─d15cba72-6aeb-11eb-2c80-65702b48e859
# ╟─01b59d8a-6637-11eb-0da0-8d3e314e23af
# ╟─264bf9ec-6af5-11eb-1ffd-79fb3466f596
# ╟─4f4000b4-6b2c-11eb-015f-d76a0adda0a0
# ╟─c50cf18c-6b11-11eb-07d3-0b8e332ec5bc
# ╟─5a36fc3f-ce74-42c8-8284-19321e0d687f
# ╟─a951b5dc-6af7-11eb-2401-5d11a14e3067
# ╟─5392d99b-70f9-48cd-90b4-58cba5fc9681
# ╟─7625d41a-dab1-4b10-947c-3667c03f85aa
# ╟─be17b97e-663a-11eb-2158-a381c19ece3f
# ╟─090347fc-6b8e-11eb-0e17-9d9d45749c0b
# ╟─124c4f94-6b99-11eb-2921-d7c2cd00b893
# ╟─2109f516-6b99-11eb-05a0-99b9ecfd0f9d
# ╟─7e7a9dc4-6ae8-11eb-128d-83544f01b78b
# ╟─512e3028-6ae9-11eb-31b4-1bc9fc66b322
# ╠═393b2f5e-6556-11eb-2119-cf7309ee7392
# ╠═abc03f64-6a11-11eb-0319-ed7cea455cb5
# ╟─1884912a-6aeb-11eb-2b4a-d14d4a321dc5
# ╟─c94ef72d-bd12-4434-935e-01e94d5a4588
# ╟─54272681-dfeb-4034-b127-7e68c19fd576
# ╠═731492c6-95c7-449f-8c19-53e22ab438b8
# ╟─93ef0431-643b-4f6f-8c09-beabff59e0c6
# ╠═1c10231d-3bea-4468-a2d8-886c05c6474c
# ╠═730b693d-cca8-46de-8382-c151b5f63352
# ╠═24152a84-523b-4027-9b5a-7e3524b9c659
# ╟─97f73ee4-3db7-43ba-93eb-17025b485f4f
# ╠═8f2d959e-5e61-48a9-bbf0-538bc1d478a8
# ╟─c866d213-4680-4b08-8ecc-1faefc8661a4
# ╠═8e9cd30f-1722-4f2f-a26b-2f558805d4a1
# ╠═13417555-be95-4662-ad11-eca4dece81b5
# ╟─9c7e922e-c82f-41a4-9513-4462a0559c3f
# ╟─4f4b580d-507c-4ad0-b1d5-5967c8ed829e
# ╠═5eb607c7-172b-4a6c-a815-acbc195108f0
# ╟─992c86a2-6b13-11eb-1e00-95bdff2736d0
# ╠═1fe4d112-6a11-11eb-37a6-bf95fbe032b1
# ╟─5cc812af-4f8e-444b-9652-fb063cc6be06
# ╠═7c2af794-a0ed-4827-bc8e-1f15ad205eca
# ╟─a08d1b29-fa5e-4312-977c-39b050de4516
# ╠═0181a0b1-ec6c-4325-8b3a-851a3fe33846
# ╟─247c2416-6c67-11eb-01df-8dac01cdbf8f
# ╠═36ebe112-6c66-11eb-11f0-7fdb946865e4
# ╟─8e105fae-6aec-11eb-1471-83aebb776241
# ╟─97c5020c-6aec-11eb-024b-513b1e603d98
# ╠═caf4a4a2-6aec-11eb-2765-49d67afa47dd
# ╠═7c03fcbe-6a11-11eb-1b7b-cbad863156a6
# ╟─fc194672-6aed-11eb-0a06-2d967ec094b1
# ╠═14db987c-6aee-11eb-06cf-a11987b98f1e
# ╠═7f83b838-6a11-11eb-3652-bdff24f3473e
# ╟─9020fa5d-7408-4161-a52f-df37b3c2e6f5
# ╠═3184d209-1cc9-40ed-a9ec-9f094c5c94b5
# ╠═b38cd229-64e2-4ca4-a78b-8881ec166b09
# ╟─663a4cae-658a-11eb-382f-cf256c08c9d1
# ╟─594e7534-6aeb-11eb-1254-3b92b71877ed
# ╟─32dce844-6aee-11eb-3cf2-3ba420d311d3
# ╠═605f06ae-6a11-11eb-0b2f-eb81f6526829
# ╟─a1f56b0a-6aeb-11eb-0a44-556fad58f368
# ╠═5be2e5d2-6a11-11eb-1421-492f5af16f9c
# ╟─58617378-6aee-11eb-23e8-c13d89b4c57f
# ╠═6668c49a-6a11-11eb-2abf-5feecaee8972
# ╟─88d941e6-658a-11eb-08a2-0f021e5ae3a4
# ╠═e9fc3d44-6559-11eb-2da7-314e8fc76ee9
# ╟─9d1f24cc-6a0f-11eb-3b16-35f89aff5d4a
# ╠═78e664d0-6618-11eb-135b-5574bb05ddef
# ╟─b0b17206-6a0f-11eb-2f5e-5fc8fa06cd36
# ╠═5694f1da-6636-11eb-0fed-9fee5c48b114
# ╟─0f960a8a-6a0f-11eb-04e2-b543192f6354
# ╟─2eb890a4-658c-11eb-1fc4-af645d74109d
# ╠═ca2d2f12-6a1a-11eb-13ca-1f93df2b8e4a
# ╟─9be41e32-6af0-11eb-0904-d1cf3c288cab
# ╠═67b538f6-6a1b-11eb-3004-2d89c2f941e8
# ╟─213e8ffa-6a0f-11eb-357e-638146193c5d
# ╟─babcb896-6af0-11eb-194a-15922bc2df83
# ╠═47426056-6af2-11eb-17f8-6d27d35003ca
# ╟─799680d0-6af1-11eb-321d-b7758a40f931
# ╠═1f0384de-659b-11eb-043e-5b86fcdd36e6
# ╟─54501a10-6b9c-11eb-29de-77afc3772fb7
# ╠═5ed3ab64-6b9c-11eb-149e-43a1ef12ac7d
# ╟─8643b03c-6af1-11eb-0aa7-67acee28d2c0
# ╠═a74998b4-659c-11eb-354d-09ff62710b87

max_stars_repo_name: Guillawme/julia-curve-fitting

max_issues_repo_name: Guillawme/julia-curve-fitting

max_forks_repo_name: Guillawme/julia-curve-fitting


text:
######################
# Include Statements #
######################
    include("quarray.jl")
    include("arraymath.jl")
    include("constructors.jl")
    include("ladderops.jl")
    include("spinops.jl")

ext: jl


text:
using NeuralNetworks

net = Serial(Float64, (10, 2, 1), (6, 1))
add_layer!(net, View(Float64, (10, 2, 1), (20, 1)))
add_layer!(net, Linear(Float64, 1, 20, 6))
add_layer!(net, Sigmoid(Float64, 6, 1))
validate(net)

initWeights!(net)

input = reshape(linspace(1.0, 10.0, 20), (10, 2, 1))
output = forward!(net, input)

println(output)

gradOutput = reshape(linspace(1.0, 10.0, 6), (6, 1))
gradInput = backward!(net, input, gradOutput)

println(gradInput)


text:
immutable State
    intercept::Float64
    slope::Float64
end

immutable Datum
    x::Float64
    y::Float64
end

immutable LinregMHSampler <: Sampler
    num_steps::Int
    prior_std::Float64
    noise_std::Float64
    data::Array{Datum,1}
    function LinregMHSampler(x::Array{Float64,1}, y::Array{Float64,1},
                            prior_std::Float64, noise_std::Float64, num_steps::Int)
        @assert length(x) == length(y)
        data = Array{Datum,1}()
        for (xi, yi) in zip(x, y)
            push!(data, Datum(xi, yi))
        end
        new(num_steps, prior_std, noise_std, data)
    end
end

function prior_sample(sampler::LinregMHSampler)
    intercept = rand(Normal(0., sampler.prior_std))
    slope = rand(Normal(0., sampler.prior_std))
    State(intercept, slope)
end

function log_prior(sampler::LinregMHSampler, state::State)
    lp = logpdf(Normal(0., sampler.prior_std), state.intercept)
    lp += logpdf(Normal(0., sampler.prior_std), state.slope)
    lp
end

function log_likelihood(sampler::LinregMHSampler, state::State)
    ll = 0.0
    for datum in sampler.data
        y_expected = state.intercept + state.slope * datum.x
        ll += logpdf(Normal(y_expected, sampler.noise_std), datum.y)
    end
    ll
end

function log_unnormalized_posterior(sampler::LinregMHSampler, state::State)
    log_prior(sampler, state) + log_likelihood(sampler, state)
end

function imh_step(sampler::LinregMHSampler, state::State)
    # perform one step of independent MH
    # propose from prior
    state_proposed = prior_sample(sampler)
    if log(rand()) <= log_likelihood(sampler, state_proposed) - log_likelihood(sampler, state)
        # accept
        new_state = state_proposed
    else
        # reject
        new_state = state
    end
    new_state
end

function run_markov_chain(sampler::LinregMHSampler, state::State, steps::Int)
    for step in 1:steps
        state = imh_step(sampler, state)
    end
    state
end

function run_inference(sampler::LinregMHSampler)
    init_state = prior_sample(sampler)
    state = run_markov_chain(sampler, init_state, sampler.num_steps)
    # returns the initial and final state of the Markov chain
    # the initial state is all we need from the trace
    return (init_state, state)
end

function log_weight(sampler::LinregMHSampler, init_state::State, output::State)
    lw = 0.
    lw += log_prior(sampler, init_state)
    lw += log_unnormalized_posterior(sampler, output)
    lw -= log_unnormalized_posterior(sampler, init_state)
    lw
end

function simulate(sampler::LinregMHSampler)
    (init_state, state) = run_inference(sampler)
    ([state.intercept, state.slope], log_weight(sampler, init_state, state))
end

function regenerate(sampler::LinregMHSampler, intercept_and_slope::Array{Float64,1})
    @assert length(intercept_and_slope) == 2
    intercept, slope = intercept_and_slope
    state = State(intercept, slope)
    init_state = run_markov_chain(sampler, state, sampler.num_steps)
    log_weight(sampler, init_state, state)
end

immutable PriorSampler <: Sampler
    x::Array{Float64,1}
    y::Array{Float64,1}
    prior_std::Float64
    noise_std::Float64
end

function simulate(sampler::PriorSampler)
    dist = Normal(0., sampler.prior_std)
    intercept = rand(dist)
    slope = rand(dist)
    ([intercept, slope], logpdf(dist, intercept) + logpdf(dist, slope))
end

function regenerate(sampler::PriorSampler, intercept_and_slope::Array{Float64,1})
    (intercept, slope) = intercept_and_slope
    dist = Normal(0., sampler.prior_std)
    logpdf(dist, intercept) + logpdf(dist, slope)
end


text:
include("euler/euler.jl")
using .Primes: prime_numbers
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(n::Int)::Int
    result = 1
    for prime ∈ prime_numbers(trunc(Int, √n))
        result *= prime
        if result * prime > n
            return result
        end
    end
end

compute(10)

compute(1_000_000)

@benchmark compute(1_000_000)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
"""
    ObservableChain

A collection of simplices that acts like an observable. When you want to plot a chain that
changes over time, wrap it in this type. Updating it with `setindex` will also update the
plot.

# Example

julia> using Ripserer

julia> simplices = [Simplex{2}(1, 1), Simplex{2}(2, 2)];

julia> data = [(0, 0), (0, 1), (1, 0), (1, 1)];

julia> oc = ObservableChain(simplices, data)
ObservableChain(nv=4, ne=5, nt=2)

julia> oc[] = []; # or ()

julia> oc
ObservableChain(nv=0, ne=0, nt=0)

julia> oc[] = simplices;

julia> oc
ObservableChain(nv=4, ne=5, nt=2)

"""
struct ObservableChain{V,E,T,D}
    vertices::Observable{V}
    edges::Observable{E}
    triangles::Observable{T}
    data::D
end

function ObservableChain(chain, data)
    data = collect(data)
    return ObservableChain(
        Observable(convert_arguments(Scatter, chain, data)[1]),
        Observable(convert_arguments(LineSegments, chain, data)[1]),
        Observable(convert_arguments(Mesh, chain, data)[1]),
        convert_arguments(Scatter, data)[1],
    )
end
function ObservableChain(data)
    return ObservableChain(Simplex{1,Int,Int}[], data)
end

function Base.show(io::IO, ch::ObservableChain)
    nv = length(ch.vertices[])
    ne = length(ch.edges[]) ÷ 2
    nt = length(ch.triangles[])
    print(io, "ObservableChain(nv=$nv, ne=$ne, nt=$nt)")
    return nothing
end

function Base.setindex!(oc::ObservableChain, chain)
    if isempty(chain)
        oc.vertices[] = eltype(oc.vertices)[]
        oc.edges[] = eltype(oc.edges)[]
        oc.triangles[] = GeometryBasics.Mesh(oc.data, typeof(GLTriangleFace(1, 2, 3))[])
    else
        oc.vertices[] = convert_arguments(Scatter, chain, oc.data)[1]
        oc.edges[] = convert_arguments(LineSegments, chain, oc.data)[1]
        oc.triangles[] = convert_arguments(Mesh, chain, oc.data)[1]
    end
end

function AbstractPlotting.default_theme(scene::SceneLike, ::Type{<:Plot(ObservableChain)})
    return Theme(; palette=DEFAULT_PALETTE, color=1)
end

function AbstractPlotting.plot!(p::Plot(ObservableChain))
    oc = to_value(p[1])
    color = get_color(p, p[:color])
    dim_is_1 = @lift length($(oc.triangles)) == 0
    edgecolor = @lift $dim_is_1 ? $color : RGB(0.0, 0.0, 0.0)
    linewidth = @lift $dim_is_1 ? 5 : 1

    mesh!(p, oc.triangles; color, shading=false)
    linesegments!(p, oc.edges; color=edgecolor, linewidth, shading=false)
    return scatter!(p, oc.vertices; color)
end


text:
using Laplacians

fileName = "docs/API/wholeAPI.md"

function extraInfo(x)
    fh = IOBuffer()
    if isa(x,Function)
        mt = methods(x)
        println(fh,"\n```julia")
        loc = " "
        firstit = true
        for meth in mt
            str = string(meth)
            ind = rsearchindex(str," at ")
            println(fh, str[1:(ind-1)])
            if firstit
                ind2 = rsearchindex(str,"/")+1
                loc = str[ind2:end]
                firstit = false
            end
        end
        println(fh,"```\n")

        println(fh, loc)

    elseif isa(x,DataType)
        str = Docs.typesummary(x)
        writemime(fh, "text/plain", str)
    end

    str = takebuf_string(fh)
end

noDocString = "No documentation found."

nms = names(Laplacians)

n = length(nms)

docstrs = Array(AbstractString,n)
linenums = zeros(Int64,n)
fileins = Array(AbstractString,n)
extras = Array(AbstractString,n)
strings = Array(AbstractString,n)

for i in 1:length(nms)
    sym = nms[i]
    x = eval(sym)

    strings[i] = string(sym)

    try

        if isa(x,Function)
            mt = methods(x)
            file, line = functionloc(mt.defs)
            fileins[i] = rsplit(file,"/")[end]
            linenums[i] = line
        end

        docmd = @doc string(x)

        docstr = stringmime("text/plain", docmd )
        if (length(docstr) < 23) || (docstr[1:23] != noDocString)
            docstrs[i] = docstr;
        end

        extras[i] = extraInfo(x)

    catch
        println("didn't generate docs for ", string(sym))
    end

end

ind = zeros(Bool,length(fileins))
for i in 1:length(fileins)
    ind[i] = isdefined(fileins,i)
end

u = unique(fileins[ind]);

fileins[!ind] = ""

for j in 1:length(u)
    infile = find(fileins .== u[j])
    infile = infile[sortperm(linenums[infile])]

    fn = split(u[j],'.')[1]

    fileName = "docs/API/" * fn * "API.md"
    fh = open(fileName, "w")

    println(fh, "# ", fn)

    for i in infile
        if isdefined(docstrs,i)
            println(fh, "### ", strings[i])
            println(fh, docstrs[i])
            println(fh, extras[i])
            println(fh,"\n")
        end
    end

    close(fh)

end


text:
using Distributions
using Gadfly
using RDatasets

function spiral(radius, theta, sequence)
    rs = map(radius, sequence)
    ts = map(theta, sequence)
    xs = map((r, t) -> r * sin(t), rs, ts)
    ys = map((r, t) -> r * cos(t), rs, ts)

    return [xs ys]
end

typealias Weight Float64
typealias Value Float64
typealias Point Vector{Value}
typealias Class Matrix{Value}
typealias DataSet Vector{Class}

sigmoid = n -> 1 / (1 + exp(-n))

unvec = (x) -> reshape(x, 1, length(x))

function computeOutput(weights, input)
    output = Matrix{Value}[input]
    for w in weights
        push!(output, map(sigmoid, last(output) * w))
    end
    return output
end

function classify(weights, input)
    output = input
    for w in weights
        output = map(sigmoid, output * w)
    end
    return indmax(output)
end

function mse(target, output)
    diff = target - output
    return sum(diff * diff) / 2
end

function learn(weights, input, target, learningRate = 0.8)
    n = size(weights, 1)
    outputs = computeOutput(weights, input)
    deltas = fill(Matrix{Value}(1, 1), n)

    # tabulate the output layer first
    o = last(outputs)
    deltas[n] = (target - o) .* o .* (1 - o)

    # tabulate from output to input
    for i in (n-1):-1:1
        o = outputs[i+1]
        w = weights[i+1]
        d = deltas[i+1]
        sums = (w * d')'
        deltas[i] = sums .* o .* (1 - o)
    end

    for i in 1:n
        delta = [d * o for o = outputs[i], d = deltas[i]]
        if (all(abs(delta) .< 0.0001))
            weights[i] = weights[i] + delta * learningRate * 1000
        else
            weights[i] = weights[i] + delta * learningRate
        end
    end

    return Dict(
        "weights" => weights,
        "error" => sum(map(mse, target, last(outputs)))
    )
end

function train(weights, dataset, iteration = 1000)
    n = size(dataset, 1)
    totalSize = sum([size(class, 1) for class = dataset])
    errors = []
    for i in 1:iteration
        accumulatedError = 0
        for (j, class) in enumerate(dataset)    # for each class
            target = zeros((1, n))
            target[j] = 1
            m = size(class, 1)
            for k in 1:m
                result = learn(weights, [class[k,:] 1], target)    # augmented
                weights = result["weights"]
                accumulatedError = accumulatedError + result["error"]
            end
        end
        push!(errors, accumulatedError)
        if (accumulatedError < 0.005 * totalSize)
            break
        end
        if (i % 100 == 0)
            @show (i, last(errors))
        end
    end
    return Dict(
        "weights" => weights,
        "error" => errors
    )
end

function generateRandomWeights(layers)
    n = size(layers, 1)
    w = Array{Weight}[]
    for i in 1:(n-1)
        s = layers[i] * layers[i+1]
        push!(w, reshape(rand(s), layers[i], layers[i+1]))
    end
    return w
end

function generateDataSet(dataset)
    result = DataFrame(x = [], y = [], class = [])
    for (i, class) in enumerate(dataset)
        result = [result; DataFrame(x = class[:,1], y = class[:,2], class = i)]
    end
    return result
end

function generateGrid(seq, weights)
    n = length(seq)
    xs = repeat(seq, inner=[n])
    ys = repeat(seq, outer=[n])
    rs = map((x, y) -> classify(weights, [x y 1]), xs, ys)
    return DataFrame(x = xs, y = ys, class = rs)
end

function drawPoint(data)
    return layer(data, x = "x", y = "y", color = "class", Geom.point,
        Theme(default_point_size = 3px, highlight_width = 1.5px))
end

function drawGrid(data)
    return layer(data, x = "x", y = "y", color = "class", Geom.rectbin)
end

function drawError(data)
    x = collect(1:length(data))
    return layer(x = x, y = data, Geom.line)
end

################################################################################
#   2
################################################################################

# PDF
spiralPDFS = Cairo.CairoPDFSurface("spiral.pdf", 600, 600)
spiralPDFC = Cairo.CairoContext(spiralPDFS)
spiralPDF = Compose.CAIROSURFACE(spiralPDFS)

# data points
spiral0 = spiral(i -> 6.5 * (104 - i) / 104,
    i -> pi * i / 16,
    40:96)

spiral1 = spiral(i -> 6.5 * (104 - i) / -104,
    i -> pi * i / 16,
    40:96)
spiralDataPoints = Class[spiral0, spiral1]

# weights
spiralWeights = generateRandomWeights([3, 80, 20, 10, 2])
spiralResult = train(spiralWeights, spiralDataPoints, 10000)
spiralWeights = spiralResult["weights"]
spiralErrors = spiralResult["error"]

# plot spirals
spiralDataSet = generateDataSet(spiralDataPoints)
spiralDataSetLayer = drawPoint(spiralDataSet)

# plot the grid
spiralGrid = generateGrid(collect(-8:0.1:8), spiralWeights)
spiralGridLayer = drawGrid(spiralGrid)

spiralPlot = plot(spiralDataSetLayer, spiralGridLayer,
    Coord.cartesian(fixed = true),
    Guide.title("The Two Spirals Problem"),
    Scale.color_discrete_manual(colorant"coral", colorant"cornflowerblue"))

spiralErrorsLayer = drawError(spiralErrors)
spiralErrorsPlot = plot(spiralErrorsLayer,
    Coord.Cartesian(ymin = 0),
    Guide.ylabel("Error"),
    Guide.xlabel("Iteration"))

draw(spiralPDF, spiralPlot)
Cairo.show_page(spiralPDFC)
draw(spiralPDF, spiralErrorsPlot)
Cairo.finish(spiralPDFS)

################################################################################
#   2
################################################################################

moonsPDFS = Cairo.CairoPDFSurface("moons.pdf", 600, 600)
moonsPDFC = Cairo.CairoContext(moonsPDFS)
moonsPDF = Compose.CAIROSURFACE(moonsPDFS)

# the double moon problem
n = 200
thetas = linspace(-180, 180, n)*pi/360;
r = 8
x0 = -5 + r * sin(thetas)' + randn(1, n)
y0 =      r * cos(thetas)' + randn(1, n)
x1 =  5 + r * sin(thetas)' + randn(1, n)
y1 =     -r * cos(thetas)' + randn(1, n)
moons = Class[[vec(x0) vec(y0)], [vec(x1) vec(y1)]]

#
moonsWeights = generateRandomWeights([3, 8, 8, 2])
moonsResult = train(moonsWeights, moons, 10000)
moonsWeights = moonsResult["weights"]
moonsErrors = moonsResult["error"]

# plot moons
moonsDS = generateDataSet(moons)
moonsDSLayer = drawPoint(moonsDS)

# plot the grid
moonsGrid = generateGrid(collect(-18:0.5:18), moonsWeights)
moonsGridLayer = drawGrid(moonsGrid)

moonsPlot = plot(moonsDSLayer, moonsGridLayer,
    Coord.cartesian(fixed = true),
    Guide.title("The Double Moons Problem"),
    Scale.color_discrete_manual(colorant"coral", colorant"cornflowerblue"))

moonsErrorsLayer = drawError(moonsResult["error"])
moonsErrorsPlot = plot(moonsErrorsLayer,
    Guide.ylabel("Error"),
    Guide.xlabel("Iteration"))

draw(moonsPDF, moonsPlot)
Cairo.show_page(moonsPDFC)
draw(moonsPDF, moonsErrorsPlot)
Cairo.finish(moonsPDFS)

################################################################################
#   3
################################################################################

# PDF
mvrPDFS = Cairo.CairoPDFSurface("mvrs.pdf", 600, 600)
mvrPDFC = Cairo.CairoContext(mvrPDFS)
mvrPDF = Compose.CAIROSURFACE(mvrPDFS)

mvr0 = rand(MvNormal([0.0, 0.0], [1.0 0.0; 0.0 1.0]), 100)'
mvr1 = rand(MvNormal([10.0, 0.0], [1.0 0.0; 0.0 4.0]), 100)'
mvr2 = rand(MvNormal([5.0, 10.0], [4.0 0.0; 0.0 1.0]), 100)'
mvr3 = rand(MvNormal([5.0, 5.0], [1.0 0.0; 0.0 1.0]), 100)'
mvrDataPoints = Class[mvr0, mvr1, mvr2, mvr3]

# weights
mvrWeights = generateRandomWeights([3, 10, 10, 4])
mvrResult = train(mvrWeights, mvrDataPoints, 10000)
mvrWeights = mvrResult["weights"]
mvrErrors = mvrResult["error"]

# plot distributions
mvrDataSet = generateDataSet(mvrDataPoints)
mvrDataSetLayer = drawPoint(mvrDataSet)

# plot the grid
mvrGrid = generateGrid(collect(-8:0.5:18), mvrWeights)
mvrGridLayer = drawGrid(mvrGrid)

mvrPlot = plot(mvrDataSetLayer, mvrGridLayer,
    Coord.cartesian(fixed = true),
    Guide.title("4 Multivariate Distributions"))
    # Scale.color_discrete_manual(colorant"coral", colorant"cornflowerblue"))

mvrErrorsLayer = drawError(mvrErrors)
mvrErrorsPlot = plot(mvrErrorsLayer,
    Coord.Cartesian(ymin = 0),
    Guide.ylabel("Error"),
    Guide.xlabel("Iteration"))

draw(mvrPDF, mvrPlot)
Cairo.show_page(mvrPDFC)
draw(mvrPDF, mvrErrorsPlot)
Cairo.finish(mvrPDFS)

max_stars_repo_name: banacorn/neuron

max_issues_repo_name: banacorn/neuron

max_forks_repo_name: banacorn/neuron


text:
module DiffEqML

using DiffEqBase, Flux, DiffResults, DiffEqSensitivity, ForwardDiff

include("Flux/layers.jl")

export diffeq_fd, diffeq_rd, diffeq_adjoint
end


text:
using Test, GeodesicTracer

struct TestMetricParams{T} <: GeodesicTracer.AbstractMetricParams{T}
    test::T
end

GeodesicTracer.constrain(m::TestMetricParams{T}, u, v; μ = μ) where {T} = sum(v)

SVector = GeodesicTracer.SVector

@testset "Constraints" begin
    tmp = TestMetricParams(1.0)

    @testset "Vectors" begin
        u = Float64[1.0, 1.0, 1.0, 1.0]
        v = Float64[1.0, 2.0, -1.0, -2.0]
        expected = [0.0, 2.0, -1.0, -2.0]

        res_v = GeodesicTracer.constrain_all(tmp, u, v, 0.0)
        @test res_v isa Vector{Float64}
        @test all(res_v .== expected)

        u_static = SVector(u...)
        v_static = SVector(1.0, 2.0, -1.0, -2.0)
        res_v_static = GeodesicTracer.constrain_all(tmp, u_static, v_static, 0.0)
        @test all(res_v_static .== expected)
    end

    @testset "Vector of Vectors" begin
        us = Vector{Float64}[
            [1.0, 1.0, 1.0, 1.0],
            [1.0, 1.0, 1.0, 1.0],
            [1.0, 1.0, 1.0, 1.0]
        ]
        vs = Vector{Float64}[
            [1.0, 2.0, -1.0, -2.0],
            [2.0, 2.0, -1.0, -2.0],
            [3.0, 2.0, -1.0, -2.0]
        ]

        res_vs = GeodesicTracer.constrain_all(tmp, us, vs, 0.0)
        for (i, v) in enumerate(res_vs)
            @test v[1] .≈ i - 1.0
        end
    end

end


text:
module DocumentationGenerator
using Pkg
using JSON, Downloads

include("utils/misc.jl")
include("utils/pkgtools.jl")
include("utils/runners.jl")
include("builders.jl")
include("metadata.jl")
include("registry.jl")

## when job doesn't log anything to stdout
const RUNNER_TIMEOUT = 60*60
## maxtime a job is allowed to run
const RUNNER_MAX_TIMEOUT = 3*60*60
## graceful termination timeout
const RUNNER_KILL_TIMEOUT = 60

function try_install_package(packagespec, envdir)
    @assert ispath(envdir)
    success = false
    try
        Pkg.activate(envdir)
        Pkg.add(packagespec)
        success = true
    catch err
        @error("Could not `Pkg.add($(packagespec))`.", exception=err)
    end
    return success
end

function try_use_package(packagespec)
    pkg_sym = Symbol(packagespec.name)

    pkg_module = try
        @eval(Main, (using $pkg_sym; $pkg_sym))
    catch err
        @error("`using $(pkg_sym) did not succeed.`", exception=err)
        nothing
    end

    return pkg_module
end

function build_package_docs(packagespec::Pkg.Types.PackageSpec, buildpath, registry; src_prefix="", href_prefix="")
    type, uri = doctype(packagespec, registry)

    @info("$(packagespec.name) specifies docs of type $(type).")
    out = try
        if type == "hosted"
            build_hosted_docs(packagespec, buildpath, uri)
        elseif type == "git-repo"
            build_git_docs(packagespec, buildpath, uri; src_prefix=src_prefix, href_prefix=href_prefix)
        elseif type == "vendored"
            build_local_docs(packagespec, buildpath, uri; src_prefix=src_prefix, href_prefix=href_prefix)
        else
            @error("Invalid doctype specified: $(type).")
            Dict(
                "success" => false
            )
        end
    catch err
        @error("Error while generating docs.", exception=(err, catch_backtrace()))
        Dict(
            "success" => false
        )
    end

    return out
end

function build_documentation(
        packages;
        processes::Int = 8,
        sleeptime = 0.5,
        juliacmd = first(Base.julia_cmd()),
        basepath = joinpath(@__DIR__, ".."),
        envpath = normpath(joinpath(@__DIR__, "..")),
        filter_versions = last,
        sync_registry = true,
        deployment_url = "pkg.julialang.org/docs",
        update_only = false,
        registry = joinpath(homedir(), ".julia/registries/General"),
        timeout = RUNNER_TIMEOUT,
        max_timeout = RUNNER_MAX_TIMEOUT,
        kill_timeout = RUNNER_KILL_TIMEOUT
    )

    has_xvfb = try
        success(`xvfb-run --help`)
    catch err
        @warn("No `xvfb` installed. Running without it.")
        false
    end

    regpath = get_registry(basepath, sync = sync_registry)
    process_queue = []

    # make sure registry is updated *before* we start multiple processes that might try that at the same time
    Pkg.Registry.update()

    envmod = []

    local x_server_proc
    if has_xvfb
        display_server = string(':', find_free_x_servernum())

        @info("Running Xvfb on display $(display_server).")

        x_server_proc = run(`Xvfb $(display_server)`, wait=false)
        envmod = ["DISPLAY" => display_server]
    end

    withenv(envmod...) do
        for package in packages
                # make sure we're not queueing new processes over the limit
                while length(process_queue) >= processes
                    filter!(process_running, process_queue)
                    sleep(sleeptime)
                end

                # separate process for each version of a package
                for version in vcat(filter_versions(package.versions))
                    proc = start_builder(package, version;
                                           basepath = basepath,
                                           juliacmd = juliacmd,
                                           registry_path = regpath,
                                           deployment_url = deployment_url,
                                           update_only = update_only,
                                           timeout = timeout,
                                           max_timeout = max_timeout,
                                           kill_timeout = kill_timeout)
                    push!(process_queue, proc)
                end
        end

        # wait for all queued processes to finish
        for proc in process_queue
            wait(proc)
        end
    end

    if has_xvfb
        kill(x_server_proc)
    end

    # record dependency relations specified in registry
    generate_dependency_list(packages, basepath = basepath, registry = registry, filter_versions = filter_versions)
end

function http_get_request(url)
    output = IOBuffer()
    response = Downloads.request(url, output = output, method = "GET")
    (response, String(take!(output)))
end

function get_pkg_eval_data()
    pkg_eval = Dict()

    resp, bodystring = try
        http_get_request("https://raw.githubusercontent.com/JuliaCI/NanosoldierReports/master/pkgeval/by_date/latest")
    catch ex
        @warn "Failed to fetch latest Nanosoldier report" ex
        return pkg_eval
    end

    if resp.status == 200
        latest_date = bodystring

        if occursin(r"\d{4}\-\d{2}/\d{2}", latest_date)
            last_db_url = "https://raw.githubusercontent.com/JuliaCI/NanosoldierReports/master/pkgeval/by_date/$(latest_date)/db.json"
            resp, body_string = try
                 http_get_request(last_db_url)
            catch ex
                @warn "Failed to fetch Nanosoldier report" ex
                nothing
            end
            if resp != nothing && resp.status == 200
                pkg_eval = JSON.parse(body_string)
            end
        end
    end

    return pkg_eval
end

function generate_dependency_list(packages;
        basepath = joinpath(@__DIR__, ".."),
        registry = joinpath(homedir(), ".julia/registries/General"),
        filter_versions = last
    )
    @info "Generating deps info"
    pkg_eval_data = get_pkg_eval_data()
    deps = dependencies_per_package(registry)
    rdeps = reverse_dependencies_per_package(deps)
    for package in packages
        for version in vcat(filter_versions(package.versions))
            try
                builddir = joinpath(basepath, "build", get_docs_dir(package.name, package.uuid), string(version))
                metatoml = joinpath(builddir, "meta.toml")
                isfile(metatoml) || continue

                meta = Pkg.TOML.parsefile(metatoml)
                 if haskey(pkg_eval_data, "tests") && haskey(pkg_eval_data["tests"], package.uuid)
                    meta["pkgeval"] = pkg_eval_data["tests"][package.uuid]
                end
                meta["deps"] = collect(alldeps(package.uuid, string(version), deps))
                meta["reversedeps"] = collect(allreversedeps(package.uuid, string(version), rdeps))
                open(metatoml, "w") do io
                    if VERSION < v"1.6"
                        Pkg.TOML.print(io, meta)
                    else
                        Pkg.TOML.print(x -> string(x), io, meta)
                    end
                end
                open(joinpath(builddir, "pkg.json"), "w") do f
                    readme = joinpath(builddir, "_readme", "readme.html")
                    isfile(readme) && (meta["readme"] = read(readme, String))
                    print(f, JSON.json(meta))
                end
            catch err
                @error(exception=(err, catch_backtrace()))
            end
        end
    end
end

function start_builder(package, version;
        basepath = error("`basepath` is a required argument."),
        juliacmd = error("`juliacmd` is a required argument."),
        registry_path = error("`registry_path` is a required argument."),
        deployment_url = error("`deployment_url` is a required argument."),
        update_only = error("`update_only` is a required argument."),
        src_prefix = nothing,
        href_prefix = nothing,
        timeout = RUNNER_TIMEOUT,
        max_timeout = RUNNER_MAX_TIMEOUT,
        kill_timeout = RUNNER_KILL_TIMEOUT
    )

    workerfile = joinpath(@__DIR__, "workerfile.jl")
    buildpath = joinpath(basepath, "build")
    logpath = joinpath(basepath, "logs")

    isdir(buildpath) || mkpath(buildpath)
    isdir(logpath) || mkpath(logpath)

    name = package.name
    uuid = package.uuid
    url = package.url
    src_prefix  = haskey(package, :src_prefix) ? package.src_prefix : string("/docs/", get_docs_dir(name, uuid), '/', string(version), "/_packagesource/")
    href_prefix = haskey(package, :href_prefix) ? package.href_prefix : string("/ui/Code/docs/", get_docs_dir(name, uuid), '/', string(version), "/_packagesource/")

    builddir = joinpath(buildpath, get_docs_dir(name, uuid), string(version))
    isdir(builddir) || mkpath(builddir)

    logfile = joinpath(builddir, "..", "$(version).log")

    thisproject = Base.active_project()

    cmd = ```
        $(juliacmd)
            --project="$(thisproject)"
            --color=no
            --compiled-modules=no
            -O0
            $workerfile
            $uuid
            $name
            $url
            $version
            $builddir
            $registry_path
            $deployment_url
            $src_prefix
            $href_prefix
            $(update_only ? "update" : "build")
    ```
    process, task = run_with_timeout(cmd,
                                     log=logfile,
                                     name = string("docs build for ", name, "@", version, " (", uuid, ")"),
                                     timeout = timeout,
                                     max_timeout = max_timeout,
                                     kill_timeout = kill_timeout)
    return process
end
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/optimize/optimize_convex.jl
# Contains the solve_local_nlp! routine which computes the optimal value
# of a convex function. This is used to compute the upper bound in the
# branch and bound routine. A number of utility function required for
# solve_local_nlp! are also included.
#############################################################################

"""

Shifts the resulting local nlp objective value `f*` by `(1.0 + relative_tolerance/100.0)*f* + absolute_tolerance/100.0`.
This assumes that the local solvers relative tolerance and absolute tolerance is significantly lower than the global
tolerance (local problem is minimum).
"""
function stored_adjusted_upper_bound!(d::Optimizer, v::Float64)
    adj_atol = d._parameters.absolute_tolerance/100.0
    adj_rtol = d._parameters.relative_tolerance/100.0
    if v > 0.0
        d._upper_objective_value = v*(1.0 + adj_rtol) + adj_atol
    else
        d._upper_objective_value = v*(1.0 - adj_rtol) + adj_atol
    end

    return nothing
end

revert_adjusted_upper_bound!(t::ExtensionType, d::Optimizer) = nothing

function revert_adjusted_upper_bound!(t::DefaultExt, d::Optimizer)

    adj_atol = d._parameters.absolute_tolerance/100.0
    adj_rtol = d._parameters.relative_tolerance/100.0

    adj_objective_value = d._global_upper_bound
    adj_objective_value -= adj_atol
    if adj_objective_value > 0.0
        adj_objective_value /= (1.0 + adj_rtol)
    else
        adj_objective_value /= (1.0 - adj_rtol)
    end
    d._global_upper_bound = adj_objective_value

    return nothing
end

# translates quadratic cone
function add_soc_constraints_as_quad!(m::Optimizer, opt::T) where T

    for (func, set) in m._input_problem._conic_second_order
        # quadratic cone implies variable[1] >= 0.0, bounds contracted accordingly in initial_parse!
        quad_terms = SQT[SQT((), func.variables[i], func.variables[i]) for i = 1:length(func.variables)]
        sqf = SQF(SQT[], SAF[], 0.0)
        MOI.add_constraint(opt, sqf, LT_ZERO)
    end

    return nothing
end

"""

Constructs and solves the problem locally on on node `y` updated the upper
solution informaton in the optimizer.
"""
function solve_local_nlp!(m::Optimizer)

    upper_optimizer = m.upper_optimizer
    MOI.empty!(upper_optimizer)

    upper_variables = m._upper_variables
    for i = 1:m._working_problem._variable_count
        @inbounds upper_variables[i] = MOI.add_variable(upper_optimizer)
    end

    n = m._current_node
    sol_to_branch_map = m._sol_to_branch_map
    lower_variable_bounds = n.lower_variable_bounds
    upper_variable_bounds = n.upper_variable_bounds
    variable_info = m._input_problem._variable_info

    lvb = 0.0
    uvb = 0.0
    x0 = 0.0

    for i = 1:m._input_problem._variable_count
        vinfo = @inbounds variable_info[i]
        single_variable = MOI.SingleVariable(@inbounds upper_variables[i])

        if vinfo.branch_on === BRANCH
            if vinfo.is_integer
            else
                indx = @inbounds sol_to_branch_map[i]
                lvb  = @inbounds lower_variable_bounds[indx]
                uvb  = @inbounds upper_variable_bounds[indx]
                if vinfo.is_fixed
                    MOI.add_constraint(upper_optimizer, single_variable, ET(lvb))

                elseif vinfo.has_lower_bound
                    if vinfo.has_upper_bound
                        MOI.add_constraint(upper_optimizer, single_variable, LT(uvb))
                        MOI.add_constraint(upper_optimizer, single_variable, GT(lvb))

                    else
                        MOI.add_constraint(upper_optimizer, single_variable, GT(lvb))

                    end
                elseif vinfo.has_upper_bound
                    MOI.add_constraint(upper_optimizer, single_variable, LT(uvb))

                end
            end
            x0 = 0.5*(lvb + uvb)
            upper_variable_index = @inbounds upper_variables[i]
            MOI.set(upper_optimizer, MOI.VariablePrimalStart(), upper_variable_index, x0)

        else
            # not branch variable
            if vinfo.is_integer
            else
                lvb  = vinfo.lower_bound
                uvb  = vinfo.upper_bound
                if vinfo.is_fixed
                    MOI.add_constraint(upper_optimizer, single_variable, ET(lvb))

                elseif vinfo.has_lower_bound
                    if vinfo.has_upper_bound
                        MOI.add_constraint(upper_optimizer, single_variable, LT(uvb))
                        MOI.add_constraint(upper_optimizer, single_variable, GT(lvb))

                    else
                        MOI.add_constraint(upper_optimizer, single_variable, GT(lvb))

                    end
                elseif vinfo.has_upper_bound
                    MOI.add_constraint(upper_optimizer, single_variable, LT(uvb))
                end
                x0 = 0.5*(lvb + uvb)
                upper_variable_index = @inbounds upper_variables[i]
                MOI.set(upper_optimizer, MOI.VariablePrimalStart(), upper_variable_index, x0)
            end
        end
    end

    # Add linear and quadratic constraints to model
    add_linear_constraints!(m, upper_optimizer)

    for (func, set) in m._input_problem._quadratic_leq_constraints
        MOI.add_constraint(upper_optimizer, func, set)
    end
    for (func, set) in m._input_problem._quadratic_geq_constraints
        MOI.add_constraint(upper_optimizer, func, set)
    end
    for (func, set) in m._input_problem._quadratic_eq_constraints
        MOI.add_constraint(upper_optimizer, func, set)
    end

    if MOI.supports_constraint(upper_optimizer, VECOFVAR, SOC)
        add_soc_constraints!(m, upper_optimizer)
    else
        add_soc_constraints_as_quad!(m, upper_optimizer)
    end

    # Add nonlinear evaluation block
    MOI.set(upper_optimizer, MOI.NLPBlock(), m._working_problem._nlp_data)
    MOI.set(upper_optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    # set objective as NECESSARY
    add_sv_or_aff_obj!(m, upper_optimizer)
    if m._input_problem._objective_type === SCALAR_QUADRATIC
        MOI.set(upper_optimizer, MOI.ObjectiveFunction{SQF}(), m._input_problem._objective_sqf)
    end

    # Optimizes the object
    MOI.optimize!(upper_optimizer)

    # Process output info and save to CurrentUpperInfo object
    m._upper_termination_status = MOI.get(upper_optimizer, MOI.TerminationStatus())
    m._upper_result_status = MOI.get(upper_optimizer, MOI.PrimalStatus())

    if is_feasible_solution(m._upper_termination_status, m._upper_result_status)
        m._upper_feasibility = true
        value = MOI.get(upper_optimizer, MOI.ObjectiveValue())
        stored_adjusted_upper_bound!(m, value)
        m._best_upper_value = min(value, m._best_upper_value)
        m._upper_solution .= MOI.get(upper_optimizer, MOI.VariablePrimal(), upper_variables)

    else
        m._upper_feasibility = false
        m._upper_objective_value = Inf

    end

    return nothing
end

function optimize!(::Val{DIFF_CVX}, m::Optimizer)

    solve_local_nlp!(m)

    return nothing
end


text:
function fourfilt(x,delt,tmax,tmin)
#  FOURFILT  Fourier low, high, or bandpass filter.
#
#     [filtdat]=fourfilt(x,delt,tmax,tmin)
#
#     where:   x:      data series to be filtered
#              delt:   sampling interval
#              tmax:   maximum period filter cutoff
#              tmin:   minimum period filter cutoff
#
#     usage:  lowpassdata=fourfilt(data,0.5,2000,20)
#
#               gives lowpass filter with cutoff at 20.0 sec
#               tmax set > (length(x)*delt) for no cutoff at low freq end
#
#     usage:  highpassdata=fourfilt(x,0.5,20,0.9)
#
#               gives highpass filter with cutoff at 20.0 sec
#               tmin set < (2*delt) for no cutoff at high freq end
#
#     usage:  bandpassdata=fourfilt(x,0.5,20,10)
#
#               gives bandpass filter passing 10-20 sec. band
#
#     REQUIREMENTS: Julia "FFTW" package
#
# Reference:
# Walters, R. A.  and Heston, C., 1982. Removing the tidal-period variations from time-series
# data using low-pass digital filters. Journal of Physical Oeanography, 12 112-115 .
#
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Julia version that has been translated from MATLAB version
# by William Pringle on 08/21/2019
# MATLAB verion is available from Rich Signell's "RPSstuff" github
#
# Version 1.0  (12/4/96)  Jeff List (jlist@usgs.gov)
# Version 1.1  (1/8/97)  Rich Signell (rsignell@usgs.gov)
#     removed argument for number of points and add error trapping for matrices
# Version 1.1b (12/1/2005) Rich Signell (rsignell@usgs.gov)
#     added reference
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global f

if ndims(x) != 1
  display("fourfilt cannot handle matrices");
  return
end
if sum(isnan.(x)) > 0
  display("fourfilt will not work with nans in timeseries")
  return
end

npts = length(x);
nby2 = convert(Int,floor(npts/2));
tfund = npts*delt;
ffund = 1.0/tfund;

#  remove the mean from data:
datamean = mean(x);
x = x .- datamean;

# fourier transform data:
coeffs = fft(x);

#  filter coefficients:
f = ffund
for i = 2:nby2+1
  t = 1/f
  if t > tmax || t < tmin
     coeffs[i] = 0.0*coeffs[i];
  end
  f += ffund
end

#  calculate the remaining coefficients:
for i = 2:nby2
   coeffs[npts+2-i] = conj(coeffs[i])
end

#  backtransform data and take real part:
backcoeffs = ifft(coeffs);
filtdat    = real(backcoeffs);

# add back the mean:
return filtdat .+ datamean

end


text:
"""
    calc_boundary_state!(::NumericalFluxGradient, ::Impenetrable{FreeSlip}, ::ModelSetup)
apply free slip boundary condition for velocity
sets non-reflective ghost point
"""
function calc_boundary_state!(
    ::NumericalFluxGradient,
    ::Impenetrable{FreeSlip},
    ::Union{ModelSetup},
    ::Union{AbstractDiffusion, Nothing},
    state⁺,
    aux⁺,
    n⁻,
    state⁻,
    aux⁻,
    t,
    _...,
)
    state⁺.ρ = state⁻.ρ

    ρu⁻ = state⁻.ρu
    state⁺.ρu = ρu⁻ - n⁻ ⋅ ρu⁻ .* SVector(n⁻)

    return nothing
end

"""
    calc_boundary_state!(::NumericalFluxGradient, ::Impenetrable{NoSlip}, ::ModelSetup)
apply no slip boundary condition for velocity
set numerical flux to zero for U
"""
@inline function calc_boundary_state!(
    ::NumericalFluxGradient,
    ::Impenetrable{NoSlip},
    ::Union{ModelSetup},
    ::Union{AbstractDiffusion, Nothing},
    state⁺,
    aux⁺,
    n⁻,
    state⁻,
    aux⁻,
    t,
    _...,
)
    FT = eltype(state⁺)
    state⁺.ρu = @SVector zeros(FT, 3)

    return nothing
end

"""
    calc_boundary_state!(::NumericalFluxGradient, ::Insulating, ::HBModel)

apply insulating boundary condition for temperature
sets transmissive ghost point
"""
@inline function calc_boundary_state!(
    ::Union{NumericalFluxGradient},
    ::Union{Insulating, CoupledPrimaryBoundary},
    ::Union{ModelSetup},
    ::Union{AbstractDiffusion, Nothing},
    state⁺,
    aux⁺,
    n⁻,
    state⁻,
    aux⁻,
    t,
    _...,
)
    state⁺.ρθ = state⁻.ρθ

    return nothing
end

calc_boundary_state!(
    ::Union{NumericalFluxGradient},
    ::Union{Insulating, CoupledSecondaryBoundary},
    ::Union{SlabLandModelSetup},
    _...,
) = nothing

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
# StackedLSTM.jl
#
# Layers for stacked LSTM

mutable struct StackedLSTMCell{A,S}
	chain::A
	state::S
end

"""
	StackedLSTM(in, out, hiddensize, layers)

Stacked LSTM network. Feeds the data through a chain of LSTM layers, where the hidden state
of the previous layer gets fed to the next one. The first layer corresponds to
`LSTM(in, hiddensize)`, the hidden layers to `LSTM(hiddensize, hiddensize)`, and the final
layer to `LSTM(hiddensize, out)`. Takes the keyword argument `init` for the initialization
of the layers.
"""
function StackedLSTM(in::Integer, out::Integer, hiddensize::Integer, layers::Integer;
			init=Flux.glorot_uniform)
	if layers == 1
		chain = Chain(LSTM(in, out; init=init))
	elseif layers == 2
		chain = Chain(LSTM(in, hiddensize; init=init),
					  LSTM(hiddensize, out; init=init))
	else
		chain_vec = [LSTM(in, hiddensize; init=init)]
		for i = 1:layers - 2
			push!(chain_vec, LSTM(hiddensize, hiddensize; init=init))
		end
		chain = Chain(chain_vec..., LSTM(hiddensize, out; init=init))
	end
	return StackedLSTMCell(chain,  zeros(Float32, out))
end

function (m::StackedLSTMCell)(x)
	out = m.chain(x)
	m.state = out
	return out
end

Flux.@functor StackedLSTMCell
Flux.trainable(m::StackedLSTMCell) = (m.chain)

function Base.show(io::IO, l::StackedLSTMCell)
	if length(l.chain) == 1
		print(io, l.chain[1])
	else
		print(io, "StackedLSTM(", size(l.chain[1].cell.Wi, 2), ", ", size(l.chain[end].cell.Wi, 1) ÷ 4)
		print(io, ", ", size(l.chain[end].cell.Wi, 2))
		print(io, ", ", length(l.chain), ")")
	end
end

# Initialize forget gate bias with 1
function initialize_bias!(l::StackedLSTMCell)
	for i = 1:length(l.chain)
		l.chain[i].cell.b .= 1
	end
	return nothing
end

max_stars_repo_name: sdobber/FluxArchitectures

max_issues_repo_name: sdobber/FluxArchitectures

max_forks_repo_name: sdobber/FluxArchitectures


text:
function test_cashflow_irr()
    daycount = ActualActualISDA()
    evaluationdate = Date(2014, 6, 4)
    payments = [
        Payment(evaluationdate, -996.90)
        Payment(Date(2014, 12, 29), 531.5)
        Payment(Date(2015, 12, 29), 15.75)
        Payment(Date(2016, 12, 29), 15.75)
        Payment(Date(2017, 12, 29), 15.75)
        Payment(Date(2017, 12, 29), 500.0)
    ]
    cf = CashFlow(payments, daycount, evaluationdate)
    @test irr(cf) ≈ 0.040142375624945

    daycount = ActualActualISDA()
    evaluationdate = Date(2014, 12, 5)
    payments = [
        Payment(evaluationdate, -1007.572602739726)
        Payment(Date(2015, 12, 3), 31.5)
        Payment(Date(2016, 12, 3), 31.5)
        Payment(Date(2017, 12, 4), 31.5)
        Payment(Date(2018, 12, 3), 31.5)
        Payment(Date(2018, 12, 3), 1000.0)
    ]
    cf = CashFlow(payments, daycount, evaluationdate)
    @test irr(cf) ≈ 0.029507633486583
end

function test_cashflow_all()
    print(rpad("Finance.CashFlow...", 50, ' '))

    test_cashflow_irr()

    println("PASS")
end


text:
abstract type AbstractSymmetricMatrix{T,N,L} <: AbstractMatrix{T} end
abstract type MutableSymmetricMatrix{T,N,L} <: AbstractSymmetricMatrix{T,N,L} end

"""
Actually Hermitian, because adjoint returns itself.
"""
struct StaticSymmetricMatrix{T,N,L} <: AbstractSymmetricMatrix{T,N,L}
    data::NTuple{L,T}
end
mutable struct SymmetricMatrix{T,N,L} <: MutableSymmetricMatrix{T,N,L}
    data::NTuple{L,T}
    @generated SymmetricMatrix{T,N}() where {T,N} = :(SymmetricMatrix{$T,$N,$(N*N)}())
    function SymmetricMatrix{T,N,L}() where {T,N,L}
        isbits(T) || error("Can only construct mutable isbits matrices.")
        new{T,N,L}()
    end
    function SymmetricMatrix{T, N, L}(d::NTuple{L,T}) where {T,N,L}
        isbits(T) || error("Can only construct mutable isbits matrices.")
        new{T,N,L}(d)
    end
end
struct PointerSymmetricMatrix{T,N,L} <: MutableSymmetricMatrix{T,N,L}
    data::Ptr{T}
end
Base.size(::AbstractSymmetricMatrix{T,N}) where {T,N} = (N,N)

function StaticSymmetricMatrix(A::AbstractMatrix)
    m, n = size(A)
    @assert m == n
    StaticSymmetricMatrix(A, Val(n))
end
@generated function StaticSymmetricMatrix(A::AbstractMatrix{T}, ::Val{N}) where {T,N}
    q, L = upper_triangle_quote(N, T)
    push!(q.args, :(StaticSymmetricMatrix{$T,$N,$L}( @ntuple $L A )))
    q
end
function SymmetricMatrix(A::AbstractMatrix)
    m, n = size(A)
    @assert m == n
    SymmetricMatrix(A, Val(n))
end
@generated function SymmetricMatrix(A::AbstractMatrix{T}, ::Val{N}) where {T,N}
    q, L = upper_triangle_quote(N, T)
    push!(q.args, :(SymmetricMatrix{$T,$N,$L}( data )))
    q
end
function StaticSymmetricMatrix(data::AbstractVector)
    L = length(data)
    N = inv_triangle(L)
    StaticSymmetricMatrix(data, Val(N), Val(L))
end
function StaticSymmetricMatrix(data::AbstractVector{T}, ::Val{N}, ::Val{L}) where {T,N,L}
    StaticSymmetricMatrix{T,N,L}(ntuple(i -> data[i], L))
end
StaticSymmetricMatrix(data::NTuple{L,T}) where {T,L} = StaticSymmetricMatrix(data,ValI(Val{L}()))
SymmetricMatrix(data::NTuple{L,T}) where {T,L} = SymmetricMatrix(data,ValI(Val{L}()))
StaticSymmetricMatrix(data::SVector{L,T}) where {T,L} = StaticSymmetricMatrix(data.data,ValI(Val{L}()))
SymmetricMatrix(data::SVector{L,T}) where {T,L} = SymmetricMatrix(data.data,ValI(Val{L}()))
StaticSymmetricMatrix(data::NTuple{L,T}, ::Val{N}) where {T,N,L} = StaticSymmetricMatrix{T,N,L}(data)
SymmetricMatrix(data::NTuple{L,T}, ::Val{N}) where {T,N,L} = SymmetricMatrix{T,N,L}(data)
StaticSymmetricMatrix(data::SymmetricMatrix{T,N,L}) where {T,N,L} = StaticSymmetricMatrix{T,N,L}(data.data)
SymmetricMatrix(data::StaticSymmetricMatrix{T,N,L}) where {T,N,L} = SymmetricMatrix{T,N,L}(data.data)

point(A::MutableSymmetricMatrix{T}) where T = Base.unsafe_convert(Ptr{T}, pointer_from_objref(A))
point(A::PointerSymmetricMatrix) = A.data

@inline Base.getindex(A::StaticSymmetricMatrix, i::Integer) = A.data[i]
@inline Base.getindex(A::SymmetricMatrix, i::Integer) = A.data[i]
@inline Base.getindex(A::PointerSymmetricMatrix, i::Integer) = unsafe_load(A.data, i)
@inline function Base.getindex(A::StaticSymmetricMatrix{T,N,L}, i::Integer, j::Integer) where {T,N,L}
    @boundscheck begin
        max(i,j) > N && throw("BoundsError: index ($i, $j) is out of bounds.")
    end
    i, j = minmax(i, j)
    @inbounds out = A.data[small_triangle(j)+i]
    out
end

# @inline function Base.getindex(A::MutableSymmetricMatrix{T,N,L}, i::Integer, j::Integer) where {T,N,L}
#     @boundscheck begin
#         if (max(i,j) > 2) || (min(i,j) < 1)
#             throw(BoundsError("When recursing only indicies of 1 and 2 are supported. Received $i,$j."))
#         end
#     end
#     i, j = minmax(i, j)
#     unsafe_load(point(A.data), small_triangle(j)+i)
# end

# for N ∈ 1:cutoff
#     L = big_triangle(N)

#     @eval Base.size(A::MutableSymmetricMatrix{T,$N,$L}) = ($N,$N)

#     @eval @inline function Base.getindex(A::MutableSymmetricMatrix{T,$N,$L}, i::Integer, j::Integer) where {T,$N,$L}
#         @boundscheck begin
#             max(i,j) > $N && throw("BoundsError: index ($i, $j) is out of bounds.")
#         end
#         i, j = minmax(i, j)
#         unsafe_load(point(A.data), small_triangle(j)+i)
#     end

# end

"""
Note that Cartesian indexing is relatively slow, and results in branches even with `@inbounds`.
This means that `@simd`, for example, will not work while using Cartesian indices.
While linear indexing is fast, note that the underlying data is neither column nor row-major, so
be careful that your computations are actually doing what you intend.

You can use triangle_sub2ind (which incurs branching) to explore the behavior of
Cartesian -> Linear indices for Symmetric and Triangular matrices.
"""
@inline function Base.getindex(A::MutableSymmetricMatrix{T,N}, i::Integer, j::Integer) where {T,N}
    @boundscheck begin
        max(i,j) > N && throw("BoundsError: index ($i, $j) is out of bounds.")
    end
    i, j = minmax(i, j) # j >= i
    A[triangle_sub2ind(Val{N}(), i, j)]
end
@generated function Base.getindex(A::MutableSymmetricMatrix{T,N}, ::Val{1}, ::Val{1}) where {T, N}
    Nhalf = cld(N,2)
    triangle_size = big_triangle(Nhalf)
    :(PointerSymmetricMatrix{$T,$Nhalf,$triangle_size}(point(A)))
end
@generated function Base.getindex(A::MutableSymmetricMatrix{T,N}, ::Val{2}, ::Val{1}) where {T, N}
    Nremain, r = divrem(N,2)
    Nhalf = Nremain + r
    L = Nhalf * Nremain
    pointer_offset = big_triangle(Nhalf) * sizeof(T)
    :(LinearAlgebra.adjoint(PointerMatrix{$T,$Nhalf,$Nremain,$L}(point(A) + $pointer_offset )))
end
@generated function Base.getindex(A::MutableSymmetricMatrix{T,N}, ::Val{1}, ::Val{2}) where {T, N}
    Nremain, r = divrem(N,2)
    Nhalf = Nremain + r
    L = Nhalf * Nremain
    pointer_offset = big_triangle(Nhalf) * sizeof(T)
    :(PointerMatrix{$T,$Nhalf,$Nremain,$L}(point(A) + $pointer_offset ))
end
@generated function Base.getindex(A::MutableSymmetricMatrix{T,N,L}, ::Val{2}, ::Val{2}) where {T, N, L}
    Nremain = N ÷ 2
    triangle_size = big_triangle(Nremain)
    pointer_offset = (L - triangle_size) * sizeof(T)
    :(PointerSymmetricMatrix{$T,$Nremain,$triangle_size}(point(A) + $pointer_offset ))
end

@inline function Base.setindex!(A::MutableSymmetricMatrix{T,N,L}, val, i::Integer) where  {T,N,L}
    @boundscheck begin
        i > L && throw(BoundsError())
    end
    unsafe_store!(point(A), convert(T, val), i)
    return val
end
@inline function Base.setindex!(A::SymmetricMatrix{T,N,L}, val, i::Integer, j::Integer) where  {T,N,L}
    @boundscheck begin
        max(i,j) > N && throw(BoundsError())
    end
    i, j = minmax(i, j)
    unsafe_store!(point(A), convert(T, val), triangle_sub2ind(Val{N}(), i, j))
    return val
end

LinearAlgebra.transpose(x::SymmetricMatrix) = x
@static if VERSION >= v"0.7.0-DEV"
    LinearAlgebra.adjoint(x::SymmetricMatrix) = x
end


text:
import Base: AbstractSet, length, in, isempty, issubset, ==, push!, intersect, intersect!, union, union!, iterate, copy

mutable struct CustomSet{T} <: AbstractSet{T}
    values::AbstractArray{T}
    CustomSet(values::AbstractArray{T}) where {T} = new{T}(unique(values))
end

length(s::CustomSet) = length(s.values)

in(n, s::CustomSet) = in(n, s.values)

isempty(s::CustomSet) = isempty(s.values)

intersect(a::CustomSet, b::CustomSet) = CustomSet([x for x in a.values if x in b.values])

issubset(a::CustomSet, b::CustomSet) = intersect(a, b) == a

disjoint(a::CustomSet, b::CustomSet) = isempty(intersect(a, b))

==(a::CustomSet, b::CustomSet) = sort(a.values) == sort(b.values)

function push!(s::CustomSet, x)
    if !(x in s)
        push!(s.values, x)
    end
end

function intersect!(a::CustomSet, b::CustomSet)
    a.values = intersect(a, b).values
end

complement(a::CustomSet, b::CustomSet) = CustomSet([x for x in a.values if !(x in b.values)])

function complement!(a::CustomSet, b::CustomSet)
    a.values = complement(a, b).values
end

union(a::CustomSet, b::CustomSet) = CustomSet([a.values; b.values])

function union!(a::CustomSet, b::CustomSet)
    a.values = unique([a.values; b.values])
end

iterate(s::CustomSet, i...) = iterate(s.values, i...)

copy(s::CustomSet) = CustomSet([s.values...])

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
Pkg.add("JuMP")
Pkg.add("Gurobi")
Pkg.add("Complementarity")
Pkg.add("Taro")
Pkg.add("Plots")
Pkg.add("DataFrames")


text:
using Random
using DPMMSubClusters
using Glob
import YAML
using Pandas
include("run_experiments.jl")

# load config:
cfg = YAML.load_file("/home/vlad/projects/DPMMSubClusters.jl/benchmarks/configs/basic.yaml")

results_dir = ARGS[1]

# Extract params:
Ns = cfg["data"]["Ns"]
Ds = cfg["data"]["Ds"]
Ks = cfg["data"]["Ks"]
SEED = cfg["data"]["seed"]
DS_PATH = cfg["data"]["path"]
data_path = DS_PATH * "seed_" * string(SEED)
real_datasets = cfg["data"]["real_datasets"]

REPEATS = cfg["models"]["repeats"]
BURNOUT = cfg["models"]["burnout"]
ITERS = cfg["models"]["iterations"]
ALPHA = cfg["models"]["alpha"]

Random.seed!(SEED)

function load_data(path)
    data = HDF5.h5read(path, "data")
    gt = HDF5.h5read(path, "labels")
    return data, gt
end

# # Synthetic Data:

# iterate over number of points:
for l in 1:length(Ns)
    N = Ns[l]

    # iterate over dimensions:
    for i in 1:length(Ds)
        D = Ds[i]

        # iterate over number of clusters
        for j in 1:length(Ks)
            K = Ks[j]

            filename = "*D=" * string(D) * "*K=" * string(K) *"*.h5"

            # get path:
            datapath = glob(filename, data_path)[1]
            print("loading from: ")
            println(datapath)

            # load data:
            x, y = load_data(datapath)

            # run methods:

            # Random:
            # check if results exist:
            csv_results_random = "DPMM_D=" * string(D) * "_K=" * string(K) * "_results_random.csv"
            if length(glob(csv_results_random, data_path * "/results/")) > 0
                println("Found existing results for: " * csv_results_random)
            else
                println("No existing results for: " * csv_results_random * ", running inference:")

                res_random = repeat_experiment_xy(x, y, ALPHA, REPEATS,"none", BURNOUT, ITERS)
                # save results
                df_random = DataFrame(res_random)
                df_random["iter"] = 1:length(df_random)
                df_random["method"] = "random"

                save_path =  data_path * "/results/" * csv_results_random
                println("Saving results at: " * save_path)
                to_csv(df_random, save_path)
            end

            # kmeans
            csv_results_kmeans = "DPMM_D=" * string(D) * "_K=" * string(K) * "_results_kmeans.csv"
            if length(glob(csv_results_kmeans, data_path * "/results/")) > 0
                println("Found existing results for: " * csv_results_kmeans)
            else
                println("No existing results for: " * csv_results_kmeans * ", running inference:")

                res_kmeans = repeat_experiment_xy(x, y, ALPHA, REPEATS,"kmeans", BURNOUT, ITERS)
                # save results

                df_kmeans = DataFrame(res_kmeans)
                df_kmeans["iter"] = 1:length(df_kmeans)
                df_kmeans["method"] = "kmeans"

                save_path =  data_path * "/results/" * csv_results_kmeans
                println("Saving results at: " * save_path)
                to_csv(df_kmeans, save_path)
            end

            # Splitnet
            if D == 2
                res_splitnet = repeat_experiment_xy(x, y, ALPHA, REPEATS, "splitnet_2d", BURNOUT, ITERS)
            elseif D == 10
                res_splitnet = repeat_experiment_xy(x, y, ALPHA, REPEATS, "splitnet_10d", BURNOUT, ITERS)
            elseif D == 128
                res_splitnet = repeat_experiment_xy(x, y, ALPHA, REPEATS, "splitnet_128d", BURNOUT, ITERS)
            end

            df_splitnet = DataFrame(res_splitnet)
            df_splitnet["iter"] = 1:length(df_splitnet)
            df_splitnet["method"] = "splitnet"

            csv_save_path = results_dir * "/DPMM_D=" * string(D) * "_K=" * string(K) * "_results_splitnet.csv"
            to_csv(df_splitnet, csv_save_path)
            print("Julia: Saved results to: ")
            println(csv_save_path)

        end
    end
end

# "Real" Datasets:
for i in 1:length(real_datasets)
    ds_type = real_datasets[i]
    ds_name = ds_type["ds_name"]
    ds_path_r = DS_PATH * "/real_datasets/" * ds_type["path"]

    # load data:
    x, y = load_data(ds_path_r)

    # Random:
    # check if results exist:
    csv_results_random = "DPMM_" * string(ds_name) * "_results_random.csv"
    if length(glob(csv_results_random, DS_PATH * "/real_datasets/results/")) > 0
        println("Found existing results for: " * csv_results_random)
    else
        println("No existing results for: " * csv_results_random * ", running inference:")

        res_random = repeat_experiment_xy(x, y, ALPHA, 1,"none", 4, 10)
        # res_random = repeat_experiment_xy(x, y, ALPHA, REPEATS,"none", BURNOUT, ITERS)
        # save results
        df_random = DataFrame(res_random)
        df_random["iter"] = 1:length(df_random)
        df_random["method"] = "random"

        save_path =  DS_PATH * "/real_datasets/results/" * csv_results_random
        println("Saving results at: " * save_path)
        to_csv(df_random, save_path)
    end

    # kmeans
    csv_results_kmeans = "DPMM_" * string(ds_name) * "_results_kmeans.csv"
    if length(glob(csv_results_kmeans, DS_PATH * "/real_datasets/results/")) > 0
        println("Found existing results for: " * csv_results_kmeans)
    else
        println("No existing results for: " * csv_results_kmeans * ", running inference:")

        res_kmeans = repeat_experiment_xy(x, y, ALPHA, 1,"none", 4, 10)
        # res_kmeans = repeat_experiment_xy(x, y, ALPHA, REPEATS,"kmeans", BURNOUT, ITERS)

        # save results
        df_kmeans = DataFrame(res_kmeans)
        df_kmeans["iter"] = 1:length(df_kmeans)
        df_kmeans["method"] = "kmeans"

        save_path =  DS_PATH * "/real_datasets/results/" * csv_results_kmeans
        println("Saving results at: " * save_path)
        to_csv(df_kmeans, save_path)
    end

    # it's always D=128
    res_splitnet = repeat_experiment_xy(x, y, ALPHA, 1,"splitnet_128d", 4, 10)

    # save results
    df_splitnet = DataFrame(res_splitnet)
    df_splitnet["iter"] = 1:length(df_splitnet)
    df_splitnet["method"] = "splitnet"

    save_path =  results_dir * "/DPMM_" * string(ds_name) * "_results_splitnet.csv"
    println("Saving results at: " * save_path)
    to_csv(df_splitnet, save_path)

end

# REPEATS=1
# BURNOUT=5

# ITERS=40
# # ITERS=5
# Ks = 5
# x,y,clusters = generate_gaussian_data(10^5,2,Ks,80.0)

# run(`rm -rf $dir`)

# run(`mkdir $dir`)
# # labels,clusters,weights = fit(x,50.0, iters = ITERS, seed = SEED, gt=y, init_type="splitnet-gpu", burnout=BURNOUT)
# gauss2d_splitnet = repeat_experiment_xy(x,y,REPEATS,"splitnet_2d",BURNOUT, ITERS)
# # run(`convert -delay 50 -loop 0 ./test_dir/clusters_*.png ./compare_methods_gifs/clusters_K-"$Ks"_B-"$BURNOUT"_seed_"$SEED"_splitnet.gif`)
# run(`rm -rf $dir`)

# run(`mkdir $dir`)
# gauss2d_random = repeat_experiment_xy(x,y,REPEATS,"none",BURNOUT, ITERS)
# # labels,clusters,weights = DPMMSubClusters.fit(x,100.0, iters = ITERS, seed = SEED, gt=y, init_type="none", burnout=BURNOUT)
# # run(`convert -delay 50 -loop 0 ./test_dir/clusters_*.png ./compare_methods_gifs/clusters_K-"$Ks"_B-"$BURNOUT"_seed_"$SEED"_random.gif`)
# run(`rm -rf $dir`)

# run(`mkdir $dir`)
# # # labels,clusters,weights = fit(x,50.0, iters = ITERS, seed = SEED, gt=y, init_type="kmeans", burnout=BURNOUT)
# gauss2d_kmeans = repeat_experiment_xy(x,y,REPEATS,"kmeans",BURNOUT, ITERS)
# # run(`convert -delay 50 -loop 0 ./test_dir/clusters_*.png ./compare_methods_gifs/clusters_K-"$Ks"_B-"$BURNOUT"_seed_"$SEED"_kmeans.gif`)
# run(`rm -rf $dir`)

# all_results = Dict(
#     "rabdom" => gauss2d_random,
#     "kmeans" => gauss2d_kmeans,
#     "splitnet_gpu" => gauss2d_splitnet,
# )

# title = "K-" * string(Ks) * "_Burnout-" * string(BURNOUT) * "seed-" * string(SEED) * "_model_"*string(model)

# fig = plot_nmi_for_iters(all_results, title)
# savefig("./compare_methods_gifs/result_"* string(title) *".png")

# ITERS=40
# Ks = 10
# x,y,clusters = generate_gaussian_data(10^5,2,Ks,80.0)

# run(`mkdir $dir`)
# # labels,clusters,weights = fit(x,50.0, iters = ITERS, seed = SEED, gt=y, init_type="splitnet-gpu", burnout=BURNOUT)
# gauss2d_splitnet = repeat_experiment_xy(x,y,REPEATS,"splitnet-gpu",BURNOUT, ITERS)
# run(`convert -delay 50 -loop 0 ./test_dir/clusters_*.png ./compare_methods_gifs/clusters_K-"$Ks"_B-"$BURNOUT"_seed_"$SEED"_splitnet.gif`)
# run(`rm -rf $dir`)

# run(`mkdir $dir`)
# gauss2d_random = repeat_experiment_xy(x,y,REPEATS,"none",BURNOUT, ITERS)
# # labels,clusters,weights = DPMMSubClusters.fit(x,100.0, iters = ITERS, seed = SEED, gt=y, init_type="none", burnout=BURNOUT)
# run(`convert -delay 50 -loop 0 ./test_dir/clusters_*.png ./compare_methods_gifs/clusters_K-"$Ks"_B-"$BURNOUT"_seed_"$SEED"_random.gif`)
# run(`rm -rf $dir`)

# run(`mkdir $dir`)
# # # labels,clusters,weights = fit(x,50.0, iters = ITERS, seed = SEED, gt=y, init_type="kmeans", burnout=BURNOUT)
# gauss2d_kmeans = repeat_experiment_xy(x,y,REPEATS,"kmeans",BURNOUT, ITERS)
# run(`convert -delay 50 -loop 0 ./test_dir/clusters_*.png ./compare_methods_gifs/clusters_K-"$Ks"_B-"$BURNOUT"_seed_"$SEED"_kmeans.gif`)
# run(`rm -rf $dir`)

# all_results = Dict(
#     "rabdom" => gauss2d_random,
#     "kmeans" => gauss2d_kmeans,
#     "splitnet_gpu" => gauss2d_splitnet,
# )

# title = "K-" * string(Ks) * "_Burnout-" * string(BURNOUT) * "seed-" * string(SEED) * "_model_"*string(model)

# fig = plot_nmi_for_iters(all_results, title)
# savefig("./compare_methods_gifs/result_"* string(title) *".png")

# ITERS=100
# Ks = 40
# x,y,clusters = generate_gaussian_data(2*10^5,2,Ks,160.0)

# run(`mkdir $dir`)
# # labels,clusters,weights = fit(x,50.0, iters = ITERS, seed = SEED, gt=y, init_type="splitnet-gpu", burnout=BURNOUT)
# gauss2d_splitnet = repeat_experiment_xy(x,y,REPEATS,"splitnet-gpu",BURNOUT, ITERS)
# run(`convert -delay 50 -loop 0 ./test_dir/clusters_*.png ./compare_methods_gifs/clusters_K-"$Ks"_B-"$BURNOUT"_seed_"$SEED"_splitnet.gif`)
# run(`rm -rf $dir`)

# run(`mkdir $dir`)
# gauss2d_random = repeat_experiment_xy(x,y,REPEATS,"none",BURNOUT, ITERS)
# # labels,clusters,weights = DPMMSubClusters.fit(x,100.0, iters = ITERS, seed = SEED, gt=y, init_type="none", burnout=BURNOUT)
# run(`convert -delay 50 -loop 0 ./test_dir/clusters_*.png ./compare_methods_gifs/clusters_K-"$Ks"_B-"$BURNOUT"_seed_"$SEED"_random.gif`)
# run(`rm -rf $dir`)

# run(`mkdir $dir`)
# # # labels,clusters,weights = fit(x,50.0, iters = ITERS, seed = SEED, gt=y, init_type="kmeans", burnout=BURNOUT)
# gauss2d_kmeans = repeat_experiment_xy(x,y,REPEATS,"kmeans",BURNOUT, ITERS)
# run(`convert -delay 50 -loop 0 ./test_dir/clusters_*.png ./compare_methods_gifs/clusters_K-"$Ks"_B-"$BURNOUT"_seed_"$SEED"_kmeans.gif`)
# run(`rm -rf $dir`)

# all_results = Dict(
#     "rabdom" => gauss2d_random,
#     "kmeans" => gauss2d_kmeans,
#     "splitnet_gpu" => gauss2d_splitnet,
# )

# title = "K-" * string(Ks) * "_Burnout-" * string(BURNOUT) * "seed-" * string(SEED) * "_model_"*string(model)

# fig = plot_nmi_for_iters(all_results, title)
# savefig("./compare_methods_gifs/result_"* string(title) *".png")


text:
@testset "Controller: Implicit Dynamics" begin
	T = Float64
	s = get_simulation("particle", "flat_3D_lc", "flat_lc")

	model = s.model
	env = s.env

	nq = model.dim.q
	nu = model.dim.u
	H = 10
	h = 0.0001
	κ = 1e-6
	nq = model.dim.q

	q0 = SVector{nq,T}([0.0, 0.0, 1.0])
	q1 = SVector{nq,T}([0.0, 0.0, 1.0])

	ip_opts = ContactImplicitMPC.InteriorPointOptions(
		κ_tol = 1.0e-8, r_tol = 1.0e-10)

	sim = ContactImplicitMPC.simulator(s, q0, q1, h, H;
		p = ContactImplicitMPC.no_policy(model),
		d = ContactImplicitMPC.no_disturbances(model),
		ip_opts = ip_opts,
		sim_opts = ContactImplicitMPC.SimulatorOptions{T}())

	ContactImplicitMPC.simulate!(sim, verbose = false)
	ref_traj = deepcopy(sim.traj)
	ref_traj.κ .= κ # sets κ_tol for the im_traj IP solvers
	traj0 = deepcopy(ref_traj)

	im_traj = ContactImplicitMPC.ImplicitTraj(ref_traj, s, mode = :configuration)
	ContactImplicitMPC.implicit_dynamics!(im_traj, s, traj0)

	# Implicit dynamics contraint violation is ≈ 0
	for t = 1:H
		@test norm(im_traj.d[t], Inf) < 1.0e-5
	end

	# Quadruped reference trajectory
	s = ContactImplicitMPC.get_simulation("quadruped", "flat_2D_lc", "flat")
	model = s.model
	env = s.env

	ref_traj = deepcopy(ContactImplicitMPC.get_trajectory(model, env,
	    joinpath(module_dir(), "src/dynamics/quadruped/gaits/gait2.jld2"),
	    load_type = :split_traj_alt))

	ip_opts = eval(ContactImplicitMPC.interior_point_options(:interior_point))(
				κ_tol = 2.0 * 1.0e-4,
				r_tol = 1.0e-8,
				diff_sol = true,
				solver = :empty_solver)

	im_traj = ContactImplicitMPC.ImplicitTraj(ref_traj, s,
		ip_type = :interior_point,
		κ = 1.0e-4,
		mode = :configuration,
		opts=ip_opts)

	ContactImplicitMPC.implicit_dynamics!(im_traj, s, ref_traj)

	@test all([norm(im_traj.dq2[t], Inf) < 1.0e-2 for t = 1:ref_traj.H-1])
end


text:
function hamiltonian_monte_carlo(sel::K, cl::C; L=10, eps=0.1) where {K <: UnconstrainedSelection, C <: CallSite}
    local u_cl

    p_mod_score = get_score(cl)
    sel_values, choice_grads = get_choice_gradients(sel, cl, 1.0)
    vals = array(sel_values, Float64)
    grads = array(choice_grads, Float64)
    d = MvNormal(length(vals), 1.0)
    mom = rand(d)
    p_mom_score = logpdf(d, mom)
    for step in 1 : L
        mom += (eps / 2) * grads
        vals += eps * mom
        sel_values = selection(sel_values, vals)
        ret, u_cl, w, _ = update(sel_values, cl)
        _, choice_grads = get_choice_gradients(sel, u_cl, 1.0)
        grads = array(choice_grads, Float64)
        mom += (eps / 2) * grads
    end
    n_mod_score = get_score(u_cl)
    n_mom_score = logpdf(d, -mom)
    alpha = get_score(u_cl) - p_mod_score + n_mom_score - p_mom_score

    # Accept or reject.
    if log(rand()) < alpha
        (u_cl, true)
    else
        (cl, false)
    end
end


text:
function circles(h,k,r)
	# see this for more information: https://discourse.julialang.org/t/plot-a-circle-with-a-given-radius-with-plots-jl/23295/5
theta = collect(range(0,2*π,length = 500));
return h .+ r*sin.(theta), k .+ r*cos.(theta);
end


text:
# read in arc data

type Arc
    initNode::Int
    termNode::Int
    capacity::Float64
    freeflowtime::Float64
    flow::Float64
    flow_car::Float64
    flow_truck::Float64
end

Arc(initNode::Int, termNode::Int, capacity::Float64,freeflowtime::Float64) =
    Arc(initNode, termNode, capacity, freeflowtime, 0., 0., 0.)

function arcData(arc_file)
    arcs = Dict()
    file = open(arc_file)
    inHeader=true
    for line in eachline(file)
        if inHeader
            inHeader = !contains(line, "Init node")
            continue
        end
        vals = split(line, )
        arcs[(parse(Int, vals[1]), parse(Int, vals[2]))] = Arc(parse(Int, vals[1]), parse(Int, vals[2]),
        parse(Float64, vals[3]), parse(Float64, vals[5]))
    end
    close(file)
    return arcs
end

# add flow data to arcs
function observFlow(arc_file, tapFlowDic)
    arcs = arcData(arc_file)
    ix = 0
    for key in keys(arcs)
        arcs[key].flow = 1.0 * tapFlowDic["car"][key] + 2.0 * tapFlowDic["truck"][key]
        arcs[key].flow_car = tapFlowDic["car"][key]
        arcs[key].flow_truck = tapFlowDic["truck"][key]
    end
    return arcs
end

# obtain important parameters of the network

include("../../Julia_files/load_network_uni_class.jl")

function paraNetwork(nameNetwork)
    ta_data = load_ta_network(nameNetwork)
    numNodes = maximum(map(pair->pair[1], keys(demandsDict[0]["car"])))
    start_node = ta_data.start_node
    capacity = ta_data.capacity
    free_flow_time = ta_data.free_flow_time
    numLinks = size(start_node)[1]
    numODpairs = numNodes * (numNodes - 1)
    return numNodes, numLinks, numODpairs, capacity, free_flow_time, ta_data
end

function tapFlowVecToLinkCostDict(tapFlowVec, fcoeffsInvVI)
    linkCostVecCar = 1.0 * BPR(1.0 * tapFlowVec[1,:] + 2.0 * tapFlowVec[2,:], fcoeffsInvVI)
    linkCostVecTruck = 1.1 * BPR(1.0 * tapFlowVec[1,:] + 2.0 * tapFlowVec[2,:], fcoeffsInvVI)

    temp_dict_car = Dict{}()
    temp_dict_truck = Dict{}()
    for i = 1:length(linkCostVecCar)
        temp_dict_car["$(i-1)"] = linkCostVecCar[i]
        temp_dict_truck["$(i-1)"] = linkCostVecTruck[i]
    end

    linkCostDict = Dict{}()
    linkCostDict["car"] = temp_dict_car
    linkCostDict["truck"] = temp_dict_truck

    return linkCostDict
end

max_stars_repo_name: jingzbu/InverseVITraffic

max_issues_repo_name: jingzbu/InverseVITraffic

max_forks_repo_name: jingzbu/InverseVITraffic


text:
### struct for defining the Muller-Brown potential
mutable struct CustomModel <: AbstractSinglePotential
  dim::Int64
  pfunc::Union{Function,nothing}
  ffunc::Union{Function,nothing}
  f::Vector{Float64}
  function CustomModel(dim::Int64,pfunc::Union{Function,nothing},ffunc::Union{Function,nothing})
    new(dim,pfunc,ffunc2,zeros(Float64,dim))
  end
end
CustomModel(dim::Int64;pfunc::Union{Function,nothing}=nothing,ffunc::Union{Function,nothing}=nothing) = CustomModel(dim,pfunc,ffunc)

function getdimensionality(model::CustomModel)
  return model.dim
end

"""
Potential functions
"""
function potential(system::AbstractSystem, model::CustomModel)
  if model.pfunc != nothing
    return model.pfunc(system)
  end
end

"""
Force functions with return value
"""
function force(system::AbstractSystem, model::CustomModel)
  if model.ffunc != nothing
    return model.ffunc(system)
  else
    return zeros(Float64,model.dim)
  end
end

"""
Force functions with in-place changing
"""
function force!(system::AbstractSystem, model::CustomModel)
  if model.ffunc != nothing
    model.ffunc(system,model.f)
  else
    modelf .= 0.0
  end
end


text:
using Dance.Router

route("/../files/image.jpg", output_file_as_string; method=GET, endpoint=EP_STATIC)
route("/files/image.jpg", output_file_as_string; method=GET, endpoint=EP_STATIC)
static_dir("/static", "../sample/static")


text:
using DeformableRegistration: ImageProcessing, Distance, Transformation, Interpolation, Examples, Distance, regOptions
 using Base.Test
 using Logging

@testset "Parametric Registration" begin
@testset "SSD/NGF shifted and scaled"  for D=[ssdDistance, ngfDistance]
##
    testimage = dirname(Base.source_path()) * "/testdata/luebeck.jpg"
    refImg = loadImage(testimage)
    # define a cell centered grid, transform it and create a template image
    centeredGrid = getCellCenteredGrid(refImg)
    affineParametersInitial = [1.2,0.0,-50,0.0,1.2,-0]
    transformedGrid = transformGridAffine(centeredGrid,affineParametersInitial)
    temImg = interpolateImage(refImg,transformedGrid,interpolationScheme=InterpLinearFast)
    temImg = createImage(temImg)
    options = regOptions()
    options.levels = [7,6,5,4]
    options.parametricOnly = true;
    options.matrixFree = true;
    options.useEdgeParameterInNumerator

    affineParameters = registerImagesParametric(refImg,temImg, options, measureDistance=D)
    @test affineParameters[1] ≈ 1/affineParametersInitial[1] atol=0.1
    @test affineParameters[2] ≈ 0  atol=0.1
    @test affineParameters[3] ≈ -affineParametersInitial[3]/affineParametersInitial[1]  atol=2
    @test affineParameters[4] ≈ 0  atol=0.1
    @test affineParameters[5] ≈ 1/affineParametersInitial[5]  atol=1

end
end


text:
InputFolder = './Images/Fibers/';
OutputFolder = './Results/Images/Fibers/';
Fill = -1;

@iA = '*.tif';		% Image filter

@fxg_gmLocOrientation[iA] > [O, M];
params.GRad = 1;
params.CellSize = 8;
params.MinMag = 75;
params.Vecs = 2;
/endf

/show iA > M;
/keep O > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER