text:
module CryptoSignatures

using CryptoGroups
using Random

rng() = RandomDevice()

# const _default_rng = Ref{RandomDevice}()
# function __init__()
#     _default_rng[] = RandomDevice()
# end

# default_rng() = _default_rng[]

abstract type AbstractSignature end

id(s::AbstractSignature) = s.pubkey ### One can overwrite this as one wishes

"""
Never use this for deployed application!!!
"""
#verify(data,hashnum::UInt64) = hash("$data")==hashnum
verify(data,hash) = errror("Must be implemented by hash type.")
verify(data,s::AbstractSignature,G::AbstractGroup) = verify(s,G) && verify(data,s.hash)

abstract type AbstractSigner end
#id(s::AbstractSigner) = s.pubkey

abstract type AbstractEncryptionKey end
abstract type AbstractDecryptionKey end

# The singer does store group to prevent the private key to be accidentally missused with different group which would lower the security. And the singature would be useless in such case.
struct Signer{T} <: AbstractSigner where T<:Integer
    privkey::T
    pubkey::T
    G::AbstractGroup
end

import Base.Dict
function Dict(signer::Signer)
    dict = Dict()
    dict["priv"] = string(signer.privkey,base=16)
    dict["pub"] = string(signer.pubkey,base=16)
    return dict
end

function Signer{BigInt}(dict::Dict,G::AbstractGroup)
    priv = parse(BigInt,dict["priv"],base=16)
    pub = parse(BigInt,dict["pub"],base=16)
    Signer(priv,pub,G)
end

function Signer(G::AbstractGroup;rng=rng())
    x = rand(1:order(G))
    y = value(G^x)
    Signer(x,y,G)
end

import Base.==
==(x::Signer,y::Signer) = x.privkey==y.privkey && x.pubkey==y.pubkey && x.G==y.G

### There are many different ways one could sing stuff.
# include("rsasignatures.jl")
include("dsasignatures.jl")

export verify, Signer, DSASignature

end # module


text:
abstract type DynamicInverterComponent <: DynamicComponent end
abstract type Converter <: DynamicInverterComponent end
abstract type DCSource <: DynamicInverterComponent end
abstract type Filter <: DynamicInverterComponent end
abstract type FrequencyEstimator <: DynamicInverterComponent end
abstract type InnerControl <: DynamicInverterComponent end

abstract type ActivePowerControl <: DeviceParameter end
abstract type ReactivePowerControl <: DeviceParameter end

ext: jl


text:
function unspliced(τ, u₀, α, β)
    expᵘ = exp(-β * τ)
    u₀ * expᵘ + α / β * (1 - expᵘ)
end

function spliced(τ, s₀, u₀, α, β, γ)
    c = (α - u₀ * β) / (γ - β)
    expᵘ = exp(-β * τ)
    expˢ = exp(-γ * τ)
    s₀ * expˢ + α / γ * (1 - expˢ) + c * (expˢ - expᵘ)
end

function mRNA(τ, u₀, s₀, α, β, γ)
    expᵘ, expˢ = exp(-β * τ), exp(-γ * τ)
    expᵘˢ = (α - u₀ * β) * inv(γ - β) * (expˢ - expᵘ)
    u = u₀ * expᵘ + α / β * (1 - expᵘ)
    s = s₀ * expˢ + α / γ * (1 - expˢ) + expᵘˢ
    u, s
end


text:
using Oceananigans.Advection: AbstractAdvectionScheme
using Oceananigans.Operators: ℑxᶠᵃᵃ, ℑxᶜᵃᵃ, ℑyᵃᶠᵃ, ℑyᵃᶜᵃ, ℑzᵃᵃᶠ, ℑzᵃᵃᶜ
using Oceananigans.TurbulenceClosures: AbstractTurbulenceClosure, AbstractTimeDiscretization

const ATC = AbstractTurbulenceClosure
const ATD = AbstractTimeDiscretization

struct RasterDepthMask end

struct GridFittedBoundary{M, S} <: AbstractImmersedBoundary
    mask :: S
    mask_type :: M
end

GridFittedBoundary(mask; mask_type=nothing) = GridFittedBoundary(mask, mask_type)

@inline is_immersed(i, j, k, underlying_grid, ib::GridFittedBoundary) = ib.mask(node(c, c, c, i, j, k, underlying_grid)...)
@inline is_immersed(i, j, k, underlying_grid, ib::GridFittedBoundary{<:RasterDepthMask}) = ib.mask(i, j)

const IBG = ImmersedBoundaryGrid

const c = Center()
const f = Face()

#####
##### GridFittedImmersedBoundaryGrid
#####

const GFIBG = ImmersedBoundaryGrid{FT, TX, TY, TZ, G, <:GridFittedBoundary} where {FT, TX, TY, TZ, G}

@inline solid_cell(i, j, k, ibg) = is_immersed(i, j, k, ibg.grid, ibg.immersed_boundary)

@inline solid_node(LX, LY, LZ, i, j, k, ibg) = solid_cell(i, j, k, ibg) # fallback (for Center or Nothing LX, LY, LZ)

@inline solid_node(::Face, LX, LY, i, j, k, ibg) = solid_cell(i, j, k, ibg) || solid_cell(i-1, j, k, ibg)
@inline solid_node(LX, ::Face, LZ, i, j, k, ibg) = solid_cell(i, j, k, ibg) || solid_cell(i, j-1, k, ibg)
@inline solid_node(LX, LY, ::Face, i, j, k, ibg) = solid_cell(i, j, k, ibg) || solid_cell(i, j, k-1, ibg)

@inline solid_node(::Face, ::Face, LZ, i, j, k, ibg) = solid_node(c, f, c, i, j, k, ibg) || solid_node(c, f, c, i-1, j, k, ibg)
@inline solid_node(::Face, LY, ::Face, i, j, k, ibg) = solid_node(c, c, f, i, j, k, ibg) || solid_node(c, c, f, i-1, j, k, ibg)
@inline solid_node(LX, ::Face, ::Face, i, j, k, ibg) = solid_node(c, f, c, i, j, k, ibg) || solid_node(c, f, c, i, j, k-1, ibg)

@inline solid_node(::Face, ::Face, ::Face, i, j, k, ibg) = solid_node(c, f, f, i, j, k, ibg) || solid_node(c, f, f, i-1, j, k, ibg)

@inline conditional_flux_ccc(i, j, k, ibg::IBG{FT}, flux, args...) where FT = ifelse(solid_node(c, c, c, i, j, k, ibg), zero(FT), flux(i, j, k, ibg, args...))
@inline conditional_flux_ffc(i, j, k, ibg::IBG{FT}, flux, args...) where FT = ifelse(solid_node(f, f, c, i, j, k, ibg), zero(FT), flux(i, j, k, ibg, args...))
@inline conditional_flux_fcf(i, j, k, ibg::IBG{FT}, flux, args...) where FT = ifelse(solid_node(f, c, f, i, j, k, ibg), zero(FT), flux(i, j, k, ibg, args...))
@inline conditional_flux_cff(i, j, k, ibg::IBG{FT}, flux, args...) where FT = ifelse(solid_node(c, f, f, i, j, k, ibg), zero(FT), flux(i, j, k, ibg, args...))

@inline conditional_flux_fcc(i, j, k, ibg::IBG{FT}, flux, args...) where FT = ifelse(solid_node(f, c, c, i, j, k, ibg), zero(FT), flux(i, j, k, ibg, args...))
@inline conditional_flux_cfc(i, j, k, ibg::IBG{FT}, flux, args...) where FT = ifelse(solid_node(c, f, c, i, j, k, ibg), zero(FT), flux(i, j, k, ibg, args...))
@inline conditional_flux_ccf(i, j, k, ibg::IBG{FT}, flux, args...) where FT = ifelse(solid_node(c, c, f, i, j, k, ibg), zero(FT), flux(i, j, k, ibg, args...))

#####
##### Advective fluxes
#####

# ccc, ffc, fcf
 @inline _viscous_flux_ux(i, j, k, ibg::GFIBG, args...) = conditional_flux_ccc(i, j, k, ibg, viscous_flux_ux, args...)
 @inline _viscous_flux_uy(i, j, k, ibg::GFIBG, args...) = conditional_flux_ffc(i, j, k, ibg, viscous_flux_uy, args...)
 @inline _viscous_flux_uz(i, j, k, ibg::GFIBG, args...) = conditional_flux_fcf(i, j, k, ibg, viscous_flux_uz, args...)

 # ffc, ccc, cff
 @inline _viscous_flux_vx(i, j, k, ibg::GFIBG, args...) = conditional_flux_ffc(i, j, k, ibg, viscous_flux_vx, args...)
 @inline _viscous_flux_vy(i, j, k, ibg::GFIBG, args...) = conditional_flux_ccc(i, j, k, ibg, viscous_flux_vy, args...)
 @inline _viscous_flux_vz(i, j, k, ibg::GFIBG, args...) = conditional_flux_cff(i, j, k, ibg, viscous_flux_vz, args...)

 # fcf, cff, ccc
 @inline _viscous_flux_wx(i, j, k, ibg::GFIBG, args...) = conditional_flux_fcf(i, j, k, ibg, viscous_flux_wx, args...)
 @inline _viscous_flux_wy(i, j, k, ibg::GFIBG, args...) = conditional_flux_cff(i, j, k, ibg, viscous_flux_wy, args...)
 @inline _viscous_flux_wz(i, j, k, ibg::GFIBG, args...) = conditional_flux_ccc(i, j, k, ibg, viscous_flux_wz, args...)

# fcc, cfc, ccf
@inline _diffusive_flux_x(i, j, k, ibg::GFIBG, args...) = conditional_flux_fcc(i, j, k, ibg, diffusive_flux_x, args...)
@inline _diffusive_flux_y(i, j, k, ibg::GFIBG, args...) = conditional_flux_cfc(i, j, k, ibg, diffusive_flux_y, args...)
@inline _diffusive_flux_z(i, j, k, ibg::GFIBG, args...) = conditional_flux_ccf(i, j, k, ibg, diffusive_flux_z, args...)

#####
##### Advective fluxes
#####

# dx(uu), dy(vu), dz(wu)
# ccc,    ffc,    fcf
@inline _advective_momentum_flux_Uu(i, j, k, ibg::GFIBG, args...) = conditional_flux_ccc(i, j, k, ibg, advective_momentum_flux_Uu, args...)
@inline _advective_momentum_flux_Vu(i, j, k, ibg::GFIBG, args...) = conditional_flux_ffc(i, j, k, ibg, advective_momentum_flux_Vu, args...)
@inline _advective_momentum_flux_Wu(i, j, k, ibg::GFIBG, args...) = conditional_flux_fcf(i, j, k, ibg, advective_momentum_flux_Wu, args...)

# dx(uv), dy(vv), dz(wv)
# ffc,    ccc,    cff
@inline _advective_momentum_flux_Uv(i, j, k, ibg::GFIBG, args...) = conditional_flux_ffc(i, j, k, ibg, advective_momentum_flux_Uv, args...)
@inline _advective_momentum_flux_Vv(i, j, k, ibg::GFIBG, args...) = conditional_flux_ccc(i, j, k, ibg, advective_momentum_flux_Vv, args...)
@inline _advective_momentum_flux_Wv(i, j, k, ibg::GFIBG, args...) = conditional_flux_cff(i, j, k, ibg, advective_momentum_flux_Wv, args...)

# dx(uw), dy(vw), dz(ww)
# fcf,    cff,    ccc
@inline _advective_momentum_flux_Uw(i, j, k, ibg::GFIBG, args...) = conditional_flux_fcf(i, j, k, ibg, advective_momentum_flux_Uw, args...)
@inline _advective_momentum_flux_Vw(i, j, k, ibg::GFIBG, args...) = conditional_flux_cff(i, j, k, ibg, advective_momentum_flux_Vw, args...)
@inline _advective_momentum_flux_Ww(i, j, k, ibg::GFIBG, args...) = conditional_flux_ccc(i, j, k, ibg, advective_momentum_flux_Ww, args...)

   @inline _advective_tracer_flux_x(i, j, k, ibg::GFIBG, args...) = conditional_flux_fcc(i, j, k, ibg, advective_tracer_flux_x, args...)
   @inline _advective_tracer_flux_y(i, j, k, ibg::GFIBG, args...) = conditional_flux_cfc(i, j, k, ibg, advective_tracer_flux_y, args...)
   @inline _advective_tracer_flux_z(i, j, k, ibg::GFIBG, args...) = conditional_flux_ccf(i, j, k, ibg, advective_tracer_flux_z, args...)

#####
##### "Boundary-aware" interpolation
#####
##### Don't interpolate dead cells.
#####

@inline near_x_boundary(i, j, k, ibg, ::AbstractAdvectionScheme{0}) = false
@inline near_y_boundary(i, j, k, ibg, ::AbstractAdvectionScheme{0}) = false
@inline near_z_boundary(i, j, k, ibg, ::AbstractAdvectionScheme{0}) = false

@inline near_x_boundary(i, j, k, ibg, ::AbstractAdvectionScheme{1}) = solid_cell(i - 1, j, k, ibg) || solid_cell(i, j, k, ibg) || solid_cell(i + 1, j, k, ibg)
@inline near_y_boundary(i, j, k, ibg, ::AbstractAdvectionScheme{1}) = solid_cell(i, j - 1, k, ibg) || solid_cell(i, j, k, ibg) || solid_cell(i, j + 1, k, ibg)
@inline near_z_boundary(i, j, k, ibg, ::AbstractAdvectionScheme{1}) = solid_cell(i, j, k - 1, ibg) || solid_cell(i, j, k, ibg) || solid_cell(i, j, k + 1, ibg)

