text:
#Combined Gradient Descent with RMSprop element
#http://www.cs.toronto.edu/~tijmen/csc321/slides/lectureslideslec6.pdf

struct RMSprop{Tf,Tg,Th,Ti, IL,L} <: FirstOrderOptimizer
    mu::Tf
    alpha::Tg
	beta::Th
	epsilon::Ti
    alphaguess!::IL
    linesearch!::L
    manifold::Manifold
end

Base.summary(::RMSprop) = "Root Mean Square propagation"

"""
# Root Mean Square propagation (RMSprop)
## Constructor
```julia
RMSprop(; mu::Real = 0.9,
 alpha::Real = 0.00000001,
 beta::Real = 0.99,
 epsilon::Real = 0.000001,
 alphaguess = LineSearches.InitialPrevious(),
 linesearch = LineSearches.HagerZhang(),
```

## Description
RMSProp (for Root Mean Square Propagation) is a method in which the learning rate(here alpha) is adapted for
each of the parameters. The idea is to divide the learning rate for a weight by a running average of the
magnitudes of recent gradients for that function.

In this case, RMSprop is combined with simple gradient descent algorithm, which uses a line search step
to compute the next step.

## Formula
y = y_previous * beta + ( 1 - beta ) * g ^ 2

x = x_previous + mu * a * g + alpha * g / (sqrt(y) + epsilon)

parameters:

a= alpha calculated of line search

g= gradient of recent iteration

## References
 - RMSprop is an unpublished method proposed by Geoff Hinton in
  http://www.cs.toronto.edu/~tijmen/csc321/slides/lecture_slides_lec6.pdf
"""
function RMSprop(; mu::Real = 0.9,
								 alpha::Real = 0.00000001, # TODO: investigate good defaults
								 beta::Real = 0.99, # TODO: investigate good defaults
								 epsilon::Real = 0.000001, # TODO: investigate good defaults
                                 alphaguess = LineSearches.InitialPrevious(), # TODO: investigate good defaults
                                 linesearch = LineSearches.HagerZhang(),        # TODO: investigate good defaults
                                 manifold::Manifold=Flat())
    RMSprop(mu,alpha, beta,epsilon, alphaguess, linesearch, manifold)
end

mutable struct RMSpropState{Tx, T} <: AbstractOptimizerState
    x::Tx
    x_previous::Tx
	y::Tx
    f_x_previous::T
    s::Tx
    @add_linesearch_fields()
end

function initial_state(method::RMSprop, options, d, initial_x)
    T = eltype(initial_x)
    initial_x = copy(initial_x)
	initial_y = copy(initial_x)
	fill!(initial_y,0.0)
    retract!(method.manifold, initial_x)

    value_gradient!!(d, initial_x)

    project_tangent!(method.manifold, gradient(d), initial_x)

    RMSpropState(initial_x, # Maintain current state in state.x
                                 copy(initial_x), # Maintain previous state in state.x_previous
								 initial_y,
                                 real(T)(NaN), # Store previous f in state.f_x_previous
                                 similar(initial_x), # Maintain current search direction in state.s
                                 @initial_linesearch()...)
end

function update_state!(d, state::RMSpropState, method::RMSprop)
    project_tangent!(method.manifold, gradient(d), state.x)
    # Search direction is always the negative gradient
    state.s .= .-gradient(d)

    # Update moving average of squared gradients
    state.y .= state.y.*method.beta .+(1-method.beta).*NaNMath.pow.(state.s,2)

    # Determine the distance of movement along the search line
    lssuccess = perform_linesearch!(state, method, ManifoldObjective(method.manifold, d))

    # Update current position with alpha step + element calculated of RMSprop algorithm
    state.x .+= method.mu.*state.alpha.*state.s .+method.alpha.*state.s./(Base.sqrt.(state.y).+method.epsilon)
    retract!(method.manifold, state.x)
    lssuccess == false # break on linesearch error
end

function trace!(tr, d, state, iteration, method::RMSprop, options, curr_time=time())
  common_trace!(tr, d, state, iteration, method, options, curr_time)
end

function default_options(method::RMSprop)
    Dict(:allow_f_increases => true)
end


