text:
using NLopt
using ForwardDiff

function w_next(c, h, states, agg_state, 𝔼V, params, hh::Owner, s::NoState=hh.state)
  @unpack p, r = agg_state
  @unpack δ = params
  @unpack z, w = states

  w_next = w + z - c - p * h * (r + δ)
end

function w_next(c, h_next, states, agg_state, 𝔼V, params, hh::Owner, s::IsState=hh.state)
  @unpack p, r = agg_state
  @unpack δ = params
  @unpack z, w, h = states

  w_next = w + z - c - p * h_next * (r + δ) #+ 0.1 * p * h * (h_next != h)
end

function m_next(c, h, states, agg_state, 𝔼V, params, hh::Owner)
  @unpack p = agg_state
  @unpack z, w = states

  m = p * h + c - z - w
end

function objective0(c, h, states, agg_state, 𝔼V, params, hh::Owner, ::NoState)
  @unpack β = params

  w_next_ = w_next(c, h, states, agg_state, 𝔼V, params, hh, hh.state)

  u(c,h) + β * 𝔼V(w_next_)
end

function objective0(c, h, states, agg_state, 𝔼V, params, hh::Owner, ::IsState)
  @unpack β = params

  w_next_ = w_next(c, h, states, agg_state, 𝔼V, params, hh, hh.state)

  u(c,h) + β * 𝔼V(w_next_, h)
end

# function objective0(c, h, states, agg_state, 𝔼V, params, hh::Owner{<:Aiyagari.Conditional})
#   @unpack β = params
#
#   w_next_ = w_next(c, h, states, agg_state, 𝔼V, params, hh)
#
#   u(c,h) + β * 𝔼V([w_next_, w_next_, w_next_ - 0.8 * agg_state.p * h])
# end

function constraint0(c, h, states, agg_state, 𝔼V, params, hh::Owner)
  @unpack p, r = agg_state
  @unpack β, θ, δ = params

  m = m_next(c, h, states, agg_state, 𝔼V, params, hh::Owner)

  #(1+r) * m <= p * h * (1-δ) * θ
  (1+r) * m - (p * h * (1-δ) * θ)
end

function objective_nlopt(x::Vector, grad::Vector, args...)
  if length(grad) > 0
    ForwardDiff.gradient!(grad, XX -> objective0(XX[1], XX[2], args...), x)
  end
  objective0(x[1], x[2], args...)
end

function constraint_nlopt(x::Vector, grad::Vector, args...)
  if length(grad) > 0
    ForwardDiff.gradient!(grad, xx -> constraint0(xx[1], xx[2], args...), x)
  end
  constraint0(x[1], x[2], args...)
end

function Aiyagari.get_optimum(states, agg_state, 𝔼V, params, a_grid, hh::Owner)
  @unpack h_thres = params

  # 1. check if feasible set is non-empty
  h_max = let
    @unpack z, w = states
    @unpack δ, θ = params
    @unpack p, r = agg_state

    (w + z) / (p * (1 - (1-δ) * θ / (1+r)))
  end

  if h_max < h_thres
    conv = true
    pol_full = (c=NaN, h=NaN, m=NaN, w_next=NaN, ret=:infeasible, conv=conv, count=0)
    val = -Inf

  else
    opt = Opt(:LD_MMA, 2)
    #opt = Opt(:LD_SLSQP, 2)
    lower_bounds!(opt, [eps(), params.h_thres])

    xtol_rel!(opt, 1e-10)
    ftol_rel!(opt, 1e-10)

    maxeval!(opt, 250)

    max_objective!(opt, (x,g) -> objective_nlopt(x, g, states, agg_state, 𝔼V, params, hh, hh.state))
    inequality_constraint!(opt, (x,g) -> constraint_nlopt(x, g, states, agg_state, 𝔼V, params, hh), 1e-12)

    guess = (states.w + states.z)/2
    #guesses = [guess, guess / agg_state.p]
    guesses = [guess/10, max(guess / agg_state.p, params.h_thres)]
    (max_f, max_x, ret) = NLopt.optimize(opt, guesses)

    val = max_f
    c, h = max_x
    m = m_next(c, h, states, agg_state, 𝔼V, params, hh::Owner)
    w_ = w_next(c, h, states, agg_state, 𝔼V, params, hh::Owner)

    conv = ret in [:FTOL_REACHED, :XTOL_REACHED, :SUCCESS, :LOCALLY_SOLVED]

    pol_full = (c=c, h=h, m=m, w_next=w_, ret=ret, conv=conv, count= opt.numevals)

  end
  pol = pol_full.w_next, pol_full.h

  (pol=pol, pol_full=pol_full, val=val, conv= conv)
end

#state1 =  (a = 0.01, z = 0.5)
#grad1 = [0.0, 0.0]
#objective_nlopt([sum(state1)/2, sum(state1)/(2 * agg_state.p)], grad1, state1, agg_state, log, param)
#constraint_nlopt([sum(state1)/2, sum(state1)/(2 * agg_state.p)], grad1, state1, agg_state, log, param)
#@show grad1


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

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

# ╔═╡ 0b4c1052-9180-11ec-2cdd-1b61bde8bdc6
begin
using Plots
using Printf
using Unitful
using DataFrames
using CSV
using Images
end

# ╔═╡ 2d15cff8-622e-4f60-a16f-77a4d16ccc0d
begin
using Interpolations
using QuadGK
using UnitfulEquivalences
end

# ╔═╡ ffc498a9-3522-4309-9ced-48955cab1759
begin
using Markdown
using InteractiveUtils
using PlutoUI
end

# ╔═╡ 91aa2620-9728-4de7-9c71-f01a6c765350
begin
using LsqFit
using Distributions
end

# ╔═╡ 9d66c472-b742-49b3-9495-1c3be998108b
import Unitful:
    nm, μm, mm, cm, m, km,
    mg, g, kg,
    ps, ns, μs, ms, s, minute, hr, d, yr, Hz, kHz, MHz, GHz,
    eV,
    μJ, mJ, J,
	μW, mW, W,
    A, N, mol, mmol, V, L, M

# ╔═╡ c3b89fa7-67f8-4cc3-b3ba-fcb304f44669
import PhysicalConstants.CODATA2018: N_A

# ╔═╡ bd9eed27-fc29-41ff-bf12-ec0cd3881ea2
function ingredients(path::String)
	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object
	name = Symbol(basename(path))
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	m
end

# ╔═╡ ca6f7233-a1d1-403d-9416-835a51d360db
lfi = ingredients("../src/LaserLab.jl")

# ╔═╡ 4c6c9627-d80d-41ac-8141-f935df893b73
md"# Observing fluorescence in a ML of Ru-sil molecules (RUSL)"

# ╔═╡ 34b3d66e-3b39-43dc-922b-9ba6c15c1bac
md"## The RuSL molecules

RuSL molecules are a type of phosphorescente molecules which can be used to calibrate the TOPATU laser setup.

The molecules emit ligh peaked in the red, with a lifetime in the range of 500 ns.
"

# ╔═╡ 7e150a69-0070-4371-b5c0-02b7ad70d813
md"### Fluorescence cross section and quantum yield"

# ╔═╡ 16cdbadf-4deb-4a66-8ab7-84437a4fe3d4
load("../notebooks/img/RuSlAbs.png")

# ╔═╡ 9cc3dfcb-b95d-4086-a228-ed4753f6ca0d
begin
	ϵabs = 8801.5/(M*cm)
	Q    = 0.9
	λexc = 485.0nm
	λEM  = 690.0nm
	println("")
end

# ╔═╡ 4a0e30fe-398e-4d80-86a2-bfc384cbc6e6
md"
- The absorption cross section at 469 nm is ϵ = $ϵabs
- Measured quantum yield is Q= $Q
- The excitation wavelength of the laser is $λexc
"

# ╔═╡ 591f9fcb-7ad1-400e-af0a-29686a4914ec
md"### Emission spectrum on ML

- ML of RuSL, at a nominal packing of 1 molecule per nm
- Peak emission around $λEM
"

# ╔═╡ 768feb2a-855f-4db6-9ae0-6abec23707c1
begin
	spG = lfi.LaserLab.CsvG(';',',')  # spanish: decimals represented with ',' delimited with ';'
	enG = lfi.LaserLab.CsvG(',','.')  # english
	enG2 = lfi.LaserLab.CsvG('\t','.')  # english
	println("")
end

# ╔═╡ 9cc092e9-40d2-4c5a-9b55-4e2adccb3382
begin
	dfrusl = lfi.LaserLab.load_df_from_csv("/Users/jj/JuliaProjects/LaserLab/data/fluorimeter", "RuSL_quartz.csv", spG)
	dfrusl[!, :λ] = convert.(Float64, dfrusl[:, :λ])
	first(dfrusl,10)
	println("")
end

# ╔═╡ d75a8426-c97b-422e-8c69-3f3ce94c5370
begin
plot(dfrusl.λ, dfrusl.QUARTZ_Rusilatrane_A, lw=2, label="RuSL on quartz")
xlabel!("λ (nm)")
ylabel!("I (a.u.)")
end

# ╔═╡ a7d330cf-5093-490f-adc8-e482e3084806
md"## Temporal dependence of the phosoprescence"

# ╔═╡ 280a17c4-6d5e-4ff3-a5ab-dbbcb7199bb2
begin
	dfruslt = lfi.LaserLab.load_df_from_csv("/Users/jj/JuliaProjects/LaserLab/data/fluorimeter", "Ru_SL_time.csv", spG)
	dfruslt[!, :λ] = convert.(Float64, dfruslt[:, :λ])
	first(dfruslt,10)
	println("")
end

# ╔═╡ a230e061-b6ab-49d4-bdad-730d75e20e9c
md"### The data fits (but not too well) to a single exponential

- At large times the fluorimeter may be measuring a constant pedestal
  "

# ╔═╡ 6d5e295c-1c57-4a92-aee2-3ffcbb3306df
begin
	mexp(t, p) = p[1] * exp.(-t/p[2])
	pa0 = [3500.0, 500.0]
	tdata = dfruslt.λ
	ydata= dfruslt.Ru_Quartz
	fit = curve_fit(mexp, tdata, ydata, pa0)
	cofe = coef(fit)
	stder = stderror(fit)
	println("")
end

# ╔═╡ 80447dc7-1103-4019-bc64-283d4a9368a0
@info "fit coefficients" cofe

# ╔═╡ 08ec5b69-946f-407e-ac8d-bf8814cc0121
@info "coefficient errors (std)" stder

# ╔═╡ 58b227ac-627f-49ba-bc58-75039c65733b
md"## The RUSL experiment

The goal of the RUSL experiment is to use the (modified) TOPATU setup to measure the spectrum (color) and temporal response of the RuSL ML on quartz. This experiment has already been performed in the fluorimeter (see results below), and the main purpose of RUSL is to repeat it with the laser setup to calibrate the system (including temporal response)

The experiment ingredients are:

- A ML of RuSL on quartz.
- A pulsed laser of 480 nm (EPL 480 from Edimburgh).
- The topatu setup, modified to allow measurement of spectrum (with filters) and time response.
- The measurement must include spectral response (number of photons observed as a function of the wavelength) and time response (which requires time-stamps for the observed photons)
"

# ╔═╡ 4223f9f8-8ee6-4ea9-a7bc-6379c81c48c0
load("../notebooks/img/Nuevo_SET_UP_RuSl.png")

# ╔═╡ fabba61a-f28d-4535-b0e3-e94f5098bb0a
md"## The Monolayers

A Monolayer (ML) of fluorophores is characterised by the packing of the fluorophores, the absorption cross section of the individual fluorophores and their emission spectra. These properties may also be correlated or can be modified when moving from solution to ML. For example colective effects of the interference with substrate can modify the naive assumption that the total fluorescence is the product of the number of fluorophores and the fluorescence per fluorophore. However, it is useful to start with the simplest assumption that all cross sections measured in solution hold in solid/gas interface

"

# ╔═╡ 7a1fcff4-6948-4333-b971-8fbf31c7d3ee
md"#### Define the fluorophores at nominal excitation ($λexc)"

# ╔═╡ a917927d-d471-4465-abd1-30d959af0b45
frusl = lfi.LaserLab.Fluorophore(λexc, λEM, ϵabs, Q)

# ╔═╡ aa4283d7-37d4-4aaf-a5b7-b515466e545d
begin
	λepl    = 485.0nm
	Pkepl   = 35.0mW
	P200khz = 0.6μW
	fepl    = 200.0kHz
	wepl    = 140.0ps
	dc      = uconvert(μs, 1.0/fepl)
end

# ╔═╡ 93b78d1a-ee07-4d66-b5a4-c7a109a24f81
md"## Laser
- The experiment requires a VUV pulsed laser. For the experiment we will use a repetition rate of $fepl (pulsed each 1 μs)
-
- The nominal laser for the experiment is the EPL485

| λepl   | Pkepl      | fepl  | wepl | dc| P|
|:-------| ---------- |:-----:|:-----:|:-----:|:-----:|
| $λepl  | $Pkepl| $fepl | $wepl |$dc| $P200khz

"

# ╔═╡ 19b7d7fb-adab-4f29-82af-d059525921ee
md"#### Define laser"

# ╔═╡ 102a2054-2fd5-406e-8bb4-20ecb47c278f
epl485 = lfi.LaserLab.PulsedLaser(λepl, Pkepl, fepl, P200khz, wepl)

# ╔═╡ e97474ef-145f-45c8-96f6-213acf4a3b41
md"### Beam shape
- The beam has an oval shape, with a long axis of 3.5 mm and short axis of 1.5 mm
- The entrance iris diameter of the objective is 2.5 mm. The beam overfills one of the dimensions and not quite the other.
- We will make the approch that the beam fills the entrance iris
"

# ╔═╡ 106f7063-719e-4eff-857d-3566d6e6d4c8
load("../notebooks/img/beamspot.png")

# ╔═╡ 177762d5-b9f2-449b-abba-256f3d4318ad
md"## Objective

The objective directs the laser light into the ML. For the first round of experiments, it is convenient to focus the laser to the smaller possible spot (diffraction limit). This is done by filling the entrance pupil of the objective with the laser.

Let's assume a setup in which the back lens of the objective is filled up with a laser beam (assumed to be gaussian). The waist of the beam, assuming $z_r >> f$ (where $f$ is the focal distance of the objective and $z_r$ is the depth of focus of the gaussian beam) is then $w_0 = d/2$, where $d$ is the diameter of the back lens of the objective

The experiment must be conducted in a dry atmosphere, to avoid quenching the phosporescence. Thus the sample must be in a box at vacuum or filled with an inert gas (e.g, argon, N2). The objective may be inside the box (if working distance is small) or outside (if working distance is large)