@inline near_x_boundary(i, j, k, ibg, ::AbstractAdvectionScheme{2}) = solid_cell(i - 2, j, k, ibg) || solid_cell(i - 1, j, k, ibg) || solid_cell(i, j, k, ibg) || solid_cell(i + 1, j, k, ibg) || solid_cell(i + 2, j, k, ibg)
@inline near_y_boundary(i, j, k, ibg, ::AbstractAdvectionScheme{2}) = solid_cell(i, j - 2, k, ibg) || solid_cell(i, j - 1, k, ibg) || solid_cell(i, j, k, ibg) || solid_cell(i, j + 1, k, ibg) || solid_cell(i, j + 2, k, ibg)
@inline near_z_boundary(i, j, k, ibg, ::AbstractAdvectionScheme{2}) = solid_cell(i, j, k - 2, ibg) || solid_cell(i, j, k - 1, ibg) || solid_cell(i, j, k, ibg) || solid_cell(i, j, k + 1, ibg) || solid_cell(i, j, k + 2, ibg)

# Takes forever to compile, but works.
# @inline near_x_boundary(i, j, k, ibg, ::AbstractAdvectionScheme{buffer}) where buffer = any(ntuple(δ -> solid_cell(i - buffer - 1 + δ, j, k, ibg), Val(2buffer + 1)))
# @inline near_y_boundary(i, j, k, ibg, ::AbstractAdvectionScheme{buffer}) where buffer = any(ntuple(δ -> solid_cell(i, j - buffer - 1 + δ, k, ibg), Val(2buffer + 1)))
# @inline near_z_boundary(i, j, k, ibg, ::AbstractAdvectionScheme{buffer}) where buffer = any(ntuple(δ -> solid_cell(i, j, k - buffer - 1 + δ, ibg), Val(2buffer + 1)))

for bias in (:symmetric, :left_biased, :right_biased)
    for (d, ξ) in enumerate((:x, :y, :z))

        code = [:ᵃ, :ᵃ, :ᵃ]

        for loc in (:ᶜ, :ᶠ)
            code[d] = loc
            second_order_interp = Symbol(:ℑ, ξ, code...)
            interp = Symbol(bias, :_interpolate_, ξ, code...)
            alt_interp = Symbol(:_, interp)

            near_boundary = Symbol(:near_, ξ, :_boundary)

            # Conditional high-order interpolation in Bounded directions
            @eval begin
                import Oceananigans.Advection: $alt_interp
                using Oceananigans.Advection: $interp

                @inline $alt_interp(i, j, k, ibg::ImmersedBoundaryGrid, scheme, ψ) =
                    ifelse($near_boundary(i, j, k, ibg, scheme),
                           $second_order_interp(i, j, k, ibg.grid, ψ),
                           $interp(i, j, k, ibg.grid, scheme, ψ))

                # @inline $alt_interp(i, j, k, ibg::IBG, scheme, ψ) = $interp(i, j, k, ibg.grid, scheme, ψ)
            end
        end
    end
end

#####
##### Masking for GridFittedBoundary
#####

@inline function scalar_mask(i, j, k, grid, ::GridFittedBoundary, LX, LY, LZ, value, field)
    return @inbounds ifelse(solid_node(LX, LY, LZ, i, j, k, grid),
                            value,
                            field[i, j, k])
end

mask_immersed_velocities!(U, arch, grid::GFIBG) = Tuple(mask_immersed_field!(q) for q in U)


text:
# Run this cell to load the graphics packages

using Plots; gr()
using Interact

# ------------------------------------------------------------------------------------------
# ## Adding a function parameter
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# In the last notebook, we saw an example of adding **parameters** to functions. These are
# values that control the behavior of a function. Let's look at that in some more detail.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Let's go back to our original version of the σ function:
# ------------------------------------------------------------------------------------------

σ(x) = 1 / (1 + exp(-x))

# ------------------------------------------------------------------------------------------
# Instead of working with a single function, we can work with a whole class (set) of
# functions that look similar but differ in the value of a **parameter**. Let's make a new
# function that uses the previous $\sigma$ function, but also has a parameter, $w$.
# Mathematically, we could write
#
# $$f_w(x) = f(x; w) = \sigma(w \, x).$$
#
# (Here, $w$ and $x$ are multiplied in the argument of $\sigma$; we could write $w \times x$
# or $w * x$, or even $w \cdot x$, but usually the symbols are not written.)
#
# Mathematically speaking, we can think of $f_w$ as a different function for each different
# value of the parameter $w$.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# In Julia, we write this as follows:
# ------------------------------------------------------------------------------------------

f(x, w) = σ(w * x)

# ------------------------------------------------------------------------------------------
# Note that Julia just treats parameters as additional *arguments* to the function; the
# function `f` has two arguments, the value of `x` and the value of `w` that we want to use.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We can now investigate the effect of $w$ interactively. To do so, we need a way of writing
# in Julia "the function of one variable $x$ that we obtain when we fix the value of $w$".
# We write this as an "anonymous function", as we saw in the notebook on functions:
#
#     x -> f(x, w)
#
# We can read this as "the function that maps $x$ to the value of $f(x, w)$, for a value of
# $w$ that was previously given".
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Now we are ready to draw the function. For each plot, we *fix* a value of the parameter
# $w$ and draw the resulting function as a function of $x$. However, `Interact.jl` then
# allows us to modify interactively the value of $w$, and plot the new function that comes
# out:
# ------------------------------------------------------------------------------------------

@manipulate for w in -2:0.01:2

    plot(x->f(x, w), -5, 5, ylims=(0,1), label="sigmoid")
    plot!(x->(x>0), -5,5, label="Square Wave")
end

# ------------------------------------------------------------------------------------------
# #### Exercise 1
#
# Try writing your own function that takes a parameter. Start by copying and executing
#
# ```julia
# square(x) = x^2
# ```
#
# Then use `square` to declare a new function `square_and_scale` that takes two inputs, `a`
# and `x` such that
#
# $$\mathrm{square\_and\_scale}(x; a) := a \cdot x^2$$
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 2
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Once you have declared `square_and_scale`, uncomment the code below and see how the
# parameter `a` scales the function `square` :
# ------------------------------------------------------------------------------------------

# x = -10:10
# @manipulate for a in 0:0.01:10
#     plot(x, square.(x), label="x^2")
#     plot!(x, square_and_scale.(x, a), ls=:dash, label="ax^2")
# end

# ------------------------------------------------------------------------------------------
# ## Fitting a function to data
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# As we saw in the previous notebook, what we would like to do is use the fact that we now
# have a parameter in our function in order to do something useful! Namely, we want to model
# data with it.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# So suppose that we are given a single data point $(x_0, y_0) = (2, 0.8)$. We can try to
# "fit" the function $f_w$ by adjusting the parameter $w$ until the function passes through
# the data.
#
# **Game**: Move the slider until the graph of the function hits the data point. Which value
# of $w$ does that correspond to?
# ------------------------------------------------------------------------------------------

x0, y0 = 2, 0.8

@manipulate for w in -2:0.01:2
    plot(x->f(x, w), -5, 5, ylims=(0, 1), label="f")
    scatter!([x0], [y0], label="data")
end

# ------------------------------------------------------------------------------------------
# ## Quantifying how far we are from the goal: the *loss function*
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We can see visually when the graph of the function passes through the data point. But the
# goal is to be able to automate this process so that the computer can do it unaided.
#
# So we will need a more precise way of deciding and quantifying (i.e. measuring with a
# number) *how far away we are from the goal*; here, the goal means hitting the data point
# with the function.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 3
#
# Can you think of a way of measuring how far away our function is from the data point?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Defining the loss function
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We need to measure how far away the curve is from the data point when we choose a
# particular value of $w$.
# One way to do this is by finding the vertical distance $d$ from the curve to the data
# point.
#
# Instead of just taking the distance, it is common to take the *square* of the distance,
# $d^2$.
#
# Since we are taking the vertical distance, we need the distance at the given value of
# $x_0$ where the data point lies. For a given value of the parameter $w$, the height of the
# point on the curve with that value of $x_0$ is $f(x_0, w)$.
#
# So we take
# $$d := y_0 - f(x_0, w)$$
#
# and
# $$d^2 = [y_0 - f(x_0, w)]^2.$$
#
# This is our measure of distance. It will change when $w$ changes -- in other words, it is
# itself a *function of $w$*; we will denote this function by $L(w)$, and call it the **loss
# function**:
#
# $$L(w) := [y_0 - f(x_0, w)]^2.$$
#
# So the goal is to find the value $w^*$ of $w$ where the loss function is *least*; in other
# words, we need to *minimize* the loss function!
#
# (Another name for a loss function is a *cost function*.)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 4
#
# (a) Define the loss function `L(w)` in Julia.
#
# (b) Draw the data point and the function `x -> f(x, w)`. Also draw a vertical line from
# the data point to the function `x -> f(x, w)`.
#
# (c) Make the plot interactive.
#
# (d) Add as the plot title the value of the loss function for the current value of $w$.
#
# (e) Use the slider to find the value $w^*$ of $w$ for which the loss function reaches its
# minimum value. What is $w^*$? What is the value of the loss function there, $L(w^*)$?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## What does the loss function look like?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# The loss function $L(w)$ tells us how far away the function $f_w$ is from the data when
# the parameter value is $w$, represented visually as the vertical line in the previous
# plot. When the data are fixed, this is a function only of the parameter $w$. What does
# this function look like as a function of $w$? Let's draw it!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 5
#
# Draw the function $L(w)$ as a function of $w$.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Features of the loss function
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# This graph quantifies how far we are from the data point for a given value of $w$.
# What features can we see from the graph?
#
# Firstly, we see that $L(w)$ is always bigger than $0$, for any value of $w$. This is
# because we want $L$ to be some kind of measure of *distance*, and distances cannot be
# negative.
#
# Secondly, we see that there is a special value $w^*$ of $w$ where the function $L$ reaches
# its minimum value. In this particular case, it actually reaches all the way down to $0$!
# This means that the original function $f$ (the one we manipulated above) passes exactly
# through the data point $(x_0, y_0)$.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 6
#
# Draw a zoomed-in version of the graph to find the place $w^*$ where the function hits $0$
# more precisely.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### A different way of defining the loss function
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# **Why did we use such a complicated function $L$ with those squares inside?** We could
# instead just have used the absolute distance, instead of the distance squared, using the
# *absolute value* function, written mathematically as $| \cdot |$, and in Julia as `abs`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 7
#
# Define a new loss function, `L_abs`, using the absolute value, and see what it looks like.
# ------------------------------------------------------------------------------------------

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/intro-to-ml/06. Tools - Adding a function parameter.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
#=
P33 (*) Determine whether two positive integer numbers are coprime.
Two numbers are coprime if their greatest common divisor equals 1.
Example:
* (coprime 35 64)
T
=#

# From P32 but more terse
function gcd(a, b)
    b == 0 ? a : gcd(b, mod(a, b))
end

function coprime(a, b)
    gcd(a, b) == 1
end

@assert coprime(35, 64)
@assert coprime(4, 9)

println("Tests passed: JL-33.jl")


text:
#################################################################################
#Type Alias for long type signatures
const MinMax = NamedTuple{(:min, :max), NTuple{2, Float64}}
const NamedMinMax = Tuple{String, MinMax}
const UpDown = NamedTuple{(:up, :down), NTuple{2, Float64}}
const InOut = NamedTuple{(:in, :out), NTuple{2, Float64}}
const StartUpStages = NamedTuple{(:hot, :warm, :cold), NTuple{3, Float64}}

# Type Alias From other Packages
const PM = PowerModels
const PSY = PowerSystems
const PSI = PowerSimulations
const IS = InfrastructureSystems
const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities
const PJ = ParameterJuMP
const MOPFM = MOI.FileFormats.Model
const TS = TimeSeries

const BUILD_SIMULATION_TIMER = TimerOutputs.TimerOutput()
const RUN_SIMULATION_TIMER = TimerOutputs.TimerOutput()

#Type Alias for JuMP and PJ containers
const JuMPExpressionMatrix = Matrix{<:JuMP.AbstractJuMPScalar}
const PGAE = PJ.ParametrizedGenericAffExpr{Float64, JuMP.VariableRef}
const GAE = JuMP.GenericAffExpr{Float64, JuMP.VariableRef}
const JuMPAffineExpressionArray = Matrix{GAE}
const JuMPAffineExpressionVector = Vector{GAE}
const JuMPConstraintArray = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}
const JuMPVariableArray = JuMP.Containers.DenseAxisArray{JuMP.VariableRef}
const JuMPParamArray = JuMP.Containers.DenseAxisArray{PJ.ParameterRef}
const DenseAxisArrayContainer = Dict{Symbol, JuMP.Containers.DenseAxisArray}

@enum BUILD_STATUS begin
    BUILT = 1
    IN_PROGRESS = -1
    EMPTY = 0
end

@enum SOS_STATUS_VARIABLE begin
    NO_VARIABLE = 1
    PARAMETER = 2
    VARIABLE = 3
end

@enum STAGE_STATUS begin
    SUCESSFUL_RUN = 0
    FAILED_RUN = 1
end

# Settings constants
const UNSET_HORIZON = 0
const UNSET_INI_TIME = Dates.DateTime(0)

# Tolerance of comparisons
# MIP gap tolerances in most solvers are set to 1e-4
const ABSOLUTE_TOLERANCE = 1.0e-3
const BALANCE_SLACK_COST = 1e6
const SERVICES_SLACK_COST = 1e5
const COST_EPSILON = 1e-3
const MISSING_INITIAL_CONDITIONS_TIME_COUNT = 999.0
const SECONDS_IN_MINUTE = 60.0
const MINUTES_IN_HOUR = 60.0
const SECONDS_IN_HOUR = 3600.0
const MAX_START_STAGES = 3
const OBJECTIVE_FUNCTION_POSITIVE = 1.0
const OBJECTIVE_FUNCTION_NEGATIVE = -1.0
# The DEFAULT_RESERVE_COST value is used to avoid degeneracy of the solutions, reseve cost isn't provided.
const DEFAULT_RESERVE_COST = 1.0e-4