text:
"""
    AbstractSource

The abstract supertype for custom data sources. Implement a subtype of `AbstractSource` if
`Source{S}` is not sufficient for your source requirements (e.g. your source has additional
information besides the path, such as encoding parameters, compression indicator, etc, that
needs to be associated with each instance).

# Extended help

All subtypes of `AbstractSource` **must**:

- have a `path` field or extend the `sourcepath` function.
- accept a single string, the absolute path of the source file, in the constructor.

Minimum method definitions necessary for complete functionality of a subtype of
`AbstractSource` (i.e. `Source{S}` or otherwise):
- [`readsource`](@ref)
- [`sourcepath`](@ref) (if defining a new subtype of `AbstractSource`)

Additional methods to improve user experience and/or additional functionality
- [`readsegment`](@ref)
- [`generatesource`](@ref) (if enabling `requiresource!` generation)
- [`dependencies`](@ref) (if defining a `generatesource` method)
- [`srcext`](@ref)
- [`srcname_default`](@ref)

## Example

```julia
struct MySource <: AbstractSource
    path::String
end

DatasetManager.readsource(src::MySource) = customreadfunc(sourcepath(src))
DatasetManager.readsegment(seg::Segment{MySource}) = readsource(seg.source)[start:finish]
```
"""
abstract type AbstractSource end

(::Type{S})() where S <: AbstractSource = S(tempname()*srcext(S))

"""
    Source{S}(path)

A basic source, where `S` can be any singleton or existing type
"""
struct Source{S} <: AbstractSource
    path::String
end

Source{S}() where S = Source{S}(tempname()*srcext(Source{S}))

"""
    sourcepath(src) -> String

Return the absolute path to the `src`.
"""
sourcepath(src::AbstractSource) = src.path::String

"""
    srcext(src::Union{S,Type{S}})

Actual file extension or default file extension for a `src` or `src` type
"""
srcext(src::AbstractSource) = splitext(sourcepath(src))[2]
srcext(::Type{<:AbstractSource}) = ""

"""
    readsource(src::S; kwargs...) where {S <: AbstractSource}

Read the source data from file.
"""
function readsource end

struct UnknownDeps end

"""
    dependencies(src::Union{S,Type{S}})

Get the sources that `src` depends on to be generated by `generatesource`.
"""
dependencies(::S) where S <: AbstractSource = dependencies(S)
dependencies(::Type{<:AbstractSource}) = UnknownDeps()

"Default name for a source in `Trial.sources`"
srcname_default(::S) where S <: AbstractSource = srcname_default(S)
srcname_default(s) = string(s)

struct MissingSourceError <: Exception
    msg::String
end

function Base.showerror(io::IO, e::MissingSourceError)
    print(io, "MissingSourceError: ")
    print(io, e.msg)
end

"""
    requiresource!(trial, src::Union{S,Type{S}}; force=false, deps, kwargs...) -> nothing
    requiresource!(trial, name::String => src::Union{S,Type{S}}; force=false, deps, kwargs...)

Require a source `src` in `trial`, generate `src` if not present, or throw an error if `src`
is not present and cannot be generated. Returns `nothing`. Keyword argument `deps` defaults to
`dependencies(src)`, but can be manually specified as a `name => src::AbstractSource` pair
per dependency if there exist multiple sources that would match a dependent source type.
`kwargs...` are passed on to the `generatesource` function for `src`.

# Keyword arguments
- `force=false`: Force generating `src`, even if it already exists
- `deps=dependencies(src)`: Manually specify particular dependencies

# Examples
```julia
julia> dependencies(Source{Events})
(Source{C3DFile},)

julia> requiresource!(trial, Source{Events})

julia> requiresource!(trial, Source{Events}; force=true, deps=("mainc3d" => Source{C3DFile}))

```
"""
function requiresource!(trial, src::AbstractSource, parent=nothing; kwargs...)
    requiresource!(trial, srcname_default(src) => src, parent; kwargs...)
end

function requiresource!(trial, src::Type{<:AbstractSource}, parent=nothing; kwargs...)
    requiresource!(trial, srcname_default(src) => src(), parent; kwargs...)
end