We will use a reflection objective, the MM40XF-VUV, characterized by a large working distance and large (for an air coupled) NA.
"

# ╔═╡ 4e4ab4e5-bacc-4c58-8000-10602a1f8465
#load("../notebooks/img/LMM40XVUV.png")

# ╔═╡ e35f8650-dc2e-46f3-9220-26754e8860e0
#md"### Characteristics of the MM40xVUV

#| Feature   | Value
#|:-------| ---------- |
#| Entrance pupil diameter  | $dd|
#| Focal length  | $fl|
#|NA  | $NA|
#|M (magnification)  | $MM|
#|working distance  | $wd|
#|Tranmission (250-1000 nm)  | $T|
#|Damage threshold  | $dth|
#"

# ╔═╡ 13827ddc-bebd-44d5-929f-f4ac6f43b093
#begin
#	dd    = 5.1mm
#	fl    = 5.0mm
#	NA    = 0.5
#	dth   = 0.3J/cm^2
#	wd    = 7.8mm
#	T     = 0.85
#	MM     = 40.0
#end

# ╔═╡ 2a07dd38-bac4-410c-8135-5c4e6f851df6
#lmm40xf_uvv  = lfi.LaserLab.Objective("LMM40XF-UVVV", fl, dd, MM)

# ╔═╡ ebe4e491-1c61-45cb-a559-64fa3f5fbb9d
load("../notebooks/img/NikonMUE31900.png")

# ╔═╡ d6d6dc94-5b99-4de6-a5c3-ab2be2b49d32
begin
	dd    = 2.4mm
	fl    = 2.0mm
	NA    = 0.6
	wd    = 10.0mm
	MM    = 100.0
end

# ╔═╡ a3703591-17d4-4049-8c1b-21a4c8329ffd
md"### Characteristics of the NikonMUE31900

| Feature   | Value
|:-------| ---------- |
| Entrance pupil diameter  | $dd|
| Focal length  | $fl|
|NA  | $NA|
|M (magnification)  | $MM|
|working distance  | $wd|
"

# ╔═╡ b0112287-c957-4ce3-ba24-9c47c8128b2d
tobj  = lfi.LaserLab.Objective("Nikon", fl, dd, MM)

# ╔═╡ 417fecfd-b316-4e63-8e68-303d4c568434
md"## The laser beam as a Gaussian Laser"

# ╔═╡ bbd09112-d6cd-4f0f-9a87-7b85be983e09
md"### Focusing the beam"

# ╔═╡ 8c405cc2-e8b5-4125-b21d-f7a22f94fb59
md"The beam is now focused in a narrow spot by the objective. "

# ╔═╡ ebe77978-7d7f-407d-9a73-4be3568265ef
gepl485 = lfi.LaserLab.propagate_paralell_beam(epl485, tobj)

# ╔═╡ e7e7b1ca-b27f-406c-9df0-5ffea702719f
md"### Spot size and Depth of focus

The spot size is $(round(lfi.LaserLab.spot_size(gepl485)/nm,digits=1)) nm, while the depth of focus is $(round(lfi.LaserLab.depth_of_focus(gepl485)/nm, digits=1)) nm. "

# ╔═╡ f48eb494-f9f5-4e4f-9fc0-ffa20312155a
i0mWcm2f = round(uconvert(W/cm^2, gepl485.I0)/(W*cm^-2), sigdigits=2);

# ╔═╡ fc5b2f60-2c04-4835-a222-9971189ac54e
ng0 = round(gepl485.γ0/(Hz*cm^-2), sigdigits=2);

# ╔═╡ 4eeabd93-cd82-4051-b2d0-bee3db1723e4
md"### Power density

- The power density in the spot is now much larger: $(i0mWcm2f ) W/cm2
- Or in term of photon density: $(ng0) Hz/cm2
"

# ╔═╡ e40a68c6-5885-49d1-9259-fdd3be09fee4
md"In a gaussian laser, the intensity as a function of the radial direction ($\rho$) and the direction of propagation (z) is:

$I(\rho, z) = I_0 ( W_0 / W(z))^2 \exp{-2 \rho^2/W^2(z)}$
"

# ╔═╡ f79fb47b-b4f7-442a-927f-0c186e673b80
md"Since both the spot size and the depth of focus are small, we can approximate the intensity that will illuminate the molecules of the mono layer with I(0,0)"

# ╔═╡ 8da11b18-5934-4207-941b-795969173f21
fI = lfi.LaserLab.I(gepl485) ;

# ╔═╡ ce47be17-7773-4797-beb7-202d3c02555a
begin()
zl=-5.0:0.01:5.0
p1 = plot(zl, fI.(0.0*μm, zl*μm)/(mW * cm^-2), label="I(0,z)")
xlabel!("z (μm)")
ylabel!("I(0,z)")

rl=-1.0:0.01:1.0
p2 = plot(rl, fI.(rl*μm, 0*μm)/(mW * cm^-2), label="I(ρ,0)")
xlabel!("ρ (μm)")
ylabel!("I(0,z)")

plot(p1,p2, layout = (1, 2), legend=false, fmt = :png)
end

# ╔═╡ 498d145b-ce5f-46a5-b77e-31f412e04eb9
md"### Packing of the monolayer

- One can define the packing of the ML in terms of the 'pitch' or distance separating two molecules. This is a crucial parameter of the experiment.

"

# ╔═╡ 0d4a0651-4ce4-497e-8764-2bcbbef83cf1
md"##### molecular pitch (in nm)"

# ╔═╡ d437cc55-8c19-4d0b-86d8-760da3956895
@bind mp NumberField(1.0:10.0^3; default=1.0)

# ╔═╡ 7dedd16b-4279-4f50-8120-d9ef394f3e13
pitch = mp*nm;

# ╔═╡ 649fa1b3-0dd2-487d-9bdf-7db97a0ec178
ml = lfi.LaserLab.Monolayer(pitch);

# ╔═╡ 92a5b436-1d8e-4436-8dda-8b1d3518bdea
md"This corresponds to $(uconvert(cm^-2, ml.σ)) molecules"

# ╔═╡ 8d4516fd-8838-4e5f-a61d-9dc1f65b31ad
aspot = π * gepl485.w0^2;

# ╔═╡ c5b0405a-a991-45a5-aa04-d09020b0c7f0
md"### Spot area
- The area of the spot iluminated by the beam is $(round(uconvert(μm^2,aspot)/μm^2, digits=1)) μm2"

# ╔═╡ 03e15d19-ea63-413e-8ec5-4d50e28558ac
nmol = uconvert(μm^2, aspot) * ml.σ;

# ╔═╡ 4a1cf5fc-b44f-4b19-ac8e-d610da0a17cb
md"### Number of molecules in the spot
- The number of molecules in the spot illuminated by the beam is: $(round(nmol))"

# ╔═╡ 39a7f7b4-058b-4283-b35c-4409ea9e478a
md"### Fluorescence per molecule
- The fluorescence of each molecule is the product of the beam density, the fluorescence cross section and the quantum yield:

$f = I_0 \cdot \sigma \cdot Q$
"

# ╔═╡ abe34c91-9824-45e4-8865-7b3f73ff8758
fmrs = lfi.LaserLab.fluorescence(frusl, gepl485.γ0);

# ╔═╡ 52aeae45-0dfc-48f2-a362-20d1add0ff7f
md"### Fluorescence per molecule for free and chelated species
- The fluorescence per molecule for RuSL is: $(round(fmrs/Hz)) Hz,
"

# ╔═╡ 9d97d16a-8f90-4ff7-ac0b-ee610c20ee32
sfrs = fmrs * nmol;

# ╔═╡ 960ba1ef-495d-4bb2-8aff-73cb68ae440e
md"### Total fluorescence in the spot

The total fluorescence in the spot is the product of the number of molecules and the fluorescence per molecule:

- Fluorescence in spot for RuSL = $(round(sfrs/Hz, sigdigits=1)) Hz

"

# ╔═╡ 8815a3e2-e559-4857-b20e-14aa5f91d341
md"## Filters and dichroics"

# ╔═╡ 89aa75b3-2263-450e-9f2e-cd7c875a7919
md"### Dichroic DMLP 567
- Transmits > 98% of the light above 600 nm.
"

# ╔═╡ b5045344-fb52-4b2e-88b9-1c1d4d1d50f6
begin
	dmlp567 = lfi.LaserLab.load_df_from_csv("/Users/jj/JuliaProjects/LaserLab/data/Filters", "dmlp567.csv", spG)
	dmlp567[!, :λ] = convert.(Float64, dmlp567[:, :λ])
	first(dmlp567,10)
	println("")
end

# ╔═╡ f6fefa8e-8490-4c0d-b707-c9a95011cbb1
begin
	dp1 = plot(dmlp567.λ, dmlp567.T, lw=2, label = "T", fmt = :png)
	plot(dp1, dmlp567.λ, dmlp567.R, lw=2, label = "R", fmt = :png)
	xlabel!("λ (nm)")
	ylabel!("T (R) (%)")
end

# ╔═╡ 3555a982-9879-460f-babf-6f31c11c6f7f
md"### High band pass filter FGL550"

# ╔═╡ 954d503d-19f2-4575-84d7-2d1722a28a7f
begin
	fgl550 = lfi.LaserLab.load_df_from_csv("/Users/jj/JuliaProjects/LaserLab/data/Filters", "fgl550.csv", spG)
	fgl550[!, :λ] = convert.(Float64, fgl550[:, :λ])
	first(fgl550,10)
	println("")
end

# ╔═╡ 41eb41f7-bc6e-4dce-8d39-68142fd329db
begin
	plot(fgl550.λ, fgl550.T, lw=2, label = "T", fmt = :png)
	xlabel!("λ (nm)")
	ylabel!("T  (%)")
end

# ╔═╡ bc6ba394-d1b8-4c1b-8a01-c23bcd29178c
md"### Color filters FF01"

# ╔═╡ 7140a27e-4cde-4c7d-a794-9a624e540677
begin
	ff550 = lfi.LaserLab.load_df_from_csv("/Users/jj/JuliaProjects/LaserLab/data/Filters", "FF01-550_49_Spectrum.csv", enG2)
	ff600 = lfi.LaserLab.load_df_from_csv("/Users/jj/JuliaProjects/LaserLab/data/Filters", "FF01-600_52-25.csv", enG2)
	ff650 = lfi.LaserLab.load_df_from_csv("/Users/jj/JuliaProjects/LaserLab/data/Filters", "FF01-650_54_Spectrum.csv", enG2)
	ff692 = lfi.LaserLab.load_df_from_csv("/Users/jj/JuliaProjects/LaserLab/data/Filters", "FF01-692_40_Spectrum.csv", enG2)
	ff732 = lfi.LaserLab.load_df_from_csv("/Users/jj/JuliaProjects/LaserLab/data/Filters", "FF01-732_68_Spectrum.csv", enG2)
	ff810 = lfi.LaserLab.load_df_from_csv("/Users/jj/JuliaProjects/LaserLab/data/Filters", "FF01-810_Spectrum.csv", spG)
	ff810[!, :T] = ff810[!, :T] * 0.01
	first(ff550,10)
	println("")
end

# ╔═╡ b0bf01be-6a53-4c42-9d9a-752bf4652986
begin
	ff550s = filter(df -> df.T >= 0.01, ff550)
	ff600s = filter(df -> df.T >= 0.01, ff600)
	ff650s = filter(df -> df.T >= 0.01, ff650)
	ff692s = filter(df -> df.T >= 0.01, ff692)
	ff732s = filter(df -> df.T >= 0.01, ff732)
	ff810s = filter(df -> df.T >= 0.01, ff810)
	println("")
end

# ╔═╡ 0ee9bd36-7192-4922-abaf-7d7d08605915
begin
	p550 = plot(ff550s.λ, ff550s.T, lw=2, label = "FF01-550", fmt = :png)
	p600 = plot(p550, ff600s.λ, ff600s.T, lw=2, label = "FF01-600", fmt = :png)
	p650 = plot(p600, ff650s.λ, ff650s.T, lw=2, label = "FF01-650", fmt = :png)
	p692 = plot(p650, ff692s.λ, ff692s.T, lw=2, label = "FF01-692", fmt = :png)
	p732 = plot(p692, ff732s.λ, ff732s.T, lw=2, label = "FF01-732", fmt = :png)
	p810 = plot(p732, ff810s.λ, ff810s.T, lw=2, label = "FF01-810", fmt = :png)
	xlabel!("λ (nm)")
	ylabel!("T  ")
end

# ╔═╡ df1ff4c2-8c84-455a-9cc4-7fdc34e2cb83
begin
	λm = [550.0, 600.0, 650.0, 692.0, 732.0] * nm
	Im = [0.23, 0.43, 3.44, 1.7, 1.22] * MHz
	wm = [size(ff550s)[1], size(ff600s)[1], size(ff650s)[1], size(ff692s)[1], size(ff732s)[1]] * nm
	Ilm = uconvert.(MHz*μm^-1, Im ./λm)
	Ilmx = round.(Ilm ./ (MHz*μm^-1),  sigdigits=2)
	It = sum(Im)
end

# ╔═╡ d5ef9a03-f716-454b-a5de-0c67ee935679
md"## Summary of measurements with TOPATU setup and RuSL ML
- P = $(epl485.P)
- f = $(epl485.f)

| λ filter (nm)   | FF550      | FF600  | FF650 | FF692| FF732|
|:-------| ---------- |:-----:|:-----:|:-----:|:-----:|
| λ0 (nm)  | $(λm[1])| $(λm[2]) | $(λm[3]) |$(λm[4])| $(λm[5])
| w (nm)  | $(wm[1])| $(wm[2]) | $(wm[3]) |$(wm[4])| $(wm[5])
| I (MHz)  | $(Im[1])| $(Im[2]) | $(Im[3]) |$(Im[4])| $(Im[5])
| Ilm (MHz/nm)  | $(Ilmx[1])| $(Ilmx[2]) | $(Ilmx[3]) |$(Ilmx[4])| $(Ilmx[5])

Total observed light = $It

"

# ╔═╡ 6396d4af-e668-4a73-b8be-b73b8b41267c
begin
ps2 = scatter(λm/nm, Ilmx, markersize=3,
		color = :black,
	    legend=false,
		fmt = :png)
plot(ps2, λm/nm, Ilmx, lw=2)
xlabel!("λ (nm)")
ylabel!("Il (MHz/μm)")
end

# ╔═╡ 6445e2c6-4dcc-44dd-bdda-564e4c9b3911
effCCD = lfi.LaserLab.ccd()