# Interface limitations
const OPERATIONS_ACCEPTED_KWARGS = [
    :horizon,
    :initial_time,
    :use_forecast_data,
    :PTDF,
    :use_parameters,
    :optimizer,
    :warm_start,
    :balance_slack_variables,
    :services_slack_variables,
    :system_to_file,
    :constraint_duals,
    :export_pwl_vars,
]

const OPERATIONS_SOLVE_KWARGS = [:optimizer, :save_path]

const STAGE_ACCEPTED_KWARGS = [
    :PTDF,
    :warm_start,
    :balance_slack_variables,
    :services_slack_variables,
    :constraint_duals,
    :system_to_file,
    :export_pwl_vars,
    :allow_fails,
]

const UNSUPPORTED_POWERMODELS =
    [PM.SOCBFPowerModel, PM.SOCBFConicPowerModel, PM.IVRPowerModel]

const PSI_NAME_DELIMITER = "__"

const M_VALUE = 1e6


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

using Markdown
using InteractiveUtils

# ╔═╡ 2444e25e-76e3-11eb-021a-e9e976b267bb
begin
	using PlutoUI
	using NativeSVG
end

# ╔═╡ d77c9e58-76e2-11eb-27a0-23914499f29f
md"""# Strings

Strings are not like integers, floats, and booleans. A string is a *sequence*, which means it is an ordered collection of other values. In this chapter you’ll see how to access the characters that make up a string, and you’ll learn about some of the string helper functions provided by Julia.
"""

# ╔═╡ 3eeafd32-76e3-11eb-01dc-53af6f26eaf7
md"""## Characters

English language speakers are familiar with characters such as the letters of the alphabet (A, B, C, ...), numerals, and common punctuation. These characters are standar‐ dized and mapped to integer values between 0 and 127 by the *ASCII standard* (American Standard Code for Information Interchange).

There are, of course, many other characters used in non-English languages, including variants of the ASCII characters with accents and other modifications, related scripts such as Cyrillic and Greek, and scripts completely unrelated to ASCII and English, including Arabic, Chinese, Hebrew, Hindi, Japanese, and Korean.

The *Unicode standard* tackles the complexities of what exactly a character is, and is generally accepted as the definitive standard addressing this problem. It provides a unique number for every character on a worldwide scale.

A `Char` value represents a single character and is surrounded by single quotes:

```julia
julia> 'x'
'x': ASCII/Unicode U+0078 (category Ll: Letter, lowercase)
julia> typeof('x')
Char
julia> '🍌'
'🍌': Unicode U+1F34C (category So: Symbol, other)
```

Even emojis are part of the Unicode standard (**`\:banana: TAB`**).
"""

# ╔═╡ b29b3ca6-76e3-11eb-0970-1bee030205f2
md"""## A String Is a Sequence

A string is a sequence of characters. You can access the characters one at a time with the bracket operator (`[]`):

```julia
julia> fruit = "banana"
"banana"
julia> letter = fruit[1]
'b': ASCII/Unicode U+0062 (category Ll: Letter, lowercase)
```

The second statement selects character number 1 from `fruit` and assigns it to `letter`.

The expression in brackets is called an *index*. The index indicates which character in the sequence you want (hence the name).

All indexing in Julia is 1-based—the first element of any integer-indexed object is found at index `1` or `begin` and the last element at index `end`:

```julia
julia> fruit[begin]
'b': ASCII/Unicode U+0062 (category Ll: Letter, lowercase)
julia> fruit[end]
'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)
```

As an index, you can use an expression that contains variables and operators:

```julia
julia> i = 1
1
julia> fruit[i+1]
'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)
julia> fruit[end-1]
'n': ASCII/Unicode U+006e (category Ll: Letter, lowercase)
```

But the value of the index has to be an integer. Otherwise, you get:

```julia
julia> letter = fruit[1.5]
ERROR: MethodError: no method matching getindex(::String, ::Float64)
```
"""

# ╔═╡ ade32e66-76e4-11eb-165e-13867b2cdc2d
md"""## `length`

`length` is a built-in function that returns the number of characters in a string:

```julia
julia> fruits = "🍌 🍎 🍐"
"🍌 🍎 🍐"
julia> len = length(fruits)
5
```

To get the last letter of a string, you might be tempted to try something like this:

```julia
julia> last = fruits[len]
' ': ASCII/Unicode U+0020 (category Zs: Separator, space)
```

But you might not get what you expect.

Strings are encoded using *UTF-8 encoding*. UTF-8 is a variable-width encoding, meaning that not all characters are encoded in the same number of bytes.

The function `sizeof` gives the number of bytes in a string:

```julia
julia> sizeof("🍌")
4
```

Because an emoji is encoded in 4 bytes and string indexing is byte-based, the fifth element of fruits is a `SPACE`.

This also means that not every byte index into a UTF-8 string is necessarily a valid index for a character. If you index into a string at an invalid byte index, an error is thrown:

```julia
julia> fruits[2]
ERROR: StringIndexError("🍌 🍎 🍐", 2)
```

In the case of fruits, the character `🍌` is a 4-byte character, so the indices 2, 3, and 4 are invalid and the next character’s index is 5; this next valid index can be computed by `nextind(fruits, 1)`, the next index after that by `nextind(fruits, 5)` and so on.
"""

# ╔═╡ ac377fee-76e5-11eb-1736-9190b5c92214
md"""## Traversal

A lot of computations involve processing a string one character at a time. Often they start at the beginning, select each character in turn, do something to it, and continue until the end. This pattern of processing is called a *traversal*. One way to write a traversal is with a while loop:

```julia
index = firstindex(fruits)
while index <= sizeof(fruits)
	letter = fruits[index] println(letter)
	index = nextind(fruits, index)
end
```

This loop traverses the string and displays each letter on a line by itself. The loop con‐ dition is `index <= sizeof(fruit)`, so when `index` is larger than the number of bytes in the string, the condition is `false` and the body of the loop doesn’t run.

The function `firstindex` returns the first valid byte index.
"""

# ╔═╡ 20c034f2-76e6-11eb-1f6d-891ef3f93b29
md"""#### Exercise 8-1

Write a function that takes a string as an argument and displays the letters backward, one per line.
"""

# ╔═╡ 2c1d4452-76e6-11eb-1cc4-0be827cdb64d
md"""Another way to write a traversal is with a `for` loop:

```julia
for letter in fruits
	println(letter)
end
```

Each time through the loop, the next character in the string is assigned to the variable `letter`. The loop continues until no characters are left.

The following example shows how to use concatenation (string multiplication) and a for loop to generate an abecedarian series (i.e., in alphabetical order). In Robert McCloskey’s book *Make Way for Ducklings* (Puffin), the names of the ducklings are Jack, Kack, Lack, Mack, Nack, Ouack, Pack, and Quack. This loop outputs these names in order:

```julia
prefixes = "JKLMNOPQ"
suffix = "ack"
for letter in prefixes
	println(letter * suffix)
end
```

Although the output isn’t quite right, because “Ouack” and “Quack” are misspelled:

```
Jack
Kack
Lack
Mack
Nack
Oack
Pack
Qack
```
"""

# ╔═╡ 84220d54-76e6-11eb-24d3-b7085b8766b8
md"""#### Exercise 8-2

Modify the program to fix this error.
"""

# ╔═╡ c170dc56-76e6-11eb-2bd1-050911e1915f
md"""## String Slices

A segment of a string is called a *slice*. Selecting a slice is similar to selecting a character:

```julia
julia> str = "Julius Caesar";

julia> str[1:6]
"Julius"
```

!!! tip
    A semicolon in REPL mode not only allows you to put multiple statements on one line but also hides the output.

The operator `[n:m]` returns the part of the string from the `n`th byte to the `m`th byte, so the same caution is needed as for simple indexing.

The `begin` / `end` keyword can be used to indicate the first / last byte of the string:

```julia
julia> str[begin:6]
"Julius"
julia> str[8:end]
"Caesar"
```

If the first index is greater than the second the result is an empty string, represented by two quotation marks:

```julia
julia> str[8:7]
""
```

An empty string contains no characters and has length 0, but other than that, it is the same as any other string.
"""

# ╔═╡ 559c8a42-76e7-11eb-371c-391c5deba644
md"""#### Exercise 8-3

Continuing this example, what do you think `str[:]` means? Try it and see."""

# ╔═╡ 67110424-76e7-11eb-14c3-e195649e6d41
md"""## Strings Are Immutable

It is tempting to use the `[]` operator on the left side of an assignment, with the intention of changing a character in a string. For example:

```julia
julia> greeting = "Hello, world!"
"Hello, world!"
julia> greeting[1] = 'J'
ERROR: MethodError: no method matching setindex!(::String, ::Char, ::Int64)
```

The reason for the error is that strings are *immutable*, which means you can’t change an existing string. The best you can do is create a new string that is a variation on the original:

```julia
julia> greeting = "J" * greeting[2:end]
"Jello, world!"
```

This example concatenates a new first letter onto a slice of `greeting`. It has no effect on the original string.
"""

# ╔═╡ e6a4882e-79fc-11eb-3957-95357985303f
md"""## String Interpolation

Constructing strings using concatenation can become a bit cumbersome. To reduce the need for these verbose calls to string or repeated multiplications, Julia allows *string interpolation* using `$`:

```julia
julia> greet = "Hello"
"Hello"
julia> whom = "World"
"World"
julia> "$greet, $(whom)!"
"Hello, World!"
```

This is more readable and convenient than string concatenation:

```julia
greet * ", " * whom * "!"
```

The shortest complete expression after the `$` is taken as the expression whose value is to be interpolated into the string. Thus, you can interpolate any expression into a string using parentheses:

```julia
julia> "1 + 2 = $(1 + 2)"
"1 + 2 = 3"
```
"""

# ╔═╡ 36ac2624-79fd-11eb-1099-910ca7cea1c9
md"""## Searching
What does the following function do?

```julia
function find(word, letter)
	index = firstindex(word)
	while index <= sizeof(word)
		if word[index] == letter
			return index
		end
		index = nextind(word, index)
	end
	-1
end
```

In a sense, `find` is the inverse of the `[]` operator. Instead of taking an index and extracting the corresponding character, it takes a character and finds the index where that character appears. If the character is not found, the function returns `-1`.

This is the first example we have seen of a `return` statement inside a loop. If `word[index] == letter`, the function breaks out of the loop and returns immediately.

If the character doesn’t appear in the string, the program exits the loop normally and returns `-1`.

This pattern of computation—traversing a sequence and returning when we find what we are looking for—is called a *search*.
"""

# ╔═╡ e04750c8-79fd-11eb-0249-ed3fb7e031b2
md"""#### Exercise 8-4

Modify `find` so that it has a third parameter, the index in word where it should start looking.
"""

# ╔═╡ 019a3968-79fe-11eb-206a-3d92d9352524
md"""## Looping and Counting

The following program counts the number of times the letter a appears in a string:

```julia
word = "banana"
counter = 0
for letter in word
	if letter == 'a'
		counter = counter + 1
	end
end
println(counter)
```

This program demonstrates another pattern of computation called a *counter*. The variable `counter` is initialized to `0` and then incremented each time an `a` is found. When the loop exits, counter contains the result—the total number of `a`’s.
"""

# ╔═╡ 5af18f32-79fe-11eb-1860-0f08e3d25a3f
md"""#### Exercise 8-5

Encapsulate this code in a function named `count`, and generalize it so that it accepts the string and the letter as arguments.

Then rewrite the function so that instead of traversing the string, it uses the three parameter version of `find` from the previous section.
"""

# ╔═╡ 8d5f176e-79fe-11eb-29fc-813db53b35a6
md"""## String Library

Julia provides functions that perform a variety of useful operations on strings. For example, the function `uppercase` takes a string and returns a new string with all uppercase letters:

```julia
julia> uppercase("Hello, World!")
"HELLO, WORLD!"
```

As it turns out, there is a function named `findfirst` that is remarkably similar to the `find` function we wrote:

```julia
julia> findfirst("a", "banana")
2:2
```

Actually, the `findfirst` function is more general than our function; it can find sub‐ strings, not just characters:

```julia
julia> findfirst("na", "banana")
3:4
```

By default, `findfirst` starts at the beginning of the string, but the function `findnext` takes a third argument, the index where it should start looking:

```julia
julia> findnext("na", "banana", 4)
5:6
```
"""

# ╔═╡ e80e049a-79fe-11eb-2ade-83a99e44e5a5
md"""## The `∈` Operator

The operator `∈` (**`\in TAB`**) is a Boolean operator that takes a character and a string and returns `true` if the first appears in the second:

```julia
julia> 'a' ∈ "banana" # 'a' in "banana"
true
```

For example, the following function prints all the letters from `word1` that also appear in `word2`:

```julia
function inboth(word1, word2)
	for letter in word1
		if letter ∈ word2
			print(letter, " ")
		end
	end
end
```

With well-chosen variable names, Julia sometimes reads like English. You could read this loop as “for (each) letter in (the first) word, if (the) letter is an element of (the second) word, print (the) letter.”

Here’s what you get if you compare `"apples"` and `"oranges"`:

```julia
julia> inboth("apples", "oranges")
a e s
```
"""

# ╔═╡ 5487357e-79ff-11eb-2492-d35dfed491cc
md"""## String Comparison

The relational operators work on strings. To see if two strings are equal, use `==`:

```julia
word = "Pineapple"
if word == "banana"
	println("All right, bananas.")
end
```

Other relational operations are useful for putting words in alphabetical order:

```julia
if word < "banana"
	println("Your word, $word, comes before banana.")
elseif word > "banana"
	println("Your word, $word, comes after banana.")
else
	println("All right, bananas.")
end
```

Julia does not handle uppercase and lowercase letters the same way people do. All the uppercase letters come before all the lowercase letters, so:

```
Your word, Pineapple, comes before banana.
```

!!! tip
    A common way to address this problem is to convert strings to a standard format, such as all lowercase, before performing the comparison.
"""

