text:
"""
This module defines the data types related to seismic data:
* `TD` : time domain representation

It also provides methods that apply source and receiver filters onto
seismic data.
"""
module Data

using Misfits
using LinearAlgebra
using Statistics
import GeoPhyInv.Interpolation
import GeoPhyInv.Acquisition
import GeoPhyInv.Utils
import GeoPhyInv.Coupling
using DSP
using Random

"""
Time domain representation of Seismic Data.

# Fields

* `d::Array{Array{Float64,2},2}` : data
* `fields::Vector{Symbol}` : components recorded at each receiver
* `tgrid` : grid to represent time
* `acqgeom::Acquisition.Geom` : acquisition geometry used to generate the data
"""
mutable struct TD
	d::Matrix{Matrix{Float64}}
	fields::Vector{Symbol}
	tgrid::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
	acqgeom::Acquisition.Geom
	"adding conditions that are to be false while construction"
	TD(d, fields, tgrid, acqgeom) =
		any([
		  any([fields[ifield] ∉ [:P, :Vx, :Vz] for ifield in 1:length(fields)]),
		  length(fields) == 0,
		  broadcast(size,d) != [(length(tgrid),acqgeom.nr[iss]) for iss=1:acqgeom.nss, ifield=1:length(fields)]
		  ]) ?
		error("error in TD construction") : new(d, fields, tgrid, acqgeom)

end

"""
Method to resample data in time.

# Arguments

* `data::TD` : input data of type `TD`
* `tgrid` : resampling in time according to this time grid

# Return

* data after resampling as `TD`
"""
function interp(data::TD,
		tgrid::StepRangeLen
		)
	nss = data.acqgeom.nss
	nr = data.acqgeom.nr
	dataout = TD(
	      [zeros(length(tgrid),data.acqgeom.nr[iss]) for iss=1:nss, ifield=1:length(data.fields)],
	      data.fields,tgrid,data.acqgeom)
	interp_spray!(data, dataout)
	return dataout
end

function taper!(data::TD, perc=0.0; bperc=perc,eperc=perc)
	nr = data.acqgeom.nr;	nss = data.acqgeom.nss;	nt = length(data.tgrid);
	for ifield = 1:length(data.fields), iss = 1:nss
		dd=data.d[iss, ifield]
		Utils.taper!(dd,bperc=bperc,eperc=eperc)
	end
	return data
end

"""
Method to resample data in time.
Can reduce allocations =========

# Arguments

* `data::TD` : input data of type `TD`
* `dataout::TD` : preallocated data of type `TD` that is modified
"""
function interp_spray!(data::TD, dataout::TD, attrib=:interp, Battrib=:B1; pa=nothing)
	# check if datasets are similar
	nss = data.acqgeom.nss
	nr = data.acqgeom.nr
	xin=data.tgrid
	xout=dataout.tgrid
	if(pa===nothing)
		pa=Interpolation.Kernel([xin], [xout], :B1)
	end
	for ifield = 1:length(data.fields), iss = 1:nss
		dat=data.d[iss,ifield]
		dato=dataout.d[iss,ifield]
		for ir = 1:nr[iss]
			din=view(dat,:,ir)
			dout=view(dato,:,ir)
	 		Interpolation.interp_spray!(din, dout, pa, attrib)
		end
	end
	return dataout
end

"""
Method used to preallocate `TD` with zeros.

# Arguments

* `fields::Vector{Symbol}` : number of components
* `tgrid` : time domain grid
* `acqgeom::Acquisition.Geom` : acquisition geometry

# Return

* data with zeros as `TD`
"""
function TD_zeros(fields::Vector{Symbol}, tgrid::StepRangeLen, acqgeom::Acquisition.Geom)
	return TD([zeros(length(tgrid),acqgeom.nr[iss]) for iss=1:acqgeom.nss, ifield=1:length(fields)],fields,
	   deepcopy(tgrid),deepcopy(acqgeom))
end
function Base.fill!(data::TD, k::Float64)
	for iss=1:data.acqgeom.nss, ifield=1:length(data.fields)
		fill!(data.d[iss,ifield],k)
	end
end
function TD_zeros(d::TD)
	return TD([zeros(length(d.tgrid),d.acqgeom.nr[iss]) for iss=1:d.acqgeom.nss, ifield=1:length(d.fields)],d.fields,
	   deepcopy(d.tgrid),deepcopy(d.acqgeom))
end

"Same as `TD_zeros`, except for returning ones"
function TD_ones(fields::Vector{Symbol}, tgrid::StepRangeLen, acqgeom::Acquisition.Geom)
	return TD([ones(length(tgrid),acqgeom.nr[iss]) for iss=1:acqgeom.nss, ifield=1:length(fields)],
	   fields,deepcopy(tgrid),deepcopy(acqgeom))
end

"""
Time reverse the records of each receiver in `TD`

# Arguments

* `data::TD` : input data that is modified
"""
function TD_tr!(data::TD)
	data.d = copy([reverse(data.d[i,j],dims=1) for i in 1:data.acqgeom.nss, j in 1:length(data.fields)]);
end

function addnoise!(dataN::TD, data::TD, SNR)

	σx=Statistics.var(data)

	σxN=sqrt(σx^2*inv(10^(SNR/10.)))

	# factor to be multiplied to each scalar
	α=sqrt(σxN)
	for ifield = 1:length(data.fields), iss = 1:data.acqgeom.nss
		for ir = 1:data.acqgeom.nr[iss], it = 1:length(data.tgrid)

			dataN.d[iss, ifield][it, ir] = dataN.d[iss, ifield][it, ir] + α*Random.randn()
		end
	end
end

# include rest of the files
for file in ["base", "statistics", "processing", "misfits", "weights"]
	fn=joinpath(@__DIR__, string(file,".jl"))
	include(fn)
	#Revise.track(@__MODULE__,fn)
end

#=
function DDecon(d::TD, wav::AbstractVector{Float64}, ϵ=1e-2)

	dout=TD_zeros(d)
	ntd=length(dout.tgrid)
	wavv=deepcopy(wav);

	paD=Decon.ParamD(ntd=ntd,nts=length(wav), s=wavv)

	paD.ϵ=ϵ

	DDecon!(dout, d, paD)
	return dout
end

function DDecon!(dataout::TD, data::TD, paD)
	nr = data.acqgeom.nr;	nss = data.acqgeom.nss;	nt = length(data.tgrid);
	for ifield = 1:length(data.fields), iss = 1:nss
		dd=data.d[iss, ifield]
		ddo=dataout.d[iss, ifield]
		for ir = 1:nr[iss]
			for it in 1:nt
				paD.d[it]=dd[it,ir]
			end
			Decon.mod!(paD)
			for it in 1:nt
				ddo[it,ir]=paD.g[it]
			end
		end
	end
	return dataout
end

=#

end # module


text:
module Attenuations

using AxisArrays
using HTTP
using Unitful

import Unitful: g, cm, eV, keV, MeV
export eV, keV, MeV

export μ,
    μᵨ,
    Element,
    Compound,
    Mixture,
    Elements,
    Materials,
    PhotoelectricAbsorption,
    Coherent,
    Incoherent,
    InNuclearField,
    InElectronField,
    WithCoherent,
    WithoutCoherent,
    val

abstract type Attenuation end

"""
PhotoelectricAbsorption is one of the principal forms of interaction of x-ray
and gamma photons with matter. A low energy photon interacts with an electron
in the atom and removes it from its shell.

https://en.wikipedia.org/wiki/Photoelectric_effect
https://radiopaedia.org/articles/photoelectric-effect
"""
struct PhotoelectricAbsorption <: Attenuation end

abstract type Scattering <: Attenuation end
"""
Coherent scattering (also known as unmodified, Rayleigh, classical or elastic
scattering) is one of three forms of photon interaction which occurs when the
energy of the x-ray or gamma photon is small in relation to the ionization
energy of the atom. It therefore occurs with low energy radiation.

Upon interacting with the attenuating medium, the photon does not have enough
energy to liberate the electron from its bound state (i.e. the photon energy is
well below the binding energy of the electron) so no energy transfer occurs.
There is no energy deposition and thus no dose resulting from coherent
scattering. The only change is a change of direction (scatter) of the photon,
hence 'unmodified' scatter. Coherent scattering is not a major interaction
process encountered in radiography at the energies normally used.

Coherent scattering varies with the atomic number of the absorber (Z) and
incident photon energy (E) by Z/E².

https://radiopaedia.org/articles/coherent-scattering
"""
struct Coherent <: Scattering end

"""
Incoherent scatter (Compton effect or Compton scatter) is one of principle forms
of photon interaction. It is the main cause of scattered radiation in a
material. It occurs due to the interaction of the photon (x-ray or gamma) with
free electrons (unattached to atoms) or loosely bound valence shell
(outer shell) electrons. The resultant incident photon is scattered
(changes direction) and imparts energy to the electron (recoil electron). The
scattered photon will have a different wavelength (observed phenomenon) and thus
a different energy (E=hc/λ). Energy and momentum are conserved in this process.
The Compton effect is a partial absorption process and as the original photon
has lost energy, known as Compton shift (i.e. a shift of wavelength/frequency).
The wavelength change of the scattered photon can be determined by
0.024*(1- cosθ), where θ is scattered photon angle. Thus, the energy of the
scattered photon decreases with increasing scattered photon angle.

https://en.wikipedia.org/wiki/Compton_scattering
https://radiopaedia.org/articles/compton-effect?lang=us
"""
struct Incoherent <: Scattering end

abstract type PairProduction <: Attenuation end
struct InNuclearField <: PairProduction end
struct InElectronField <: PairProduction end
abstract type Total <: Attenuation end
struct WithCoherent <: Total end
struct WithoutCoherent <: Total end

const DefaultAttenuation = WithCoherent

abstract type Matter end

μᵨ(
    m::Matter,
    energy::T,
    a::Type{A} = DefaultAttenuation,
) where {T<:Unitful.Energy,A<:Attenuation} = μᵨ(m, [energy], a)[1]

μ(
    m::Matter,
    energy::T,
    a::Type{A} = DefaultAttenuation,
) where {T<:Unitful.Energy,A<:Attenuation} = μ(m, [energy], a)[1]

