text:
##
## Turkey!
##
"""
    function turkey()

Turkey!

# Example
```jldoctest
julia> cowsay("Gobble, gobble", cow=Cowsay.turkey)
 ________________
< Gobble, gobble >
 ----------------
  \\                                  ,+*^^*+___+++_
   \\                           ,*^^^^              )
    \\                       _+*                     ^**+_
     \\                    +^       _ _++*+_+++_,         )
              _+^^*+_    (     ,+*^ ^          \\+_        )
             {       )  (    ,(    ,_+--+--,      ^)      ^\\
            { (@)    } f   ,(  ,+-^ __*_*_  ^^\\_   ^\\       )
           {:;-/    (_+*-+^^^^^+*+*<_ _++_)_    )    )      /
          ( /  (    (        ,___    ^*+_+* )   <    <      \\
           U _/     )    *--<  ) ^\\-----++__)   )    )       )
            (      )  _(^)^^))  )  )\\^^^^^))^*+/    /       /
          (      /  (_))_^)) )  )  ))^^^^^))^^^)__/     +^^
         (     ,/    (^))^))  )  ) ))^^^^^^^))^^)       _)
          *+__+*       (_))^)  ) ) ))^^^^^^))^^^^^)____*^
          \\             \\_)^)_)) ))^^^^^^^^^^))^^^^)
           (_             ^\\__^^^^^^^^^^^^))^^^^^^^)
             ^\\___            ^\\__^^^^^^))^^^^^^^^)\\\\
                  ^^^^^\\uuu/^^\\uuu/^^^^\\^\\^\\^\\^\\^\\^\\^\\
                     ___) >____) >___   ^\\_\\_\\_\\_\\_\\_\\)
                    ^^^//\\\\_^^//\\\\_^       ^(\\_\\_\\_\\)
                      ^^^ ^^ ^^^ ^
```
"""
function turkey(;eyes="oo", tongue="  ", thoughts="\\")
the_cow =
"""
  $thoughts                                  ,+*^^*+___+++_
   $thoughts                           ,*^^^^              )
    $thoughts                       _+*                     ^**+_
     $thoughts                    +^       _ _++*+_+++_,         )
              _+^^*+_    (     ,+*^ ^          \\+_        )
             {       )  (    ,(    ,_+--+--,      ^)      ^\\
            { (@)    } f   ,(  ,+-^ __*_*_  ^^\\_   ^\\       )
           {:;-/    (_+*-+^^^^^+*+*<_ _++_)_    )    )      /
          ( /  (    (        ,___    ^*+_+* )   <    <      \\
           U _/     )    *--<  ) ^\\-----++__)   )    )       )
            (      )  _(^)^^))  )  )\\^^^^^))^*+/    /       /
          (      /  (_))_^)) )  )  ))^^^^^))^^^)__/     +^^
         (     ,/    (^))^))  )  ) ))^^^^^^^))^^)       _)
          *+__+*       (_))^)  ) ) ))^^^^^^))^^^^^)____*^
          \\             \\_)^)_)) ))^^^^^^^^^^))^^^^)
           (_             ^\\__^^^^^^^^^^^^))^^^^^^^)
             ^\\___            ^\\__^^^^^^))^^^^^^^^)\\\\
                  ^^^^^\\uuu/^^\\uuu/^^^^\\^\\^\\^\\^\\^\\^\\^\\
                     ___) >____) >___   ^\\_\\_\\_\\_\\_\\_\\)
                    ^^^//\\\\_^^//\\\\_^       ^(\\_\\_\\_\\)
                      ^^^ ^^ ^^^ ^
"""
return the_cow
end


text:
__precompile__(true)
@doc raw"""
# ElectricalEngineering

This is a Julia package on electrical engineering based on
[Unitful](https://github.com/PainterQubits/Unitful.jl) and
[PyPlot](https://github.com/JuliaPy/PyPlot.jl). The package
ElectricalEngineering.jl is tested with Julia 1.

## Plotting

- `arrowaxes` Creates a plot with a horizontal and a vertical axis instead of a
  frame

- `removeaxes` Removes the axis of the actual plot

- `lengthdimension` Draws length dimension with label based on rectangular x-
  and y-coordinates. The length dimension arrow can be shifted parallel to the
  specified  coordinates by means the single parameter `par`. Additionally,
  auxiliary lines  from the specified coordinates to the length dimension are
  drawn.

## Complex Phasors

- `pol` Creates a complex quantity with length `r` and angle `phi`

- `∠` Creates a complex quantity with length 1 and angle `phi`

- `phasor` Draws a phasor from a starting point `origin` and end point
  `origin`+`c`. The phasor consists of a shaft and an arrow head. Each phasor
  `c`  is plotted as a relative quantity, i.e., `c/ref` is actually shown the
  plot figure. This concept of plotting a per unit phasor is used to be able to
  plot phasor with different quantities, e.g., voltage and current phasors. It
  is important that the variables `c`, `origin` and `ref` have the same units
  (defined through Unitful).

- `phasorsine` Draws a phasor with magnitude between 0 and 1 on the left subplot
  of a figure and a sine diagram on the right subplot of the figure. Such graph
  is used in electrical engineering to explain the relationship between phasors
  and time domain wave forms.

- `angulardimension` Draws an arrowed arc, intended to label the angle between
  phasors. The arc is drawn from angle `phi1` (begin) to `phi2` (end). The begin
  and end arrow shapes can be set. The arc can be labeled and optionally a dot
  maker can be used to indicate right angles (90°).

- `phasordimension` Draws length dimension with label based on rectangular x-
  and y-coordinates. The length dimension arrow can be shifted parallel to the
  specified  coordinates by means the single parameter `par`. Additionally,
  auxiliary lines  from the specified coordinates to the length dimension are
  drawn.

## Complex Circuit Calculations

- `∥` Calculates the parallel connection of impedances. Optionally, the
  calculation can be performed using units from the module Unitful.

## Input and Output

- `printuln` Prints a real or complex variable in an optionally specified unit

- `usprint` Returns a string containing the number with six digits plus unit
  without space in between, such that the string output can be copied and
  pasted.

- `save2fig` Save the actual figure in the file formats `png` and `pdf`
    including subdirectories.  These two graphics format are relevant when
    processing figures in LaTeX, LyX, LibreOffice or other word processors

- `save3fig` Save the actual figure in the file formats `png`, `eps` and `pdf`
  including subdirectories.  These three graphics format are relevant when
  processing figures in LaTeX and LyX

## See also `?` on

- `ElectricalEngineering.Gray`

- `ElectricalEngineering.Tab20bc`

# Other Packages related with Electrical Engineering

- [CircuitComponentRounding.jl](https://github.com/KronosTheLate/CircuitComponentRounding.jl) provides functions to round given values to the nearest standardized value for circuit components

- [Modia.jl](https://github.com/ModiaSim/Modia.jl) is a package for modeling and simulation of multidomain engineering systems described by differential equations, algebraic equations, and (space-discretized) partial differential equations

"""
module ElectricalEngineering

    # Implementation of ElectricalEngineering is based on
    # - PyPlot              Plotting phasors
    # - Unitful             Handling of quantity units
    using Printf, Base, Statistics, PyPlot, Unitful

    # Default colors for plotting (with PyPlot)
    include("plot.jl")
    include("Tab20bc.jl")
    include("Gray.jl")

    # Complex phasor calculation and plotting
    include("phasor.jl")

    # Calculations in ElectricalEngineering circuits
    include("circuit.jl")

    # Winding calculations and winding plot routines
    include("winding.jl")

    # Input and output functions to print phasors, figures, etc.
    include("io.jl")

    # Unitful additional functions
    include("unitful.jl")
end

ext: jl


text:
#all OpenGL Types
typealias GLvoid                                                      Void
typealias GLint                                                       Cint
typealias GLsizeiptr                                                  Cint
typealias Pointer                                                     Ptr{Void}
typealias GLhalfNV                                                    Cushort
typealias GLshort                                                     Cshort
typealias GLdouble                                                    Cdouble
typealias GLushort                                                    Cushort
typealias GLuint                                                      Cuint
typealias GLint64                                                     Clonglong
typealias GLbyte                                                      Cchar
typealias GLchar                                                      Cuchar
typealias GLubyte                                                     Cuchar
typealias GLintptr                                                    Cptrdiff_t
typealias GLfloat                                                     Cfloat
typealias GLhalfARB                                                   Cushort
typealias GLhalf                                                      Cushort
typealias GLenum                                                      Cuint
typealias GLboolean                                                   Cuchar
typealias GLclampf                                                    Cfloat
typealias GLsizei                                                     Cint
typealias GLsync                                                      Ptr{Void}
typealias GLuint64                                                    Culonglong
typealias GLclampd                                                    Cdouble
typealias GLbitfield                                                  Cuint
export GLvoid
export GLint
export GLsizeiptr
export Pointer
export GLhalfNV
export GLshort
export GLdouble
export GLushort
export GLuint
export GLint64
export GLbyte
export GLchar
export GLubyte
export GLintptr
export GLfloat
export GLhalfARB
export GLhalf
export GLenum
export GLboolean
export GLclampf
export GLsizei
export GLsync
export GLuint64
export GLclampd
export GLbitfield


text:
function knotHash(d)
    d = vcat(map(Int, collect(d)), [17,31,73,47,23])
    l = collect(0:255)
    i = 0
    skip = 0
    for _ in 1:64
        for j in d
            l = circshift(l, -i)
            reverse!(l, 1, j)
            l = circshift(l, i)
            i += j + skip
            skip += 1
        end
    end
    l = mapslices(x -> reduce(⊻, x), reshape(l, (16,16)), 1)

    return mapreduce(x -> hex(x, 2), *, l)
end

input = "ljoxqyyw"