# ╔═╡ a05f7e2a-79ff-11eb-16f9-db358e178910
md"""## Debugging

When you use indices to traverse the values in a sequence, it is tricky to get the beginning and end of the traversal right. Here is a function that is supposed to compare two words and return `true` if one of the words is the reverse of the other, but it contains two errors:

```julia
function isreverse(word1, word2)
	if length(word1) != length(word2)
		return false
    end
	i = firstindex(word1)
	j = lastindex(word2)
	while j >= 0
		j = prevind(word2, j)
		if word1[i] != word2[j]
			return false
		end
		i = nextind(word1, i)
	end
	true
end
```

The first `if` statement checks whether the words are the same length. If not, we can return `false` immediately. Otherwise, for the rest of the function, we can assume that the words are the same length. This is an example of the guardian pattern; see “Checking Types”.

`i` and `j` are indices: `i` traverses `word1` forward while `j` traverses `word2` backward. If we find two letters that don’t match, we can return false immediately. If we get through the whole loop and all the letters match, we return `true`.

The function `lastindex` returns the last valid byte index of a string and `prevind` finds the previous valid index of a character.

If we test this function with the words `“pots”` and `“stop”`, we expect the return value `true`, but we get `false`:

```julia
julia> isreverse("pots", "stop")
false
```

For debugging this kind of error, my first move is to print the values of the indices:

```julia
while j >= 0
	j = prevind(word2, j)
	@show i j
	if word1[i] != word2[j]
```

Now when I run the program again, I get more information:

```julia
julia> isreverse("pots", "stop")
i = 1
j = 3
false
```

The first time through the loop, the value of `j` is `3`, but it has to be `4`. This can be fixed by moving `j = prevind(word2, j)` to the end of the `while` loop.

If I fix that error and run the program again, I get:

```julia
julia> isreverse("pots", "stop")
i = 1
j = 4
i = 2
j = 3
i = 3
j = 2
i = 4
j = 1
i = 5
j = 0
ERROR: BoundsError: attempt to access "pots"
at index [5]
```

This time a `BoundsError` has been thrown. The value of `i` is `5`, which is out of range for the string `"pots"`.
"""

# ╔═╡ c657df0e-7a00-11eb-06e7-b1839a01a3b1
md"""#### Exercise 8-6

Run the program on paper, changing the values of `i` and `j` during each iteration. Find and fix the second error in this function.
"""

# ╔═╡ e8baffa4-7a00-11eb-3373-3d1dc1e6a952
md"""## Glossary

*sequence*:
An ordered collection of values where each value is identified by an integer index.

*ASCII standard*:
A character encoding standard for electronic communication specifying 128 characters.

*Unicode standard*:
A computing industry standard for the consistent encoding, representation, and handling of text expressed in most of the world’s writing systems.

*index*:
An integer value used to select an item in a sequence, such as a character in a string. In Julia indices start from 1.

*UTF-8 encoding*:
A variable-width character encoding capable of encoding all 1,112,064 valid code points in Unicode using one to four 8-bit bytes.

*traverse*:
To iterate through the items in a sequence, performing a similar operation on each.
slice

*empty string*:
A string with no characters and length 0, represented by two quotation marks.

*immutable*:
The property of a sequence whose items cannot be changed.

*string interpolation*:
The process of evaluating a string containing one or more placeholders, yielding a result in which the placeholders are replaced with their corresponding values.

*search*:
A pattern of traversal that stops when it finds what it is looking for.

*counter*:
A variable used to count something, usually initialized to zero and then incremented.
"""

# ╔═╡ 55115388-7a01-11eb-2a9e-a769127bc90d
md"""## Exercises

#### Exercise 8-7

Read the [documentation of the string functions](https://docs.julialang.org/en/v1/manual/strings/). You might want to experiment with some of them to make sure you understand how they work. `strip` and `replace` are particularly useful.

The documentation uses a syntax that might be confusing. For example, in `search(string::AbstractString, chars::Chars, [start::Integer])`, the brackets indicate optional arguments. So, `string` and `chars` are required, but `start` is optional.
"""

# ╔═╡ bcb82726-7a01-11eb-1b75-7136ad9adfb6
md"""#### Exercise 8-8

There is a built-in function called count that is similar to the function in “Looping and Counting”. Read the documentation of this function and use it to count the number of `a`’s in `"banana"`.
"""

# ╔═╡ d33a25a0-7a01-11eb-35cc-573b861aadbb
md"""#### Exercise 8-9

A string slice can take a third index. The first specifies the start, the third the end, and the second the “step size”; that is, the number of spaces between successive characters. A step size of 2 means every other character; 3 means every third, etc. For example:

```julia
julia> fruit = "banana"
"banana"
julia> fruit[1:2:6]
"bnn"
```

A step size of `-1` goes through the word backward, so the slice `[end:-1:begin]` generates a reversed string.

Use this idiom to write a one-line version of `ispalindrome` from “Exercise 6-6”.
"""

# ╔═╡ 0d718e70-7a02-11eb-27e0-0f145b163d55
md"""#### Exercise 8-10

The following functions are all *intended* to check whether a string contains any lowercase letters, but at least some of them are wrong. For each function, describe what the function actually does (assuming that the parameter is a string).

```julia
function anylowercase1(s)
	for c in s
		if islowercase(c)
			return true
		else
			return false
		end
	end
end

function anylowercase2(s)
	for c in s
		if islowercase('c')
			return "true"
		else
			return "false"
		end
	end
end

function anylowercase3(s)
	for c in s
		flag = islowercase(c)
	end
	flag
end

function anylowercase4(s)
	flag = false
	for c in s
		flag = flag || islowercase(c)
	end
	flag
end

function anylowercase5(s)
	for c in s
		if !islowercase(c)
			return false
		end
	end
	true
end
```
"""

# ╔═╡ 7e1acfec-7a02-11eb-2b50-ff3f9e93545e
md"""#### Exercise 8-11

A Caesar cypher is a weak form of encryption that involves “rotating” each letter by a fixed number of places. To rotate a letter means to shift it through the alphabet, wrap‐ ping around to the beginning if necessary, so *`A`* rotated by 3 is *`D`* and *`Z`* rotated by 1 is *`A`*.

To rotate a word, rotate each letter by the same amount. For example, `"cheer"` rota‐ ted by 7 is `"jolly"` and `"melon"` rotated by –10 is `"cubed"`. In the movie *2001: A Space Odyssey*, the ship’s computer is called “HAL,” which is “IBM” rotated by –1.

Write a function called `rotateword` that takes a string and an integer as parameters, and returns a new string that contains the letters from the original string rotated by the given amount.

!!! tip
    You might want to use the built-in functions `Int`, which converts a character to a numeric code, and `Char`, which converts numeric codes to characters. Letters of the alphabet are encoded in alphabetical order, so, for example:

    ```julia
    julia> Int('c') - Int('a')
    2
    ```

    because *c* is the third letter of the alphabet. But beware—the numeric codes for uppercase letters are different:

    ```julia
    julia> Char(Int('A') + 32)
    'a': ASCII/Unicode U+0061 (category Ll: Letter, lower case)
    ```

Potentially offensive jokes on the internet are sometimes encoded in ROT13, which is a Caesar cypher with rotation 13. If you are not easily offended, find and decode some of them.
"""

# ╔═╡ Cell order:
# ╟─2444e25e-76e3-11eb-021a-e9e976b267bb
# ╟─d77c9e58-76e2-11eb-27a0-23914499f29f
# ╟─3eeafd32-76e3-11eb-01dc-53af6f26eaf7
# ╟─b29b3ca6-76e3-11eb-0970-1bee030205f2
# ╟─ade32e66-76e4-11eb-165e-13867b2cdc2d
# ╟─ac377fee-76e5-11eb-1736-9190b5c92214
# ╟─20c034f2-76e6-11eb-1f6d-891ef3f93b29
# ╟─2c1d4452-76e6-11eb-1cc4-0be827cdb64d
# ╟─84220d54-76e6-11eb-24d3-b7085b8766b8
# ╟─c170dc56-76e6-11eb-2bd1-050911e1915f
# ╟─559c8a42-76e7-11eb-371c-391c5deba644
# ╟─67110424-76e7-11eb-14c3-e195649e6d41
# ╟─e6a4882e-79fc-11eb-3957-95357985303f
# ╟─36ac2624-79fd-11eb-1099-910ca7cea1c9
# ╟─e04750c8-79fd-11eb-0249-ed3fb7e031b2
# ╟─019a3968-79fe-11eb-206a-3d92d9352524
# ╟─5af18f32-79fe-11eb-1860-0f08e3d25a3f
# ╟─8d5f176e-79fe-11eb-29fc-813db53b35a6
# ╟─e80e049a-79fe-11eb-2ade-83a99e44e5a5
# ╟─5487357e-79ff-11eb-2492-d35dfed491cc
# ╟─a05f7e2a-79ff-11eb-16f9-db358e178910
# ╟─c657df0e-7a00-11eb-06e7-b1839a01a3b1
# ╟─e8baffa4-7a00-11eb-3373-3d1dc1e6a952
# ╟─55115388-7a01-11eb-2a9e-a769127bc90d
# ╟─bcb82726-7a01-11eb-1b75-7136ad9adfb6
# ╟─d33a25a0-7a01-11eb-35cc-573b861aadbb
# ╟─0d718e70-7a02-11eb-27e0-0f145b163d55
# ╟─7e1acfec-7a02-11eb-2b50-ff3f9e93545e


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests related to orbit representations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/orbit/representations/rv.jl
# =======================================

# Functions: kepler_to_rv and rv_to_kepler
# ----------------------------------------

################################################################################
#                                 TEST RESULTS
################################################################################
#
# Scenario 01
# ===========
#
# Example 2-5: Finding position and velocity vectors (COE2RV Test Case) [1, p.
# 119-120].
#
# Cartesian representation:
#
#     r = 6525.344    I + 6861.535    J + 6449.125    K km
#     v =    4.902276 I +    5.533124 J -    1.975709 K km
#
# Orbit elements
#
#     ╔═════════════════╦══════════════╗
#     ║    Parameter    ║    Values    ║
#     ╠═════════════════╬══════════════╣
#     ║ p               ║ 11067.790 km ║
#     ║ Eccentricity    ║ 0.83285      ║
#     ║ Inclination     ║ 87.87°       ║
#     ║ RAAN            ║ 227.89°      ║
#     ║ Arg. of Perigee ║ 53.38°       ║
#     ║ True Anomaly    ║ 92.335°      ║
#     ╚═════════════════╩══════════════╝
#
################################################################################

@testset "Function kepler_to_rv and rv_to_kepler" begin

    ## kepler_to_rv
    ## ============

    p    = 11067.790*1000
    e    = 0.83285
    i    = 87.87*pi/180
    RAAN = 227.89*pi/180
    w    = 53.38*pi/180
    f    = 92.335*pi/180
    a    = p/(1-e^2)

    (r_i, v_i) = kepler_to_rv(a, e, i, RAAN, w, f)

    @test r_i[1]/1000 ≈ +6525.344 atol=5e-2
    @test r_i[2]/1000 ≈ +6861.535 atol=5e-2
    @test r_i[3]/1000 ≈ +6449.125 atol=5e-2

    @test v_i[1]/1000 ≈ +4.902276 atol=1e-4
    @test v_i[2]/1000 ≈ +5.533124 atol=1e-4
    @test v_i[3]/1000 ≈ -1.975709 atol=1e-4

    ## rv_to_kepler
    ## ============

    r_i = [6525.344; 6861.535; 6449.125]*1000
    v_i = [4.902276; 5.533124; -1.975709]*1000

    oe = rv_to_kepler(r_i..., v_i...)

    a, e, i, RAAN, w, f = oe.a, oe.e, oe.i, oe.Ω, oe.ω, oe.f
    p = a*(1-e^2)

    @test p/1000      ≈ 11067.790 atol=5e-2
    @test e           ≈ 0.83285   atol=1e-5
    @test i*180/pi    ≈ 87.87     atol=1e-2
    @test RAAN*180/pi ≈ 227.89    atol=1e-2
    @test w*180/pi    ≈ 53.38     atol=1e-2
    @test f*180/pi    ≈ 92.335    atol=1e-3
end