val(a::AbstractArray{T}) where {T<:Unitful.AbstractQuantity} =
    [i.val for i in a.data]

val(a::T) where {T<:Unitful.AbstractQuantity} = a.val

"""
    Element

An element of the periodic table.
"""
struct Element{T,S} <: Matter where {T<:Unitful.Energy,S<:Unitful.Density}
    Z::Int
    symbol::String
    name::String
    ZAratio::Float64
    I::T
    ρ::S
end

Base.show(io::IO, e::Element) = print(
    io,
    "$(e.Z) $(e.symbol) $(e.name) Z/A=$(e.ZAratio) I=$(e.I.val)eV ρ=$(e.ρ.val)g/cm³",
)

function μᵨ(
    e::Element,
    energies::AbstractArray{<:Unitful.Energy},
    a::Type{<:Attenuation} = DefaultAttenuation,
)
    body = Dict{String,String}(
        "Method" => "1",
        "ZNum" => "$(e.Z)",
        bodykey(a) => "on",
        "Energies" => formatenergies(energies),
    )

    μᵨ = XCOM(body) * cm^2 ./ g
    AxisArray(μᵨ, Axis{:energy}(energies))
end

μ(
    e::Element,
    energies::AbstractArray{<:Unitful.Energy},
    a::Type{<:Attenuation} = DefaultAttenuation,
) = AxisArray(e.ρ * μᵨ(e, energies, a), Axis{:energy}(energies))

"""
    Compound
"""
struct Compound <: Matter
    formula::String
end

function μᵨ(
    c::Compound,
    energies::AbstractArray{<:Unitful.Energy},
    a::Type{<:Attenuation} = DefaultAttenuation,
)
    body = Dict{String,String}(
        "Method" => "2",
        "Formula" => c.formula,
        bodykey(a) => "on",
        "Energies" => formatenergies(energies),
    )

    μᵨ = XCOM(body) * cm^2 ./ g
    AxisArray(μᵨ, Axis{:energy}(energies))
end

"""
    Mixture
"""
struct Mixture{T} <: Matter where {T<:AbstractFloat}
    formulae::Dict{String,T}
end

function μᵨ(
    m::Mixture,
    energies::AbstractArray{<:Unitful.Energy},
    a::Type{<:Attenuation},
)
    body = Dict{String,String}(
        "Method" => "3",
        "Formulae" => join(["$k $v" for (k, v) in m.formulae], '\n'),
        bodykey(a) => "on",
        "Energies" => formatenergies(energies),
    )

    μᵨ = XCOM(body) * cm^2 ./ g
    AxisArray(μᵨ, Axis{:energy}(energies))
end

struct Material{T,S} <: Matter where {T<:Unitful.Energy,S<:Unitful.Density}
    name::String
    ZAratio::Float64
    I::T
    ρ::S
    composition::Dict{Int,Float64}
end

μᵨ(
    m::Material,
    energies::AbstractArray{<:Unitful.Energy},
    a::Type{<:Attenuation} = DefaultAttenuation,
) = μᵨ(
    Mixture(Dict([(Elements[k].symbol, v) for (k, v) in m.composition])),
    energies,
    a,
)

μ(
    m::Material,
    energies::AbstractArray{<:Unitful.Energy},
    a::Type{<:Attenuation} = DefaultAttenuation,
) = AxisArray(m.ρ * μᵨ(m, energies, a), Axis{:energy}(energies))

Base.show(io::IO, m::Material) = print(
    io,
    "$(m.name) Z/A=$(m.ZAratio) I=$(m.I) ρ=$(m.ρ)\r\n",
    join(["$k: $v" for (k, v) in m.composition], "\r\n"),
)

include("xcom.jl")

include("elements.jl")
include("materials.jl")

end # module

ext: jl


text:
# Example for HINT layer (Kruse et al, 2020)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using InvertibleNetworks, LinearAlgebra, Test

# Input
nx = 32
ny = 32
n_channel = 32
n_hidden = 64
batchsize = 2

# Input image
X = glorot_uniform(nx, ny, n_channel, batchsize)
Y = glorot_uniform(nx, ny, n_channel, batchsize)

# Conditional HINT layer
CH = ConditionalLayerHINT(nx, ny, n_channel, n_hidden, batchsize)

# Forward/inverse
Zx, Zy, logdet = CH.forward(X, Y)
X_, Y_ = CH.inverse(Zx, Zy)

@test isapprox(norm(X - X_)/norm(X), 0f0; atol=1f-4)
@test isapprox(norm(Y - Y_)/norm(Y), 0f0; atol=1f-4)

# Forward/inverse Y-lane only
Zy = CH.forward_Y(Y)
Y_ = CH.inverse_Y(Zy)

@test isapprox(norm(Y - Y_)/norm(Y), 0f0; atol=1f-4)


text:
using Quaternions
const Q = Quaternions

@enum ProjectionEnum Perspective Orthographic

abstract type AbstractCamera end

@reactivecomposed type BasicCamera <: AbstractCamera
    Area
    Projection
    View
    ProjectionView
end

@composed type PerspectiveCamera <: AbstractCamera
    <: BasicCamera
    ProjectionType
    Translation
    Rotation
    RotationSpeed
    TranslationSpeed
    EyePosition
    LookAt
    UpVector
    Fov
    Near
    Far
end

default(::Type{Rotation}, ::Partial{<: AbstractCamera}) = Vec3f0(0, 0, 0)

function add_translation(cam, scene, key, button)
    local last_mousepos::Vec2f0 = Vec2f0(0, 0)
    lift_node(scene[:mousedrag]) do drag
        if ispressed(scene, key) && ispressed(scene, button)
            if drag == Mouse.down
                #just started pressing, nothing to do yet
                last_mousepos = Vec2f0(to_value(scene, :mouseposition))
            elseif drag == Mouse.pressed
                 # we need the difference, although I'm wondering if absolute wouldn't be better.
                 # TODO FIND OUT! Definitely would have more precision problems
                mousepos = Vec2f0(scene[:mouseposition])
                diff = (last_mousepos - mousepos) * to_value(cam, :translationspeed)
                last_mousepos = mousepos
                cam[:translation] = Vec3f0(0f0, diff[1], - diff[2])
            end
        end
        return
    end
    lift_node(scene, Mouse.Scroll) do scroll
        if ispressed(scene, button)
            cam[:translation] = Vec3f0(scroll[2], 0f0, 0f0)
        end
        return
    end
end
function addrotation(cam, ::Type{Rotation}, canvas, key, button)
    local last_mousepos::Vec2f0 = Vec2f0(0, 0)
    on(canvas, Mouse.Drag) do drag
        if ispressed(canvas, key) && ispressed(canvas, button)
            if drag == Mouse.down
                last_mousepos = Vec2f0(canvas[Mouse.Position])
            elseif drag == Mouse.pressed
                mousepos = Vec2f0(canvas[Mouse.Position])
                mp = (last_mousepos - mousepos) * cam[RotationSpeed]
                last_mousepos = mousepos
                cam[Rotation] = Vec3f0(mp[1], -mp[2], 0f0)
            end
        end
        return
    end
end

function add!(cam, ::Type{EyePosition}, ::Type{LookAt})
    on(cam, Translation) do translation
        translation == Vec3f0(0) && return

        lookat = cam[LookAt]; eyepos = cam[EyePosition]; up = cam[UpVector]
        projview = cam[ProjectionView]; area = cam[Area]; prjt = cam[ProjectionType]

        dir = eyepos - lookat
        dir_len = norm(dir)
        cam_res = Vec2f0(widths(area))
        zoom, x, y = translation
        zoom *= 0.1f0 * dir_len

        if prjt != Perspective
            x, y = to_worldspace(Vec2f0(x, y), projview, cam_res)
        else
            x, y = (Vec2f0(x, y) ./ cam_res) .* dir_len
        end
        dir_norm = normalize(dir)
        right = normalize(cross(dir_norm, up))
        zoom_trans = dir_norm * zoom

        side_trans = right * (-x) + normalize(up) * y
        newpos = eyepos + side_trans + zoom_trans
        cam[EyePosition] = newpos
        cam[LookAt] = lookat + side_trans
        return
    end
    on(cam, Rotation) do theta_v
        theta_v == Vec3f0(0) && return #nothing to do!
        eyepos_v = cam[EyePosition]; lookat_v = cam[LookAt]; up_v = cam[UpVector]

        dir = normalize(eyepos_v - lookat_v)
        right_v = normalize(cross(up_v, dir))
        up_v  = normalize(cross(dir, right_v))
        rotation = rotate_cam(theta_v, right_v, Vec3f0(0, 0, 1), dir)
        r_eyepos = lookat_v + rotation * (eyepos_v - lookat_v)
        r_up = normalize(rotation * up_v)
        cam[EyePosition] = r_eyepos
        cam[UpVector] = r_up
        return
    end
end

function add!(cam::PerspectiveCamera, ::Type{ProjectionView})
    on(cam, Area, Fov, Near, ProjectionType, LookAt, EyePosition, UpVector) do area, fov, near, projectiontype, lookatv, eyeposition, upvector
        zoom = norm(lookatv - eyeposition)
        # TODO this means you can't set FarClip... SAD!
        far = max(zoom * 5f0, 30f0)
        proj = projection_switch(area, fov, near, far, projectiontype, zoom)
        view = lookat(eyeposition, lookatv, upvector)
        cam[Projection] = proj
        cam[View] = view
        cam[ProjectionView] = proj * view
    end
end

function projection_switch{T <: Real}(
        wh::SimpleRectangle,
        fov::T, near::T, far::T,
        projection::ProjectionEnum, zoom::T
    )
    aspect = T(wh.w / wh.h)
    h = T(tan(fov / 360.0 * pi) * near)
    w = T(h * aspect)
    projection == Perspective && return frustum(-w, w, -h, h, near, far)
    h, w = h * zoom, w * zoom
    orthographicprojection(-w, w, -h, h, near, far)
end

function rotate_cam{T}(
        theta::Vec{3, T},
        cam_right::Vec{3, T}, cam_up::Vec{3, T}, cam_dir::Vec{3, T}
    )
    rotation = one(Q.Quaternion{T})
    # first the rotation around up axis, since the other rotation should be relative to that rotation
    if theta[1] != 0
        rotation *= Q.qrotation(cam_up, theta[1])
    end
    # then right rotation
    if theta[2] != 0
        rotation *= Q.qrotation(cam_right, theta[2])
    end
    # last rotation around camera axis
    if theta[3] != 0
        rotation *= Q.qrotation(cam_dir, theta[3])
    end
    rotation