# Part 1
hashes = map(x -> knotHash("$input-$x"), 0:127)
grid = map(h -> vcat(map(x -> reverse!(digits(x, 2, 8)), hex2bytes(h))...)', hashes)
println(mapreduce(sum, +, grid))

# Part 2
n = length(grid)
function visit(i, j)
    if !(1 <= i <= n && 1 <= j <= n) || grid[i][j] != 1
        return 0
    end
    grid[i][j] = 0
    for d in [[1, 0], [0, -1], [-1, 0], [0, 1]]
        visit(i + d[1], j + d[2])
    end
    return 1
end
println(sum(visit(i, j) for i in 1:n, j in 1:n))


text:
using MultidimensionalSpectroscopy, PyPlot, QuantumOptics, LinearAlgebra, FFTW, Colors,
        Printf, DelimitedFiles

## some theory
# chem.libretexts.org/Bookshelves/Physical_and_Theoretical_Chemistry_Textbook_Maps/Book%3A_Time_Dependent_Quantum_Mechanics_and_Spectroscopy_(Tokmakoff)/13%3A_Coupling_of_Electronic_and_Nuclear_Motion/13.01%3A_The_Displaced_Harmonic_Oscillator_Model
# www.scm.com/doc/ADF/Input/Vibrationally_resolved_el_spectra.html

# make sure to set script directory as pwd()
cd(@__DIR__)

pygui(true)

cmp = create_colormap("bright");

calc_2d = true

## Huang-Rhys factor
D_HR = .25
d = sqrt(D_HR)
D_HR2 = .55
d2 = sqrt(D_HR2)

b_tls = NLevelBasis(2)          # Hilbert-space of system                   Basis: {|ui⟩}
b_vib = FockBasis(3)            # Hilbert-space of oscillator               Basis: {|vi⟩}
b     = b_tls ⊗ b_vib          # combined basis

j21 = transition(b_tls,2,1)                     # |e⟩⟨g|
j12 = dagger(j21)                               # |g⟩⟨e|
at  = create(b_vib)                             # ...
a   = dagger(at)
D   = displace(b_vib,d)
D2  = displace(b_vib,d2)

# ground states
Eg  = 0
ωg  = 0.10415
ωg2 = 0.16862
Hg_el = Eg * j12 * j21     # == HG = g * Eg ⊗ dagger(g)
#Hg_vib = ωg * (at*a + a*at * 1/2)
Hg_vib = ωg * (number(b_vib) + identityoperator(b_vib) * 1/2)
Hg_vib2 = ωg2 * (number(b_vib) + identityoperator(b_vib) * 1/2)
#Hg_vib = ωg * (at*a + identityoperator(b_vib) * 1/2)
#Hg_vib2 = ωg2 * (at*a + identityoperator(b_vib) * 1/2)
Hg = Hg_el ⊗ one(b_vib) ⊗ one(b_vib) + j12*j21 ⊗ (1 * Hg_vib) ⊗ one(b_vib) + j12*j21 ⊗ one(b_vib) ⊗ (1 * Hg_vib2)
#Hg = Hg_el ⊗ Hg_vib

# HAVE TO created excited states w/o displaced potential and then use displacement
# operator to act on transition dipole moment ... or calculate FC factors and use
# this as transition dipole operator ... #CHECK

# excited state
Ee  = 1.65
ωe  = .9 * 0.10415
ωe2 = .9 * 0.16862
He_el = Ee * j21 * j12
#He_vib = dagger(D) * Hg_vib * D # issue with D and dagger(D)
#He_vib = D' * Hg_vib * D # issue with D and dagger(D)
#He_vib2 = D' * Hg_vib * D
#He_vib = Hg_vib  # TODO figure out what's going on here
He_vib  = ( ωe  * (number(b_vib) + identityoperator(b_vib) * 1/2))
He_vib2 = ( ωe2 * (number(b_vib) + identityoperator(b_vib) * 1/2))
He = He_el ⊗ one(b_vib) ⊗ one(b_vib) + j21*j12 ⊗ (1 * He_vib) ⊗ one(b_vib) + j21*j12 ⊗ one(b_vib) ⊗ (1 * He_vib2)
#He = He_el ⊗ He_vib

## full system Hamiltonian
Het = (one(b_tls) ⊗ D ⊗ one(b_vib))' * He * (one(b_tls) ⊗ D ⊗ one(b_vib)) +
      (one(b_tls) ⊗ one(b_vib) ⊗ D2)' * He * (one(b_tls) ⊗ one(b_vib) ⊗ D2)

#He = (one(b_tls) ⊗ D ⊗ D2)' * He * (one(b_tls) ⊗ D ⊗ D2)

H = Hg + He
H = (H + dagger(H)) / 2

## diagonalized Hamiltonian
states   = eigvecs(dense(H).data)
energies = eigvals(dense(H).data)

energies, states = eigenstates(dense(H))

# and GS
energies_g, states_g = eigenstates(dense(Hg_vib))
energies_g2, states_g2 = eigenstates(dense(Hg_vib2))
# and ES
energies_e, states_e = eigenstates(dense(D' * He_vib * D))
energies_e2, states_e2 = eigenstates(dense(D2' * He_vib2 * D2))

## plot energy levels
N = length(energies)
println(Int8(ceil(N/2)))
figure(figsize=(6,4))
subplot(121)
scatter(zeros(length(energies_g)), energies_g)
scatter(zeros(length(energies_g2)), energies_g2)

x = [-1:0.01:2*d;]

# TODO: transform to position basis to plot PECs
#xx = PositionBasis(0,10,200)
#xxx = position(xx)

m = 50                        # arbitrary value for m
f = m * ωg^2
f2 = m * ωg2^2
# GS
plot(x,f*x.^2)
plot(x,f2*x.^2)
# ES
f = m * ωe^2
f2 = m * ωe2^2
scatter(zeros(length(energies_e)) .+ d, energies_e .+ Ee)
scatter(zeros(length(energies_e2)) .+ d2, energies_e2 .+ Ee)
plot(x,f * (x.-d).^2 .+ Ee)
plot(x,f2 * (x.-d2).^2 .+ Ee)
ylim(0, 3.5)

#for i=1:length(states_g)
#    plot(xpoints, abs2.(states_g[i].data).*40 .+ energies_g[i],color="k",linewidth=1)
#    plot(xpoints, abs2.(states_e[i].data).*40 .+ energies_e[i],color="g",linewidth=1)
#end

## calculate Franck-Condon factors
# |⟨v|i⟩² , where |v⟩ is the vibrational level on the excited electronic state and |i⟩ the initial state
FC = complex(zeros(length(energies_g),length(energies_e)))
for ii = 1:length(energies_g)
    for jj = 1:length(energies_e)
        #FC[ii,jj] = abs(conj(states_g[1]' * (D*states_e[jj])) * states_g[1]' * (D*states_e[ii]))
        #FC[ii,jj] = abs(conj((D*states_e[jj])' * states_g[1]) * (D*states_e[ii])' * states_g[1])
        FC[ii,jj] = abs(states_e[jj]' * states_g[1])^2
        #FC[ii,jj] = abs(conj(states_g[1]' * states_e[jj]) * (states_g[1]' * states_e[ii]))
    end
end

FC2 = complex(zeros(length(energies_g2),length(energies_e2)))
for ii = 1:length(energies_g2)  # only for lowest vib level on ground state ... length(energies_g2)
    for jj = 1:length(energies_e2)
        #FC2[ii,jj] = abs(conj(states_g2[1]' * (D2*states_e2[jj])) * states_g2[1]' * (D2*states_e2[ii]))
        #FC2[ii,jj] = abs(conj((D2*states_e2[jj])' * states_g2[1]) * (D2*states_e2[ii])' * states_g2[1])
        FC2[ii,jj] = abs(states_e2[jj]' * states_g2[1])^2
    end
end

## RAMAN overlaps
## calculate Franck-Condon factors
# ⟨f|v⟩⟨v|i⟩
FR = complex(zeros(length(energies_g),length(energies_e)))
for ii = 1:length(energies_g)
    for jj = 1:length(energies_e)
        FR[ii,jj] = abs(conj(states_g[2]' * D*states_e[jj]) * ( (D*states_e[ii])' * states_g[1]))^2
    end
end

## RAMAN overlaps
## calculate Franck-Condon factors
# ⟨f|v⟩⟨v|i⟩
FR2 = complex(zeros(length(energies_g),length(energies_e)))
for ii = 1:length(energies_g)
    for jj = 1:length(energies_e)
        FR2[ii,jj] = abs(conj(states_g2[2]' * D*states_e2[jj]) * ( (D*states_e2[ii])' * states_g2[1]))^2
    end
end

# the previous algorithm corresponds to calculating the FC-factors
l = complex(zeros(length(energies_g),length(energies_e)))
for ii = 1:length(energies_g)
    for jj = 1:length(energies_e)
        l[ii,jj] = exp(-d^2) * d^(2*jj) / factorial(jj)
    end
end

ln = zeros(5,)
for n = 1:5
    ln[n] = exp(-d^2) * d^(2*n) / factorial(n)
end

μ12 = j12 ⊗ D ⊗ D2# + j21 ⊗ D' ⊗ D2'
#μ12 = j12 ⊗ D ⊗ one(b_vib) + j12 ⊗ one(b_vib) ⊗ D2
μ12 = μ12 + dagger(μ12)
#μ12 = μ ⊗ D ⊗ one(b_vib) + μ ⊗ one(b_vib) ⊗ D2
#μ12 = μ ⊗ (a*at) ⊗ (a*at)  # ???
#μ12 = rho0 * μ12 + μ12 * rho0
#μ12 = rho0 * (μ ⊗ D ⊗ D2)
#μ12 = μ12 + μ12'
#μ12 = μ ⊗ (a + at) ⊗ (a + at)
μ23 = copy(μ12)
fill!(μ23.data,0)

# initialize in groud state
Psi0 = nlevelstate(b_tls,1) ⊗ fockstate(b_vib,0) ⊗ fockstate(b_vib,0)
rho0 = Psi0 ⊗ dagger(Psi0)  # alt.: rho0 = dm(Psi0)

Psi1 = nlevelstate(b_tls,2) ⊗ fockstate(b_vib,0) ⊗ fockstate(b_vib,0)
rho1 = dm(Psi1)

 ## TODO used thermally populated ground state instead
#T = 0.01                               # Temperature T in what unit ?
#rho0 = thermalstate(H,T)

Γ = [sqrt(0.0001), sqrt(0.0001), sqrt(0.045)]
L = Γ .* [one(b_tls) ⊗ a ⊗ one(b_vib), one(b_tls) ⊗ one(b_vib) ⊗ a, j21*j12 ⊗ one(b_vib) ⊗ one(b_vib)]

tlist = [0:0.1:30;]*2*π

corr = timecorrelations.correlation(tlist, rho0, H, L, μ12, μ12)

# calculate ⟨f|i(t)⟩ from https://pubs.acs.org/doi/pdf/10.1021/ar960240c eq.10
#corr = timecorrelations.correlation(tlist, rho0, H, L, μ12, j12⊗a⊗one(b_vib))
#figure()
#plot(tlist,corr)

subplot(222)
plot(tlist,real.(corr))

zp = 11
corr    = zeropad(corr,zp)
tnew, ~ = interpt(tlist,zp)

ω_abs, spec_abs = timecorrelations.correlation2spectrum(tnew, corr;
                                                        normalize_spec=true)

corr = []
corr = timecorrelations.correlation(tlist, rho1, H, L, μ12, μ12)
#tout, rhot = timeevolution.master(tlist,rho1*μ12,H,L)
#corr = real(expect(μ12,rhot))
zp = 11
corr    = zeropad(corr,zp)
tnew, ~ = interpt(tlist,zp)

ω_em, spec_em = timecorrelations.correlation2spectrum(tnew, corr;
                                                        normalize_spec=true)

subplot(224)
plot(ω_abs,spec_abs)
plot(-ω_em,spec_em,"r:")

legend(["absorption", "emission"])

γ = 0.02 # γ has an effect on the intensity envelop
abs_crossSection = real(ω_abs .* sum([γ * FC[1,i]./((i * ωe + Ee - ωg .- ω_abs).^2 .- (im * γ)^2) for i in 1:length(energies_e)]))
plot(-ω_abs,abs_crossSection./maximum(abs_crossSection)./2)

abs_crossSection2 = real(ω_abs .* sum([γ * FC2[1,i]./((i * ωe2 + Ee - ωg2 .- ω_abs).^2 .- (im * γ)^2) for i in 1:length(energies_e)]))
plot(-ω_abs,abs_crossSection2./maximum(abs_crossSection2)./2)

abs_crossSec = (abs_crossSection + abs_crossSection2) ./maximum(abs_crossSection + abs_crossSection2)

plot(-ω_abs,abs_crossSec)

#temp = FC
#temp = temp ./ maximum(real(temp))
#for jj = 1:length(energies_g)
#    plot(-([Ee-Eg, Ee-Eg] .+ (jj-1)*ωg), [0, real(temp[jj,jj])], color="r",linewidth=2)
#end

#temp = FC2
#temp = temp ./ maximum(real(temp))
#for jj = 1:length(energies_g)
#    plot(-([Ee-Eg, Ee-Eg] .+ (jj-1)*ωg2), [0, real(temp[jj,jj])], color="b",linewidth=1)
#end
#=
temp = FR
temp = temp ./ maximum(real(temp))
for jj = 1:length(energies_g)
    plot(-([Ee-Eg, Ee-Eg] .+ (jj-1)*ωg) .+ 0.005, [0, real(temp[jj,jj])], color="y",linewidth=2,linestyle=":")
end

temp = FR2
temp = temp ./ maximum(real(temp))
for jj = 1:length(energies_g)
    plot(-([Ee-Eg, Ee-Eg] .+ (jj-1)*ωg2) .+ 0.005, [0, real(temp[jj,jj])], color="g",linewidth=2,linestyle=":")
end
=#
tight_layout()

F = L # works with method "lindblad"

if calc_2d

    ## zeropad up to 10^zp
    zp = 11

    ## calculate 2d spectra at
    T = [0.]

    spectra2d = Array{out2d}(undef, length(T))

    # use multithreding...
    Threads.@threads for i = 1:length(T)
    # ... or not
    #for i = 1:length(T)
        spectra2d[i] = make2Dspectra(tlist,rho0,H,F,μ12,μ23,T[i],
                                            "lindblad";debug=true,use_sub=false,
                                                zp=zp, t2coh=false);
    end

    ## crop 2D data and increase dw
    spectra2d = [crop2d(spectra2d[i],1;w_max=4,step=1) for i = 1:length(T)]

    ## simulate effect of laser spectrum
    laserSpec =  .3 * exp.(-(spectra2d[1].ω.-1.35).^2/(2*(.19)^2)) +
                  exp.(-(spectra2d[1].ω.-1.6).^2/(2*(.11)^2)) +
                   .3 * exp.(-(spectra2d[1].ω.-1.85).^2/(2*(.05)^2))

    ## plot laser spectrum
    #figure()
    #plot(spectra2d[1].ω,laserSpec)

    ## 1D convoltion (in ω_excitation)
    IRF = (laserSpec.^2*(ones(length(laserSpec)))')
    ## 2D convolution
    IRF = (laserSpec.^2*(laserSpec.^2)')

    #for i in 1:length(out2d)
        conv2d = spectra2d[1].full2d .* IRF
    #end
end

if calc_2d

    ## plot 2D spectra for each(?) T
    # what to plot
    rep  = "absorptive"
    scal = "lin"

    ## make  subplot layout
    nplots = length(T);
    ncols  = Int32(ceil(sqrt(nplots)));
    nrows  = Int32(ceil(nplots / ncols));

    # determine maximum value in dataset out2d[:].full2d[:,:]
    maxi = maximum([maximum(real(spectra2d[i].full2d)) for i in 1:length(spectra2d)])

    # plot 2D spectra
    fig, ax = subplots(nrows,ncols,sharex=true,sharey=true,figsize=(ncols*3.2,nrows*3))
    if nplots == 1
        ax = [ax]
    end

    fig.suptitle(rep * " 2D spectrum (" * scal * ". scaling)")
    k = 0
    for i = 1:ncols
        for j = 1:nrows
            global k += 1
            if k > nplots
                continue
            end
            sca(ax[i,j])
            ax[i,j].set_aspect="equal"
            plot2d(spectra2d[k].ω,round.(spectra2d[k].full2d,digits=1);repr=rep,scaling=scal,norm=maxi)
            #plot2d(spectra2d[k].ω,conv2d;repr=rep,scaling=scal,norm=0)
            title("2D spectrum at $(T[k]) fs")
        end
    end
    tight_layout()
    subplots_adjust(top=0.8)

    ω = spectra2d[1].ω

    #DELETE or improve
    #plot([Ee-Eg, Ee-Eg] .+ 0*ωg,[ω[1], ω[end]],"k--",linewidth=0.5)
    #plot([Ee-Eg, Ee-Eg] .+ 2*ωg,[ω[1], ω[end]],"k--",linewidth=0.5)
    #plot([Ee-Eg, Ee-Eg] .+ 4*ωg,[ω[1], ω[end]],"k--",linewidth=0.5)
    #plot([ω[1], ω[end]],[Ee-Eg, Ee-Eg] .+ 0*ωg,"k--",linewidth=0.5)
    #plot([ω[1], ω[end]],[Ee-Eg, Ee-Eg] .+ 2*ωg,"k--",linewidth=0.5)
    #plot([ω[1], ω[end]],[Ee-Eg, Ee-Eg] .+ 4*ωg,"k--",linewidth=0.5)

end

max_stars_repo_name: lbmnky/CMDS

max_issues_repo_name: lbmnky/CMDS

max_forks_repo_name: lbmnky/CMDS


text:
__precompile__()
module TrackedDistributions

using Reexport

export
    TDiagonalNormal,
    TMVDiagonalNormal,
    kl_q_p,
    data,
    logσ

using LinearAlgebra
using Random
@reexport using Distributions

include("./distributions.jl")

end


text:
tests = ["lreg",
         "whiten",
         "pca",
         "cca",
         "cmds",
         "lda",
         "mclda",
         "ica",
         "ppca",
         "kpca",
         "fa",
         "facrot"]

for test in tests
    include(test*".jl")
end


text:
using jInv.Mesh;
using Helmholtz

const plotting = false;
if plotting
	using PyPlot;
	close("all")
end

# m = readdlm("SEGmodel2Dsalt.dat"); m = m';
m = ones(257,129);
m = m*1e-3;
Minv = getRegularMesh([0.0,13.5,0.0,4.2],collect(size(m)) .- 1);

pad = 25;

m = ones(size(m));
m = 1.0./m.^2

if plotting
	figure();
	imshow(1.0./sqrt.(Matrix(m')))
end
f = 1.0;
w = 2*pi*f

maxOmega = getMaximalFrequency(m,Minv);
ABLamp = maxOmega;

println("omega*h:");
println(w*Minv.h*sqrt(maximum(m)));
pad = pad*ones(Int64,Minv.dim);
H = GetHelmholtzOperator(Minv,m,w,ones(size(m))*0.01,true,pad,ABLamp,true)[1];
SH = H + GetHelmholtzShiftOP(m, w,0.1);

n_nodes = Minv.n.+1;
nnodes_tup = tuple(n_nodes...);

q,src = getAcousticPointSource(Minv,ComplexF64);

s = H\q[:];
s = real(reshape(s,nnodes_tup));
s[loc2cs(n_nodes,src)] = 0.0;

if plotting
	figure();
	imshow(s');title("Helmholtz Solution");
end

sh = SH\q[:];
sh = real(reshape(sh,nnodes_tup));
sh[loc2cs(n_nodes,src)]= 0.0;

if plotting
	figure();
	imshow(sh');title("Shifted Helmholtz Solution");
end


text:
using Documenter, DiffusionMCMC

makedocs(;
    modules=[DiffusionMCMC],
    format=Documenter.HTML(
        mathengine = Documenter.MathJax(
            Dict(
                :TeX => Dict(
                    :equationNumbers => Dict(
                        :autoNumber => "AMS"
                    ),
                    :Macros => Dict(
                        :dd => "{\\textrm d}",
                        :RR => "\\mathbb{R}",
                        :wt => ["\\widetilde{#1}", 1]
                    ),
                ),
            ),
        ),
        collapselevel = 1,
    ),
    pages=[
        "Home" => "index.md",
        "Get started" => joinpath("get_started", "overview.md"),
        "User manual" => Any[
            "Extension of ExtensibleMCMC.jl" => joinpath("manual", "how_extend.md"),
            "Updates & Decorators" => joinpath("manual", "updates_and_decorators.md"),
            "(TODO) Blocking" => joinpath("manual", "blocking.md"),
            "(TODO) Callbacks" => joinpath("manual", "callbacks.md"),
            "Diffusion Workspaces" => joinpath("manual", "workspaces.md"),
            "Conjugate updates" => joinpath("manual", "conjugate_updates.md"),
        ],
        "How to..." => Any[
            "(TODO) Do smoothing" => joinpath("how_to_guides", "smoothing.md"),
            "(TODO) Do blocking" => joinpath("how_to_guides", "blocking.md"),
        ],
        "Tutorials" => Any[
            "Inference & smoothing for the FitzHugh–Nagumo model" => joinpath("tutorials", "simple_inference.md"),
            "Conjugate updates for the FitzHugh–Nagumo model" => joinpath("tutorials", "conjugate_updates.md"),
            "(TODO) Inference from first-passage time observations" => joinpath("tutorials", "inference_fpt.md"),
            "(TODO) Mixed-effect models" => joinpath("tutorials", "mixed_effects_models.md"),
        ],
        "Index" => "module_index.md",
    ],
    repo="https://github.com/JuliaDiffusionBayes/DiffusionMCMC.jl/blob/{commit}{path}#L{line}",
    sitename="DiffusionMCMC.jl",
    authors="Sebastiano Grazzi, Frank van der Meulen, Marcin Mider, Moritz Schauer",
    assets=String[],
)

deploydocs(;
    repo="github.com/JuliaDiffusionBayes/DiffusionMCMC.jl",
)


text:
module StochasticSemiDiscretizationMethod

import InteractiveUtils
using Reexport
@reexport using LinearAlgebra
@reexport using SparseArrays
@reexport using StaticArrays
@reexport using Arpack
using QuadGK
using Lazy: iterated, take

import SemiDiscretizationMethod
import SemiDiscretizationMethod:
AbstractLDDEProblem, AbstractResult,
MatrixOrFunction, ArrayOrFunction, VectorOrFunction, RealOrFunction, CyclicVector,
DiscretizationMethod, SemiDiscretization, NumericSD, methodorder, lagr_el0,
Coefficients, CoefficientMatrix, AdditiveVector,
ProportionalMX, DelayMX, Additive,
Delay,
subArray, SubMX, SubV,
calculate_Aavgs,
# calculateResults,
addSubmatrixToResult!,addSubvectorToResults!,
DiscreteMapping, DiscreteMappingSteps,
subMxRange, rOfDelay, nStepOfLength, prodl,
reduce_additive

calculateDetResults! = SemiDiscretizationMethod.calculateResults!

include("structures_input.jl")
include("structures_method.jl")
include("structures_result.jl")

include("functions_method.jl")
include("functions_stoch_utilities.jl")
include("functions_discretization.jl")

export  SemiDiscretization, NumericSD,
ProportionalMX,
Delay,DelayMX,
stCoeffMX,
Additive, stAdditive,
LDDEProblem,
DiscreteMapping_M1, DiscreteMapping_M2,
MxToCovVec, VecToCovMx,
# DiscreteMapping_M1_1step, DiscreteMapping_M2_1step,
fixPointOfMapping, spectralRadiusOfMapping

end # module


text:
is_call(ex::Expr) = Meta.isexpr(ex, :call)
is_call(ex::Expr, op::Symbol) = Meta.isexpr(ex, :call) && ex.args[1] == op
is_call(::Any) = false
is_call(::Any, ::Any) = false

"""
    is_formula(ex::Expr)

Detect formula expressions (either two-sided or one-sided).

"""
is_formula(ex::Expr) = is_call(ex, :~) || is_call(ex) && is_call(ex.args[2], :~)
is_formula(::Any) = false


text:
import Base: mapreduce, sum
import Statistics: mean

@eval begin
    @explicit_intercepts(
        mapreduce,
        Tuple{Function, $plustype, AbstractArray{<:∇Scalar}},
        [false, false, true],
        (dims=:,),
    )
    function ∇(
        ::typeof(mapreduce),
        ::Type{Arg{3}},
        p, y, ȳ, f,
        ::$plustype,
        A::AbstractArray{<:∇Scalar};
        dims=:,
    )
        return broadcast((An, ȳn)->ȳn * ForwardDiff.derivative(f, An), A, ȳ)
    end
end

@explicit_intercepts(
    sum,
    Tuple{Function, AbstractArray{<:∇Scalar}},
    [false, true],
    (dims=:,),
)
function ∇(
    ::typeof(sum),
    ::Type{Arg{2}},
    p, y, ȳ,
    f::Function,
    A::AbstractArray{<:∇Scalar};
    dims=:,
)
    # Just pass through to mapreduce
    return ∇(mapreduce, Arg{3}, p, y, ȳ, f, Base.add_sum, A; dims=dims)
end
# sum(abs2, xs) is in ChainRules, but it results in method ambiguties with the
# version that accepts any function above
function ∇(
    ::typeof(sum),
    ::Type{Arg{2}},
    p, y, ȳ,
    ::typeof(abs2),
    A::AbstractArray{<:Real};
    dims=:,
)
    return 2ȳ .* A
end

@explicit_intercepts(
    mean,
    Tuple{Function, AbstractArray{<:∇Scalar}},
    [false, true],
    #(dims=:,)  # https://github.com/JuliaLang/julia/issues/31412
)

_denom(x, dims::Colon) = length(x)
_denom(x, dims::Integer) = size(x, dims)
_denom(x, dims) = mapreduce(i->size(x, i), Base.mul_prod, unique(dims), init=1)

function ∇(
    ::typeof(mean),
    ::Type{Arg{2}},
    p, y, ȳ,
    f::Function,
    x::AbstractArray{<:∇Scalar},
)
    return ∇(sum, Arg{2}, p, y, ȳ, f, x; dims=:) / _denom(x, :)
end


text:
facts("ss2float") do

  context("parses as expected") do
      @fact isempty("") => true
  end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ f1fb47b0-63ff-11ec-350c-7d4755d75c7b
begin
	using Pkg; Pkg.activate(@__DIR__); Pkg.instantiate()
	Pkg.precompile()

	using PlutoUI
	using Clustering
	using VegaLite
	using VegaDatasets
	using DataFrames
	using Statistics
	using JSON
	using CSV
	using Distances
end;

# ╔═╡ cd0d79a3-2e64-4991-ade0-dac7422e323c
PlutoUI.TableOfContents(aside=true, indent=true, depth=3)

# ╔═╡ 55d15235-24d0-41c8-85e6-67f807a3b759
html"""

"""

# ╔═╡ 0ea01882-d01b-477f-bbb1-fc4f3042cc50
md"""
# 🗃 Clustering

**Cluster analysis** or **clustering** is the task of grouping a set of objects in such a way that objects in the same group (called a cluster) are more similar (in some sense) to each other than to those in other groups (clusters). It is a main task of exploratory data analysis, and a common technique for statistical data analysis, used in many fields.

Put simply, the task of clustering is to place observations that seem similar (in _feature space_) within the same cluster. Clustering is commonly used in two dimensional data where the goal is to create clusters based on coordinates. Here, we will use something similar. We will cluster houses based on their latitude-longitude locations using several different clustering methods.

[YouTube Link](https://www.youtube.com/watch?v=cwurgt7cn5s)
"""

# ╔═╡ fbfaeaf4-8e40-4c4e-b618-01a0bb16f9b9
md"""
## Load data

We will start off by getting some data. We will use data of 20,000+ California houses dataset. We will then learn whether house prices directly correlate with map location.
"""

# ╔═╡ 678ccf3f-aee0-45e9-b410-1efe156e2ef5
begin
	download("https://raw.githubusercontent.com/alexhegit/handson-ml2/master/datasets/housing/housing.csv", "newhouses.csv")
	houses = CSV.read("newhouses.csv", DataFrame)
end

# ╔═╡ d6fe53b6-a8bc-4bb7-aa2f-8e53bead237a
md"""
Let's check the available variables...
"""

# ╔═╡ 298b9aad-ed34-4cba-a74c-d103ffe1e347
names(houses)

# ╔═╡ 17d0f736-1faf-40e6-aa6a-17ef2e74bdb4
md"""
## Plot maps

We will use the [VegaLite.jl](https://github.com/queryverse/VegaLite.jl) package here for plotting. This package makes it very easy to plot information on a map. All you need is a JSON file of the map you intend to draw. Here, we will use the California counties JSON file and plot each house on the map and color code it via a heatmap of the price. This is done by the following line:

```julia
color="median_house_value:q"
```
"""

# ╔═╡ 6a3627e2-9995-4846-b11b-7631b5a2501e
begin
	# importing California counties shape
	cali_shape = JSON.parsefile("data/california-counties.json")
	VV = VegaDatasets.VegaJSONDataset(cali_shape, "california-counties.json")
end

# ╔═╡ 2f1c19e5-3719-4bc7-ab7a-80f1f6c53fff
begin
	# figure layout
	@vlplot(width=500, height=300) +
	# california counties shape
	@vlplot(mark={:geoshape,
				 fill=:black,
				 stroke=:white},
		   data={values=VV,
				 format={type=:topojson,
				 		 feature=:cb_2015_california_county_20m}},
		   projection={type=:albersUsa},) +
	# houses colored by values
	@vlplot(:circle,
			data=houses,
			projection={type=:albersUsa},
			longitude="longitude:q",
			latitude="latitude:q",
			size={value=12},
			color="median_house_value:q")
end

# ╔═╡ df1ab399-d07b-462f-ae51-598c2c9f7e05
md"""
One thing we will try and explore in this notebook is if clustering the houses has any direct relationship with their prices, so we will bucket the houses into intervals of $50,000 and re-perform the color codes based on each bucket.
"""

# ╔═╡ b6de785c-58c4-4c17-be08-f269791c4316
begin
	# insert bucket price column
	bucketprice = Int.(div.(houses[!,:median_house_value], 50000))
	insertcols!(houses, 3, :cprice => bucketprice)
end

# ╔═╡ 513e85d6-3480-4502-a6a8-cd7c9ed0077f
begin
	# figure layout
	@vlplot(width=500, height=300) +
	# california counties shape
	@vlplot(mark={:geoshape,
				 fill=:black,
				 stroke=:white},
		   data={values=VV,
				 format={type=:topojson,
				 		 feature=:cb_2015_california_county_20m}},
		   projection={type=:albersUsa},) +
	# houses colored by bucket price
	@vlplot(:circle,
			data=houses,
			projection={type=:albersUsa},
			longitude="longitude:q",
			latitude="latitude:q",
			size={value=12},
			color="cprice:n")
end

# ╔═╡ d2477d5e-ccc1-4e91-aec7-d6446961db1b
md"""
## k-Means Clustering

**k-Means Clustering** is a method of vector quantization, originally from signal processing, that aims to partition n observations into k clusters in which each observation belongs to the cluster with the nearest mean (cluster centers or cluster centroid), serving as a prototype of the cluster. This results in a partitioning of the data space into Voronoi cells.

k-Means Clustering minimizes within-cluster variances (squared Euclidean distances), but not regular Euclidean distances, which would be the more difficult Weber problem: the mean optimizes squared errors, whereas only the geometric median minimizes Euclidean distances. For instance, better Euclidean solutions can be found using k-medians and k-medoids.

**Note:** for more math details about k-Means, check this [link](https://en.wikipedia.org/wiki/K-means_clustering).

First, we will apply k-Means to cluster coordinates into 10 cluesters and then we will insert a new column `cluster_means`.

**Note:** k-Means requires that observations must be represented as columns (not rows).
"""

# ╔═╡ c8e16e6c-b99e-405d-b0e4-0148e09cc818
begin
	# DataFrame of coordinates
	X = houses[!, [:latitude, :longitude]]
	# k-means model
	KM = kmeans(Matrix(X)', 10)
end

# ╔═╡ ece5be38-b8ba-4615-8755-6c05c9e6e77f
insertcols!(houses, 3, :cluster_means => KM.assignments)

# ╔═╡ d8171788-16dd-4f0b-95da-86b7a35fab92
md"""
Now we can visualize the data colored by the new column `cluster_means`...
"""

# ╔═╡ 0073e46e-0350-4f9c-975c-31355d9d3040
begin
	# figure layout
	@vlplot(width=500, height=300) +
	# california counties shape
	@vlplot(mark={:geoshape,
				 fill=:black,
				 stroke=:white},
		   data={values=VV,
				 format={type=:topojson,
				 		 feature=:cb_2015_california_county_20m}},
		   projection={type=:albersUsa},) +
	# houses colored by cluster_means column
	@vlplot(:circle,
			data=houses,
			projection={type=:albersUsa},
			longitude="longitude:q",
			latitude="latitude:q",
			size={value=12},
			color="cluster_means:n")
end

# ╔═╡ 65ece652-2caa-4c9f-a607-b19a7087a065
md"""
As we can see, location indeed affects price of house. we may want to know if this relationship remains true for location as in proximity to water, downtown or bus stops.
"""

# ╔═╡ f9ab7250-a2db-4824-9fda-b68f3fcb582e
md"""
## k-Medoids Clustering

The **k-Medoids** problem is a clustering problem similar to k-Means. The name was coined by Leonard Kaufman and Peter J. Rousseeuw with their PAM algorithm. Both the k-Means and k-Medoids algorithms are partitional (breaking the dataset up into groups) and attempt to minimize the distance between points labeled to be in a cluster and a point designated as the center of that cluster.

In contrast to the k-Means algorithm, k-Medoids chooses actual data points as centers (medoids or exemplars), and thereby allows for greater interpretability of the cluster centers than in k-Means, where the center of a cluster is not necessarily one of the input data points (it is the average between the points in the cluster).

Furthermore, k-Medoids can be used with arbitrary dissimilarity measures, whereas k-Means generally requires Euclidean distance for efficient solutions. Because k-Medoids minimizes a sum of pairwise dissimilarities instead of a sum of squared Euclidean distances, it is more robust to noise and outliers than k-Means.

**Note:** for more math details about k-Medoids, check this [link](https://en.wikipedia.org/wiki/K-medoids).

For k-Medoids, we need to build a distance matrix as we did in the last notebook (UMAP for dimensionality reduction). We will use the [Distances.jl](https://github.com/JuliaStats/Distances.jl) package for this purpose and compute the pairwise Euclidean distances.
"""

# ╔═╡ 2d933b98-83c3-4f9e-9858-cd832cb6b457
begin
	# transpose matrix
	Xᵗ = Matrix(X)'

	# calculating Euclidean distance matrix
	D = pairwise(Euclidean(), Xᵗ, Xᵗ, dims=2)
end

# ╔═╡ 6a42ec90-b400-4638-91aa-4b982c2da069
# k-Medoids model
KMD = kmedoids(D, 10)

# ╔═╡ 8f55e8e6-37b2-4d1d-bb0c-db2c83f854be
# inserting k-Medoids column
insertcols!(houses, 3, :cluster_medoids => KMD.assignments)

# ╔═╡ 1dd518f3-3c63-4610-9736-7369f980392a
md"""
Now we can visualize the data colored by the new column `cluster_medoids`...
"""

# ╔═╡ 0ea39069-dbbf-4bd9-becb-560cdbeb9e8e
begin
	# figure layout
	@vlplot(width=500, height=300) +
	# california counties shape
	@vlplot(mark={:geoshape,
				 fill=:black,
				 stroke=:white},
		   data={values=VV,
				 format={type=:topojson,
				 		 feature=:cb_2015_california_county_20m}},
		   projection={type=:albersUsa},) +
	# houses colored by cluster_medoids column
	@vlplot(:circle,
			data=houses,
			projection={type=:albersUsa},
			longitude="longitude:q",
			latitude="latitude:q",
			size={value=12},
			color="cluster_medoids:n")
end

# ╔═╡ 33eed20d-9961-4e96-9997-e2866c20e097
md"""
## Hierarchical Clustering

**Hierarchical Clustering (HCA)** is a method of cluster analysis which seeks to build a hierarchy of clusters. Strategies for hierarchical clustering generally fall into two types:

- **Agglomerative**: This is a "bottom-up" approach: each observation starts in its own cluster, and pairs of clusters are merged as one moves up the hierarchy.

- **Divisive**: This is a "top-down" approach: all observations start in one cluster, and splits are performed recursively as one moves down the hierarchy.

In general, the merges and splits are determined in a greedy manner. The results of Hierarchical Clustering are usually presented in a dendrogram.

**Note:** for more math details about Hierarchical Clustering, check this [link](https://en.wikipedia.org/wiki/Hierarchical_clustering).

First, we will input the distance matrix `D` to `hclust` function and then generate 10 clusters using the `cutree` function. Finally, we will add the new `cluster_hca` column to our data.
"""

# ╔═╡ 177f94ef-6b45-4324-9b8a-6b07dc853e54
# HCA
HCA = hclust(D)

# ╔═╡ a5a2a363-9ae7-44d8-89cc-3aaaf7aa3d9e
# generating clusters
hca_clusters = cutree(HCA; k=10)

# ╔═╡ 56f2d129-6a8a-4750-b292-1682f7387165
# inserting k-Medoids column
insertcols!(houses, 3, :cluster_hca => hca_clusters)

# ╔═╡ bdc23b6c-b686-4971-984d-63c2154e661e
md"""
Now we can visualize the data colored by the new column `cluster_hca`...
"""

# ╔═╡ f4fc1f80-601f-4075-8c6d-16cb6e9f0223
begin
	# figure layout
	@vlplot(width=500, height=300) +
	# california counties shape
	@vlplot(mark={:geoshape,
				 fill=:black,
				 stroke=:white},
		   data={values=VV,
				 format={type=:topojson,
				 		 feature=:cb_2015_california_county_20m}},
		   projection={type=:albersUsa},) +
	# houses colored by cluster_hca column
	@vlplot(:circle,
			data=houses,
			projection={type=:albersUsa},
			longitude="longitude:q",
			latitude="latitude:q",
			size={value=12},
			color="cluster_hca:n")
end

# ╔═╡ b6091ca6-873c-4cc2-8ee2-30d884db514f
md"""
## DBSCAN

**Density-based Spatial Clustering of Applications with Noise (DBSCAN)** is a data clustering algorithm proposed by Martin Ester, Hans-Peter Kriegel, Jörg Sander and Xiaowei Xu in 1996. It is a density-based clustering non-parametric algorithm: given a set of points in some space, it groups together points that are closely packed together (points with many nearby neighbors), marking as outliers points that lie alone in low-density regions (whose nearest neighbors are too far away).

DBSCAN is one of the most common clustering algorithms and also most cited in scientific literature.

**Note:** for more math details about DBSCAN, check this [link](https://en.wikipedia.org/wiki/DBSCAN).

We just need to pass the distance matrix `D` to perform DBSCAN. We may also set both `eps` (the radius of a point neighborhood) and `minpts` (the minimum number of neighboring points to qualify a point as a density point) arguments.

**Note:** we cannot pass the number of clusters when using DBSCAN.
"""

# ╔═╡ 35702907-ecec-4441-a6c6-3497b02bfd51
# creating a new distance matrix (Square of Euclidean distance)
D² = pairwise(SqEuclidean(), Xᵗ, dims=2)

# ╔═╡ 431648b8-406b-4064-b42f-1bb7c0721026
# DBSCAN
DBS = dbscan(D², 0.05, 10)

# ╔═╡ 5f4aba99-22c4-43b9-be6e-ccc45f349158
# inserting DBSCAN column
insertcols!(houses, 3, :cluster_dbscan => DBS.assignments)

# ╔═╡ 7db3222d-eab6-4ddd-8d5e-7990f437a85e
md"""
Now we can visualize the data colored by the new column `cluster_dbscan`...
"""

# ╔═╡ b632e115-6278-4382-9d1b-1d708d40ac9a
begin
	# figure layout
	@vlplot(width=500, height=300) +
	# california counties shape
	@vlplot(mark={:geoshape,
				 fill=:black,
				 stroke=:white},
		   data={values=VV,
				 format={type=:topojson,
				 		 feature=:cb_2015_california_county_20m}},
		   projection={type=:albersUsa},) +
	# houses colored by cluster_dbscan column
	@vlplot(:circle,
			data=houses,
			projection={type=:albersUsa},
			longitude="longitude:q",
			latitude="latitude:q",
			size={value=12},
			color="cluster_dbscan:n")
end

# ╔═╡ Cell order:
# ╟─f1fb47b0-63ff-11ec-350c-7d4755d75c7b
# ╟─cd0d79a3-2e64-4991-ade0-dac7422e323c
# ╟─55d15235-24d0-41c8-85e6-67f807a3b759
# ╟─0ea01882-d01b-477f-bbb1-fc4f3042cc50
# ╟─fbfaeaf4-8e40-4c4e-b618-01a0bb16f9b9
# ╠═678ccf3f-aee0-45e9-b410-1efe156e2ef5
# ╟─d6fe53b6-a8bc-4bb7-aa2f-8e53bead237a
# ╠═298b9aad-ed34-4cba-a74c-d103ffe1e347
# ╟─17d0f736-1faf-40e6-aa6a-17ef2e74bdb4
# ╠═6a3627e2-9995-4846-b11b-7631b5a2501e
# ╠═2f1c19e5-3719-4bc7-ab7a-80f1f6c53fff
# ╟─df1ab399-d07b-462f-ae51-598c2c9f7e05
# ╠═b6de785c-58c4-4c17-be08-f269791c4316
# ╠═513e85d6-3480-4502-a6a8-cd7c9ed0077f
# ╟─d2477d5e-ccc1-4e91-aec7-d6446961db1b
# ╠═c8e16e6c-b99e-405d-b0e4-0148e09cc818
# ╠═ece5be38-b8ba-4615-8755-6c05c9e6e77f
# ╟─d8171788-16dd-4f0b-95da-86b7a35fab92
# ╠═0073e46e-0350-4f9c-975c-31355d9d3040
# ╟─65ece652-2caa-4c9f-a607-b19a7087a065
# ╟─f9ab7250-a2db-4824-9fda-b68f3fcb582e
# ╠═2d933b98-83c3-4f9e-9858-cd832cb6b457
# ╠═6a42ec90-b400-4638-91aa-4b982c2da069
# ╠═8f55e8e6-37b2-4d1d-bb0c-db2c83f854be
# ╟─1dd518f3-3c63-4610-9736-7369f980392a
# ╠═0ea39069-dbbf-4bd9-becb-560cdbeb9e8e
# ╟─33eed20d-9961-4e96-9997-e2866c20e097
# ╠═177f94ef-6b45-4324-9b8a-6b07dc853e54
# ╠═a5a2a363-9ae7-44d8-89cc-3aaaf7aa3d9e
# ╠═56f2d129-6a8a-4750-b292-1682f7387165
# ╟─bdc23b6c-b686-4971-984d-63c2154e661e
# ╠═f4fc1f80-601f-4075-8c6d-16cb6e9f0223
# ╟─b6091ca6-873c-4cc2-8ee2-30d884db514f
# ╠═35702907-ecec-4441-a6c6-3497b02bfd51
# ╠═431648b8-406b-4064-b42f-1bb7c0721026
# ╠═5f4aba99-22c4-43b9-be6e-ccc45f349158
# ╟─7db3222d-eab6-4ddd-8d5e-7990f437a85e
# ╠═b632e115-6278-4382-9d1b-1d708d40ac9a

max_stars_repo_name: fnaghetini/DataScience

max_issues_repo_name: fnaghetini/DataScience

max_forks_repo_name: fnaghetini/DataScience


text:
#####
##### Tuning MCTS parameters
#####

using AlphaZero
using AlphaZero: Log, Util
using Formatting
using Base: @kwdef

const NUM_GAMES = 1000
const RESET_MCTS_EVERY = 10
const LOGGER = Log.Logger()

#####
##### MCTS parameters to tune
#####

@kwdef struct Config
  n :: Int      # num simulations
  w :: Int      # num workers
  c :: Float64  # cpuct
  τ :: Float64
  ϵ :: Float64
  α :: Float64
  p :: Float64 # prior temperature
end

Util.generate_update_constructor(Config) |> eval

function Base.string(conf::Config)
  c = fmt(".1f", conf.c)
  τ = fmt(".2f", conf.τ)
  ϵ = fmt(".2f", conf.ϵ)
  α = fmt(".2f", conf.α)
  p = fmt(".2f", conf.p)
  return "$(conf.n) simulations, $(conf.w) workers; " *
    "c=$c, τ=$τ, ϵ=$ϵ, α=$α, p=$p"
end

#####
##### Experiment Report
#####

@kwdef struct Report
  reward :: Float64
  time_per_game :: Float64
  inference_time_ratio :: Float64
  average_exploration_depth :: Float64
end

function print_report(report, legend)
  wr   = Log.ColType(10, x -> format("{:.2f}", x))
  time = Log.ColType(10, x -> format("{:.1f} min", x * 100 / 60))
  itr  = Log.ColType(5,  x -> format("{}%", round(Int, x * 100)))
  expd = Log.ColType(5,  x -> format("{:.1f}", x))
  tab  = Log.Table([
    ("AVGR", wr,   r -> r.reward),
    ("T100", time, r -> r.time_per_game),
    ("ITR",  itr,  r -> r.inference_time_ratio),
    ("EXPD", expd, r -> r.average_exploration_depth)])
  Log.table_row(LOGGER, tab, report, [legend])
end

#####
##### Measurement Code
#####

using Juno

function evaluate_player(network, params, baseline, n)
  player = AlphaZero.MctsPlayer(network, params)
  Rtot = 0.0
  time = @elapsed Juno.@progress for i in 1:n
    if i % RESET_MCTS_EVERY == 0
      MCTS.reset!(player.mcts)
    end
    Rtot += AlphaZero.play_game(player, baseline)
  end
  return Report(
    reward=(Rtot / n),
    time_per_game=(time / n),
    inference_time_ratio=MCTS.inference_time_ratio(player.mcts),
    average_exploration_depth=MCTS.average_exploration_depth(player.mcts))
end

function run_experiments(network, baseline, configs)
  for config in configs
    params = MctsParams(
      use_gpu=true,
      num_workers=config.w,
      num_iters_per_turn=config.n,
      cpuct=config.c,
      temperature=ConstSchedule(config.τ),
      dirichlet_noise_ϵ=config.ϵ,
      dirichlet_noise_α=config.α,
      prior_temperature=config.p)
    evaluate_player(network, params, baseline, 1) # Compilation
    rep = evaluate_player(network, params, baseline, NUM_GAMES)
    print_report(rep, string(config))
  end
end

#####
##### Experiments
#####

include("../games.jl")
const GAME = "connect-four"
const SelectedGame = GAME_MODULE[GAME]
using .SelectedGame: Game, Training

const SESSION_DIR = "sessions/connect-four"

session = Session(
  Game, Training.Network{Game},
  Training.params, Training.netparams, benchmark=Training.benchmark,
  dir=SESSION_DIR, autosave=false, save_intermediate=false)

configs = [
  Config(n=n, w=w, c=c, τ=τ, ϵ=ϵ, α=α, p=p)
  for τ in [1.0, 0.5, 0.1]
  for n in [400, 600]
  for w in [8, 16, 32]
  for ϵ in [0.1, 0.2]
  for c in [3.0]
  for α in [1.0]
  for p in [1.0]
]

network = session.env.bestnn
baseline = MinMax.Player{Game}(;depth=5, τ=0.2)
run_experiments(network, baseline, configs)


text:
@testset "Duality-Based Bound Tightening" begin

    ylower = Float64[1.0, 1.0, 1.0, 1.0]
    yupper = Float64[4.0, 4.0, 4.0, 4.0]
    ymult_lo = Float64[50.0, 0.0, 1.0, 0.0]
    ymult_hi = Float64[0.0, 0.0, 0.8, 3.0]

    n = EAGO.NodeBB(ylower, yupper, -Inf, Inf, 3, 2)
    @inferred EAGO.variable_dbbt!(n, ymult_lo, ymult_hi, 1.0, 3.0, 4)
    lvb = n.lower_variable_bounds
    uvb = n.upper_variable_bounds

    @test lvb[1] == 1.0; @test uvb[1] == 1.04
    @test lvb[2] == 1.0; @test uvb[2] == 4.0
    @test lvb[3] == 1.0; @test uvb[3] == 3.0
    # lvb[3] doesn't tighten since dbbt assumes variables aren't fixed...
    @test isapprox(lvb[4], 3.33333333, atol= 1E-5); @test uvb[4] == 4.0
end

#=
@testset "Classify Quadratic Types" begin

    @test EAGO.get_value(MOI.LessThan{Float64}(1.1)) == 1.1
    @test EAGO.get_value(MOI.GreaterThan{Float64}(2.1)) == 2.1
    @test EAGO.get_value(MOI.EqualTo{Float64}(1.3)) == 1.3

    opt1 = EAGO.Optimizer()
    opt1._parameters.verbosity = 0
    x = MOI.add_variables(opt1, 3)

    MOI.add_constraint(opt1,MOI.SingleVariable(x[1]), MOI.LessThan(4.0))
    MOI.add_constraint(opt1,MOI.SingleVariable(x[1]), MOI.GreaterThan(-2.0))
    MOI.add_constraint(opt1,MOI.SingleVariable(x[2]), MOI.LessThan(4.0))
    MOI.add_constraint(opt1,MOI.SingleVariable(x[2]), MOI.GreaterThan(0.0))
    MOI.add_constraint(opt1,MOI.SingleVariable(x[3]), MOI.LessThan(1.0))
    MOI.add_constraint(opt1,MOI.SingleVariable(x[3]), MOI.GreaterThan(-1.0))

    # adds four univariate quadratic constraints
    qterm1a = [MOI.ScalarQuadraticTerm{Float64}(1.5, MOI.VariableIndex(1), MOI.VariableIndex(1))]
    qterm2a = [MOI.ScalarQuadraticTerm{Float64}(1.5, MOI.VariableIndex(1), MOI.VariableIndex(1))]
    qterm3a = [MOI.ScalarQuadraticTerm{Float64}(1.5, MOI.VariableIndex(2), MOI.VariableIndex(2))]
    qterm4a = [MOI.ScalarQuadraticTerm{Float64}(1.5, MOI.VariableIndex(2), MOI.VariableIndex(2))]
    MOI.add_constraint(opt1,MOI.ScalarQuadraticFunction{Float64}([MOI.ScalarAffineTerm(1.1, x[1])], qterm1a, 2.0), MOI.GreaterThan(2.0))
    MOI.add_constraint(opt1,MOI.ScalarQuadraticFunction{Float64}([MOI.ScalarAffineTerm(1.2, x[1])], qterm2a, 1.0), MOI.LessThan(0.0))
    MOI.add_constraint(opt1,MOI.ScalarQuadraticFunction{Float64}([MOI.ScalarAffineTerm(1.3, x[2])], qterm3a, -2.0), MOI.GreaterThan(2.0))
    MOI.add_constraint(opt1,MOI.ScalarQuadraticFunction{Float64}([MOI.ScalarAffineTerm(1.4, x[2])], qterm4a, -1.0), MOI.EqualTo(0.0))

    # adds four bivariate quadratic constraints
    qterm1b = [MOI.ScalarQuadraticTerm{Float64}(1.5, MOI.VariableIndex(1), MOI.VariableIndex(2))]
    qterm2b = [MOI.ScalarQuadraticTerm{Float64}(1.5, MOI.VariableIndex(2), MOI.VariableIndex(2))]
    qterm3b = [MOI.ScalarQuadraticTerm{Float64}(1.5, MOI.VariableIndex(1), MOI.VariableIndex(1))]
    qterm4b = [MOI.ScalarQuadraticTerm{Float64}(1.5, MOI.VariableIndex(1), MOI.VariableIndex(2))]
    MOI.add_constraint(opt1,MOI.ScalarQuadraticFunction{Float64}([MOI.ScalarAffineTerm(2.1, x[1])], qterm1b, -2.0), MOI.GreaterThan(2.0))
    MOI.add_constraint(opt1,MOI.ScalarQuadraticFunction{Float64}([MOI.ScalarAffineTerm(2.2, x[1])], qterm2b, -1.0), MOI.LessThan(0.0))
    MOI.add_constraint(opt1,MOI.ScalarQuadraticFunction{Float64}([MOI.ScalarAffineTerm(2.3, x[2])], qterm3b, -2.0), MOI.GreaterThan(2.0))
    MOI.add_constraint(opt1,MOI.ScalarQuadraticFunction{Float64}([MOI.ScalarAffineTerm(2.4, x[2])], qterm4b, -1.0), MOI.EqualTo(0.0))
end
=#

#=
@testset "Quadratic Domain Reduction (Univariate)" begin
    m = EAGO.Optimizer(verbosity = 0)

    m._current_node = NodeBB([-10.0, -10.0, -10.0], [10.0, 10.0, 10.0], -Inf, Inf, 1, 1)
    a1, b1, c1 = 2.0, -4.0, 3.0
    a2, b2, c2 = 2.0, -4.0, 3.0
    a3, b3, c3 = 2.0, -4.0, 3.0
    push!(m._univariate_quadratic_geq_constraints,(a1,b1,c1,1))
    push!(m._univariate_quadratic_leq_constraints,(a2,b2,c2,2))
    push!(m._univariate_quadratic_eq_constraints,(a3,b3,c3,3))
    feas = EAGO.univariate_quadratic(m)

    @test feas == true
    @test isapprox(m._current_node.lower_variable_bounds[1],-10.0,atol=1E-3)
    @test isapprox(m._current_node.lower_variable_bounds[2],-0.58113,atol=1E-3)
    @test isapprox(m._current_node.lower_variable_bounds[3],-0.58113,atol=1E-3)
    @test isapprox(m._current_node.upper_variable_bounds[1],10.0,atol=1E-3)
    @test isapprox(m._current_node.upper_variable_bounds[2],2.5811,atol=1E-3)
    @test isapprox(m._current_node.upper_variable_bounds[3],2.5811,atol=1E-3)
end
=#

#=
@testset "Quadratic Domain Reduction (Bivariate)" begin
end

@testset "Interval CSP" begin
    # Vigerske 2017 example
    opt1 = EAGO.Optimizer()
    x = MOI.add_variable(opt1)
    y = MOI.add_variable(opt1)
    z = MOI.add_variable(opt1)

    MOI.add_constraint(opt1,MOI.SingleVariable(x), MOI.LessThan(Inf))
    MOI.add_constraint(opt1,MOI.SingleVariable(x), MOI.GreaterThan(-Inf))
    MOI.add_constraint(opt1,MOI.SingleVariable(y), MOI.LessThan(Inf))
    MOI.add_constraint(opt1,MOI.SingleVariable(y), MOI.GreaterThan(-Inf))
    MOI.add_constraint(opt1,MOI.SingleVariable(z), MOI.LessThan(Inf))
    MOI.add_constraint(opt1,MOI.SingleVariable(z), MOI.GreaterThan(-Inf))

    MOI.add_constraint(opt1,MOI.ScalarAffineFunction{Float64}([MOI.ScalarAffineTerm(1.0, x[1]), MOI.ScalarAffineTerm(1.0, x[2])], 0.0), MOI.GreaterThan(4.0))
    MOI.add_constraint(opt1,MOI.ScalarAffineFunction{Float64}([MOI.ScalarAffineTerm(1.0, x[1]), MOI.ScalarAffineTerm(1.0, x[2])], 0.0), MOI.LessThan(1.0))

    n1 = NodeData()
    n1.lower_variable_bounds = [-Inf, -Inf, 0.0]
    n1.upper_variable_bounds = [Inf, Inf, 1.0]
    n1.lower_bound = -Inf
    n1.upper_bound = Inf
    n1.depth = 3

    # C1 = exp(x^2 + y^2) - z <= 0
    # C2 = z*sqrt(x^2+y^2) <= 1

    feas1 = EAGO.IntervalCSP(opt1,n1)
    @test (n1.lower_variable_bounds[1] == 0.0)
    @test (n1.lower_variable_bounds[2] == 0.0)
    @test (n1.lower_variable_bounds[3] == 1.0)
    @test (n1.upper_variable_bounds[1] == 0.0)
    @test (n1.upper_variable_bounds[2] == 0.0)
    @test (n1.upper_variable_bounds[3] == 1.0)
    @test feas1 == true
end
=#
#=
@testset "Optimization-Based Bound Tightening (Linear)" begin
    m = Model(optimizer_with_attributes(EAGO.Optimizer, "verbosity" => 0))
    xL = [-4.0; -2.0]
    xU = [4.0; 2.0]

    @variable(m, xL[i] <= x[i=1:2] <= xU[i])

    @objective(m, Min, x[2] - x[1])
    @constraint(m, o1, x[1] >= -1.0)
    @constraint(m, o2, x[2] <= 1.0)
    @constraint(m, e1, -0.1*x[1] + x[2] >= -1.6)
    @constraint(m, e2, -3.7*x[1] + x[2] >= -12.8)
    @constraint(m, e3, 0.1*x[1] - x[2] >= -1.6)
    @constraint(m, e4, 3.7*x[1] - x[2] >= -12.8)

    m.moi_backend.optimizer.model.optimizer._global_lower_bound = Inf
    m.moi_backend.optimizer.model.optimizer._global_upper_bound = -Inf

    optimize!(m)
    opt = m.moi_backend.optimizer.model.optimizer
    m.moi_backend.optimizer.model.optimizer._branch_variables[1] = true
    m.moi_backend.optimizer.model.optimizer._branch_variables[2] = true
    m.moi_backend.optimizer.model.optimizer._global_lower_bound = -12.0
    m.moi_backend.optimizer.model.optimizer._global_upper_bound = 0.0
    opt._current_node = NodeBB(xL, xU,-12.0,0.0,1,1)
    opt.obbt_variable_values = [true; true]
    feas = EAGO.obbt!(opt)

    @test feas
    @test isapprox(opt._current_node.lower_variable_bounds[1], -1.0, atol = 1E-6)
    @test isapprox(opt._current_node.lower_variable_bounds[2], -1.7, atol = 1E-6)
    @test isapprox(opt._current_node.upper_variable_bounds[1], 3.729729, atol = 1E-6)
    @test isapprox(opt._current_node.upper_variable_bounds[2], 1.0, atol = 1E-6)
end

@testset "Optimization-Based Bound Tightening (Nonlinear)" begin
    m = EAGO.Optimizer()
    m._parameters.verbosity = 0

    xL = [0.0; -2.0]
    xU = [4.0; 4.0]

    @variable(m, xL[i] <= x[i=1:2] <= xU[i])
    @objective(m, Min, x[2] - x[1])
    @NLconstraint(m, e1, exp(x[1]/3) - x[2] == 0)
    @NLconstraint(m, e3, x[2] <= 3*log(x[1]+4))

    m.moi_backend.optimizer.model.optimizer._global_lower_bound = Inf
    m.moi_backend.optimizer.model.optimizer._global_upper_bound = -Inf

    optimize!(m)
    opt = m.moi_backend.optimizer.model.optimizer
    m.moi_backend.optimizer.model.optimizer.branch_variable[1] = true
    m.moi_backend.optimizer.model.optimizer.branch_variable[2] = true
    m.moi_backend.optimizer.model.optimizer._global_lower_bound = -10.0
    m.moi_backend.optimizer.model.optimizer._global_upper_bound = 10.0
    oldn = pop!(opt._stack)
    opt._current_node = NodeBB(oldn.lower_variable_bounds, oldn.upper_variable_bounds,-10.0,10.0,1,1)
    opt.obbt_variable_values[1] = true
    opt.obbt_variable_values[2] = true
    feas = EAGO.obbt(opt)

    @test feas
    @test opt._current_node.lower_variable_bounds[1] == 0.0
    @test isapprox(opt._current_node.upper_variable_bounds[1], 4.0, atol = 1E-6)
    @test isapprox(opt._current_node.lower_variable_bounds[2], 0.6492446803515586, atol = 1E-6)
    @test isapprox(opt._current_node.upper_variable_bounds[2], 3.7936678946831783, atol = 1E-6)
end
=#


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.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
module A

function f(x)
    a = x^2
    eval(:(a + 1))
end

end

A.f(3)

# %%
module A1

function f(x)
    global a = x^2
    eval(:(a + 1))
end

end

A1.f(3)

# %%
module B

function f(x)
    a = x^2
    eval(:($a + 1))
end

end

B.f(3)

# %%
@code_warntype B.f(3)

# %%
module C

function f(x)
    a = x^2
    eval(Expr(:call, :+, a, 1))
end

end

C.f(3)

# %%
@code_warntype C.f(3)

# %%
module D

function f(x)
    a = x^2
    a + 1
end

end

D.f(3)

# %%
@code_warntype D.f(3)

# %%
z = 9
:($z + 1) == Expr(:call, :+, z, 1)

# %%

max_stars_repo_path: 0013/eval in function.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0013/eval in function.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0013/eval in function.jl

max_forks_repo_name: genkuroki/public


text:
module OceanStateEstimation

using Artifacts, LazyArtifacts, Downloads, Tar, CodecZlib
using Statistics, FortranFiles, MeshArrays, MITgcmTools

export dataverse_lists, get_from_dataverse
export get_ecco_variable_if_needed, get_ecco_velocity_if_needed
export get_occa_variable_if_needed, get_occa_velocity_if_needed
export ECCOclim_path, OCCAclim_path
export MITPROFclim_path, CBIOMESclim_path, ECCOdiags_path

##

p=dirname(pathof(OceanStateEstimation))

artifact_toml = joinpath(p, "../Artifacts.toml")

ECCOclim_hash = artifact_hash("ECCOclim", artifact_toml)
ECCOclim_path = artifact_path(ECCOclim_hash)*"/"

OCCAclim_hash = artifact_hash("OCCAclim", artifact_toml)
OCCAclim_path = artifact_path(OCCAclim_hash)*"/"

MITPROFclim_hash = artifact_hash("MITPROFclim", artifact_toml)
MITPROFclim_path = artifact_path(MITPROFclim_hash)*"/"

"""
    MITPROFclim_download()

Download lazy artifact to `MITPROFclim_path`.
"""
MITPROFclim_download() = artifact"MITPROFclim"

CBIOMESclim_hash = artifact_hash("CBIOMESclim", artifact_toml)
CBIOMESclim_path = artifact_path(CBIOMESclim_hash)*"/"

"""
    CBIOMESclim_download()

Download lazy artifact to `CBIOMESclim_path`.
"""
CBIOMESclim_download() = artifact"CBIOMESclim"

ECCOdiags_hash = artifact_hash("ECCOdiags", artifact_toml)
ECCOdiags_path = artifact_path(ECCOdiags_hash)*"/"

"""
    ECCOdiags_download()

Download lazy artifact to `ECCOdiags_path`.
"""
ECCOdiags_download() = artifact"ECCOdiags"

## Dataverse Donwloads

"""
    get_from_dataverse(lst::String,nam::String,pth::String)

```
using OceanStateEstimation, CSV, DataFrames
pth=dirname(pathof(OceanStateEstimation))
lst=joinpath(pth,"../examples/OCCA_climatology.csv")
nams=CSV.read(lst,DataFrame)
[get_from_dataverse(lst,nam,OCCAclim_path) for nam in nams.name[:]]
```
"""
function get_from_dataverse(lst::String,nam::String,pth::String)
    lists=dataverse_lists(lst)
    ii = findall([occursin("$nam", lists.name[i]) for i=1:length(lists.ID)])
    for i in ii
        nam1=Downloads.download(lists.URL[i])
        if length(ii)>1
            !isdir(joinpath(pth,nam)) ? mkdir(joinpath(pth,nam)) : nothing
            nam2=joinpath(pth,nam,lists.name[i])
            mv(nam1,nam2)
        else
            nam2=joinpath(pth,lists.name[i])
            mv(nam1,nam2)
        end
    end
end

get_from_dataverse(nam::String,pth::String) = get_from_dataverse("../examples/nctiles_climatology.csv",nam,pth)

"""
    dataverse_lists(lst::String)

Read and derive lists (ID,name,URL) from csv file (ID,name) and return as tuple
```
lists=dataverse_lists(lst)
```
"""
function dataverse_lists(lst::String)
    tmp=readlines(lst)
    ID=[parse(Int,tmp[j][1:findfirst(isequal(','),tmp[j])-1]) for j=2:length(tmp)]
    name=[tmp[j][findfirst(isequal(','),tmp[j])+1:end] for j=2:length(tmp)]
    tmp="https://dataverse.harvard.edu/api/access/datafile/"
    URL=[tmp*"$(ID[j])" for j=1:length(ID)]
    return (ID=ID,name=name,URL=URL)
end

"""
    get_ecco_files(γ::gcmgrid,v::String,t=1)

```
using MeshArrays, OceanStateEstimation
γ=GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
tmp=get_ecco_files(γ,"oceQnet")
```
"""
function get_ecco_files(γ::gcmgrid,v::String,t=1)
    get_ecco_variable_if_needed(v)
    #return read_nctiles(ECCOclim_path*"$v/$v","$v",γ,I=(:,:,:,t))
    return read_nctiles(ECCOclim_path*"$v/$v","$v",γ,I=(:,:,t))
end

"""
    get_ecco_variable_if_needed(v::String)

Download ECCO output for variable `v` to `ECCOclim_path` if needed
"""
function get_ecco_variable_if_needed(v::String)
    p=dirname(pathof(OceanStateEstimation))
    lst=joinpath(p,"../examples/nctiles_climatology.csv")
    pth=ECCOclim_path
    !isdir(pth*v) ? get_from_dataverse(lst,v,pth) : nothing
end

"""
    get_ecco_velocity_if_needed()

Download ECCO output for `u,v,w` to `ECCOclim_path` if needed
"""
function get_ecco_velocity_if_needed()
    get_ecco_variable_if_needed("UVELMASS")
    get_ecco_variable_if_needed("VVELMASS")
    get_ecco_variable_if_needed("WVELMASS")
end

"""
    get_occa_variable_if_needed(v::String)

Download OCCA output for variable `v` to `OCCAclim_path` if needed
"""
function get_occa_variable_if_needed(v::String)
    p=dirname(pathof(OceanStateEstimation))
    lst=joinpath(p,"../examples/OCCA_climatology.csv")
    fil=joinpath(OCCAclim_path,v*".0406clim.nc")
    !isfile(fil) ? get_from_dataverse(lst,v,OCCAclim_path) : nothing
end

"""
    get_occa_velocity_if_needed()

Download OCCA output for `u,v,w` to `OCCAclim_path` if needed
"""
function get_occa_velocity_if_needed()
    nams = ("DDuvel","DDvvel","DDwvel","DDtheta","DDsalt")
    [get_occa_variable_if_needed(nam) for nam in nams]
    "done"
end

## zenodo.org Downloads

"""
    ECCOdiags_add(nam::String)

Add data to the ECCOdiags_path folder. Known options for `nam` include
"release3", "release3", "release4", and "interp_coeffs". Note that
"release2" is the estimate that's readily donwloaded by ECCOdiags_download().
"""
function ECCOdiags_add(nam::String)
    if nam=="release1"
        url="https://zenodo.org/record/5787104/files/ECCOv4r1_analysis.tar.gz?download=1"
        fil="ECCOv4r1_analysis.tar.gz"
    elseif nam=="release3"
        url="https://zenodo.org/record/5786213/files/ECCOv4r3_analysis.tar.gz?download=1"
        fil="ECCOv4r3_analysis.tar.gz"
    elseif nam=="release4"
        url="https://zenodo.org/record/5786229/files/ECCOv4r4_analysis.tar.gz?download=1"
        fil="ECCOv4r4_analysis.tar.gz"
    elseif nam=="interp_coeffs"
        url="https://zenodo.org/record/5784905/files/interp_coeffs_halfdeg.jld2?download=1"
        fil="interp_coeffs_halfdeg.jld2"
    else
        println("unknown release name")
        url=missing
        fil=missing
    end
    if (!ismissing(url))&&(nam=="interp_coeffs")
        Downloads.download(url,joinpath(ECCOdiags_path,fil))
    elseif (!ismissing(url))&&(!isdir(joinpath(ECCOdiags_path,fil)[1:end-7]))
        println("downloading "*nam*" ... started")
        Downloads.download(url,joinpath(ECCOdiags_path,fil))
        tmp_path=open(joinpath(ECCOdiags_path,fil)) do io
            Tar.extract(CodecZlib.GzipDecompressorStream(io))
        end
        mv(joinpath(tmp_path,fil[1:end-7]),joinpath(ECCOdiags_path,fil[1:end-7]))
        rm(joinpath(ECCOdiags_path,fil))
        println("downloading "*nam*" ... done")
    end
end

end # module


text:
using JuLIP, JuLIP.Potentials
using JuLIPMaterials: Vec3
import JuLIPMaterials
const CLE = JuLIPMaterials.CLE

using JuLIPMaterials.CLE: grad
using JuLIPMaterials.Testing
using Einsum
using GaussQuadrature: legendre
using LinearAlgebra, ForwardDiff, Printf, LinearAlgebra

println("------------------------------------------------------------")
println(" Testing the 3D Anisotropic Green's Function Implementation")
println("------------------------------------------------------------")

print("check conversion between spherical and euclidean: ")
maxerr = 0.0
for n = 1:10
   x = randvec3(1.0, 1.0)
   global maxerr = max(maxerr, norm(x - CLE.euclidean(CLE.spherical(x)...)))
end
println(@test maxerr < 1e-14)
println("maxerr = $maxerr")

print("Test agreement between anisotopic and isotropic Gr fcn: ")
λ = 1.0 + rand()
μ = 1.0 + rand()
Ciso = CLE.isotropic_moduli(λ, μ)
G = CLE.GreenFunction(Ciso, Nquad = 4)
Giso = CLE.IsoGreenFcn3D(λ, μ)
maxerr = 0.0
maxerr_g = 0.0
for n = 1:10
   x = randvec3()
   global maxerr = max( maxerr, norm(G(x) - Giso(x), Inf) )
   global maxerr_g = max(maxerr_g, norm(grad(G, x) - grad(Giso, x), Inf) )
end
println("maxerr = $maxerr, maxerr_g = $maxerr_g")
println(@test maxerr < 1e-12)
println(@test maxerr_g < 1e-12)

Crand = randmoduli()

for (G, id, C) in [ (CLE.IsoGreenFcn3D(λ, μ), "IsoGreenFcn3D", Ciso),
            #  The next one SHOULD fail, so it is commented out
            # (CLE.GreenFunction(Crand, Nquad = 6), "GreenFunction(6)", Crand),
           (CLE.GreenFunction(Crand, Nquad = 30), "GreenFunction(30)", Crand) ]
   print("G = $id : test that ∇G is consistent with G: ")
   maxerr = 0.0
   for n = 1:10
      a, x = randvec3(), randvec3()
      u = x_ -> (a' * G(x_))[:]
      ∂u = x_ -> reshape(a' *  reshape(grad(G, x_), 3, 9), 3, 3)
      ∂uad = x_ -> ForwardDiff.jacobian(u, x_)
      maxerr = max( maxerr, norm(∂u(x) - ∂uad(x), Inf) )
   end
   println("maxerr = $maxerr")
   println(@test maxerr < 1e-12)
end

for (G, id, C) in [
            (CLE.IsoGreenFcn3D(λ, μ), "IsoGreenFcn3D", Ciso),
            (CLE.GreenFunction(Crand, Nquad = 30), "GreenFunction(30)", Crand) ]
   print("G = $id: test that G satisfies the PDE: ")
   maxerr = 0.0
   for n = 1:10
      a, x = randvec3(), randvec3()
      u = x_ -> G(x_) * a
      maxerr = max( norm(cleforce(Vec3(x), u, C), Inf), maxerr )
   end
   println("maxerr = $maxerr")
   println(@test maxerr < 1e-10)
end

for (G, id, C) in [(CLE.IsoGreenFcn3D(λ, μ), "IsoGreenFcn3D", Ciso),
                   (CLE.GreenFunction(Crand, Nquad = 30), "GreenFunction(30)", Crand) ]
   print("G = $id : test normalisation of G: ")
   err = 0.0

   # Test normal derivative integral over sphere via Gaussian quadrature
   # (Could use Lebedev, but no obvious Julia package)
   n = 30;
   c, w = legendre(n)
   I = zeros(3,3)
   DGnu = zeros(3,3)
   for ω in range(0.0, step=pi/n, length=2*n), i=1:n
      x = [sqrt(1-c[i]^2)*cos(ω),sqrt(1-c[i]^2)*sin(ω),c[i]]
      # TODO: @einsum seems to have changed behaviour or have a bug
      # @einsum DGnu[a,b]  = C[a,β,γ,δ] * CLE.grad(G,x)[b,γ,δ] * x[β]
      for a = 1:3, b = 1:3
         DGnu[a,b] = sum( C[a,β,γ,δ] * CLE.grad(G,x)[b,γ,δ] * x[β]
                          for β=1:3,γ=1:3,δ=1:3 )
      end
      I -= DGnu*w[i]
   end
   I = I*pi/n
   maxerr = norm(I  - [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0])
   println("maxerr = $maxerr")
   println(@test maxerr < 1e-12)
end

println("Convergence of G with random C (please test this visually!): ")
println(" nquad |    err    |   err_g")
println("-------|-----------|-----------")
C = Crand
xtest = [ randvec3() for n = 1:10 ]
G0 = CLE.GreenFunction(C, Nquad = 30)
for nquad in [2, 4, 6, 8, 10, 12, 14]
   global xtest, G0, G
   G = CLE.GreenFunction(C, Nquad = nquad)
   err = maximum( norm(G0(x) - G(x), Inf)   for x in xtest )
   err_g = maximum( norm(grad(G0, x) - grad(G, x), Inf)   for x in xtest )
   @printf("   %2d  | %.3e | %.3e  \n", nquad, err, err_g)
end


text:
# GPU memory management and pooling

using Printf
using Logging
using TimerOutputs

include("pool/utils.jl")
using .PoolUtils

## allocation statistics

mutable struct AllocStats
  # pool allocation requests
  pool_nalloc::Int
  pool_nfree::Int
  ## in bytes
  pool_alloc::Int

  # actual CUDA allocations
  actual_nalloc::Int
  actual_nfree::Int
  ## in bytes
  actual_alloc::Int
  actual_free::Int

  pool_time::Float64
  actual_time::Float64
end

const alloc_stats = AllocStats(0, 0, 0, 0, 0, 0, 0, 0, 0)

Base.copy(alloc_stats::AllocStats) =
  AllocStats((getfield(alloc_stats, field) for field in fieldnames(AllocStats))...)

AllocStats(b::AllocStats, a::AllocStats) =
  AllocStats(
    b.pool_nalloc - a.pool_nalloc,
    b.pool_nfree - a.pool_nfree,
    b.pool_alloc - a.pool_alloc,
    b.actual_nalloc - a.actual_nalloc,
    b.actual_nfree - a.actual_nfree,
    b.actual_alloc - a.actual_alloc,
    b.actual_free - a.actual_free,
    b.pool_time - a.pool_time,
    b.actual_time - a.actual_time)

## CUDA allocator

const alloc_to = TimerOutput()

"""
    alloc_timings()

Show the timings of the CUDA allocator. Assumes [`CUDA.enable_timings()`](@ref) has been
called.
"""
alloc_timings() = (show(alloc_to; allocations=false, sortby=:name); println())

const usage = PerDevice{Threads.Atomic{Int}}() do dev
  Threads.Atomic{Int}(0)
end
const usage_limit = PerDevice{Int}() do dev
  if haskey(ENV, "JULIA_CUDA_MEMORY_LIMIT")
    parse(Int, ENV["JULIA_CUDA_MEMORY_LIMIT"])
  elseif haskey(ENV, "CUARRAYS_MEMORY_LIMIT")
    Base.depwarn("The CUARRAYS_MEMORY_LIMIT environment flag is deprecated, please use JULIA_CUDA_MEMORY_LIMIT instead.", :__init_pool__)
    parse(Int, ENV["CUARRAYS_MEMORY_LIMIT"])
  else
    typemax(Int)
  end
end

@memoize function allocatable_memory(dev::CuDevice)
  # NOTE: this function queries available memory, which obviously changes after we allocate,
  #       so we memoize it to ensure only the first value is ever returned.
  device!(dev) do
    Base.min(available_memory(), usage_limit[dev])
  end
end

function reserved_memory(dev::CuDevice)
  # taken from TensorFlow's `MinSystemMemory`
  #
  # if the available memory is < 2GiB, we allocate 225MiB to system memory.
  # otherwise, depending on the capability version assign
  #  500MiB (for cuda_compute_capability <= 6.x) or
  # 1050MiB (for cuda_compute_capability <= 7.x) or
  # 1536MiB (for cuda_compute_capability >= 8.x)
  available = allocatable_memory(dev)
  if available <= 1<<31
    225 * 1024 * 1024
  else
    cap = capability(dev)
    if cap <= v"6"
      500 * 1024 * 1024
    elseif cap <= v"7"
      1050 * 1024 * 1024
    else
      1536 * 1024 * 1024
    end
  end
end

# This limit depends on the actually available memory, which might be an underestimation
# (e.g. when the GPU was in use initially). The hard limit does not rely on such heuristics.
@memoize function soft_limit(dev::CuDevice)
  available = allocatable_memory(dev)
  reserve = reserved_memory(dev)
  return available - reserve
end

function hard_limit(dev::CuDevice)
  # ignore the available memory heuristic, and even allow to eat in to the reserve
  usage_limit[dev]
end

function actual_alloc(dev::CuDevice, bytes::Integer, last_resort::Bool=false)
  buf = @device! dev begin
    # check the memory allocation limit
    if usage[dev][] + bytes > (last_resort ? hard_limit(dev) : soft_limit(dev))
      return nothing
    end

    # try the actual allocation
    try
      time = Base.@elapsed begin
        @timeit_debug alloc_to "alloc" begin
          buf = Mem.alloc(Mem.Device, bytes; async=true,
                          pool = (async_alloc[] ? pool() : nothing))
          # we only need a memory pool when we'll be using the async allocator.
          # this avoids a needless warning when running under cuda-memcheck,
          # which doesn't support the stream-ordered memory allocator.
        end
      end

      Threads.atomic_add!(usage[dev], bytes)
      alloc_stats.actual_time += time
      alloc_stats.actual_nalloc += 1
      alloc_stats.actual_alloc += bytes

      buf
    catch err
      (isa(err, CuError) && err.code == ERROR_OUT_OF_MEMORY) || rethrow()
      return nothing
    end
  end

  return Block(buf, bytes; state=AVAILABLE)
end

function actual_free(dev::CuDevice, block::Block)
  @assert iswhole(block) "Cannot free $block: block is not whole"
  @assert block.off == 0
  @assert block.state == AVAILABLE "Cannot free $block: block is not available"

  @device! dev begin
    # free the memory
    @timeit_debug alloc_to "free" begin
      time = Base.@elapsed begin
        Mem.free(block.buf; async=true)
      end
      block.state = INVALID

      Threads.atomic_sub!(usage[dev], sizeof(block.buf))
      alloc_stats.actual_time += time
      alloc_stats.actual_nfree += 1
      alloc_stats.actual_free += sizeof(block.buf)
    end
  end

  return
end

## memory pools

"""
    pool_timings()

Show the timings of the currently active memory pool. Assumes
[`CUDA.enable_timings()`](@ref) has been called.
"""
pool_timings() = (show(PoolUtils.to; allocations=false, sortby=:name); println())

# pool API:
# - init()
# - alloc(::CuDevice, sz)::Block
# - free(::CuDevice, ::Block)
# - reclaim(::CuDevice, nb::Int=typemax(Int))::Int
# - cached_memory()

module Pool
@enum MemoryPool None Simple Binned Split
end
const active_pool = Ref{Pool.MemoryPool}()
const async_alloc = Ref{Bool}()

macro pooled(ex)
    @assert Meta.isexpr(ex, :call)
    f, args... = ex.args
    quote
      if active_pool[] == Pool.None
        NoPool.$(f)($(map(esc, args)...))
      elseif active_pool[] == Pool.Simple
        SimplePool.$(f)($(map(esc, args)...))
      elseif active_pool[] == Pool.Binned
        BinnedPool.$(f)($(map(esc, args)...))
      elseif active_pool[] == Pool.Split
        SplitPool.$(f)($(map(esc, args)...))
      else
        error("unreachable")
      end
    end
end

include("pool/none.jl")
include("pool/simple.jl")
include("pool/binned.jl")
include("pool/split.jl")

## interface

export OutOfGPUMemoryError

const allocated_lock = NonReentrantLock()
const allocated = PerDevice{Dict{CuPtr,Tuple{Block,Int}}}() do dev
    Dict{CuPtr,Tuple{Block,Int}}()
end

const requested_lock = NonReentrantLock()
const requested = PerDevice{Dict{CuPtr{Nothing},Vector}}() do dev
  Dict{CuPtr{Nothing},Vector}()
end

"""
    OutOfGPUMemoryError()

An operation allocated too much GPU memory for either the system or the memory pool to
handle properly.
"""
struct OutOfGPUMemoryError <: Exception
  sz::Int

  OutOfGPUMemoryError(sz::Integer=0) = new(sz)
end

function Base.showerror(io::IO, err::OutOfGPUMemoryError)
    print(io, "Out of GPU memory")
    if err.sz > 0
      print(io, " trying to allocate $(Base.format_bytes(err.sz))")
    end
    println(io)
    memory_status(io)
end

"""
    alloc(sz)

Allocate a number of bytes `sz` from the memory pool. Returns a `CuPtr{Nothing}`; may throw
a [`OutOfGPUMemoryError`](@ref) if the allocation request cannot be satisfied.
"""
@inline function alloc(sz)
  # 0-byte allocations shouldn't hit the pool
  sz == 0 && return CU_NULL

  dev = device()

  time = Base.@elapsed begin
    @pool_timeit "pooled alloc" block = @pooled alloc(dev, sz)
  end
  block::Union{Nothing,Block}
  block === nothing && throw(OutOfGPUMemoryError(sz))

  # record the memory block
  ptr = pointer(block)
  @lock allocated_lock begin
      @assert !haskey(allocated[dev], ptr)
      allocated[dev][ptr] = block, 1
  end

  # record the allocation site
  if Base.JLOptions().debug_level >= 2
    bt = backtrace()
    @lock requested_lock begin
      @assert !haskey(requested[dev], ptr)
      requested[dev][ptr] = bt
    end
  end

  alloc_stats.pool_time += time
  alloc_stats.pool_nalloc += 1
  alloc_stats.pool_alloc += sz

  if MEMDEBUG && ptr == CuPtr{Cvoid}(0xbbbbbbbbbbbbbbbb)
    error("Allocated a scrubbed pointer")
  end

  return ptr
end

"""
    alias(ptr)

Increase the reference count of a buffer pointed to by `ptr`. As a result, it will require
multiple calls to `free` before this buffer is put back into the memory pool.
"""
@inline function alias(ptr::CuPtr{Nothing})
  # 0-byte allocations shouldn't hit the pool
  ptr == CU_NULL && return

  dev = device()

  # look up the memory block
  @spinlock allocated_lock begin
    block, refcount = allocated[dev][ptr]
    allocated[dev][ptr] = block, refcount+1
  end

  return
end

"""
    free(ptr)

Releases a buffer pointed to by `ptr` to the memory pool.
"""
@inline function free(ptr::CuPtr{Nothing})
  # 0-byte allocations shouldn't hit the pool
  ptr == CU_NULL && return

  dev = device()
  last_use[dev] = time()

  if MEMDEBUG && ptr == CuPtr{Cvoid}(0xbbbbbbbbbbbbbbbb)
    Core.println("Freeing a scrubbed pointer!")
  end

  # this function is typically called from a finalizer, where we can't switch tasks,
  # so perform our own error handling.
  try
    # look up the memory block, and bail out if its refcount isn't 1
    block = @spinlock allocated_lock begin
        block, refcount = allocated[dev][ptr]
        if refcount == 1
          delete!(allocated[dev], ptr)
          block
        else
          # we can't actually free this block yet, so decrease its refcount and return
          allocated[dev][ptr] = block, refcount-1
          return
        end
    end

    # look up the allocation site
    if Base.JLOptions().debug_level >= 2
      @lock requested_lock begin
        @assert haskey(requested[dev], ptr)
        delete!(requested[dev], ptr)
      end
    end

    time = Base.@elapsed begin
      @pool_timeit "pooled free" @pooled free(dev, block)
    end

    alloc_stats.pool_time += time
    alloc_stats.pool_nfree += 1
  catch ex
    Base.showerror_nostdio(ex, "WARNING: Error while freeing $ptr")
    Base.show_backtrace(Core.stdout, catch_backtrace())
    Core.println()
  end

  return
end

"""
    reclaim([sz=typemax(Int)])

Reclaims `sz` bytes of cached memory. Use this to free GPU memory before calling into
functionality that does not use the CUDA memory pool. Returns the number of bytes
actually reclaimed.
"""
function reclaim(sz::Int=typemax(Int))
  dev = device()
  @pooled reclaim(dev, sz)
end

"""
    @retry_reclaim isfailed(ret) ex

Run a block of code `ex` repeatedly until it successfully allocates the memory it needs.
Retries are only attempted when calling `isfailed` with the current return value is true.
At each try, more and more memory is freed from the CUDA memory pool. When that is not
possible anymore, the latest returned value will be returned.

This macro is intended for use with CUDA APIs, which sometimes allocate (outside of the
CUDA memory pool) and return a specific error code when failing to.
"""
macro retry_reclaim(isfailed, ex)
  quote
    ret = nothing
    for phase in 1:4
      ret = $(esc(ex))
      $(esc(isfailed))(ret) || break

      # incrementally more costly reclaim of cached memory
      if phase == 1
        reclaim()
      elseif phase == 2
        GC.gc(false)
        reclaim()
      elseif phase == 3
        GC.gc(true)
        reclaim()
      elseif phase == 4 && async_alloc[]
        # this phase is unique to retry_reclaim, as regular allocations come from the pool
        # so are assumed to never need to trim its contents.
        trim(pool())
      end
    end
    ret
  end
end

# XXX: function version for use in CUDAdrv where we haven't loaded pool.jl yet
function retry_reclaim(f, check)
  @retry_reclaim check f()
end

## management

const last_use = PerDevice{Union{Nothing,Float64}}() do dev
  nothing
end

# reclaim unused pool memory after a certain time
function pool_cleanup()
  while true
    t1 = time()
    @pool_timeit "cleanup" for dev in devices()
      t0 = last_use[dev]
      t0 === nothing && continue

      if t1-t0 > 300
        # the pool hasn't been used for a while, so reclaim unused buffers
        @pooled reclaim(dev)
      end
    end

    sleep(60)
  end
end

## utilities

used_memory(dev=device()) = @lock allocated_lock begin
    mapreduce(sizeof∘first, +, values(allocated[dev]); init=0)
end

"""
    @allocated

A macro to evaluate an expression, discarding the resulting value, instead returning the
total number of bytes allocated during evaluation of the expression.
"""
macro allocated(ex)
    quote
        let
            local f
            function f()
                b0 = alloc_stats.pool_alloc
                $(esc(ex))
                alloc_stats.pool_alloc - b0
            end
            f()
        end
    end
end

"""
    @time ex

Run expression `ex` and report on execution time and GPU/CPU memory behavior. The GPU is
synchronized right before and after executing `ex` to exclude any external effects.

"""
macro time(ex)
    quote
        local val, cpu_time,
        cpu_alloc_size, cpu_gc_time, cpu_mem_stats,
        gpu_alloc_size, gpu_gc_time, gpu_mem_stats = @timed $(esc(ex))

        local cpu_alloc_count = Base.gc_alloc_count(cpu_mem_stats)
        local gpu_alloc_count = gpu_mem_stats.pool_nalloc

        local gpu_lib_time = gpu_mem_stats.actual_time

        Printf.@printf("%10.6f seconds", cpu_time)
        for (typ, gctime, libtime, bytes, allocs) in
            (("CPU", cpu_gc_time, 0, cpu_alloc_size, cpu_alloc_count),
             ("GPU", gpu_gc_time, gpu_lib_time, gpu_alloc_size, gpu_alloc_count))
          if bytes != 0 || allocs != 0
              allocs, ma = Base.prettyprint_getunits(allocs, length(Base._cnt_units), Int64(1000))
              if ma == 1
                  Printf.@printf(" (%d%s %s allocation%s: ", allocs, Base._cnt_units[ma], typ, allocs==1 ? "" : "s")
              else
                  Printf.@printf(" (%.2f%s %s allocations: ", allocs, Base._cnt_units[ma], typ)
              end
              print(Base.format_bytes(bytes))
              if gctime > 0
                  Printf.@printf(", %.2f%% gc time", 100*gctime/cpu_time)
                if libtime > 0
                    Printf.@printf(" of which %.2f%% spent allocating", 100*libtime/gctime)
                end
              end
              print(")")
          elseif gctime > 0
              Printf.@printf(", %.2f%% %s gc time", 100*gctime/cpu_time, typ)
          end
        end
        println()

        val
    end
end

macro timed(ex)
    quote
        while false; end # compiler heuristic: compile this block (alter this if the heuristic changes)

        # @time(d) might surround an application, so be sure to initialize CUDA before that
        CUDA.initialize_cuda_context()

        # coarse synchronization to exclude effects from previously-executed code
        synchronize()

        local gpu_mem_stats0 = copy(alloc_stats)
        local cpu_mem_stats0 = Base.gc_num()
        local cpu_time0 = time_ns()

        # fine-grained synchronization of the code under analysis
        local val = @sync blocking=false $(esc(ex))

        local cpu_time1 = time_ns()
        local cpu_mem_stats1 = Base.gc_num()
        local gpu_mem_stats1 = copy(alloc_stats)

        local cpu_time = (cpu_time1 - cpu_time0) / 1e9

        local cpu_mem_stats = Base.GC_Diff(cpu_mem_stats1, cpu_mem_stats0)
        local gpu_mem_stats = AllocStats(gpu_mem_stats1, gpu_mem_stats0)

        (value=val, time=cpu_time,
         cpu_bytes=cpu_mem_stats.allocd, cpu_gctime=cpu_mem_stats.total_time / 1e9, cpu_gcstats=cpu_mem_stats,
         gpu_bytes=gpu_mem_stats.pool_alloc, gpu_gctime=gpu_mem_stats.pool_time, gpu_gcstate=gpu_mem_stats)
    end
end

cached_memory() = @pooled cached_memory()

"""
    memory_status([io=stdout])

Report to `io` on the memory status of the current GPU and the active memory pool.
"""
function memory_status(io::IO=stdout)
  dev = device()

  free_bytes, total_bytes = Mem.info()
  used_bytes = total_bytes - free_bytes
  used_ratio = used_bytes / total_bytes
  @printf(io, "Effective GPU memory usage: %.2f%% (%s/%s)\n",
              100*used_ratio, Base.format_bytes(used_bytes),
              Base.format_bytes(total_bytes))

  @printf(io, "CUDA allocator usage: %s", Base.format_bytes(usage[dev][]))
  if usage_limit[dev] !== typemax(Int)
    @printf(io, " (capped at %s)", Base.format_bytes(usage_limit[dev]))
  end
  println(io)

  alloc_used_bytes = used_memory()
  alloc_cached_bytes = cached_memory()
  alloc_total_bytes = alloc_used_bytes + alloc_cached_bytes
  @printf(io, "Memory pool '%s' usage: %s (%s allocated, %s cached)\n", string(active_pool[]),
              Base.format_bytes(alloc_total_bytes), Base.format_bytes(alloc_used_bytes),
              Base.format_bytes(alloc_cached_bytes))

  # check if the memory usage as counted by the CUDA allocator wrapper
  # matches what is reported by the pool implementation
  discrepancy = Base.abs(usage[dev][] - alloc_total_bytes)
  if discrepancy != 0
    println(io, "Discrepancy of $(Base.format_bytes(discrepancy)) between memory pool and allocator!")
  end

  if Base.JLOptions().debug_level >= 2
    requested′, allocated′ = @lock requested_lock begin
      copy(requested[dev]), copy(allocated[dev])
    end
    for (ptr, bt) in requested′
      block = allocated′[ptr]
      @printf(io, "\nOutstanding memory allocation of %s at %p",
              Base.format_bytes(sizeof(block)), Int(ptr))
      stack = stacktrace(bt, false)
      StackTraces.remove_frames!(stack, :alloc)
      Base.show_backtrace(io, stack)
      println(io)
    end
  end
end

## init

function __init_pool__()
  # usage
  initialize!(usage, ndevices())
  initialize!(last_use, ndevices())
  initialize!(usage_limit, ndevices())

  # allocation tracking
  initialize!(allocated, ndevices())
  initialize!(requested, ndevices())

  # memory pool configuration
  default_pool = version() >= v"11.2" ? "cuda" : "binned"
  pool_name = get(ENV, "JULIA_CUDA_MEMORY_POOL", default_pool)
  active_pool[], async_alloc[] = if pool_name == "none"
      Pool.None, false
  elseif pool_name == "simple"
      Pool.Simple, false
  elseif pool_name == "binned"
      Pool.Binned, false
  elseif pool_name == "split"
      Pool.Split, false
  elseif pool_name == "cuda"
      @assert version() >= v"11.2" "The CUDA memory pool is only supported on CUDA 11.2+"
      Pool.None, true
  else
      error("Invalid memory pool '$pool_name'")
  end
  @pooled init()

  TimerOutputs.reset_timer!(alloc_to)
  TimerOutputs.reset_timer!(PoolUtils.to)

  if isinteractive()
    @async @pooled pool_cleanup()
  end
end


text:
include("LS-method-7s-init.jl");

begin
    @load "LS_pureIOU_2000_saved.jld2" LS_pureIOU_2000_saved;
    r_list=LS_pureIOU_2000_saved[1]
    a_list=LS_pureIOU_2000_saved[2]
    Λ1=LS_pureIOU_2000_saved[3]
    tbl2=LS_pureIOU_2000_saved[4]
end;

## Two important checkpoints!
sum(Λ1) # should print 1
[sum(Λ1[i,:]) for i=1:nos]-BSP # should print a list of zeros

begin
    A=hcat(r_list,a_list)
    Aplus=A[A[:,2].>0,:]
    Aminus=A[A[:,2].<=0,:]
    iloc1=findmin(Aplus[:,2])[2]
    iloc2=findmax(Aminus[:,2])[2]
end;

#the nearest interest rates
Aplus[iloc1,1],Aminus[iloc2,1]

#the respective aggregate demand
Aplus[iloc1,2],Aminus[iloc2,2]

#the gap in the demand a(r)
Aplus[iloc1,2]-Aminus[iloc2,2]

#the gap in the interest rate is
Aplus[iloc1,1]-Aminus[iloc2,1]

begin
    scatter(a_list,r_list,label="",color=:black,border=true)
    plot!((0.02:0.0001:0.0305)*0,0.02:0.0001:0.0305,color=:black,label="")
    xlabel!("average assets")
    ylabel!("interest rate")
end

begin
    r=r_list[end]
    b=3.0;
    minkap=-b;
    if r<=0.0
        minkap = -b;
    else
        minkap=-min(b, wage*hours[1]/r);
    end;
    maxkap=16.0;
    Lg=2000; # The only place to set the grid size.
    grda=[minkap+i*(maxkap-minkap)/(Lg-1) for i=0:(Lg-1)];
    β=exp(-ρ);
end;

begin
    bond_price=ι/(1+r)
    cons=[(1+r)*grda[i]+wage*hours[k]-grda[Int64(tbl2[k,i])] for k=1:nos, i=1:Lg];
    aassets=[grda[Int64(tbl2[k,i])] for k=1:nos, i=1:Lg];
    shares=[grda[Int64(tbl2[k,i])]/bond_price for k=1:nos, i=1:Lg];
end;

begin
    SL=[sum(Λ1[state,:]) for state=1:nos];
    LF=[[Λ1[state,1]/SL[state]] for state=1:nos];
    for i=2:Lg
        for state=1:nos
            LF[state]=vcat(LF[state],LF[state][end]+Λ1[state,i]/SL[state])
        end
    end
end;

begin
    plot(aassets[1,2:920],LF[1][2:920],label="(LS)",color=:black)
    plot!(aassets[2,2:920],LF[2][2:920],label="",color=:black)
    plot!(aassets[3,2:920],LF[3][2:920],label="",color=:black)
    plot!(aassets[4,2:920],LF[4][2:920],label="",color=:black)
    plot!(aassets[5,2:920],LF[5][2:920],label="",color=:black)
    plot!(aassets[6,2:920],LF[6][2:920],label="",color=:black)
    plot!(aassets[7,2:920],LF[7][2:920],label="",color=:black)
end


text:
# Define Noise type which the solvers accepts
function default_seed()
    return rand(UInt)
end

function default_rng(seed)
    return MersenneTwister(seed)
end

"""
Arguments
T : Noise Type
m : Specifies the number of noise channels

Keyword arguments
seed: There is another more convenient constructor that given a seed instantiates a RNG using the deault_rng function.
"""
struct Noise
    m::Int
    rng::AbstractRNG
    Noise(m::Int, rng::AbstractRNG) = new(m, rng)
    Noise(m::Int; seed::Integer=default_seed()) = new(m, default_rng(seed))
end

W0(noise::Noise) = zeros(Float64, noise.m)

import Base:randn!

function randn!(noise::Noise, ΔW)
    randn!(noise.rng, ΔW)
end

max_stars_repo_name: onoderat/SDEs

max_issues_repo_name: onoderat/SDEs

max_forks_repo_name: onoderat/SDEs


text:
# Helper functions to read Harwell-Boeing and Rutherford-Boeing data.

function decode_int_fmt(fmt :: AbstractString)
  if fmt[1] == '('
    fmt = uppercase(fmt[2:end-1])
  end
  return map(s -> isempty(s) ? 1 : parse(Int, s), split(fmt, 'I'))
end

function decode_real_fmt(fmt :: AbstractString)
  fmt = join(split(fmt))  # Remove all white spaces.
  if fmt[1] == '('
    fmt = uppercase(fmt[2:end-1])
  end
  scale = "0"
  if (',' in fmt)  # Process scale factor, e.g., 1P,5D16.9
    scale, fmt = split(fmt, ',')
    scale, _ = split(scale, 'P')
  elseif ('P' in fmt)
    scale, fmt = split(fmt, 'P')
  end
  scale = parse(Int, scale)

  fmt1 = split(fmt, '.')[1]
  if occursin('E', fmt1)
    (npl, len) = map(s -> isempty(s) ? 1 : parse(Int, s), split(fmt1, 'E'))
  elseif occursin('D', fmt1)
    (npl, len) = map(s -> isempty(s) ? 1 : parse(Int, s), split(fmt1, 'D'))
  elseif occursin('F', fmt1)
    (npl, len) = map(s -> isempty(s) ? 1 : parse(Int, s), split(fmt1, 'F'))
  else
    error("Malformed real format")
  end

  return (npl, len, scale)
end

function standardize_real(number_as_str :: AbstractString)
  s = join(split(number_as_str))  # for numbers in the form "0.24555165E 00".
  # change ".16000000+006" to ".16000000e+006". The first char could be +/-.
  if !any(occursin.(['E', 'D', 'e', 'd'], s))
    if occursin('+', s[2:end])
      s = s[1:1] * join(split(s[2:end], '+'), "e+")
    elseif occursin('-', s[2:end])
      s = s[1:1] * join(split(s[2:end], '-'), "e-")
    end
  end
  return s
end

function read_array(io :: IO, n :: Int, fmt :: AbstractString; is_complex :: Bool=false)
  if 'I' in fmt
    scale = 0
    (npl, len) = decode_int_fmt(fmt)
    conv = s -> parse(Int, s)
    typ = Int
  else
    (npl, len, scale) = decode_real_fmt(fmt)
    conv = s -> parse(Float64, s)
    typ = Float64
    if is_complex
      n *= 2
    end
  end

  x = zeros(typ, n)
  for j = 1 : div(n, npl)
    if typ == Float64
      line = join(split(uppercase(readline(io)), 'D'), 'e')
    else
      line = readline(io)
    end
    chunk = [line[len*(i-1)+1:len*i] for i = 1 : npl]
    if typ == Float64
      chunk = map(standardize_real, chunk)
    end
    x[npl * (j-1) + 1 : npl * j] = map(conv, chunk)
  end
  rem = mod(n, npl)
  if rem > 0
    if typ == Float64
      line = join(split(uppercase(readline(io)), 'D'), 'e')
    else
      line = readline(io)
    end
    chunk = [line[len*(i-1)+1:len*i] for i = 1 : rem]
    if typ == Float64
      chunk = map(standardize_real, chunk)
    end
    x[end-rem+1 : end] = map(conv, chunk)
  end
  if scale != 0
    x /= 10.0^scale
  end
  return is_complex ? [ComplexF64(x[i], x[i+1]) for i = 1 : 2 : n-1] : x
end

function sortsparse!(colptr :: Vector{Ti}, rowind :: Vector{Ti}, values :: Vector{Tf}) where {Ti <: Integer, Tf <: Number}
  # ensure row indices are sorted in each column
  ncol = length(colptr) - 1
  for col = 1 : ncol
    colbeg = colptr[col]
    colend = colptr[col + 1] - 1
    rows = rowind[colbeg:colend]
    if !issorted(rows)
      p = sortperm(rows)
      rowind[colbeg:colend] = rows[p]
      values[colbeg:colend] = values[colbeg:colend][p]
    end
  end
end