@testset "Issue #25 - Special cases in kepler_to_rv" begin
    # Circular and equatorial
    # =======================

	orb_i = KeplerianElements(0.0, 42164e3, 0, 0, 0, 0, 0 )
    orb_o = sv_to_kepler( kepler_to_sv( orb_i ) )
	orb_e = KeplerianElements(0.0, 42164e3, 0, 0, 0, 0, 0 )

    @test orb_o.a ≈ orb_e.a   atol = 1e-7
    @test orb_o.e ≈ orb_e.e   atol = 1e-7
    @test orb_o.i ≈ orb_e.i   atol = 1e-7
    @test orb_o.Ω ≈ orb_e.Ω   atol = 1e-7
    @test orb_o.ω ≈ orb_e.ω   atol = 1e-7
    @test orb_o.f ≈ orb_e.f   atol = 1e-7

	orb_i = KeplerianElements(0.0, 42164e3, 0, 0, 10*π/180, 11*π/180, 12*π/180 )
    orb_o = sv_to_kepler( kepler_to_sv( orb_i ) )
    orb_e = KeplerianElements(0.0, 42164e3, 0, 0, 0, 0, (10+11+12)*π/180 )

    @test orb_o.a ≈ orb_e.a   atol = 1e-7
    @test orb_o.e ≈ orb_e.e   atol = 1e-7
    @test orb_o.i ≈ orb_e.i   atol = 1e-7
    @test orb_o.Ω ≈ orb_e.Ω   atol = 1e-7
    @test orb_o.ω ≈ orb_e.ω   atol = 1e-7
    @test orb_o.f ≈ orb_e.f   atol = 1e-7

    # Elliptical equatorial
    # =====================

    orb_i = KeplerianElements(0.0, 42164e3, 0.05, 0, 10*π/180, 11*π/180, 12*π/180 )
    orb_o = sv_to_kepler( kepler_to_sv( orb_i ) )
    orb_e = KeplerianElements(0.0, 42164e3, 0.05, 0, 0, (10+11)*π/180, 12*π/180 )

    @test orb_o.a ≈ orb_e.a   atol = 1e-7
    @test orb_o.e ≈ orb_e.e   atol = 1e-7
    @test orb_o.i ≈ orb_e.i   atol = 1e-7
    @test orb_o.Ω ≈ orb_e.Ω   atol = 1e-7
    @test orb_o.ω ≈ orb_e.ω   atol = 1e-7
    @test orb_o.f ≈ orb_e.f   atol = 1e-7

    # Circular inclined
    # =================

    orb_i = KeplerianElements(0.0, 42164e3, 0, π/4, 10*π/180, 11*π/180, 12*π/180 )
    orb_o = sv_to_kepler( kepler_to_sv( orb_i ) )
    orb_e = KeplerianElements(0.0, 42164e3, 0, π/4, 10*π/180, 0, (11+12)*π/180 )

    @test orb_o.a ≈ orb_e.a   atol = 1e-7
    @test orb_o.e ≈ orb_e.e   atol = 1e-7
    @test orb_o.i ≈ orb_e.i   atol = 1e-7
    @test orb_o.Ω ≈ orb_e.Ω   atol = 1e-7
    @test orb_o.ω ≈ orb_e.ω   atol = 1e-7
    @test orb_o.f ≈ orb_e.f   atol = 1e-7
end


text:
"""
    nfeasible(x::LogicalCombo; feasible=true)

    Returns a count of the number of feasible outcomes in x.

* feasible:

#### Example
```
julia> nfeasible(logicalparse("a,b in 1:3; a!=b", verbose = false))
6
julia> nfeasible(logicalparse("a,b in 1:3; a!=b", verbose = false), feasible=false)
3
```
"""
nfeasible(x::LogicalCombo; feasible::Bool=true) =
  (feasible ? sum(x.logical) : size(x,1) - sum(x.logical))

nfeasible(x::LogicalCombo, feasible::Bool) = nfeasible(x::LogicalCombo; feasible=feasible)

variablerange(x::LogicalCombo; feasible=true) = [unique(x[:,i,feasible]) for i in 1:size(x,2)]

"""
    percievedfeasible(x::LogicalCombo)

    Returns a count of the number of feasible outcomes that could be infered if
    one were just looking at the number of possible values each variable can take.

* feasible:

#### Example
```
julia> percievedfeasible(logicalparse("a,b in 1:3; a=b", verbose = false))
9
julia> percievedfeasible(logicalparse("a,b in 1:3; a=1", verbose = false))
3
```
"""
percievedfeasible(x::LogicalCombo) = prod(length.(variablerange(x)))


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

# Name functions
Jan = Dates.DateTime(2013,1,1) #Tuesday
Feb = Dates.DateTime(2013,2,2) #Saturday
Mar = Dates.DateTime(2013,3,3) #Sunday
Apr = Dates.DateTime(2013,4,4) #Thursday
May = Dates.DateTime(2013,5,5) #Sunday
Jun = Dates.DateTime(2013,6,7) #Friday
Jul = Dates.DateTime(2013,7,7) #Sunday
Aug = Dates.DateTime(2013,8,8) #Thursday
Sep = Dates.DateTime(2013,9,9) #Monday
Oct = Dates.DateTime(2013,10,10) #Thursday
Nov = Dates.DateTime(2013,11,11) #Monday
Dec = Dates.DateTime(2013,12,11) #Wednesday
monthnames = ["January","February","March","April",
                "May","June","July","August","September",
                "October","November","December"]
daysofweek = [Dates.Tue,Dates.Sat,Dates.Sun,Dates.Thu,Dates.Sun,Dates.Fri,
              Dates.Sun,Dates.Thu,Dates.Mon,Dates.Thu,Dates.Mon,Dates.Wed]
dows = ["Tuesday","Saturday","Sunday","Thursday","Sunday","Friday",
        "Sunday","Thursday","Monday","Thursday","Monday","Wednesday"]
daysinmonth = [31,28,31,30,31,30,31,31,30,31,30,31]
for (i,dt) in enumerate([Jan,Feb,Mar,Apr,May,Jun,Jul,Aug,Sep,Oct,Nov,Dec])
    @test Dates.month(dt) == i
    @test Dates.monthname(dt) == monthnames[i]
    @test Dates.monthname(i) == monthnames[i]
    @test Dates.monthabbr(dt) == monthnames[i][1:3]
    @test Dates.monthabbr(i) == monthnames[i][1:3]
    @test Dates.dayofweek(dt) == daysofweek[i]
    @test Dates.dayname(dt) == dows[i]
    @test Dates.dayname(Dates.dayofweek(dt)) == dows[i]
    @test Dates.dayabbr(dt) == dows[i][1:3]
    @test Dates.dayabbr(Dates.dayofweek(dt)) == dows[i][1:3]
    @test Dates.daysinmonth(dt) == daysinmonth[i]
end

# Customizing locale
const french_daysofweek = Dict(1=>"Lundi",2=>"Mardi",3=>"Mercredi",4=>"Jeudi",
                               5=>"Vendredi",6=>"Samedi",7=>"Dimanche")
Dates.VALUETODAYOFWEEK["french"] = french_daysofweek
@test Dates.dayname(Nov;locale="french") == "Lundi"
@test Dates.dayname(Jan;locale="french") == "Mardi"
@test Dates.dayname(Dec;locale="french") == "Mercredi"
@test Dates.dayname(Apr;locale="french") == "Jeudi"
@test Dates.dayname(Jun;locale="french") == "Vendredi"
@test Dates.dayname(Feb;locale="french") == "Samedi"
@test Dates.dayname(May;locale="french") == "Dimanche"

const french_months = Dict(1=>"janvier",2=>"février",3=>"mars",4=>"avril",5=>"mai",6=>"juin",
                           7=>"juillet",8=>"août",9=>"septembre",10=>"octobre",11=>"novembre",12=>"décembre")
Dates.VALUETOMONTH["french"] = french_months
@test Dates.monthname(Jan;locale="french") == "janvier"
@test Dates.monthname(Feb;locale="french") == "février"
@test Dates.monthname(Mar;locale="french") == "mars"
@test Dates.monthname(Apr;locale="french") == "avril"
@test Dates.monthname(May;locale="french") == "mai"
@test Dates.monthname(Jun;locale="french") == "juin"
@test Dates.monthname(Jul;locale="french") == "juillet"
@test Dates.monthname(Aug;locale="french") == "août"
@test Dates.monthname(Sep;locale="french") == "septembre"
@test Dates.monthname(Oct;locale="french") == "octobre"
@test Dates.monthname(Nov;locale="french") == "novembre"
@test Dates.monthname(Dec;locale="french") == "décembre"

@test Dates.isleapyear(Dates.DateTime(1900)) == false
@test Dates.isleapyear(Dates.DateTime(2000)) == true
@test Dates.isleapyear(Dates.DateTime(2004)) == true
@test Dates.isleapyear(Dates.DateTime(2008)) == true
@test Dates.isleapyear(Dates.DateTime(0)) == true
@test Dates.isleapyear(Dates.DateTime(1)) == false
@test Dates.isleapyear(Dates.DateTime(-1)) == false
@test Dates.isleapyear(Dates.DateTime(4)) == true
@test Dates.isleapyear(Dates.DateTime(-4)) == true

@test Dates.daysinyear(2000) == 366
@test Dates.daysinyear(2001) == 365
@test Dates.daysinyear(2000) == 366
@test Dates.daysinyear(2001) == 365

@test Dates.daysinyear(Dates.Date(2000)) == 366
@test Dates.daysinyear(Dates.Date(2001)) == 365
@test Dates.daysinyear(Dates.DateTime(2000)) == 366
@test Dates.daysinyear(Dates.DateTime(2001)) == 365

# Days of week from Monday = 1 to Sunday = 7
@test Dates.dayofweek(Dates.DateTime(2013,12,22)) == 7
@test Dates.dayofweek(Dates.DateTime(2013,12,23)) == 1
@test Dates.dayofweek(Dates.DateTime(2013,12,24)) == 2
@test Dates.dayofweek(Dates.DateTime(2013,12,25)) == 3
@test Dates.dayofweek(Dates.DateTime(2013,12,26)) == 4
@test Dates.dayofweek(Dates.DateTime(2013,12,27)) == 5
@test Dates.dayofweek(Dates.DateTime(2013,12,28)) == 6
@test Dates.dayofweek(Dates.DateTime(2013,12,29)) == 7

# There are 5 Sundays in December, 2013
@test Dates.daysofweekinmonth(Dates.DateTime(2013,12,1)) == 5
# There are 4 Sundays in November, 2013
@test Dates.daysofweekinmonth(Dates.DateTime(2013,11,24)) == 4

@test Dates.dayofweekofmonth(Dates.DateTime(2013,12,1)) == 1
@test Dates.dayofweekofmonth(Dates.DateTime(2013,12,8)) == 2
@test Dates.dayofweekofmonth(Dates.DateTime(2013,12,15)) == 3
@test Dates.dayofweekofmonth(Dates.DateTime(2013,12,22)) == 4
@test Dates.dayofweekofmonth(Dates.DateTime(2013,12,29)) == 5

@test Dates.dayofyear(2000,1,1) == 1
@test Dates.dayofyear(2004,1,1) == 1
@test Dates.dayofyear(20013,1,1) == 1
# Leap year
@test Dates.dayofyear(2000,12,31) == 366
# Non-leap year
@test Dates.dayofyear(2001,12,31) == 365

@test Dates.dayofyear(Dates.DateTime(2000,1,1)) == 1
@test Dates.dayofyear(Dates.DateTime(2004,1,1)) == 1
@test Dates.dayofyear(Dates.DateTime(20013,1,1)) == 1
# Leap year
@test Dates.dayofyear(Dates.DateTime(2000,12,31)) == 366
# Non-leap year
@test Dates.dayofyear(Dates.DateTime(2001,12,31)) == 365
# Test every day of a year
dt = Dates.DateTime(2000,1,1)
for i = 1:366
    @test Dates.dayofyear(dt) == i
    dt += Dates.Day(1)
end
dt = Dates.DateTime(2001,1,1)
for i = 1:365
    @test Dates.dayofyear(dt) == i
    dt += Dates.Day(1)
end

@test Dates.quarterofyear(Dates.Date(2000,1,1))  == 1
@test Dates.quarterofyear(Dates.Date(2000,1,31))  == 1
@test Dates.quarterofyear(Dates.Date(2000,2,1))  == 1
@test Dates.quarterofyear(Dates.Date(2000,2,29))  == 1
@test Dates.quarterofyear(Dates.Date(2000,3,1))  == 1
@test Dates.quarterofyear(Dates.Date(2000,3,31))  == 1
@test Dates.quarterofyear(Dates.Date(2000,4,1)) == 2
@test Dates.quarterofyear(Dates.Date(2000,4,30)) == 2
@test Dates.quarterofyear(Dates.Date(2000,5,1)) == 2
@test Dates.quarterofyear(Dates.Date(2000,5,31)) == 2
@test Dates.quarterofyear(Dates.Date(2000,6,1)) == 2
@test Dates.quarterofyear(Dates.Date(2000,6,30)) == 2
@test Dates.quarterofyear(Dates.Date(2000,7,1)) == 3
@test Dates.quarterofyear(Dates.Date(2000,7,31)) == 3
@test Dates.quarterofyear(Dates.Date(2000,8,1)) == 3
@test Dates.quarterofyear(Dates.Date(2000,8,31)) == 3
@test Dates.quarterofyear(Dates.Date(2000,9,1)) == 3
@test Dates.quarterofyear(Dates.Date(2000,9,30)) == 3
@test Dates.quarterofyear(Dates.Date(2000,10,1)) == 4
@test Dates.quarterofyear(Dates.Date(2000,10,31)) == 4
@test Dates.quarterofyear(Dates.Date(2000,11,1)) == 4
@test Dates.quarterofyear(Dates.Date(2000,11,30)) == 4
@test Dates.quarterofyear(Dates.Date(2000,12,1)) == 4
@test Dates.quarterofyear(Dates.Date(2000,12,31)) == 4