# ╔═╡ 7e0b0616-5dd0-44d3-beab-4fa32521d3ff
begin
	wl = 350.0:10.0:1000.0
	eccd = effCCD.(collect(wl))
	plot(wl, eccd, lw=2)
	xlabel!("λ (nm)")
	ylabel!("ϵ")
end

# ╔═╡ b49c15a7-d9de-4942-b09c-7f2ed9b4550e
function ccd_eff(lmn,wmn)
	effx = zeros(1,5)
	for i in 1:5
		effx[i] = (effCCD(lmn[i] - wmn[i]/2.0) + effCCD(lmn[i] + wmn[i]/2.0))/2.0
	end
	#@info "effx" effx
	sum(effx)/length(effx)
end

# ╔═╡ 356b7b35-6794-40d0-8c88-b8e066f086a6
begin
	lmn = λm/nm
	wmn = wm/nm
	ϵobj = 0.95
	ϵd = 0.95
	ϵf = 0.85
	ϵPMT = 0.1
	ϵNA      = lfi.LaserLab.transmission(tobj)
	efccd = ccd_eff(lmn,wmn)
	ϵT =  ϵobj^2 * ϵd^2 * ϵf * ϵPMT * ϵNA
	ϵTc =  ϵobj^2 * ϵd^2 * ϵf * efccd * ϵNA
	qf = 0.1
	println("")
end

# ╔═╡ f38318dd-f7bb-4bf1-9bc9-d1f7ed8a8397
md"## Detected Light

The detected light is the product of the emitted fluorescence and the detection efficiency, which in turns includes:

- Transmission efficiency of the objective : ϵ_obj_vuv =$(ϵobj^2)
- Tranmission due to the NA ϵ_NA = $(round(ϵNA, digits=2))
- Transmission due to the dichroic ϵd =$(ϵd^2)
- Tranmission due to the filters ϵf = $ϵf
- Transmission due to PMT ϵPMT = $ϵPMT
- Transmission due to CCD ϵCCD = $(round(efccd, digits=2))
- Quenching factor of phosphorescence due to oxygen = $qf
- The total tranmission (with PMT) is $(round(ϵT, digits=3))
- The total tranmission (with CCD) is $(round(ϵTc, digits=3))
"

# ╔═╡ f615fc39-bfd9-45ce-84fe-a28921bde525

# ╔═╡ 43e8a5f1-512c-46b7-91f6-89d3c7e81368
begin
	osf  = sfrs * ϵT * qf
	osfc = sfrs * ϵTc * qf
end

# ╔═╡ d49e7e9b-6487-407b-aef7-2884461879a0
md" ### Expected observed light in the PMT

Thus the total expected light in the PMT for a ML of ~ $(uconvert(cm^-2, ml.σ)) molecules is:

- Expected observed light:
  -  with PMT= $(round(osf/Hz, sigdigits=1)) Hz
  -  with CCD= $(round(osfc/Hz, sigdigits=1)) Hz

"

# ╔═╡ ae5f17d2-3c31-4a9c-85e4-f2f22625d86b
begin
	P0 = 1.5μW
	f0  = 500kHz
	println("")
end

# ╔═╡ c3952378-79bb-4605-8ffb-828c1e3e3321
#load("../notebooks/img/RuSlCCD250222.png")

# ╔═╡ 0f7ff8a4-43e4-4829-a4b4-78594664cee2
md"## TCSP experiments"

# ╔═╡ 131f4e45-85f5-43bc-8c32-d59f8803bfd6
md"### Typical setup"

# ╔═╡ a4066207-4c50-4360-b43f-218d5355ff3e
load("../notebooks/img/tcsp_setup.png")

# ╔═╡ 218fc76e-1060-40ff-a52f-d884441380e2
md"### Photon counting technique"

# ╔═╡ 58b61eac-e1de-48f3-a37b-ecc1f8d5f8ab
load("../notebooks/img/tcspc.png")

# ╔═╡ 320f9c75-5047-4a4c-890b-ffdc70767634
md" The phothon counting technique requires that one photon is recorded on average per pulse (in this case in the period 1-5 mus). This is to avoid inefficiencies which can bias the measurement as ilustrated below."

# ╔═╡ b97cd0a9-ee46-4520-a408-f60150bbea74
load("../notebooks/img/tcspc_deadtime.png")

# ╔═╡ b96961fd-6a9a-40bb-b582-2b3586f56edb
md"However, most TCSPC techniques are designed for short interval times. The dead time is typically a few nanoseconds, thus no loss for dead time is expected here. Nevertheless the system is designed to record just one photon per pulse. This implies that the average number of photons must be reduced to 1 per 1 (5) μs. This can be done by:
- attenuating the laser light
- reducing the laser power
- spacing the molecules in the ML
"

# ╔═╡ eba3554f-ae6f-4a3d-8131-43ffa2743977
md"# Appendix"

# ╔═╡ f57ca807-9b3a-4f77-9607-74ee3a411990
md"## Fitting"

# ╔═╡ 7a6a3b8a-8b1c-41d5-ab25-7d294f1bee3d
md"### *func1dfit* is a light wrapper to curve_fit"

# ╔═╡ 92698c05-edf3-4b9e-a10a-1da9ed0dc82a
"""
    func1dfit(ffit::Function, x::Vector{<:Real},
              y::Vector{<:Real}, p0::Vector{<:Real},
              lb::Vector{Float64}, ub::Vector{Float64})

Fit a function to the data x, y with start prediction p0
and return coefficients and errors.
"""
func1dfit(ffit::Function, x::Vector{<:Real},
          y::Vector{<:Real}, p0::Vector{<:Real},
          lb::Vector{<:Real}, ub::Vector{<:Real}) = curve_fit(ffit, x, y, p0,
			                                                 lower=lb, upper=ub)

# ╔═╡ 0dd448d1-9952-438c-9284-e0c320c955aa
md"### Example: fit to a polynomial"

# ╔═╡ f8b61a8e-2679-48cd-842c-17d6e6ee760e

pol3(x, a, b, c, d) = a + b*x + c*x^2 + d*x^3

# ╔═╡ bfc6e346-3525-47d7-a436-ffa02dab11b9
begin
	err_sigma = 0.04
	x=collect(LinRange(0., 10., 100))
	p0 = [10.0, 1.0, 0.7, 0.5]
	y = pol3.(x, p0...)
    y += rand(Normal(0, err_sigma), length(y))
    lb = fill( 0.0, length(p0))
    ub = fill(20.0, length(p0))
    pol3_fit = @. pol(x, p) = p[1] + p[2] * x + p[3] * x^2 + p[4] * x^3
    fq = func1dfit(pol3_fit, x, y, p0, lb, ub)
	cfq = coef(fq)
	sfq = stderror(fq)
end

# ╔═╡ 2ff6f551-2662-43a5-9776-70951ff40364
@info "fit coefficients" cfq

# ╔═╡ ec3cd40c-6a52-48fd-a3ce-9e091250e981
yf = pol3.(x, cfq...);

# ╔═╡ 525b0316-8374-4882-aa5b-84f2bf33c5c3
@info "coefficient errors (std)" sfq

# ╔═╡ 840fa50e-bf0f-42f9-8e7f-ec3cc4bdb9af
sfq

# ╔═╡ 5738e27c-5307-4623-bd43-67f66b7b97d2
@info "margin_of_error (90%)" margin_error(fq, 0.1)

# ╔═╡ 14360d77-340e-488d-bed0-00f408ef1dd4
all(isapprox.(cfq, p0; atol=err_sigma))

# ╔═╡ 685fd840-64b1-427d-83f3-217664ea9798
begin
pp1 = scatter(x, y,
	          label="p3",
			  markersize=2,
			  color = :black,
	          legend=false,
		   	  fmt = :png)
pp2 = plot(pp1, x, yf, lw=2)

xlabel!("x")
ylabel!("p3(x)")
end

# ╔═╡ 8ecc8206-275b-4d55-9b01-e82e7f2df5dc
md"### Fit an exponential"

# ╔═╡ aac7a640-97fe-46c0-89b3-7e76978baf0d
exp(1.0)

# ╔═╡ 89ef5f08-7a1e-42a5-8fcf-b7efa63a8a68
expo(t, N, λ) = N*exp(-t/λ)

# ╔═╡ 426f8190-7865-4bd1-bc23-0adb5fc1892c
begin
	err_sigma2 = 0.05
	t=collect(LinRange(0.0, 5000.0, 1000))
	p0t = [3500.0, 500.0]
	yt = expo.(t, p0t...)
	yts = yt + rand(Normal(0, err_sigma2), length(yt)) .* yt
    lbt = [0.0, 0.0]
    ubt = [50000.0, 50000.0]
    expo_fit = @. expo(t, p) = p[1]*exp(-t/p[2])
    fqe = func1dfit(expo_fit, t, yts, p0t, lbt, ubt)
	cfqe = coef(fqe)
	sfqe = stderror(fqe)
end

# ╔═╡ 2fdad658-0b71-4d42-8591-9284fee3aeb6

begin
	tft = expo.(tdata, coef(fit)...);
	println("")
end

# ╔═╡ 5a0cc7d8-a3c3-45b9-8b91-9528e4938fb0
begin
ps1 = scatter(tdata, ydata, markersize=1,
		color = :black,
	    legend=false,
		fmt = :png)
pp = plot(ps1, tdata, tft, lw=2,fmt = :png)
xlabel!("t (ns)")
ylabel!("I (a.u.)")
end

# ╔═╡ f1e6c6b6-49fa-44bb-bb9a-07986f6a1d13
expo(0,3000.,500.)

# ╔═╡ c54205d5-ea1a-4416-b4b0-3bdb168dae61
fqe.converged

# ╔═╡ b83d4ff4-cc73-4cf7-abad-78da291eb404
@info "fit coefficients" cfqe

# ╔═╡ 13ad7221-0d2d-4a77-9258-9edace85fde0
@info "coefficient errors (std)" sfqe

# ╔═╡ f2b2cc4d-54ed-4f2f-80bc-bc3bf82bb2e8
begin
pp3 = scatter(t, yt,
	          label="expo",
			  markersize=2,
			  color = :black,
	          legend=false,
		   	  fmt = :png)
#pp2 = plot(pp1, x, yf, lw=2)

xlabel!("t")
ylabel!("expo(t)")
end