end

@defaults function camera3d(backend, scene, kw_args)
    rotationspeed = to_float(1)
    translationspeed = to_float(1)
    eyeposition = Vec3f0(3)
    lookat = Vec3f0(0)
    upvector = Vec3f0((0, 0, 1))
    fov = to_float(45f0)
    near = to_float(0.01f0)
    far = to_float(100f0)
    projectiontype = to_projection(Perspective)
end


text:
using PkgBenchmark, CSV, WeakRefStrings
if !is_windows()
    using DecFP
end

prep(io::IO, ::Type{Int64}) = write(io, "10")
prep(io::IO, ::Type{Float64}) = write(io, "10.0")
prep(io::IO, ::Type{WeakRefString{UInt8}}) = write(io, "hey there sailor")
prep(io::IO, ::Type{String}) = write(io, "hey there sailor")
prep(io::IO, ::Type{Date}) = write(io, "2016-09-28")
prep(io::IO, ::Type{DateTime}) = write(io, "2016-09-28T03:21:00")
if !is_windows()
    prep(io::IO, ::Type{Dec64}) = write(io, "10.0")
end

function prep{T}(::Type{IOBuffer}, ::Type{T})
    io = IOBuffer()
    prep(io, T)
    seekstart(io)
    return io, ()->return
end
function prep{T}(::Type{IOStream}, ::Type{T})
    t = tempname()
    io = open(t, "w")
    prep(io, T)
    close(io)
    io = open(t, "r")
    return io, ()->rm(t)
end

TYPES = !is_windows() ? (Int, Float64, WeakRefString{UInt8}, String, Date, DateTime, Dec64) : (Int, Float64, WeakRefString{UInt8}, String, Date, DateTime)

@benchgroup "CSV" begin
    @benchgroup "CSV.parsefield" begin
        opts = CSV.Options()
        row = col = 1
        state = Ref{CSV.ParsingState}(CSV.None)
        for I in (IOBuffer, IOStream)
            for T in TYPES
                io, f = prep(I, T)
                @bench "$I - $T" CSV.parsefield($io, $T, opts, row, col, state)
                f()
            end
        end
    end
    FILE = joinpath(dirname(@__FILE__), "randoms_small.csv")

    @benchgroup "CSV.read" begin
        @bench "CSV.read" CSV.read(FILE)
    end

    @benchgroup "CSV.write" begin
        df = CSV.read(FILE)
        t = tempname()
        @bench "CSV.write" CSV.write(t, df)
    end

end

# generate single column files w/ 1M rows for each type
using WeakRefStrings

val = "hey"
for i in (1001, 100.1, WeakRefString{UInt8}(pointer(val), 3, 0), Date(2008, 1, 3), DateTime(2008, 3, 4))
    open("/Users/jacobquinn/Downloads/randoms_$(typeof(i)).csv", "w") do f
        for j = 1:1_000_000
            write(f, string(i))
            write(f, "\n")
        end
    end
end

using CSV, TextParse
for T in (Int, Float64, WeakRefStrings.WeakRefString{UInt8}, Date, DateTime)
    println("comparing for T = $T...")
    # T == WeakRefStrings.WeakRefString{UInt8} && continue
    @time CSV.read("/Users/jacobquinn/Downloads/randoms_$(T).csv");
    # @time TextParse.csvread("/Users/jacobquinn/Downloads/randoms_$T.csv");
end

for T in (Int, Float64, WeakRefStrings.WeakRefString{UInt8}, Date, DateTime)
    println("comparing for T = $T...")
    # T == WeakRefStrings.WeakRefString{UInt8} && continue
    # @time CSV.read("/Users/jacobquinn/Downloads/randoms_$(T).csv");
    @time TextParse.csvread("/Users/jacobquinn/Downloads/randoms_$T.csv");
end

@time CSV.read("/Users/jacobquinn/Downloads/yellow_tripdata_2015-01.csv");
@time TextParse.csvread("/Users/jacobquinn/Downloads/yellow_tripdata_2015-01.csv");

for T in ('Int64', 'Float64', 'WeakRefString{UInt8}', 'Date', 'DateTime'):
    start = time.time()
    delim = ','
    table = pandas.read_csv("/Users/jacobquinn/Downloads/randoms_" + T + ".csv", delimiter=delim)
    end = time.time()
    print(end - start)

start = time.time()
delim = ','
table = pandas.read_csv("/Users/jacobquinn/Downloads/yellow_tripdata_2015-01.csv", delimiter=delim)
end = time.time()
print(end - start)
@time df = CSV.read("/Users/jacobquinn/Downloads/file.txt"; delim=' ');
@time TextParse.csvread("/Users/jacobquinn/Downloads/randoms_$(T).csv")
# julia> for T in (Int, Float64, WeakRefStrings.WeakRefString{UInt8}, Date, DateTime)
#            println("comparing for T = $T...")
#            @time CSV.read("/Users/jacobquinn/Downloads/randoms_$(T).csv");
#            @time TextParse.csvread("/Users/jacobquinn/Downloads/randoms_$(T).csv");
#        end
# comparing for T = Int64...
# pre-allocating DataFrame w/ rows = 999999
#   0.043684 seconds (1.00 M allocations: 22.929 MiB, 31.61% gc time)
#   0.045556 seconds (460 allocations: 15.575 MiB, 3.20% gc time)
# comparing for T = Float64...
# pre-allocating DataFrame w/ rows = 999999
#   0.080026 seconds (1.00 M allocations: 22.974 MiB, 23.80% gc time)
#   0.082530 seconds (457 allocations: 16.528 MiB)
# comparing for T = WeakRefString{UInt8}...
# pre-allocating DataFrame w/ rows = 999999
#   0.058446 seconds (1.89 k allocations: 22.986 MiB, 8.53% gc time)
#   0.069034 seconds (595 allocations: 5.188 MiB)
# comparing for T = Date...
# pre-allocating DataFrame w/ rows = 999999
#   0.125229 seconds (2.00 M allocations: 53.504 MiB, 20.94% gc time)
#   0.120472 seconds (1.00 M allocations: 51.846 MiB, 6.73% gc time)
# comparing for T = DateTime...
# pre-allocating DataFrame w/ rows = 999999
#   0.175855 seconds (2.00 M allocations: 53.504 MiB, 23.30% gc time)
#   0.187619 seconds (1.00 M allocations: 60.516 MiB, 4.40% gc time)

T = Int64
@time source = CSV.Source("/Users/jacobquinn/Downloads/randoms_$(T).csv";)
@time source = CSV.Source("/Users/jacobquinn/Downloads/randoms_small.csv"; allowmissing=:auto)
@time source = CSV.Source("/Users/jacobquinn/Downloads/randoms_small.csv"; allowmissing=:none)
# source.schema = DataStreams.Data.Schema(DataStreams.Data.header(source.schema), (Int, String, String, Float64, Float64, Date, DateTime), 9)
# @time df = CSV.read(source, NamedTuple);
sink = Si = NamedTuple
transforms = Dict{Int,Function}(1=>x->x-1)
append = false
args = kwargs = ()
source_schema = DataStreams.Data.schema(source)
sink_schema, transforms2 = DataStreams.Data.transform(source_schema, transforms, true);
sinkstreamtype = DataStreams.Data.Field
sink = Si(sink_schema, sinkstreamtype, append, args...; kwargs...);
columns = []
filter = x->true
@code_warntype DataStreams.Data.stream!(source, sinkstreamtype, sink, source_schema, sink_schema, transforms2, filter, columns)
@time DataStreams.Data.stream!(source, sinkstreamtype, sink, source_schema, sink_schema, transforms2, filter, columns)

function testt(t)
    a = getfield(t, 1)
    b = getfield(t, 2)
    c = getfield(t, 3)
    d = getfield(t, 4)
    e = getfield(t, 5)
    f = getfield(t, 6)
    g = getfield(t, 7)
    return (a, b, c, d, e, f, g)
end
@code_warntype testt((i1=(?Int)[], i2=(?String)[], i3=(?String)[], i4=(?Float64)[], i5=(?Float64)[], i6=(?Date)[], i7=(?DateTime)[]))

@code_llvm DataStreams.Data.stream!(source, sinkstreamtype, sink, source_schema, sink_schema, transforms2, filter, columns)
@time DataStreams.Data.stream!(source, sinkstreamtype, sink, source_schema, sink_schema, transforms2, filter, columns)

@code_warntype @time CSV.parsefield(IOBuffer(), ?Int, CSV.Options(), 0, 0, CSV.STATE)

t = Vector{Int}(1000000)

# having CSV.parsefield(io, T) where T !>: Missing decreases allocations by 1.00M
# inlining CSV.parsefield also dropped allocations
# making CSV.Options not have a type parameter also sped things up
#

using BenchmarkTools

g(x) = x < 5 ? x : -1
A = [i for i = 1:10]
function get_then_set(A)
    @simd for i = 1:10
        @inbounds A[i] = g(i)
    end
    return A
end
@code_warntype g(1)
@code_warntype get_then_set(A)
@benchmark get_then_set(A) # 20ns

@inline g3(x) = g2(x)
@inline function g2(x)
    if x < 20
        return x * 20
    end

    if x < 15
        return nothing
    end

    if x < 12
        return 2x
    end

    if x * 20 / 4 % 2 == 0
        return 1
    end

    if x < 0
        return nothing
    end
    return nothing
end

A = Union{Int, Void}[i for i = 1:10]
@inline function get_then_set2(A)
    @simd for i = 1:10
        # Base.arrayset(A, g2(i), i)
        val = g3(i)
        if val isa Void
            @inbounds A[i] = val#::Union{Int, Void}
        else
            @inbounds A[i] = val#::Union{Int, Void}
        end
    end
    return A
end
function run_lots(N)
    A = Union{Int, Void}[i for i = 1:10]
    for i = 1:N
        get_then_set2(A)
    end
    return
end

@code_warntype g2(1)
@code_warntype get_then_set2(A)
@code_llvm get_then_set2(A)
@benchmark get_then_set2(A) # 155ns