@test Dates.quarterofyear(Dates.DateTime(2000,1,1))  == 1
@test Dates.quarterofyear(Dates.DateTime(2000,1,31))  == 1
@test Dates.quarterofyear(Dates.DateTime(2000,2,1))  == 1
@test Dates.quarterofyear(Dates.DateTime(2000,2,29))  == 1
@test Dates.quarterofyear(Dates.DateTime(2000,3,1))  == 1
@test Dates.quarterofyear(Dates.DateTime(2000,3,31))  == 1
@test Dates.quarterofyear(Dates.DateTime(2000,4,1)) == 2
@test Dates.quarterofyear(Dates.DateTime(2000,4,30)) == 2
@test Dates.quarterofyear(Dates.DateTime(2000,5,1)) == 2
@test Dates.quarterofyear(Dates.DateTime(2000,5,31)) == 2
@test Dates.quarterofyear(Dates.DateTime(2000,6,1)) == 2
@test Dates.quarterofyear(Dates.DateTime(2000,6,30)) == 2
@test Dates.quarterofyear(Dates.DateTime(2000,7,1)) == 3
@test Dates.quarterofyear(Dates.DateTime(2000,7,31)) == 3
@test Dates.quarterofyear(Dates.DateTime(2000,8,1)) == 3
@test Dates.quarterofyear(Dates.DateTime(2000,8,31)) == 3
@test Dates.quarterofyear(Dates.DateTime(2000,9,1)) == 3
@test Dates.quarterofyear(Dates.DateTime(2000,9,30)) == 3
@test Dates.quarterofyear(Dates.DateTime(2000,10,1)) == 4
@test Dates.quarterofyear(Dates.DateTime(2000,10,31)) == 4
@test Dates.quarterofyear(Dates.DateTime(2000,11,1)) == 4
@test Dates.quarterofyear(Dates.DateTime(2000,11,30)) == 4
@test Dates.quarterofyear(Dates.DateTime(2000,12,1)) == 4
@test Dates.quarterofyear(Dates.DateTime(2000,12,31)) == 4

@test Dates.dayofquarter(Dates.Date(2014,1,1)) == 1
@test Dates.dayofquarter(Dates.Date(2014,4,1)) == 1
@test Dates.dayofquarter(Dates.Date(2014,7,1)) == 1
@test Dates.dayofquarter(Dates.Date(2014,10,1)) == 1
@test Dates.dayofquarter(Dates.Date(2014,3,31)) == 90
@test Dates.dayofquarter(Dates.Date(2014,6,30)) == 91
@test Dates.dayofquarter(Dates.Date(2014,9,30)) == 92
@test Dates.dayofquarter(Dates.Date(2014,12,31)) == 92
@test Dates.dayofquarter(Dates.DateTime(2014,1,1)) == 1
@test Dates.dayofquarter(Dates.DateTime(2014,4,1)) == 1
@test Dates.dayofquarter(Dates.DateTime(2014,7,1)) == 1
@test Dates.dayofquarter(Dates.DateTime(2014,10,1)) == 1
@test Dates.dayofquarter(Dates.DateTime(2014,3,31)) == 90
@test Dates.dayofquarter(Dates.DateTime(2014,6,30)) == 91
@test Dates.dayofquarter(Dates.DateTime(2014,9,30)) == 92
@test Dates.dayofquarter(Dates.DateTime(2014,12,31)) == 92

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
# add Blink   # in REPL package mode
using Blink

w = Window()
body!(w, "Hello from Blink!")

loadurl(w, "https://google.com")
load!(w, "ui/app.css")
load!(w, "ui/frameworks/jquery-3.3.1.js")


text:
module ECMMethods

# stdlib stuff
using Random, LinearAlgebra, Printf

# Packages we need
using BasisMatrices, Optim, QuantEcon, Parameters
using BasisMatrices: Degree, Derivative

# Export functions that need to be exposed
export NeoclassicalGrowth, ValueCoeffs, f, df, u, du, duinv, expendables_t,
       SolutionMethod, IterateOnPolicy, VFI_ECM, VFI_EGM, VFI, PFI_ECM, PFI,
       dVFI_ECM, EulEq, build_V, build_k, solve, copy, simulate,
       EulerEquation!, ee_residuals, compute_EV, env_condition_kp, update_v!

# Include our other files
include("model.jl")
include("solutionmethods.jl")
include("sim_resids.jl")

#=

shocks = vec(readcsv("../EE_SHOCKS.csv"))
ncgm = NeoclassicalGrowth()
@time begin
vp = ValueCoeffs(ncgm, Val{2}, IterateOnPolicy())
solve(ncgm, vp, 1e-7, 5000, 1.0, 200, false)
vp_vfi = copy(vp, VFI_ECM())
Profile.clear()
solve(ncgm, vp_vfi, 1e-7, 5000, 1.0, 50, false)
end  # time

=#

end  # Module

max_stars_repo_name: sglyon/CLMMJuliaPythonMatlab

max_issues_repo_name: sglyon/CLMMJuliaPythonMatlab

max_forks_repo_name: sglyon/CLMMJuliaPythonMatlab


text:
using LaserTypes, Test

GCL  = setup_laser(GaussLaser, :atomic; profile = ConstantProfile())
GGL  = setup_laser(GaussLaser, :atomic; profile = GaussProfile(z₀ = 0.0, t₀ = 0.0, τ = 744.145))
G2L  = setup_laser(GaussLaser, :atomic; profile = Cos²Profile(z₀ = 0.0, t₀ = 0.0, τ = 744.145))
GQL  = setup_laser(GaussLaser, :atomic; profile = QuasiRectangularProfile(z₀ = 0.0, t₀ = 0.0, τ = 744.145, Δt = 7441.45))

LGCL = setup_laser(LaguerreGaussLaser, :atomic; profile = ConstantProfile(), p = 1, m = 1)
LGCD = setup_laser(LaguerreGaussLaser, :atomic; profile = ConstantProfile(), p = 1, m = 1, ξx = 1/√2, ξy =  1im/√2)
LGCS = setup_laser(LaguerreGaussLaser, :atomic; profile = ConstantProfile(), p = 1, m = 1, ξx = 1/√2, ξy = -1im/√2)

vec = [GCL, LGCL]

# DO NOT DISTURB THE STRINGS UNLESS YOU KNOW WHAT YOU ARE DOING OR YOU HAVE TIME TO TEST THINGS OUT!

GCL_str =
"""GaussLaser with atomic units
a₀ = 1.0
ω  = 0.0569542
ϕ₀ = 0.0
w₀ = 1.09604e6
with linear polarization
and temporal ConstantProfile"""

GGL_str =
"""GaussLaser with atomic units
a₀ = 1.0
ω  = 0.0569542
ϕ₀ = 0.0
w₀ = 1.09604e6
with linear polarization
and temporal Gaussian profile
centered in z₀ = 0.0 and t₀ = 0.0
with duration of pulse (FWHM) τ = 744.145"""

G2L_str =
"""GaussLaser with atomic units
a₀ = 1.0
ω  = 0.0569542
ϕ₀ = 0.0
w₀ = 1.09604e6
with linear polarization
and temporal Cos² profile
centered in z₀ = 0.0 and t₀ = 0.0
with duration of pulse τ = 744.145"""

GQL_str =
"""GaussLaser with atomic units
a₀ = 1.0
ω  = 0.0569542
ϕ₀ = 0.0
w₀ = 1.09604e6
with linear polarization
and temporal Quasi rectangular profile
centered in z₀ = 0.0 and t₀ = 0.0
with exponential decay width τ = 744.145
length of the plateau Δt = 7441.45"""

LGCL_str =
"""LaguerreGaussLaser with atomic units
a₀ = 1.0
ω  = 0.0569542
ϕ₀ = 0.0
w₀ = 1.09604e6
p = 1 and m = 1
with linear polarization
and temporal ConstantProfile"""

LGCD_str =
"""LaguerreGaussLaser with atomic units
a₀ = 1.0
ω  = 0.0569542
ϕ₀ = 0.0
w₀ = 1.09604e6
p = 1 and m = 1
with right-handed polarization
and temporal ConstantProfile"""

LGCS_str =
"""LaguerreGaussLaser with atomic units
a₀ = 1.0
ω  = 0.0569542
ϕ₀ = 0.0
w₀ = 1.09604e6
p = 1 and m = 1
with left-handed polarization
and temporal ConstantProfile"""

vec_str =
"""
2-element Vector{LaserTypes.AbstractLaser}:
 Gaussian laser: a₀ = 1.0, ω = 0.0569542, ϕ₀ = 0.0, w₀ = 1.09604e6, temporal ConstantProfile
 Laguerre-Gauss laser: a₀ = 1.0, ω = 0.0569542, ϕ₀ = 0.0, w₀ = 1.09604e6, p = 1, m = 1, temporal ConstantProfile"""

@testset "Printing/Show" begin
    io = IOBuffer()

    # Gauss Tests
    # Constant Profile
    show(io, MIME"text/plain"(), GCL)
    @test String(take!(io)) == GCL_str

    # Gaussian Profile
    show(io, MIME"text/plain"(), GGL)
    @test String(take!(io)) == GGL_str

    # Cos² Profile
    show(io, MIME"text/plain"(), G2L)
    @test String(take!(io)) == G2L_str

    # Quasi Rectangular Profile
    show(io, MIME"text/plain"(), GQL)
    @test String(take!(io)) == GQL_str

    # Laguerre Gauss Laser
    # Linearly Polarized
    show(io, MIME"text/plain"(), LGCL)
    @test String(take!(io)) == LGCL_str

    # Right Circularly Polarized
    show(io, MIME"text/plain"(), LGCD)
    @test String(take!(io)) == LGCD_str

    # Left Circularly Polarized
    show(io, MIME"text/plain"(), LGCS)
    @test String(take!(io)) == LGCS_str

    # Brief Printing Test
    # Vector Test
    show(io, MIME"text/plain"(), vec)
    @test String(take!(io)) == vec_str
end


text:
module Theta0
function f(q::Vector{T}, q_dot::Vector{T}, xi::T) where T
    l1 = q[1]
    l2 = q[2]
    l3 = q[3]
    l1_dot = q_dot[1]
    l2_dot = q_dot[2]
    l3_dot = q_dot[3]

    [8.8412689977681e-5*xi.^2.*(0.5*l1 - 5.0e-5).*(l2 - 0.002) + 3.31547587416304e-5*xi.^2.*(0.5*l1 - 5.0e-5).*(l3 - 0.0015) - 3.03918621798279e-5*xi.^2.*(l1 - 0.0001).^2 - 2.41125518120948e-5*xi.^2.*(0.5*l2 - 0.001).*(l3 - 0.0015) - 1.60750345413965e-5*xi.^2.*(l2 - 0.002).^2 + 0.5*(l3 - 0.0015).^2.*(-4.52110346476778e-6*xi.^2 - 1.04773789644241e-9) + 1.0 -2.32023137992695e-5*xi.^2.*(0.5*l1 - 5.0e-5).*(l2 - 0.002) - 8.70086767472609e-6*xi.^2.*(0.5*l1 - 5.0e-5).*(l3 - 0.0015) + 7.97579536849891e-6*xi.^2.*(l1 - 0.0001).^2 + 6.32790376343716e-6*xi.^2.*(0.5*l2 - 0.001).*(l3 - 0.0015) + 4.2186025089581e-6*xi.^2.*(l2 - 0.002).^2 + 5.93240977822234e-7*xi.^2.*(l3 - 0.0015).^2 2.96701739469302*xi.*(0.5*l1 - 5.0e-5).*(l2 - 0.002) - 4.79372586530627*xi.*(0.5*l1 - 5.0e-5).*(l3 - 0.0015) + 0.456677117653311*xi.*(l1 - 0.0001).^2 - 0.00779639175257732*xi.*(l1 - 0.0001) - 5.26092039536614*xi.*(0.5*l2 - 0.001).*(l3 - 0.0015) + 0.57347575016828*xi.*(l2 - 0.002).^2 + 0.00567010309278351*xi.*(l2 - 0.002) + 2.5136615651681*xi.*(l3 - 0.0015).^2 + 0.00212628865979382*xi.*(l3 - 0.0015);
8.8412689977681e-5*xi.^2.*(0.5*l1 - 5.0e-5).*(l2 - 0.002) + 3.31547587416304e-5*xi.^2.*(0.5*l1 - 5.0e-5).*(l3 - 0.0015) - 3.03918621798279e-5*xi.^2.*(l1 - 0.0001).^2 - 2.41125518120948e-5*xi.^2.*(0.5*l2 - 0.001).*(l3 - 0.0015) - 1.60750345413965e-5*xi.^2.*(l2 - 0.002).^2 + 0.5*(l3 - 0.0015).^2.*(-4.52110346476778e-6*xi.^2 - 1.04773789644241e-9) + 1.0 -2.32023137992695e-5*xi.^2.*(0.5*l1 - 5.0e-5).*(l2 - 0.002) - 8.70086767472609e-6*xi.^2.*(0.5*l1 - 5.0e-5).*(l3 - 0.0015) + 7.97579536849891e-6*xi.^2.*(l1 - 0.0001).^2 + 6.32790376343716e-6*xi.^2.*(0.5*l2 - 0.001).*(l3 - 0.0015) + 4.2186025089581e-6*xi.^2.*(l2 - 0.002).^2 + 5.93240977822234e-7*xi.^2.*(l3 - 0.0015).^2 2.96701739469302*xi.*(0.5*l1 - 5.0e-5).*(l2 - 0.002) - 4.79372586530627*xi.*(0.5*l1 - 5.0e-5).*(l3 - 0.0015) + 0.456677117653311*xi.*(l1 - 0.0001).^2 - 0.00779639175257732*xi.*(l1 - 0.0001) - 5.26092039536614*xi.*(0.5*l2 - 0.001).*(l3 - 0.0015) + 0.57347575016828*xi.*(l2 - 0.002).^2 + 0.00567010309278351*xi.*(l2 - 0.002) + 2.5136615651681*xi.*(l3 - 0.0015).^2 + 0.00212628865979382*xi.*(l3 - 0.0015);
8.8412689977681e-5*xi.^2.*(0.5*l1 - 5.0e-5).*(l2 - 0.002) + 3.31547587416304e-5*xi.^2.*(0.5*l1 - 5.0e-5).*(l3 - 0.0015) - 3.03918621798279e-5*xi.^2.*(l1 - 0.0001).^2 - 2.41125518120948e-5*xi.^2.*(0.5*l2 - 0.001).*(l3 - 0.0015) - 1.60750345413965e-5*xi.^2.*(l2 - 0.002).^2 + 0.5*(l3 - 0.0015).^2.*(-4.52110346476778e-6*xi.^2 - 1.04773789644241e-9) + 1.0 -2.32023137992695e-5*xi.^2.*(0.5*l1 - 5.0e-5).*(l2 - 0.002) - 8.70086767472609e-6*xi.^2.*(0.5*l1 - 5.0e-5).*(l3 - 0.0015) + 7.97579536849891e-6*xi.^2.*(l1 - 0.0001).^2 + 6.32790376343716e-6*xi.^2.*(0.5*l2 - 0.001).*(l3 - 0.0015) + 4.2186025089581e-6*xi.^2.*(l2 - 0.002).^2 + 5.93240977822234e-7*xi.^2.*(l3 - 0.0015).^2 2.96701739469302*xi.*(0.5*l1 - 5.0e-5).*(l2 - 0.002) - 4.79372586530627*xi.*(0.5*l1 - 5.0e-5).*(l3 - 0.0015) + 0.456677117653311*xi.*(l1 - 0.0001).^2 - 0.00779639175257732*xi.*(l1 - 0.0001) - 5.26092039536614*xi.*(0.5*l2 - 0.001).*(l3 - 0.0015) + 0.57347575016828*xi.*(l2 - 0.002).^2 + 0.00567010309278351*xi.*(l2 - 0.002) + 2.5136615651681*xi.*(l3 - 0.0015).^2 + 0.00212628865979382*xi.*(l3 - 0.0015)]
end
end