# ╔═╡ Cell order:
# ╠═0b4c1052-9180-11ec-2cdd-1b61bde8bdc6
# ╠═2d15cff8-622e-4f60-a16f-77a4d16ccc0d
# ╠═ffc498a9-3522-4309-9ced-48955cab1759
# ╠═91aa2620-9728-4de7-9c71-f01a6c765350
# ╠═9d66c472-b742-49b3-9495-1c3be998108b
# ╠═c3b89fa7-67f8-4cc3-b3ba-fcb304f44669
# ╠═bd9eed27-fc29-41ff-bf12-ec0cd3881ea2
# ╠═ca6f7233-a1d1-403d-9416-835a51d360db
# ╟─4c6c9627-d80d-41ac-8141-f935df893b73
# ╟─34b3d66e-3b39-43dc-922b-9ba6c15c1bac
# ╟─7e150a69-0070-4371-b5c0-02b7ad70d813
# ╠═16cdbadf-4deb-4a66-8ab7-84437a4fe3d4
# ╟─4a0e30fe-398e-4d80-86a2-bfc384cbc6e6
# ╠═9cc3dfcb-b95d-4086-a228-ed4753f6ca0d
# ╟─591f9fcb-7ad1-400e-af0a-29686a4914ec
# ╠═768feb2a-855f-4db6-9ae0-6abec23707c1
# ╠═9cc092e9-40d2-4c5a-9b55-4e2adccb3382
# ╠═d75a8426-c97b-422e-8c69-3f3ce94c5370
# ╟─a7d330cf-5093-490f-adc8-e482e3084806
# ╟─280a17c4-6d5e-4ff3-a5ab-dbbcb7199bb2
# ╟─a230e061-b6ab-49d4-bdad-730d75e20e9c
# ╠═6d5e295c-1c57-4a92-aee2-3ffcbb3306df
# ╠═80447dc7-1103-4019-bc64-283d4a9368a0
# ╠═08ec5b69-946f-407e-ac8d-bf8814cc0121
# ╠═2fdad658-0b71-4d42-8591-9284fee3aeb6
# ╠═5a0cc7d8-a3c3-45b9-8b91-9528e4938fb0
# ╟─58b227ac-627f-49ba-bc58-75039c65733b
# ╠═4223f9f8-8ee6-4ea9-a7bc-6379c81c48c0
# ╟─fabba61a-f28d-4535-b0e3-e94f5098bb0a
# ╟─7a1fcff4-6948-4333-b971-8fbf31c7d3ee
# ╠═a917927d-d471-4465-abd1-30d959af0b45
# ╠═93b78d1a-ee07-4d66-b5a4-c7a109a24f81
# ╟─aa4283d7-37d4-4aaf-a5b7-b515466e545d
# ╟─19b7d7fb-adab-4f29-82af-d059525921ee
# ╠═102a2054-2fd5-406e-8bb4-20ecb47c278f
# ╟─e97474ef-145f-45c8-96f6-213acf4a3b41
# ╟─106f7063-719e-4eff-857d-3566d6e6d4c8
# ╟─177762d5-b9f2-449b-abba-256f3d4318ad
# ╟─4e4ab4e5-bacc-4c58-8000-10602a1f8465
# ╟─e35f8650-dc2e-46f3-9220-26754e8860e0
# ╟─13827ddc-bebd-44d5-929f-f4ac6f43b093
# ╟─2a07dd38-bac4-410c-8135-5c4e6f851df6
# ╟─a3703591-17d4-4049-8c1b-21a4c8329ffd
# ╟─ebe4e491-1c61-45cb-a559-64fa3f5fbb9d
# ╟─d6d6dc94-5b99-4de6-a5c3-ab2be2b49d32
# ╠═b0112287-c957-4ce3-ba24-9c47c8128b2d
# ╟─417fecfd-b316-4e63-8e68-303d4c568434
# ╟─bbd09112-d6cd-4f0f-9a87-7b85be983e09
# ╟─8c405cc2-e8b5-4125-b21d-f7a22f94fb59
# ╟─ebe77978-7d7f-407d-9a73-4be3568265ef
# ╟─e7e7b1ca-b27f-406c-9df0-5ffea702719f
# ╟─f48eb494-f9f5-4e4f-9fc0-ffa20312155a
# ╟─4eeabd93-cd82-4051-b2d0-bee3db1723e4
# ╟─fc5b2f60-2c04-4835-a222-9971189ac54e
# ╟─e40a68c6-5885-49d1-9259-fdd3be09fee4
# ╟─f79fb47b-b4f7-442a-927f-0c186e673b80
# ╠═8da11b18-5934-4207-941b-795969173f21
# ╟─ce47be17-7773-4797-beb7-202d3c02555a
# ╟─498d145b-ce5f-46a5-b77e-31f412e04eb9
# ╟─0d4a0651-4ce4-497e-8764-2bcbbef83cf1
# ╟─d437cc55-8c19-4d0b-86d8-760da3956895
# ╟─7dedd16b-4279-4f50-8120-d9ef394f3e13
# ╟─649fa1b3-0dd2-487d-9bdf-7db97a0ec178
# ╟─92a5b436-1d8e-4436-8dda-8b1d3518bdea
# ╟─c5b0405a-a991-45a5-aa04-d09020b0c7f0
# ╟─8d4516fd-8838-4e5f-a61d-9dc1f65b31ad
# ╟─4a1cf5fc-b44f-4b19-ac8e-d610da0a17cb
# ╟─03e15d19-ea63-413e-8ec5-4d50e28558ac
# ╟─39a7f7b4-058b-4283-b35c-4409ea9e478a
# ╟─52aeae45-0dfc-48f2-a362-20d1add0ff7f
# ╟─abe34c91-9824-45e4-8865-7b3f73ff8758
# ╟─960ba1ef-495d-4bb2-8aff-73cb68ae440e
# ╟─9d97d16a-8f90-4ff7-ac0b-ee610c20ee32
# ╟─8815a3e2-e559-4857-b20e-14aa5f91d341
# ╟─89aa75b3-2263-450e-9f2e-cd7c875a7919
# ╟─b5045344-fb52-4b2e-88b9-1c1d4d1d50f6
# ╟─f6fefa8e-8490-4c0d-b707-c9a95011cbb1
# ╟─3555a982-9879-460f-babf-6f31c11c6f7f
# ╟─954d503d-19f2-4575-84d7-2d1722a28a7f
# ╟─41eb41f7-bc6e-4dce-8d39-68142fd329db
# ╟─bc6ba394-d1b8-4c1b-8a01-c23bcd29178c
# ╟─7140a27e-4cde-4c7d-a794-9a624e540677
# ╠═b0bf01be-6a53-4c42-9d9a-752bf4652986
# ╟─0ee9bd36-7192-4922-abaf-7d7d08605915
# ╠═d5ef9a03-f716-454b-a5de-0c67ee935679
# ╠═df1ff4c2-8c84-455a-9cc4-7fdc34e2cb83
# ╠═6396d4af-e668-4a73-b8be-b73b8b41267c
# ╠═f38318dd-f7bb-4bf1-9bc9-d1f7ed8a8397
# ╠═356b7b35-6794-40d0-8c88-b8e066f086a6
# ╠═6445e2c6-4dcc-44dd-bdda-564e4c9b3911
# ╠═7e0b0616-5dd0-44d3-beab-4fa32521d3ff
# ╠═b49c15a7-d9de-4942-b09c-7f2ed9b4550e
# ╠═f615fc39-bfd9-45ce-84fe-a28921bde525
# ╠═d49e7e9b-6487-407b-aef7-2884461879a0
# ╠═43e8a5f1-512c-46b7-91f6-89d3c7e81368
# ╟─ae5f17d2-3c31-4a9c-85e4-f2f22625d86b
# ╠═c3952378-79bb-4605-8ffb-828c1e3e3321
# ╟─0f7ff8a4-43e4-4829-a4b4-78594664cee2
# ╟─131f4e45-85f5-43bc-8c32-d59f8803bfd6
# ╠═a4066207-4c50-4360-b43f-218d5355ff3e
# ╟─218fc76e-1060-40ff-a52f-d884441380e2
# ╟─58b61eac-e1de-48f3-a37b-ecc1f8d5f8ab
# ╠═320f9c75-5047-4a4c-890b-ffdc70767634
# ╠═b97cd0a9-ee46-4520-a408-f60150bbea74
# ╠═b96961fd-6a9a-40bb-b582-2b3586f56edb
# ╠═eba3554f-ae6f-4a3d-8131-43ffa2743977
# ╠═f57ca807-9b3a-4f77-9607-74ee3a411990
# ╠═7a6a3b8a-8b1c-41d5-ab25-7d294f1bee3d
# ╠═92698c05-edf3-4b9e-a10a-1da9ed0dc82a
# ╠═0dd448d1-9952-438c-9284-e0c320c955aa
# ╠═f8b61a8e-2679-48cd-842c-17d6e6ee760e
# ╠═bfc6e346-3525-47d7-a436-ffa02dab11b9
# ╠═2ff6f551-2662-43a5-9776-70951ff40364
# ╠═ec3cd40c-6a52-48fd-a3ce-9e091250e981
# ╠═525b0316-8374-4882-aa5b-84f2bf33c5c3
# ╠═840fa50e-bf0f-42f9-8e7f-ec3cc4bdb9af
# ╠═5738e27c-5307-4623-bd43-67f66b7b97d2
# ╠═14360d77-340e-488d-bed0-00f408ef1dd4
# ╠═685fd840-64b1-427d-83f3-217664ea9798
# ╠═8ecc8206-275b-4d55-9b01-e82e7f2df5dc
# ╠═aac7a640-97fe-46c0-89b3-7e76978baf0d
# ╠═89ef5f08-7a1e-42a5-8fcf-b7efa63a8a68
# ╠═f1e6c6b6-49fa-44bb-bb9a-07986f6a1d13
# ╠═426f8190-7865-4bd1-bc23-0adb5fc1892c
# ╠═c54205d5-ea1a-4416-b4b0-3bdb168dae61
# ╠═b83d4ff4-cc73-4cf7-abad-78da291eb404
# ╠═13ad7221-0d2d-4a77-9258-9edace85fde0
# ╠═f2b2cc4d-54ed-4f2f-80bc-bc3bf82bb2e8

ext: jl

max_stars_repo_name: proscrite/LaserLab

max_issues_repo_name: proscrite/LaserLab

max_forks_repo_name: proscrite/LaserLab


text:
@testset "Test Stopping user structure" begin

  ###############################################################################
  #
  # The Stopping structure eases the implementation of algorithms and the
  # stopping criterion.
  # We illustrate here the basic features of Stopping.
  #
  # Stopping has a mysterious attribute "stopping_user_struct", we illustrate
  # in this example how convenient it can be when specializing your own object.
  # stop! also calls a function _user_check! that can also be specialized by the
  # user.
  # These user-dependent attribute/function allow specializing a new Stopping
  # without too much of copy/paste.
  #
  ###############################################################################
  #using LinearAlgebra, NLPModels, Stopping, Test

  #Let us consider an NLPStopping:
  f, x0 = x -> x[1], ones(6)
  meta = NLPModelMeta(
    6,
    x0 = x0,
    lvar = -Inf * ones(6),
    uvar = Inf * ones(6),
    ncon = 1,
    y0 = [0.0],
    lcon = [-Inf],
    ucon = [6.0],
  )
  nlp = ADNLPModel(meta, Counters(), ADNLPModels.ForwardDiffAD(1, 1, f, x0), f, x -> [sum(x)])

  structtest = Dict(:feasible => true, :rho => 1e-1)
  stop_bnd = NLPStopping(nlp, user_struct = structtest)

  io = IOBuffer()
  show(io, stop_bnd)

  @test stop_bnd.stopping_user_struct[:feasible] == true
  @test stop_bnd.stopping_user_struct[:rho] == 1e-1

  ##############################################################################
  function uss_func(stp::NLPStopping, start::Bool)
    x = stp.current_state.x
    cx = stp.current_state.cx
    feas =
      max.(stp.pb.meta.lcon - cx, cx - stp.pb.meta.ucon, stp.pb.meta.lvar - x, x - stp.pb.meta.uvar)
    tol = max(stp.meta.atol, stp.stopping_user_struct[:rho])
    stp.stopping_user_struct[:feasible] = norm(feas, Inf) <= tol
  end

  stop_bnd.meta.user_check_func! = uss_func

  #=
  import Stopping._user_check!
  #We now redefine _user_check! to verify the feasibility
  function _user_check!(stp :: NLPStopping, x :: T) where T
   cx   = stp.current_state.cx
   feas = max.( stp.pb.meta.lcon - cx,
                cx - stp.pb.meta.ucon,
                stp.pb.meta.lvar - x,
                x - stp.pb.meta.uvar )
   tol = max(stp.meta.atol, stp.stopping_user_struct.rho)
   stp.stopping_user_struct.feasible = norm(feas, Inf) <= tol
  end

  #remove the _user_check!(stp :: NLPStopping, x :: Stopping.Iterate) from the workspace
  #Base.delete_method(which(_user_check!, (NLPStopping,Union{Number, AbstractVector},)))
  =#
  ##############################################################################

  sol = 2 * ones(6)
  fill_in!(stop_bnd, sol)
  stop!(stop_bnd)
  @test stop_bnd.stopping_user_struct[:feasible] == false

  sol2 = ones(6)
  sol2[1] = 1.0 + 1e-1
  reinit!(stop_bnd, rstate = true)
  fill_in!(stop_bnd, sol2)
  stop!(stop_bnd)
  #since violation of the constraints smaller than rho.
  @test stop_bnd.stopping_user_struct[:feasible] == true
end


text:
using IGC
using Test
# import IGC: IGCDocument, Abstract_IGC_record
import IGC: IGCTimeFormat, IGCDateFormat, EmptyTime, EmptyDate
import IGC: IGCLatitude, IGCLongitude, FixValidity, IGCPressureAltitude, IGCGpsAltitude, TZ
import IGC: flc, A_record, B_record, D_record, GpsQualifier
import IGC: Abstract_C_record, C_record_task_info, C_record_waypoint_info
import IGC: E_record, F_record, G_record, IGCExtension, I_record, J_record
import IGC: Abstract_H_record, update, DEFAULT_SOURCE
import IGC: H_record_FiXAccuracy, H_record_DaTE, H_record_PiLoT, H_record_Copilot, H_record_GliderType, H_record_GliderRegistration
import IGC: H_record_GpsDatum, H_record_FirmwareRevision, H_record_HardwareRevision, H_record_ManufacturerModel
import IGC: H_record_GpsReceiver, MaxAlt, H_record_PressureAltitudeSensor
import IGC: H_record_CompetitionId, H_record_CompetitionClass
import IGC: H_record_TimeZoneOffset, H_record_MeansOfPropulsion, H_record_Site, H_record_UnitsOfMeasure
import IGC: K_record, L_record
import IGC: update!, Header, EOL_DEFAULT

using Dates
using TimeZones

include("test_reader.jl")
include("test_writer.jl")
include("test_filename.jl")


text:
struct PlotContext{T₁ <: AbstractControl,   T₂ <: AbstractModel,  T₃ <: AbstractDisplayProperties} <: AbstractPlotContext
    control::T₁
    model::T₂
    display_properties::T₃
end

function (context::PlotContext{<: PlotlinesControl,   <: PlotlinesModel,  <: PlotlinesDisplayProperties})()
        control = context.control
        model = context.model
        display_properties = context.display_properties
        id = get_id(display_properties)
        captioner = get_captioner(display_properties)
        caption = captioner(nothing)
        is_new_window(display_properties) ? CImGui.Begin(caption*id,C_NULL, CImGui.ImGuiWindowFlags_NoBringToFrontOnFocus) : nothing
        isenabled(control) ? control(model, display_properties) : nothing
        is_new_window(display_properties) ? CImGui.End() : nothing
end


text:
a  = read_SSR(joinpath(dirname(@__FILE__), "..", "data", "test_Hz19.5-testing.bdf"))

@test a.processing["Side"] == "Bilateral"
@test a.processing["Name"] == "P2"
@test a.processing["Amplitude"] == 70.0
@test a.processing["Carrier_Frequency"] == 1000.0
@test modulationrate(a) == 40.0390625


text:
"""
    lcp(str1::T, str2::T)::T where T <: AbstractString

Get the longest common prefix.
"""
function lcp(str1::T, str2::T)::T where T <: AbstractString
    minl::Int64 = min(length(str1), length(str2))
    outs::T = ""
    minl == 0 && return outs
    for i in 1:minl
      str1[i] == str2[i] ? outs *= str1[i] : return outs
    end # for
    return outs
end


text:
@Base.deprecate parseWKT(geom::String) readgeom(geom)
@Base.deprecate geomFromWKT(geom::String) _readgeom(geom)
@Base.deprecate geomToWKT(geom::Ptr{GEOSGeometry}) writegeom(geom)


text:
module JuliaCommunityStatistics

using GitHub
using ProgressMeter
using Dates
using DataFrames
import GitHub: name

export jlrepo, auth
const auth = authenticate(ENV["GH_AUTH"])
const jlrepo = repo("JuliaLang/julia"; auth=auth)

export get_all_prs
function get_all_prs(;state="all")
    prs = PullRequest[]
    @showprogress 1 "Fetching... " for iP = 1:1000
        thisprs, _ = pull_requests(jlrepo; auth=auth, page_limit=1, params = Dict("state"=>state, "page"=>iP))
        !isempty(thisprs) || break;
        prs = append!(prs, thisprs)
    end
    prs
end

function get_pr_info_df(pr)
    prdf = pr_to_df(pr)
    cs, _ = commits(jlrepo, pr;auth=auth)
    cfs = pull_request_files(jlrepo, pr;auth=auth)
    cdfs = commit_to_df.(cs)
    cfdfs = changed_file_to_df.(cfs)
    insertcols!(prdf, :commits=>[cdfs], :changed_files=>[cfdfs])
end

export get_all_pr_info_df
function get_all_pr_info_df(prs)
    df = DataFrame()
    @showprogress 1 "Fetching... " for pr in prs
        try
            prdf = get_pr_info_df(pr)
            df = vcat(df, prdf)
        catch
            sleep_till_reset()
            prdf = get_pr_info_df(pr)
            df = vcat(df, prdf)
        end
    end
    df
end

export get_commits
function get_commits(prs)
    prd = Dict()
    @showprogress 1 "Fetching... " for iPR in prs
        try
            prd[iPR], _ = commits(jlrepo, pr;auth=auth)
        catch
            sleep_till_reset()
            prd[iPR], _ = commits(jlrepo, pr;auth=auth)
        end
    end
    prd
end

export get_changed_files
function get_changed_files(prs)
    prd = Dict()
    @showprogress 1 "Fetching... " for iPR in prs
        try
            prd[iPR] = pull_request_files(jlrepo, iPR;auth=auth)
        catch
            sleep_till_reset()
            prd[iPR] = pull_request_files(jlrepo, iPR;auth=auth)
        end
    end
    prd
end

