text:
module Benchmarks

using BenchmarkTools
using Unitful
using StaticArrays

using ..OpticSim
using ..OpticSim: Sphere, Ray, replprint, trace, Cylinder, AcceleratedParametricSurface, newton, Infinity, RayOrigin, NullInterface, IntervalPoint, intervalintersection, LensTrace, FresnelInterface, wavelength, mᵢandmₜ, refractedray, direction, origin, pathlength, LensAssembly, NoPower, power, snell, fresnel, reflectedray, BoundingBox
using ..OpticSim.TestData
using ..OpticSim.Examples

rayz() = Ray([0.0, 0.0, 10.0], [0.0, 0.0, -1.0])
perturbrayz() = Ray([0.0, 0.0, 10.0], [0.001, 0.001, -1.0])
rayx() = Ray([10.0, 0.0, 0.0], [-1.0, 0.0, 0.0])
bezierray() = Ray([0.5, 0.5, 10.0], [0.0, 0.0, -1.0])
optray() = OpticalRay(origin(rayz()), direction(rayz()), 1.0, 0.55)
perturboptray() = OpticalRay(origin(perturbrayz()), direction(perturbrayz()), 1.0, 0.55)

#############################
### BENCHMARK DEFINITIONS ###
#############################

# benchmarks are listed in the form: benchmark() = function, (args...)
# e.g. surfaceintersection, (surface, ray)

double_convex() = trace, (TestData.doubleconvex(), optray())
double_concave() = trace, (TestData.doubleconcave(), optray())
zernike_lens() = trace, (TestData.zernikesystem(), perturboptray())
aspheric_lens() = trace, (TestData.asphericsystem(), perturboptray())
cooke_triplet() = trace, (TestData.cooketriplet(), optray())
chebyshev_lens() = trace, (TestData.chebyshevsystem(), perturboptray())
gridsag_lens() = trace, (TestData.gridsagsystem(), perturboptray())
multi_hoe() = trace, (TestData.multiHOE(), OpticalRay(SVector(0.0, 3.0, 3.0), SVector(0.0, -1.0, -1.0), 1.0, 0.55))

system_benchmarks() = double_concave, double_convex, aspheric_lens, zernike_lens, chebyshev_lens, gridsag_lens, cooke_triplet, multi_hoe

triangle() = surfaceintersection, (Triangle(SVector(-1.0, -1.0, 0.0), SVector(1.0, 0.0, 0.0), SVector(0.0, 1.0, 0.0)), rayz())
sphericalcap() = surfaceintersection, (SphericalCap(1.0, π / 3), rayz())
sphere_outside() = surfaceintersection, (Sphere(0.5), rayz())
sphere_inside() = surfaceintersection, (Sphere(0.5), Ray([0.0, 0.0, 0.0], [0.0, 0.0, -1.0]))
plane() = surfaceintersection, (Plane(0.0, 0.0, 1.0, 0.0, 0.0, 0.0), rayz())
rectangle() = surfaceintersection, (Rectangle(1.0, 1.0, SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, 0.0)), rayz())
ellipse() = surfaceintersection, (Ellipse(1.0, 1.0, SVector(0.0, 0.0, 1.0), SVector(0.0, 0.0, 0.0)), rayz())
cylinder_parallel() = surfaceintersection, (Cylinder(1.0), rayz())
cylinder_perp() = surfaceintersection, (Cylinder(1.0), rayx())
bbox_parallel() = surfaceintersection, (BoundingBox(-1.0, 1.0, -1.0, 1.0, -5.0, 5.0), rayz())
bbox_perp() = surfaceintersection, (BoundingBox(-1.0, 1.0, -1.0, 1.0, -5.0, 5.0), rayx())

simple_surface_benchmarks() = plane, rectangle, ellipse, triangle, sphere_outside, sphere_inside, sphericalcap, cylinder_parallel, cylinder_perp, bbox_parallel, bbox_perp

zernike_surface1() = surfaceintersection, (AcceleratedParametricSurface(TestData.zernikesurface1(), 30), perturbrayz())
zernike_surface2() = surfaceintersection, (AcceleratedParametricSurface(TestData.zernikesurface3(), 30), perturbrayz())
qtype_surface() = surfaceintersection, (AcceleratedParametricSurface(TestData.qtypesurface1(), 30), perturbrayz())
bezier_surface() = surfaceintersection, (AcceleratedParametricSurface(TestData.beziersurface(), 30), bezierray())
chebyshev_surface1() = surfaceintersection, (AcceleratedParametricSurface(TestData.chebyshevsurface1(), 30), perturbrayz())
chebyshev_surface2() = surfaceintersection, (AcceleratedParametricSurface(TestData.chebyshevsurface2(), 30), perturbrayz())
gridsag_surface() = surfaceintersection, (AcceleratedParametricSurface(TestData.gridsagsurfacebicubic(), 30), perturbrayz())

accel_surface_benchmarks() = zernike_surface1, zernike_surface2, qtype_surface, bezier_surface, chebyshev_surface1, chebyshev_surface2, gridsag_surface

all_benchmarks() = (simple_surface_benchmarks()..., accel_surface_benchmarks()..., system_benchmarks()...)

#############################

function runbenchmark(b; kwargs...)
    f, args = b()
    qkwargs = [:($(keys(kwargs)[i]) = $(values(kwargs)[i])) for i in 1:length(kwargs)]
    if f === trace
        return @eval (@benchmark($f(($args)..., test = true), $(qkwargs...)))
    else
        return @eval (@benchmark($f(($args)...), setup = (OpticSim.emptyintervalpool!()), $(qkwargs...)))
    end
end

timings(f::Function) = timings(f()...)
function timings(functions::Vararg{Function})
    for func in functions
        println("timing $func")
        replprint(runbenchmark(func))
        println("\n\n***************\n")
    end
end

function runforpipeline(ismaster::Bool)
    io = open(ismaster ? "timings_master.csv" : "timings_pr.csv", "w")
    write(io, "Function, Memory, Allocs, Min Time, Mean Time, Max Time\n")
    for f in all_benchmarks()
        b = runbenchmark(f)
        write(io, "$f, $(BenchmarkTools.prettymemory(b.memory)), $(b.allocs), $(BenchmarkTools.prettytime(minimum(b.times))), $(BenchmarkTools.prettytime(BenchmarkTools.mean(b.times))), $(BenchmarkTools.prettytime(maximum(b.times)))\n")
    end
    close(io)
end

########################################################

end #module

export Benchmarks


text:
## Generic and Fallbacks
"""
    discount(yc, to)
    discount(yc, from,to)

The discount factor for the yield curve `yc` for times `from` through `to`.
"""
discount(yc::T, from, to) where {T<:AbstractYieldCurve}= discount(yc, to) / discount(yc, from)

"""
    forward(yc, from, to, CompoundingFrequency=Periodic(1))

The forward `Rate` implied by the yield curve `yc` between times `from` and `to`.
"""
function forward(yc::T, from, to) where {T<:AbstractYieldCurve}
    return forward(yc, from, to, DEFAULT_COMPOUNDING)
end

function forward(yc::T, from, to, cf::CompoundingFrequency) where {T<:AbstractYieldCurve}
    r = Periodic((accumulation(yc, to) / accumulation(yc, from))^(1 / (to - from)) - 1, 1)
    return convert(cf, r)
end

function forward(yc::T, from) where {T<:AbstractYieldCurve}
    to = from + 1
    return forward(yc, from, to)
end

function CompoundingFrequency(curve::T) where {T<:AbstractYieldCurve}
    return DEFAULT_COMPOUNDING
end

"""
    par(curve,time;frequency=2)

Calculate the par yield for maturity `time` for the given `curve` and `frequency`. Returns a `Rate` object with periodicity corresponding to the `frequency`. The exception to this is if `time` is less than what the payments allowed by frequency (e.g. a time `0.5` but with frequency `1`) will effectively assume frequency equal to 1 over `time`.

# Examples

julia> c = Yields.Constant(0.04);

julia> Yields.par(c,4)
Yields.Rate{Float64, Yields.Periodic}(0.03960780543711406, Yields.Periodic(2))

julia> Yields.par(c,4;frequency=1)
Yields.Rate{Float64, Yields.Periodic}(0.040000000000000036, Yields.Periodic(1))

julia> Yields.par(c,0.6;frequency=4)
Yields.Rate{Float64, Yields.Periodic}(0.039413626195875295, Yields.Periodic(4))

julia> Yields.par(c,0.2;frequency=4)
Yields.Rate{Float64, Yields.Periodic}(0.039374942589460726, Yields.Periodic(5))

julia> Yields.par(c,2.5)
Yields.Rate{Float64, Yields.Periodic}(0.03960780543711406, Yields.Periodic(2))

"""
function par(curve, time; frequency=2)
    mat_disc = discount(curve, time)
    coup_times = coupon_times(time,frequency)
    coupon_pv = sum(discount(curve,t) for t in coup_times)
    Δt = step(coup_times)
    r = (1-mat_disc) / coupon_pv
    cfs = [t == last(coup_times) ? 1+r : r for t in coup_times]
    # `sign(r)`` is used instead of `1` because there are times when the coupons are negative so we want to flip the sign
    cfs = [-1;cfs]
    r = internal_rate_of_return(cfs,[0;coup_times])
    frequency_inner = min(1/Δt,max(1 / Δt, frequency))
    r = convert(Periodic(frequency_inner),r)
    return r
end

"""
    zero(curve,time)
    zero(curve,time,CompoundingFrequency)

Return the zero rate for the curve at the given time.
"""
function Base.zero(c::YC, time) where {YC<:AbstractYieldCurve}
     zero(c, time, CompoundingFrequency(c))
end

function Base.zero(c::YC, time, cf::C) where {YC<:AbstractYieldCurve,C<:CompoundingFrequency}
    df = discount(c, time)
    r = -log(df)/time
    return convert(cf, Continuous(r)) # c.zero is a curve of continuous rates represented as floats. explicitly wrap in continuous before converting
end

"""
    accumulation(yc, from, to)

The accumulation factor for the yield curve `yc` for times `from` through `to`.
"""
function accumulation(yc::AbstractYieldCurve, time)
    return 1 ./ discount(yc, time)
end

function accumulation(yc::AbstractYieldCurve, from, to)
    return 1 ./ discount(yc, from, to)
end

"""
    Par(rates, maturities=eachindex(rates)
    Par(p::YieldCurveFitParameters, rates, maturities=eachindex(rates)

Construct a curve given a set of bond equivalent yields and the corresponding maturities. Assumes that maturities <= 1 year do not pay coupons and that after one year, pays coupons with frequency equal to the CompoundingFrequency of the corresponding rate (normally the default for a `Rate` is `1`, but when constructed via `Par` the default compounding Frequency is `2`).

See [`bootstrap`](@ref) for more on the `interpolation` parameter, which is set to `QuadraticSpline()` by default.

# Examples

```julia-repl

julia> par = [6.,8.,9.5,10.5,11.0,11.25,11.38,11.44,11.48,11.5] ./ 100
julia> maturities = [t for t in 1:10]
julia> curve = Par(par,maturities);
julia> zero(curve,1)
Rate(0.06000000000000005, Periodic(1))

```
"""
function Yields.Par(rates,maturities=eachindex(rates))
    # bump to a constant yield if only given one rate
    length(rates) == 1 && return Constant(first(rates))
    return Yields.Par(Bootstrap(),rates,maturities)
end

"""
    Forward(rates,maturities=eachindex(rates))
    Forward(p::YieldCurveFitParameters, rates,maturities=eachindex(rates))

Takes a vector of 1-period forward rates and constructs a discount curve. The method of fitting the curve to the data is determined by the [`YieldCurveFitParameters`](@ref) object `p`, which is a `Boostrap(QuadraticSpline())` by default.

If `rates` is a vector of floating point number instead of a vector `Rate`s, see the [`YieldCurveFitParameters`](@ref) for how the rate will be interpreted.

# Examples

```julia-repl
julia> Yields.Forward( [0.01,0.02,0.03] );

julia> Yields.Forward( Yields.Continuous.([0.01,0.02,0.03]) );

```
"""
function Yields.Forward(rates,maturities=eachindex(rates))
    # bump to a constant yield if only given one rate
    length(rates) == 1 && return Constant(first(rates))
    return Yields.Forward(Bootstrap(),rates,maturities)
end