function requiresource!(trial, namesrc::Pair, parent=nothing;
    force=false, deps=dependencies(namesrc.second), kwargs...
)
    name, src = namesrc
    if !force
        if hassource(trial, name) || hassource(trial, src) ||
            (parent !== nothing && hassource(trial, typeof(src)))
            return nothing
        elseif isfile(sourcepath(src))
            sources(trial)[name] = src
            return nothing
        end
    end

    if deps isa UnknownDeps
        throw(MissingSourceError("unable to generate missing source $src for $trial"))
    else
        foreach(deps) do reqsrc
            requiresource!(trial, reqsrc, src; force=false)
        end
    end
    _src = generatesource(trial, src, deps; kwargs...)
    isfile(sourcepath(_src)) ||
        throw(MissingSourceError("failed to generate source $name => $_src for $trial"))

    sources(trial)[name] = _src

    return nothing
end

"""
    generatesource(trial, src::Union{S,Type{S}}, deps; kwargs...) where S <: AbstractSource -> newsrc::typeof(src)

Generate source `src` using dependent sources `deps` from `trial`. Returns a source of the same
type as `src`, but is not required to be exactly equal to `src` (i.e. a different
`sourcepath(newsrc)` is acceptable).
"""
function generatesource end

generatesource(trial, ::Type{S}, deps; kwargs...) where S <: AbstractSource =
    generatesource(trial, S(), deps; kwargs...)

ext: jl


text:
"""

$(TYPEDEF)

Auxiliary structure used to update counters.

$(TYPEDFIELDS)

"""
mutable struct DminMol
    d::Float64
    jmol::Int
    iat::Int
    jat::Int
end


text:
function uhlig04()
#
#  test polynomial used by F. Uhlig
#
    r = [(3/11)*ones(12); 11/3; 11/3; (2*im/7)*ones(4); (2.5+im/4)*ones(2); 1/8]
    p = poly(r);
    z = [3/11 12; 11/3 2; 2*im/7 4; 2.5+im/4 2; 1/8 1];

    p, PolyZeros(z)
end


text:
@testset "orthogonal_matrix" begin
    U = rand(5, 3)
    S = Diagonal(rand(3))

    @test_throws ArgumentError Orthogonal(U, S)

    U, S, _ = svd(rand(5,3))
    H = Orthogonal(U, Diagonal(S))

    @test size(H) == size(U)
    @test H.S == Diagonal(S)
    @test H.U == U
    @test H == U * sqrt(Diagonal(S))
end
@info "Ran orthogonal_matrix tests."


text:
include("pristine_graphene.jl")

"""
    ImpurityState(ϵ::Float64, coupling::Vector{Tuple{Float64,GrapheneCoord}})

An impurity state of energy `ϵ` (in eV) coupled to the graphene system. The
tuples in the `coupling` field contain all the coupling energies (in eV) and the
 corresponding [`GrapheneCoord`](@ref)'s.
"""
struct ImpurityState
    ϵ::Float64
    coupling::Vector{Tuple{Float64,GrapheneCoord}}
end

"""
    GrapheneSystem(
        μ::Float64,
        T::Float64,
        Δ::Array{Float64,2},
        V::Array{Float64,2},
        scattering_atoms::Vector{GrapheneCoord},
        imps::Vector{Float64},
    )

A structure describing the perturbed graphene system.

See also [`mkGrapheneSystem`](@ref).
"""
struct GrapheneSystem
    μ::Float64                              # Chemical potential
    T::Float64                              # Temperature
    Δ::Array{Float64,2}                     # Δ matrix
    V::Array{Float64,2}                     # V Matrix
    scattering_atoms::Vector{GrapheneCoord} # List of all perturbed atoms
    imps::Vector{Float64}                   # Impurity energies
end