export sleep_till_reset
function sleep_till_reset()
    rate_lim = rate_limit(auth=auth)["rate"]
    if rate_lim["remaining"] > 0
        return
    end
    reset_time = Dates.unix2datetime(rate_lim["reset"])
    sleeptime = Dates.Millisecond(reset_time - now()) + Dates.Millisecond(2000)
    println("Sleeping for $(sleeptime)")
    sleep(sleeptime)
end

export pr_to_df
function pr_to_df(pr)
    DataFrame(
        number = pr.number,
        state = pr.state,
        owner = name(pr.user),
        created_at = pr.created_at,
        closed_at = pr.closed_at,
        updated_at = pr.updated_at,
        merged_at = pr.merged_at,
        base = name(pr.base),
        head = name(pr.head),
        merge_commit = pr.merge_commit_sha
    )
end

export commit_to_df
function commit_to_df(c)
    DataFrame(
        sha = c.sha,
        author = name(c.author)
    )
end

export changed_file_to_df
function changed_file_to_df(cf)
    DataFrame(
        file = cf.filename,
        status = cf.status,
        changes = cf.changes,
        additions = cf.additions,
        deletions = cf.deletions
    )
end
end


text:
module Runtests

include("TestFunctions.jl")

end


text:
module GLMakie

using ModernGL, FixedPointNumbers, Colors, GeometryBasics
using Makie, FileIO

using Makie: @key_str, Key, broadcast_foreach, to_ndim, NativeFont
using Makie: Scene, Lines, Text, Image, Heatmap, Scatter
using Makie: convert_attribute, @extractvalue, LineSegments
using Makie: @get_attribute, to_value, to_colormap, extrema_nan
using Makie: ClosedInterval, (..)
using Makie: inline!
using Makie: spaces, is_data_space, is_pixel_space, is_relative_space, is_clip_space
using ShaderAbstractions
using FreeTypeAbstraction
using GeometryBasics: StaticVector

using Base: RefValue
import Base: push!, isopen, show
using Base.Iterators: repeated, drop

using LinearAlgebra

# re-export Makie, including deprecated names
for name in names(Makie, all=true)
    if Base.isexported(Makie, name)
        @eval using Makie: $(name)
        @eval export $(name)
    end
end

export inline!
import ShaderAbstractions: Sampler, Buffer
export Sampler, Buffer

struct GLBackend <: Makie.AbstractBackend
end

loadshader(name) = normpath(joinpath(@__DIR__, "..", "assets", "shader", name))

# don't put this into try catch, to not mess with normal errors
include("gl_backend.jl")

function activate!(use_display=true)
    b = GLBackend()
    Makie.register_backend!(b)
    Makie.set_glyph_resolution!(Makie.High)
    Makie.current_backend[] = b
    Makie.inline!(!use_display)
end

function __init__()
    activate!()
end

export set_window_config!

if Base.VERSION >= v"1.4.2"
    include("precompile.jl")
    _precompile_()
end

end


text:
"""
    MultivariateUniform(lower, upper)

Uniform distribution in `d` dimensions, where `d == length(lower) == length(upper)`.
# Example
```julia-repl
julia> d = MultivariateUniform([-2., -3], [4, 1])
MultivariateUniform{Float64} in 2 dimensions

julia> logpdf(d, [-3, 0])
-Inf

julia> logpdf(d, [-1, 0])
-3.1780538303479453
```
"""
struct MultivariateUniform{S} <: ContinuousMultivariateDistribution
    lower::Vector{S}
    upper::Vector{S}
    range::Vector{S}
    function MultivariateUniform(lower, upper)
        length(lower) != length(upper) && error("Lower and upper bounds must have the same length.")
        (|)((lower .> upper)...) && error("Upper bounds must not be smaller than lower bounds")
        new{eltype(lower)}(lower, upper, upper .- lower)
    end
end
Base.length(d::MultivariateUniform) = length(d.lower)
lower(d::MultivariateUniform) = d.lower
upper(d::MultivariateUniform) = d.upper
function Base.show(io::IO, ::MIME"text/plain", d::MultivariateUniform{S}) where S
    println(io, "MultivariateUniform{$S} in $(length(d.lower)) dimensions")
end
function Distributions._rand!(rng::Random.AbstractRNG,
                              d::MultivariateUniform,
                              x::AbstractVector{<:Real})
#     length(x) != length(d) && error("length(x) = $(length(x)), but length(distribution) = $(length(d))")
    for i in eachindex(d.lower)
        x[i] = d.lower[i] + rand(rng) * d.range[i]
    end
    x
end
function Distributions._logpdf(d::MultivariateUniform, x::AbstractVector{<:Real})
    for i in eachindex(x)
        (d.lower[i] > x[i] || d.upper[i] < x[i]) && return -Inf
    end
    -sum(log.(d.range))
end
Distributions.logpdf(d::Nothing, ::Any) = 0.

"""
    TruncatedMultivariateNormal(mvnormal, lower, upper)

Multivariate normal distribution `mvnormal` truncated to a box
given by `lower` and `upper` bounds.
Simple rejection sampling is implemented.
IMPORTANT: `logpdf` is not properly normalized.
# Example
```julia-repl
julia> using Distributions

julia> d = TruncatedMultivariateNormal(MultivariateNormal([2., 3.], .3*I), [0, 0], [4, 5]);

julia> logpdf(d, [2, 3])
-0.6339042620834094

julia> logpdf(d, [2, 3]) == logpdf(d.mvnormal, [2, 3])
true

julia> logpdf(d, [-1, 1])
-Inf
```
"""
struct TruncatedMultivariateNormal{B, T} <: ContinuousMultivariateDistribution
    mvnormal::B
    lower::Vector{T}
    upper::Vector{T}
end
"""
    TruncatedMultivariateNormal(m, cov; lower, upper)

Constructs a `TruncatedMultivariateNormal` with mean `m` covariance matrix `cov`
and bounds `lower` and `upper`.
# Example
```julia-repl
julia> d = TruncatedMultivariateNormal([0, 0], 3I;
                                       lower = [-1, -4], upper = [4, 5])
```
"""
function TruncatedMultivariateNormal(m, sig; lower, upper)
    TruncatedMultivariateNormal(MultivariateNormal(m, sig), lower, upper)
end
Base.length(d::TruncatedMultivariateNormal) = length(d.lower)
lower(d::TruncatedMultivariateNormal) = d.lower
upper(d::TruncatedMultivariateNormal) = d.upper
function Distributions._rand!(rng::Random.AbstractRNG,
                              d::TruncatedMultivariateNormal,
                              x::AbstractVector{<:Real})
    # To do this more efficiently
    # see e.g. https://github.com/BrianNaughton/TruncatedMVN.jl/blob/master/src/multivariate.jl
    while true
        Distributions._rand!(rng, d.mvnormal, x)
        for i in eachindex(x)
            (d.lower[i] > x[i] || d.upper[i] < x[i]) && break
            if i == length(x)
                return x
            end
        end
    end
end
function Distributions._logpdf(d::TruncatedMultivariateNormal,
                               x::AbstractVector{<:Real})
    # This is wrong; but it doesn't matter here.
    # A constant is missing due to truncation.
    for i in eachindex(x)
        (d.lower[i] > x[i] || d.upper[i] < x[i]) && return -Inf
    end
    Distributions._logpdf(d.mvnormal, x)
end

"""
    MultiParticleNormal(particles, weights;
                        lower = nothing, upper = nothing,
                        diagonal = false,
                        scaling = 2, regularization = 0.)

Mixture of `k` multivariate Gaussians with means given by an array of `particles`
(`length(particles) == k`), weighted by `weights` and with covariance matrix
given by `scaling * cov(particles) + regularization * I` for each Gaussian,
if `diagonal == false`; otherwise `scaling * std(particles) + regularization`.
Sampling truncated Gaussians with diagonal covariance matrix can be more efficient.
The mean of the distribution is the weighted average of the particles.

If `lower` and `upper` are not `nothing` the Gaussians are truncated to the box
confined by the `lower` and `upper` bounds.
Note that the logpdf of the truncated version is not properly normalized
(see also [`TruncatedMultivariateNormal`](@ref)).
# Example
```julia-repl
julia> d = MultiParticleNormal([[0, 0], [2, 3], [4, 0]], [1/3, 1/3, 1/3]);

julia> mean(d)
2-element Array{Float64,1}:
 2.0
 1.0
```
"""
mutable struct MultiParticleNormal{P, S, B} <: ContinuousMultivariateDistribution
    particles::Vector{P}
    weights::Vector{Float64}
    cov::S
    scaling::Float64
    regularization::Float64
    lower::B
    upper::B
end
"""
    MultiParticleNormal(; kwargs...)

Returns a anonymous function that accepts `particles` as input and creates a
`MultiParticleNormal` distributions with uniform weights. It accepts the same
keyword arguments `kwargs...` as the normal `MultiParticleNormal` constructor.
# Example
```julia-repl
julia> constructor = MultiParticleNormal(scaling = 3, lower = [-1, 0], upper = [5, 5]);

julia> d = constructor([[0, 0], [2, 3], [4, 0], [5, 1]])
MultiParticleNormal(4 particles in 2 dimensions, scaling = 3.0, regularization = 0.0, with bounds)

julia> d.weights
4-element Array{Float64,1}:
 0.25
 0.25
 0.25
 0.25
```
"""
function MultiParticleNormal(; kwargs...)
    particles -> MultiParticleNormal(particles,
                                     fill(1/length(particles), length(particles));
                                     kwargs...)
end
function cov_matrix(::Type{<:PDMats.AbstractPDMat}, particles, weights, scaling, regularization)
    m = scaling .* cov(particles[weights .> 0]) + regularization * I
    for _ in 1:10
        try
            return PDMats.PDMat(m)
        catch err
            if typeof(err)!=LinearAlgebra.PosDefException
                throw(err)
            end
            n = size(m, 1)
            ϵ = 1e-6 * tr(m)/n
            @inbounds for i in 1:n
                m[i, i] += ϵ
            end
        end
    end
end
function cov_matrix(::Type{<:AbstractVector}, particles, weights, scaling, regularization)
    scaling .* std(particles[weights .> 0]) .+ regularization
end
function MultiParticleNormal(particles, weights;
                             lower = nothing, upper = nothing,
                             diagonal = false,
                             scaling = 2, regularization = 0.)
    sum(weights) ≈ 1 || error("Particle weights must be normalized.")
    MultiParticleNormal(particles, weights,
                        cov_matrix(diagonal ? Vector : PDMats.PDMat,
                                   particles, weights, scaling, regularization),
                        float(scaling), float(regularization),
                        lower, upper)
end
Base.length(d::MultiParticleNormal) = size(d.cov, 1)
lower(d::MultiParticleNormal{<:Any, <:Any, <:AbstractVector}) = d.lower
upper(d::MultiParticleNormal{<:Any, <:Any, <:AbstractVector}) = d.upper
function Base.show(io::IO, ::MIME"text/plain", d::MultiParticleNormal{P,S,B}) where {P, S, B}
    println(io, "MultiParticleNormal($(length(d.particles)) particles in $(size(d.cov, 1)) dimensions, scaling = $(round(d.scaling, sigdigits = 3)), regularization = $(round(d.regularization, sigdigits = 3)), $(B <: Nothing ? "without" : "with") bounds)")
end
function Distributions._rand!(rng::Random.AbstractRNG,
                              d::MultiParticleNormal,
                              x::AbstractVector{<:Real})
    idx = wsample(d.weights)
    Distributions._rand!(rng, d, x, idx)
end
function Distributions._rand!(rng::Random.AbstractRNG,
                              d::MultiParticleNormal{<:Any, <:Any, <:Nothing},
                              x::AbstractVector{<:Real},
                              idx::Int)
    Distributions.add!(unwhiten!(d.cov, randn!(rng, x)), d.particles[idx])
end
function Distributions._rand!(rng::Random.AbstractRNG,
                              d::MultiParticleNormal,
                              x::AbstractVector{<:Real},
                              idx::Int)
    while true
        Distributions.add!(unwhiten!(d.cov, randn!(rng, x)), d.particles[idx])
        for i in eachindex(x)
            (d.lower[i] > x[i] || d.upper[i] < x[i]) && break
            if i == length(x)
                return x
            end
        end
    end
        # To do this more efficiently
        # see e.g. https://github.com/BrianNaughton/TruncatedMVN.jl/blob/master/src/multivariate.jl
        # return gibbsTMVN(d.particles[idx], d.cov, d.lower, d.upper, diagm(ones(length(d.upper))), 1)
end
function Distributions._rand!(rng::Random.AbstractRNG,
                              d::MultiParticleNormal{<:Any, <:AbstractVector},
                              x::AbstractVector{<:Real},
                              idx::Int)
    particle = d.particles[idx]
    for i in eachindex(x)
        a = d.lower[i]
        b = d.upper[i]
        r = b - a
        m = particle[i]
        sig = d.cov[i]
        if r/sig <= 2.5066282746310002
            denom = 1/(2*sig^2)
            while true
                tmp = rand() * r + a
                if exp(-(tmp-m)^2*denom) > rand()
                    x[i] = tmp
                    break
                end
            end
        else
            while true
                tmp = m + randn() * sig
                if a <= tmp && b >= tmp
                    x[i] = tmp
                    break
                end
            end
        end
    end
    x
end
function _lpdf(d, x)
    log(sum(d.weights[i] * exp(logpdf(MultivariateNormal(d.particles[i], d.cov), x))
            for i in eachindex(d.weights)))
end
function Distributions._logpdf(d::MultiParticleNormal{<:Any, <:Any, <:Nothing},
                               x::AbstractVector{<:Real})
    _lpdf(d, x)
end
function Distributions._logpdf(d::MultiParticleNormal, x::AbstractVector{<:Real})
    # This is wrong if upper and lower are finite; but it doesn't matter here.
    # A constant is missing due to truncation.
    for i in eachindex(x)
        (d.lower[i] > x[i] || d.upper[i] < x[i]) && return -Inf
    end
    _lpdf(d, x)
end
function Distributions.mean(d::MultiParticleNormal)
    sum(d.weights[i] .* d.particles[i] for i in eachindex(d.weights))
end
function update!(d::MultiParticleNormal{<:Any, S}, particles, weights; log_weights = true) where S
    log_weights ? weights!(d.weights, weights) : d.weights .= weights
    for i in eachindex(particles)
        d.particles[i] .= particles[i]
    end
    d.cov = cov_matrix(S, d.particles, d.weights, d.scaling, d.regularization)
    d
end

lower(d) = fill(-Inf, length(d))
upper(d) = fill(Inf, length(d))


text:
export const_controller

