text:
@testset "8.7 Zeta function" begin
    (a, b, s, x, ) = @variables a b s x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integration*Problems*Involving*Zeta*Functions=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*Zeta(2, a+b*x)=#

    @test_int [x^2*Zeta(2, a + b*x), x, 4, -((2*x*LogGamma(a + b*x))/b^2) + (2*PolyGamma(-2, a + b*x))/b^3 + (x^2*PolyGamma(0, a + b*x))/b]
    @test_int [x^1*Zeta(2, a + b*x), x, 3, -(LogGamma(a + b*x)/b^2) + (x*PolyGamma(0, a + b*x))/b]
    @test_int [x^0*Zeta(2, a + b*x), x, 2, PolyGamma(0, a + b*x)/b]
    @test_int [Zeta(2, a + b*x)/x^1, x, 1, Unintegrable(PolyGamma(1, a + b*x)/x, x)]
    @test_int [Zeta(2, a + b*x)/x^2, x, 2, b*Unintegrable(PolyGamma(2, a + b*x)/x, x) - PolyGamma(1, a + b*x)/x]
    @test_int [Zeta(2, a + b*x)/x^3, x, 3, (1/2)*b^2*Unintegrable(PolyGamma(3, a + b*x)/x, x) - PolyGamma(1, a + b*x)/(2*x^2) - (b*PolyGamma(2, a + b*x))/(2*x)]

    @test_int [Zeta(2, a + b*x)/x^2 - b*(PolyGamma(2, a + b*x)/x), x, 3, -(PolyGamma(1, a + b*x)/x)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*Zeta(s, a+b*x)=#

    @test_int [x^2*Zeta(s, a + b*x), x, 3, If(13>=8, (2*Zeta(-3 + s, a + b*x))/(b^3*(1 - s)*(2 - s)*(3 - s)) - (2*x*Zeta(-2 + s, a + b*x))/(b^2*(1 - s)*(2 - s)) + (x^2*Zeta(-1 + s, a + b*x))/(b*(1 - s)), (2*Zeta(-3 + s, a + b*x))/(b^3*(1 - s)*(6 - 5*s + s^2)) - (2*x*Zeta(-2 + s, a + b*x))/(b^2*(1 - s)*(2 - s)) + (x^2*Zeta(-1 + s, a + b*x))/(b*(1 - s)))]
    @test_int [x^1*Zeta(s, a + b*x), x, 2, -(Zeta(-2 + s, a + b*x)/(b^2*(1 - s)*(2 - s))) + (x*Zeta(-1 + s, a + b*x))/(b*(1 - s))]
    @test_int [x^0*Zeta(s, a + b*x), x, 1, Zeta(-1 + s, a + b*x)/(b*(1 - s))]
    @test_int [Zeta(s, a + b*x)/x^1, x, 0, CannotIntegrate(Zeta(s, a + b*x)/x, x)]
    @test_int [Zeta(s, a + b*x)/x^2, x, 1, (-b)*s*CannotIntegrate(Zeta(1 + s, a + b*x)/x, x) - Zeta(s, a + b*x)/x]
    @test_int [Zeta(s, a + b*x)/x^3, x, 2, (1/2)*b^2*s*(1 + s)*CannotIntegrate(Zeta(2 + s, a + b*x)/x, x) - Zeta(s, a + b*x)/(2*x^2) + (b*s*Zeta(1 + s, a + b*x))/(2*x)]

    @test_int [Zeta(s, a + b*x)/x^2 + b*s*(Zeta(1 + s, a + b*x)/x), x, 2, -(Zeta(s, a + b*x)/x)]
    end

max_stars_repo_path: test/8 Special functions/8.7 Zeta function.jl


text:
# Responses
printstyled("  InstrumentResponse subtypes\n", color=:light_green)
redirect_stdout(out) do
  v = 1.0 + 1.0*im
  X = rand(12,3)
  Y = rand(12,3)
  R = GenResp("", X, Y)
  @test R == GenResp(complex.(X,Y))
  @test hash(R) == hash(GenResp(complex.(X,Y)))
  @test sizeof(R) > sizeof(X)+sizeof(Y)
  @test R[10] == getindex(R.resp, 10) == getindex(R, 10)
  @test R[11,2] == getindex(R.resp, 11, 2) == getindex(R, 11, 2)
  R[3] = v
  R[4,2] = 1.0
  @test getindex(R.resp, 3) == v == getindex(R, 3, 1)
  @test real(getindex(R.resp, 4, 2)) == 1.0
  show(stdout, R)
  repr(R, context=:compact=>true)

  for T in (CoeffResp, GenResp, MultiStageResp, PZResp, PZResp64)
    R = T()
    @test isempty(R) == true
    show(stdout, R)
    repr(R, context=:compact=>true)
    repr(R, context=:compact=>false)

    R2 = T()
    @test R == R2
    @test hash(R) == hash(R2)
    @test code2resptyp(resptyp2code(R)) == T
    @test sizeof(R) == sizeof(R2)
  end

  nr = 255
  R = MultiStageResp(nr)
  for f in (:fs, :gain, :fg, :delay, :corr, :fac, :os)
    setfield!(R, f, rand(eltype(getfield(R, f)), nr))
  end
  i = Array{String,1}(undef,nr)
  o = Array{String,1}(undef,nr)
  o[1] = "m/s"
  i[1] = "{counts}"
  for j = 2:nr
    i[j] = randstring(2^rand(2:6))
    o[j] = i[j-1]
  end
  R.i = i
  R.o = o
  R.stage[1] = RandSeis.randResp()
  for i = 2:nr
    R.stage[i] = CoeffResp(b = rand(rand(1:1200)))
  end
  show(R)
end

printstyled("    resp codes\n", color=:light_green)
codes = (0x00, 0x01, 0x02, 0x03, 0x04, 0xff)
types = (GenResp, PZResp, PZResp64, CoeffResp, MultiStageResp, Nothing)
for c in codes
  T = code2resptyp(c)()
  @test c == resptyp2code(T)
end
for (i, T) in enumerate(types)
  resp = T()
  @test resptyp2code(resp) == codes[i]
end

ext: jl


text:
using ExoplanetsSysSim

function set_test_param(sim_param_closure::SimParam)
    p_lim_arr_num = [0.5, 1.25, 2.5, 5., 10., 20., 40., 80., 160., 320.]
    r_lim_arr_num = [0.5, 0.75, 1., 1.25, 1.5, 1.75, 2., 2.5, 3., 4., 6., 8., 12., 16.]
    rate_tab_init = [0.004 0.031 0.124 0.075 0.162 0.175 0.328 0.503 0.591 ;
                     0.004 0.012 0.082 0.165 0.234 0.407 0.321 1.223 2.381 ;
                     0.001 0.007 0.057 0.103 0.186 0.234 0.764 0.700 1.915 ;
                     0.009 0.012 0.127 0.229 0.483 0.822 0.910 1.463 1.795 ;
                     0.004 0.029 0.136 0.599 1.409 1.889 2.279 3.366 4.987 ;
                     0.007 0.042 0.218 0.932 2.270 2.832 3.500 3.135 3.166 ;
                     0.010 0.059 0.289 1.148 2.527 3.658 3.550 3.616 6.367 ;
                     0.036 0.057 0.380 0.796 1.552 1.181 1.959 2.420 3.505 ;
                     0.047 0.126 0.738 1.089 1.909 1.924 2.254 2.675 2.199 ;
                     0.053 0.126 0.630 1.608 2.405 2.411 2.272 4.370 7.607 ;
                     0.088 0.180 0.524 1.523 2.801 3.055 4.071 3.216 8.540 ;
                     0.100 0.295 1.022 3.156 3.559 7.929 0.0852 2.722 10.4 ;
                     0.039 0.404 2.615 5.095 5.378 1.969 7.989  0.0   0.0  ] * 0.01
    #rate_tab_init = reshape(rate_tab_init, (length(p_lim_arr_num)-1)*(length(r_lim_arr_num)-1))
    add_param_fixed(sim_param_closure, "p_lim_arr", p_lim_arr_num)
    add_param_fixed(sim_param_closure, "r_lim_arr", r_lim_arr_num*ExoplanetsSysSim.earth_radius)
    add_param_fixed(sim_param_closure, "obs_par", rate_tab_init)

  add_param_active(sim_param_closure,"sigma_hk",0.03)
  add_param_active(sim_param_closure,"sigma_incl",1.0)   # degrees; 0 = coplanar w/ generate_kepler_target_simple; ignored by generate_planetary_system_uncorrelated_incl

    return sim_param_closure
end

#=
function set_test_param(sim_param_closure::SimParam)
    return sim_param_closure
    add_param_fixed(sim_param_closure,"num_targets_sim_pass_one",5000)  # For faster simulated catalogs
    p_lim_arr_num = [10., 20.]
    r_lim_arr_num = [1.0, 1.25]
    rate_tab_init = reshape([2.5]*0.01, (1,1))
    add_param_fixed(sim_param_closure, "p_lim_arr", p_lim_arr_num)
    add_param_fixed(sim_param_closure, "r_lim_arr", r_lim_arr_num*ExoplanetsSysSim.earth_radius)
    add_param_active(sim_param_closure, "obs_par", rate_tab_init)
    return sim_param_closure
end
=#


text:
@testset "Dynamic" begin
    @testset "Check Dimensions" begin
        real_k = 3
        real_n = 3
        real_p = [1, 2, 3]
        real_nreps = 4
        real_T = 10

        Y = ones(Float64, real_T, real_n * real_nreps)
        F = [ones(Float64, real_n, real_p[j]) for j = 1:real_k]
        G = [ones(Float64, real_p[j], real_p[j]) for j = 1:real_k]

        n, nreps, p, T, k, index_map = DynamicMixtures.check_dimensions(Y, F, G)

        @test k == real_k
        @test n == real_n
        @test p == real_p
        @test nreps == real_nreps
        @test T == real_T

        @test index_map[1] == 1:3
        @test index_map[2] == 4:6
        @test index_map[3] == 7:9
        @test index_map[4] == 10:12

        #TODO: Add tests checking the throws
    end
end


text:
export free_cell_graph

"""
`free_cell_graph()` returns an `ImplicitGraph` whose vertices
are all possible FreeCell game positions in which there is an
edge from `S` to `T` if there's a game move from `S` to `T`.
"""
function free_cell_graph()
    vcheck(::Tableau) = true
    return ImplicitGraph{Tableau}(vcheck, move_maker)
end


text:
abstract type AbstractRadialBasisFunction <: InterpolationMethod end
abstract type RadialBasisFunction <: AbstractRadialBasisFunction end
abstract type GeneralizedRadialBasisFunction <: AbstractRadialBasisFunction end

Base.iterate(x::AbstractRadialBasisFunction) = (x, nothing)
Base.iterate(x::AbstractRadialBasisFunction, ::Any) = nothing

export  Gaussian,
        Multiquadratic,
        InverseQuadratic,
        InverseMultiquadratic,
        Polyharmonic,
        ThinPlate,
        GeneralizedMultiquadratic,
        GeneralizedPolyharmonic

# Define types for the different kinds of radial basis functions
for rbf in (:Gaussian,
            :Multiquadratic,
            :InverseQuadratic,
            :InverseMultiquadratic,
            :Polyharmonic)
    @eval begin
        struct $rbf{T} <: RadialBasisFunction where T <: Real
            ɛ::T
        end

        # Define default constructors
        $rbf() = $rbf(1)
    end
end

# Generalized RBF:s
struct GeneralizedMultiquadratic{T<:Real, S<:Real, U<:Integer} <: GeneralizedRadialBasisFunction
    ε::T
    β::S
    degree::U
end

struct GeneralizedPolyharmonic{S<:Integer, U<:Integer} <: GeneralizedRadialBasisFunction
    k::S
    degree::U
end

@doc "
    Gaussian(ɛ = 1)

Define a Gaussian Radial Basis Function

```math
ϕ(r) = e^{-(ɛr)^2}
```
" Gaussian
(rbf::Gaussian)(r) = exp(-(rbf.ɛ*r)^2)

@doc "
    Multiquadratic(ɛ = 1)

Define a Multiquadratic Radial Basis Function

```math
ϕ(r) = \\sqrt{1 + (ɛr)^2}
```
" Multiquadratic
(rbf::Multiquadratic)(r) = (sqrt(1 + (rbf.ɛ*r)^2))

@doc "
    InverseQuadratic(ɛ = 1)

Define an Inverse Quadratic Radial Basis Function

```math
ϕ(r) = \\frac{1}{1 + (ɛr)^2}
```
" InverseQuadratic
(rbf::InverseQuadratic)(r) = (1/(1 + (rbf.ɛ*r)^2))

@doc "
    InverseMultiquadratic(ɛ = 1)

Define an Inverse Multiquadratic Radial Basis Function

```math
ϕ(r) = \\frac{1}{\\sqrt{1 + (ɛr)^2}}
```
" InverseMultiquadratic
(rbf::InverseMultiquadratic)(r) = (1/sqrt(1 + (rbf.ɛ*r)^2))

@doc "
    Polyharmonic(k = 1)

Define a Polyharmonic Spline Radial Basis Function

```math
ϕ(r) = r^k, k = 1, 3, 5, ...
\\\\
ϕ(r) = r^k ln(r), k = 2, 4, 6, ...
```
" Polyharmonic
function (rbf::Polyharmonic{T})(r) where T <: Integer
    @assert rbf.ɛ > 0

    # Distinguish odd and even cases
    expr = if rbf.ɛ % 2 == 0
        (r > 0 ? r^rbf.ɛ*log(r) : 0.0)
    else
        (r^rbf.ɛ)
    end

    expr
end

@doc "
    ThinPlate()

Define a Thin Plate Spline Radial Basis Function

```math
ϕ(r) = r^2 ln(r)
```

This is a shorthand for `Polyharmonic(2)`.
" ThinPlate
ThinPlate() = Polyharmonic(2)

@doc "
    GeneralizedMultiquadratic(ɛ, β, degree)

Define a generalized Multiquadratic Radial Basis Function

```math
ϕ(r) = (1 + (ɛ*r)^2)^β
```
Results in a positive definite system for a 'degree' of ⌈β⌉ or higher.

" GeneralizedMultiquadratic
(rbf::GeneralizedMultiquadratic)(r) = (1 + (rbf.ε*r)^2)^rbf.β

@doc "
    GeneralizedPolyharmonic(k, degree)

Define a generalized Polyharmonic Radial Basis Function

```math
ϕ(r) = r^k, k = 1, 3, 5, ...
\\\\
ϕ(r) = r^k ln(r), k = 2, 4, 6, ...
```
Results in a positive definite system for a 'degree' of ⌈k/2⌉ or higher for k = 1, 3, 5, ...
and of exactly k + 1 for k = 2, 4, 6, ...
" GeneralizedPolyharmonic
function (rbf::GeneralizedPolyharmonic)(r)
    @assert rbf.k > 0

    # Distinguish odd and even cases
    expr = if rbf.k % 2 == 0
        (r > 0 ? r^rbf.k*log(r) : 0.0)
    else
        (r^rbf.k)
    end

    expr
end

abstract type RadialBasisInterpolant <: ScatteredInterpolant end

struct RBFInterpolant{T1, T2, F, M} <: RadialBasisInterpolant where {T1 <: AbstractArray, T2 <: AbstractMatrix{<:Real}}

    w::T1
    points::T2
    rbf::F
    metric::M
end

struct GeneralizedRBFInterpolant{T1, T2, F, M} <: RadialBasisInterpolant where {T1 <: AbstractArray, T2 <: AbstractMatrix{<:Real}}

    w::T1
    λ::T1
    points::T2
    rbf::F
    metric::M
end

function interpolate(rbf::Union{T, AbstractVector{T}} where T <: AbstractRadialBasisFunction,
                     points::AbstractArray{<:Real,2},
                     samples::AbstractArray{<:Number,N};
                     metric = Euclidean(), returnRBFmatrix::Bool = false,
                     smooth::Union{S, AbstractVector{S}} = false) where {N} where {S<:Number}

    #hinder smooth from being set to true and interpreted as the value 1
    @assert smooth != true "set the smoothing value as a number or vector of numbers"

    # Compute pairwise distances, apply the Radial Basis Function
    # and optional smoothing (ridge regression)
    A = pairwise(metric, points;dims=2)

    A = evaluateRBF!(A, rbf, smooth)

    # Solve for the weights
    itp = solveForWeights(A, points, samples, rbf, metric)

    # Create and return an interpolation object
    if returnRBFmatrix    # Return matrix A
        return itp, A
    else
        return itp
    end

end

@inline function evaluateRBF!(A, rbf, symmetric::Bool = false)
    A .= rbf.(A)

    if symmetric
        A = Symmetric(A)
    end

    return A
end
@inline function evaluateRBF!(A, rbfs::AbstractVector{<:AbstractRadialBasisFunction}, symmetric::Bool = false)
    for (j, rbf) in enumerate(rbfs)
        A[:,j] .= rbf.(@view A[:,j])
    end
    return A
end
@inline function evaluateRBF!(A, rbf, smooth)
    A = evaluateRBF!(A, rbf)
    A = addSmoothing!(A, smooth)
    return A
end

@inline function addSmoothing!(A, smooth::Vector{T}) where T <: Number
    for i = 1:size(A, 1)
        A[i,i] += smooth[i]
    end
    return A
end
@inline function addSmoothing!(A, smooth::T) where T <: Number
    for i = 1:size(A, 1)
        A[i,i] += smooth
    end
    return A
end

@inline function solveForWeights(A, points, samples,
                                    rbf::Union{T, AbstractVector{T}} where T <: RadialBasisFunction,
                                    metric)
    w = A\samples
    RBFInterpolant(w, points, rbf, metric)
end
@inline function solveForWeights(A, points, samples,
                                    rbf::Union{T, AbstractVector{T}} where T <: Union{GeneralizedRadialBasisFunction, RadialBasisFunction},
                                    metric)
    # Use the maximum degree among the generalized RBF:s
    P = getPolynomial(rbf, points)

    # Solve for the weights and polynomial coefficients
    # We end up with a blocked system, so we don't have to form the full matrix
    Af = factorize(A)
    B = -P'*(Af\P)
    E = B\(P'*(Af\samples))

    w = A\(I*samples + P*E)
    λ = -E

    GeneralizedRBFInterpolant(w, λ, points, rbf, metric)
end

function evaluate(itp::RadialBasisInterpolant, points::AbstractArray{<:Real, 2})

    # Compute distance matrix and evaluate the RBF
    A = pairwise(itp.metric, points, itp.points; dims=2)
    A = evaluateRBF!(A, itp.rbf, false)

    # Compute polynomial matrix for generalized RBF:s
    P = getPolynomial(itp.rbf, points)

    # Compute the interpolated values
    return computeInterpolatedValues(A, P, itp)
end

@inline getPolynomial(rbf::Union{T, AbstractVector{T}} where T <: RadialBasisFunction, points) = nothing
@inline function getPolynomial(rbf::Union{T, AbstractVector{T}} where T <: Union{GeneralizedRadialBasisFunction, RadialBasisFunction}, points)

    # Use the maximum degree among the generalized RBF:s
    degree = maximum(x isa GeneralizedRadialBasisFunction ? x.degree : 0 for x in rbf)
    P = generateMultivariatePolynomial(points, degree)
end

@inline computeInterpolatedValues(A, P, itp::RBFInterpolant) = A*itp.w
@inline computeInterpolatedValues(A, P, itp::GeneralizedRBFInterpolant) = A*itp.w + P*itp.λ

# Helper function to generate matrices defining complete homogenic symmetric polynomials
function generateMultivariatePolynomial(points::AbstractArray{<:Real, 2}, degree::Integer)

    # How big should the matrix be?
    nDimensions, nPoints = size(points)
    nTerms = binomial(degree + nDimensions, degree)

    P = ones(nPoints, nTerms)

    # Start with the lowest orders and work upwards
    position = 2
    while position <= nTerms
        for order = 1:degree
            for combination in with_replacement_combinations(1:nDimensions, order)
                for var in combination
                    P[:, position] .*= points[var, :]
                end
                position += 1
            end
        end
    end

    P
end


text:
"""
$(TYPEDEF)

Common options for all methods provided via `DataDrivenDiffEq`.
"""
mutable struct DataDrivenCommonOptions{T,K}
    """Maximum iterations"""
    maxiter::Int
    """Absolute tolerance"""
    abstol::T
    """Relative tolerance"""
    reltol::T

    """Show a progress"""
    progress::Bool
    """Display log - Not implemented right now"""
    verbose::Bool

    """Denoise the data using singular value decomposition"""
    denoise::Bool
    """Normalize the data"""
    normalize::Bool
    """Sample options, see `DataSampler`"""
    sampler::AbstractSampler

    """Mapping from the candidate solution of a problem to features used for pareto analysis"""
    f::Function
    """Scalarization of the features for a candidate solution"""
    g::Function
    """Significant digits for the parameters - used for rounding. Default = 10"""
    digits::Int

    """Additional kwargs"""
    kwargs::K
end

DataDrivenCommonOptions(opt::AbstractKoopmanAlgorithm, ::Type{T} = Float64, args...;
    maxiter = 100, abstol = sqrt(eps()), reltol = sqrt(eps()),
    progress = false, verbose = false,
    denoise = false, normalize = false,
    sampler = DataSampler(),
    f = F(opt), g = G(opt), digits = 10, kwargs...) where T = begin
   DataDrivenCommonOptions{eltype(T), typeof(kwargs)}(
       maxiter, abstol, reltol,
       progress, verbose, denoise, normalize,
       sampler, f, g, digits, kwargs
   )
end

DataDrivenCommonOptions(opt::AbstractOptimizer{T}, args...;
    maxiter = 100, abstol = sqrt(eps()), reltol = sqrt(eps()),
    progress = false, verbose = false,
    denoise = false, normalize = false,
    sampler = DataSampler(),
    f = F(opt), g = G(opt),
    digits = 10,
    kwargs...) where T = begin
   DataDrivenCommonOptions{eltype(T), typeof(kwargs)}(
       maxiter, abstol, reltol,
       progress, verbose, denoise, normalize,
       sampler, f, g, digits, kwargs
   )
end

## Normalization etc
function normalize_theta!(scales::AbstractVector, theta::AbstractMatrix)
    map(1:length(scales)) do i
        scales[i] = norm(theta[i,:], 2)
        theta[i, :] .= theta[i,:]./scales[i]
    end
    return
end

function rescale_xi!(xi::AbstractMatrix, scales::AbstractVector, round_::Bool)
    digs = 10
    @inbounds for i in 1:length(scales), j in 1:size(xi, 2)
        iszero(xi[i,j]) ? continue : nothing
        round_ && (xi[i,j] % 1) != zero(xi[i,j]) ? digs = round(Int64,-log10(abs(xi[i,j]) % 1))+1 : nothing
        xi[i,j] = xi[i,j] / scales[i]
        round_ ? xi[i,j] = round(xi[i,j], digits = digs) : nothing
    end
    return
end

function candidate_matrix(b::Basis, n_o::Int)
    eqs = map(x->x.rhs, equations(b))
    xs = states(b)
    ys = implicit_variables(b)

    isempty(ys) && return ones(Bool, n_o, length(eqs))

    c = zeros(Bool, length(ys), length(eqs))

    for i in 1:length(ys), j in 1:length(eqs)
        c[i,j] = is_dependent(Num(eqs[j]), Num(ys[i]))
        c[i,j] && continue
        c[i,j] = all(map(xi->is_not_dependent(Num(eqs[j]), Num(xi)), ys))
    end

    return c
end


text:
module ImageToEIS

#using ElectricalEquivalentCircuits
using ImageIO
using FileIO
using IterativeSolvers
using PyPlot
using Dates
using Printf
using Colors
using DataFrames
using CSV
using Interpolations
using Optim
using SparseArrays
using Statistics

using Base

include("physical_properties.jl")
export matrix_to_file

# general part

include("material_matrix_to_lin_sys.jl")
include("material_matrix_to_impedance.jl")
export material_matrix_to_impedance

# user requested part

include("equivalent_circuit_support.jl")
include("Z_view_export.jl")
include("image_to_EIS_interface.jl")
export image_to_EIS

# maybe independent package for proper geometry generation

include("generate_matrix.jl")
export generate_matrix
export three_column_domain_template
export three_column_domain_matrix
export chess_matrix

# macro stuff consisting of par study, evaluate results, ploting, cluster operations

include("macro_stuff.jl")
export par_study

# specific material dependences

include("temperature_dependences.jl")
export TI
export TI_2D

end # module

max_stars_repo_name: Masicko/ImageToEIS

max_issues_repo_name: Masicko/ImageToEIS

max_forks_repo_name: Masicko/ImageToEIS


text:
### interacting with the system
julia>;

shell> pwd
/Users/rahullakhanpal/learning-julia

# or alternatively, using
julia> pwd()
"/Users/rahullakhanpal/learning-julia"

# to know your home directory
julia> homedir()
"/Users/rahullakhanpal"

# to join two paths
julia> joinpath(homedir(), "Documents")
"/Users/rahullakhanpal/Documents"

# read the contents of a directory
julia> readdir(pwd())
8-element Array{String,1}:
 ".git"
 ".gitignore"
 "chapter-1"
 "chapter-2"
 "chapter-3"
 "chapter-4"
 "chapter-5"
 "chapter-6"

# stat over any random file
julia> stat(".gitignore")
StatStruct(mode=0o100644, size=21)

# check if its a file
julia> isfile(".gitignore")
true

# check if its a directory
julia> isdir(".git")
true

# knowing basename of the home directory
julia> basename(homedir())
"rahullakhanpal"

# to split the directory structure
julia> splitdir(homedir())
("/Users", "rahullakhanpal")

# to check if its a path
julia> ispath(homedir())
true

# to know absolute path of current directory
julia> abspath(".")
"/Users/rahullakhanpal/learning-julia/"

### I/O Operations

# to open a file and read its contents in Julia
julia> file = open("sample.txt")
IOStream(<file sample.txt>)
julia> readlines(file)
2-element Array{String,1}:
 "Hi there! my name is Rahul."
 "Learning Julia is so much fun."

# to write into a file in julia
julia> write("sample.txt", "Thats a line using the write() function")
39
julia> readlines("sample.txt")
1-element Array{String,1}:
 "Thats a line using the write() function"

# close a file, where,
# file = "sample.txt"
julia> close(file)

### Example to create a shell script that accepts arguments

# create a file input.txt in the current directory and put the following in it
shell> cat input.txt
Julia is a high-level, high-performance dynamic programming language for numerical computing. It provides a sophisticated compiler, distributed parallel execution, numerical accuracy, and an extensive mathematical function library. Julia’s Base library, largely written in Julia itself, also integrates mature, best-of-breed open source C and Fortran libraries for linear algebra, random number generation, signal processing, and string processing. In addition, the Julia developer community is contributing a number of external packages through Julia’s built-in package manager at a rapid pace. IJulia, a collaboration between the Jupyter and Julia communities, provides a powerful browser-based graphical notebook interface to Julia.

shell> cat sample.jl
# Arguments
# ARGS[0] is the command itself
in_file = ARGS[1]
out_file = ARGS[2]

# Keeping track using a counter
counter = 0
for line in eachline(in_file)
	for word in split(line)
		if word == "Julia"
			counter += 1
		end
	end
end

# Write the contents to the o/p file
write(out_file, "The count for the word julia is $counter")

# Finally, read the contents to inform the user
for line in readlines(out_file)
	println(line)
end

# run the script like
shell> julia sample.jl input.txt out.txt
The count for the word julia is 4

### calling C from julia

julia> ccall((:clock, :libc), Int64, ())
5818283

julia> syspath = ccall( (:getenv, :libc), Ptr{Cchar}, (Ptr{Cchar},), "SHELL")
Ptr{Int8} @0x00007ffee62cfc00

julia> unsafe_string(syspath)
"/bin/bash"

# Python
# simple hello world
julia> using PyCall
INFO: Precompiling module PyCall.

julia> py"""
         print 'hello world'
         """
hello world

# calling a function
julia> py"""
       def cube(n):
           return n*n*n
       print cube(10)
       """
1000

# imports in python
julia> @pyimport numpy as np

julia> np.sin(120)
0.5806111842123143

# Using pybuiltin, directy create a python dict object
julia> pybuiltin(:dict)(a=1,b=2)
Dict{Any,Any} with 2 entries:
  "b" => 2
  "a" => 1

julia> pycall(pybuiltin("dict"), Any, a=1, b=2)
PyObject {'a': 1, 'b': 2}

julia> d = pycall(pybuiltin("dict"), Any, a=1, b=2)
PyObject {'a': 1, 'b': 2}

# making sure of the type, still a python object
julia> typeof(d)
PyCall.PyObject

# using pyDict to convert python's object to a juia dict object
julia> julia_dictionary = PyDict{Symbol, Int64}(pycall(pybuiltin("dict"), Any, a=1, b=2))
PyCall.PyDict{Symbol,Int64,true} with 2 entries:
  :a => 1
  :b => 2

# as we see, its now a julia dict object
julia> typeof(julia_dictionary)
PyCall.PyDict{Symbol,Int64,true}

# accessing dictionary elements
julia> julia_dictionary[:a]
1

# accessing dictionary elements
julia> julia_dictionary[:b]
2

# a function that takes kw arguments
julia> f(; a=0, b=0) = [10a, b]
f (generic function with 1 method)

julia> f(;julia_dictionary...)
2-element Array{Int64,1}:
 10
  2

### Expressions and Macros

julia> code = "println(\"hello world \")"
"println(\"hello world \")"

julia> expression = parse(code)
:(println("hello world "))

julia> typeof(expression)
Expr

julia> fieldnames(expression)
3-element Array{Symbol,1}:
 :head
 :args
 :typ

julia> expression.args
2-element Array{Any,1}:
 :println
 "hello world "

julia> expression.head
:call

julia> expression.typ
Any

julia> dump(expression)
Expr
  head: Symbol call
  args: Array{Any}((2,))
    1: Symbol println
    2: String "hello world "
  typ: Any

julia> result = :hello
:hello

julia> typeof(result)
Symbol

julia> fieldnames(result)
0-element Array{Symbol,1}

julia> sample_expr = Expr(:call, +, 10, 20)
:((+)(10,20))

julia> eval(sample_expr)
30

julia> sample_expr.args
3-element Array{Any,1}:
   +
 10
 20

julia> eval(sample_expr)
20

# interpolation of variables
julia> x = 10
10

julia> y = 10
10

julia> sample_expr = Expr(:call, +, :x, :y)
:((+)(10,10))

# or even this works
julia> sample_expr = Expr(:call, +, x, y)
:((+)(10,10))

julia> eval(sample_expr)
20

# interpolation using quotes(:)
julia> quote
       30 * 100
       end
quote  # REPL[12], line 2:
    30 * 100
end

julia> eval(ans)
3000

julia> :(30 * 100)
:(30 * 100)

julia> eval(ans)
3000

# interpolation of varibales using dollar($)
julia> x = 10
10

julia> y = 10
10

julia> e = :($x + $y)
:(10 + 10)

julia> eval(e)
20

# Macros
julia> macro HELLO(name)
                  :( println("Hello! ", $name))
              end
@HELLO (macro with 1 method)

julia> @HELLO("rahul")
Hello! rahul

julia> macroexpand(:(@HELLO("rahul")))
:(println("Hello!","rahul"))

# Example of using macro
julia> for header in [:h1, :h2, :h3, :h4, :h5, :h6]
                         @eval function $(Symbol(header))(text::String)
                             println("<" * $(string(header))* ">"  * " $text " * "</" * $(string(header))* ">")
                         end
                     end

julia> h1("Hello world!")
<h1> Hello world! </h1>

julia> h2("Hello world!")
<h2> Hello world! </h2>

### Built in macros

# checking the time taken by a function to run, using @time macro
# simple function to find recursive sum
julia> function recursive_sum(n)
           if n == 0
               return 0
           else
               return n + recursive_sum(n-1)
       end
       end
recursive_sum (generic function with 1 method)

# A bit slow to run for the 1st Time, as the function gets compiled.
julia> @time recursive_sum(10000)
  0.003905 seconds (450 allocations: 25.816 KiB)
50005000

# Much much faster in the second run!
julia> @time recursive_sum(10000)
  0.000071 seconds (5 allocations: 176 bytes)
50005000

# @elapsed
julia> @elapsed average(10000000, 1000000000)
2.144e-6

julia> typeof(@elapsed average(10000000, 1000000000))
Float64

# @show
julia> @show(println("hello world"))
hello world
println("hello world") = nothing

julia> @show(:(println("hello world")))
$(Expr(:quote, :(println("hello world")))) = :(println("hello world"))
:(println("hello world"))

julia> @show(:(3*2))
$(Expr(:quote, :(3 * 2))) = :(3 * 2)
:(3 * 2)

julia> @show(3*2)
3 * 2 = 6
6

julia> @show(Int64)
Int64 = Int64
Int64

# @which
# create a function that tripples an Integer
julia> function tripple(n::Int64)
           3n
       end
tripple (generic function with 1 method)

# redefine the same function to accept Float
julia> function tripple(n::Float64)
           3n
       end
tripple (generic function with 2 methods)

# check the methods available for this function
julia> methods(tripple)
# 2 methods for generic function "tripple":
tripple(n::Float64) in Main at REPL[22]:2
tripple(n::Int64) in Main at REPL[21]:2

# Get the correct method , when 'n' is an Int64
julia> @which tripple(10)
tripple(n::Int64) in Main at REPL[21]:2

# Get the correct method , when 'n' is Float64
julia> @which tripple(10.0)
tripple(n::Float64) in Main at REPL[22]:2

# @task
julia> say_hello() = println("hello world")
say_hello (generic function with 1 method)

julia> say_hello_task = @task say_hello()
Task (runnable) @0x000000010dcdfa90

julia> istaskstarted(say_hello_task)
false

julia> schedule(say_hello_task)
hello world
Task (queued) @0x000000010dcdfa90

julia> yield()

julia> istaskdone(say_hello_task)
true

# @code_lowered, @code_typed, @code_llvm and @code_native
julia> function fibonacci(n::Int64)
           if n < 2
               n
           else
               fibonacci(n-1) + fibonacci(n-2)
           end
       end
fibonacci (generic function with 1 method)

# OR, can also define it this way
julia> fibonacci(n::Int64) = n < 2 ? n : fibonacci(n-1) + fibonacci(n-2)
fibonacci (generic function with 1 method)

julia> fibonacci(10)
55

julia> @code_lowered fibonacci(10)
CodeInfo(:(begin
        nothing
        unless n < 2 goto 4
        return n
        4:
        return (Main.fibonacci)(n - 1) + (Main.fibonacci)(n - 2)
    end))

julia> @code_typed fibonacci(10)
CodeInfo(:(begin
        unless (Base.slt_int)(n, 2)::Bool goto 3
        return n
        3:
        SSAValue(1) = $(Expr(:invoke, MethodInstance for fibonacci(::Int64), :(Main.fibonacci), :((Base.sub_int)(n, 1)::Int64)))
        SSAValue(0) = $(Expr(:invoke, MethodInstance for fibonacci(::Int64), :(Main.fibonacci), :((Base.sub_int)(n, 2)::Int64)))
        return (Base.add_int)(SSAValue(1), SSAValue(0))::Int64
    end))=>Int64

julia> @code_warntype fibonacci(10)
Variables:
  #self#::#fibonacci
  n::Int64

Body:
  begin
      unless (Base.slt_int)(n::Int64, 2)::Bool goto 3
      return n::Int64
      3:
      SSAValue(1) = $(Expr(:invoke, MethodInstance for fibonacci(::Int64), :(Main.fibonacci), :((Base.sub_int)(n, 1)::Int64)))
      SSAValue(0) = $(Expr(:invoke, MethodInstance for fibonacci(::Int64), :(Main.fibonacci), :((Base.sub_int)(n, 2)::Int64)))
      return (Base.add_int)(SSAValue(1), SSAValue(0))::Int64
  end::Int64

julia> @code_llvm fibonacci(10)

define i64 @julia_fibonacci_61143.2(i64) #0 !dbg !5 {
top:
  %1 = icmp sgt i64 %0, 1
  br i1 %1, label %L3, label %if

if:                                               ; preds = %top
  ret i64 %0

L3:                                               ; preds = %top
  %2 = add i64 %0, -1
  %3 = call i64 @julia_fibonacci_61143(i64 %2)
  %4 = add i64 %0, -2
  %5 = call i64 @julia_fibonacci_61143(i64 %4)
  %6 = add i64 %5, %3
  ret i64 %6
}

julia> @code_native fibonacci(10)
  .section  __TEXT,__text,regular,pure_instructions
Filename: REPL[50]
  pushq %rbp
  movq  %rsp, %rbp
  pushq %r15
  pushq %r14
  pushq %rbx
  pushq %rax
  movq  %rdi, %rbx
Source line: 1
  cmpq  $1, %rbx
  jle L63
  leaq  -1(%rbx), %rdi
  movabsq $fibonacci, %r15
  callq *%r15
  movq  %rax, %r14
  addq  $-2, %rbx
  movq  %rbx, %rdi
  callq *%r15
  addq  %r14, %rax
  addq  $8, %rsp
  popq  %rbx
  popq  %r14
  popq  %r15
  popq  %rbp
  retq
L63:
  movq  %rbx, %rax
  addq  $8, %rsp
  popq  %rbx
  popq  %r14
  popq  %r15
  popq  %rbp
  retq
  nopl  (%rax)

## Type introspection
# making use of typeof

julia> type Student
           name::String
           age::Int64
       end

julia> alpha = Student("alpha",24)
Student("alpha", 24)

julia> beta = Student("beta",25)
Student("beta", 25)

julia> typeof(alpha)
Student

# similar to isinstance in python
julia> isa(alpha, Student)
true

## Reflection
# the first method tries to take in all integer values
julia> function calculate_quad(a::Int64,b::Int64,c::Int64,x::Int64)
           return a*x^2 + b*x + c
       end
calculate_quad (generic function with 2 methods)

julia> calculate_quad(1,2,3,4)
27

# the second method takes all but x as integer values
julia> function calculate_quad(a::Int64,b::Int64,c::Int64,x::Float64)
           return a*x^2 + b*x + c
       end
calculate_quad (generic function with 3 methods)

julia> calculate_quad(1,2,3,4.75)
35.0625

# to know what all methods does the function supports
# which as we can see that there are 2 currently
julia> methods(calculate_quad)
# 3 methods for generic function "calculate_quad":
calculate_quad(a::Int64, b::Int64, c::Int64, x::Float64) in Main at REPL[31]:2
calculate_quad(a::Int64, b::Int64, c::Int64, x::Int64) in Main at REPL[29]:2
calculate_quad(a, b, c, x) in Main at REPL[27]:2

# from the already declared Student class
julia> fieldnames(Student)
2-element Array{Symbol,1}:
 :name
 :age

julia> Student.types
svec(String, Int64)

max_stars_repo_name: PacktPublishing/Learning-Julia

max_issues_repo_name: PacktPublishing/Learning-Julia

max_forks_repo_name: PacktPublishing/Learning-Julia


text:
export DefiniteIntegral,DefiniteLineIntegral

abstract CalculusFunctional{S,T} <: Operator{T}

@functional CalculusFunctional

##TODO: Add ConcreteOp

macro calculus_functional(Op)
    ConcOp=parse("Concrete"*string(Op))
    WrappOp=parse(string(Op)*"Wrapper")
    return esc(quote
        abstract $Op{SSS,TTT} <: CalculusFunctional{SSS,TTT}
        immutable $ConcOp{S,T} <: $Op{S,T}
            domainspace::S
        end
        immutable $WrappOp{BT<:Operator,S<:Space,T} <: $Op{S,T}
            func::BT
        end

        # We expect the operator to be real/complex if the basis is real/complex
        $ConcOp(dsp::Space) = $ConcOp{typeof(dsp),eltype(dsp)}(dsp)

        $Op() = $Op(UnsetSpace())
        $Op(dsp) = $ConcOp(dsp)
        $Op(d::Domain) = $Op(Space(d))

        promotedomainspace(::$Op,sp::Space) = $Op(sp)

        Base.convert{T}(::Type{Operator{T}},Σ::$ConcOp) =
            T==eltype(Σ)?Σ:$ConcOp{typeof(Σ.domainspace),T}(Σ.domainspace)

        domain(Σ::$ConcOp) = domain(Σ.domainspace)
        domainspace(Σ::$ConcOp) = Σ.domainspace

        getindex(::$ConcOp{UnsetSpace},kr::Range) =
            error("Spaces cannot be inferred for operator")

        $WrappOp(op::Operator) =
            $WrappOp{typeof(op),typeof(domainspace(op)),eltype(op)}(op)

        Base.convert{T}(::Type{Operator{T}},Σ::$WrappOp) =
            T==eltype(Σ)?Σ:$WrappOp(convert(Operator{T},Σ.func))

        #Wrapper just adds the operator it wraps
        getindex(D::$WrappOp,k::Range) = D.func[k]
        getindex(D::$WrappOp,k::Integer) = D.func[k]
        domainspace(D::$WrappOp) = domainspace(D.func)
        bandinds(D::$WrappOp) = bandinds(D.func)
    end)
end

@calculus_functional(DefiniteIntegral)
@calculus_functional(DefiniteLineIntegral)

#default implementation
function getindex(B::ConcreteDefiniteIntegral,kr::Range)
    S=domainspace(B)
    Q=Integral(S)
    A=(Evaluation(S,true)-Evaluation(S,false))*Q
    A[kr]
end

function getindex(B::ConcreteDefiniteIntegral,kr::Integer)
    S=domainspace(B)
    Q=Integral(S)
    A=(Evaluation(S,true)-Evaluation(S,false))*Q
    A[kr]
end

function DefiniteIntegral(sp::Space)
    if typeof(canonicaldomain(sp)).name==typeof(domain(sp)).name
        ConcreteDefiniteIntegral{typeof(sp),eltype(sp)}(sp)
    else
        M = Multiplication(fromcanonicalD(sp),setcanonicaldomain(sp))
        Op = DefiniteIntegral(rangespace(M))*M
        DefiniteIntegralWrapper(SpaceOperator(Op,sp,rangespace(Op)))
    end
end

function DefiniteLineIntegral(sp::Space)
    if typeof(canonicaldomain(sp)).name==typeof(domain(sp)).name
        ConcreteDefiniteLineIntegral{typeof(sp),eltype(sp)}(sp)
    else
        M = Multiplication(abs(fromcanonicalD(sp)),setcanonicaldomain(sp))
        Op = DefiniteLineIntegral(rangespace(M))*M
        DefiniteLineIntegralWrapper(SpaceOperator(Op,sp,rangespace(Op)))
    end
end

#TODO: Remove SPECIALOPS reimplement
# *{T,D<:DefiniteIntegral,M<:Multiplication}(A::TimesFunctional{T,D,M},b::Fun) = bilinearform(A.op.f,b)
# *{T,D<:DefiniteLineIntegral,M<:Multiplication}(A::TimesFunctional{T,D,M},b::Fun) = linebilinearform(A.op.f,b)
# *{T,D<:Union{DefiniteIntegral,DefiniteLineIntegral},
#   M<:Multiplication,V}(A::FunctionalOperator{TimesFunctional{T,D,M},V},b::Fun) =
#     Fun(A.func*b)


text:
export elementary, bi_elementary
export elementary_waypoints, bi_elementary_waypoints

function elementary((x0, y0, θ0)::StaticVector{3,T0}, (xf, yf, θf)::StaticVector{3,Tf}; v=1, λ=1//2) where {T0,Tf}
    ψ=atan(yf - y0, xf - x0)
    @assert adiff(ψ, θ0) ≈ adiff(θf, ψ) "$((x0, y0, θ0)) and $((xf, yf, θf)) are not symmetric for elementary path"
    @assert abs(adiff(ψ, θ0)) < π/2     "$((xf, yf, θf)) is not in front of $((x0, y0, θ0)) for elementary path"
    T = promote_type(T0, Tf)
    α = adiff(θf, θ0)
    d = hypot(xf - x0, yf - y0)
    α ≈ 0 && return (cost=d/v, controls=SVector(zero(VelocityCurvRateStep{T}),
                                                VelocityCurvRateStep{T}(d/v, VelocityCurvRateControl{T}(v, 0)),
                                                zero(VelocityCurvRateStep{T})))
    L = d/D(α, λ)
    σ = 4*α/(L*L*(1 - λ*λ))
    ctrl = SVector(VelocityCurvRateStep{T}(L*(1 - λ)/2, VelocityCurvRateControl{T}(1, σ)),
                   VelocityCurvRateStep{T}(L*λ, VelocityCurvRateControl{T}(1, 0)),
                   VelocityCurvRateStep{T}(L*(1 - λ)/2, VelocityCurvRateControl{T}(1, -σ)))
    (cost=L/v, controls=scalespeed.(ctrl, v))
end

function bi_elementary(q0::StaticVector{3,T0}, qf::StaticVector{3,Tf}; v=1, γ=1//2, λ=1//2) where {T0,Tf}
    T = promote_type(T0, Tf)
    x0, y0, θ0 = q0
    xf, yf, θf = qf
    @assert adiff(θ0, θf) ≈ 0 "angles of $q0 and $qf must be equal for bi-elementary path"
    q0 ≈ qf && return (cost=T(0), controls=zeros(SVector{6,VelocityCurvRateStep{T}}))
    dx = xf - x0
    dy = yf - y0
    ψ  = atan(dy, dx)
    θi = mod2piF(2*ψ - θ0)
    qi = SE2State(x0 + γ*dx, y0 + γ*dy, θi)
    ctrl = [elementary(q0, qi, v=v, λ=λ).controls; elementary(qi, qf, v=v, λ=λ).controls]
    (cost=duration(ctrl), controls=ctrl)
end

function elementary_waypoints(q0::StaticVector, qf::StaticVector, dt_or_N; v=1, λ=1//2)
    waypoints(SimpleCarDynamics{0,1}(), SE2κState(q0), elementary(q0, qf, v=v, λ=λ).controls, dt_or_N)
end

function bi_elementary_waypoints(q0::StaticVector, qf::StaticVector, dt_or_N; v=1, γ=1//2, λ=1//2)
    waypoints(SimpleCarDynamics{0,1}(), SE2κState(q0), bi_elementary(q0, qf, v=v, γ=γ, λ=λ).controls, dt_or_N)
end

@inline function D(α::T, λ) where {T}
    α = abs(α)
    c = α/(1 + λ)
    s = sqrt(α/(T(π)*(1 - λ*λ)))
    sin(c*λ)/c + (cos(α/2)*fresnelC(s*(1 - λ)) + sin(α/2)*fresnelS(s*(1 - λ)))/s
end


text:
module NeuronModel

export izhikevicmodel_fire
export izhikevicmodel_step
export thresh

const N = 1000
const Ne = 800
const a = [(i<=Ne) ? 0.02 : 0.1 for i in 1:N]
const d = [(i<=Ne) ? 8.0 : 2.0 for i in 1:N]
const thresh = 30

function izhikevicmodel_step(v::Array{Float64},u::Array{Float64},I::Array{Float64}) # Module neuron_model
    v .= @. v+0.5*((0.04*v+5)*v+140-u+I)
    v .= @. v+0.5*((0.04*v+5)*v+140-u+I)
    u .= @. u+a*(0.2*v-u)
    return v,u
end

function izhikevicmodel_fire(v::Array{Float64},u::Array{Float64},fired::Array{Int64}) # Module neuron_model
    if length(fired) != 0
        v[fired] .= -65.0
        u[fired] .= u[fired] .+ d[fired]
    end
    return v,u
end

end #end of module

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
using IntSets
using Base.Test

include("intset.jl")


text:
module MyPackage

"""
Prints a greetin message.
"""
greet() = print("Hello World!")

"""
    addone(x)

Adds one to `x`.
"""
addone(x) = x + 1

export greet, addone

end # module


text:
module MATT

include("./Component.jl")
include("./React.jl")
inlcude("./Server.jl")

greet() = print("Hello World!")

end # module


text: is_push_required(gl::GitLink) = !_is_stage_up_to_day(gl)


text:
using CPUTime

function y(x_start, dx, x_end)
    y = []
    for xi in x_start:dx:x_end
        push!(y, xi^2 + xi + 1)
        println(xi)
    end
    return y
end

@CPUtime ylist = y(0, 1, 10)
println(ylist)


text:
@testset "get_mesh_surface.jl" begin

end


text:
#=
--- Day 3: Toboggan Trajectory ---

With the toboggan login problems resolved, you set off toward the airport. While travel by toboggan might be easy, it's certainly not safe: there's very minimal steering and the area is covered in trees. You'll need to see which angles will take you near the fewest trees.

Due to the local geology, trees in this area only grow on exact integer coordinates in a grid. You make a map (your puzzle input) of the open squares (.) and trees (#) you can see. For example:

..##.......
#...#...#..
.#....#..#.
..#.#...#.#
.#...##..#.
..#.##.....
.#.#.#....#
.#........#
#.##...#...
#...##....#
.#..#...#.#

These aren't the only trees, though; due to something you read about once involving arboreal genetics and biome stability, the same pattern repeats to the right many times:

..##.........##.........##.........##.........##.........##.......  --->
#...#...#..#...#...#..#...#...#..#...#...#..#...#...#..#...#...#..
.#....#..#..#....#..#..#....#..#..#....#..#..#....#..#..#....#..#.
..#.#...#.#..#.#...#.#..#.#...#.#..#.#...#.#..#.#...#.#..#.#...#.#
.#...##..#..#...##..#..#...##..#..#...##..#..#...##..#..#...##..#.
..#.##.......#.##.......#.##.......#.##.......#.##.......#.##.....  --->
.#.#.#....#.#.#.#....#.#.#.#....#.#.#.#....#.#.#.#....#.#.#.#....#
.#........#.#........#.#........#.#........#.#........#.#........#
#.##...#...#.##...#...#.##...#...#.##...#...#.##...#...#.##...#...
#...##....##...##....##...##....##...##....##...##....##...##....#
.#..#...#.#.#..#...#.#.#..#...#.#.#..#...#.#.#..#...#.#.#..#...#.#  --->

You start on the open square (.) in the top-left corner and need to reach the bottom (below the bottom-most row on your map).

The toboggan can only follow a few specific slopes (you opted for a cheaper model that prefers rational numbers); start by counting all the trees you would encounter for the slope right 3, down 1:

From your starting position at the top-left, check the position that is right 3 and down 1. Then, check the position that is right 3 and down 1 from there, and so on until you go past the bottom of the map.

The locations you'd check in the above example are marked here with O where there was an open square and X where there was a tree:

..##.........##.........##.........##.........##.........##.......  --->
#..O#...#..#...#...#..#...#...#..#...#...#..#...#...#..#...#...#..
.#....X..#..#....#..#..#....#..#..#....#..#..#....#..#..#....#..#.
..#.#...#O#..#.#...#.#..#.#...#.#..#.#...#.#..#.#...#.#..#.#...#.#
.#...##..#..X...##..#..#...##..#..#...##..#..#...##..#..#...##..#.
..#.##.......#.X#.......#.##.......#.##.......#.##.......#.##.....  --->
.#.#.#....#.#.#.#.O..#.#.#.#....#.#.#.#....#.#.#.#....#.#.#.#....#
.#........#.#........X.#........#.#........#.#........#.#........#
#.##...#...#.##...#...#.X#...#...#.##...#...#.##...#...#.##...#...
#...##....##...##....##...#X....##...##....##...##....##...##....#
.#..#...#.#.#..#...#.#.#..#...X.#.#..#...#.#.#..#...#.#.#..#...#.#  --->

In this example, traversing the map using this slope would cause you to encounter 7 trees.

Starting at the top-left corner of your map and following a slope of right 3 and down 1, how many trees would you encounter?
=#

function read_input()::Array{String}

    input::Array{String} = []
    open("./input", "r") do f
        while ! eof(f)
            s::String = readline(f)
            push!(input,s)
        end
    end

    return input
end
input = read_input()

function count_trees(input)

    trees = 0
    displacement = 1
    pattern_len = length(input[1])

    for i in input[2:end]
        displacement+=3

        # clock arithmetic
        index = mod(displacement,pattern_len)
        if index==0
            index = pattern_len
        end

        encounter = i[index]
        if encounter == '#'
            trees+=1
        end
    end
    return trees
end

println("Trees found in choosen trajectory: ", count_trees(input))

#=
--- Part Two ---

Time to check the rest of the slopes - you need to minimize the probability of a sudden arboreal stop, after all.

Determine the number of trees you would encounter if, for each of the following slopes, you start at the top-left corner and traverse the map all the way to the bottom:

    Right 1, down 1.
    Right 3, down 1. (This is the slope you already checked.)
    Right 5, down 1.
    Right 7, down 1.
    Right 1, down 2.

In the above example, these slopes would find 2, 7, 3, 4, and 2 tree(s) respectively; multiplied together, these produce the answer 336.

What do you get if you multiply together the number of trees encountered on each of the listed slopes?
=#

function generic_count_trees(input,move_right,move_down)

    trees = 0
    displacement = 1
    pattern_len = length(input[1])

    for i in input[move_down+1:move_down:end]
        displacement+=move_right

        # clock arithmetic
        index = mod(displacement,pattern_len)
        if index==0
            index = pattern_len
        end

        encounter = i[index]
        if encounter == '#'
            trees+=1
        end
    end
    return trees
end

function compute_trees_for_each_slope(input)

    multiplication = 1
    slopes = [(1,1), (3,1), (5,1), (7,1), (1,2)]
    for slope in slopes
        trees = generic_count_trees(input,slope[1],slope[2])
        println("Trees with slope ", slope, ": ", trees)
        multiplication *= trees
    end
    println("Solution: ", multiplication)
end

println()
compute_trees_for_each_slope(input)


text:
type ScatterPlot <: Graph
	dataGroups::Vector{DataGroup}
end

ScatterPlot(dg::DataGroup) = ScatterPlot([dg])

function plotData(p::PlotFrame, g::ScatterPlot, showLegend::Bool)
	global dgCount

	for x in 1:length(g.dataGroups)
		dg = g.dataGroups[x]
		if dg.lineColor == "blank"
			dg.lineColor = goodColors[mod(x, length(goodColors))]
		end

		if dg.markerColor == "blank"
			dg.markerColor = goodColors[mod(x, length(goodColors))]
		end

		if dg.markerType == "blank"
			dg.markerType = goodSymbols[mod(x, length(goodSymbols))]
		end

		plotPoints(p, dg)

		if showLegend
			dgCount += 1
			if dg.legend.label == "Data Group "
				dg.legend.label = dg.legend.label*string(dgCount)
			end
		end
	end
end

function scatterplot(dgs::Vector{DataGroup})
	plot(PlotFrame(ScatterPlot(dgs)))
end

function scatterplot(dg::DataGroup)
	plot(PlotFrame(ScatterPlot(dg)))
end

function scatterplot(dfs::Vector{DataFrame})
	dgs = DataGroup[]
	for df in dfs
		push!(dgs, DataGroup(df))
	end
	scatterplot(dgs)
end

function scatterplot(df::DataFrame)
	scatterplot(DataGroup(df))
end

function scatterplot(ys::Array{Float64, 2})
	dfs = DataFrame[]
	for r = 1:size(y, 1)
		push!(dfs, DataFrame(x = collect(1:length(y)), y = y[r,:]))
	end
	scatterplot(dfs)
end

function scatterplot(xs::Array{Float64, 2}, ys::Array{Float64, 2})
	dfs = DataFrame[]
	for r = 1:size(y, 1)
		push!(dfs, DataFrame(x = x[r,:], y = y[r,:]))
	end
	scatterplot(dfs)
end

function scatterplot(xs::Array{Float64, 2}, ys::Array{Float64, 2}, yerrs::Array{Float64, 2})
	dfs = DataFrame[]
	for r = 1:size(ys, 1)
		push!(dfs, DataFrame(x = vec(xs[r,:]), y = vec(ys[r,:]), yerr = vec(yerrs[r,:])))
	end
	scatterplot(dfs)
end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{typeof(SentinelArrays.defaultsentinel),Type{T} where T})
    Base.precompile(Tuple{typeof(empty),Dict{Any,Any},Type{Symbol},Type{SentinelArray{Int64,1,Int64,Missing,Array{Int64,1}}}})
    Base.precompile(Tuple{typeof(empty),Dict{Any,Any},Type{Symbol},Type{SentinelArray{Int64,1,Int64,Missing,ChainedVector{Int64,Array{Int64,1}}}}})
    Base.precompile(Tuple{typeof(empty),Dict{Symbol,SentinelArray{Int64,1,Int64,Missing,ChainedVector{Int64,Array{Int64,1}}}},Type{Symbol},Type{SentinelArray{T,1,S,Missing,A} where A<:AbstractArray{T,1} where S where T}})
    Base.precompile(Tuple{typeof(merge!),Dict{Symbol,SentinelArray{T,1,S,Missing,A} where A<:AbstractArray{T,1} where S where T},Dict{Symbol,SentinelArray{Int64,1,Int64,Missing,Array{Int64,1}}}})
    Base.precompile(Tuple{typeof(merge!),Dict{Symbol,SentinelArray{T,1,S,Missing,A} where A<:AbstractArray{T,1} where S where T},Dict{Symbol,SentinelArray{Int64,1,Int64,Missing,ChainedVector{Int64,Array{Int64,1}}}}})
    Base.precompile(Tuple{typeof(resize!),SentinelArray{Dates.Date,1,Dates.Date,Missing,Array{Dates.Date,1}},Int64})
    Base.precompile(Tuple{typeof(resize!),SentinelArray{Dates.DateTime,1,Dates.DateTime,Missing,Array{Dates.DateTime,1}},Int64})
    Base.precompile(Tuple{typeof(resize!),SentinelArray{Dates.Time,1,Dates.Time,Missing,Array{Dates.Time,1}},Int64})
    Base.precompile(Tuple{typeof(resize!),SentinelArray{Float64,1,Float64,Missing,Array{Float64,1}},Int64})
    Base.precompile(Tuple{typeof(resize!),SentinelArray{Int64,1,Int64,Missing,Array{Int64,1}},Int64})
    Base.precompile(Tuple{typeof(setindex!),Dict{Symbol,AbstractArray{T,1} where T},ChainedVector{Union{Missing, Bool},Array{Union{Missing, Bool},1}},Symbol})
    Base.precompile(Tuple{typeof(setindex!),Dict{Symbol,AbstractArray{T,1} where T},SentinelArray{Dates.Date,1,Dates.Date,Missing,Array{Dates.Date,1}},Symbol})
    Base.precompile(Tuple{typeof(setindex!),Dict{Symbol,AbstractArray{T,1} where T},SentinelArray{Dates.Date,1,Dates.Date,Missing,ChainedVector{Dates.Date,Array{Dates.Date,1}}},Symbol})
    Base.precompile(Tuple{typeof(setindex!),Dict{Symbol,AbstractArray{T,1} where T},SentinelArray{Dates.DateTime,1,Dates.DateTime,Missing,Array{Dates.DateTime,1}},Symbol})
    Base.precompile(Tuple{typeof(setindex!),Dict{Symbol,AbstractArray{T,1} where T},SentinelArray{Dates.DateTime,1,Dates.DateTime,Missing,ChainedVector{Dates.DateTime,Array{Dates.DateTime,1}}},Symbol})
    Base.precompile(Tuple{typeof(setindex!),Dict{Symbol,AbstractArray{T,1} where T},SentinelArray{Dates.Time,1,Dates.Time,Missing,Array{Dates.Time,1}},Symbol})
    Base.precompile(Tuple{typeof(setindex!),Dict{Symbol,AbstractArray{T,1} where T},SentinelArray{Dates.Time,1,Dates.Time,Missing,ChainedVector{Dates.Time,Array{Dates.Time,1}}},Symbol})
    Base.precompile(Tuple{typeof(setindex!),Dict{Symbol,AbstractArray{T,1} where T},SentinelArray{Float64,1,Float64,Missing,Array{Float64,1}},Symbol})
    Base.precompile(Tuple{typeof(setindex!),Dict{Symbol,AbstractArray{T,1} where T},SentinelArray{Float64,1,Float64,Missing,ChainedVector{Float64,Array{Float64,1}}},Symbol})
    Base.precompile(Tuple{typeof(setindex!),Dict{Symbol,AbstractArray{T,1} where T},SentinelArray{Int64,1,Int64,Missing,Array{Int64,1}},Symbol})
    Base.precompile(Tuple{typeof(setindex!),Dict{Symbol,AbstractArray{T,1} where T},SentinelArray{Int64,1,Int64,Missing,ChainedVector{Int64,Array{Int64,1}}},Symbol})
    Base.precompile(Tuple{typeof(setindex!),Dict{Symbol,AbstractArray{T,1} where T},SentinelArray{String,1,UndefInitializer,Missing,Array{String,1}},Symbol})
    Base.precompile(Tuple{typeof(setindex!),Dict{Symbol,AbstractArray{T,1} where T},SentinelArray{String,1,UndefInitializer,Missing,ChainedVector{String,Array{String,1}}},Symbol})
    Base.precompile(Tuple{typeof(setindex!),Dict{Symbol,SentinelArray{Int64,1,Int64,Missing,Array{Int64,1}}},SentinelArray{Int64,1,Int64,Missing,Array{Int64,1}},Symbol})
    Base.precompile(Tuple{typeof(setindex!),Dict{Symbol,SentinelArray{Int64,1,Int64,Missing,ChainedVector{Int64,Array{Int64,1}}}},SentinelArray{Int64,1,Int64,Missing,ChainedVector{Int64,Array{Int64,1}}},Symbol})
    Base.precompile(Tuple{typeof(setindex!),Dict{Symbol,SentinelArray{T,1,S,Missing,A} where A<:AbstractArray{T,1} where S where T},SentinelArray{Float64,1,Float64,Missing,Array{Float64,1}},Symbol})
    Base.precompile(Tuple{typeof(setindex!),Dict{Symbol,SentinelArray{T,1,S,Missing,A} where A<:AbstractArray{T,1} where S where T},SentinelArray{Float64,1,Float64,Missing,ChainedVector{Float64,Array{Float64,1}}},Symbol})
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 727ca485-8350-4adc-9aa0-59bbb84a9205
using Manopt, Manifolds, Random, Colors, PlutoUI

# ╔═╡ 6bf76330-ad0e-11ec-0c00-894872624127
md"""
# Get started: Optimize!

This example illustrates how to set up and solve optimization problems and how
to further get data from the algorithm using debug output and record data.
We will use the Riemannian mean and median as simple examples.

To start from the quite general case: A __Solver__ is an algorithm that aims
to solve

```math
\operatorname*{argmin}_{x∈\mathcal M} f(x)
```

where ``\mathcal M`` is a [Manifold](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.Manifold) and ``f:\mathcal M → ℝ`` is the cost function.

In `Manopt.jl` a __Solver__ is an algorithm that requires a [`Problem`](https://manoptjl.org/stable/plans/index.html#Manopt.Problem)
`p` and [`Options`](https://manoptjl.org/stable/plans/index.html#Manopt.Options) `o`. While former contains __static__ data,
most prominently the manifold ``\mathcal M`` (usually as `p.M`) and the cost
function ``f`` (usually as `x->get_cost(p, x)`), the latter contains __dynamic__
data, i.e. things that usually change during the algorithm, are allowed to
change, or specify the details of the algorithm to use. Together they form a
__plan__. A __plan__ uniquely determines the algorithm to use and provide all
necessary information to run the algorithm.
"""

# ╔═╡ 94dee66e-2f37-4cc0-8451-c0bbb5eae2c9
md"""
# Setup

Let‘s first set up a few variables
"""

# ╔═╡ 4235a1ba-3cf2-49dc-9a26-32fafc7a7008
begin
    localpath = join(splitpath(@__FILE__)[1:(end - 1)], "/") # files folder
    image_prefix = localpath * "/optimize"
    @info image_prefix
    render_asy = false # on CI or when you do not have asymptote, this should be false
end;

# ╔═╡ 65489070-9066-46bb-b5b1-52732dbe9bc7
md"""
# Example

A gradient plan consists of a [`GradientProblem`](https://manoptjl.org/stable/plans/index.html#Manopt.GradientProblem) with the fields `M`,
cost function ``f`` as well as `gradient` storing the gradient function
 corresponding to ``f``. Accessing both functions can be done directly but should
be encapsulated using [`get_cost`](https://manoptjl.org/stable/plans/index.html#Manopt.get_cost)`(p,x)` and [`get_gradient`](https://manoptjl.org/stable/plans/index.html#Manopt.get_gradient)`(p,x)`,
where in both cases `x` is a point on the [`Manifold`](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.Manifold) `M`.
Second, the [`GradientDescentOptions`](https://manoptjl.org/stable/solvers/gradient_descent.html#Manopt.GradientDescentOptions) specify that the algorithm to solve
the `GradientProblem` is the [gradient
descent](https://en.wikipedia.org/wiki/Gradient_descent) algorithm.
It requires an initial value `o.x0`, a `StoppingCriterion` `o.stop`, a
`Stepsize` `o.stepsize` and a retraction `o.retraction`.
Internally is stores the last evaluation of the gradient at `o.gradient` for convenience.
The only mandatory parameter is the initial value `x0`, though the defaults for
both the stopping criterion ([`StopAfterIteration`](@ref)`(100)`) as well as the
stepsize ([`ConstantStepsize`](@ref)`(1.)` are quite conservative, but are
chosen to be as simple as possible.

With these two at hand, running the algorithm just requires to call `x_opt = solve(p,o)`.

In the following two examples we will see, how to use a higher level interface
that allows to more easily activate for example a debug output or record values during the iterations.
"""

# ╔═╡ 177cc292-94d3-4344-857e-30483f592a55
md"""
Let‘s load a few colors from [Paul Tol](https://personal.sron.nl/~pault/)
"""

# ╔═╡ 0b405c42-19a5-480d-b1dc-0fb8811a48fa
begin
    black = RGBA{Float64}(colorant"#000000")
    TolVibrantOrange = RGBA{Float64}(colorant"#EE7733")
    TolVibrantBlue = RGBA{Float64}(colorant"#0077BB")
    TolVibrantTeal = RGBA{Float64}(colorant"#009988")
    TolVibrantMagenta = RGBA{Float64}(colorant"#EE3377")
end;

# ╔═╡ 803fc640-bbed-4700-8a1e-f414c6446eea
md"""
# The data set
"""

# ╔═╡ 278dbd3c-aa91-4d9d-ad49-c3b4b336efe2
md"""
We take a look at a srandom set of points.
"""

# ╔═╡ 7e0301fb-7465-410c-b47c-04686bf44ab1
begin
    n = 100
    σ = π / 8
    M = Sphere(2)
    x = 1 / sqrt(2) * [1.0, 0.0, 1.0]
    Random.seed!(42)
    data = [exp(M, x, random_tangent(M, x, Val(:Gaussian), σ)) for i in 1:n]
end

# ╔═╡ 9b130a57-293d-429d-88b5-78bfacbf836f
asymptote_export_S2_signals(
    image_prefix * "/startDataAndCenter.asy";
    points=[[x], data],
    colors=Dict(:points => [TolVibrantBlue, TolVibrantTeal]),
    dot_size=3.5,
    camera_position=(1.0, 0.5, 0.5),
)

# ╔═╡ d2bc2e63-6ae7-4f54-b176-f74e66365b1a
render_asy && render_asymptote(image_prefix * "/startDataAndCenter.asy"; render=2)

# ╔═╡ 6ea36d50-96f0-46e1-8d90-529f0b23120d
PlutoUI.LocalResource(image_prefix * "/startDataAndCenter.png")

# ╔═╡ e21d6d03-4c61-457b-a9c7-fad5b4f369db
md"""
## Computing the Mean

To compute the mean on the manifold we use the characterization, that the
Euclidean mean minimizes the sum of squared distances, and end up with the
following cost function. Its minimizer is called
[Riemannian Center of Mass](https://arxiv.org/abs/1407.2087).

> **Note.**
> There are more sophisticated methods tailored for the specific manifolds available in [Manifolds.jl](https://juliamanifolds.github.io/Manifolds.jl/) see [mean](https://juliamanifolds.github.io/Manifolds.jl/stable/features/statistics.html#Statistics.mean-Tuple{Manifold,AbstractArray{T,1}%20where%20T,AbstractArray{T,1}%20where%20T,ExtrinsicEstimation}).
"""

# ╔═╡ 515eaa1d-8307-45ba-ae63-be070dc2ff1c
F(M, y) = sum(1 / (2 * n) * distance.(Ref(M), Ref(y), data) .^ 2)

# ╔═╡ 00bc7b5b-7cc9-43bc-bf96-b5d2b85ddb9a
gradF(M, y) = sum(1 / n * grad_distance.(Ref(M), data, Ref(y)))

# ╔═╡ 3cc7ce2b-8ae6-4ba9-8298-2050fe7081eb
md"""
Note that the [grad_distance](https://manoptjl.org/stable/functions/gradients.html#Manopt.grad_distance) defaults to the case `p=2`, i.e. the
gradient of the squared distance. For details on convergence of the gradient
descent for this problem, see [^AfsariTronVidal2013]

The easiest way to call the gradient descent is now to call
[gradient_descent](https://manoptjl.org/stable/solvers/gradient_descent.html#Manopt.gradient_descent).
"""

# ╔═╡ c067f4d1-b54b-4228-85cf-70cdbbdc948b
gradient_descent(M, F, gradF, data[1])

# ╔═╡ 7aee0626-14ba-431e-af84-79c2dfc021da
md"""
In order to get more details, we further add the `debug=` options, which
act as a [decorator pattern](https://en.wikipedia.org/wiki/Decorator_pattern).

The following debug prints

```juliaREPL
# i | Last Change: | F(x): | x:
```
as well as the reason why the algorithm stopped at the end.

The elements passed to `debug=` are postprocessed, there are specifiic symbols and formats we can use. For example `:Iteration` to plot the iteratin, but we can also add a format for the print using `(:Change, "| Last change: %e3.8")`, i.e. the second string is a format for printf.

Note that here we use `PlutoUI` to see the output also here in the notebook
"""

# ╔═╡ 51daafb5-84a0-47fd-ac40-4d53888cd914
with_terminal() do
    global xMean = gradient_descent(
        M,
        F,
        gradF,
        data[1];
        debug=[
            :Iteration,
            (:Change, "change: %1.9f | "),
            (:Cost, " F(x): %1.11f | "),
            "\n",
            :Stop,
        ],
    )
end

# ╔═╡ 863bf8b8-272c-40d6-985f-0a7cf9454756
md"""
A way to get better performance and for convex and coercive costs a guaranteed convergence is to switch the default
[`ConstantStepsize`](@ref)(1.0) with a step size that performs better, for
example the [`ArmijoLinesearch`](https://manoptjl.org/stable/plans/index.html#Manopt.ArmijoLinesearch).
We can tweak the default values for the `contractionFactor` and the `sufficientDecrease`  beyond constant step size which is already quite fast. We get
"""

# ╔═╡ 38df2fb3-f742-4652-857c-baa403985ff8
with_terminal() do
    global xMean2 = gradient_descent(
        M,
        F,
        gradF,
        data[1];
        stepsize=ArmijoLinesearch(1.0, ExponentialRetraction(), 0.99, 0.5),
        debug=[
            :Iteration, (:Change, "change: %1.9f | "), :Cost, (:x, " | %s"), "\n", :Stop
        ],
    )
end

# ╔═╡ 244eb6ea-0bdb-4443-8dc1-40419966198a
md"""
which finishes in 5 steaps, just slightly better than the previous computation.
"""

# ╔═╡ eedfedbc-f305-48ff-8aa7-78b6aa6c4d02
F(M, xMean) - F(M, xMean2)

# ╔═╡ 887bd166-c278-48d3-bfb1-31c7026433b1
md"""
Note that other optimization tasks may have other speedup opportunities.

For even more precision, we can further require a smaller gradient norm.
This is done by changing the `StoppingCriterion` used, where several
criteria can be combined using `&` and/or `|`.  If we want to decrease the final
gradient (from less that 1e-8) norm but keep the maximal number of iterations
to be 200, we can run
"""

# ╔═╡ a99a5603-6ef5-43e8-a082-54dd20226956
with_terminal() do
    global xMean3 = gradient_descent(
        M,
        F,
        gradF,
        data[1];
        stepsize=ArmijoLinesearch(1.0, ExponentialRetraction(), 0.99, 0.5),
        stopping_criterion=StopAfterIteration(200) | StopWhenGradientNormLess(1e-15),
        debug=[
            :Iteration,
            (:Change, "change: %1.9f | "),
            (:Cost, "F(x): %1.9f"),
            (:x, " | %s"),
            "\n",
            :Stop,
        ],
    )
end

# ╔═╡ 98028747-31dd-4bf8-b4b5-0959d5afb75c
md"""
Let‘s add this point to out data image
"""

# ╔═╡ fb07943f-54b4-4cb3-b1fd-f3ab06b4d033
asymptote_export_S2_signals(
    image_prefix * "/startDataCenterMean.asy";
    points=[[x], data, [xMean3]],
    colors=Dict(:points => [TolVibrantBlue, TolVibrantTeal, TolVibrantOrange]),
    dot_size=3.5,
    camera_position=(1.0, 0.5, 0.5),
)

# ╔═╡ 5e1e6db2-39da-4857-8745-eda6ae510fa8
render_asy && render_asymptote(image_prefix * "/startDataCenterMean.asy"; render=2)

# ╔═╡ c9e09455-1af2-40aa-aa05-7fa329b5eec7
PlutoUI.LocalResource(image_prefix * "/startDataCenterMean.png")

# ╔═╡ a17ddd91-0839-48fe-ab78-445226ee4ff9
md"""
## Computing the Median

> **Note.**
> There are more sophisticated methods tailored for the specific manifolds available in
> [Manifolds.jl](https://juliamanifolds.github.io/Manifolds.jl/) see [`median`](https://juliamanifolds.github.io/Manifolds.jl/stable/features/statistics.html#Statistics.median-Tuple{Manifold,AbstractArray{T,1}%20where%20T,AbstractArray{T,1}%20where%20T,CyclicProximalPointEstimation}).

Similar to the mean you can also define the median as the minimizer of the
distances, see for example [^Bačák2014], but since
this problem is not differentiable, we employ the Cyclic Proximal Point (CPP)
algorithm, described in the same reference. We define
"""

# ╔═╡ 439b700d-dce7-43c4-bb1e-1f263a3f54a4
F2(M, y) = sum(1 / (2 * n) * distance.(Ref(M), Ref(y), data))

# ╔═╡ 79921477-db82-424f-8f8e-9bade1cd63c3
proxes = Function[(M, λ, y) -> prox_distance(M, λ / n, di, y, 1) for di in data]

# ╔═╡ dc236508-283c-4047-a216-c570e35bc791
md"""
So we call the cyclic proximal point algorithm this time with a recording and activate the return of the complete options to access the recorded values. We further increase the display of the cost function to more digits.
"""

# ╔═╡ d2a8250e-7796-454b-a0bf-9970b1b9a2aa
with_terminal() do
    global o = cyclic_proximal_point(
        M,
        F2,
        proxes,
        data[1];
        debug=[
            :Iteration,
            " | ",
            :Change,
            " | ",
            (:Cost, "F(x): %1.12f"),
            " | ",
            :x,
            "\n",
            50,
            :Stop,
        ],
        record=[:Iteration, :Change, :Cost, :x],
        return_options=true,
    )
end

# ╔═╡ 90414c9b-b164-490c-b939-2472d8284887
xMedian = get_solver_result(o)

# ╔═╡ aa88f273-2a01-4d93-857b-8b9c19fcdd1a
md"""
where the differences to `gradient_descent` are as follows

* the third parameter is now an Array of proximal maps
* debug is reduces to only every 50th iteration
* we further activated a `RecordAction` using the `record=` optional
  parameter. These work very similar to those in debug, but they
  collect their data in an array. The high level interface then returns two
  variables; the `values` do contain an array of recorded
  datum per iteration. Here a Tuple containing the iteration, last change and
  cost respectively; see [RedordOptions](https://manoptjl.org/stable/plans/index.html#RecordOptions-1) for details.

We can access the recorded values using `get_record` and contains of a tuple per iteration and contains the iteration number, the change and the cost.
"""

# ╔═╡ c835b5ec-085e-4c9d-b777-76036515bcd1
values = get_record(o)

# ╔═╡ 9a19d3bb-4487-4469-a856-b1a0f3f540a7
asymptote_export_S2_signals(
    image_prefix * "/startDataCenterMedianAndMean.asy";
    points=[[x], data, [xMean], [xMedian]],
    colors=Dict(
        :points => [TolVibrantBlue, TolVibrantTeal, TolVibrantOrange, TolVibrantMagenta]
    ),
    dot_size=3.5,
    camera_position=(1.0, 0.5, 0.5),
)

# ╔═╡ e6e8cc14-a67b-4e6a-a83b-73f9ba1e1dcb
render_asy && render_asymptote(image_prefix * "/startDataCenterMedianAndMean.asy"; render=2)

# ╔═╡ 6cd0aaab-e0bf-428e-a84c-f456a32f8e36
md"""
In the following image the mean (orange), median (magenta) are shown.
"""

# ╔═╡ 97e1bec2-3577-423a-a012-ca1b5413b29a
PlutoUI.LocalResource(image_prefix * "/startDataCenterMedianAndMean.png")

# ╔═╡ 0ba01385-6502-4e66-afee-4c4af391c9b9
md"""
## Literature

[^Bačák2014]:
	> Bačák, M: __Computing Medians and Means in Hadamard Spaces.__
    > SIAM Journal on Optimization, Volume 24, Number 3, pp. 1542–1566,
    > doi: [10.1137/140953393](https://doi.org/10.1137/140953393),
    > arxiv: [1210.2145](https://arxiv.org/abs/1210.2145)
[^AfsariTronVidal2013]:
	> Afsari, B; Tron, R.; Vidal, R.: __On the Convergence of Gradient Descent for Finding the Riemannian Center of Mass__,
	> SIAM Journal on Control and Optimization, Volume 51, Issue 3, pp. 2230–2260.
	> doi: [10.1137/12086282X](https://doi.org/10.1137/12086282X),
	> arxiv: [1201.0925](https://arxiv.org/abs/1201.0925).
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
Manopt = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[compat]
Colors = "~0.12.8"
Manifolds = "~0.7.7"
Manopt = "~0.3.20"
PlutoUI = "~0.7.37"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.2"
manifest_format = "2.0"

[[deps.AbstractFFTs]]
deps = ["ChainRulesCore", "LinearAlgebra"]
git-tree-sha1 = "6f1d9bc1c08f9f4a8fa92e3ea3cb50153a1b40d4"
uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
version = "1.1.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.3"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.ArnoldiMethod]]
deps = ["LinearAlgebra", "Random", "StaticArrays"]
git-tree-sha1 = "62e51b39331de8911e4a7ff6f5aaf38a5f4cc0ae"
uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
version = "0.2.0"

[[deps.ArrayInterface]]
deps = ["Compat", "IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
git-tree-sha1 = "6e8fada11bb015ecf9263f64b156f98b546918c7"
uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
version = "5.0.5"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.Calculus]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f641eb0a4f00c343bbc32346e1217b86f3ce9dad"
uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
version = "0.5.1"

[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "9950387274246d08af38f6eef8cb5480862a435f"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.14.0"

[[deps.ChangesOfVariables]]
deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.2"

[[deps.ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random"]
git-tree-sha1 = "12fc73e5e0af68ad3137b886e3f7c1eacfca2640"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.17.1"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "96b0bc6c52df76506efc8a441c6cf1adcb1babc4"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.42.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.CovarianceEstimation]]
deps = ["LinearAlgebra", "Statistics", "StatsBase"]
git-tree-sha1 = "a3e070133acab996660d31dcf479ea42849e368f"
uuid = "587fd27a-f159-11e8-2dae-1979310e6154"
version = "0.2.7"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.DensityInterface]]
deps = ["InverseFunctions", "Test"]
git-tree-sha1 = "80c3e8639e3353e5d2912fb3a1916b8455e2494b"
uuid = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
version = "0.4.0"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.Distributions]]
deps = ["ChainRulesCore", "DensityInterface", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "Test"]
git-tree-sha1 = "c43e992f186abaf9965cc45e372f4693b7754b22"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.25.52"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.DualNumbers]]
deps = ["Calculus", "NaNMath", "SpecialFunctions"]
git-tree-sha1 = "5837a837389fccf076445fce071c8ddaea35a566"
uuid = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
version = "0.6.8"

[[deps.Einsum]]
deps = ["Compat"]
git-tree-sha1 = "4a6b3eee0161c89700b6c1949feae8b851da5494"
uuid = "b7d42ee7-0b51-5a75-98ca-779d3107e4c0"
version = "0.4.1"

[[deps.FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "246621d23d1f43e3b9c368bf3b72b2331a27c286"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.13.2"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Graphs]]
deps = ["ArnoldiMethod", "Compat", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
git-tree-sha1 = "57c021de207e234108a6f1454003120a1bf350c4"
uuid = "86223c79-3864-5bf0-83f7-82e725a168b6"
version = "1.6.0"

[[deps.HybridArrays]]
deps = ["LinearAlgebra", "Requires", "StaticArrays"]
git-tree-sha1 = "eb6b23460f5544c5d09efae0818b86736cefcd3d"
uuid = "1baab800-613f-4b0a-84e4-9cd3431bfbb9"
version = "0.4.10"

[[deps.HypergeometricFunctions]]
deps = ["DualNumbers", "LinearAlgebra", "SpecialFunctions", "Test"]
git-tree-sha1 = "65e4589030ef3c44d3b90bdc5aac462b4bb05567"
uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
version = "0.3.8"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.IfElse]]
git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1"
uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
version = "0.1.1"

[[deps.Inflate]]
git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c"
uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
version = "0.1.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "91b5dcf362c5add98049e6c29ee756910b03051d"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.3"

[[deps.IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.4.1"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

[[deps.Kronecker]]
deps = ["LinearAlgebra", "NamedDims", "SparseArrays", "StatsBase"]
git-tree-sha1 = "a51f46415c844dee694cb8b20a3fcbe6dba342c2"
uuid = "2c470bb0-bcc8-11e8-3dad-c9649493f05e"
version = "0.5.1"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "58f25e56b706f95125dcb796f39e1fb01d913a71"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.10"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

[[deps.Manifolds]]
deps = ["Colors", "Distributions", "Einsum", "Graphs", "HybridArrays", "Kronecker", "LinearAlgebra", "ManifoldsBase", "Markdown", "Random", "RecipesBase", "RecursiveArrayTools", "Requires", "SimpleWeightedGraphs", "SpecialFunctions", "StaticArrays", "Statistics", "StatsBase"]
git-tree-sha1 = "e8fb7c69d6e67d377152fd882a20334535db050f"
uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
version = "0.7.7"

[[deps.ManifoldsBase]]
deps = ["LinearAlgebra", "Markdown"]
git-tree-sha1 = "058ba95cf4a41d4c7b88879f5b961352880ec919"
uuid = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb"
version = "0.12.12"

[[deps.Manopt]]
deps = ["ColorSchemes", "ColorTypes", "Colors", "DataStructures", "Dates", "LinearAlgebra", "ManifoldsBase", "Markdown", "Printf", "Random", "Requires", "SparseArrays", "StaticArrays", "Statistics", "Test"]
git-tree-sha1 = "bcd147170ad2699518155ba9b12ad33f835c5a4d"
uuid = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5"
version = "0.3.20"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NaNMath]]
git-tree-sha1 = "737a5957f387b17e74d4ad2f440eb330b39a62c5"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "1.0.0"

[[deps.NamedDims]]
deps = ["AbstractFFTs", "ChainRulesCore", "CovarianceEstimation", "LinearAlgebra", "Pkg", "Requires", "Statistics"]
git-tree-sha1 = "0856b62716585eb90cc1dada226ac9eab5f69aa5"
uuid = "356022a1-0364-5f58-8944-0da4b18d706f"
version = "0.2.47"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[deps.OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.PDMats]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
git-tree-sha1 = "e8185b83b9fc56eb6456200e873ce598ebc7f262"
uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.11.7"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "85b5da0fa43588c75bb1ff986493443f821c70b7"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.3"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "bf0a1121af131d9974241ba53f601211e9303a9e"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.37"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "d3538e7f8a790dc8903519090857ef8e1283eecd"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.5"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "78aadffb3efd2155af139781b8a8df1ef279ea39"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.2"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.RecipesBase]]
git-tree-sha1 = "6bf3f380ff52ce0832ddd3a2a7b9538ed1bcca7d"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.2.1"

[[deps.RecursiveArrayTools]]
deps = ["Adapt", "ArrayInterface", "ChainRulesCore", "DocStringExtensions", "FillArrays", "LinearAlgebra", "RecipesBase", "Requires", "StaticArrays", "Statistics", "ZygoteRules"]
git-tree-sha1 = "f5dd036acee4462949cc10c55544cc2bee2545d6"
uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
version = "2.25.1"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

[[deps.Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
version = "0.7.0"

[[deps.Rmath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
version = "0.3.0+0"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.SimpleTraits]]
deps = ["InteractiveUtils", "MacroTools"]
git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231"
uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
version = "0.9.4"

[[deps.SimpleWeightedGraphs]]
deps = ["Graphs", "LinearAlgebra", "Markdown", "SparseArrays", "Test"]
git-tree-sha1 = "a6f404cc44d3d3b28c793ec0eb59af709d827e4e"
uuid = "47aef6b3-ad0c-573a-a1e2-d07658019622"
version = "1.2.1"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "5ba658aeecaaf96923dce0da9e703bd1fe7666f9"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "2.1.4"

[[deps.Static]]
deps = ["IfElse"]
git-tree-sha1 = "87e9954dfa33fd145694e42337bdd3d5b07021a6"
uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
version = "0.6.0"

[[deps.StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "4f6ec5d99a28e1a749559ef7dd518663c5eca3d5"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.4.3"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.StatsAPI]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "c3d8ba7f3fa0625b062b82853a7d5229cb728b6b"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.2.1"

[[deps.StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "8977b17906b0a1cc74ab2e3a05faa16cf08a8291"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.16"

[[deps.StatsFuns]]
deps = ["ChainRulesCore", "HypergeometricFunctions", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "25405d7016a47cf2bd6cd91e66f4de437fd54a07"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.16"

[[deps.SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.ZygoteRules]]
deps = ["MacroTools"]
git-tree-sha1 = "8c1a8e4dfacb1fd631745552c8db35d0deb09ea0"
uuid = "700de1a5-db45-46bc-99cf-38207098b444"
version = "0.2.2"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─6bf76330-ad0e-11ec-0c00-894872624127
# ╟─94dee66e-2f37-4cc0-8451-c0bbb5eae2c9
# ╠═4235a1ba-3cf2-49dc-9a26-32fafc7a7008
# ╟─65489070-9066-46bb-b5b1-52732dbe9bc7
# ╠═727ca485-8350-4adc-9aa0-59bbb84a9205
# ╟─177cc292-94d3-4344-857e-30483f592a55
# ╠═0b405c42-19a5-480d-b1dc-0fb8811a48fa
# ╟─803fc640-bbed-4700-8a1e-f414c6446eea
# ╟─278dbd3c-aa91-4d9d-ad49-c3b4b336efe2
# ╠═7e0301fb-7465-410c-b47c-04686bf44ab1
# ╠═9b130a57-293d-429d-88b5-78bfacbf836f
# ╠═d2bc2e63-6ae7-4f54-b176-f74e66365b1a
# ╟─6ea36d50-96f0-46e1-8d90-529f0b23120d
# ╟─e21d6d03-4c61-457b-a9c7-fad5b4f369db
# ╠═515eaa1d-8307-45ba-ae63-be070dc2ff1c
# ╠═00bc7b5b-7cc9-43bc-bf96-b5d2b85ddb9a
# ╟─3cc7ce2b-8ae6-4ba9-8298-2050fe7081eb
# ╠═c067f4d1-b54b-4228-85cf-70cdbbdc948b
# ╟─7aee0626-14ba-431e-af84-79c2dfc021da
# ╠═51daafb5-84a0-47fd-ac40-4d53888cd914
# ╟─863bf8b8-272c-40d6-985f-0a7cf9454756
# ╠═38df2fb3-f742-4652-857c-baa403985ff8
# ╟─244eb6ea-0bdb-4443-8dc1-40419966198a
# ╠═eedfedbc-f305-48ff-8aa7-78b6aa6c4d02
# ╟─887bd166-c278-48d3-bfb1-31c7026433b1
# ╠═a99a5603-6ef5-43e8-a082-54dd20226956
# ╟─98028747-31dd-4bf8-b4b5-0959d5afb75c
# ╠═fb07943f-54b4-4cb3-b1fd-f3ab06b4d033
# ╠═5e1e6db2-39da-4857-8745-eda6ae510fa8
# ╟─c9e09455-1af2-40aa-aa05-7fa329b5eec7
# ╟─a17ddd91-0839-48fe-ab78-445226ee4ff9
# ╠═439b700d-dce7-43c4-bb1e-1f263a3f54a4
# ╠═79921477-db82-424f-8f8e-9bade1cd63c3
# ╟─dc236508-283c-4047-a216-c570e35bc791
# ╠═d2a8250e-7796-454b-a0bf-9970b1b9a2aa
# ╠═90414c9b-b164-490c-b939-2472d8284887
# ╟─aa88f273-2a01-4d93-857b-8b9c19fcdd1a
# ╠═c835b5ec-085e-4c9d-b777-76036515bcd1
# ╠═9a19d3bb-4487-4469-a856-b1a0f3f540a7
# ╠═e6e8cc14-a67b-4e6a-a83b-73f9ba1e1dcb
# ╠═6cd0aaab-e0bf-428e-a84c-f456a32f8e36
# ╟─97e1bec2-3577-423a-a012-ca1b5413b29a
# ╟─0ba01385-6502-4e66-afee-4c4af391c9b9
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
include("tools.jl")
using tools
using Base.Test

@testset "HOSVD tools" begin
    @testset "Matricisations" begin
        a = collect(1:8);
        a = reshape(a, (2,2,2));

        @test unfolding(a, 1) == [1 3 5 7;
                                 2 4 6 8]

        @test unfolding(a, 2) == [1 2 5 6;
                                 3 4 7 8]

        @test unfolding(a, 3) == [1 2 3 4;
                                 5 6 7 8]

        @test a == folding(unfolding(a,1),1,size(a))
        @test a == folding(unfolding(a,2),2,size(a))
        @test a == folding(unfolding(a,3),3,size(a))

    end

    @testset "Matricisation of function, mode $n" for n in 1:3
        # the unfolding for arrays is assumed to work now, as it
        # passed the tests earlier
        dims = (10,10,10);
        a = rhs_sin(10);
        b = get_rhs_sin(10);

        b_u = unfolding(b, n);
        a_u = unfolding_fun(a, n, dims);

        a_u_arr = [a_u(i,j) for i=1:dims[1],j=1:prod(dims[1:end .!= n])];

        @test a_u_arr ≈ b_u
    end

    @testset "HOSVD" begin
        # set the random number generator to a defined state for
        # reproducibility
        srand(0);
        b = rand(100,100,100);
        ϵ = 1e-4;
        c = hosvd(b, ϵ)

        # for the comparison need to multiply ϵ with norm of b, as the error
        # in this hosvd implementation is relative
        @test norm((fullten(c)-b)[:],2) <= ϵ * norm(b[:],2)

        # use the built in isapprox function, it chooses a default tolerance
        @test fullten(c) ≈ b
    end
end;


text:
"""
    γ = opnorm1(op)

Compute `γ`, the induced `1`-norm of the linear operator `op`, as the maximum of `1`-norm of the
columns of the associated `m x n` matrix ``M`` `= Matrix(op)`:
```math
\\gamma = \\|op\\|_1 := \\max_{1 ≤ j ≤ n} \\|M_j\\|_1
```
with ``M_j`` the `j`-th column of ``M``. This function is not recommended to be used for large order operators.

# Examples
```jldoctest
julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
3×3 Array{Float64,2}:
 -6.0  -2.0   1.0
  5.0   1.0  -1.0
 -4.0  -2.0  -1.0

julia> opnorm1(lyapop(A))
30.0

julia> opnorm1(invlyapop(A))
3.7666666666666706
```
"""
function opnorm1(op::AbstractLinearOperator)
  (m, n) = size(op)
  T = eltype(op)
  Tnorm = real(T)
  Tsum = promote_type(Float64, Tnorm)
  nrm::Tsum = 0
  for j = 1 : n
      ej = zeros(Tsum, n)
      ej[j] = 1
      try
         nrm = max(nrm,norm(op*ej,1))
      catch err
         # if isnothing(findfirst("SingularException",string(err))) &&
         #    isnothing(findfirst("LAPACKException",string(err)))
         findfirst("SingularException",string(err)) === nothing &&
         findfirst("LAPACKException",string(err)) === nothing ? rethrow() : (return Inf)
      end
  end
  return convert(Tnorm, nrm)
end
"""
    γ = opnorm1est(op)

Compute `γ`, a lower bound of the `1`-norm of the square linear operator `op`, using
reverse communication based computations to evaluate `op * x` and `op' * x`.
It is expected that in most cases ``\\gamma > \\|op\\|_1/10``, which is usually
acceptable for estimating the condition numbers of linear operators.

# Examples
```jldoctest
julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
3×3 Array{Float64,2}:
 -6.0  -2.0   1.0
  5.0   1.0  -1.0
 -4.0  -2.0  -1.0

julia> opnorm1est(lyapop(A))
18.0

julia> opnorm1est(invlyapop(A))
3.76666666666667
```
"""
function opnorm1est(op::AbstractLinearOperator)
  m, n = size(op)
  m == n || throw(DimensionMismatch("The operator op must be square"))
  T = promote_type(Float64,eltype(op))
  TR = real(T)
  TR == Float64 ? SMLNUM = reinterpret(Float64, 0x2000000000000000) : SMLNUM = reinterpret(Float32, 0x20000000)
  BIGNUM = 2*eps(TR) / SMLNUM
  cmplx = T<:Complex
  V = Array{T,1}(undef,n)
  X = Array{T,1}(undef,n)
  cmplx ? ISGN = Array{Int,1}(undef,1) : ISGN = Array{Int,1}(undef,n)
  ISAVE = Array{Int,1}(undef,3)
  ANORM = zero(TR)
  KASE = 0
  finish = false
  while !finish
     if cmplx
        ANORM, KASE = LapackUtil.lacn2!(V, X, ANORM, KASE, ISAVE )
     else
        ANORM, KASE = LapackUtil.lacn2!(V, X, ISGN, ANORM, KASE, ISAVE )
     end
     (isfinite(ANORM) && ANORM < BIGNUM) || (return Inf)
     if KASE != 0
        try
           KASE == 1 ? X = op*X : X = op'*X
        catch err
         #   if isnothing(findfirst("SingularException",string(err))) &&
         #      isnothing(findfirst("LAPACKException",string(err)))
           findfirst("SingularException",string(err)) === nothing &&
           findfirst("LAPACKException",string(err)) === nothing ? rethrow() : (return Inf)
        end
      else
        finish = true
     end
  end
  return ANORM
end

"""
    sep = opsepest(opinv; exact = false)

Compute `sep`, an estimation of the `1`-norm separation of a linear operator
`op`, where `opinv` is the inverse operator `inv(op)`. The estimate is computed as
``\\text{sep}  = 1 / \\|opinv\\|_1`` , using an estimate of the `1`-norm, if `exact = false`, or
the computed exact value of the `1`-norm, if `exact = true`.
The `exact = true` option is not recommended for large order operators.

The separation of the operator `op` is defined as
```math
\\text{sep} = \\displaystyle\\min_{X\\neq 0} \\frac{\\|op*X\\|}{\\|X\\|}.
```

An estimate of the reciprocal condition number of `op` can be computed as ``\\text{sep}/\\|op\\|_1``.

# Example
```jldoctest
julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
3×3 Array{Float64,2}:
 -6.0  -2.0   1.0
  5.0   1.0  -1.0
 -4.0  -2.0  -1.0

julia> opsepest(invlyapop(A))
0.26548672566371656

julia> 1/opnorm1est(invlyapop(A))
0.26548672566371656

julia> opsepest(invlyapop(A),exact = true)
0.26548672566371656

julia> 1/opnorm1(invlyapop(A))
0.26548672566371656
```
"""
function opsepest(opinv::AbstractLinearOperator; exact = false)
   T = promote_type(Float64,eltype(opinv))
   TR = real(T)
   TR == Float64 ? SMLNUM = reinterpret(Float64, 0x2000000000000000) : SMLNUM = reinterpret(Float32, 0x20000000)
   BIGNUM = 2*eps(TR) / SMLNUM
   ZERO = zero(TR)
   exact ? opinvnrm1 = opnorm1(opinv) : opinvnrm1 = opnorm1est(opinv)
   if opinvnrm1 >= BIGNUM
      return ZERO
   end
   return one(TR)/opinvnrm1
end
"""
    rcond = oprcondest(op, opinv; exact = false)

Compute `rcond`, an estimation of the `1`-norm reciprocal condition number
of a linear operator `op`, where `opinv` is the inverse operator `inv(op)`. The estimate is computed as
``\\text{rcond} = 1 / (\\|op\\|_1\\|opinv\\|_1)``, using estimates of the `1`-norm, if `exact = false`, or
computed exact values of the `1`-norm, if `exact = true`.
The `exact = true` option is not recommended for large order operators.

_Note:_ No check is performed to verify that `opinv = inv(op)`.

# Examples
```jldoctest
julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
3×3 Array{Float64,2}:
 -6.0  -2.0   1.0
  5.0   1.0  -1.0
 -4.0  -2.0  -1.0

julia> oprcondest(lyapop(A),invlyapop(A))
0.014749262536873142

julia> 1/opnorm1est(lyapop(A))/opnorm1est(invlyapop(A))
0.014749262536873142

julia> oprcondest(lyapop(A),invlyapop(A),exact = true)
0.008849557522123885

julia> 1/opnorm1(lyapop(A))/opnorm1(invlyapop(A))
0.008849557522123885
```
"""
function oprcondest(op:: LinearOperator, opinv::LinearOperator; exact = false)
   return opsepest(op, exact = exact)*opsepest(opinv, exact = exact)
end