g4(x::Int) = 1
g4(x::Void) = 0

A = [i for i = 1:10]
function get_sum(A)
    s = 0
    for a in A
        s += g4(a)
    end
    return s
end
@code_warntype get_sum(A)
@code_llvm get_sum(A)
@benchmark get_sum(A) # 24ns

A = Union{Int, Void}[i for i = 1:10]
A[[3, 5, 7]] = nothing
function get_sum2(A)
    s = 0
    for a in A
        s += g4(a)
    end
    return s
end
@code_warntype get_sum2(A)
@code_llvm get_sum(A)
@benchmark get_sum2(A) # 100ns

function getstatic{T}(t::T)
    return t[1]
end


text:
function update(b::HyperRectangle{N, T}, v::Vec{N, T2}) where {N, T, T2}
    update(b, Vec{N, T}(v))
end
function update(b::HyperRectangle{N, T}, v::Vec{N, T}) where {N, T}
    m = min.(minimum(b), v)
    maxi = maximum(b)
    mm = if isnan(maxi)
        v-m
    else
        max.(v, maxi) - m
    end
    HyperRectangle{N, T}(m, mm)
end

# Min maximum distance functions between hrectangle and point for a given dimension
@inline function min_dist_dim(rect::HyperRectangle{N, T}, p::Vec{N, T}, dim::Int) where {N, T}
    max(zero(T), max(minimum(rect)[dim] - p[dim], p[dim] - maximum(rect)[dim]))
end

@inline function max_dist_dim(rect::HyperRectangle{N, T}, p::Vec{N, T}, dim::Int) where {N, T}
    max(maximum(rect)[dim] - p[dim], p[dim] - minimum(rect)[dim])
end

@inline function min_dist_dim(rect1::HyperRectangle{N, T},
                              rect2::HyperRectangle{N, T},
                              dim::Int) where {N, T}
    max(zero(T), max(
        minimum(rect1)[dim] - maximum(rect2)[dim],
        minimum(rect2)[dim] - maximum(rect1)[dim]
    ))
end

@inline function max_dist_dim(rect1::HyperRectangle{N, T},
                              rect2::HyperRectangle{N, T},
                              dim::Int) where {N, T}
    max(
        maximum(rect1)[dim] - minimum(rect2)[dim],
        maximum(rect2)[dim] - minimum(rect1)[dim]
    )
end

# Total minimum maximum distance functions
@inline function min_euclideansq(rect::HyperRectangle{N, T},
                                 p::Union{Vec{N, T},
                                 HyperRectangle{N, T}}) where {N, T}
    minimum_dist = T(0.0)
    for dim in 1:length(p)
        d = min_dist_dim(rect, p, dim)
        minimum_dist += d*d
    end
    return minimum_dist
end

# could add an @symmetric macro, which defines min_euclidean(pt, s) from
# min_euclidean(s, pt) automatically etc.
@inline min_euclidean(pt1::Vec, pt2::Vec) = norm(pt1-pt2)
min_euclidean(pt::Vec, s::Simplex) = √(sqdist(pt, s))
min_euclidean(s::Simplex, pt::Vec) = min_euclidean(pt, s)
min_euclidean(r1::HyperRectangle, r2::Vec) = sqrt(min_euclideansq(r1, r2))
min_euclidean(r1::Vec, r2::HyperRectangle) = sqrt(min_euclideansq(r1, r2))
min_euclidean(r1::HyperRectangle, r2::HyperRectangle) = sqrt(min_euclideansq(r1, r2))
min_euclidean(c1::AbstractConvexHull, pt::Vec) = gjk(c1,pt)
min_euclidean(pt::Vec, c2::AbstractConvexHull) = gjk(pt,c2)
min_euclidean(c1::AbstractConvexHull, c2::AbstractConvexHull) = gjk(c1,c2)

@inline function max_euclideansq(rect::HyperRectangle{N, T}, p::Union{Vec{N, T}, HyperRectangle{N, T}}) where {N, T}
    maximum_dist = T(0.0)
    for dim in 1:length(p)
        d = max_dist_dim(rect, p, dim)
        maximum_dist += d*d
    end
    return maximum_dist
end
function max_euclidean(rect::HyperRectangle{N, T}, p::Union{Vec{N, T}, HyperRectangle{N, T}}) where {N, T}
    sqrt(max_euclideansq(rect, p))
end

# Functions that return both minimum and maximum for convenience
@inline function minmax_dist_dim(rect::HyperRectangle{N, T}, p::Union{Vec{N, T}, HyperRectangle{N, T}}, dim::Int) where {N, T}
    minimum_d = min_dist_dim(rect, p, dim)
    maximum_d = max_dist_dim(rect, p, dim)
    return minimum_d, maximum_d
end

@inline function minmax_euclideansq(rect::HyperRectangle{N, T}, p::Union{Vec{N, T}, HyperRectangle{N, T}}) where {N, T}
    minimum_dist = min_euclideansq(rect, p)
    maximum_dist = max_euclideansq(rect, p)
    return minimum_dist, maximum_dist
end

function minmax_euclidean(rect::HyperRectangle{N, T}, p::Union{Vec{N, T}, HyperRectangle{N, T}}) where {N, T}
    minimumsq, maximumsq = minmax_euclideansq(rect, p)
    return sqrt(minimumsq), sqrt(maximumsq)
end


text:
@testset "Storage data misspecification" begin
    # See https://discourse.julialang.org/t/how-to-use-test-warn/15557/5 about testing for warning throwing
    warn_message = "The data doesn't include devices of type GenericBattery, consider changing the device models"
    model = DeviceModel(GenericBattery, BookKeeping)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    op_problem = OperationsProblem(MockOperationProblem, DCPPowerModel, c_sys5)
    @test_logs (:info,) (:warn, warn_message) match_mode = :any mock_construct_device!(
        op_problem,
        model,
    )
    op_problem = OperationsProblem(MockOperationProblem, DCPPowerModel, c_sys14)
    @test_logs (:info,) (:warn, warn_message) match_mode = :any mock_construct_device!(
        op_problem,
        model,
    )
end

@testset "Storage Basic Storage With DC - PF" begin
    model = DeviceModel(GenericBattery, BookKeeping)
    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat")
    op_problem = OperationsProblem(MockOperationProblem, DCPPowerModel, c_sys5_bat)
    mock_construct_device!(op_problem, model)
    moi_tests(op_problem, false, 72, 0, 72, 72, 24, false)
    psi_checkobjfun_test(op_problem, GAEVF)
end

@testset "Storage Basic Storage With AC - PF" begin
    model = DeviceModel(GenericBattery, BookKeeping)
    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat")
    op_problem = OperationsProblem(MockOperationProblem, ACPPowerModel, c_sys5_bat)
    mock_construct_device!(op_problem, model)
    moi_tests(op_problem, false, 96, 0, 96, 96, 24, false)
    psi_checkobjfun_test(op_problem, GAEVF)
end

@testset "Storage with Reservation DC - PF" begin
    model = DeviceModel(GenericBattery, BookKeepingwReservation)
    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat")
    op_problem = OperationsProblem(MockOperationProblem, DCPPowerModel, c_sys5_bat)
    mock_construct_device!(op_problem, model)
    moi_tests(op_problem, false, 96, 0, 72, 72, 24, true)
    psi_checkobjfun_test(op_problem, GAEVF)
end

@testset "Storage with Reservation With AC - PF" begin
    model = DeviceModel(GenericBattery, BookKeepingwReservation)
    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat")
    op_problem = OperationsProblem(MockOperationProblem, ACPPowerModel, c_sys5_bat)
    mock_construct_device!(op_problem, model)
    moi_tests(op_problem, false, 120, 0, 96, 96, 24, true)
    psi_checkobjfun_test(op_problem, GAEVF)
end

@testset "BatteryEMS with EndOfPeriodEnergyTarget with DC - PF" begin
    model = DeviceModel(BatteryEMS, EndOfPeriodEnergyTarget)
    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat_ems")
    op_problem = OperationsProblem(MockOperationProblem, DCPPowerModel, c_sys5_bat)
    mock_construct_device!(op_problem, model)
    moi_tests(op_problem, false, 72, 0, 72, 72, 24, false)
    psi_checkobjfun_test(op_problem, GAEVF)
end

@testset "BatteryEMS with EndOfPeriodEnergyTarget With AC - PF" begin
    model = DeviceModel(BatteryEMS, EndOfPeriodEnergyTarget)
    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat_ems")
    op_problem = OperationsProblem(MockOperationProblem, ACPPowerModel, c_sys5_bat)
    mock_construct_device!(op_problem, model)
    moi_tests(op_problem, false, 96, 0, 96, 96, 24, false)
    psi_checkobjfun_test(op_problem, GAEVF)
end


text:
using SpeedyWeather
using Test

include("utility_functions.jl")
include("spectral_transform.jl")
include("diffusion.jl")
include("time_stepping.jl")
include("initialize_from_rest.jl")


text:
using Test
using JAC, ..Defaults, ..TestFrames