function const_controller(cfg::Dict; c::Float64=1.0)
    nouts = output_lengths(cfg)
    cell_division(x::Array{Float64}) = Bool(round(c))
    new_cell_params(x::Array{Float64}) = ones(nouts[2]) .* c
    cell_death(x::Array{Float64}) = Bool(round(c))
    cell_state_update(x::Array{Float64}) = ones(nouts[4]) .* c
    cell_param_update(x::Array{Float64}) = ones(nouts[5]) .* c
    connect(x::Array{Float64}) = Bool(round(c))
    new_conn_params(x::Array{Float64}) = ones(nouts[7]) .* c
    disconnect(x::Array{Float64}) = Bool(round(c))
    conn_state_update(x::Array{Float64}) = ones(nouts[9]) .* c
    conn_param_update(x::Array{Float64}) = ones(nouts[10]) .* c

    Controller(cell_division, new_cell_params, cell_death,
               cell_state_update, cell_param_update,
               connect, new_conn_params, disconnect,
               conn_state_update, conn_param_update)
end


text:
module UtilsTests

using TensorPolynomialBases: _gradient_type
using Test
using StaticArrays
using TensorValues

@test _gradient_type(SVector{3,Int},SVector{2,Int}) == SMatrix{2,3,Int,6}

@test _gradient_type(SMatrix{3,4,Int},SVector{2,Int}) == SArray{Tuple{2,3,4},Int,3,24}

@test _gradient_type(VectorValue{3,Int},VectorValue{2,Int}) == MultiValue{Tuple{2,3},Int,2,6}

@test _gradient_type(TensorValue{3,Int,9},VectorValue{2,Int}) == MultiValue{Tuple{2,3,3},Int,3,18}

end # module UtilsTests


text:
using CSV, DataFrames, MLDataUtils, KernelFunctions, MVDA, StableRNGs
using LinearAlgebra, Statistics

# using MKL
BLAS.set_num_threads(8)

run = function(dir, example, data, sparse2dense::Bool=false; at::Real=0.8, nfolds::Int=5, kwargs...)
    # Create MVDAProblem instance w/o kernel and construct grids.
    @info "Creating MVDAProblem for $(example)"
    problem = MVDAProblem(data..., intercept=true, kernel=nothing)
    n, p, c = MVDA.probdims(problem)
    n_train = round(Int, n * at * (nfolds-1)/nfolds)
    n_validate = round(Int, n * at * 1/nfolds)
    n_test = round(Int, n * (1-at))
    fill!(problem.coeff.all, 1/(p+1))
    fill!(problem.coeff_prev.all, 1/(p+1))
    ϵ_grid = [MVDA.maximal_deadzone(problem)]
    s_grid = sort!([1-k/p for k in p:-1:0], rev=sparse2dense)
    grids = (ϵ_grid, s_grid)

    # Collect data for cross-validation replicates.
    @info "CV split: $(n_train) Train / $(n_validate) Validate / $(n_test) Test"
    subsets = (n_train, n_validate, n_test)
    rng = StableRNG(1903)
    replicates = MVDA.cv_estimation(MMSVD(), problem, grids;
        at=at,                  # propagate CV / Test split
        nfolds=nfolds,          # propagate number of folds
        rng=rng,                # random number generator for reproducibility
        nouter=10^2,            # outer iterations
        ninner=10^6,            # inner iterations
        gtol=1e-3,              # tolerance on gradient for convergence of inner problem
        dtol=1e-3,              # tolerance on distance for convergence of outer problem
        rtol=0.0,               # use strict distance criteria
        lambda=1e-3,            # regularization level
        tol=1e-6,               # convergence tolerance in initialization
        maxiter=10^4,           # maximum iterations in initialization
        nesterov_threshold=100, # delay on Nesterov acceleration
        show_progress=true,     # display progress over replicates
        kwargs...               # propagate other keywords
    )

    # Write results to disk.
    @info "Writing to disk"
    traversal = sparse2dense ? "S2D" : "D2S"
    filename = joinpath(dir, "$(example)-L-path=$(traversal).dat")
    kernel = "none"
    problem_info = (example, subsets, p, c, sparse2dense, kernel)
    MVDA.save_cv_results(filename, problem_info, replicates, grids)

    return nothing
end

# Examples ordered from easiest to hardest
examples = (
    ("colon", 3, 0.8),
    ("srbctA", 3, 0.8),
    ("leukemiaA", 3, 0.8),
    ("lymphomaA", 3, 0.8),
    ("brain", 3, 0.8),
    ("prostate", 3, 0.8),
)

dir = ARGS[1]
@info "Output directory: $(dir)"

for (example, nfolds, split) in examples
    tmp = CSV.read("/home/alanderos/Desktop/data/$(example).DAT", DataFrame, header=false)
    df = shuffleobs(tmp, rng=StableRNG(1234))
    data = (Vector(df[!,end]), Matrix{Float64}(df[!,1:end-1]))
    run(dir, example, data, false;
        at=split,           # CV set / Test set split
        nfolds=nfolds,      # number of folds
        nreplicates=50,     # number of CV replicates
    )
end


text:
using BenchmarkTools, Compat, DataFrames, IntervalRootFinding, IntervalArithmetic, StaticArrays

function randVec(n::Int)
    a = randn(n)
    A = Interval.(a)
    mA = MVector{n}(A)
    sA = SVector{n}(A)
    return A, mA, sA
end

function randMat(n::Int)
    a = randn(n, n)
    A = Interval.(a)
    mA = MMatrix{n, n}(A)
    sA = SMatrix{n, n}(A)
    return A, mA, sA
end

function benchmark(max=10)
    df = DataFrame()
    df[:Method] = ["Array", "MArray", "SArray", "Contractor", "ContractorMArray", "ContractorSArray"]
    for n in 1:max
        A, mA, sA = randMat(n)
        b, mb, sb = randVec(n)
        t1 = @belapsed gauss_seidel_interval($A, $b)
        t2 = @belapsed gauss_seidel_interval($mA, $mb)
        t3 = @belapsed gauss_seidel_interval($sA, $sb)
        t4 = @belapsed gauss_seidel_contractor($A, $b)
        t5 = @belapsed gauss_seidel_contractor($mA, $mb)
        t6 = @belapsed gauss_seidel_contractor($sA, $sb)
        df[Symbol("$n")] = [t1, t2, t3, t4, t5, t6]
    end
    a = []
    for i in 1:max
        push!(a, Symbol("$i"))
    end
    df1 = stack(df, a)
    dfnew = unstack(df1, :variable, :Method, :value)
    dfnew = rename(dfnew, :variable => :n)
    println(dfnew)
    dfnew
end

function benchmark_elimination(max=10)
    df = DataFrame()
    df[:Method] = ["Gauss Elimination", "Base.\\"]
    for n in 1:max
        A, mA, sA = randMat(n)
        b, mb, sb = randVec(n)
        t1 = @belapsed gauss_elimination_interval($A, $b)
        t2 = @belapsed gauss_elimination_interval1($A, $b)
        df[Symbol("$n")] = [t1, t2]
    end
    a = []
    for i in 1:max
        push!(a, Symbol("$i"))
    end
    df1 = stack(df, a)
    dfnew = unstack(df1, :variable, :Method, :value)
    dfnew = rename(dfnew, :variable => :n)
    println(dfnew)
    dfnew
end


text:
function solve()
    n, x = [parse(Int, x) for x in split(readline())]
    as = [parse(Int, x) for x in split(readline())]

    join([a for a in as if a != x], " ")
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
# Julia wrapper for header: /usr/include/scip/scip_sepa.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeSepa(scip, name, desc, priority, freq, maxbounddist, usessubscip, delay, sepacopy, sepafree, sepainit, sepaexit, sepainitsol, sepaexitsol, sepaexeclp, sepaexecsol, sepadata)
    ccall((:SCIPincludeSepa, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Cstring, Cint, Cint, Cdouble, UInt32, UInt32, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{SCIP_SEPADATA}), scip, name, desc, priority, freq, maxbounddist, usessubscip, delay, sepacopy, sepafree, sepainit, sepaexit, sepainitsol, sepaexitsol, sepaexeclp, sepaexecsol, sepadata)
end

function SCIPincludeSepaBasic(scip, sepa, name, desc, priority, freq, maxbounddist, usessubscip, delay, sepaexeclp, sepaexecsol, sepadata)
    ccall((:SCIPincludeSepaBasic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_SEPA}}, Cstring, Cstring, Cint, Cint, Cdouble, UInt32, UInt32, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{SCIP_SEPADATA}), scip, sepa, name, desc, priority, freq, maxbounddist, usessubscip, delay, sepaexeclp, sepaexecsol, sepadata)
end

function SCIPsetSepaCopy(scip, sepa, sepacopy)
    ccall((:SCIPsetSepaCopy, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_SEPA}, Ptr{Cvoid}), scip, sepa, sepacopy)
end

function SCIPsetSepaFree(scip, sepa, sepafree)
    ccall((:SCIPsetSepaFree, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_SEPA}, Ptr{Cvoid}), scip, sepa, sepafree)
end

function SCIPsetSepaInit(scip, sepa, sepainit)
    ccall((:SCIPsetSepaInit, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_SEPA}, Ptr{Cvoid}), scip, sepa, sepainit)
end

function SCIPsetSepaExit(scip, sepa, sepaexit)
    ccall((:SCIPsetSepaExit, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_SEPA}, Ptr{Cvoid}), scip, sepa, sepaexit)
end

function SCIPsetSepaInitsol(scip, sepa, sepainitsol)
    ccall((:SCIPsetSepaInitsol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_SEPA}, Ptr{Cvoid}), scip, sepa, sepainitsol)
end

function SCIPsetSepaExitsol(scip, sepa, sepaexitsol)
    ccall((:SCIPsetSepaExitsol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_SEPA}, Ptr{Cvoid}), scip, sepa, sepaexitsol)
end

function SCIPfindSepa(scip, name)
    ccall((:SCIPfindSepa, libscip), Ptr{SCIP_SEPA}, (Ptr{SCIP_}, Cstring), scip, name)
end

function SCIPgetSepas(scip)
    ccall((:SCIPgetSepas, libscip), Ptr{Ptr{SCIP_SEPA}}, (Ptr{SCIP_},), scip)
end