"""
    mkGrapheneSystem(
        μ::Float64,
        T::Float64,
        imps::Vector{ImpurityState},
        pert::Vector{Tuple{GrapheneCoord,GrapheneCoord,Float64}},
    )

Construct [`GrapheneSystem`](@ref).

# Arguments
* `μ`: chemical potential
* `T`: temperature
* `imps`: a list of [`ImpurityState`]'s
* `pert`: a list of 3-tuples describing the coupling between carbon atoms

When supplying `pert`, ensure that each coordinate pair appears only once as
repeated pairs with different couplings will cause eariler values to be
overwritten. The order of the coordinates does not matter. A tuple with the same
coordinate given twice creates and on-site potential for the coordinate.

The function constructs a [`GrapheneSystem`](@ref) with the `Δ` and `V` matrix
computed from `imps` and `pert`. In addition, a list of all
[`GrapheneCoord`](@ref) that are perturbed (`scattering_atoms` field in
[`GrapheneCoord`](@ref)) and a list of
impurity energies (`imps` field in [`GrapheneCoord`](@ref)) are included.
"""
function mkGrapheneSystem(
    μ::Float64,
    T::Float64,
    imps::Vector{ImpurityState},
    pert::Vector{Tuple{GrapheneCoord,GrapheneCoord,Float64}},
)
    # Create a coupling dictionary
    pert_dict =
        reduce(
            vcat,
            map(x -> [((x[1], x[2]), x[3]), ((x[2], x[1]), x[3])], pert),
        ) |> Dict
    # Get the coordinates of all the directly-perturbed atoms
    perturbed_atoms = map(x -> x[1], pert_dict |> keys |> collect) |> unique
    # Get the coordinates of atoms coupled to impurities
    coupled_atoms = map(
        y -> y[2],
        reduce(
            vcat,
            map(x -> x.coupling, imps),
            init = Vector{Tuple{Float64,GrapheneCoord}}[],
        ),
    )
    # Combine the two sets of coordinates and keep only unique entries
    all_atoms = vcat(perturbed_atoms, coupled_atoms) |> unique |> sort
    # Assemble Δ
    all_atoms_M = repeat(all_atoms, 1, length(all_atoms))
    all_atoms_M_T = permutedims(all_atoms_M)
    Δ = map((x, y) -> get!(pert_dict, (x, y), 0.0), all_atoms_M, all_atoms_M_T)
    # Assemble V
    V_array = map(
        imp -> map(
            atom -> sum(
                map(c -> ((atom == c[2]) * c[1]), imp.coupling),
            )::Float64,
            all_atoms,
        ),
        imps,
    )
    if !isempty(V_array)
        V = reduce(hcat, V_array)
    else
        V = Array{GrapheneCoord}(undef, 0, 0)
    end
    imp_energies = map(x -> x.ϵ, imps)
    return GrapheneSystem(μ, T, Δ, V, all_atoms, imp_energies)
end


text:
using Test
using ParallelMapReduce
using Distributed

@testset "Pmapreduce" begin

    if nprocs() == 1
        addprocs(2)
    elseif nworkers() > 2
        error("Too many processes")
    end

    @everywhere using Pkg
    @everywhere pkg"activate ."
    @everywhere pkg"precompile"
    @everywhere using ParallelMapReduce
    function test()
        @everywhere function f(n)
            sleep(n)
            n
        end

        args = [5, ones(Int, 9)...]

        val_even, t_even  = @timed pmapreduce(f, +, args)
        val_red_local, t_red_local  = @timed pmapreduce(
            f, +, args; algorithm = :reduction_local
        )
        val_red_master, t_red_master  = @timed pmapreduce(
            f, +, args; algorithm = :reduction_master
        )

        @test val_even ==  val_red_master == val_red_local
        @test t_red_local < t_even
        @test t_red_master < t_even
    end

    test()

end


text:
const LOG2 = log(2)
@inline vmullog2(x::T) where {T<:Number} = vmul(T(LOG2), x)
@inline vmullog2(v::V) where {W,T,V<:AbstractSIMDVector{W,T}} = vmul(vbroadcast(V, T(LOG2)), v)
@inline vmullog2add!(y, x::T) where {T<:Number} = vfmadd(T(LOG2), x, y)
@inline vmullog2add!(y, v::V) where {W,T,V<:AbstractSIMDVector{W,T}} = vfmadd(vbroadcast(V, T(LOG2)), v, y)

const LOG10 = log(10)
@inline vmullog10(x::T) where {T<:Number} = vmul(T(LOG10), x)
@inline vmullog10(v::V) where {W,T,V<:AbstractSIMDVector{W,T}} = vmul(vbroadcast(V, T(LOG10)), v)
@inline vmullog10add!(y, x::T) where {T<:Number} = vfmadd(T(LOG10), x, y)
@inline vmullog10add!(y, v::V) where {W,T,V<:AbstractSIMDVector{W,T}} = vfmadd(vbroadcast(V, T(LOG10)), v, y)