@testset "Name" begin
    ##x global short = true
    printstyled("\nPerform tests on the JAC program; this may take a while .... \n", color=:cyan)
    ## Defaults.Constants.define("print test: open", pwd() * "/runtests.report")

    ##x include("inc-halfintegers.jl")
    ##x include("inc-wignersymbols.jl")
    ##x include("inc-racahsum.jl")

    ##x redirect_stdout(streamDummy) do
    @testset "JAC methods" begin
        @test TestFrames.testMethod_Wigner_3j()
    end

    @testset "JAC evaluations" begin
        @test TestFrames.testEvaluation_Wigner_3j_specialValues()
        @test TestFrames.testEvaluation_Wigner_6j_specialValues()
        @test TestFrames.testEvaluation_Wigner_9j_specialValues()
        @test TestFrames.testEvaluation_sumRulesForOneWnj()
        @test TestFrames.testEvaluation_sumRulesForTwoWnj()
        @test RacahAlgebra.testSpecialValuesW3j()
        @test RacahAlgebra.testSpecialValuesW6j()
        @test RacahAlgebra.testSpecialValuesW9j()
        @test RacahAlgebra.testSumRules()
    end

    @testset "JAC representations" begin
        @test TestFrames.testRepresentation_MeanFieldBasis_CiExpansion()
        @test TestFrames.testRepresentation_RasExpansion()
        @test TestFrames.testRepresentation_GreenExpansion()
    end

    @testset "JAC amplitudes" begin
        @test TestFrames.testModule_MultipoleMoment()
        ## @test TestFrames.testModule_ParityNonConservation()
    end
    ##x close(streamDummy)

    @testset "JAC properties" begin
        @test TestFrames.testModule_Einstein()
        ## @test TestFrames.testModule_Hfs()
        @test TestFrames.testModule_LandeZeeman()
        @test TestFrames.testModule_IsotopeShift()
        @test TestFrames.testModule_AlphaVariation()
        @test TestFrames.testModule_FormFactor()
        @test TestFrames.testModule_DecayYield()
        @test TestFrames.testModule_MultipolePolarizibility()
        @test TestFrames.testModule_PlasmaShift()
    end

    @testset "JAC processes" begin
        @test TestFrames.testModule_PhotoEmission()
        @test TestFrames.testModule_PhotoExcitation()
        @test TestFrames.testModule_PhotoIonization()
        @test TestFrames.testModule_PhotoRecombination()
        @test TestFrames.testModule_AutoIonization()
        @test TestFrames.testModule_Dielectronic()
        ## @test TestFrames.testModule_PhotoExcitationFluores()
        ## @test TestFrames.testModule_PhotoExcitationAutoion()
        ## @test TestFrames.testModule_RayleighCompton()
        ## @test TestFrames.testModule_MultiPhotonDeExcitation()
        ## @test TestFrames.testModule_CoulombExcitation()
    end

    @testset "JAC cascades" begin
        @test TestFrames.testModule_Cascade_StepwiseDecay()
        @test TestFrames.testModule_Cascade_PhotonIonization()
        @test TestFrames.testModule_Cascade_PhotonExcitation()
        ## @test TestFrames.testModule_Cascade_Simulation()
    end

    ## Defaults.setDefaults("print test: close", "")
end


text:
module LincombTests

using Test
using Gridap.Arrays
using Gridap.Fields
using Gridap.Fields: MockField, MockBasis
using FillArrays
using Gridap.TensorValues

np = 4
p = Point(1,2)
x = fill(p,np)
v = 2.0
d = 2
ndof = 8
wi = 3.0
w = fill(wi,ndof)
f = MockBasis{d}(v,ndof)
g = lincomb(f,w)
fx = evaluate(f,x)
∇fx = evaluate(∇(f),x)
gx = fx*w
∇gx = ∇fx*w
test_field(g,x,gx,grad=∇gx)

l = 10
af = Fill(f,l)
ax = fill(x,l)
aw = fill(w,l)
ag = lincomb(af,aw)
agx = fill(gx,l)
a∇gx = fill(∇gx,l)
test_array_of_fields(ag,ax,agx,grad=a∇gx)

l = 0
af = Fill(f,l)
ax = fill(x,l)
aw = fill(w,l)
ag = lincomb(af,aw)
agx = fill(gx,l)
a∇gx = fill(∇gx,l)
test_array_of_fields(ag,ax,agx,grad=a∇gx)

end #module


text:
import LBFGSB

struct LBFGSB_Result
    # TODO: get rid of this data structure, once we move everything the
    # info-hook needs into wrk
    f::Float64
    x::Vector{Float64}
    x_previous::Vector{Float64}
    message_in::String
    message_out::String
end

function run_optimizer(optimizer::LBFGSB.L_BFGS_B, wrk, fg!, info_hook, check_convergence!)

    m = 10
    factr = 1e7
    pgtol = 1e-5
    iprint = -1 # TODO: set this based on the value of some verbosity flag
    x = wrk.pulsevals
    x_previous = copy(x)
    n = length(x)
    obj = optimizer
    bounds = zeros(3, n) # TODO: allow the problem to specify constraints
    f = 0.0
    message = "n/a"
    # clean up
    fill!(obj.task, Cuchar(' '))
    fill!(obj.csave, Cuchar(' '))
    fill!(obj.lsave, zero(Cint))
    fill!(obj.isave, zero(Cint))
    fill!(obj.dsave, zero(Cdouble))
    fill!(obj.wa, zero(Cdouble))
    fill!(obj.iwa, zero(Cint))
    fill!(obj.g, zero(Cdouble))
    fill!(obj.nbd, zero(Cint))
    fill!(obj.l, zero(Cdouble))
    fill!(obj.u, zero(Cdouble))
    # set bounds
    for i = 1:n
        obj.nbd[i] = bounds[1, i]
        obj.l[i] = bounds[2, i]
        obj.u[i] = bounds[3, i]
    end
    # start
    obj.task[1:5] = b"START"
    message_in = "START"
    message_out = ""
    # iprint = 100  # enable for L-BFGS-B internal trace printing

    while true # task loop

        message_in = strip(String(copy(obj.task)))
        # println("- start of task loop: $message_in")
        # println("Calling setulb with f=$f, |g|=$(norm(obj.g)), |x|=$(norm(x))")
        LBFGSB.setulb(
            n,
            m,
            x,
            obj.l,
            obj.u,
            obj.nbd,
            f,
            obj.g,
            factr,
            pgtol,
            obj.wa,
            obj.iwa,
            obj.task,
            iprint,
            obj.csave,
            obj.lsave,
            obj.isave,
            obj.dsave
        )
        message_out = strip(String(copy(obj.task)))
        # println("  task -> $message_out")

        if obj.task[1:2] == b"FG" # FG_LNSRCH or FG_START
            f = fg!(f, obj.g, x)
            # println("calling fg! for |x|=$(norm(x)) -> f=$f, |g|=$(norm(obj.g))")
            if (message_in == "START") || (message_in == "NEW_X")
                # new iteration
                n0 = obj.isave[13]
                wrk.searchdirection .= obj.wa[n0:n0+n-1]
            end
            if obj.task[1:5] == b"FG_ST" # FG_START
                # x is the guess for the 0 iteration
                copyto!(wrk.gradient, obj.g)
                iter_res = LBFGSB_Result(f, x, x, message_in, message_out)
                update_result!(wrk, iter_res, obj, 0)
                wrk.result.J_T = f # TODO: don't mutate result outside of update_result!
                #update_hook!(...) # TODO
                info_tuple = info_hook(wrk, iter_res, obj, 0)
                wrk.fg_count .= 0
                (info_tuple !== nothing) && push!(wrk.result.records, info_tuple)
            end
        elseif obj.task[1:5] == b"NEW_X"
            # x is the optimized pulses for the current iteration
            iter = wrk.result.iter_start + obj.isave[30]
            iter_res = LBFGSB_Result(f, x, x_previous, message_in, message_out)
            update_result!(wrk, iter_res, obj, iter)
            #update_hook!(...) # TODO
            info_tuple = info_hook(wrk, iter_res, obj, wrk.result.iter)
            wrk.fg_count .= 0
            (info_tuple !== nothing) && push!(wrk.result.records, info_tuple)
            check_convergence!(wrk.result)
            if wrk.result.converged
                fill!(obj.task, Cuchar(' '))
                obj.task[1:24] = b"STOP: NEW_X -> CONVERGED"
            else # prepare for next iteration
                copyto!(x_previous, x)
                copyto!(wrk.gradient, obj.g)
            end
        else
            if wrk.result.message == "in progress"
                wrk.result.message = String(copy(obj.task))
            end
            fill!(obj.task, Cuchar(' '))
            obj.task[1:10] = b"STOP: DONE"
            # print_lbfgsb_trace(wrk, obj, message_in, message_out)  # enable for trace-debugging
            break
        end
        # print_lbfgsb_trace(wrk, obj, message_in, message_out)  # enable for trace-debugging

    end # task loop

    return LBFGSB_Result(f, x, x_previous, message_in, message_out)

end

function update_result!(
    wrk::GrapeWrk,
    iter_res::LBFGSB_Result,
    optimizer::LBFGSB.L_BFGS_B,
    i::Int64
)
    # TODO: make this depend only on wrk. Should not be backend-dependent
    res = wrk.result
    for (k, Ψ) in enumerate(wrk.fw_states)
        copyto!(res.states[k], Ψ)
    end
    res.J_T_prev = res.J_T
    res.J_T = iter_res.f
    (i > 0) && (res.iter = i)
    if i >= res.iter_stop
        res.converged = true
        res.message = "Reached maximum number of iterations"
        # Note: other convergence checks are done in user-supplied
        # check_convergence routine
    end
    prev_time = res.end_local_time
    res.end_local_time = now()
    res.secs = Dates.toms(res.end_local_time - prev_time) / 1000.0
end

function finalize_result!(wrk::GrapeWrk, optim_res::LBFGSB_Result)
    L = length(wrk.controls)
    res = wrk.result
    res.end_local_time = now()
    ϵ_opt = reshape(optim_res.x, L, :)
    for l = 1:L
        res.optimized_controls[l] = discretize(ϵ_opt[l, :], res.tlist)
    end
    res.optim_res = optim_res
end