function SCIPgetNSepas(scip)
    ccall((:SCIPgetNSepas, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPsetSepaPriority(scip, sepa, priority)
    ccall((:SCIPsetSepaPriority, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_SEPA}, Cint), scip, sepa, priority)
end

function SCIPgetSepaMinEfficacy(scip)
    ccall((:SCIPgetSepaMinEfficacy, libscip), Cdouble, (Ptr{SCIP_},), scip)
end


text:
# Description:
#   Declare the Euler data object
#   Includes all the files for the Euler module

module EulerEquationMod

using PDESolver  # setup LOAD_PATH to find all PDESolver components
using SolverCommon
#using ArrayViews
using ODLCommonTools  # abstract type definitions + common functions
using SummationByParts
using PdePumiInterface  # common mesh interface implementation - pumi
#using NonlinearSolvers
using Jacobian
using LinearSolvers
using Utils
import ODLCommonTools.sview
using MPI
using Input  # input file processing
using PETSc2
#using PETSc
# using FreeFormDeformation
# using MeshMovement

# the AbstractEquation type is declared in ODLCommonTools
# every equation will have to declare a new type that is a subtype of AbstractEquation

export AbstractEulerData, EulerData, EulerData_, run_euler, eval_dJdaoa

# add a layer of abstraction
@doc """
### EulerEquationMod.AbstractEulerData{Tsol, Tres}

  This abstract type should be the supertype of *all* solution data objects
  that are related to the Euler equations.

  It should be used for specifying the type of a function argument only when
  the function does no operations on the solution data object itself, it just
  passes it onto other functions that do the work (thus AbstractEulerData
  should be used for only the highest level functions).

  Another way of saying say it is that this type should only be used when
  the function only needs to ensure that it is solving the Euler equations,
  but does not care even a little bit about how.
"""->
abstract type AbstractEulerData{Tsol, Tres} <: AbstractSolutionData{Tsol, Tres} end
#=
Use this type to leverage multiple dispatch.
This type holds any data, including large arrays of solution values,
  that are specific to the equation for the Euler equations.
  This includes the solution variables q, the fluxes in the xi and eta
  direction, and the result of the calculation the inverse mass matrix is
  stored here as well (not sure if that fits better here or in the mesh object)
  things like the coordinate field, the jacobian etc. are stored
  in the mesh object

The aux_vars array holds all auxiliary variables that are stored over
  the entire mesh.
Although it can be accessed directly, the preferred method is to use the macros
  defined in the euler_macros.jl file.
These macros return either a scalar or an ArrayView of the specified indices,
  depending on if the quantity requested is scalar or vector.
Every time a new variable is added to the array, the size must be updated
  and a new macro must be created.

The uses of aux_vars should mirror that of eqn.q, in that entire columns
  should be passed to low level functions and the low level functions
  use the macros to access individual variables.
The advantages of macros vs functions for access to variables remains unclear
  if aux_vars is a fixed size.
If it is variable sized then macros give the advantage of doing location lookup
  at compile time.
=#

#=
  Initial Condition:
    All initial condition functions use conservative variables, and are
    converted to the variables used to solve the equation.
=#

#=
  Boundary conditions:
    All boundary condition calculations are done in conservative variables.
    If using some other variables, they are converted before being passed
    to the boundary condition functions.

    The functors that compute boundary conditions are gathered from BCDict
    during init() and stored in mesh.bndry_funcs.  The functors compute the
    flux from the boundary condition at a node.  The boundary fluxes for all
    the nodes on the boundary are stored in mesh.bndryflux, to be integrated
    later using boundaryintegrate!.
=#

#=
  Variable Conversion:
    If solving the equations is some set of variables other than conservative,
    (lets call them the v variables, and the conservative variables the q
     variables), the following functions must be defined:

    # convert to conservative variables at a node
    convertToConservative(params::ParamType, qe::AbstractArray{Tsol, 1},
                          qc::AbstractArray{Tsol, 1}

    # convert from conservative to the other variables at a node
    convertToEntropy(params::ParamType, qe::AbstractArray{Tsol, 1},
                     qc::AbstractArray{Tsol, 1})

    # convert from conservative variables to the variables the equation
    # is being solved in
    convertFromConservativeToWorkingVars(params::ParamType{Tdim, :var_type,
                        qe::AbstractArray{Tsol, 1}, qc::AbstractArray{Tsol, 1})

    # calculate the coefficient matrix of the time derivate (dq/dv) at a node
    calcA0(params::ParamType, q::AbstractArray{Tsol, 1},
           A0::AbstractArray{Tsol, 2})
    # calculate inv(dq/dv)
    calcA0inv(params::ParamType, q::AbstractArray{Tsol, 1},
              A0::AbstractArray{Tsol, 2})

    # multiply a 3D array by inv(A0)
    matVecA0Inv(mesh, sbp, eqn, opts, arr::AbstractArray{Tsol, 3})

    # calculate the Euler flux at a node
    calcEulerFlux(params::ParamType, q::AbstractArray{Tsol, 1},
                  aux_vars::AbstractArray{Tsol, 1}, dir::AbstractArray{Tmsh},
                  F::AbstractArray{Tsol, 1})

=#

#=
  Calculating Other Quantities:
    To calculate physical quantities from the solution variables (for which
    the calculation is different for the different sets of variables) at a node,
    a separate function should be used, usnig the params object to dispatch
    to the right function for the variables being used.

    Some examples:

    calcPressure(params::ParamType, q::AbstractArray{Tsol, 1})
    calcEntropy(params::ParamType, q::AbstractArray{Tsol, 1})
    calcSpeedOfSound(params::ParamType, q::AbstractArray{Tsol, 1})
=#

@doc """
### EulerEquationMod.EulerData

  This type, although abstract, is the type functions should use for their
  input arguments if they do any operations on the solution data object.
  It stores all data used in evaluting the Euler Equations.

  It is paramaterized on the types Tsol, the type of the
  conservative variables q, and Tdim, the dimension of the equation

  It should have the following fields:

    * res_type : datatype of residual (depreciated)
    * q  : 3D array holding conservative variables
    * q_vec  : vector to assemble q into
    * aux_vars : 3D array holding auxiliary variables
    * flux_parametric : 4D array [ndof per node, nnodes per element, nelements, Tdim]
             holding the Euler flux in the xi and eta directions
    * res  : 3D array holding residual
    * res_vec   : vector form of res
    * edgestab_alpha : paramater used for edge stabilization, 4d array
    * bndryflux : 3D array holding boundary flux data
    * stabscale : 2D array holding edge stabilization scale factor
    * M : vector holding the mass matrix
    * Minv :  vector holding inverse mass matrix
    # Minv3D :  3D array holding inverse mass matrix for application to res (not res_vec)
"""->
abstract type EulerData{Tsol, Tres, Tdim, var_type} <: AbstractEulerData{Tsol, Tres} end

"""
  Functor type for faceElementIntegrals.  These integrals operate on a face,
  but require data from the entirety of the elements that make up the
  face, rather than data interpolated to the face

  These functors have the signature:

  ```
    calcFaceElementIntegral(obj::FaceElementIntegralType,
        params::AbstractParamType, sbpface::AbstractFace, iface::Interface,
        qL::AbstractMatrix, qR::AbstractMatrix, aux_vars::AbstractMatrix,
        nrm_face::AbstractMatrix, functor::FluxType,
        resL::AbstractMatrix, resR::AbstractMatrix)
  ```

  where `qL` and `qR` contains the solutions at the volume nodes of the left
  and right elements (`numDofPerNode` x `numNodesPerElement`), `aux_vars`
  contains the auxiliary variables for `qL`, `nrm_face` contains the normal
  vector at each face quadracture note, `functor` is the flux function used
  for the entropy conservative integrals, and `resL` and `resR` are the
  residual for the left and right elements, to be updated (same shape as
  `qL`, `qR`).

  Some of the `FaceElementIntegralType`s compute the entropy conservative
  integrals, others apply an entropy dissipative penalty of one kind
  or another.

  For use with explicit jacobian calculation, functors should also define
  a method for:

  ```
    calcFaceElementIntegral_diff(obj::FaceElementIntegralType,
        params::AbstractParamType, sbpface::AbstractFace, iface::Interface,
        qL::AbstractMatrix, qR::AbstractMatrix, aux_vars::AbstractMatrix,
        nrm_face::AbstractMatrix, functor::FluxType_diff,
        jacLL::AbstractArray{T,4}, resLR::AbstractArray{T,4},
        jacRL::AbstractArray{T,4}, resRR::AbstractArray{T,4})
  ```

  Most of the arguments are the same as the regular version, the output
  arguments are

   **Inputs/Outputs**

    * jacLL: jacobian of `resL` wrt `qL`
    * jacLR: jacobian of `resL` wrt `qR`
    * jacRL: jacobian of `resR` wrt `qL`
    * jacRR: jacobian of `resR` wrt `qR`

"""
abstract type FaceElementIntegralType end

"""
  Abstract type for all kernel operations used with entropy penalty functions
  (ie. given the state at the interface, apply a symmetric semi-definite
  operation).  Every type must extend this function with a new method:

  ```
    applyEntropyKernel(obj::AbstractEntropyKernel, params::ParamType,
                        q_avg::AbstractVector, delta_w::AbstractVector,
                        nrm_in::AbstractVector, flux::AbstractVector)
  ```

  specializing the first argument.  Here `q_avg` is the state at the interface
  in conservative variables, `delta_w` is the vector to multiply against,
  `nrm_in` is the normal vector at the node, and `flux` is the vector to
  overwrite with the result.

  Each `AbstractEntropyKernel` should have an outer constructor with signaturea

  ```
    function MyEntropyKernel(mesh::AbstractMesh, eqn::EulerData)
  ```

  For use in computing the Jacobian, a second function should be extended with
  a new method

  ```
    function applyEntropyKernel_diff(obj::AbstractEntropyKernel, params::ParamType,
                          q_avg::AbstractVector, q_avg_dot::AbstractMatrix,
                          delta_w::AbstractVector, delta_w_dot::AbstractMatrix,
                          nrm::AbstractVector,
                          flux::AbstractVector, flux_dot::AbstractMatrix)

  ```

  **Inputs**

   * obj: `AbstractEntropyKernel` implementation
   * params: ParamType
   * q_avg: same as regular method
   * q_avg_dot: `numDofPerNode` x `nd` matrix containing `nd` dual vectors for
                `q_avg`
   * delta_w: same as regular method
   * delta_w_dot: `numDofPerNode` x `nd` matrix containing `nd` dual vectors for
                  `delta_w`
   * nrm: same as regular method
   * flux: vector to be overwritten with the flux
   * flux_dot: `numDofPerNode` x `nd` array the derivative values will be added
               to.

  Note that `nd` must be less than or equal to the `nd` argument of the
  `AbstractEntropyKernel` constructor.

  For reverse mode with respect to the solution, a function should be defined:

  ```
    function applyEntropyKernel_revq(obj::LFKernel, params::ParamType,
                            q_avg::AbstractVector, q_bar::AbstractVector,
                            delta_w::AbstractVector, delta_w_bar::AbstractVector,
                            nrm::AbstractVector, flux::AbstractVector,
                            flux_bar::AbstractVector)
  ```

  **Inputs**

   * obj: `AbstractEntropyKernel` implementation
   * params: ParamType
   * q_avg: same as regular method
   * delta_w: same as regular method
   * nrm: same as regular method
   * flux_bar: seed vector for dual part of `flux`

  **Inputs/Outputs**

   * q_bar: vector to accumulate the result for `q_avg` into (not overwritten)
   * delta_w_bar: vector to accumulate result for `delta_w` into (not overwritten)
   * flux: vector, length `numDofPerNode` to overwrite with flux

  Similarly for reverse mode with respect to the metrics:

  ```
    function applyEntropyKernel_revm(obj::AbstractEntropyKernel, params::ParamType,
                            q_avg::AbstractVector, delta_w::AbstractVector,
                            nrm::AbstractVector, nrm_bar::AbstractVector,
                            flux::AbstractVector,
                            flux_bar::AbstractVector)
  ```

  **Inputs**

   * obj: `AbstractEntropyKernel` implementation
   * params: ParamType
   * q_avg: same as regular method
   * delta_w: same as regular method
   * nrm: same as regular method
   * flux_bar: seed vector for dual part of `flux`

  **Inputs/Outputs**

   * nrm_bar: vector to accumulate result into (not overwritten)
   * flux: vector, length `numDofPerNode` to overwrite with flux

"""
abstract type AbstractEntropyKernel end

"""
  Functionals that use entropy penalties
"""
abstract type EntropyPenaltyFunctional{Topt} <: AbstractIntegralFunctional{Topt} end

"""
  Abstract type for Local Discontinuous Galerkin flux functions
"""
abstract type AbstractLDGFlux end

# high level functions should take in an AbstractEulerData, remaining
# agnostic to the dimensionality of the equation
# Mid level function should take in an EulerData{Tsol, Tdim}, so they
# know the dimensionality of the equation, but have a single method that
# can handle all dimensions

# low level functions should take in EulerData{Tsol, 2} or EulerData{Tsol, 3}
# this allows them to have different methods for different dimension equations.

# TEMPORARY
const Abstract4DArray{T} = AbstractArray{T, 4}
const Abstract5DArray{T} = AbstractArray{T, 5}

import PDESolver: getShockSensor, setShockSensor, setShockSensorAlpha

include("types.jl")  # type definitions
include(joinpath(Pkg.dir("PDESolver"), "src/solver/debug.jl"))  # debug macro
include("euler_macros.jl")
include("common_funcs.jl")
include("euler_funcs.jl")
include("conversion.jl")
include("euler.jl")
include("ic.jl")
include("bc.jl")
include("stabilization.jl")
include("faceElementIntegrals.jl")
include("flux.jl")
# include("artificialViscosity.jl")
# include("constant_diff.jl")
include("GLS2.jl")
include("functionals.jl")
include("boundary_functional.jl")
include("functional_deriv.jl")
include("evalFunctional.jl")
include("source.jl")
include("PressureMod.jl")
include("entropy_flux.jl")
include("eigensystem.jl")
include("check_options.jl")
include("eqn_deepcopy.jl")
include("startup_func.jl")  # function for invoking the solver
include("evaldRdm.jl")
include("evaldRdq.jl")
include("homotopy.jl")
include("shock_capturing_mesh.jl")
include("shock_diffusion.jl")
include("shock_capturing.jl")
include("shock_sensors.jl")
include("projection_shock_capturing.jl")
include("volume_shock_capturing.jl")
include("ldg_shock_capturing.jl")
include("sbp_sat_shock_capturing.jl")
include("sbp_sat_reduced_sc.jl")
include("sbp_sat_reduced2_sc.jl")
include("sbp_cartesian.jl")
include("lp_stab.jl")
include("face_functional_integrate.jl")

# Jacobian calculation
include("evalJacobian.jl")
include("euler_funcs_diff.jl")
include("flux_diff.jl")
include("bc_solvers_diff.jl")
include("bc_diff.jl")
include("homotopy_diff.jl")
include("shock_sensors_diff.jl")
include("shock_capturing_diff.jl")
include("sbp_cartesian_diff.jl")
include("sbp_sat_shock_capturing_diff.jl")
include("sbp_sat_reduced2_sc_diff.jl")

"""
  This physics is named `Euler`
"""
global const PhysicsName = "Euler"

register_physics(PhysicsName, EulerEquationMod, createObjects, checkOptions)

end # end module


text:
using AutomotiveDrivingModels
using Parameters

const MAIN_LANE_ID = 1
const MERGE_LANE_ID = 2

"""
    MergingEnvironment

A road network with a main lane and a merging lane. The geometry can be modified by
passing the parameters as keyword arguments in the constructor

# Parameters
- `lane_width::Float64 = 3.0`
- `main_lane_vmax::Float64 = 15.0`
- `merge_lane_vmax::Float64 = 10.0`
- `main_lane_length::Float64 = 100.0`
- `main_lane_angle::Float64 = float(pi)/4`
- `merge_lane_angle::Float64 = float(pi)/4`
- `merge_lane_length::Float64 = 50.0`
- `after_merge_length::Float64 = 50.0`

# Internals

- `roadway::Roadway{Float64}` contains all the road segment and lane information
- `merge_point::VecSE2{Float64}` coordinate of the merge point in cartesian frame (0.0, 0.0, 0.0) by default
- `merge_proj::RoadProjection{Int64, Float64}` projection of the merge point on the roadway
- `merge_index::RoadIndex`

"""
@with_kw struct MergingEnvironment
    lane_width::Float64 = 3.0
    main_lane_vmax::Float64 = 15.0
    merge_lane_vmax::Float64 = 10.0
    main_lane_length::Float64 = 100.0
    main_lane_angle::Float64 = float(pi)/4
    merge_lane_angle::Float64 = float(pi)/4
    merge_lane_length::Float64 = 50.0
    after_merge_length::Float64 = 50.0

    # internals
    roadway::Roadway{Float64} = generate_merging_roadway(lane_width,
                                                main_lane_vmax,
                                                merge_lane_vmax,
                                                main_lane_length,
                                                merge_lane_length,
                                                after_merge_length,
                                                main_lane_angle,
                                                merge_lane_angle)
    merge_point::VecSE2{Float64} = VecSE2(0.0, 0.0, 0.0)
    merge_proj::RoadProjection{Int64, Float64} = proj(merge_point, roadway)
    merge_index::RoadIndex{Int64, Float64} = RoadIndex(merge_proj.curveproj.ind, merge_proj.tag)
end

function append_to_curve!(target::Curve{T}, newstuff::Curve{T}) where T <: Real
    s_end = target[end].s
    for c in newstuff
        push!(target, CurvePt{T}(c.pos, c.s+s_end, c.k, c.kd))
    end
    return target
end

"""
    generate_merging_roadway(lane_width::Float64 = 3.0, main_lane_vmax::Float64 = 20.0, merge_lane_vmax::Float64 = 15.0, main_lane_length::Float64 = 20.0, merge_lane_length::Float64 = 20.0, after_merge_length::Float64 = 20.0, main_lane_angle::Float64 = float(pi)/4, merge_lane_angle::Float64 = float(pi)/4)

Generate a `Roadway` object representing a merging scenario.
The merge point is defined at (0., 0.) by default.
"""
function generate_merging_roadway(lane_width::Float64 = 3.0,
                                   main_lane_vmax::Float64 = 20.0,
                                   merge_lane_vmax::Float64 = 15.0,
                                   main_lane_length::Float64 = 20.0,
                                   merge_lane_length::Float64 = 20.0,
                                   after_merge_length::Float64 = 20.0,
                                   main_lane_angle::Float64 = float(pi)/4,
                                   merge_lane_angle::Float64 = float(pi)/4)
    # init empty roadway
    roadway = Roadway()
    n_pts = 2 # sample points for the roadway, only two needed each time, since all segments are straight
    main_tag = LaneTag(MAIN_LANE_ID, 1)
    merge_tag = LaneTag(MERGE_LANE_ID, 1)
    # after_merge_tag = LaneTag(AFTER_MERGE_LANE_ID, 1)

    # define curves
    merge_point = VecE2(0.0, 0.0)
    main_lane_startpt = merge_point + polar(main_lane_length, -float(pi) - main_lane_angle)
    main_curve = gen_straight_curve(main_lane_startpt, merge_point, n_pts)
    merge_index = curveindex_end(main_curve)
    append_to_curve!(main_curve, gen_straight_curve(merge_point, merge_point + polar(after_merge_length, 0.0), n_pts)[2:end])
    merge_lane_startpt = merge_point + polar(merge_lane_length, float(pi) + merge_lane_angle)
    merge_curve = gen_straight_curve(merge_lane_startpt, merge_point, n_pts)

    # define lanes with connections
    main_lane = Lane(main_tag, main_curve, width = lane_width, speed_limit=SpeedLimit(0.,main_lane_vmax))
    merge_lane = Lane(merge_tag, merge_curve, width = lane_width,speed_limit=SpeedLimit(0.,merge_lane_vmax),
                        next=RoadIndex(merge_index, main_tag))

    # add segments to roadway
    push!(roadway.segments, RoadSegment(MAIN_LANE_ID, [main_lane]))
    push!(roadway.segments, RoadSegment(MERGE_LANE_ID, [merge_lane]))

    return roadway
end

## convenience
"""
    main_lane(env::MergingEnvironment)
returns the main lane of the merging scenario
"""
main_lane(env::MergingEnvironment) = env.roadway[LaneTag(MAIN_LANE_ID, 1)]

"""
    merge_lane(env::MergingEnvironment)
returns the merging lane of the merging scenario
"""
merge_lane(env::MergingEnvironment) = env.roadway[LaneTag(MERGE_LANE_ID, 1)]


text:
function test_header(hdStr)
    println("-----  Testing $hdStr")
end

function test_file_dir()
    return joinpath(@__DIR__, "test_files")
end

if !isdir(test_file_dir())
    mkdir(test_file_dir())
end

# ------------


text:
function manifoldupdate!(integ, residualf; maxiters=100, ϵ₁=1e-25, ϵ₂=1e-15)
    @unpack x, SolProj = integ.cache
    f(m) = residualf(SolProj * m)
    return copy!(x, manifoldupdate(x, f; maxiters=maxiters, ϵ₁=ϵ₁, ϵ₂=ϵ₂))
end

function manifoldupdate(x, f; maxiters=100, ϵ₁=1e-25, ϵ₂=1e-15)
    m, C = x
    m_i = copy(m)

    local m_i_new, C_i_new
    for i in 1:maxiters
        z = f(m_i)
        J = ForwardDiff.jacobian(f, m_i)
        S = X_A_Xt(C, J)

        m_i_new, C_i_new = update(x, Gaussian(z .+ (J * (m - m_i)), S), J)

        if norm(m_i_new .- m_i) < ϵ₁ && norm(z) < ϵ₂
            break
        end
        m_i = m_i_new
    end
    return Gaussian(m_i_new, C_i_new)
end

function ManifoldUpdate(residualf, args...; maxiters=100, ϵ₁=1e-25, ϵ₂=1e-15, kwargs...)
    condition(u, t, integ) = true
    affect!(integ) = manifoldupdate!(integ, residualf; maxiters=maxiters, ϵ₁=ϵ₁, ϵ₂=ϵ₂)
    return DiscreteCallback(condition, affect!, args...; kwargs...)
end


text:
include("src/base.jl")
include("src/envs/toytext_mdp.jl")
include("src/policies/value_table_policy.jl")
include("src/solvers/value_table_policy_solver.jl")

gym_env = ToyTextMDP(gym.make("Taxi-v2"))
#gym_env = ToyTextMDP(gym.make("FrozenLake-v0"))

#solver = ValueTablePolicySolver(5000; max_frame_iterations = 250, epsilon = 0.01, learning_rate = 0.85, discount = 0.99)
solver = ValueTablePolicySolver(5000; max_frame_iterations = 100, epsilon = 0.25, epsilon_discount = 0.999, learning_rate = 0.15, discount = 0.99)
policy = solve(solver, gym_env, verbose = true)

# println(policy.value_map)

max_stars_repo_name: dmitrijsc/practical-rl

max_issues_repo_name: dmitrijsc/practical-rl

max_forks_repo_name: dmitrijsc/practical-rl


text:
# main simulation module for evaluating dynamics and creating the world
module sim

dir_path = @__DIR__
include_dependency(dir_path * "/../vis/vis.jl")

using LinearAlgebra
push!(LOAD_PATH, dir_path * "/../vis")
using vis

include("rk4.jl")
include("sim_figure8.jl")
include("sim_types.jl")

export Road, Agent, World
export advance!, make_figure8_path, diagnostic

# Custom Rendering ############################################################
function diagnostic(agent::Agent, world::World, t::Float64=0.0)
  (x, y) = sp2xy(agent.x[1], agent.x[3], world.road.path)
  (sx, sy) = sp2sxsy(agent.x[1], agent.x[3], world.road.path)

  dx = fill(0.0, 3)
  u = fill(0.0, 2)
  agent.dynamics!(dx, agent.x, Pair(agent, world), t)
  agent.controller!(u, agent.x, dx, Pair(agent, world), t)

  return (x, y, sx, sy, dx, u)
end

# Dynamics and Control ########################################################
function default_controller!(u::AbstractArray{Float64},
                             x::AbstractArray{Float64},
                             dx::AbstractArray{Float64},
                             agent_world::Pair{Agent, World}, t::Float64)
  u[1] = -0.1 * (x[2] - 120 / 3.6)
  u[2] = -0.01 * x[3] / dx[1] # counteract lateral movement

  return
end

function default_dynamics!(dx::AbstractArray{Float64},
                           x::AbstractArray{Float64},
                           agent_world::Pair{Agent, World}, t::Float64)
  agent = agent_world.first
  world = agent_world.second

  # extract states
  s = x[1]
  ds = x[2]
  p = x[3]

  # get local interpolated path curvature
  path = world.road.path

  s = mod(s, path.S[end])
  idx = binary_search(path.S, s)
  if idx + 1 > length(path.R)
    display(path.S[end])
    display(s)
  end
  R = path.R[idx]
  R1 = path.R[idx + 1]
  S = path.S[idx]
  S1 = path.S[idx + 1]
  r = (R1 - R) / (S1 - S) * (s - S) + R

  # compute rates
  ds *= isnan(r / (r - p)) ? 1 : r / (r - p)
  dds = 0.0
  dds += abs(p) > world.road.width / 2 ? -1e1 * ds : 0.0
  dp = isnan(ds^2 / r) ? 0.0 : -2e-1 * ds^2 / r

  dx[1] = ds
  dx[2] = dds
  dx[3] = dp

  # controller acts
  u = fill(0.0, 2)
  agent.controller!(u, x, dx, agent_world, t)

  dx[2] += u[1]
  dx[3] += isnan(u[2]) || isinf(u[2]) ? 0.0 : ds * u[2]

  # done
  return
end

# Utility Functions ###########################################################
function sp2xy(s::Float64, p::Float64, path)
  s = mod(s, path.S[end])
  idx = binary_search(path.S, s)

  X = path.X[idx]
  X1 = path.X[idx + 1]

  Y = path.Y[idx]
  Y1 = path.Y[idx + 1]

  S = path.S[idx]
  S1 = path.S[idx + 1]

  Px = path.Px[idx]
  Px1 = path.Px[idx + 1]

  Py = path.Py[idx]
  Py1 = path.Py[idx + 1]

  x = (X1 - X) / (S1 - S) * (s - S) + X
  y = (Y1 - Y) / (S1 - S) * (s - S) + Y

  x += ((Px1 - Px) / (S1 - S) * (s - S) + Px) * p
  y += ((Py1 - Py) / (S1 - S) * (s - S) + Py) * p

  return (x, y)
end

function sp2sxsy(s::Float64, p::Float64, path)
  s = mod(s, path.S[end])
  idx = binary_search(path.S, s)

  dX = path.dX[idx]
  dX1 = path.dX[idx + 1]

  dY = path.dY[idx]
  dY1 = path.dY[idx + 1]

  S = path.S[idx]
  S1 = path.S[idx + 1]

  sx = (dX1 - dX) / (S1 - S) * (s - S) + dX
  sy = (dY1 - dY) / (S1 - S) * (s - S) + dY

  return (sx, sy)
end

function binary_search(sorted_x::Array{T}, el::T) where T
  a = 1
  len = length(sorted_x)
  b = len
  if el > sorted_x[b]
    return b
  end

  m = div(a + b, 2)
  while b - a > 1
    if sorted_x[m] > el
      b = m
    else
      a = m
    end
    m = div(a + b, 2)
  end
  m = m == len ? m - 1 : m;

  return m
end

end


text:
####
#### Defines list of tutorials given `generated_dir`
####

generate_tutorials = true

tutorials = Any[]

# Allow flag to skip generated
# tutorials since this is by
# far the slowest part of the
# docs build.
if generate_tutorials

    tutorials_dir = joinpath(@__DIR__, "..", "tutorials")      # julia src files

    # generate tutorials
    import Literate

    tutorials_jl = [
        joinpath(root, f)
        for (root, dirs, files) in Base.Filesystem.walkdir(tutorials_dir)
        for f in files
    ]
    filter!(x -> endswith(x, ".jl"), tutorials_jl) # only grab .jl files

    filter!(x -> !occursin("topo.jl", x), tutorials_jl)                       # currently broken, TODO: Fix me!
    filter!(x -> !occursin("dry_rayleigh_benard.jl", x), tutorials_jl)        # currently broken, TODO: Fix me!
    filter!(x -> !occursin("ex_001_periodic_advection.jl", x), tutorials_jl)  # currently broken, TODO: Fix me!
    filter!(x -> !occursin("ex_002_solid_body_rotation.jl", x), tutorials_jl) # currently broken, TODO: Fix me!
    filter!(x -> !occursin("ex_003_acoustic_wave.jl", x), tutorials_jl)       # currently broken, TODO: Fix me!
    filter!(x -> !occursin("ex_004_nonnegative.jl", x), tutorials_jl)         # currently broken, TODO: Fix me!
    filter!(x -> !occursin("KinematicModel.jl", x), tutorials_jl)             # currently broken, TODO: Fix me!
    filter!(x -> !occursin("ex_1_saturation_adjustment.jl", x), tutorials_jl) # currently broken, TODO: Fix me!
    filter!(x -> !occursin("ex_2_Kessler.jl", x), tutorials_jl)               # currently broken, TODO: Fix me!

    println("Building literate tutorials:")
    for tutorial in tutorials_jl
        println("    $(tutorial)")
    end

    for tutorial in tutorials_jl
        gen_dir =
            joinpath(generated_dir, relpath(dirname(tutorial), tutorials_dir))
        input = abspath(tutorial)
        script = Literate.script(input, gen_dir)
        code = strip(read(script, String))
        mdpost(str) = replace(str, "@__CODE__" => code)
        Literate.markdown(input, gen_dir, postprocess = mdpost)
        # Literate.notebook(input, gen_dir, execute = true)
    end

    # TODO: Should we use AutoPages.jl?

    # These files mirror the .jl files in `CLIMA/tutorials/`:
    tutorials = Any[
        "Atmos" => Any[
            "Dry Idealized GCM" => "generated/Atmos/heldsuarez.md",
            "Rising Bubble LES" => "ExtendingCLIMA/Atmos/Model/risingbubble.md",
        ],
        "Ocean" => Any[],
        "Numerics" => Any[
            "LinearSolvers" => Any["Conjugate Gradient" => "generated/Numerics/LinearSolvers/cg.md",],
            "Contributing" => Any["Notes on Literate" => "generated/literate_markdown.md",],
        ],
    ]

end


text:
cd("lesson-modules")

using Pkg
Pkg.add("SpecialFunctions")
#%%
using SpecialFunctions

gamma(3)

sinint(5)

#%% restart the REPL befor running this piece of code
using SpecialFunctions: gamma, sinint

gamma(3)

sinint(5)

cosint(5)

#%%
function gamma(x)
    println("I am another 'gamma' function")
    return x^2
end

using SpecialFunctions

gamma(3)

SpecialFunctions.gamma(3)

#%% a better solution
import SpecialFunctions

gamma(3)

SpecialFunctions.gamma(3)

function gamma(x)
    println("I am another 'gamma' function")
    return x^2
end

gamma(3)

#%%

module MyModule
export func2

a = 42
function func1(x)
    return x^2
end

function func2(x)
    return func1(x) + a
end

end #end of module

#%%

using Main.MyModule

func2(3)
func1(3)
MyModule.func1(3)

#%% using MyBigModule
include("big-module.jl")

using .MyBigModule

func2big(3)

#%%

max_stars_repo_name: JeffreySarnoff/techytok-examples

max_issues_repo_name: JeffreySarnoff/techytok-examples

max_forks_repo_name: JeffreySarnoff/techytok-examples


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# https://discourse.julialang.org/t/how-to-read-only-the-last-line-of-a-file-txt/68005/6

# %%
str = """
some header option
header other option
random metadata
A,B,C,X,Y,Z
1,1,1,2.0,0.0,102.0
1,1,2,2.0,0.0,202.0
1,1,3,2.0,0.0,302.0
1,2,1,3.0,1.0,103.0
1,2,2,3.0,1.0,203.0
1,2,3,3.0,1.0,303.0
1,3,1,4.0,2.0,104.0
1,3,2,4.0,2.0,204.0
1,3,3,4.0,2.0,304.0
1,4,1,5.0,3.0,105.0
1,4,2,5.0,3.0,205.0
1,4,3,5.0,3.0,305.0
1,5,1,6.0,4.0,106.0
1,5,2,6.0,4.0,206.0
1,5,3,6.0,4.0,306.0
1,6,1,7.0,5.0,107.0
1,6,2,7.0,5.0,207.0
1,6,3,7.0,5.0,307.0
"""

write("test.txt", str)

# %%
function read_last(file)
    open(file) do io
        seekend(io)
        seek(io, position(io) - 1)
        while Char(peek(io)) != '\n'
            seek(io, position(io) - 1)
        end
        read(io, Char)
        read(io, String)
    end
end

# %%
function read_last2(file)
    open(file) do io
        seekend(io)
        seek(io, position(io) - 1)
        p = position(io)
        while Char(peek(io)) != '\n'
            seek(io, position(io) - 1)
        end
        if position(io) == p
            seek(io, position(io) - 1)
            while Char(peek(io)) != '\n'
                seek(io, position(io) - 1)
            end
        end
        read(io, Char)
        read(io, String)
    end |> chomp
end

# %%
read_last("test.txt")

# %%
read_last2("test.txt")

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public