"""
    Zero(rates, maturities=eachindex(rates))
    Zero(p::YieldCurveFitParameters,rates, maturities=eachindex(rates))

Construct a yield curve with given zero-coupon spot `rates` at the given `maturities`. The method of fitting the curve to the data is determined by the [`YieldCurveFitParameters`](@ref) object `p`, which is a `Boostrap(QuadraticSpline())` by default.

If `rates` is a vector of floating point number instead of a vector `Rate`s, see the [`YieldCurveFitParameters`](@ref) for how the rate will be interpreted.

# Examples

```julia-repl
julia> Yields.Zero([0.01,0.02,0.04,0.05],[1,2,5,10])

              ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀Yield Curve (Yields.BootstrapCurve)
              ┌────────────────────────────────────────────────────────────┐
         0.05 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠤⠒⠒⠒⠒⠒⠤⠤⠤⢄⣀⣀⣀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ Zero rates
              │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠖⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠉⠉⠉⠒⠒⠒⠢⠤⠤⠤⣄⣀⣀⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
              │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠖⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠉⠉⠑⠒⠒⠒⠦⠤⠤⠤⣀⣀⣀⡀⠀│
              │⠀⠀⠀⠀⠀⠀⠀⠀⠀⡔⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠉│
              │⠀⠀⠀⠀⠀⠀⠀⢠⠎⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
              │⠀⠀⠀⠀⠀⠀⢠⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
              │⠀⠀⠀⠀⠀⢠⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
   Continuous │⠀⠀⠀⠀⢀⠇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
              │⠀⠀⠀⠀⡜⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
              │⠀⠀⠀⡸⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
              │⠀⠀⢠⠇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
              │⠀⠀⡜⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
              │⠒⠒⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
              │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
            0 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
              └────────────────────────────────────────────────────────────┘
              ⠀0⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀time⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀30

"""
function Yields.Zero(rates,maturities=eachindex(rates))
    # bump to a constant yield if only given one rate
    length(rates) == 1 && return Constant(first(rates))
    return Yields.Zero(Bootstrap(),rates,maturities)
end

"""
    Yields.CMT(rates, maturities; interpolation=QuadraticSpline())
    Yields.CMT(p::YieldCurveFitParameters,rates, maturities)

Takes constant maturity (treasury) yields (bond equivalent), and assumes that instruments <= one year maturity pay no coupons and that the rest pay semi-annual.

The method of fitting the curve to the data is determined by the [`YieldCurveFitParameters`](@ref) object `p`, which is a `Boostrap(QuadraticSpline())` by default.

# Examples

```
# 2021-03-31 rates from Treasury.gov
rates =[0.01, 0.01, 0.03, 0.05, 0.07, 0.16, 0.35, 0.92, 1.40, 1.74, 2.31, 2.41] ./ 100
mats = [1/12, 2/12, 3/12, 6/12, 1, 2, 3, 5, 7, 10, 20, 30]

Yields.CMT(rates,mats)
```
"""
function Yields.CMT(rates,maturities=eachindex(rates))
    # bump to a constant yield if only given one rate
    length(rates) == 1 && return Constant(first(rates))
    return Yields.CMT(Bootstrap(),rates,maturities)
end

"""
    OIS(rates, maturities)
    OIS(p::YieldCurveFitParameters, rates, maturities)

Takes Overnight Index Swap rates, and assumes that instruments <= one year maturity are settled once and other agreements are settled quarterly with a corresponding CompoundingFrequency.

The method of fitting the curve to the data is determined by the [`YieldCurveFitParameters`](@ref) object `p`, which is a `Boostrap(QuadraticSpline())` by default.

# Examples
```
julia> ois = [1.8, 2.0, 2.2, 2.5, 3.0, 4.0] ./ 100;
julia> mats = [1 / 12, 1 / 4, 1 / 2, 1, 2, 5];
julia> curve = Yields.OIS(ois, mats)

              ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀Yield Curve (Yields.BootstrapCurve)
              ┌────────────────────────────────────────────────────────────┐
          0.1 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ Zero rates
              │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⣀⡤│
              │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⣀⡤⠔⠒⠋⠁⠀⠀│
              │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⣀⡤⠔⠒⠋⠁⠀⠀⠀⠀⠀⠀⠀⠀│
              │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⣀⡤⠔⠒⠋⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
              │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⣀⡤⠔⠒⠋⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
              │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⣀⡤⠔⠒⠋⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
   Continuous │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⣀⡤⠔⠒⠋⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
              │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⣀⡤⠔⠒⠋⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
              │⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⣀⡤⠔⠒⠋⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
              │⠀⠀⠀⠀⠀⠀⡠⠔⠊⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
              │⠀⠀⠀⣀⠔⠋⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
              │⠀⡠⠎⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
              │⠜⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
         0.01 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
              └────────────────────────────────────────────────────────────┘
              ⠀0⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀time⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀30
```
"""
function Yields.OIS(rates,maturities=eachindex(rates))
    # bump to a constant yield if only given one rate
    length(rates) == 1 && return Constant(first(rates))
    return Yields.OIS(Bootstrap(),rates,maturities)
end

ext: jl


text:
flatten = A -> collect(Iterators.flatten(A));

function volumesphere(n)
    V = 1;
    if n==1
        V = 2;
    elseif n>1
        V = (2*π/n)*volumesphere(n-2);
    elseif n < 0
        throw(DomainError(n,"dimension argument must be non-negative"));
    end
    return V
end

function surfaceareasphere(n)
    return (n+1)*volumesphere(n+1)
end

function sphericalquadrature(n,deg)
    x,wx = Float64[1,-1],Float64[1,1];
    if n > 1
        z,wz = gaussjacobi(deg,(n-3)/2,(n-3)/2);
        y,wy = sphericalquadrature(n-1,deg);
        xA = map(x->flatten([x[1],sqrt(1-x[1]^2).*x[2]]),Base.product(z,y))
        x = xA[:];
        wxA = map(prod,Base.product(wz,wy))
        wx = wxA[:];
    elseif n < 1
        throw(DomainError(n,"dimension argument must be positive"));
    end
    return x,wx
end


text:
struct TGV{T <: Real, Int <: Integer} <: FluidFlowProblem
    q::Quadrature
    ρ_0::T
    u_max::T
    u_0::T
    τ::T
    ν::T
    NX::Int
    NY::Int
    static::Bool
    domain_size::Tuple{T, T}
end
function TGV(
    q::Quadrature,
    τ, # Effective relaxation time
    scale = 2,
    NX = 16 * scale,
    NY = NX,
    # u_max = sqrt(0.001) / scale
    u_max = 0.02 / scale,
)
    ν = (τ - 0.5) / (q.speed_of_sound_squared)
    Re = NX * u_max / ν

    TGV(q, 1.0, 1.0, u_max, τ, ν, NX, NY, false, (1.0, 1.0))
end

function density(q::Quadrature, problem::TGV, x::T, y::T, t::Real = 0.0) where {T <: Real}
    k_x = 2pi / problem.NX
    k_y = 2pi / problem.NY

    x *= problem.NX
    y *= problem.NY

    ν = problem.ν
    t_d = 1 / (ν * (k_x^2 + k_y^2))

    return problem.ρ_0 * (
        1.0 -
        q.speed_of_sound_squared *
        (problem.u_0^2 / 4) *
        ((k_y / k_x) * cos(2 * k_x * x) + (k_x / k_y) * cos(2 * k_y * y)) *
        exp(-2 * t / t_d)
    )
end

function pressure(q::Quadrature, problem::TGV, x::T, y::T, t::Real = 0.0) where {T <: Real}
    k_x = 2pi / problem.NX
    k_y = 2pi / problem.NY

    x *= problem.NX
    y *= problem.NY

    ν = problem.ν
    t_d = 1 / (ν * (k_x^2 + k_y^2))

    p_0 = problem.ρ_0
    ρ = 1.0

    return p_0 -
           q.speed_of_sound_squared *
    ρ *
    (problem.u_0^2 / 4) *
    ((k_y / k_x) * cos(2 * k_x * x) + (k_x / k_y) * cos(2 * k_y * y)) *
    exp(-2 * t / t_d)
end

function velocity(problem::TGV, x::T, y::T, t::Real = 0.0) where {T <: Real}
    k_x = 2pi / problem.NX
    k_y = 2pi / problem.NY

    x *= problem.NX
    y *= problem.NY

    ν = problem.ν
    t_d = 1 / (ν * (k_x^2 + k_y^2))

    return problem.u_0 *
    exp(-t / t_d) *
    [
        -sqrt(k_y / k_x) * cos(k_x * x) * sin(k_y * y)
        sqrt(k_x / k_y) * sin(k_x * x) * cos(k_y * y)
    ]
end

function velocity_gradient(problem::TGV, x::T, y::T, t::Real = 0.0) where {T <: Real}
    k_x = 2pi / problem.NX
    k_y = 2pi / problem.NY

    x *= problem.NX
    y *= problem.NY

    ν = problem.ν
    t_d = 1 / (ν * (k_x^2 + k_y^2))

    u_x = sqrt(k_y * k_x) * sin(k_x * x) * sin(k_y * y)
    v_y = -sqrt(k_y * k_x) * sin(k_x * x) * sin(k_y * y)

    u_y = -sqrt(k_y^3 / k_x) * cos(k_x * x) * cos(k_y * y)
    v_x = sqrt(k_x^3 / k_y) * cos(k_x * x) * cos(k_y * y)

    return exp(-t / t_d) * problem.u_0 * [u_x v_x; u_y v_y]
end

function decay(problem::TGV, x::T, y::T, t::Real) where {T <: Real}
    ν = problem.ν
    k_x = 2pi / problem.NX
    k_y = 2pi / problem.NY
    t_d = 1 / (ν * (k_x^2 + k_y^2))

    return problem.static ? 1.0 : exp(-t / t_d)
end

function force(problem::TGV, x::T, y::T, t::Real = 0.0) where {T <: Real}
    ν = problem.ν

    return problem.static ? 2 * ν * velocity(problem, x, y, 0.0) : [0.0 0.0]
end

has_external_force(problem::TGV) = problem.static

viscosity(problem::TGV) = problem.ν # (problem.τ - 0.5) / (problem.q.speed_of_sound_squared)
delta_x(problem::TGV) = 1.0
delta_t(problem::TGV) = 1.0

function decay_time(problem)
    ν = viscosity(problem)
    k_x = 2pi / problem.NX
    k_y = 2pi / problem.NY
    t_d = 1 / (ν * (k_x^2 + k_y^2))
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export dal_dashaladinlite

"""
    dal_dashaladinlite(;kwargs...)

A DashAladinLite component.
DashAladinLite
Keyword arguments:
- `id` (String; optional): The ID used to identify this component in Dash callbacks.
- `autoFov` (Bool; optional): Adjust the fov automatically for target if set to True.
- `className` (String; optional): The class of the container.
- `custom_script_calls` (Dict; optional): Custom script calls to make. The keys should match those in
``custom_script``, and the vialues are passed to
``custom_scripts``.
- `custom_scripts` (Dict; optional): Custom scripts to be used in ``custom_script_calls``.
The values are inline functions that have signature like
``function(aladin, data, props)``, where data are
passed in ``custom_script_calls``.
- `footprintClicked` (Dict; optional): Clicked footprint.
- `footprintHovered` (Dict; optional): Hovered footprint.
- `fov` (Real; optional): The FOV of the display.
- `layers` (optional): Layers to show.
Each layer can be one of types ['catalog', 'overlay', 'moc'].. layers has the following type: Array of lists containing elements 'type', 'data', 'options'.
Those elements have the following types:
  - `type` (a value equal to: "catalog"; optional)
  - `data` (optional): . data has the following type: String | Array of lists containing elements 'ra', 'dec'.
Those elements have the following types:
  - `ra` (Real; optional)
  - `dec` (Real; optional)s
  - `options` (optional): . options has the following type: lists containing elements 'show', 'name', 'color', 'sourceSize', 'shape', 'maxNbsources', 'raField', 'decField', 'displayLabel', 'labelColor', 'labelFont', 'labelColumn', 'onClick'.
Those elements have the following types:
  - `show` (Bool; optional)
  - `name` (String; optional)
  - `color` (String; optional)
  - `sourceSize` (Real; optional)
  - `shape` (a value equal to: 'circle', "plus", "rhomb", 'cross', 'triangle', 'square'; optional)
  - `maxNbsources` (Real; optional)
  - `raField` (String; optional)
  - `decField` (String; optional)
  - `displayLabel` (Bool; optional)
  - `labelColor` (String; optional)
  - `labelFont` (String; optional)
  - `labelColumn` (String; optional)
  - `onClick` (a value equal to: 'showTable', 'showPopup'; optional) | lists containing elements 'type', 'data', 'options'.
Those elements have the following types:
  - `type` (a value equal to: "overlay"; optional)
  - `data` (optional): . data has the following type: Array of lists containing elements 'type', 'ra', 'dec', 'radius', 'color'.
Those elements have the following types:
  - `type` (a value equal to: 'circle'; optional)
  - `ra` (Real; optional)
  - `dec` (Real; optional)
  - `radius` (Real; optional)
  - `color` (String; optional) | lists containing elements 'type', 'data', 'color', 'lineWidth'.
Those elements have the following types:
  - `type` (a value equal to: 'polyline'; optional)
  - `data` (Array of Arrays; optional)
  - `color` (String; optional)
  - `lineWidth` (Real; optional) | lists containing elements 'type', 'data'.
Those elements have the following types:
  - `type` (a value equal to: 'polygon'; optional)
  - `data` (Array of Arrays; optional)s
  - `options` (optional): . options has the following type: lists containing elements 'show', 'name', 'color', 'lineWidth'.
Those elements have the following types:
  - `show` (Bool; optional)
  - `name` (String; optional)
  - `color` (String; optional)
  - `lineWidth` (Real; optional) | lists containing elements 'type', 'data', 'options'.
Those elements have the following types:
  - `type` (a value equal to: 'moc'; optional)
  - `data` (String | Dict; optional)
  - `options` (optional): . options has the following type: lists containing elements 'show', 'name', 'color', 'opacity', 'adaptativeDisplay'.
Those elements have the following types:
  - `show` (Bool; optional)
  - `name` (String; optional)
  - `color` (String; optional)
  - `opacity` (Real; optional)
  - `adaptativeDisplay` (Bool; optional)s
- `objectClicked` (Dict; optional): Clicked object.
- `objectHovered` (Dict; optional): Hovered object.
- `options` (Dict; optional): Additional options passed to the aladinlite.js
- `position` (Dict; optional): Position.
- `style` (Dict; optional): The style of the container.
- `survey` (optional): The survey to display as base image layer.. survey has the following type: String | lists containing elements 'id', 'name', 'rootUrl', 'cooFrame', 'maxOrder', 'options'.
Those elements have the following types:
  - `id` (String; optional)
  - `name` (String; optional)
  - `rootUrl` (String; optional)
  - `cooFrame` (a value equal to: 'equatorial', 'galactic'; optional)
  - `maxOrder` (Real; optional)
  - `options` (optional): . options has the following type: lists containing elements 'imgFormat'.
Those elements have the following types:
  - `imgFormat` (a value equal to: 'jpeg', 'fits', 'png'; optional)
- `target` (String; required): The target to display.
"""
function dal_dashaladinlite(; kwargs...)
        available_props = Symbol[:id, :autoFov, :className, :custom_script_calls, :custom_scripts, :footprintClicked, :footprintHovered, :fov, :layers, :objectClicked, :objectHovered, :options, :position, :style, :survey, :target]
        wild_props = Symbol[]
        return Component("dal_dashaladinlite", "DashAladinLite", "dash_aladin_lite", available_props, wild_props; kwargs...)