function print_lbfgsb_trace(
    wrk,
    optimizer::LBFGSB.L_BFGS_B,
    message_in::AbstractString,
    message_out::AbstractString;
    show_details=true
)
    n = length(wrk.pulsevals)
    println("- end of task loop: $message_in -> $message_out")
    #! format: off
    if occursin("NEW_X", message_out) && show_details
        println("   lsave[1] = $(optimizer.lsave[1]):\t initial x has been replaced by its projection in the feasible set?")
        println("   lsave[2] = $(optimizer.lsave[2]):\t problem is constrained?")
        println("   lsave[3] = $(optimizer.lsave[3]):\t every variable is upper and lower bounds?")
        println("   isave[22] = $(optimizer.isave[22]):\t total number of intervals explored in the search of Cauchy points")
        println("   isave[26] = $(optimizer.isave[26]):\t the total number of skipped BFGS updates before the current iteration")
        println("   isave[30] = $(optimizer.isave[30]):\t the number of current iteration")
        println("   isave[31] = $(optimizer.isave[31]):\t the total number of BFGS updates prior the current iteration")
        println("   isave[33] = $(optimizer.isave[33]):\t the number of intervals explored in the search of Cauchy point in the current iteration")
        println("   isave[34] = $(optimizer.isave[34]):\t the total number of function and gradient evaluations")
        println("   isave[36] = $(optimizer.isave[36]):\t the number of function value or gradient evaluations in the current iteration")
        println("   isave[37] = $(optimizer.isave[37]):\t subspace argmin is (0) within / (1) beyond the box")
        println("   isave[38] = $(optimizer.isave[38]):\t the number of free variables in the current iteration")
        println("   isave[39] = $(optimizer.isave[39]):\t the number of active constraints in the current iteration")
        println("   isave[40] = $(optimizer.isave[40]):\t n+1-isave[40] = $(n+1-optimizer.isave[40]) = the number of variables leaving the set of active constraints in the current iteration")
        println("   isave[41] = $(optimizer.isave[41]):\t the number of variables entering the set of active constraints in the current iteration")
        println("   dsave[01] = $(optimizer.dsave[1]):\t current θ in the BFGS matrix")
        println("   dsave[02] = $(optimizer.dsave[2]):\t f(x) in the previous iteration")
        println("   dsave[03] = $(optimizer.dsave[3]):\t factr*epsmch")
        println("   dsave[04] = $(optimizer.dsave[4]):\t 2-norm of the line search direction vector")
        println("   dsave[05] = $(optimizer.dsave[5]):\t the machine precision epsmch generated by the code")
        println("   dsave[07] = $(optimizer.dsave[7]):\t the accumulated time spent on searching for Cauchy points")
        println("   dsave[08] = $(optimizer.dsave[8]):\t the accumulated time spent on subspace minimization")
        println("   dsave[09] = $(optimizer.dsave[9]):\t the accumulated time spent on line search")
        println("   dsave[11] = $(optimizer.dsave[11]):\t the slope of the line search function at the current point of line search")
        println("   dsave[12] = $(optimizer.dsave[12]):\t the maximum relative step length imposed in line search")
        println("   dsave[13] = $(optimizer.dsave[13]):\t the infinity norm of the projected gradient")
        println("   dsave[14] = $(optimizer.dsave[14]):\t the relative step length in the line search")
        println("   dsave[15] = $(optimizer.dsave[15]):\t the slope of the line search function at the starting point of the line search")
        println("   dsave[16] = $(optimizer.dsave[16]):\t the square of the 2-norm of the line search direction vector")
    end
    #! format: on
end

function print_table(
    wrk,
    iter_res::LBFGSB_Result,
    optimizer::LBFGSB.L_BFGS_B,
    iteration,
    args...
)
    # TODO: make this depend only on wrk. Should not be backend-dependent
    J_T = wrk.result.J_T
    ΔJ_T = J_T - wrk.result.J_T_prev
    secs = wrk.result.secs

    iter_stop = "$(get(wrk.kwargs, :iter_stop, 5000))"
    widths = [max(length("$iter_stop"), 6), 11, 11, 11, 8, 8]

    if iteration == 0
        header = ["iter.", "J_T", "|∇J_T|", "ΔJ_T", "FG(F)", "secs"]
        for (header, w) in zip(header, widths)
            print(lpad(header, w))
        end
        print("\n")
    end

    strs = (
        "$iteration",
        @sprintf("%.2e", J_T),
        @sprintf("%.2e", norm(wrk.gradient)),
        (iteration > 0) ? @sprintf("%.2e", ΔJ_T) : "n/a",
        @sprintf("%d(%d)", wrk.fg_count[1], wrk.fg_count[2]),
        @sprintf("%.1f", secs),
    )
    for (str, w) in zip(strs, widths)
        print(lpad(str, w))
    end
    print("\n")
    flush(stdout)
end


text:
module introspection
using InteractiveUtils

include("__text_utils.jl")

import Term: highlight, theme

import ..consoles: console
import ..panel: Panel, TextBox
import ..layout: Spacer, hLine
include("_inspect.jl")

export inspect

# ---------------------------------------------------------------------------- #
#                                   TYPEINFO                                   #
# ---------------------------------------------------------------------------- #

"""
    TypeInfo

Stores metadata about a DataType
"""
struct TypeInfo
    name::String
    supertypes::Union{Nothing,Tuple}
    subtypes::Union{Nothing,Vector}
    fields::Union{Nothing,Dict}
    constructors::Vector
    methods::Vector  # functions using the target type
    docs::Union{Nothing,Vector{Docs.DocStr}}
    docstring::String
end

"""
    TypeInfo(type::DataType)

Extract information from a DataType and store it as a `TypeInfo` object.
"""
function TypeInfo(type::DataType)
    # get super/sub types
    super = length(supertypes(type)) > 0 ? (supertypes(type)) : nothing
    sub = length(subtypes(type)) > 0 ? subtypes(type) : nothing

    # get docstring
    doc, docstring = get_docstring(Symbol(type))

    # get fields
    if !isabstracttype(type) && length(fieldnames(type)) > 0
        fields = Dict("names" => fieldnames(type), "types" => fieldtypes(type))
    else
        fields = nothing
    end

    # get constructors and methods
    _constructors = split_lines(string(methods(type)))
    constructors = length(_constructors) > 1 ? _constructors[2:end] : []

    _methods = methodswith(type)

    return TypeInfo(
        string(type),
        super,
        sub,
        fields,
        constructors,
        _methods,
        doc,
        docstring,
    )
end

"""
    TypeInfo(fun::Function)

Exctract information from a function object
"""
function TypeInfo(fun::Function)
    # get docstring
    doc, docstring = get_docstring(fun)

    # get methods with same name
    _methods = split_lines(string(methods(fun)))
    _methods = length(_methods) > 1 ? _methods[2:end] : []

    return TypeInfo(string(fun), nothing, nothing, nothing, [], _methods, doc, docstring)
end

# ---------------------------------------------------------------------------- #
#                                    INSPECT                                   #
# ---------------------------------------------------------------------------- #

"""
    inspect(type::DataType; width::Int=120)

Introspect a  type.

Extract  info like docstring, fields, types etc. and show it in a structured
terminal output.
"""
function inspect(
    type::DataType;
    width::Union{Nothing,Int} = nothing,
    max_n_methods::Int = 3,
)
    width = isnothing(width) ? min(console.width, 92) - 4 : width - 4
    # extract type info
    info = TypeInfo(type)

    # ------------------------------ types hierarchy ----------------------------- #
    hierarchy = TextBox(
        "",
        style_super_types(info),
        "",
        style_sub_types(info),
        width = width,
        title = "Types hierarchy",
        title_style = "bold underline yellow",
    )

    # ----------------------------------- docs ----------------------------------- #
    # @info "making docstirng box" info.docstring width
    docs = TextBox(
        info.docstring,
        title = "Docstring",
        title_style = "bold underline yellow",
        width = width,
    )

    # ---------------------------------- fields ---------------------------------- #
    if !isnothing(info.fields)
        formatted_fields::Vector{AbstractString} = []
        if !isnothing(info.fields)
            for (name, type) in zip(info.fields["names"], info.fields["types"])
                push!(
                    formatted_fields,
                    "[bold white]$(string(name))[/bold white]" *
                    highlight("::$(type)", theme, :type),
                )
            end
        end

        fields_panel = Panel(
            isnothing(formatted_fields) ? "[dim]No arguments[/dim]" : formatted_fields,
            title = "Arguments",
            title_style = "bold yellow",
            style = "dim yellow",
            width = width - 2,
        )
        insights_panel = (docs / Spacer(width - 2, 2) / fields_panel)
    else
        insights_panel = docs
    end

    # ------------------------------- constructors ------------------------------- #
    constructors = do_by_line((x) -> style_method_line(x; trim = true), info.constructors)
    n_constructors =
        length(split_lines(constructors)) > 1 ? Int(length(split_lines(constructors)) / 2) :
        0
    if n_constructors > max_n_methods
        constructors =
            join_lines(split_lines(constructors)[1:max_n_methods*2]) *
            "\n\n[grey53]( additional constructors not shown... )[/grey53]"
    end
    constructors =
        n_constructors > 1 ? constructors : "[dim]No constructors          [/dim]"

    constructors_panel = TextBox(
        constructors,
        title = "Constructors[dim]($n_constructors)",
        title_style = "bold underline yellow",
        width = width,
    )

    # ---------------------------------- methods --------------------------------- #
    if length(info.methods) > 0
        methods = do_by_line(style_method_line, info.methods)
        n_methods =
            length(split_lines(methods)) > 1 ? Int(length(split_lines(methods)) / 2) : 1

        if n_methods > max_n_methods
            methods =
                join_lines(split_lines(methods)[1:max_n_methods*2]) *
                "\n\n[grey53]( additional methods not shown... )[/grey53]"
        end
    else
        methods = "[dim]No methods          [/dim]"
        n_methods = 0
    end

    methods_panel = TextBox(
        methods,
        title = "Methods[dim]($n_methods)",
        title_style = "bold underline yellow",
        width = width,
    )

    # ------------------------------- CREATE PANEL ------------------------------- #
    _title = isabstracttype(type) ? " [dim](Abstract)[/dim]" : ""
    panel = Panel(
        Spacer(width - 2, 1),
        hierarchy,
        hLine(width - 6; style = "blue dim"),
        insights_panel,
        hLine(width - 6; style = "blue dim"),
        constructors_panel,
        hLine(width - 6; style = "blue dim"),
        methods_panel,
        title = "$(typeof(type)): [bold]$(info.name)" * _title,
        title_style = "red",
        style = "blue",
        width = width + 4,
    )

    println(panel)

end