text:
"""
    array(p::TensorProto)

Return `p` as an `Array` of the correct type. Second argument can be used to change type of the returned array
"""
function array(p::TensorProto, wrap=Array)
    # Copy pasted from jl
    # Can probably be cleaned up a bit
    # TODO: Add missing datatypes...
    if p.data_type === TensorProto_DataType.INT64
        if hasproperty(p, :int64_data) && !isempty(p.int64_data)
            return reshape(reinterpret(Int64, p.int64_data), reverse(p.dims)...) |> wrap
        end
        return reshape(reinterpret(Int64, p.raw_data), reverse(p.dims)...) |> wrap
    end

    if p.data_type === TensorProto_DataType.INT32
        if hasproperty(p, :int32_data) && !isempty(p.int32_data)
            return reshape(p.int32_data , reverse(p.dims)...) |> wrap
        end
        return reshape(reinterpret(Int32, p.raw_data), reverse(p.dims)...) |> wrap
    end

    if p.data_type === TensorProto_DataType.INT8
        return reshape(reinterpret(Int8, p.raw_data), reverse(p.dims)...) |> wrap
    end

    if p.data_type === TensorProto_DataType.DOUBLE
        if hasproperty(p, :double_data) && !isempty(p.double_data)
            return reshape(p.double_data , reverse(p.dims)...) |> wrap
        end
        return reshape(reinterpret(Float64, p.raw_data), reverse(p.dims)...) |> wrap
    end

    if p.data_type === TensorProto_DataType.FLOAT
        if hasproperty(p,:float_data) && !isempty(p.float_data)
            return reshape(reinterpret(Float32, p.float_data), reverse(p.dims)...) |> wrap
        end
        return reshape(reinterpret(Float32, p.raw_data), reverse(p.dims)...) |> wrap
    end

    if p.data_type === TensorProto_DataType.FLOAT16
        return reshape(reinterpret(Float16, p.raw_data), reverse(p.dims)...) |> wrap
    end
end

Base.size(vip::ValueInfoProto) = size(vip._type)
Base.size(tp::TypeProto) = size(tp.tensor_type)
Base.size(tp::TensorProto) = tp.dims
Base.size(tp_t::TypeProto_Tensor) = hasproperty(tp_t, :shape) ? size(tp_t.shape) : missing
Base.size(tsp::TensorShapeProto) = size.(Tuple(reverse(tsp.dim)))
Base.size(tsp_d::TensorShapeProto_Dimension) = hasproperty(tsp_d, :dim_value) ? tsp_d.dim_value : missing

"""
    attribute(p::AttributeProto)

Return attribute in `p` as a name => value pair.
"""
function attribute(p::AttributeProto)
    # Copy paste from ONNX.jl
    if (p._type != 0)
        field = [:f, :i, :s, :t, :g, :floats, :ints, :strings, :tensors, :graphs][p._type]
        if field === :s
            return Symbol(p.name) => String(getproperty(p, field))
        elseif  field === :strings
            return Symbol(p.name) => String.(getproperty(p, field))
        end
        return Symbol(p.name) => getproperty(p, field)
    end
end

Base.Dict(pa::AbstractVector{AttributeProto}) = Dict(attribute(p) for p in pa)


text:
push!(LOAD_PATH,"../src/")
using Documenter, IntervalAnalysis

makedocs(sitename="Interval Analysis",doctest=true,modules=[IntervalAnalysis])


text:
"""
add_Y_diag!(Y,Ubar,Vbar,u,v,x,d)

Inplace computation of `Y[i,:] = Vbar'*u + Ubar'*v + d[i]*x`.
"""
function add_Y_diag!(Y,Ubar,Vbar,u,v,x,d)
    @inbounds for k = 1:size(Vbar,2)
        Y[k] = d*x[k]
        @inbounds for i = 1:size(Vbar,1)
            Y[k] += Vbar[i,k]*u[i] + Ubar[i,k]*v[i]
        end
    end
end
"""
add_Y_tril_diag!(Y,w,Ubar,x,d)

Inplace computation of `Y[i,:] = w'*Ubar + ds[i]*x`.
"""
function add_Y_tri_diag!(Y,w,Ubar,x,d)
    @inbounds for k = 1:size(Ubar,2)
        Y[k] = d*x[k]
        @inbounds for i = 1:size(Ubar,1)
            Y[k] += Ubar[i,k]*w[i]
        end
    end
end
"""
add_Y_tri!(Y,Wbar,u)

Inplace computation of `Y[i,:] = Wbar'*u`.
"""
function add_Y_tri!(Y,Wbar,u)
    for k = 1:size(Wbar,2)
        Y[k] = 0.0
        for i = 1:size(Wbar,1)
            Y[k] += Wbar[i,k]*u[i]
        end
    end
end

"""
add_Y!(Y,Ubar,Vbar,u,v)

Inplace computation of `Y[i,:] = Vbar'*u + Ubar'*v`.
"""
function add_Y!(Y,Ubar,Vbar,u,v)
    @inbounds for k = 1:size(Vbar,2)
        Y[k] = 0.0
        @inbounds for i = 1:size(Vbar,1)
            Y[k] += Vbar[i,k]*u[i] + Ubar[i,k]*v[i]
        end
    end
end

"""
add_inner_minus!(Ubar,u,x)

Inplace computation of `Ubar -= u .* x`.
"""
function add_inner_minus!(Ubar,u,x)
    for j=1:size(Ubar,2),i = 1:size(Ubar,1)
        Ubar[i,j] = Ubar[i,j] - u[i]*x[j]
    end
end
"""
add_inner_plus!(Vbar,v,x)

Inplace computation of `Vbar += v .* x`.
"""
function add_inner_plus!(Vbar,v,x)
    for j=1:size(Vbar,2),i = 1:size(Vbar,1)
        Vbar[i,j] = Vbar[i,j] + v[i]*x[j]
    end
end
"""
    add_outer_product!(wTw,w)

Inplace computation of `wTw = wTw + w*w'`
"""
function add_outer_product!(wTw,w)
    for j=1:size(wTw,2),i = 1:size(wTw,1)
        wTw[i,j] = wTw[i,j] + w[i]*w[j]
    end
end
"""
add_product!(Wbar,w,x)

Inplace computation of `Wbar  += w*x'`
"""
function add_product!(Wbar,w,x)
    for j=1:size(Wbar,2),i = 1:size(Wbar,1)
        Wbar[i,j] = Wbar[i,j] + w[i]*x[j]
    end
end


text:
import MultivariatePolynomials: pairzip, tuplezip

@testset "zip tests" begin
    @testset "pairzip" begin
        Mod.@polyvar x y z

        @test @inferred(pairzip((x, y), (1, 2))) == (x=>1, y=>2)
        @test_throws ArgumentError pairzip((x, y, z), (1, 2))
        @test_throws ArgumentError pairzip((x, y), (1, 2, 3))
        @test @inferred(pairzip((1, :x, r"x"), ("w", 1.0, +))) == (1=>"w", :x=>1.0, r"x"=>+)
        @test @inferred(pairzip((1, :x, r"x")=>("w", 1.0, +))) == (1=>"w", :x=>1.0, r"x"=>+)
    end

    @testset "tuplezip" begin
        Mod.@polyvar x y z

        @test @inferred(tuplezip((x, y), (1, 2))) == ((x, 1), (y, 2))
        @test_throws ArgumentError tuplezip((x, y, z), (1, 2))
        @test_throws ArgumentError tuplezip((x, y), (1, 2, 3))
        @test @inferred(tuplezip((1, :x, r"x"), ("w", 1.0, +))) == ((1, "w"), (:x, 1.0), (r"x", +))
    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule GeographicLib_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("GeographicLib")
JLLWrappers.@generate_main_file("GeographicLib", UUID("a660ed4b-f19b-5b77-bc0d-0bc97335ea30"))
end  # module GeographicLib_jll


text:
using DualNumbers, SpecialFunctions
using Test
using LinearAlgebra
import DualNumbers: value
import NaNMath

x = Dual(2, 1)
y = x^3

@test value(y) ≈ 2.0^3
@test epsilon(y) ≈ 3.0*2^2

y = x^3.0

@test value(y) ≈ 2.0^3
@test epsilon(y) ≈ 3.0*2^2

y = sin(x)+exp(x)
@test value(y) ≈ sin(2)+exp(2)
@test epsilon(y) ≈ cos(2)+exp(2)

@test x > 1
@test dual(1) < dual(2.0)
@test dual(1.0) < dual(2.0)
y = abs(-x)
@test value(y) ≈ 2.0
@test epsilon(y) ≈ 1.0

@test isequal(1.0,Dual(1.0))
@test Dual{Float32}(3) === Dual{Float32}(3.0f0, 0.0f0)

y = 1/x
@test value(y) ≈ 1/2
@test epsilon(y) ≈ -1/2^2

Q = [1.0 0.1; 0.1 1.0]
x = dual.([1.0,2.0])
x[1] = Dual(1.0,1.0)
y = (1/2)*dot(x,Q*x)
@test value(y) ≈ 2.7
@test epsilon(y) ≈ 1.2

function squareroot(x)
    it = x
    while abs(it*it - x) > 1e-13
        it = (it+x/it)/2
    end
    return it
end

@test epsilon(squareroot(Dual(10000.0,1.0))) ≈ 0.005

@test epsilon(exp(1)^Dual(1.0,1.0)) ≈ exp(1)
@test epsilon(NaNMath.pow(exp(1),Dual(1.0,1.0))) ≈ exp(1)
@test epsilon(NaNMath.sin(Dual(1.0,1.0))) ≈ cos(1)

@test Dual(1.0,3) == Dual(1.0,3.0)
x = Dual(1.0,1.0)
@test eps(x) == eps(1.0)
@test eps(Dual{Float64}) == eps(Float64)
@test one(x) == Dual(1.0,0.0)
@test one(Dual{Float64}) == Dual(1.0,0.0)
@test convert(Dual{Float64}, Inf) == convert(Float64, Inf)
@test isnan(convert(Dual{Float64}, NaN))

@test convert(Dual{Float64},Dual(1,2)) == Dual(1.0,2.0)
@test convert(Float64, Dual(10.0,0.0)) == 10.0
@test convert(Dual{Int}, Dual(10.0,0.0)) == Dual(10,0)

x = Dual(1.2,1.0)
@test floor(x) === 1.0
@test ceil(x)  === 2.0
@test trunc(x) === 1.0
@test round(x) === 1.0
@test floor(Int, x) === 1
@test ceil(Int, x)  === 2
@test trunc(Int, x) === 1
@test round(Int, x) === 1

# test Dual{Complex}

z = Dual(1.0+1.0im,1.0)
f = exp(z)
@test value(f) == exp(value(z))
@test epsilon(f) == epsilon(z)*exp(value(z))

g = sinpi(z)
@test value(g) == sinpi(value(z))
@test epsilon(g) == epsilon(z)*cospi(value(z))*π

h = z^4
@test value(h) == value(z)^4
@test epsilon(h) == 4epsilon(z)*value(z)^3

a = abs2(z)
@test value(a) == abs2(value(z))
@test epsilon(a) == conj(epsilon(z))*value(z)+conj(value(z))*epsilon(z)

l = log(z)
@test value(l) == log(value(z))
@test epsilon(l) == epsilon(z)/value(z)

s = sign(z)
@test value(s) == value(z)/abs(value(z))

a = angle(z)
@test value(a) == angle(value(z))

@test angle(Dual(0.0+im,0.0+im)) == π/2

# check bug in inv
@test inv(dual(1.0+1.0im,1.0)) == 1/dual(1.0+1.0im,1.0) == dual(1.0+1.0im,1.0)^(-1)

#
# Tests limit definition. Let z = a + b ɛ, where a and b ∈ C.
#
# The dual of |z| is lim_{h→0} (|a + bɛh| - |a|)/h
#
# and it depends on the direction (i.e. the complex value of epsilon(z)).
#

z = Dual(1.0+1.0im,1.0)
@test abs(z) ≡ sqrt(2) + 1/sqrt(2)*ɛ
z = Dual(1.0+1.0im,cis(π/4))
@test abs(z) ≡ sqrt(2) + 2/sqrt(2)^2*ɛ
z = Dual(1.0+1.0im,cis(π/2))
@test abs(z) ≡ sqrt(2) + 1/sqrt(2)*ɛ

# tests vectorized methods
const zv = dual.(collect(1.0:10.0), ones(10))