end


text:
## sugar.jl - (hopefully) temporary copy of needed parts of Sugar.jl package

# jlhome() = ccall(:jl_get_julia_home, Any, ())
jlhome() = Sys.BINDIR

function juliabasepath(file)
    srcdir = joinpath(jlhome(),"..","..","base")
    releasedir = joinpath(jlhome(),"..","share","julia","base")
    normpath(joinpath(isdir(srcdir) ? srcdir : releasedir, file))
end

function get_source_file(path::AbstractString, ln)
    isfile(path) && return path
    # if not a file, it might be in julia base
    file = juliabasepath(path)
    if !isfile(file)
        throw(LoadError(path, Int(ln), ErrorException("file $path not found")))
    end
    file
end

function get_method(f, types::Type)
    get_method(f, (types.parameters...,))
end
function get_method(ftype::DataType, types::Tuple)
    world = typemax(UInt)
    if !isclosure(ftype)
        ftype = Type{ftype}
    end
    tt = Tuple{ftype, to_tuple(types)...}
    (ti, env, meth) = Base._methods_by_ftype(tt, 1, world)[1]
    Base.func_for_method_checked(meth, tt)
end
function get_method(f, types::Tuple)
    if !all(isconcretetype, types)
        error("Not all types are concrete: $types")
    end
    # make sure there is a specialization with precompile
    # TODO, figure out a better way, since this method is not very reliable.
    # (I think, e.g. anonymous functions don't work)
    precompile(f, (types...,))
    x = methods(f, types)
    if isempty(x)
        throw(NoMethodError(f, types))
    elseif length(x) != 1
        error("
            More than one method found for signature $f $types.
            Please use more specific types!
        ")
    end
    first(x)
end

"""
Looks up the source of `method` in the file path found in `method`.
Returns the AST and source string, might throw an LoadError if file not found.
"""
function get_source_at(file, linestart)
    file = get_source_file(file, linestart)
    code, str = open(file) do io
        line = ""
        for i=1:linestart-1
            line = readline(io)
        end
        try # lines can be one off, which will result in a parse error
            Meta.parse(line)
        catch e
            line = readline(io)
        end
        while !eof(io)
            line = line * "\n" * readline(io)
            e = Base.parse_input_line(line; filename=file)
            if !(isa(e,Expr) && e.head === :incomplete)
                return e, line
            end
        end
    end
    code, str
end


text:
"""
    An external product using Einstein summation.
"""
extern_prod(x, y) = @tullio threads = false z[i, j] := x[i] * y[j]

"""
    Pairwise multiplication when `x` and `a` both have two dimensions.
"""
pairwise_mul!(x, a) = @tullio threads = false x[i, j] *= a[i, j]

"""
    Pairwise difference when `x` and `y` are vectors.
"""
pairwise_diff(x, y) = @tullio threads = false z[i] := x[i] - y[i]

"""
    Computes a pairwise multiplication between `y`, `z` and the second dimension from
`x`. Then, it uses a sum reduction of the first dimension from `x`.

It is equivalent to

```julia
sum(@. x * y * z; dims=1)
```
"""
prod_reduction(x, y, z) = @tullio threads = false w[1, j] := x[i, j] * y[i] * z[i]

prod_reduction(x, y) = @tullio threads = false w[1, j] := x[i, j] * y[i]

max_stars_repo_name: edwinb-ai/LeastSquaresSVM

max_issues_repo_name: edwinb-ai/Elysivm

max_forks_repo_name: edwinb-ai/LeastSquaresSVM


text:
# MobilityTheories.jl

"""
    polaronmobility(Trange,
                    ε_Inf, ε_S, freq, effectivemass;
                    verbose::Bool=false)

    Solves the Feynman polaron problem variationally with finite temperature
    Osaka energies.  From the resulting v, and w parameters, calculates polaron
    structure (wave function size, etc.).  Uses FHIP, Kadanoff (Boltzmann
    relaxation time) and Hellwarth direct contour integration to predict
    a temperature-dependent mobility for the material system.
    Input is a temperature range (e.g. 10:50:1000),
    reduced dielectric constants (e.g. 5, 20),
    characteristic dielectric phonon frequency (e.g. 2.25E12) - units Hertz
    bare-band effective-mass (e.g. 012) - units electron mass.

    Returns a structure of type Polaron, containing arrays of useful
    information.  Also prints a lot of information to the standard out - which
    may be more useful if you're just inquiring as to a particular data point,
    rather than plotting a temperature-dependent parameter.

    As an example, to calculate the electron polaron in MAPI at 300 K:
    polaronmobility(300, 4.5, 24.1, 2.25E12, 0.12)
"""
function polaronmobility(Trange, ε_Inf, ε_S, freq, effectivemass; verbose::Bool=false)
    println("\n\nPolaron mobility for system ε_Inf=$ε_Inf, ε_S=$ε_S, freq=$freq,
                 effectivemass=$effectivemass; with Trange $Trange ...")

    # Internally we we use 'mb' for the 'band mass' in SI units, of the effecitve-mass of the electron
    mb=effectivemass*MassElectron
    ω = (2*pi)*freq # angular-frequency, of the phonon mode

    α=frohlichalpha(ε_Inf, ε_S,  freq,    effectivemass)
    #α=2.395939683378253 # Hard coded; from MAPI params, 4.5, 24.1, 2.25THz, 0.12me

    v=7.1 # starting guess for v,w variational parameters
    w=6.5

    @printf("Polaron mobility input parameters: ε_Inf=%f ε_S=%f freq=%g α=%f \n",ε_Inf, ε_S, freq, α )
    @printf("Derived params in SI: ω =%g mb=%g \n",ω ,mb)

    # Empty struct for storing data
    # A slightly better way of doing this ф_ф ...
    p=Polaron()
    # populate data structure with (athermal) parameters supplied...
    append!(p.α,α) # appending so as not to mess with type immutability
    append!(p.mb,mb)
    append!(p.ω,ω)

    # We define βred as the subsuming the energy of the phonon; i.e. kbT c.f. ħω
    for T in Trange
        β=1/(kB*T)
        βred=ħ*ω*β
        append!(p.βred,βred)
        @printf("T: %f β: %.3g βred: %.3g ħω  = %.3g meV\t",T,β,βred, 1E3* ħ*ω  / q)

        if T==0
            v,w=feynmanvw(α, v=v,w=w)
        else
            v,w=feynmanvw(α,βred, v=v,w=w)
        end

        @printf("\n Polaron Parameters:  v= %.4f  w= %.4f",v,w)

        # From 1962 Feynman, definition of v and w in terms of the coupled Mass and spring-constant
        # See Page 1007, just after equation (18)
        # Units of M appear to be 'electron masses'
        # Unsure of units for k, spring coupling constant
        k=(v^2-w^2)
        M=(v^2-w^2)/w^2
        append!(p.k,k)
        append!(p.M,M)

        @printf("  ||   M=%f  k=%f\t",M,k)
        @printf("\n Bare-band effective mass: %f Polaron effective mass: %f Polaron mass enhancement: %f%%",effectivemass,effectivemass*(1+M),M*100)

        @printf("\n Polaron frequency (SI) v= %.2g Hz  w= %.2g Hz",
                v*ω /(2*pi),w*ω /(2*pi))

		# (46) in Feynman1955
		meSmallAlpha(α )=α /6 + 0.025*α ^2
		# (47) In Feynman1955
		meLargeAlpha(α )=16*α ^4 / (81*π ^4)
		#meLargeAlpha(α )=202*(α /10)^4
        if (verbose) # asymptotic solutions - not that interesting when you have the actual ones!
		    @printf("\n Feynman1955(46,47): meSmallAlpha(α)= %.3f meLargeAlpha(α)= %.3f",
                    meSmallAlpha(α),meLargeAlpha(α))
            @printf("\n Feynman1962: Approximate ~ Large alpha limit, v/w = %.2f  =~approx~= alpha^2 = %.2f ",
                    v/w,α ^2)
        end

        # POLARON SIZE
        @printf("\n Polaron size (rf), following Schultz1959. (s.d. of Gaussian polaron ψ )")
        # Schultz1959 - rather nicely he actually specifies everything down into units!
        # just before (2.4) in Schultz1959
        mu=((v^2-w^2)/v^2)
        # (2.4)
        rf=sqrt(3/(2*mu*v))
        # (2.4) SI scaling inferred from units in (2.5a) and Table II
        rfsi=rf*sqrt(2*me*ω )
        @printf("\n\t Schultz1959(2.4): rf= %g (int units) = %g m [SI]",rf,rfsi )
        append!(p.rfsi,rfsi)

        if (verbose)
        scale=sqrt(2*mb*ω) # Note we're using mb;
        #band effective-mass in SI units (i.e. meff*melectron)

        rfa=(3/(0.44*α ))^0.5 # As given in Schultz1959(2.5a), but that 0.44 is actually 4/9
        @printf("\n\t Schultz1959(2.5a) with 0.44: Feynman α→0 expansion: rfa= %g (int units) = %g m [SI]",rfa,scale*rfa )
        rfa=(3/((4/9)*α ))^0.5 # Rederived from Feynman1955, 8-8-2017; Yellow 2017.B Notebook pp.33-34
        @printf("\n\t Schultz1959(2.5a) with 4/9 re-derivation: Feynman α→0 expansion: rfa= %g (int units) = %g m [SI]",rfa,scale*rfa )
        append!(p.rfsmallalpha,scale*rfa)

        rfb=3*(pi/2)^0.5 * α
        @printf("\n\t Schultz1959(2.5b): Feynman α→∞ expansion: rf= %g (int units) = %g m [SI]",rfb,scale*rfb )

        # Schultz1959 - Between (5.7) and (5.8) - resonance of Feynman SHM system
        phononfreq=sqrt(k/M)
        @printf("\n\t Schultz1959(5.7-5.8): fixed-e: phononfreq= %g (int units) = %g [SI, Hz] = %g [meV]",
            phononfreq,phononfreq*ω /(2*pi), phononfreq*hbar*ω *1000/q)

        phononfreq=sqrt(k/mu) # reduced mass
        @printf("\n\t Schultz1959(5.7-5.8): reducd mass: phononfreq= %g (int units) = %g [SI, Hz] = %g [meV]",
            phononfreq,phononfreq*ω /(2*pi), phononfreq*hbar*ω *1000/q)

        @printf("\n\t Schultz1959: electronfreq= %g (int units) = %g [SI, Hz] = %g [meV]",
            sqrt(k/1),sqrt(k/1)*ω /(2*pi), sqrt(k/1)*hbar*ω *1000/q)
        @printf("\n\t Schultz1959: combinedfreq= %g (int units) = %g [SI, Hz] = %g [meV]",
            sqrt(k/(1+M)),sqrt(k/(1+M))*ω /(2*pi), sqrt(k/(1+M))*hbar*ω *1000/q)

        # Devreese1972: 10.1103/PhysRevB.5.2367
        # p.2371, RHS.
        @printf("\n Devreese1972: (Large Alpha) Franck-Condon frequency = %.2f", 4*α ^2/(9*pi))
        end

        # F(v,w,β,α)=-(A(v,w,β)+B(v,w,β,α)+C(v,w,β)) #(62a) - Hellwarth 1999
        if (T>0)
            @printf("\n Polaron Free Energy: A= %f B= %f C= %f F= %f",A(v,w,βred),B(v,w,βred,α),C(v,w,βred),F(v,w,βred,α))
            @printf("\t = %f meV",1000.0 * F(v,w,βred,α) * ħ*ω  / q)
            append!(p.A,A(v,w,βred))
            append!(p.B,B(v,w,βred,α))
            append!(p.C,C(v,w,βred))
            append!(p.F,F(v,w,βred,α))
        else # Athermal case; Enthalpy
            @printf("\n Polaron Enthalpy: F= %f = %f meV \n",F(v,w,α), 1_000*F(v,w,α) * ħ*ω/q)

            return # return early, as if T=0, all mobility theories = infinite / fall over
        end

        # FHIP
        #    - low-T mobility, final result of Feynman1962
        # [1.60] in Devreese2016 page 36; 6th Edition of Frohlich polaron notes (ArXiv)
        # I believe here β is in SI (expanded) units
        @printf("\n Polaron Mobility theories:")
        μ=(w/v)^3 * (3*q)/(4*mb*ħ*ω^2*α*β) * exp(ħ*ω*β) * exp((v^2-w^2)/(w^2*v))
        @printf("\n\tμ(FHIP)= %f m^2/Vs \t= %.2f cm^2/Vs",μ,μ*100^2)
        append!(p.T,T)
        append!(p.FHIPμ,μ*100^2)

        # Kadanoff
        #     - low-T mobility, constructed around Boltzmann equation.
        #     - Adds factor of 3/(2*beta) c.f. FHIP, correcting phonon emission behaviour
        # [1.61] in Devreese2016 - Kadanoff's Boltzmann eqn derived mob
        μ=(w/v)^3 * (q)/(2*mb*ω*α) * exp(ħ*ω*β) * exp((v^2-w^2)/(w^2*v))
        @printf("\n\tμ(Kadanoff,via Devreese2016)= %f m^2/Vs \t= %.2f cm^2/Vs",μ,μ*100^2)

        append!(p.Kμ,μ*100^2)

        ######
        # OK, now deep-diving into Kadanoff1963 itself to extract Boltzmann equation components
        # Particularly right-hand-side of page 1367
        #
        # From 1963 Kadanoff, (9), define eqm. number of phonons (just from T and phonon omega)
        Nbar=(exp(βred) - 1)^-1
        @printf("\n\t\tEqm. Phonon. pop. Nbar: %f ",Nbar)
        if (verbose)
            @printf("\n\texp(Bred): %f exp(-Bred): %f exp(Bred)-1: %f",exp(βred),exp(-βred),exp(βred)-1)
        end
        Nbar=exp(-βred)
        #Note - this is only way to get Kadanoff1963 to be self-consistent with
        #FHIP, and later statements (Devreese) of the Kadanoff mobility.
        #It suggests that Kadanoff used the wrong identy for Nbar in 23(b) for
        #the Gamma0 function, and should have used a version with the -1 to
        #account for Bose / phonon statistics

        #myv=sqrt(k*(1+M)/M) # cross-check maths between different papers
        #@printf("\nv: %f myv: %f\n",v,myv)

        # Between 23 and 24 in Kadanoff 1963, for small momenta skip intergration --> Gamma0
        Gamma0=2*α * Nbar * (M+1)^(1/2) * exp(-M/v)
        Gamma0*=ω  #* (ω *hbar) # Kadanoff 1963 uses hbar=omega=mb=1 units
            # Factor of omega to get it as a rate relative to phonon frequency
            # Factor of omega*hbar to get it as a rate per energy window
        μ=q/( mb*(M+1) * Gamma0 ) #(25) Kadanoff 1963, with SI effective mass

        if (verbose) # these are cross-checks
            @printf("\n\tμ(Kadanoff1963 [Eqn. 25]) = %f m^2/Vs \t = %.2f cm^2/Vs",μ,μ*100^2)
            @printf("\n\t\t Eqm. Phonon. pop. Nbar: %f ",Nbar)
        end

        @printf("\n\t\tGamma0 = %g rad/s = %g /s ",
                Gamma0, Gamma0/(2*pi))
        @printf(" \n\t\tTau=1/Gamma0 = %g s = %f ps",
            2*pi/Gamma0, 2*pi*1E12/Gamma0)
        Eloss=hbar*ω * Gamma0/(2*pi) # Simply Energy * Rate
        @printf("\n\t\tEnergy Loss = %g J/s = %g meV/ps",Eloss,Eloss * 1E3 / (q*1E12) )
        append!(p.Tau, 2*pi*1E12/Gamma0) # Boosted into ps ?

        # Hellwarth1999 - directly do contour integration in Feynman1962, for
        # finite temperature DC mobility
        # Hellwarth1999 Eqn (2) and (1) - These are going back to the general
        # (pre low-T limit) formulas in Feynman1962.  to evaluate these, you
        # need to do the explicit contour integration to get the polaron
        # self-energy
        R=(v^2-w^2)/(w^2*v) # inline, page 300 just after Eqn (2)

        #b=R*βred/sinh(b*βred*v/2) # This self-references b! What on Earth?
        # OK! I now understand that there is a typo in Hellwarth1999 and
        # Biaggio1997. They've introduced a spurious b on the R.H.S. compared to
        # the original, Feynman1962:
        b=R*βred/sinh(βred*v/2) # Feynman1962 version; page 1010, Eqn (47b)

        a=sqrt( (βred/2)^2 + R*βred*coth(βred*v/2))
        k(u,a,b,v) = (u^2+a^2-b*cos(v*u))^(-3/2)*cos(u) # integrand in (2)
        K=quadgk(u->k(u,a,b,v),0,Inf)[1] # numerical quadrature integration of (2)

        #Right-hand-side of Eqn 1 in Hellwarth 1999 // Eqn (4) in Baggio1997
        RHS= α/(3*sqrt(π)) * βred^(5/2) / sinh(βred/2) * (v^3/w^3) * K
        μ=RHS^-1 * (q)/(ω*mb)
        @printf("\n\tμ(Hellwarth1999)= %f m^2/Vs \t= %.2f cm^2/Vs",μ,μ*100^2)
        append!(p.Hμ,μ*100^2)

        if (verbose)
        # Hellwarth1999/Biaggio1997, b=0 version... 'Setting b=0 makes less than 0.1% error'
        # So let's test this
        R=(v^2-w^2)/(w^2*v) # inline, page 300 just after Eqn (2)
        b=0
        a=sqrt( (βred/2)^2 + R*βred*coth(βred*v/2))
        #k(u,a,b,v) = (u^2+a^2-b*cos(v*u))^(-3/2)*cos(u) # integrand in (2)
        K=quadgk(u->k(u,a,b,v),0,Inf)[1] # numerical quadrature integration of (2)

        #Right-hand-side of Eqn 1 in Hellwarth 1999 // Eqn (4) in Baggio1997
        RHS= α/(3*sqrt(π)) * βred^(5/2) / sinh(βred/2) * (v^3/w^3) * K
        μ=RHS^-1 * (q)/(ω*mb)
        @printf("\n\tμ(Hellwarth1999,b=0)= %f m^2/Vs \t= %.2f cm^2/Vs",μ,μ*100^2)
        @printf("\n\tError due to b=0; %f",(100^2*μ-p.Hμ[length(p.Hμ)])/(100^2*μ))
        #append!(Hμs,μ*100^2)
        end
        @printf("\n") # blank line at end of spiel.

        # Recycle previous variation results (v,w) as next guess
        initial=[v,w] # Caution! Might cause weird sticking in local minima

        append!(p.v,v)
        append!(p.w,w)

    end

    return(p)
end

"""
    savepolaron(fileprefix, p::Polaron)

Saves data from polaron 'p' into file "fileprefix".
This is a simple space-delimited text file, with each entry a separate temperature, for plotting with Gnuplot or similar.

Structure of file is written to the header:
# Ts, βreds, Kμs, Hμs, FHIPμs, vs, ws, ks, Ms, As, Bs, Cs, Fs, Taus, rfsis
# 1    2     3    4     5      6   7   8   9  10  11  12  13    14     15
"""
function savepolaron(fileprefix, p::Polaron)
    println("Saving data to $fileprefix.dat ...")
    f=open("$fileprefix.dat","w")

    @printf(f,"# %s\n",fileprefix) # put name / material at header
    @printf(f,"# Params in SI: ω =%g mb=%g \n",p.ω[1] ,p.mb[1])
    @printf(f,"# Alpha parameter: α = %f  \n",p.α[1] )

    @printf(f,"# Ts, βreds, Kμs, Hμs, FHIPμs, vs, ws, ks, Ms, As, Bs, Cs, Fs, Taus, rfsis\n")
    @printf(f,"#  1    2     3    4     5      6   7   8   9  10  11  12  13    14     15\n") # columns for GNUPLOT etc.

    for i in 1:length(p.T)
        @printf(f,"%d %03f %g %g %g %g %g %g %g %g %g %g %g %g %g \n",
        p.T[i], p.βred[i], p.Kμ[i], p.Hμ[i], p.FHIPμ[i],
        p.v[i], p.w[i],
        p.k[i], p.M[i], p.A[i], p.B[i], p.C[i], p.F[i],
        p.Tau[i], p.rfsi[i])
    end
    close(f)
end


text:
module OnlinePackage

import Base: Generator
import HTTP
import JSON
import JSON: json
import LibGit2
import Base64: base64encode

abstract type Remote end
const user_agent = "OnlinePackage/0.0.1"

include("talk_to.jl")
include("github.jl")
include("travis.jl")
include("travis_token.jl")
include("ssh_keygen.jl")
include("configure.jl")
include("copy.jl")
include("generate.jl")

end


text:
module SugarKelp
using RecursiveArrayTools, DiffEqBase, OrdinaryDiffEq, Roots, Interpolations, DataFrames
# Input parameters
# t: time (days)
# irr: irradiance (micro mol photons / m^2 / s)
# temp: Temperature (deg C)
# ex_n: External (to kelp) nitrate concentration (micro mol / L)
# u: Current speed (relative to kelp) (m/s)

# Variables
# a: Frond area of individual kelp(dm^2)
# n: Nitrogen reserve relative to dry weight (gN/(g sw))
# c: Carbon reserve relative to dry weight (gC/(g sw))

"""
    SugarKelp.eval_μ(a,n,c,temp,λ)

Solves Equation 2, the specific growth rate.

Parameters:
- `a`: area /dm^2
- `n`: nitrate reserve /gN/gSW
- `c`: carbon reserve /gC/gSW
- `temp`: temperature /°C
- `λ`: normalised change in day length

Returns: specific growth rate

Notes:
- Names eval_μ because it looks better in SugarKelp.equations to have μ and you can't reuse it
"""
eval_μ(a,n,c,temp,λ) = f_area(a) * f_temp(temp) * f_photo(λ) * min(1 - N_min / n, 1 - C_min / c)

"""
    SugarKelp.f_area(a)

Solves Equation 3, the effect of area on growth

Parameters: `a`, area /dm^2

Returns: effect of area on growth
"""
f_area(a) = m_1 * exp(-(a / A_0)^2) + m_2

"""
    SugarKelp.f_temp(temp)

Solves Equation 4, the effect of temperature on growth

Parameters: `temp`: temperature /°C

Returns: effect of temperature on growth
"""
function f_temp(temp)
    if -1.8 <= temp < 10
        return 0.08 * temp + 0.2
    elseif 10 <= temp <= 15
        return 1
    elseif 15 < temp <= 19
        return 19 / 4 - temp / 4
    elseif temp > 19
        return 0
    else
        return 0
    end
end

"""
    SugarKelp.f_photo(λ)

Solves Equation 5, the seasonal influence on growth rate

Parameters: `λ`: normalised change in day length

Returns: seasonal influence on growth
"""
f_photo(λ) = a_1 * (1 + sign(λ) * abs(λ)^.5) + a_2

"""
    SugarKelp.ν(a)

Solves Equation 6, the specific frond erosion rate.

Parameters: `a`: area /dm^2

Returns: specific frond erosion rate
"""
ν(a) = 1e-6 * exp(ϵ * a) / (1 + 1e-6 * (exp(ϵ * a) - 1))

"""
    SugarKelp.eval_j(ex_n,n,u)

Solves Equation 8, the specific nitrate uptake rate.

Parameters:
- `ex_n`: external nitrate concentration /mmol/m^3
- `n`: nitrogen reserve /gN/gSW
- `u`: water speed /m/s

Returns: specific nitrate uptake rate

Notes:
- Names eval_j because it looks better in SugarKelp.equations to have j and you can't reuse it
"""
eval_j(ex_n,n,u) = J_max * (ex_n / (K_X + ex_n)) * ((N_max - n) / (N_max - N_min)) * (1 - exp(-u / U_0p65))

"""
    SugarKelp.p(temp,irr)

Solves Equation 10, the gross photosynthesis.

Parameters:
- `temp`: temperature /°C
- `irr`: irradiance /mol photons/m²/day

Returns: gross photosynthesis function
"""
function p(temp, irr)
    p_max = P_1 * exp(T_AP / T_P1 - T_AP / (temp + 273.15)) / (1 + exp(T_APL / (temp + 273.15) - T_APL / T_PL) + exp(T_APH / T_PH - T_APH / (temp + 273.15)))
    β_func(x) = p_max - (α * I_sat / log(1 + α / x)) * (α / (α + x)) * (x / (α + x))^(x / α)
    β = find_zero(β_func, (0, 0.1), Bisection())
    p_s = α * I_sat / log(1 + α / β)
    return p_s * (1 - exp(-α * irr / p_s)) * exp(-β * irr / p_s)
end

"""
    SugarKelp.r(temp,μ,j,resp_model)

Solves Equation 14 (2012) if resp_model=1, or Equation 2 (2013) if resp_model=2

Parameters:
- `temp`: temperature /°C
- `μ`: specific growth rate
- `j`: specific nitrate uptake rate
- `resp_model`: choice of resparation model (see description)

Returns: respiration function
"""
function r(temp, μ, j, resp_model)
    if resp_model == 1
        return R_1 * exp(T_AR / T_R1 - T_AR / (temp + 273.15)) # temperature dependent respiration
    elseif resp_model == 2
        return (R_A * (μ / μ_max + j / J_max) + R_B) * exp(T_AR / T_R1 - T_AR / (temp + 273.15))
    end
end

e(c) = 1 - exp(γ * (C_min - c))

"""
    SugarKelp.equations(t, a, n, c, u, temp, irr, ex_n, λ, resp_model, dt)

Solves the papers main equations.

Parameters:
- `a`: area /dm^2
- `n`: nitrate reserve /gN/gSW
- `c`: carbon reserve /gC/gSW
- `u`: water speed /m/s
- `temp`: temperature /°C
- `irr`: irradiance / mol photons/m²/day
- `ex_n`: external nitrate concentration /mmol m^3
- `λ`: normalised change in day length
- `resp_model`: choice of resparation model (see SugarKelp.r)
- `dt`: timestep length /days

Returns:
- `da`: results of equation 1, the area rate
- `dn`: results of equation 7, the nitrogen rate
- `dc`: results of equation 9, the carbon rate
- `j`: the specific nitrate uptake rate
"""
function equations(a, n, c, u, temp, irr, ex_n, λ, resp_model, dt)
    μ = eval_μ(a, n, c, temp, λ)
    j = eval_j(ex_n, n, u)

    # Equation 1
    da = (μ - ν(a)) * a
    # Equation 7
    dn = j / K_A - μ * (n + N_struct)
    # Equation 9
    dc = (p(temp, irr) * (1 - e(c)) - r(temp, μ, j, resp_model)) / K_A - μ * (c + C_struct)

    # "Extreme carbon limitation"
    if c + dc < C_min
        da -= (a * (C_min - c) / C_struct) * .5 * dt
        dc = (C_min - c) * .5 * dt
    end

    return da, dn, dc, j
end

"""
    SugarKelp.solver!(y, params, t)

Interface for OrdinaryDiffEq library, extracts current enviromental variable values and solves equations for each timestep.

Parameters:
- `y`: the current state of the system as a vector (area, nitrate reserve, carbon reserve).
- `params`: the variable parameters of the model:
    - `u_arr`: interpolation object of the water speed in time
    - `temp_arr`: interpolation object of the temperature in time
    - `irr_arr`: interpolation object of the irradianec in time
    - `ex_n_arr`: interpolation object of the external nitrate concentration in time
    - `λ_arr`: array of the normalised change in day length
    - `resp_model`: the choice of respiration model, 1 is the origional from the 2012 paper and 2 is the modified version in 2013
    - `dt`: the time step length, this is important as it is used in the "extreme carbon limit" part of the equations, see NB.
- `t`: the current time (with respect to the time in the interpolations)

Returns: array of da,dn,dc,j*a - the rate of a/n/c and the nitrate uptake rate

Note:
These equations **must** be solved with an algorithm with fixed time steps and known, constant sub timestep lengths.
This is because the "extreme carbon limit" element is only implientable with known and fixed timesteps as the value of a
must be changed to a particular value rather than changing the deriviative. This can only be done (within the framework of
the ODE library) by setting the derivitive to (X(next)-X(old))/dt
"""
function solver!(y::Vector{Float64}, params, t::Float64)
    a, n, c = y

    if a > 0
        u_arr, temp_arr, irr_arr, ex_n_arr, λ_arr, resp_model, dt = params

        u = u_arr(t)::Float64
        temp = temp_arr(t)::Float64
        irr = irr_arr(t)::Float64
        ex_n = ex_n_arr(t)::Float64

        d = trunc(Int, mod(floor(t), 365) + 1)
        λ = λ_arr[d]

        da, dn, dc, j = equations(a, n, c, u, temp, irr, ex_n, λ, resp_model, dt)

    else
        da, dn, dc, j = 0, 0, 0, 0
    end
    return (vcat(da, dn, dc, j * a))
end

"""
    SugarKelp.solve(t_i, nd, u, temp, irr, ex_n, lat, a_0, n_0, c_0, params="src/parameters/origional.jl", resp_model=1, dt=1, dataframe=true)

Solves the model for some set of parametetrs and returns the ODE library solution as well as a dataframe of the useful results.

Parameters:
- `t_i`: th estart time (in days since the start of the interpolation objects "day zero")
- `nd`: the number of days to run for
- `u`: interpolation object (IO) of water speed
- `temp`: IO of temperature
- `ex_n`: IO of external nitrate concentration
- `lat`: latitude, relivant for the change of day length
- `a_0`: initial area
- `n_0`: initial nitrogen reserve (gN/gSW)
- `c_0`: initial carbon reserve (gC/gSW)
- `params`: string of the path to a parameters file, defaults to the 2012 values. Also supplied is 2013 in src/parameters/2013.jl or you can copy and vary them
- `resp_model`: choice of respiration model, 1 (default) uses the 2012 version and 2 uses the modifcations from the 2013 paper
- `dt`: the time step size to use (see equations! note), default is 1 day (seems small enough)
- `dataframe`: output as a dataframe, default to true. Alternative is an array (faster)
- `λ_arr`: array of normalised change of day length, defaults to nothing which generates the default one

Returns:
- `solution`: the ODE library solution
- `results`: dataframe or array of area/nitrogen reserve/carbon reserve/total nitrate update. All others useful quantities can be easily derived.
"""
function solve(t_i, nd, u, temp, irr, ex_n, lat, a_0, n_0, c_0, params="../src/parameters/origional.jl", resp_model=1, dt=1, dataframe=true, λ_arr=nothing)
    include(params)
    if λ_arr==nothing
        λ_arr = gen_λ(lat)
    end
    params = (u, temp, irr, ex_n, λ_arr, resp_model, dt)

    y_0 = vcat(a_0, n_0, c_0, 0)

    solver = ODEProblem(solver!, y_0, (t_i, t_i + nd), params)
    solution = OrdinaryDiffEq.solve(solver, RK4(), dt=dt, adaptive=false)# Please keep the RK4 algorithm otherwise the extreme caron limit needs to be changed

    if dataframe == true
        results = DataFrame(area=[], nitrogen=[], carbon=[], gross_nitrate=[], time=[])
        for (ind, val) in enumerate(solution.u)
            push!(val, solution.t[ind])
            push!(results, (val))
        end
        return(solution, results)
    else
        return vcat(transpose.(solution.u)...)
    end
end

"""
    SugarKelp.solvegrid(t_i, nd, a_0, n_0, c_0, arr_lon, arr_lat, arr_dep, arr_time, no3, temp, u, par_data, kd_data, att = nothing, params="src/parameters/origional.jl", resp_model=1, dt=1, progress=true)
Solve the model for a (spacially) fixed grid of inputs.

Parameters:
- `t_i`: start day (w.r.t. t=0 in time)
- `nd`: number of days to run for
- `a_0`: initial area/dm^2
- `n_0`: initial nitrogen reserve/gN/gSW
- `c_0`: initial carbon reserve/gC/gSW
- `arr_lon`: longitudes for no3 and temp
- `arr_lat`: latitudes for no3 and temp
- `arr_dep`: depths for no3 and temp
- `arr_time`: time for no3 and temp
- `no3`: array of no3 concentration in lon,lat,depth,time
- `temp`: array of temp concentration in lon,lat,depth,time
- `u`: array of water speed in lon,lat,depth,time
- `par_data`: array of:
    - `par values in lon,lat,time
    - `corresponding time
    - `par fill value
- `kd_data`: array of:
    - kd values in lon,lat,time
    - corresponding time
    - kd fill value
- `att`: array of light attenuation coefficients (PAR(z)=PAR(z=0)*att) in lon,lat,depth,time. Defaults to nothing and kd is used instead
- `params`: string of the path to a parameters file, defaults to the 2012 values. Also supplied is 2013 in src/parameters/2013.jl or you can copy and vary them
- `resp_model`: choice of respiration model, 1 (default) uses the 2012 version and 2 uses the modifcations from the 2013 paper
- `dt`: the time step size to use (see equations! note), default is 1 day (seems small enough)
- `progress`: option to update progress at each level (when multithreading not accurate but useful), defaults to true

Returns: results as an array of (area/nitrogen/carbon/nitrate update, lon, lat, depth, time)

Notes:
par and kd need their own time coordinates and fill value because they come from satelite observation which are
temporally sparse so need to be checked and interpolated in time for each point. On the other hand temp and no3
(that I'm using) are from Copurnicus' models so if a point has a value at some time it will at all times.

no3,temp and u need to be of the same shape and size and with the values corresponding to the same position/time.
"""
function solvegrid(t_i::Float64, nd::Int, a_0::Float64, n_0::Float64, c_0::Float64, arr_lon, arr_lat, arr_dep, arr_time, no3::Array{Float64,4}, temp::Array{Float64,4}, u::Array{Float64,4}, par_data, kd_data, att=nothing, params::String="../src/parameters/origional.jl", resp_model::Int=1, dt=1, progress::Bool=true)
    # Would like to annotate type for the others but for some reason they making tuples of "Number" doesn't isn't satisfied
    par, par_t, par_fill = par_data;kd, kd_t, kd_fill = kd_data

    all_results::Array{Float64,5} = repeat([NaN], 4, length(arr_lon), length(arr_lat), length(arr_dep), round(Int, nd / dt + 1));
    points::Array{Vector{Int64},3} = collect.(Iterators.product([1:length(arr_lon);], [1:length(arr_lat);], [1:length(arr_dep);]));

    Threads.@threads for (i, j, k) in points
        if (progress == true) & (i == 1) & (j == 1)
            @info("At level $k")
        end

        lat = arr_lat[j];depth = arr_dep[k]

        no3_vals = no3[i,j,k,:]
        temp_vals = temp[i,j,k,:]
        u_vals = u[i,j,k,:]
        if (!isnan(no3_vals[1])) & (!isnan(temp_vals[1])) & (!isnan(u_vals[1])) & all(n > 0 for n in no3_vals)
            no3_itp = Interpolations.LinearInterpolation(arr_time, no3_vals)
            temp_itp = Interpolations.LinearInterpolation(arr_time, temp_vals)
            u_itp = Interpolations.LinearInterpolation(arr_time, u_vals)

            par_vals_raw = par[i,j,:]
            par_vals, par_t_vals = extract_valid(par_vals_raw, par_t, par_fill)

            if att == nothing
                kd_vals_raw = kd[i,j,:]
                kd_vals, kd_t_vals = extract_valid(kd_vals_raw, kd_t, kd_fill)
            else
                att_vals = att[i,j,k,:]
            end

            if (length(par_vals) > 6)# & (att!=nothing) | (length(kd_vals) > 6)
                par_itp = Interpolations.LinearInterpolation(par_t_vals, par_vals, extrapolation_bc=Flat())
                if att == nothing
                    kd_itp = Interpolations.LinearInterpolation(kd_t_vals, kd_vals, extrapolation_bc=Flat())
                    irr_itp = Interpolations.LinearInterpolation(arr_time, par_itp.(arr_time) .* exp.(-kd_itp.(arr_time) .* depth), extrapolation_bc=Flat())
                else
                    irr_itp = Interpolations.LinearInterpolation(arr_time, par_itp.(arr_time) .* att_vals, extrapolation_bc=Flat())
                end
                solution = SugarKelp.solve(t_i, nd, u_itp, temp_itp, irr_itp, no3_itp, lat, a_0, n_0, c_0, params, resp_model, dt, false);

                all_results[:,i,j,k,1:size(solution)[1]] = reshape(solution', size(solution)[2], 1, 1, 1, size(solution)[1])
    end
        end
    end
    return all_results
end

"""
    SugarKelp.extract_valid(raw,raw_time,fill)

Extracts the valid values from an array by checking against a fill value and returns the valids and corresponding time.

Parameters:
- `raw`: the array to check
- `raw_time`: the corresponding time array
- `fill`: the fill value to check against

Returns:
- `vals`: the filtered values
- `time`: corresponding times
"""
function extract_valid(raw, raw_time, fill)
    vals, times = [], []
    for (ind, val) in enumerate(raw)
        if (val != fill) & (!isnan(val))
            push!(vals, val)
            push!(times, raw_time[ind])
        end
    end
    return (vals, times)
end

"""
    SugarKelp.gen_λ(lat)

Generates λ as described in Model desctiptions/Main equations/Photoperiodic effect (page 763) in the origional paper.

Parameter: `lat`, the latitude

Returns: array of normalised change in day length
"""
function gen_λ(lat)
    θ = 0.2163108 .+ 2 .* atan.(0.9671396 .* tan.(0.00860 .* ([1:365;] .- 186)))
    δ = asin.(0.39795 .* cos.(θ))
    p = 0.8333
    d = 24 .- (24 / pi) .* acos.((sin(p * pi / 180) .+ sin(lat * pi / 180) .* sin.(δ)) ./ (cos(lat * pi / 180) .* cos.(δ)))
    λ = diff(d)
    push!(λ, λ[end])
    return λ ./ findmax(λ)[1]
end
"""
    SugarKelp.get_ind(val, list, tol)
Function that finds the index in the list with the closest value to val. Error is thrown if no result is within tollerance, tol.

Parameters:
- `val`: the value searching for
- `list`: the list to search
- `tol`: tollerance of search

Returns:
- `index of closest valuea
"""
function get_ind(val, list, tol)
    result = findmin(abs.(list .- val))
    ind = result[2]
    if result[1] > tol
        closest = list[result[2]]
        throw("No indicie could be found within tollerance for the requested value, requested was $val, closes was $closest at $ind")
    else
        return ind
    end
end

"""
    SugarKelp.interp_deps(arr, origional_depths, desired_depths, invalid_val)
Function to interpolate a 4D array in the last dimension. Useful for lineaising the depth steps of the arrays
from Copurnicus as they have increasingly corse step sizes.

Parameters:
- `arr`: the array to interpolate where the 3rd dimension is the one to interpolate along
- `origional_depths`: the values of the 3rd dimensions coordinates (the origional depths of the data)
- `desired_depths`: the desired coordinates/depths
- `invalid_val`: the fill value to replace with NaN in the array (as we are searching it anyway)

Returns: new array with interpolated 3rd dimension
"""
function interp_deps(arr, origional_depths, desired_depths, invalid_val)
    arr_size = size(arr)
    new_arr = repeat([NaN], arr_size[1], arr_size[2], length(desired_depths), arr_size[4]);

    points = collect.(Iterators.product([1:arr_size[1];], [1:arr_size[2];], [1:arr_size[4];]));
    Threads.@threads for (i, j, k) in points
        dep_arr = arr[i,j,:,k]
        usable = [];
        finished = false
        for (ind, val) in enumerate(dep_arr)
            if (finished == false) & (!isnan(val)) & (round(val) != round(invalid_val))
                push!(usable, val)
            else
                finished = true
            end
        end
        if length(usable) > 1
            deps = origional_depths[1:length(usable)]
            dep_itp = Interpolations.LinearInterpolation(deps, usable, extrapolation_bc=Flat())

            extent = findmin(abs.(desired_depths .- deps[end]))[2]
            if extent != length(desired_depths)
                search_deps = desired_depths[1:extent - 1]
            else
                search_deps = desired_depths
            end
            new_arr[i,j,1:length(search_deps),k] = dep_itp.(search_deps)
        end
    end

    return new_arr
end
end # module


text:
include("euler/euler.jl")
using .Numbers: count_divisors

function compute(n::Int)::Int
    i, triangle = 1, 1
    while count_divisors(triangle) < n
        i += 1
        triangle += i
    end
    return triangle
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
# Rational Quadratic ARD Covariance Function

"""
    RQArd <: StationaryARD{WeightedSqEuclidean}

ARD Rational Quadratic kernel (covariance)
```math
k(x,x') = σ²(1 + (x - x')ᵀL⁻²(x - x')/(2α))^{-α}
```
with length scale ``ℓ = (ℓ₁, ℓ₂, …)``, signal standard deviation ``σ``, and shape parameter
``α`` where ``L = diag(ℓ₁, ℓ₂, …)``.
"""
mutable struct RQArd{T<:Real} <: StationaryARD{WeightedSqEuclidean}
    "Inverse squared length scale"
    iℓ2::Vector{T}
    "Signal variance"
    σ2::T
    "Shape parameter"
    α::T
    "Priors for kernel parameters"
    priors::Array
end

"""
Rational Quadratic ARD Covariance Function

    RQArd(ll::Vector{Real}, lσ::Real, lα::Real)

# Arguments
  - `ll::Vector{Real}`: vector of length scales (given on log scale)
  - `lσ::Real`: signal standard deviation (given on log scale)
  - `lα::Real`: shape parameter (given on log scale)
"""
RQArd(ll::Vector{T}, lσ::T, lα::T) where T = RQArd{T}(exp.(-2 .* ll), exp(2 * lσ), exp(lα), [])

function set_params!(rq::RQArd, hyp::AbstractVector)
    length(hyp) == num_params(rq) || throw(ArgumentError("RQArd kernel has $(num_params(rq_ard)) parameters"))
    @views @. rq.iℓ2 = exp(-2 * hyp[1:(end-2)])
    rq.σ2 = exp(2 * hyp[end-1])
    rq.α = exp(hyp[end])
end

get_params(rq::RQArd) = [-log.(rq.iℓ2) / 2; log(rq.σ2) / 2; log(rq.α)]
get_param_names(rq::RQArd) = [get_param_names(rq.iℓ2, :ll); :lσ; :lα]
num_params(rq::RQArd) = length(rq.iℓ2) + 2

cov(rq::RQArd,r::Number) = rq.σ2*(1+0.5*r/rq.α)^(-rq.α)

@inline dk_dll(rq::RQArd, r::Real, wdiffp::Real) =
    rq.σ2 * wdiffp * (1 + r / (2 * rq.α))^(-rq.α - 1)
@inline function dk_dlα(rq::RQArd, r::Real)
    part = (1 + r / (2 * rq.α))
    return rq.σ2 * part^(-rq.α) * (r / (2 * part) - rq.α * log(part))
end
@inline function dKij_dθp(rq::RQArd, X1::AbstractMatrix, X2::AbstractMatrix, i::Int, j::Int, p::Int, dim::Int)
    if p <= dim
        return dk_dll(rq, distij(metric(rq),X1,X2,i,j,dim), distijk(metric(rq),X1,X2,i,j,p))
    elseif p==dim+1
        return dk_dlσ(rq, distij(metric(rq),X1,X2,i,j,dim))
    else
        return dk_dlα(rq, distij(metric(rq),X1,X2,i,j,dim))
    end
end
@inline function dKij_dθp(rq::RQArd, X1::AbstractMatrix, X2::AbstractMatrix, data::StationaryARDData,
                          i::Int, j::Int, p::Int, dim::Int)
    return dKij_dθp(rq,X1,X2,i,j,p,dim)
end


text:
using FileIO
using Printf
using PyPlot

using ParallelTreeSampling

include("utils.jl")

# Parameters to be set for plotting.
mc_data = ["data/gridworld_large_baseline_0427-2140_1000000_$(trial).jld2" for trial in 1:3]
# mc_data = []
is_data = ["data/gridworld_large_mcts_0427-2221_500000_$(trial).jld2" for trial in 1:3]
# is_data = []
max_num_mc = -1; max_num_is = -1;
delta_n_mc = 10_000; delta_n_is = 10_000;
plot_hist = true; plot_est = true;
output_dir = homedir()

# Constants to be updated as needed.
samples_key = "risks:"; weights_key = "weights:";

# Read the jld2 files.
function load_jld2(jld2_filename)
    d = load(jld2_filename)
    @assert haskey(d, samples_key)
    return (haskey(d, weights_key) ?
        (d[samples_key], d[weights_key]) : d[samples_key])
end
n_mc_samples = [load_jld2(mc_datum)[1:(max_num_mc < 0 ? end : max_num_mc)]
                for mc_datum in mc_data]
n_is_samples = [load_jld2(is_datum)[1:(max_num_is < 0 ? end : max_num_is)]
                for is_datum in is_data]

# Returns a lambda for computing metrics used for plotting.
#
# Valid metric types include :mean, :var, :cvar, :worst.
function estimate_fn(metric_type, alpha, with_weights=false)
    if with_weights
        fn = function(samples, weights)
            metrics = eval_metrics(samples; weights=exp.(weights), alpha=alpha)
            return getproperty(metrics, metric_type)
        end
    else
        fn = function(samples)
            metrics = eval_metrics(samples; weights=nothing, alpha=alpha)
            return getproperty(metrics, metric_type)
        end
    end
    return fn
end

# Plot convergence graphs.
alphas = [1e-4]
if plot_est
    figure(figsize=(9.0, 6.0))
    if !isempty(n_mc_samples)
        for alpha in alphas
            alpha_str = @sprintf("%.1e", alpha)
            plot_estimates(n_mc_samples, estimate_fn(:var, alpha),
                           delta_n_mc, "MC-VaR-$(alpha_str)")
            plot_estimates(n_mc_samples, estimate_fn(:cvar, alpha),
                           delta_n_mc, "MC-CVaR-$(alpha_str)")
        end
    end
    if !isempty(n_is_samples)
        for alpha in alphas
            alpha_str = @sprintf("%.1e", alpha)
            plot_estimates(n_is_samples, estimate_fn(:var, alpha, true),
                           delta_n_is, "IS-VaR-$(alpha_str)")
            plot_estimates(n_is_samples, estimate_fn(:cvar, alpha, true),
                           delta_n_is, "IS-CVaR-$(alpha_str)")
        end
    end
    xlabel("no. samples"); ylabel("estimates"); legend();
    if !isempty(output_dir)
        PyPlot.savefig(joinpath(output_dir, "estimates.png"))
    end
end

# Plot histograms.
if plot_hist
    figure(figsize=(9.0, 6.0))
    if !isempty(n_mc_samples)
        samples = first(n_mc_samples)
        plot_histogram(samples, bins=50, label="MC")
    end
    if !isempty(n_is_samples)
        samples, weights = first(n_is_samples)
        plot_histogram(samples, bins=50, label="IS-uw")
        plot_histogram(samples, weights, bins=50, label="IS-w")
    end
    ylim(bottom=0.0); xlabel("costs"); ylabel("density"); legend();
    if !isempty(output_dir)
        PyPlot.savefig(joinpath(output_dir, "histogram.png"))
    end
end


text:
"""
Struct for Active Brownian particle system
"""
mutable struct ActiveBrownianSystem <: AbstractActiveSystem
  dim::Int64
  v0::Float64
  t::Float64
  x::Vector{Float64}
  f::Vector{Float64}
  thermostat::ActiveLangevin
  function ActiveBrownianSystem(dim::Int64,v0::Float64,t::Float64,x::Vector{Float64},f::Vector{Float64},thermostat::ActiveLangevin)
    new(dim,v0,t,x,f,thermostat)
  end
end

function ActiveBrownianSystem(v0::Float64,model::AbstractPotential,thermostat::ActiveLangevin)
  dim::Int64 = getdimensionality(model)
  x::Vector{Float64} = zeros(dim+1)
  f::Vector{Float64} = zeros(dim+1)
  ActiveBrownianSystem(dim+1,v0,0.0,x,f,thermostat)
end

# TODO idk yet, but this needs to be fixed to allow the different gammas
function addforce!(dt::Float64,sys::ActiveBrownianSystem,model::AbstractSinglePotential)
  @. @views sys.x[1:sys.dim-1] .+= dt .* model.f[1:sys.dim-1] ./ sys.thermostat.gammat
  if model.dim == sys.dim
    sys.x[end] += dt .* model.f[end] ./ sys.thermostat.gammar
  end
end


text:
using Match

input = joinpath(@__DIR__, "input")
lines = readlines(input)

DEBUG = false

function decode(lines)
    acc = []
    for line in lines
        @match line[1:3] begin
            "inp" => push!(acc, (:inp, line[5]))
            "add" => push!(acc, (:add, line[5], line[7:end]))
            "mul" => push!(acc, (:mul, line[5], line[7:end]))
            "div" => push!(acc, (:div, line[5], line[7:end]))
            "mod" => push!(acc, (:mod, line[5], line[7:end]))
            "eql" => push!(acc, (:eql, line[5], line[7:end]))
        end
    end
    return acc
end

function adder(registers, a, b)
    rs = copy(registers)
    if first(b) in ['w', 'x', 'y', 'z']
        rs[a] = rs[a] + rs[first(b)]
    else
        rs[a] = rs[a] + parse(Int, b)
    end
    return rs
end

function multiplier(registers, a, b)
    rs = copy(registers)
    if first(b) in ['w', 'x', 'y', 'z']
        rs[a] = rs[a] * rs[first(b)]
    else
        rs[a] = rs[a] * parse(Int, b)
    end
    return rs
end

function divider(registers, a, b)
    rs = copy(registers)
    if first(b) in ['w', 'x', 'y', 'z']
        rs[a] = floor(rs[a] / rs[first(b)])
    else
        rs[a] = floor(rs[a] / parse(Int, b))
    end
    return rs
end

function modder(registers, a, b)
    rs = copy(registers)
    if first(b) in ['w', 'x', 'y', 'z']
        rs[a] = rem(rs[a], rs[first(b)])
    else
        rs[a] = rem(rs[a], parse(Int, b))
    end
    return rs
end

function equaler(registers, a, b)
    rs = copy(registers)
    if first(b) in ['w', 'x', 'y', 'z']
        if rs[a] == rs[first(b)]
            rs[a] = 1
        else
            rs[a] = 0
        end
    else
        if rs[a] == parse(Int, b)
            rs[a] = 1
        else
            rs[a] = 0
        end
    end
    return rs
end

# generic vm. slow as shit. use runf below
function run(program, registers, input, io)
    idx = io
    for int in program
        @match first(int) begin
            :inp => begin
                _, register = int
                registers[register] = parse(Int, input[idx])
                idx += 1
                DEBUG && println("$int :: $io :: $registers")
            end
            :add => begin
                _, a, b = int
                registers = adder(registers, a, b)
                DEBUG && println("$int :: $io :: $registers")
            end
            :mul => begin
                _, a, b = int
                registers = multiplier(registers, a, b)
                DEBUG && println("$int :: $io :: $registers")
            end
            :div => begin
                _, a, b = int
                registers = divider(registers, a, b)
                DEBUG && println("$int :: $io :: $registers")
            end
            :mod => begin
                _, a, b = int
                registers = modder(registers, a, b)
                DEBUG && println("$int :: $io :: $registers")
            end
            :eql => begin
                _, a, b = int
                registers = equaler(registers, a, b)
                DEBUG && println("$int :: $io :: $registers")
            end
        end
    end
    return registers
end

# optimized version of the program. doesn't take a year to finish
function runf(program, rs, input, i)
    rs = copy(rs)
    w = parse(Int, input[i])
    mode = parse(Int, last(program[5]))
    check = parse(Int, last(program[6]))
    offset = parse(Int, last(program[16]))

    z = rs['z']
    x = rem(z, 26) + check
    z = Int(floor(z / mode))
    if x != w
        rs['z'] = 26 * z + w + offset
    else
        rs['z'] = z
    end
    return rs
end

PROGRAM = decode(lines)
INITIAL = Dict([('w', 0), ('x', 0), ('y', 0), ('z', 0)])

function optimize(f)
    candidates = Dict([([], 0)])

    # 14 inputs
    for i in 1:14
        baseline = minimum(values(candidates))
        program = PROGRAM[((18 * (i - 1)) + 1):(18 * i)]
        digits = map(string, 1:9)
        definitely = Dict()
        maybe = Dict()
        for (candidate, z) in collect(candidates)
            for digit in digits
                rs = copy(INITIAL)
                rs['z'] = z
                input = string(collect(candidate)...) * digit
                result = f(program, rs, input, i)
                nz = result['z']
                # @show((input, nz, baseline, nz < baseline))
                if nz < baseline
                    definitely[input] = nz
                elseif isempty(definitely)
                    maybe[input] = nz
                end
            end
        end
        if !isempty(definitely)
            candidates = definitely
        else
            candidates = maybe
        end
    end

    filter!(kv -> last(kv) == 0, candidates)
    return extrema(map(n -> parse(Int, n), collect(keys(candidates))))
end

# pass in run instead of runf to verify correctness
p2, p1 = optimize(runf)

println("-----------------------------------------------------------------------")
println("arithmetic logic unit -- part one :: $p1")
println("arithmetic logic unit -- part two :: $p2")
println("-----------------------------------------------------------------------")

@assert(p1 == 52926995971999)
@assert(p2 == 11811951311485)


text:
include("../common/PWGrid_v02.jl")
include("../common/ortho_gram_schmidt.jl")
include("../common/wrappers_fft.jl")

include("EnergiesT.jl")
include("PotentialsT.jl")
include("gen_dr.jl")
include("init_pot_harm_3d.jl")
include("op_K.jl")
include("op_V_loc.jl")
include("op_H.jl")
include("calc_grad_v2.jl")
include("calc_rho.jl")
include("Poisson_solve.jl")
include("LDA_VWN.jl")

include("calc_grad_Haux.jl")
include("calc_Focc.jl")
include("smear_FD.jl")

function test_main( Ns )

    const LatVecs = 16.0*diagm( ones(3) )

    pw = PWGrid( Ns, LatVecs )

    const Ω  = pw.Ω
    const r  = pw.r
    const G  = pw.gvec.G
    const G2 = pw.gvec.G2
    const Npoints = prod(Ns)
    const Ngwx = pw.gvecw.Ngwx

    #
    # Generate array of distances
    #
    center = sum(LatVecs,2)/2
    dr = gen_dr( r, center )
    #
    # Setup potential
    #
    V_ionic = init_pot_harm_3d( pw, dr )
    #
    const Nstates = 5
    const Nelectrons = 4.0
    Focc = zeros(Nstates)
    Focc[1:2] = [2.0, 2.0]

    srand(1234)
    psi = randn(Ngwx,Nstates) + im*randn(Ngwx,Nstates)
    psi = ortho_gram_schmidt(psi)

    Potentials = PotentialsT( V_ionic, zeros(Npoints), zeros(Npoints) )
    rho = calc_rho( pw, Focc, psi )
    Potentials.Hartree = real( G_to_R( Ns, Poisson_solve(pw, rho) ) )
    Potentials.XC = excVWN( rho ) + rho .* excpVWN( rho )

    # Get eigenvalues
    mu = psi' * op_H(pw, Potentials, psi)
    evals = sort(real(eigvals(mu)))
    println(evals)

    kT = 0.01
    Focc, E_fermi = calc_Focc(evals, Nelectrons, kT, is_spinpol=false)

    println("Focc = ", Focc)

    g_Haux = calc_grad_Haux( pw, Potentials, Focc, evals, psi, 0.01 )

    println("sum g_Haux = ", sum(g_Haux))

end

test_main( [30,30,30] )


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using FMIFlux
using Test

@testset "FMIFlux.jl" begin
    if Sys.iswindows()
        @info "Automated testing for Windows is supported."
        @testset "ME-NeuralFMU" begin
            include("hybrid_ME.jl")
        end
        @testset "CS-NeuralFMU" begin
            include("hybrid_CS.jl")
        end
    elseif Sys.islinux()
        @warn "Test-sets are using Windows-FMUs, automated testing for Linux is currently not supported."
    elseif Sys.isapple()
        @warn "Test-sets are using Windows-FMUs, automated testing for macOS is currently not supported."
    end
end


text:
module TestYFilesWiringDiagrams

using Test
using Catlab.WiringDiagrams, Catlab.Graphics

function read_yfiles(name::String; kw...)
  filename = joinpath(@__DIR__, "data", name)
  read_yfiles_diagram(Symbol, Symbol, filename; kw...)
end

ports(n::Int) = repeat([nothing], n)
d = WiringDiagram(ports(1), ports(1))
fv = add_box!(d, Box(:f, ports(1), ports(2)))
gv = add_box!(d, Box(:g, ports(2), ports(1)))
hv = add_box!(d, Box(:h, ports(1), ports(1)))
kv = add_box!(d, Box(:k, ports(1), ports(1)))
add_wires!(d, [
  Wire(:A, (input_id(d),1) => (fv, 1)),
  Wire(:B, (fv,1) => (hv,1)), Wire(:C, (fv,2) => (kv,1)),
  Wire(:D, (hv,1) => (gv,1)), Wire(:E, (kv,1) => (gv,2)),
  Wire(:F, (gv,1) => (output_id(d),1)),
])

d_vertical = read_yfiles("yed_vertical.graphml", orientation=TopToBottom)
@test nboxes(d_vertical) == nboxes(d)
perm = sortperm(boxes(d_vertical); by=box->box.value)
@test is_permuted_equal(d_vertical, d, perm)

d_horizontal = read_yfiles("yed_horizontal.graphml", orientation=LeftToRight)
@test nboxes(d_horizontal) == nboxes(d)
perm = sortperm(boxes(d_horizontal); by=box->box.value)
@test is_permuted_equal(d_horizontal, d, perm)

end


text:
# __author__:aliyousuuf@gmail.com

# This file contains all the codes from lecture, pactise quizzes, and assesment
# of week 02.

# Assume that we've downloaded the EBola data from wiki's website,
# and save it as csv format under file name epidata.csv.

epidata = readdlm("epidata.csv", ',')

col1 = epidata[:,1]
for date = 1:54
    col1[date] = DateTime(col1[date], "d u y")
end

#println(col1)
#println(Dates.datetime2rata(col1[1]))

dayssincemar22(x) = Dates.datetime2rata(x)
epidays = Array(Int64, 54)
for i = 1:54
    epidays[i] = dayssincemar22(col1[i])
end

epidata[:,1] = epidays
writedlm("epidata_converted.csv", epidata, ',')
epidata_converted = readdlm("epidata_converted.csv",',')

# Now we're going to plot epidemic days vs all the cases.

epidays = epidata_converted[:,1]
allcases = epidata_converted[:,2]

# akkcases contains two types of data--string and Int--we've to convert them
# into one type.

for item = 1:54
    if typeof(allcases[item,1]) == SubString{String}
        allcases[item,1] = parse.(Int64,replace(allcases[item,1], ",", ""))
    end
end

using Plots
#Plots.PyPlotBackend()
#pyplot()
#using Gadfly
#plot(epidays, allcases)

# Week_02_Practice Quiz: Q.2

f(x) = 3 * x^2 + 6 * x - 9
for x = -5:5
  #println((,x, ", ", f(x), ")")
end
using Plots
#plot(f, -4, 3) # plot f over [-4,4]
#plot!(zero, -4, 3)

#Week_02_Practice Quiz: Q.3

data = [1.6800483  -1.641695388;
        0.501309281 -0.977697538;
        1.528012113 0.52771122;
        1.70012253 1.711524991;
        1.992493625 1.891000015;
        2.706075824 -0.463427794;
        2.994931927 -0.443566619;
        3.491852811 -1.275179133;
        3.501191722 -0.690499597;
        4.459924502 -5.516130799;
        4.936965851 -6.001703074;
        5.023289852 -8.36416901;
        5.04233698 -7.924477517;
        5.50739285 -10.77482371;
        5.568665171 -10.9171878]

x,y = data[:,1] , data[:,2]
Plots.PyPlotBackend()
#plot(x, y, linetype = :scatter, leg = false)

#Week_02_Practice Quiz: Q.7
n = 20
x = sort(rand(20)); y = rand(20)
#plot!(x,y, leg=false, title = "A sample plot")

# gui()

# replace all the instances of "-" and "--" into 0, and "," into "".
row,col = size(epidata_converted)
for r = 1:row
    for c = 1: col
        if typeof(epidata_converted[r,c]) == SubString{String}
            epidata_converted[r,c] = replace(epidata_converted[r,c], ",", "")
            if isnumber(string(epidata_converted[r,c])) == false
                epidata_converted[r,c] = 0
            end
        end
    end
end

# Convert all the string into integer
for r = 1:row
    for c = 1: col
       if isnumber(string(epidata_converted[r,c])) == true
           epidata_converted[r,c] = parse(Int64,string(epidata_converted[r,c]))
       end
   end
end

max_stars_repo_name: aliyyousuff/MOOC

max_issues_repo_name: aliyyousuff/MOOC

max_forks_repo_name: aliyyousuff/MOOC


text:
"""
    linearized_air_pressure(ρ, ρe_tot, ρe_pot, ρq_tot=0, ρq_liq=0, ρq_ice=0)

The air pressure, linearized around a dry rest state, from the equation of state
(ideal gas law) where:

 - `ρ` (moist-)air density
 - `ρe_tot` total energy density
 - `ρe_pot` potential energy density
and, optionally,
 - `ρq_tot` total water density
 - `ρq_liq` liquid water density
 - `ρq_ice` ice density
"""
function linearized_air_pressure(ρ::FT, ρe_tot::FT, ρe_pot::FT,
                                 ρq_tot::FT=FT(0), ρq_liq::FT=FT(0), ρq_ice::FT=FT(0)) where {FT<:Real}
  ρ*FT(R_d)*FT(T_0) + FT(R_d)/FT(cv_d)*(ρe_tot - ρe_pot - (ρq_tot - ρq_liq)*FT(e_int_v0) + ρq_ice*(FT(e_int_i0) + FT(e_int_v0)))
end

@inline function linearized_pressure(::DryModel, orientation::Orientation, state::Vars, aux::Vars)
  ρe_pot = state.ρ * gravitational_potential(orientation, aux)
  linearized_air_pressure(state.ρ, state.ρe, ρe_pot)
end
@inline function linearized_pressure(::EquilMoist, orientation::Orientation, state::Vars, aux::Vars)
  ρe_pot = state.ρ * gravitational_potential(orientation, aux)
  linearized_air_pressure(state.ρ, state.ρe, ρe_pot, state.moisture.ρq_tot)
end

abstract type AtmosLinearModel <: BalanceLaw
end

vars_state(lm::AtmosLinearModel, FT) = vars_state(lm.atmos,FT)
vars_gradient(lm::AtmosLinearModel, FT) = @vars()
vars_diffusive(lm::AtmosLinearModel, FT) = @vars()
vars_aux(lm::AtmosLinearModel, FT) = vars_aux(lm.atmos,FT)
vars_integrals(lm::AtmosLinearModel, FT) = @vars()

function update_aux!(dg::DGModel, lm::AtmosLinearModel, Q::MPIStateArray, t::Real)
  return false
end
integrate_aux!(lm::AtmosLinearModel, integ::Vars, state::Vars, aux::Vars) = nothing
flux_diffusive!(lm::AtmosLinearModel, flux::Grad, state::Vars, diffusive::Vars, aux::Vars, t::Real) = nothing
function wavespeed(lm::AtmosLinearModel, nM, state::Vars, aux::Vars, t::Real)
  ref = aux.ref_state
  return soundspeed_air(ref.T)
end

function boundary_state!(nf::Rusanov, lm::AtmosLinearModel, x...)
  atmos_boundary_state!(nf, NoFluxBC(), lm.atmos, x...)
end
function boundary_state!(nf::CentralNumericalFluxDiffusive, lm::AtmosLinearModel, x...)
  nothing
end
init_aux!(lm::AtmosLinearModel, aux::Vars, geom::LocalGeometry) = nothing
init_state!(lm::AtmosLinearModel, state::Vars, aux::Vars, coords, t) = nothing

struct AtmosAcousticLinearModel{M} <: AtmosLinearModel
  atmos::M
  function AtmosAcousticLinearModel(atmos::M) where {M}
    if atmos.ref_state === NoReferenceState()
      error("AtmosAcousticLinearModel needs a model with a reference state")
    end
    new{M}(atmos)
  end
end

function flux_nondiffusive!(lm::AtmosAcousticLinearModel, flux::Grad, state::Vars, aux::Vars, t::Real)
  FT = eltype(state)
  ref = aux.ref_state
  e_pot = gravitational_potential(lm.atmos.orientation, aux)

  flux.ρ = state.ρu
  pL = linearized_pressure(lm.atmos.moisture, lm.atmos.orientation, state, aux)
  flux.ρu += pL*I
  flux.ρe = ((ref.ρe + ref.p)/ref.ρ - e_pot)*state.ρu
  nothing
end
function source!(lm::AtmosAcousticLinearModel, source::Vars, state::Vars, aux::Vars, t::Real)
  nothing
end

struct AtmosAcousticGravityLinearModel{M} <: AtmosLinearModel
  atmos::M
  function AtmosAcousticGravityLinearModel(atmos::M) where {M}
    if atmos.ref_state === NoReferenceState()
      error("AtmosAcousticGravityLinearModel needs a model with a reference state")
    end
    new{M}(atmos)
  end
end
function flux_nondiffusive!(lm::AtmosAcousticGravityLinearModel, flux::Grad, state::Vars, aux::Vars, t::Real)
  FT = eltype(state)
  ref = aux.ref_state
  e_pot = gravitational_potential(lm.atmos.orientation, aux)

  flux.ρ = state.ρu
  pL = linearized_pressure(lm.atmos.moisture, lm.atmos.orientation, state, aux)
  flux.ρu += pL*I
  flux.ρe = ((ref.ρe + ref.p)/ref.ρ)*state.ρu
  nothing
end
function source!(lm::AtmosAcousticGravityLinearModel, source::Vars, state::Vars, aux::Vars, t::Real)
  ∇Φ = ∇gravitational_potential(lm.atmos.orientation, aux)
  source.ρu -= state.ρ * ∇Φ
  nothing
end

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
# demo
# ====

# fibonacci
# ---------

fib(n) = n ≤ 2 ? n : fib(n-1) + fib(n-2)

fib(1000)   # never terminates in ordinal execution
fib(m)      # undef var
fib("1000") # obvious type error

# language features
# -----------------

# user-defined types, macros
struct Ty{T}
    fld::T
end

function foo(a)
    v = Ty(a)
    return bar(v)
end

# macros will be expanded
@inline bar(n::T)     where {T<:Number} = n < 0 ? zero(T) : one(T)
@inline bar(v::Ty{T}) where {T<:Number} = bar(v.fdl) # typo "fdl"
@inline bar(v::Ty)                      = bar(convert(Number, v.fld))

foo(1.2)
foo("1") # `String` can't be converted to `Number`

# even staged programming
# adapted from https://github.com/JuliaLang/julia/blob/9f665c19e076ab37cbca2d0cc99283b82e99c26f/base/namedtuple.jl#L253-L264
@generated function badmerge(a::NamedTuple{an}, b::NamedTuple{bn}) where {an, bn}
    names = Base.merge_names(an, bn)
    types = Base.merge_types(names, a, b)
    vals = Any[ :(getfield($(Base.sym_in(names[n], bn) ? :b : :a), $(names[n]))) for n in 1:length(names) ] # missing quote, just ends up with under vars
    :( NamedTuple{$names,$types}(($(vals...),)) )
end

badmerge((x=1,y=2), (y=3,z=1))


text:
using Keldysh, Test

@testset "time_gf" begin

  let c = FullContour(tmax=2.0, β=5.0), nt = 21, ntau=51
    ϵ = -1.0
    f = (t1, t2) -> -1.0im * (heaviside(t1.bpoint, t2.bpoint) - fermi(ϵ, c.β)) * exp(-1.0im * (t1.bpoint.val - t2.bpoint.val) * ϵ)

    grid = FullTimeGrid(c, nt, ntau)

    G1 = GenericTimeGF(f, grid, 1, true)
    G1_ = similar(G1)
    G1_ = zero(G1)

    @test Keldysh.jump(G1) == 1.0im

    G2 = FullTimeGF(f, grid, 1, fermionic, true)
    G2_ = similar(G2)
    G2_ = zero(G2)

    @test all(map(((t1, t2),) -> isapprox(G1[t1,t2], G2[t1,t2], atol=1e-10), Iterators.product(grid, grid)))

    G3 = TimeInvariantFullTimeGF(f, grid, 1, fermionic, true)
    G3_ = similar(G3)
    G3_ = zero(G3)

    @test all(map(((t1, t2),) -> isapprox(G1[t1,t2], G3[t1,t2], atol=1e-10), Iterators.product(grid, grid)))

    # test interpolation on linear function
    lin_f = (t1, t2) -> (2.0 * t1.val - im * t2.val) + 5.0 * heaviside(t1, t2)
    for t1 in grid
      for t2 in grid
        G1[t1, t2] = lin_f(t1.bpoint, t2.bpoint)
      end
    end

    grid_fine = FullTimeGrid(c, 41, 101)

    @test all(map(Iterators.product(grid_fine, grid_fine)) do (t1, t2)
                isapprox(G1(t1.bpoint, t2.bpoint), lin_f(t1.bpoint, t2.bpoint), atol=1e-10)
              end
             )
  end

  let c = KeldyshContour(tmax=2.0), nt = 21

    ϵ = -1.0
    β = 5.0
    f = (t1, t2) -> -1.0im * (heaviside(t1.bpoint, t2.bpoint) - fermi(ϵ, β)) * exp(-1.0im * (t1.bpoint.val - t2.bpoint.val) * ϵ)

    grid = KeldyshTimeGrid(c, nt)

    G1 = GenericTimeGF(f, grid, 1, true)
    G1_ = similar(G1)
    G1_ = zero(G1)

    @test Keldysh.jump(G1) == 1.0im

    G2 = KeldyshTimeGF(f, grid, 1, fermionic, true)
    G2_ = similar(G2)
    G2_ = zero(G2)

    @test all(map(((t1, t2),) -> isapprox(G1[t1,t2], G2[t1,t2], atol=1e-10), Iterators.product(grid, grid)))

    G3 = TimeInvariantKeldyshTimeGF(f, grid, 1, fermionic, true)
    G3_ = similar(G3)
    G3_ = zero(G3)

    @test all(map(((t1, t2),) -> isapprox(G1[t1,t2], G3[t1,t2], atol=1e-10), Iterators.product(grid, grid)))

  end

  let c = ImaginaryContour(β=5.0), ntau=51
    ϵ = -1.0
    f = (t1, t2) -> -1.0im * (heaviside(t1.bpoint, t2.bpoint) - fermi(ϵ, c.β)) * exp(-1.0im * (t1.bpoint.val - t2.bpoint.val) * ϵ)
    grid = ImaginaryTimeGrid(c, ntau)

    G1 = GenericTimeGF(f, grid, 1, true)
    G1_ = similar(G1)
    G1_ = zero(G1)

    G2 = ImaginaryTimeGF(f, grid, 1, fermionic, true)
    G2_ = similar(G2)
    G2_ = zero(G2)
    @test all(map(((t1, t2),) -> isapprox(G1[t1,t2], G2[t1,t2], atol=1e-10), Iterators.product(grid, grid)))

  end

end


text:
function add(obj::SeriesTerm, arg0::jint, arg1::jint, arg2::jdouble, arg3::jdouble)
    return jcall(obj, "add", void, (jint, jint, jdouble, jdouble), arg0, arg1, arg2, arg3)
end

function build_term(::Type{SeriesTerm}, arg0::jint, arg1::jint, arg10::jint, arg11::jint, arg12::jint, arg13::jint, arg14::jint, arg2::jint, arg3::jint, arg4::jint, arg5::jint, arg6::jint, arg7::jint, arg8::jint, arg9::jint)
    return jcall(SeriesTerm, "buildTerm", SeriesTerm, (jint, jint, jint, jint, jint, jint, jint, jint, jint, jint, jint, jint, jint, jint, jint), arg0, arg1, arg10, arg11, arg12, arg13, arg14, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function derivative(obj::SeriesTerm, arg0::BodiesElements)
    return jcall(obj, "derivative", Vector{jdouble}, (BodiesElements,), arg0)
end

function derivative(obj::SeriesTerm, arg0::FieldBodiesElements)
    return jcall(obj, "derivative", Vector{RealFieldElement}, (FieldBodiesElements,), arg0)
end

function get_cos_coeff(obj::SeriesTerm, arg0::jint, arg1::jint)
    return jcall(obj, "getCosCoeff", jdouble, (jint, jint), arg0, arg1)
end

function get_degree(obj::SeriesTerm, arg0::jint)
    return jcall(obj, "getDegree", jint, (jint,), arg0)
end

function get_sin_coeff(obj::SeriesTerm, arg0::jint, arg1::jint)
    return jcall(obj, "getSinCoeff", jdouble, (jint, jint), arg0, arg1)
end

function value(obj::SeriesTerm, arg0::BodiesElements)
    return jcall(obj, "value", Vector{jdouble}, (BodiesElements,), arg0)
end

function value(obj::SeriesTerm, arg0::FieldBodiesElements)
    return jcall(obj, "value", Vector{RealFieldElement}, (FieldBodiesElements,), arg0)
end


text:
using Test
using VkComputeUtils
using Vulkan
import SwiftShader_jll

if !haskey(ENV, "VKCOMPUTEUTILS_RUN_ON_GPU")
    @set_driver :SwiftShader
end

@testset "VkComputeUtils tests" begin
    include("simple_shader.jl")
end


text:
module MedPipe3D

include(joinpath("structs","distinctColorsSaved.jl"))
include(joinpath("fromMonai","LoadFromMonai.jl"))
include(joinpath("visualizationUtils","visualizationFromHdf5.jl"))
include(joinpath("HDF5Utils","HDF5saveUtils.jl"))

end # module