const INVLOG2 = 1 / log(2)
@inline vdivlog2(x::T) where {T<:Number} = vmul(T(INVLOG2), x)
@inline vdivlog2(v::V) where {W,T,V<:AbstractSIMDVector{W,T}} = vmul(vbroadcast(V, T(INVLOG2)), v)
@inline vdivlog2add!(y, x::T) where {T<:Number} = vfmadd(T(INVLOG2), x, y)
@inline vdivlog2add!(y, v::V) where {W,T,V<:AbstractSIMDVector{W,T}} = vfmadd(vbroadcast(V, T(INVLOG2)), v, y)

const INVLOG10 = 1 / log(10)
@inline vdivlog10(x::T) where {T<:Number} = vmul(T(INVLOG10), x)
@inline vdivlog10(v::V) where {W,T,V<:AbstractSIMDVector{W,T}} = vmul(vbroadcast(V, T(INVLOG10)), v)
@inline vdivlog10add!(y, x::T) where {T<:Number} = vfmadd(T(INVLOG10), x, y)
@inline vdivlog10add!(y, v::V) where {W,T,V<:AbstractSIMDVector{W,T}} = vfmadd(vbroadcast(V, T(INVLOG10)), v, y)

@inline vfmaddaddone(x::T) where {T} = vfmadd(x, x, VectorizationBase.vone(promote_type(T)))
@inline vfmaddaddone(x::T1, y::T2) where {T1,T2} = vfmadd(x, y, VectorizationBase.vone(promote_type(T1,T2)))

@inline vadd!(y, x) = vadd(y, x)
@inline vsub!(y, x) = vsub(y, x)
@inline vmul!(y, x) = vmul(y, x)
@inline vfdiv!(y, x) = vfdiv(y, x)
@inline vfmadd!(c, a, b) = vfmadd(a, b, c)
@inline vfnmadd!(c, a, b) = vfnmadd(a, b, c)
@inline vfmsub!(c, a, b) = vfmsub(a, b, c)
@inline vfnmsub!(c, a, b) = vfnmsub(a, b, c)


text:
# preferred units are THz, ps, and m
freq = collect(linspace(150, 250, 250))THz
neff_data = collect(linspace(3.0, 1.0, 250))
neff = OpticalAttr(freq, neff_data, "n_eff")
corefrac = OpticalAttr(1.0, "Core Fraction")
Aeff = OpticalAttr(1e-12m^2, "Aeff (m²)")
loss = OpticalAttr(0.0m^-1, "Loss (1/m)")
coupling = OpticalAttr(0.6, "Coupling")

mode_TM = Mode(neff, Aeff, loss, corefrac, coupling, :TM, true, true)
mode_TE = Mode(neff, Aeff, loss, corefrac, coupling, :TE, true, true)

# 10 mm long waveguide in default (100) orientation
wg_simple = Waveguide(Silicon)
add_straight!(wg_simple, 1e-3m)
add_turn!(wg_simple, 90°, 1e-3m)
add_mode!(wg_simple, mode_TM)
add_mode!(wg_simple, mode_TE)
add_interaction!(wg_simple, mode_TM, mode_TE)
#add_medium!(wg_simple, Diamond)
#
# Note: pulse width is in 1/THz = fs
source_simple = PulsedLaser(200.0THz, 1.0W, 0.20ps)
model_NLSE = Model(source_simple, wg_simple)
# Laser answers: δ₀ can change during experiment, ω₀, and Ein? All can be
# parameters of cw laser. Solver ?s: How many dispersive orders
# should we use (keyword, default=3)? Initial conditions for NLSE? How do we
# determine initial conditions otherwise?
#solve(model_simple, DynamicLL()) (DynamicLL, DynamicIkeda, SSLL, SSIkeda, NLSE)

res_simple = CircularResonator(300e-6, SiO2)
add_mode!(res_simple, mode_TE)
source_CW = CWLaser(200THz, 0.05, 1.0)
model_LL = Model(source_CW, [res_simple, res_simple])

sol_NLSE = solve(model_NLSE, StochasticNLSE(); time_pts=2^5)


text:
push!(LOAD_PATH,"../docs/src/")
push!(LOAD_PATH,"./src")

using Documenter, AgentBasedModels

makedocs(sitename="AgentBasedModel.jl",
pages = [
    "Home" => "index.md",
    "Usage.md",
    "API.md",
    "APIdevelopers.md"
],
format = Documenter.HTML(prettyurls = false)
)

deploydocs(
    repo = "github.com/dsb-lab/AgentBasedModels.jl",
)