f = exp.(zv)
@test all(value.(f) .== exp.(value.(zv)))
@test all(epsilon.(f) .== epsilon.(zv) .* exp.(value.(zv)))

# tests norms and inequalities
@test norm(f,Inf) ≤ norm(f) ≤ norm(f,1)

# tests for constant ɛ
@test epsilon(1.0 + ɛ) == 1.0
@test epsilon(1.0 + 0.0ɛ) == 0.0
test(x, y) = x^2 + y
@test test(1.0 + ɛ, 1.0) == 2.0 + 2.0ɛ
@test test(1.0, 1.0 + ɛ) == 2.0 + 1.0ɛ

@test ɛ*im == Dual(Complex(false,false),Complex(false,true))

@test value(mod(Dual(15.23, 1), 10)) == 5.23
@test epsilon(mod(Dual(15.23, 1), 10)) == 1

@test epsilon(Dual(-2.0,1.0)^2.0) == -4
@test epsilon(Dual(-2.0,1.0)^Dual(2.0,0.0)) == -4

# test for flipsign
flipsign(Dual(1.0,1.0),2.0) == Dual(1.0,1.0)
flipsign(Dual(1.0,1.0),-2.0) == Dual(-1.0,-1.0)
flipsign(Dual(1.0,1.0),Dual(1.0,1.0)) == Dual(1.0,1.0)
flipsign(Dual(1.0,1.0),Dual(0.0,-1.0)) == Dual(-1.0,-1.0)
flipsign(-1.0,Dual(1.0,1.0)) == -1.0

# test SpecialFunctions
@test erf(dual(1.0,1.0)) == dual(erf(1.0), 2exp(-1.0^2)/sqrt(π))
@test gamma(dual(1.,1)) == dual(gamma(1.0),polygamma(0,1.0))

let x = exp10(Dual(2, 0.01))
    @test value(x) ≈ 100.0
    @test epsilon(x) ≈ log(10)
end

@test value(3) == 3
@test epsilon(44.0) ≈ 0.0


text:
using LinearAlgebra: mul!, det

export IESH

"""
    IESH{T} <: SurfaceHopping

Independent electron surface hopping.

# References
- [Shenvi, Roy, Tully, J. Chem. Phys. 130, 174107 (2009)](https://doi.org/10.1063/1.3125436)
- [Roy, Shenvi, Tully, J. Chem. Phys. 130, 174716 (2009)](https://doi.org/10.1063/1.3122989)

"""
struct IESH{T} <: SurfaceHopping
    hopping_probability::Matrix{T}
    state::Vector{Int}
    new_state::Vector{Int}
    proposed_state::Vector{Int}
    n_electrons::Int
    overlap::Matrix{Complex{T}}
    tmp::Vector{Complex{T}}
    function IESH{T}(states::Integer, n_electrons::Integer) where {T}

        hopping_probability = zeros(Int, states, n_electrons)
        state = zeros(Int, n_electrons)
        new_state = zero(state)
        proposed_state = zero(state)
        tmp = zeros(Complex{T}, states)
        overlap = zeros(Complex{T}, n_electrons, n_electrons)

        new{T}(hopping_probability, state, new_state, proposed_state,n_electrons, overlap, tmp)
    end
end

function NQCDynamics.Simulation{IESH}(atoms::Atoms{S,T}, model::Model; n_electrons, kwargs...) where {S,T}
    NQCDynamics.Simulation(atoms, model, IESH{T}(NQCModels.nstates(model), n_electrons); kwargs...)
end

function DynamicsMethods.DynamicsVariables(sim::Simulation{<:IESH}, v, r)
    ψ = zeros(NQCModels.nstates(sim.calculator.model), sim.method.n_electrons)

    for i=1:sim.method.n_electrons
        ψ[i,i] = 1
    end
    state = collect(1:sim.method.n_electrons)

    SurfaceHoppingVariables(ComponentVector(v=v, r=r, σreal=ψ, σimag=zero(ψ)), state)
end

"""
Set the acceleration due to the force from the currently occupied states.
See Eq. 12 of Shenvi, Tully JCP 2009 paper.
"""
function acceleration!(dv, v, r, sim::Simulation{<:IESH}, t, state)
    dv .= zero(eltype(dv))
    for I in eachindex(dv)
        for k in state
            # Contribution to the force from each occupied state `k`
            dv[I] -= sim.calculator.adiabatic_derivative[I][k, k]
        end
    end
    DynamicsUtils.divide_by_mass!(dv, sim.atoms.masses)
    return nothing
end

"""
Propagation of electronic wave function happens according to Eq. (14)
in the Shenvi, Tully paper (JCP 2009)

In IESH each electron is independent so we can loop through electrons and set the
derivative one at a time, in the standard way for FSSH.
"""
function DynamicsUtils.set_quantum_derivative!(dσ, v, σ, sim::Simulation{<:IESH})
    V = sim.calculator.eigen.values
    d = sim.calculator.nonadiabatic_coupling
    @views for i in axes(dσ, 2)
        set_single_electron_derivative!(dσ[:,i], σ[:,i], V, v, d, sim.method.tmp)
    end
end

"""
```math
\\frac{d c_k}{dt} = -iV_{kk}c_k/\\hbar - \\sum_j v d_{kj} c_j
```
"""
function set_single_electron_derivative!(dc, c, V, v, d, tmp)
    @. dc = -im*V * c
    for I in eachindex(v)
        mul!(tmp, d[I], c)
        lmul!(v[I], tmp)
        @. dc -= tmp
    end
    return nothing
end

"""
Hopping probability according to equation 21 in Shenvi, Roy, Tully 2009.
"""
function evaluate_hopping_probability!(sim::Simulation{<:IESH}, u, dt)

    ψ = DynamicsUtils.get_quantum_subsystem(u)
    v = DynamicsUtils.get_velocities(u)

    S = sim.method.overlap
    proposed_state = sim.method.proposed_state
    prob = sim.method.hopping_probability
    d = sim.calculator.nonadiabatic_coupling

    compute_overlap!(S, ψ, sim.method.state)
    det_current = det(S)
    Akk = abs2(det_current)

    fill!(prob, zero(eltype(prob)))
    for i in axes(prob, 2) # each electron
        for j in axes(prob, 1) # each basis state
            if j ∉ sim.method.state
                copyto!(proposed_state, sim.method.state)
                proposed_state[i] = j
                Akj = calculate_Akj(S, ψ, det_current, proposed_state)
                for I in eachindex(v)
                    prob[j,i] += 2 * v[I] * real(Akj/Akk) * d[I][sim.method.state[i], j] * dt
                end
            end
        end
    end

    clamp!(prob, 0, 1) # Restrict probabilities between 0 and 1

    return nothing
end

"Equation 17 in Shenvi, Roy, Tully 2009. Uses equations 19 and 20."
function calculate_Akj(S, ψ, detS, new_state)

    compute_overlap!(S, ψ, new_state)
    det_new = det(S)

    return detS*conj(det_new)
end

"Equation 20 in Shenvi, Roy, Tully 2009."
function compute_overlap!(S::Matrix, ψ, state)
    for i in axes(S, 2) # each electron
        for j in axes(S, 1) # each electron
            S[j,i] = ψ[state[j],i]
        end
    end
    return nothing
end

function select_new_state(sim::AbstractSimulation{<:IESH}, u)::Vector{Int}

    prob = sim.method.hopping_probability

    random = rand()
    cumulative = 0
    for i in axes(prob, 2) # each electron
        for j in axes(prob, 1) # each basis state
            if j ∉ sim.method.state
                cumulative += prob[j,i]
                if random < cumulative
                    copyto!(sim.method.proposed_state, sim.method.state)
                    sim.method.proposed_state[i] = j
                    return sim.method.proposed_state
                end
            end
        end
    end

    return sim.method.state
end

function Estimators.diabatic_population(sim::Simulation{<:IESH}, u)
    Calculators.evaluate_potential!(sim.calculator, DynamicsUtils.get_positions(u))
    Calculators.eigen!(sim.calculator)
    U = sim.calculator.eigen.vectors

    ψ = DynamicsUtils.get_quantum_subsystem(u)
    diabatic_ψ = zero(ψ)
    @views for i in axes(ψ, 2) # each electron
        diabatic_ψ[:,i] .= U*ψ[:,i]
    end
    diabatic_population = sum(abs2.(diabatic_ψ), dims=2)

    return diabatic_population
end

function Estimators.adiabatic_population(sim::Simulation{<:IESH}, u)
    population = zeros(NQCModels.nstates(sim.calculator.model))
    population[u.state] .= 1
    return population
end

unpack_states(sim::Simulation{<:IESH}) = symdiff(sim.method.new_state, sim.method.state)

# function impurity_summary(model::DiabaticModel, R::AbstractMatrix, state::AbstractArray, σ::AbstractArray)
#     """Calculate impurity population according to MiaoSubotnik_JChemPhys_150_041711_2019"""

#     eig_vec = zeros(model.n_states,model.n_states)

#     eival = zeros(model.n_states)
#     σad = zeros(Complex, model.n_states, model.n_states)
#     tmp = 0
#     σdia = zeros(Complex, model.n_states, model.n_states)
#     eig_array = zeros(4+model.n_states)
#     V = Hermitian(zeros(model.n_states,model.n_states))
#     dvect = zeros(model.n_states)
#     dvect[2] = 1

#     # get potential
#     potential!(model,V,R)
#     eig_vec .= eigvecs(V)
#     ieig = inv(eig_vec)
#     # Get density matrix matrix
#     for i = 1:length(state)
#         eig_array[4 + i] = norm(σ[i,:])
#         for j = 1:length(state)
#             for k in axes(σ,2)
#                 σad[i,j] = σad[i,j] + σ[i,k]*σ[j,k]
#             end
#         end
#     end

#     #Turn into diabatic matrix for impurity population
#     σdia .= eig_vec *σad * ieig
#     eig_array[4] = real(σdia[2,2])^2 + imag(σdia[2,2])^2

#     # eig_array[4] = 0
#     # for k in axes(σ,2)
#     #     for i = 1:length(state)
#     #         σad[i,i] = σ[i,k]
#     #     end
#     #     σdia .= 0
#     #     σdia .= eig_vec *σad * ieig
#     #     eig_array[4] = eig_array[4] + real(σdia[2,2])^2 + imag(σdia[2,2])^2
#     # end

#     # Get the eigenvectors and values
#     eival .= eigvals(V)

#     # save position
#     eig_array[1] = R[1]
#     for i = 1:length(state)
#         # Energy
#         eig_array[2] = eig_array[2] + eival[state[i]]
#         # Hopping prob. by hopping array
#         eig_array[3] = eig_array[3] + state[i]

#     end

#     # # over electrons. This is adiabatic.
#     # for i in axes(σ, 2)
#     #     eig_array[4] = eig_array[4] + real(σ[2,i])^2 + imag(σ[2,i])^2
#     # end

#     # Export an array of eigenvalues with last two elements being hopping prob
#     eig_array = eig_array
# end


text:
"""
    UniversalFunctions

Universal functions for the surface fluxes module.
"""
module UniversalFunctions

export Pr_0_Businger,
    a_m_Businger,
    a_h_Businger,
    Pr_0_Gryanik,
    a_m_Gryanik,
    a_h_Gryanik,
    b_m_Gryanik,
    b_h_Gryanik,
    Pr_0_Grachev,
    a_m_Grachev,
    a_h_Grachev,
    b_m_Grachev,
    b_h_Grachev,
    c_h_Grachev

#####
##### Businger model
#####

# Reference:
#
#    Nishizawa, S., and Y. Kitamura. "A Surface Flux
#    Scheme Based on the Monin‐Obukhov Similarity for
#    Finite Volume Models." Journal of Advances in
#    Modeling Earth Systems 10.12 (2018): 3159-3175.
#    Appendix A.

# Original work
#
#    Businger, Joost A., et al. "Flux-profile
#    relationships in the atmospheric surface layer."
#    Journal of the atmospheric Sciences 28.2 (1971):
#    181-189.

""" Prandtl number at neutral stratification """
function Pr_0_Businger end
""" coefficient for momentum """
function a_m_Businger end
""" coefficient for heat """
function a_h_Businger end

#####
##### Gryanik model
#####

# Reference
#
#    Gryanik, Vladimir M., et al. "New modified and extended
#    stability functions for the stable boundary layer based
#    on SHEBA and parametrizations of bulk transfer coefficients
#    for climate models." Journal of the Atmospheric Sciences
#    (2020).

""" neutral-limit turbulent Prandtl number """
function Pr_0_Gryanik end
""" empirical coefficient (a) for momentum """
function a_m_Gryanik end
""" empirical coefficient (a) for heat """
function a_h_Gryanik end
""" empirical coefficient (b) for momentum """
function b_m_Gryanik end
""" empirical coefficient (b) for heat """
function b_h_Gryanik end

#####
##### Grachev model
#####

# Reference:
#
#    Grachev, A. A., E. L. Andreas, C. W. Fairall, P. S. Guest, and
#    P. O. G. Persson, 2007a: SHEBA flux–profile relationships in
#    the stable atmospheric boundary layer. Bound.-Layer Meteor.,
#    124, 315–333, https://doi.org/10.1007/s10546-007-9177-6.

""" Prandtl number at neutral stratification """
function Pr_0_Grachev end
""" empirical coefficient (a) for momentum """
function a_m_Grachev end
""" empirical coefficient (a) for heat """
function a_h_Grachev end
""" empirical coefficient (b) for momentum """
function b_m_Grachev end
""" empirical coefficient (b) for heat """
function b_h_Grachev end
""" empirical coefficient (c) for heat """
function c_h_Grachev end

end # module UniversalFunctions

max_stars_repo_name: climate-machine/CLIMAParameters

max_issues_repo_name: climate-machine/CLIMAParameters

max_forks_repo_name: climate-machine/CLIMAParameters