"""
    inspect(fun::Function; width::Int=88, max_n_methods::Int = 7)

Inspects `Function` objects providing docstrings, and methods signatures.
"""
function inspect(fun::Function; width::Union{Nothing,Int} = nothing, max_n_methods::Int = 7)
    width = isnothing(width) ? min(console.width, 92) : width

    info = TypeInfo(fun)

    # ----------------------------- prepare contents ----------------------------- #
    docs = TextBox(
        info.docstring,
        title = "Docstring",
        title_style = "bold underline yellow",
        width = width - 2,
    )

    if length(info.methods) > 0
        methods = do_by_line(m -> style_method_line(m; trim = true), info.methods)
        n_methods =
            length(split_lines(methods)) > 1 ? Int(length(split_lines(methods)) / 2) : 1

        if n_methods > max_n_methods
            methods =
                join_lines(split_lines(methods)[1:max_n_methods*2]) *
                "\n\n[grey53]( additional methods not shown... )[/grey53]"
        end
    else
        methods = "[dim]No methods          [/dim]"
        n_methods = 0
    end

    methods_panel = TextBox(
        methods,
        title = "Methods[dim]($n_methods)",
        title_style = "bold underline yellow",
        width = width - 2,
    )

    # -------------------------------- print panel ------------------------------- #
    println(
        Panel(
            Spacer(width - 2, 1),
            docs,
            hLine(width - 6; style = "blue dim"),
            methods_panel,
            title = "Function: [bold red]$(info.name)[/bold red]",
            title_style = "red",
            style = "blue",
            width = width,
        ),
    )
end

"""
generic inspect method, dispatches to type-specific methods when they can be found
"""
function inspect(obj; kwargs...)
    if typeof(obj) <: Function
        inspect(obj; kwargs...)
    elseif typeof(typeof(obj)) == DataType
        inspect(typeof(obj); kwargs...)
    else
        throw("Cannot inspect $obj ($(typeof(obj)))")
    end
end
end


text:
using NamedArrays

a = [10, 20, 30]
n = NamedArray(a)
filter!(x -> x > 15, n)


text:
const working_dir_path = expanduser(normpath(@__DIR__, "..",".."))
"""
    savetofolder() -> fullpath::String

Save the most recent DataFrame into a new CSV. File into the corresponding folder of the
.../HMCResearchRandomWalks/out/... directory.

Name the file using the most recent seed and time when savetofolder() is called.
"""
function savetofolder()
    df = current_df()
    folderpath = current_path()
    fn = nameCSVfile()
    fullpath = abspath(string(working_dir_path, folderpath, fn))
    return savehelper(fullpath, df)
end

"""
    savedata([filename, df]) -> fullpath::String

Save a DataFrame into a new CSV file labelled by file name. File will be found in .../HMCResearchRandomWalks/out

If no filename is given, name the file using the most recent seed and time when savedata() is called.
If no df is given, save the most recent DataFrame loaded by the last getdata_...() call in this terminal.
"""
function savedata(filename::AbstractString=nameCSVfile(), df::DataFrame=current_df())
    fullpath = abspath(string(working_dir_path, "out/", filename))
    return savehelper(fullpath, df)
end

"""
    loaddata([path, relative])

Load a CSV file in .../HMCResearchRandomWalks/out as a DataFrame.

If no path is given, load the most recent CSV file saved by savedata() or savetofolder()
By default relative=false, so the given path must begin from the user system.
If relative=true, the given path must begin with "out/..." or "/out/..."
"""
function loaddata(path::AbstractString=current_file(), relative::Bool=false)::DataFrame
    if relative
        fullpath = abspath(string(working_dir_path, path))
    else
        fullpath = path
    end

    if !isfile(fullpath)
        @warn("No such file found. Please double check the path.")
    else
        csvdata = CSV.File(fullpath)
        df = DataFrame(csvdata)
        return df
    end
end

"Helper function to write the CSV file and checks for a valid file name"
function savehelper(fullpath::AbstractString, df::DataFrame)
    _, ext = splitext(fullpath)
    if ext != ".csv"
        @warn("File name must end with the .csv extension")
    else
        CSV.write(fullpath, df)
        current_file(fullpath)
        println("CSV File saved at the given path:")
        println(fullpath)
        return fullpath
    end
end

"Helper function to name using time and the seed loaded in CURRENT"
function nameCSVfile()
    rightnow = string(Dates.now())
    fntemplate = string(rightnow, "_seed_", current_seed())
    clean = replace(replace(fntemplate, ":" => "_"), "." => "_")
    fn = string(clean, ".csv")
    return fn
end

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
module LeafGasExchange

using Cropbox

include("vaporpressure.jl")
include("weather.jl")
include("diffusion.jl")
include("nitrogen.jl")
include("base.jl")
include("c3.jl")
include("c4.jl")
include("boundarylayer.jl")
include("stomata.jl")
include("intercellularspace.jl")
include("irradiance.jl")
include("energybalance.jl")

@system ModelBase(
    Weather, Nitrogen,
    BoundaryLayer, StomataBase, IntercellularSpace, Irradiance, EnergyBalance
)

@system ModelC3BB(ModelBase, StomataBallBerry, C3, Controller)
@system ModelC4BB(ModelBase, StomataBallBerry, C4, Controller)

@system ModelC3MD(ModelBase, StomataMedlyn, C3, Controller)
@system ModelC4MD(ModelBase, StomataMedlyn, C4, Controller)

export ModelC3BB, ModelC3MD, ModelC4BB, ModelC4MD

end


text:
export
    CUDNN_POOLING_MAX,
    CUDNN_POOLING_AVERAGE_COUNT_INCLUDE_PADDING,
    CUDNN_POOLING_AVERAGE_COUNT_EXCLUDE_PADDING,

    CUDNN_NOT_PROPAGATE_NAN,
    CUDNN_PROPAGATE_NAN

function pooling_desc(mode, window, padding, stride, maxpoolingNanOpt)
    N = length(window)
    p = Ptr{Void}[0]
    cudnnCreatePoolingDescriptor(p)
    cwindow = Cint[window[i] for i=N:-1:1]
    cpadding = Cint[padding[i] for i=N:-1:1]
    cstride = Cint[stride[i] for i=N:-1:1]
    cudnnSetPoolingNdDescriptor(p[1], mode, maxpoolingNanOpt, N, cwindow, cpadding, cstride)
    p[1]
end

function pooling(mode, window, padding, stride, x; maxpoolingNanOpt=CUDNN_NOT_PROPAGATE_NAN, alpha=1.0, beta=0.0)
    N = length(window)
    outdims = [(size(x,i) + 2padding[i] - window[i]) ÷ stride[i] + 1 for i=1:N]
    y = similar(x, outdims..., size(x,N+1), size(x,N+2))

    T = eltype(x)
    h = handle(x)
    poolingdesc = pooling_desc(mode, window, padding, stride, maxpoolingNanOpt)
    xdesc = tensor_desc(x)
    ydesc = tensor_desc(y)
    cudnnPoolingForward(h, poolingdesc, T[alpha], xdesc, x, T[beta], ydesc, y)

    cudnnDestroyPoolingDescriptor(poolingdesc)
    cudnnDestroyTensorDescriptor(xdesc)
    cudnnDestroyTensorDescriptor(ydesc)
    y
end

function ∇pooling!(mode, windims, padding, stride, y, dy, x, dx;
    maxpoolingNanOpt=CUDNN_NOT_PROPAGATE_NAN, alpha=1.0, beta=0.0)

    h = gethandle(device(x))
    poolingdesc = pooling_desc(mode, window, padding, stride, maxpoolingNanOpt)
    ydesc = tensor_desc(y)
    dydesc = tensor_desc(dy)
    xdesc = tensor_desc(x)
    dxdesc = tensor_desc(dx)
    cudnnPoolingBackward(h, poolingdesc, T[alpha], ydesc, y, dydesc, dy, xdesc, x,
    T[beta], dxdesc, dx)

    cudnnDestroyPoolingDescriptor(poolingdesc)
    cudnnDestroyTensorDescriptor(ydesc)
    cudnnDestroyTensorDescriptor(dydesc)
    cudnnDestroyTensorDescriptor(xdesc)
    cudnnDestroyTensorDescriptor(dxdesc)
end


text:
constructor_of(::T) where T = constructor_of(T)
constructor_of(T::Type) = getfield(T.name.module, Symbol(T.name.name))
constructor_of(T::UnionAll) = constructor_of(T.body)

struct_as_dict(st) = [(n => getfield(st, n)) for n in fieldnames(typeof(st))]

"""
    recon(thing; <keyword arguments>)

Re-construct `thing` with new field values specified by the keyword
arguments.
"""
function recon(thing; kwargs...)
    constructor = constructor_of(thing)
    return constructor(; struct_as_dict(thing)..., kwargs...)
end


text:
module NetworkDepot

import Base: download

export
# Abstract Dataset
Dataset

abstract Dataset

const DDIR = joinpath(Pkg.dir(), "NetworkDepot", "downloads")

function filepath(x::String)
   joinpath(DDIR, x)
end

include("Epinions.jl")
include("LiveJournal.jl")

end # module


text:
"""
    unitsphere(ndirection; half = false, normal = [0, 0, 1])

Create `ndirection` directions unformly distributed on the unit sphere.

If `half` is `true`, the points will be discributed on the hemisphere defined by `normal`.
"""
function unitsphere(ndirection; half = false, normal = [0, 0, 1])
    # Create Fibonacci distributed directions
    directions = create_fibonacci_sphere(ndirection)
end


text:
# Query adjacency matrix in a database table using an iterator. This example queries all the elements of a set of starting vertices, calculates their degrees, and finds the maximum degree.

# Max elements in iterator.
MaxElem = 500

# Create array of starting vertices
Nv0 = 100
v0 = string.(ceil.(Int,10000*rand(Nv0,1)))

# Set up query iterator
TadjIt = getiterator(Tadj,MaxElem)

# Start query iterator
A = TadjIt[v0str,:]

# Initialize AinDeg to accumulate output
AinDeg = Assoc("","","")

# While there are still elements in A
while nnz(A) > 0
    global AinDeg, A
    # Compute in degree
    AinDeg = AinDeg + sum(str2num(A),1)

    # Get next query
    A = TadjIt[]
end

# Get vertex with the maximum degree
AmaxInDeg = (AinDeg == maximum(getadj(AinDeg)))


text:
export SSparseMatrixCSC

"""
    SSparseMatrixCSC{Tv,Ti<:Integer, NNZ, NP} <: AbstractSparseMatrix{Tv,Ti}

static version of SparseMatrixCSC
"""
struct SSparseMatrixCSC{Tv,Ti<:Integer, NNZ, NP} <: AbstractSparseMatrix{Tv,Ti}
    m::Int                  # Number of rows
    n::Int                  # Number of columns
    colptr::SVector{NP,Ti}      # Column i is in colptr[i]:(colptr[i+1]-1)
    rowval::SVector{NNZ,Ti}      # Row values of nonzeros
    nzval::SVector{NNZ,Tv}       # Nonzero values

    function SSparseMatrixCSC{Tv,Ti, NNZ, NP}(m::Integer, n::Integer, colptr::SVector{NP, Ti}, rowval::SVector{NNZ,Ti},
                                    nzval::SVector{NNZ,Tv}) where {Tv,Ti<:Integer, NNZ, NP}
        m < 0 && throw(ArgumentError("number of rows (m) must be ≥ 0, got $m"))
        n < 0 && throw(ArgumentError("number of columns (n) must be ≥ 0, got $n"))
        new(Int(m), Int(n), colptr, rowval, nzval)
    end
end

function SSparseMatrixCSC(m::Integer, n::Integer, colptr::SVector, rowval::SVector, nzval::SVector)
    Tv = eltype(nzval)
    Ti = promote_type(eltype(colptr), eltype(rowval))
    SSparseMatrixCSC{Tv,Ti,length(nzval),n+1}(m, n, colptr, rowval, nzval)
end

function size(spm::SSparseMatrixCSC{Tv, Ti, NNZ, NP}, i::Integer) where {Tv, Ti, NNZ, NP}
    i == 1 ? spm.m : (i == 2 ? NP-1 : throw(ArgumentError("dimension out of bound!")))
end
size(spm::SSparseMatrixCSC{Tv, Ti, NNZ, NP}) where {Tv, Ti, NNZ, NP} = (spm.m, NP-1)

function getindex(ssp::SSparseMatrixCSC{Tv}, i, j) where Tv
    S = ssp.colptr[j]
    E = ssp.colptr[j+1]-1
    for ii in S:E
        if i == ssp.rowval[ii]
            return ssp.nzval[ii]
        end
    end
    return Tv(0)
end

issparse(::SSparseMatrixCSC) = true
nonzeros(M::SSparseMatrixCSC) = M.nzval
nnz(spm::SSparseMatrixCSC{Tv, Ti, NNZ}) where {Tv, Ti, NNZ} = NNZ
dropzeros!(M::SSparseMatrixCSC; trim::Bool=false) = M

isdense(::SSparseMatrixCSC) = false


text:
using OrdinaryDiffEq, Test

function time_derivative(du,u,p,t)
  du[1] = -t
end
function time_derivative_analytic(u0,p,t)
  u0 .- t.^2 ./ 2
end

ff_time_derivative = ODEFunction(time_derivative,
                                 analytic=time_derivative_analytic)
u0 = [1.0]
tspan = (0.0,1.0)
prob = ODEProblem(ff_time_derivative,u0,tspan)

sol = solve(prob,Rosenbrock32(),reltol=1e-9,abstol=1e-9)
@test sol.errors[:final] < 1e-5
sol = solve(prob,Rosenbrock23())
@test sol.errors[:final] < 1e-10
sol = solve(prob,Rodas4())
@test sol.errors[:final] < 1e-10
sol = solve(prob,Rodas5())
@test sol.errors[:final] < 1e-10
sol = solve(prob,Veldd4())
@test sol.errors[:final] < 1e-10
sol = solve(prob,KenCarp4())
@test sol.errors[:final] < 1e-10
sol = solve(prob,TRBDF2())
@test sol.errors[:final] < 1e-10

sol = solve(prob,GenericImplicitEuler(),dt=1/10)
@test sol.errors[:final] < 1e-1
sol = solve(prob,GenericTrapezoid(),dt=1/10)
@test sol.errors[:final] < 1e-12
sol = solve(prob,Euler(),dt=1/100)
@test sol.errors[:final] < 6e-3

for alg in CACHE_TEST_ALGS
  sol = solve(prob,alg,dt=1/10)
  if !(typeof(alg) <: Euler)
    @test sol.errors[:final] < 4e-14
  end
end


text:
using Distributions

distro = Binomial(9, 0.5)
likelihood = pdf(distro, 6)
println(likelihood)


text:
#=
Parallelization is hard and is not my specialty, so I'll focus on a simple example.
See more [here](https://docs.julialang.org/en/v1/manual/parallel-computing/).

**Atention**: Before we start, you should either:
- Set **"Julia: Num Threads"** in your VSCode to the number of threads that you can use; or
- Pass the number of threads to the `julia` call: `julia --threads N`.

Check the number of threads being used with
=#
Threads.nthreads()
#=

We'll compute the approximate value of $\pi$ using Monte Carlo simulation.
This example is based on [Parallel Programming in Python](https://carpentries-incubator.github.io/lesson-parallel-python/).

The idea is pretty simple, and you might have seen it before.
Consider a square with vertices $(0,0)$, $(1,0)$, $(0,1)$ and $(1,1)$,
and the sector of the circle centered at $(0,0)$ with radius $1$.
=#
using Plots # hide
pyplot() # hide
plot(leg=false, axis=false, axis_ratio=:equal, ticks=false) # hide
plot!([0,1,1,0,0], [0,0,1,1,0], c=:blue, fill=true, fillcolor=:lightblue) # hide
t = range(0, π / 2, length=100) # hide
plot!([0; cos.(t); 0], [0; sin.(t); 0], c=:red, fill=true, fillcolor=:pink) # hide
png(joinpath(@OUTPUT, "parallel-1")) # hide
#=
\fig{parallel-1}

The square's area is $1$, and the circle's area is $\pi / 4$.
The ratio being $\pi / 4$, then.

The Monte Carlo approach is to generate a number $N$ os points equally distributed on the square.
Let $M$ be the number of points that fall inside the sector.
The ratio $M / N$ should assintotically approach the ratio $\pi / 4$, therefore an approximation to $\pi$ can be obtained by computing $4M / N$.

Here's a basic implementation:
=#
function calc_pi(N)
  M = 0
  for i = 1:N
    x, y = rand(), rand()
    if x^2 + y^2 < 1
      M += 1
    end
  end
  return 4 * M / N
end
#=
Here are some histogram of the error using various $N$ various for 100 simulations each:
=#
pyplot() # hide
N_values = 10 .^ (5:7)
plt = plot(layout=grid(3, 1), size=(400,600), leg=false)
for (i, N) in enumerate(N_values)
  ϵ = π .- [calc_pi(N) for _ = 1:200]
  histogram!(plt[i], ϵ, bins=20)
  title!(plt[i], "N = $N")
end
plt
png(joinpath(@OUTPUT, "parallel-2")) # hide
#=
\fig{parallel-2}

Using `@time` we can assess the computation time of `calc_pi`:
=#
@time calc_pi(10^7)
#=
Be aware that Julia compiles the code as it is first executed.
In this case, since we ran the code before, the timing is mostly accurate.

That being said, measuring execution time of code has some randomness, so it's better to run the code a few times and compute some statistics (mean, median, minimum, maximum, etc.) on the result.
Luckily, the package `BenchmarkTools` does that for us.

You can use `@btime "Your stuff"` to compute the equivalent of `@time`, but we can see more information with `@benchmark`:
=#
using BenchmarkTools

println("π ≈ $(calc_pi(10^7))")
@benchmark calc_pi(10^7)
#=
Unlike Python, vectorizing the code won't improve the performance:
=#
function calc_pi_vectorized(N)
  pts = rand(N, 2)
  M = count(sum(pts .^ 2, dims=2) .< 1)
  return 4 * M / N
end

println("π ≈ $(calc_pi_vectorized(10^7))")
@benchmark calc_pi_vectorized(10^7)
#=
We can, however, use threads to improve it.
The code below has a few simple changes to have an improved version.
=#
using Base.Threads: @threads, threadid, nthreads

function calc_pi_threaded_naive(N)
  M = zeros(Int, nthreads())
  @threads for i = 1:N
    x, y = rand(), rand()
    if x^2 + y^2 < 1
      M[threadid()] += 1
    end
  end
  return 4 * sum(M) / N
end

println("π ≈ $(calc_pi_threaded_naive(10^7))")
@benchmark calc_pi_threaded_naive(10^7)
#=
The version above creates the a vector to store the number of points inside the circle per thread.
We could be more naive and use a single `M`, but for sharing between threads, `M` must be an array.

This version simply uses the `Base.Threads` module in the form of the macro `@threads` before the `for`.
The macro splits the for between the threads automatically.

As you saw, though, the result is **slower** than the basic version.
To fix this, we have to use thread-safe random number generators.
The code below does this and improves the execution time in comparison to the original.
=#
using Random

function calc_pi_threaded(N)
  M = zeros(Int, nthreads())
  @threads for i = 1:N
    rng = Random.THREAD_RNGs[threadid()]
    x, y = rand(rng), rand(rng)
    if x^2 + y^2 < 1
      M[threadid()] += 1
    end
  end
  return 4 * sum(M) / N
end

println("π ≈ $(calc_pi_threaded(10^7))")
@benchmark calc_pi_threaded(10^7)
#=
But wait, there's more.
We can manually split the range, and just use `@threads` over the ranges.
=#
using Distributed: splitrange

function calc_pi_proper_threaded(N)
  ranges = splitrange(1, N, nthreads())
  M = Ref(0)
  @threads for rg in ranges
    rng = MersenneTwister()
    c = 0
    for _ in rg
      c += rand(rng)^2 + rand(rng)^2 < 1
    end
    M[] += c
  end
  return 4 * M[] / N
end

println("π ≈ $(calc_pi_proper_threaded(10^7))")
@benchmark calc_pi_proper_threaded(10^7)

#=
We use `Ref(0)` now to keep a reference to a number that all threads can access.

Finally, we can use the package FLoops and change very little to obtain still some more speed.
I don't know why, though.
=#

using FLoops

function calc_pi_proper_floop(N)
  ranges = splitrange(1, N, nthreads())
  @floop for rg in ranges
    rng = MersenneTwister()
    c = 0
    for _ in rg
      c += rand(rng)^2 + rand(rng)^2 < 1
    end
    @reduce(M += c)
  end
  return 4 * M / N
end

println("π ≈ $(calc_pi_proper_floop(10^7))")
@benchmark calc_pi_proper_floop(10^7)

max_stars_repo_name: chunggnuhc/julia-tutorial-at-nlesc

max_issues_repo_name: chunggnuhc/julia-tutorial-at-nlesc

max_forks_repo_name: chunggnuhc/julia-tutorial-at-nlesc