text:
struct Beam
    λ::Real; z::Real; n::Real
    x::Real; k::Real; q::Complex
end
Beam(λ::Real, w0::Real, n0::Real=1) = Beam(λ, 0, n0, 0, 0, 1im*π*n0*w0^2/λ)

η(e::Interface) = e.η
η(e::Union{Tan,Sag}) = dη(e.e)
η(e::Element) = 1

dz(e::FreeSpace) = e.L
dz(e::Union{Tan,Sag}) = dz(e.e)
dz(e::Element) = 0

function transform(e::Element, Γ::Beam)
    M = RTM(e)
    return Beam(Γ.λ, Γ.z+dz(e), Γ.n/η(e), M*[Γ.x,Γ.k]..., /(M*[Γ.q,1]...))
end

function beamtrace(elems::Vector{<:Element}, Γ0::Beam)
    Γs = Vector{Beam}(undef, length(elems)+1); Γs[1] = Γ0
    for (ind, elem) in enumerate(elems)
        Γs[ind+1] = transform(elem, Γs[ind])
    end
    return Γs
end

spotsize(Γ::Beam) = /(-Γ.λ, π*Γ.n*imag(1/Γ.q)) |> sqrt
location(Γ::Beam) = Γ.z

"""

Return a function that calculates the beam's waist radius as a
function of the propagation distance.

Expects an optical system (e.g. a vector or collection of
[`Element`](@ref)) as first argument, and a [`Beam`](@ref) as second
argument. The optional keyword argument `outside` can be used to
define a return value for beam positions outside those covered by
the system; the default is to throw a DomainError.

"""
function waistradiusfunc(elements, beam::Beam; outside=nothing)
    beams = Vector{Beam}(undef, length(elements)+1)
    beams[1] = beam
    for (i, el) in enumerate(elements)
        beam = transform(el, beam)
        beams[i+1] = beam
    end
    return function(z)
        if z < beams[1].z || z > beams[length(beams)].z
            if isnothing(outside)
                throw(DomainError(string(
                    "system does not cover ",
                    "the requested beam position"
                )))
            end
            return outside
        end
        # find last element before requested beam position z; note
        # that this could be accelerated by a more intelligent
        # algorithm since the beams[i].z should be sorted (and this
        # algorithm already depends on this expected ordering).
        i = 0
        while beams[i+1].z < z
            i += 1
        end
        i = max(i, 1)
        # add the effect of a FreeSpace element to reach the requested
        # beam position z
        beam = transform(FreeSpace(z - beams[i].z), beams[i])
        return spotsize(beam)
    end
end

discretize(e::FreeSpace, N::Int) = fill(FreeSpace(e.L/N), N)
discretize(e::Element, N::Int) = e
discretize(els::Vector{<:Element}, N::Int) = vcat(discretize.(els,N)...)


text:
@testset "$(rpad("Lobatto-Legendre",80))" begin

    @test_throws ErrorException LobattoLegendreQuadrature(1)

    for s in 2:10
        @test LobattoLegendreQuadrature(s) ≈  LobattoLegendreQuadrature(s; fast=true)
        @test LobattoLegendreQuadrature(s) == LobattoLegendreQuadrature(Float64, s)

        @test sum(weights(LobattoLegendreQuadrature(s))) ≈ 1
    end

end

ext: jl


text:
module FDMTests

include("../test_fdm.jl")

nparts = (2,2,2)
prun(test_fdm,sequential,nparts)

nparts = 4
prun(test_fdm,sequential,nparts)

end # module


text:
# random functions that dispatch either to rocRAND or GPUArrays' generic RNG

using Random

const GPUARRAY_RNG = Ref{Union{Nothing,GPUArrays.RNG}}(nothing)

function GPUArrays.default_rng(::Type{<:ROCArray})
    if GPUARRAY_RNG[] == nothing
        agent = AMDGPU.default_device().device.agent
        p = Ref{UInt32}()
        GC.@preserve p begin
            AMDGPU.getinfo(agent, HSA.AGENT_INFO_WORKGROUP_MAX_SIZE, p) |> AMDGPU.check
            N = Int(p[])
        end
        state = ROCArray{NTuple{4, UInt32}}(undef, N)
        GPUARRAY_RNG[] = GPUArrays.RNG(state)
        Random.seed!(GPUARRAY_RNG[])
    end
    return GPUARRAY_RNG[]::GPUArrays.RNG
end

gpuarrays_rng() = GPUArrays.default_rng(ROCArray)
const rocrand_rng = librocrand !== nothing ? rocRAND.default_rng : gpuarrays_rng

# the interface is split in two levels:
# - functions that extend the Random standard library, and take an RNG as first argument,
#   will only ever dispatch to rocRAND and as a result are limited in the types they support.
# - functions that take an array will dispatch to either rocRAND or GPUArrays
# - non-exported functions are provided for constructing GPU arrays from only an eltype

# these types are supported by rocrand and we use rocrand whenever possible
const rocRANDUniformType = Union{Type{UInt32}, Type{Cuchar}, Type{UInt16}, Type{Float16}, Type{Float32}, Type{Float64}}
const rocRANDUniformArray = ROCArray{<:Union{UInt32, Cuchar, UInt16, Float16, Float32,Float64}}
const rocRANDNormalType = Union{Type{Float16}, Type{Float32}, Type{Float64}}
const rocRANDNormalArray = ROCArray{<:Union{Float16, Float32,Float64}}

function seed!(seed=Base.rand(UInt64))
    Random.seed!(rocrand_rng(), seed)
    Random.seed!(gpuarrays_rng(), seed)
end

# rocRAND in-place
Random.rand!(A::rocRANDUniformArray) = Random.rand!(rocrand_rng(), A)
Random.randn!(A::rocRANDNormalArray; kwargs...) = Random.randn!(rocrand_rng(), A; kwargs...)

# FIXME: GPUArrays rng does not support Float16
Random.rand!(rng::GPUArrays.RNG, x::ROCArray{Float16}) =
    error("GPUArrays.jl does not support Float16, use rocRAND instead.")
Random.randn!(rng::GPUArrays.RNG, x::ROCArray{Float16}) =
    error("GPUArrays.jl does not support Float16, use rocRAND instead.")

# rocRAND out-of-place
rand(T::rocRANDUniformType, dims::Dims) = Random.rand(rocrand_rng(), T, dims)
randn(T::rocRANDNormalType, dims::Dims; kwargs...) = Random.randn(rocrand_rng(), T, dims; kwargs...)

# support all dimension specifications
rand(T::rocRANDUniformType, dim1::Integer, dims::Integer...) =
    Random.rand(rocrand_rng(), T, Dims((dim1, dims...)))
randn(T::rocRANDNormalType, dim1::Integer, dims::Integer...; kwargs...) =
    Random.randn(rocrand_rng(), T, Dims((dim1, dims...)); kwargs...)

# GPUArrays in-place
Random.rand!(A::ROCArray) = Random.rand!(gpuarrays_rng(), A)
Random.randn!(A::ROCArray; kwargs...) =
    error("AMDGPU.jl does not support generating normally-distributed random numbers of type $(eltype(A))")

# GPUArrays out-of-place
rand(T::Type, dims::Dims) = Random.rand!(ROCArray{T}(undef, dims...))
randn(T::Type, dims::Dims; kwargs...) = Random.randn!(ROCArray{T}(undef, dims...); kwargs...)

# support all dimension specifications
rand(T::Type, dim1::Integer, dims::Integer...) =
    Random.rand!(ROCArray{T}(undef, dim1, dims...))
randn(T::Type, dim1::Integer, dims::Integer...; kwargs...) =
    Random.randn!(ROCArray{T}(undef, dim1, dims...); kwargs...)

# untyped out-of-place
rand(dim1::Integer, dims::Integer...) =
    Random.rand(rocrand_rng(), Dims((dim1, dims...)))
randn(dim1::Integer, dims::Integer...; kwargs...) =
    Random.randn(rocrand_rng(), Dims((dim1, dims...)); kwargs...)

# rand_logn, rand_poisson
const rand_logn = librocrand !== nothing ? rocRAND.rand_logn : (x...;kwargs...) -> error("Not supported without rocRAND.")
const rand_poisson = librocrand !== nothing ? rocRAND.rand_poisson : (x...;kwargs...) -> error("Not supported without rocRAND.")


text:
# sbf_header.jl
#
#   Read the header from an SBF binary file and return the header info and an
#		open file reference.
#
# Jeremy Rüffer
# Thünen Institut
# Institut für Agrarklimaschutz
# Junior Research Group NITROSPHERE
# Julia 1.4.0
# 25.01.2018
# Last Edit: 02.04.2020

"""# sbf_header

Load SBF (Stoffflüsse Binary Format) file header information

---

### Examples

`description,cols,units,types,conversions,lsize = sbf_header(\"K:\\\\Data\\\\Site_20171222_125500.sbf\")` Load the SBF file header information.\n\n
`fid, description,cols,units,types,conversions,lsize = sbf_header(fid)` Load the SBF file header information.\n\n

---

`sbf_load(F,verbose=false,veryverbose=false)`\n
* **F**::String = File name, including path
* **verbose**::Bool (optional) = Display information about the loading/processing, FALSE is default
* **veryverbose**::Bool (optional) = Display column names, units, and conversions for each file loaded, FALSE is default

`sbf_header(source::String)`\n
* **source**::String = Source SBF file, including path

`sbf_header(fid::IOStream)`\n
* **fid**::IOStream = File ID, already opened SBF file ID
* **description**::String = File description
* **cols**::Array{String} = Column names
* **units**::Array{String} = Units of each column
* **types**::Array{DataType} = Types of each column
* **conversions**::Array{String} = Conversions for each column
* **lsize**::Int64 = Line size
"""
function sbf_header(src::String)
	if !isfile(src)
		error(src * " doesn''t exist")
	end

	fid = open(src,"r")
	magicbytes = read!(fid,Array{Int8}(undef,2))
	if magicbytes != Int8[63,42]
		close(fid)
		error(src * " is not a CLD/QCL/Gill binary file")
	end

	fid,description,cols,units,types,conversions,lsize = sbf_header(fid) # Read the header info
	close(fid)

	return description,cols,units,types,conversions,lsize
end

function sbf_header(fid::IOStream)
	function parse_header(h::String)
		f = findall((in)('\x1e'),h)
		parsed_header = Array{String}(undef,length(f) + 1)
		string_pos = 1
		for i=1:length(f)
			parsed_header[i] = h[string_pos:f[i]-1]
			string_pos = f[i] + 1
		end
		parsed_header[end] = h[string_pos:end]

		return parsed_header
	end

	#################################
	## Load and Parse Header Info  ##
	#################################
	read!(fid,Array{Int8}(undef,2)) # \x01, start of the header
	description = String(readuntil(fid,UInt8(03))[1:end]) # Starts with \x02, ends with \x03, \x1e delimited
	column_names_str = String(readuntil(fid,UInt8(03))[2:end]) # Starts with \x02, ends with \x03, \x1e delimited
	units_str = String(readuntil(fid,UInt8(03))[2:end]) # Starts with \x02, ends with \x03, \x1e delimited
	types_str = String(readuntil(fid,UInt8(03))[2:end]) # Starts with \x02, ends with \x03, \x1e delimited
	conversions_str = String(readuntil(fid,UInt8(03))[2:end]) # Starts with \x02, ends with \x03, \x1e delimited
	read!(fid,Array{Int8}(undef,1)) # \x04, end of the header

	cols = parse_header(column_names_str)
	units = parse_header(units_str)
	types_str_list = parse_header(types_str)
	conversions = parse_header(conversions_str)
	types = Array{DataType}(undef,length(cols))
	lsize = 9 # Bytes per line, start with 9 to account for the LF character at the end and the Float64 for the timestamp at the beginning
	types[1] = DateTime
	for i=2:length(types)
		types[i] = eval(Meta.parse(types_str_list[i]))
		if occursin("64",types_str_list[i])
			lsize += 8
		elseif occursin("32",types_str_list[i])
			lsize += 4
		elseif occursin("16",types_str_list[i])
			lsize += 2
		elseif occursin("8",types_str_list[i])
			lsize += 1
		end
	end

	# Clean up column names
	for i=1:length(cols)
		cols[i] = replace(cols[i],' ' => '_')
	end

	return fid, description,cols,units,types,conversions,lsize
end


text:
using Mimi

@defcomp PAGE09Discontinuity begin

    region = Index()
    y_year = Parameter(index=[time], unit="year")
    y_year_0 = Parameter(unit="year")

    rand_discontinuity = Parameter(unit="unitless", default=.5)

    irefeqdis_eqdiscimpact = Variable(index=[region], unit="%")
    wincf_weightsfactor = Parameter(index=[region], unit="unitless", default=[1, 0.8, 0.8, 0.4, 0.8, 0.8, 0.6, 0.6])
    wdis_gdplostdisc = Parameter(unit="%", default=15.)

    igdpeqdis_eqdiscimpact = Variable(index=[time,region], unit="%")
    rgdp_per_cap_NonMarketRemainGDP = Parameter(index=[time,region], unit="\$/person")
    GDP_per_cap_focus_0_FocusRegionEU = Parameter(unit="\$/person", default=27934.244777382406)
    ipow_incomeexponent = Parameter(unit="unitless", default=-0.13333333333333333)

    igdp_realizeddiscimpact = Variable(index=[time,region], unit="%")
    occurdis_occurrencedummy = Variable(index=[time], unit="unitless")
    expfdis_discdecay = Variable(index=[time], unit="unitless")

    distau_discontinuityexponent = Parameter(unit="unitless", default=90.)

    idis_lossfromdisc = Variable(index=[time], unit="degreeC")
    tdis_tolerabilitydisc = Parameter(unit="degreeC", default=3.)
    rt_g_globaltemperature = Parameter(index=[time], unit="degreeC")
    pdis_probability = Parameter(unit="%/degreeC", default=20.)

    isatg_saturationmodification = Parameter(unit="unitless")
    isat_satdiscimpact = Variable(index=[time,region], unit="%")

    isat_per_cap_DiscImpactperCapinclSaturation = Variable(index=[time,region], unit="\$/person")
    rcons_per_cap_DiscRemainConsumption = Variable(index=[time, region], unit="\$/person")
    rcons_per_cap_NonMarketRemainConsumption = Parameter(index=[time, region], unit="\$/person")

    function run_timestep(p, v, d, t)
        wdis_gdplostdisc_2009 = quantile(TriangularDist(5, 25, 15), cdf(TriangularDist(1, 5, 3), p.wdis_gdplostdisc))
        distau_discontinuityexponent_2009 = quantile(TriangularDist(20, 200, 50), cdf(TriangularDist(10, 30, 20), p.distau_discontinuityexponent))
        tdis_tolerabilitydisc_2009 = quantile(TriangularDist(2, 4, 3), cdf(TriangularDist(1, 2, 1.5), p.tdis_tolerabilitydisc))

        v.idis_lossfromdisc[t] = max(0, p.rt_g_globaltemperature[t] - tdis_tolerabilitydisc_2009)

        if is_first(t)
            if v.idis_lossfromdisc[t] * (p.pdis_probability / 100) > p.rand_discontinuity
                v.occurdis_occurrencedummy[t] = 1
            else
                v.occurdis_occurrencedummy[t] = 0
            end
            v.expfdis_discdecay[t] = exp(-(p.y_year[t] - p.y_year_0) / distau_discontinuityexponent_2009)
        else
            if v.idis_lossfromdisc[t] * (p.pdis_probability / 100) > p.rand_discontinuity
                v.occurdis_occurrencedummy[t] = 1
            elseif v.occurdis_occurrencedummy[t - 1] == 1
                v.occurdis_occurrencedummy[t] = 1
            else
                v.occurdis_occurrencedummy[t] = 0
            end
            v.expfdis_discdecay[t] = exp(-(p.y_year[t] - p.y_year[t - 1]) / distau_discontinuityexponent_2009)
        end

        for r in d.region
            v.irefeqdis_eqdiscimpact[r] = p.wincf_weightsfactor[r] * wdis_gdplostdisc_2009

            v.igdpeqdis_eqdiscimpact[t,r] = v.irefeqdis_eqdiscimpact[r] * (p.rgdp_per_cap_NonMarketRemainGDP[t,r] / p.GDP_per_cap_focus_0_FocusRegionEU)^p.ipow_incomeexponent

            if is_first(t)
                v.igdp_realizeddiscimpact[t,r] = v.occurdis_occurrencedummy[t] * (1 - v.expfdis_discdecay[t]) * v.igdpeqdis_eqdiscimpact[t,r]
            else
                v.igdp_realizeddiscimpact[t,r] = v.igdp_realizeddiscimpact[t - 1,r] + v.occurdis_occurrencedummy[t] * (1 - v.expfdis_discdecay[t]) * (v.igdpeqdis_eqdiscimpact[t,r] - v.igdp_realizeddiscimpact[t - 1,r])
            end

            if v.igdp_realizeddiscimpact[t,r] < p.isatg_saturationmodification
                v.isat_satdiscimpact[t,r] = v.igdp_realizeddiscimpact[t,r]
            else
                v.isat_satdiscimpact[t,r] = p.isatg_saturationmodification + (100 - p.isatg_saturationmodification) * ((v.igdp_realizeddiscimpact[t,r] - p.isatg_saturationmodification) / ((100 - p.isatg_saturationmodification) + (v.igdp_realizeddiscimpact[t,r] - p.isatg_saturationmodification)))
            end
            v.isat_per_cap_DiscImpactperCapinclSaturation[t,r] = (v.isat_satdiscimpact[t,r] / 100) * p.rgdp_per_cap_NonMarketRemainGDP[t,r]
            v.rcons_per_cap_DiscRemainConsumption[t,r] = p.rcons_per_cap_NonMarketRemainConsumption[t,r] - v.isat_per_cap_DiscImpactperCapinclSaturation[t,r]
        end
    end
end


text:
########################################################
# Functions related to fitting case series to SIRX model
########################################################

import Base.summary
using DataFrames

struct CaseModelFitResult
    κ::Float64
    κ0::Float64
    IXRatio::Float64
    inputs::NamedTuple
    fit::LsqFit.LsqFitResult
end

"""
fitCaseModel(nt::Int64,cases::Vector{Float64},
    N::Float64,R0Free::Float64,TInfected::Float64,
    p0::Vector{Float64})

- nt: Integer with length of series
- cases: Vector of case number series
- N: Population size
- R0Free: User supplied R0 in freely spreading population
- TInfected: Average time of infection
- p0: Floating point vector with initial values of estimates
"""
function fitCaseModel(nt::Int64,cases::Vector{Float64},
    N::Float64,R0Free::Float64,TInfected::Float64,
    p0::Vector{Float64})
    inputs = (nt=nt,C0=cases[1],N=N,R0Free=R0Free,TInfected=TInfected,C=cases)
    model(t,p) = logCaseModelUnknown(nt,exp.(p),N,cases[1],R0Free,TInfected)
    fit = curve_fit(model,(1:nt)*1.0,log.(cases),log.(p0))
    return CaseModelFitResult( exp(fit.param[1]), exp(fit.param[2]),
    exp(fit.param[3]), inputs, fit )
end

function summary(fit::CaseModelFitResult)
    r0eff = R0Eff(getParams(fit.κ,fit.κ0,fit.inputs.R0Free,
            fit.inputs.TInfected,SIRX()))
    df = exp.(DataFrame(LsqFit.confidence_interval(fit.fit)))
    rename!(df,[:lo,:hi])
    param = [ "κ", "κ0", "IXRatio" ]
    est = [ fit.κ, fit.κ0, fit.IXRatio ]
    insertcols!(df,1,parameter=param)
    insertcols!(df,3,est=est)
    return (coef=df,R0Eff=r0eff,mse=LsqFit.mse(fit.fit))
end

function fitted(fit::CaseModelFitResult,nt::Int64=0)
    if (nt==0)
        nt = fit.inputs.nt
    end
    return exp.(logCaseModelUnknown(nt,
                exp.(fit.fit.param), fit.inputs.N, fit.inputs.C0,
                fit.inputs.R0Free,fit.inputs.TInfected))
end

function estimatedStates(fit::CaseModelFitResult,nt::Int64=0)
    if (nt==0)
        nt = fit.inputs.nt
    end
    d = getParams(exp(fit.fit.param[1]),
                  exp(fit.fit.param[2]), fit.inputs.R0Free,
                  fit.inputs.TInfected, SIRX() )
    return estimatedStates( nt, fit.inputs.N, fit.inputs.C0,
           exp(fit.fit.param[3]), d )
end


text:
using DataFrames
using Dates
using LinearAlgebra
using Statistics: std, var, mean, cov
using HighFrequencyCovariance
using Random
using Test

brownian_corr_matrix = Hermitian([1.0 0.75 0.5 0.0;
                                  0.0 1.0 0.5 0.25;
                                  0.0 0.0 1.0 0.25;
                                  0.0 0.0 0.0 1.0])
assets = [:BARC, :HSBC, :VODL, :RYAL]
twister = MersenneTwister(1)
time_period_per_unit = Dates.Hour(1)

ts0, true_covar, micro_noise, update_rates = generate_random_path(4, 4; brownian_corr_matrix = brownian_corr_matrix, assets = assets, vols = [0.02,0.03,0.04,0.05], twister = deepcopy(twister), syncronous = true, time_period_per_unit = time_period_per_unit)
ts1, true_covar, micro_noise, update_rates = generate_random_path(4, 2000; brownian_corr_matrix = brownian_corr_matrix, assets = assets, vols = [0.02,0.03,0.04,0.05], twister = deepcopy(twister), syncronous = true, time_period_per_unit = time_period_per_unit)
ts2, true_covar, micro_noise, update_rates = generate_random_path(4, 50000; brownian_corr_matrix = brownian_corr_matrix, assets = assets, vols = [0.02,0.03,0.04,0.05], twister = deepcopy(twister), syncronous = true, time_period_per_unit = time_period_per_unit)

iscloser(a,b) = (a.Correlation_error + a.Volatility_error < b.Correlation_error + b.Volatility_error)

# Simple Volatility
simple_vol = simple_volatility(ts1)
all(values(simple_vol) .< 0.1)
simple_vol = simple_volatility(ts1; use_all_obs = true)
all(values(simple_vol) .< 1)
simple_vol = simple_volatility(ts1; fixed_spacing = 20.0)
all(values(simple_vol) .< 0.1)
simple_vol = simple_volatility(ts1; fixed_spacing = Dict(assets .=> repeat([20],4)))
all(values(simple_vol) .< 0.1)

# Preav Convergence
@test_logs (:warn,"Cannot estimate the correlation matrix with 4 ticks. There are insufficient ticks for [:BARC, :HSBC, :VODL, :RYAL]") preaveraged_covariance(ts0, assets) # This will not work due to insufficient data.
preav_estimate1 = preaveraged_covariance(ts1, assets)
preav_estimate2 = preaveraged_covariance(ts2, assets; regularisation = missing)
iscloser(calculate_mean_abs_distance(preav_estimate2, true_covar), calculate_mean_abs_distance(preav_estimate1, true_covar))
valid_correlation_matrix(preav_estimate1)
valid_correlation_matrix(preav_estimate2)

# simple Convergence
simple_estimate1 = simple_covariance(ts1, assets; refresh_times = true) # Testing refresh sampling
simple_estimate1 = simple_covariance(ts1, assets; fixed_spacing = 10) # Testing fixed spacing
simple_estimate1 = simple_covariance(ts1, assets)
simple_estimate2 = simple_covariance(ts2, assets)
valid_correlation_matrix(simple_estimate1)
valid_correlation_matrix(simple_estimate2)

# bnhls Convergence
@test_logs (:warn,"Cannot estimate the correlation matrix with the bnhls method with only 4 ticks.") bnhls_covariance(ts0, assets) # This will not work due to insufficient data.
bnhls_estimate1 = bnhls_covariance(ts1, assets; kernel = quadratic_spectral)
bnhls_estimate1 = bnhls_covariance(ts1, assets; kernel = fejer)
bnhls_estimate1 = bnhls_covariance(ts1, assets; kernel = tukey_hanning)
bnhls_estimate1 = bnhls_covariance(ts1, assets; kernel = bnhls_2008)
bnhls_estimate1 = bnhls_covariance(ts1, assets)
bnhls_estimate2 = bnhls_covariance(ts2, assets)
iscloser(calculate_mean_abs_distance(bnhls_estimate2, true_covar), calculate_mean_abs_distance(bnhls_estimate1, true_covar))
valid_correlation_matrix(bnhls_estimate1) == false # This comes out negative. Not a bug just bad luck with this algo.
valid_correlation_matrix(bnhls_estimate2)

# Preav Convergence
spectral_estimate1 = spectral_covariance(ts1, assets; num_blocks = 1) # not many observations so need to reduce the number of blocks here
spectral_estimate2 = spectral_covariance(ts2, assets)
iscloser(calculate_mean_abs_distance(spectral_estimate2, true_covar), calculate_mean_abs_distance(spectral_estimate1, true_covar))
valid_correlation_matrix(spectral_estimate1)
valid_correlation_matrix(spectral_estimate2)

# two scales Convergence
two_scales_estimate1 = two_scales_covariance(ts1, assets)
two_scales_estimate2 = two_scales_covariance(ts2, assets)
valid_correlation_matrix(two_scales_estimate1)
valid_correlation_matrix(two_scales_estimate2)

#############################
# Conversion to a covariance function
dura = Dates.Minute(137)
covarr = covariance(two_scales_estimate2, dura )
corr, vols =  cov2cor_and_vol(covarr, Nanosecond(two_scales_estimate2.time_period_per_unit).value, Nanosecond(dura).value     )
corr, vols =  cov2cor_and_vol(covarr,   Nanosecond(dura).value / Nanosecond(two_scales_estimate2.time_period_per_unit).value     )
all(corr .< two_scales_estimate2.correlation) .< 10*eps()
all(vols .< two_scales_estimate2.volatility) .< 10*eps()

#############################
 # Serialisation and deserialisation

true_df = to_dataframe(true_covar, Dict([:estimation] .=> ["True Covariance Matrix"]))
reconstituted_df = dataframe_to_covariancematrix(true_df)
calculate_mean_abs_distance(true_covar, reconstituted_df).Correlation_error .< 10*eps()
calculate_mean_abs_distance(true_covar, reconstituted_df).Volatility_error .< 10*eps()
# Reshuffling rows to make sure it still works.
true_df2 = true_df[randperm(MersenneTwister(1234), nrow(true_df)),:]
reconstituted_df = dataframe_to_covariancematrix(true_df2)
calculate_mean_abs_distance(true_covar, reconstituted_df).Correlation_error .< 10*eps()
calculate_mean_abs_distance(true_covar, reconstituted_df).Volatility_error .< 10*eps()

# Other regularisation algos:
two_scales_estimate_iden = two_scales_covariance(ts2, assets; regularisation = :identity_regularisation)
two_scales_estimate_nearest_corr = two_scales_covariance(ts2, assets; regularisation = :nearest_correlation_matrix)
two_scales_estimate_nearest_psd = two_scales_covariance(ts2, assets; regularisation = :nearest_psd_matrix)
two_scales_estimate_eigen = two_scales_covariance(ts2, assets; regularisation = :eigenvalue_clean)
# Testing that these are different (due to different regularisation)
dist = calculate_mean_abs_distance(two_scales_estimate_iden, two_scales_estimate_nearest_corr)
dist.Correlation_error .> 1000*eps()
dist = calculate_mean_abs_distance(two_scales_estimate_iden, two_scales_estimate_eigen)
dist.Correlation_error .> 1000*eps()

# Running regularistation on a CovarianceMatrix's correlation matrix.
psd_mat = two_scales_covariance(ts2, assets; regularisation = :nearest_correlation_matrix)
valid_correlation_matrix(psd_mat)
reg1 = identity_regularisation(psd_mat, 0.5; apply_to_covariance = false) # Inputting a weight explicitly.
calculate_mean_abs_distance(psd_mat, reg1).Correlation_error .> 10*eps()
calculate_mean_abs_distance(psd_mat, reg1).Volatility_error .< 10*eps()

reg1_corr = identity_regularisation(psd_mat, ts2; apply_to_covariance = false)
calculate_mean_abs_distance(psd_mat, reg1_corr).Correlation_error .> 10*eps()
calculate_mean_abs_distance(psd_mat, reg1_corr).Volatility_error .< 10*eps()

reg2 = nearest_correlation_matrix(psd_mat, ts2)
calculate_mean_abs_distance(psd_mat, reg2).Correlation_error .< 10*eps() # No change as we are already psd
calculate_mean_abs_distance(psd_mat, reg2).Volatility_error .< 10*eps()

reg3 = nearest_psd_matrix(psd_mat, ts2)
calculate_mean_abs_distance(psd_mat, reg3).Correlation_error .< 10*eps() # No change as we are already psd
calculate_mean_abs_distance(psd_mat, reg3).Volatility_error .< 1000*eps()

reg3_corr = nearest_psd_matrix(psd_mat, ts2; apply_to_covariance = false)
calculate_mean_abs_distance(psd_mat, reg3_corr).Correlation_error .< 10*eps() # No change as we are already psd
calculate_mean_abs_distance(psd_mat, reg3_corr).Volatility_error .< 10*eps()

reg4 = eigenvalue_clean(psd_mat, ts2)
calculate_mean_abs_distance(psd_mat, reg4).Correlation_error .> 2*eps()
calculate_mean_abs_distance(psd_mat, reg4).Volatility_error .> 10*eps()

reg4_cov = eigenvalue_clean(psd_mat, ts2; apply_to_covariance = false)
calculate_mean_abs_distance(psd_mat, reg4_cov).Correlation_error .> 2*eps()
calculate_mean_abs_distance(psd_mat, reg4_cov).Volatility_error .< 10*eps()

# Testing the covariance Function
all((covariance(true_covar, Hour(1)) * 24 .-   covariance(true_covar, Day(1))   ) .< 10*eps())


text:
using Documenter, MORFEInvariantManifold

makedocs(sitename="MORFEInvariantManifold")

makedocs(
    sitename = "Model Order Reduction in Julia",
    format = Documenter.HTML(),
    authors = "Andrea Opreni, Alessandra Vizzaccaro",
    modules = [MORFEInvariantManifold],
    pages = Any[
        "Home" => "index.md",
        "Overview" => "overview.md",
        "Tutorials" => "tutorials.md",
        "Supported Mesh Formats" => "supported_formats.md",
        "Integrators" => "integrators.md",
        "Library" => "library_all.md",
        "Authors" => "authors.md",
        "References" => "literature.md"
    ]
)


text:
#!/usr/bin/env julia
include("secirc9717eurshort.jl")
sc9717plall(sc9717sex(2))
savefig("../../images/secirc9717eurshort_kv.svg")
close()
sc9717plall(sc9717sex(1))
savefig("../../images/secirc9717eurshort_m.svg")

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
using Pkg
Pkg.add(["Pluto", "PlutoUI", "Plots"])
Pkg.precompile()


text:
# ----------------------------------------------------------------------------------------------------------
# l-∞ norm
# ----------------------------------------------------------------------------------------------------------
"""
    Linf

An empty structure representing l-infinity norm.
"""
struct Linf <: Norm end

"""
    name(d::Linf)

Returns the full name of the l-infinity norm.
"""
name(d::Linf) = "l-infinity norm"

"""
    normtype(d::Linf)

Returns the type of l-infinity norm.
"""
normtype(d::Linf) = Inf

"""
    optimal(s::Our, d::Linf, m::DRO; kwargs...)

Returns the optimal solution of the DRO model 'm' with l-infinity norm.
"""
function optimal(s::Our, d::Linf, m::DRO)
    c      = .- m.c
    perm   = sortperm(c)
    c      = c[perm]
    p      = m.q[perm]
    i      = 1
    j      = length(c)
    δ_down = 0

    while i < j
        δ1    = min(1 - p[i], m.ε)
        δ2    = 0
        p[i] += δ1

        while δ2 < δ1
            if min(p[j], m.ε - δ_down) >= δ1 - δ2
                p[j]   -= δ1 - δ2
                δ_down += δ1 - δ2
                break
            else
                δ2     += min(p[j], m.ε - δ_down)
                p[j]   -= min(p[j], m.ε - δ_down)
                δ_down  = 0
                j      -= 1

                if i == j
                    p[i] += - δ1 + δ2
                    break
                end
            end
        end
        i += 1
    end
    return p[invperm(perm)]
end

# ----------------------------------------------------------------------------------------------------------
# l-1 norm
# ----------------------------------------------------------------------------------------------------------
"""
    Lone

An empty structure representing l-1 norm.
"""
struct Lone <: Norm end

"""
    name(d::Lone)

Returns the full name of the l-1 norm.
"""
name(d::Lone) = "l-1 norm"

"""
    normtype(d::Lone)

Returns the type of l-1 norm.
"""
normtype(d::Lone) = 1

"""
    optimal(s::Our, d::Lone, m::DRO)

Returns the optimal solution of the DRO model 'm' with l-1 norm.
"""
function optimal(s::Our, d::Lone, m::DRO; kwargs...)
    c    = .- m.c
    perm = sortperm(c)
    c    = c[perm]
    p    = m.q[perm]
    i    = 1
    j    = length(c)
    δ_up = 0

    while i < j && δ_up < m.ε/2
        δ1   = min(1 - p[i], m.ε/2-δ_up)
        δ_up += δ1
        δ2   = 0
        p[i] += δ1

        while δ2 < δ1
            if p[j] >= δ1 - δ2
                p[j] += - δ1 + δ2
                break
            else
                δ2 = δ2 + p[j]
                p[j] = 0
                j -= 1
                if i == j
                    p[i] += - δ1 + δ2
                    break
                end
            end
        end
        i += 1
    end
    return p[invperm(perm)]
end

# ----------------------------------------------------------------------------------------------------------
# l-2 norm
# ----------------------------------------------------------------------------------------------------------
"""
    Ltwo

An empty structure representing l-2 norm.
"""
struct Ltwo <: Norm end

"""
    name(d::Ltwo)

Returns the full name of the l-2 norm.
"""
name(d::Ltwo) = "l-2 norm"

"""
    normtype(d::Ltwo)

Returns the type of l-2 norm.
"""
normtype(d::Ltwo) = 2

"""
    initial(d::Ltwo, m::DRO)

Returns the initial point for finding the root of the function `h` using the newton method.
"""
function initial(d::Ltwo, m::DRO)
    f(μ)  = h(d, m, μ)
    μ     = 10
    f_val = f(μ)
    add_eval()

    while f_val >= 0
        f == 0 && return μ

        μ    *= 10
        f_val = f(μ)
        add_eval()
    end
    return μ
end

function g(d::Ltwo, m::DRO, μ::Real)
    n     = length(m.q)
    s     = sort(μ.*m.q .+ m.c)
    λ     = 0
    g     = μ
    g_old = g

    for i in n-1:-1:1
        g_old = g
        g    -= (n - i)*(s[i + 1] - s[i])
        if g <= 0
            λ = s[i] - (s[i + 1] - s[i])/(g_old - g)*g
            break
        end
    end
    if g > 0
        λ = Statistics.mean(m.c)
    end
    return λ
end

function h(d::Ltwo, m::DRO, μ::Real; λ::Real = g(d, m, μ))
    return sum((min.(λ .- m.c, μ .* m.q)).^2) - m.ε^2 * μ^2
end

function ∇h(d::Ltwo, m::DRO, μ::Real; λ::Real = g(d, m, μ))
    q_i = @view m.q[λ .- m.c .>  μ .* m.q]
    c_i = @view m.c[λ .- m.c .<= μ .* m.q]

    return 2*μ*(sum(abs2, q_i) - m.ε^2) - 2*sum(q_i)*sum(λ .- c_i)/length(c_i)
end

"""
    optimal(s::Our, d::Ltwo, m::DRO; kwargs...)

Returns the optimal solution of the DRO model 'm' with l-2 norm.
"""
function optimal(s::Our, d::Ltwo, m::DRO; kwargs...)
    Ilen = length(m.Imax)
    Isum = sum(m.q[m.Imax])
    p    = zero(m.q)
    p[m.Imax] .= m.q[m.Imax] .+ 1/Ilen .- Isum/Ilen

    if sum(abs2, p - m.q) <= m.ε^2
        return p
    else
        μ = newton(d, m; kwargs...)
        λ = g(d, m, μ)
        return max.(m.q .- (λ .- m.c)./μ, 0)
    end
end

"""
    optimal(s::Philpott, d::Philpott, m::DRO; atol::Real = 1e-10)

Returns the optimal solution given by Philpott's algorithm of the DRO model 'm' with l-2 norm.
"""
function optimal(s::Philpott, d::Ltwo, m::DRO; atol::Real = 1e-10, kwargs...)
    n      = length(m.q)
    c_mean = Statistics.mean(m.c)
    c_std  = Statistics.stdm(m.c, c_mean; corrected = false)
    p      = m.q .+ m.ε .* (m.c .- c_mean) ./ (sqrt(n) * c_std)

    isapprox(c_std, 0, atol = atol)         && return m.q
    m.ε <= sqrt(n / (n - 1)) * minimum(m.q) && return p

    I  = collect(1:n)
    j  = 0
    a1 = 0
    a2 = 0
    a3 = sqrt(n)*m.ε

    @views for I_len in n-1:-1:1

        minimum(p[I]) >= 0 && return p

        ## find index j
        J   = I[p[I] .< 0]
        ε_p = @. (a1^2 + (c_std * (I_len * m.q[J] + a1) / (m.c[J] - c_mean))^2) / I_len + a2
        j   = J[findmin(ε_p)[2]]

        deleteat!(I, searchsortedfirst(I, j))
        p[j] = 0

        ## update p
        a1 += m.q[j]
        a2 += m.q[j]^2
        a3 = sqrt(I_len*(m.ε^2 - a2) - a1^2)

        c_mean = ((I_len + 1) * c_mean - m.c[j]) / I_len
        c_std  = Statistics.stdm(m.c[I], c_mean; corrected = false)

        @. p[I] = m.q[I] + (a1 + a3 * (m.c[I] - c_mean) / c_std) / I_len

        if isapprox(c_std, 0, atol = atol)
            @. p[I] = m.q[I] + a1 / I_len
            return p
        end
    end
    p[I] = 1
    return p
end

max_stars_repo_name: VaclavMacha/Projections

max_issues_repo_name: VaclavMacha/Projections

max_forks_repo_name: VaclavMacha/Projections


text:
function _slowing_down(E::T, p::T, E0::T, p0::T;
                      th_amu = H2_amu, b_amu = th_amu,
                      Te = 1, Zeff = 1, ne = 1, ni = ne, Zi = 1) where T<:Number
    m_e = mass_u*e_amu
    m_b = mass_u*b_amu
    m_th = mass_u*th_amu

    tau_s = Te^1.5 # there are other factors as well
    E >= E0  && return 0.0
    v = sqrt(2*(e0*E*1e3)/(m_b))
    vb = sqrt(2*(e0*Eb*1e3)/(m_b))
    vc3 = 0.75*sqrt(pi)*((2*(Te*e0*1e3)/m_e)^(1.5))*((ni/ne) * Zi^2 * (m_e/(m_th)))

    beta = (m_th*Zeff)/(2*m_b)
    alpha = (beta/3)*(1 - p0^2)*log((1 + vc3/(v^3))/(1+vc3/(vb^3)))
    g = tau_s*inv(v^3 + vc3)*inv(sqrt(4*pi*alpha))*exp(-((p - p0)^2)/(4*alpha))
end

function slowing_down(energy::AbstractVector, pitch::AbstractVector, E0, p0; kwargs...)

    nenergy = length(energy)
    npitch = length(pitch)

    f_slow = zeros(nenergy,npitch)
    for i=1:nenergy, j=1:npitch
        f_slow[i,j] = _slowing_down(energy[i],pitch[j], E0, p0; kwargs...)
    end

    N, err = hcubature(x->_slowing_down(x[1],x[2], E0, p0; kwargs...), (0.0,-1.0), (1.1*E0,1.0))
    return f_slow/N
end

function _bimaxwellian(energy::S, pitch::S, T_perp::S, T_para::S;
                      vd_para=0.0, sigma=-1, th_amu = H2_amu) where S<:Number

    m_th = mass_u*th_amu

    T_perp = T_perp*1e3*e0
    T_para = T_para*1e3*e0
    Ed_para = (0.5*m_th*vd_para^2)

    E = energy*e0*1e3
    xi = sigma*pitch

    A = sqrt(E/(pi*T_para*T_perp^2))
    T1 = -(E*xi^2 + Ed_para - vd_para*xi*sqrt(2*m_th*E))/T_para
    T2 = -(1-xi^2)*E/T_perp
    return A*exp(T1 + T2)
end

function bimaxwellian(energy::AbstractVector, pitch::AbstractVector, T_perp, T_para; kwargs...)
    nenergy = length(energy)
    npitch = length(pitch)
    f_ep = zeros(nenergy, npitch)
    for i=1:nenergy, j=1:npitch
        f_ep[i,j] = _bimaxwellian(energy[i], pitch[j], T_perp, T_para; kwargs...)
    end
    T = (T_para + 2*T_perp)/3
    N, err = hcubature(x->_bimaxwellian(x[1], x[2], T_perp, T_para; kwargs...), (0.0,-1.0), (10*T,1.0))
    return f_ep/N
end

function maxwellian(energy, pitch, T; kwargs...)
    bimaxwellian(energy, pitch, T, T; kwargs...)
end


text:
using timbpMisc, Test, Dates, Distributions

@testset "timbpMisc.jl" begin
  @testset "Financial years" begin
    d = Dates.Date("2019")
    @test FinancialYear(2019).value == 2019
    @test FinancialYear(d).value == 2019
    @test FinancialYear("2019").value == 2019
    @test FinancialYear("2019-06").value == 2019
    @test FinancialYear("2019-07").value == 2020
    @test FinancialYear("2019-06-30").value == 2019
    @test FinancialYear("2019-07-01").value == 2020
    @test FinancialYear(2019) - FinancialYear(2017) == Year(2)
    @test FinancialYear(2017) + Year(2) == FinancialYear(2019)
    @test FinancialYear(2017) - Year(2) == FinancialYear(2015)
    @test show(FinancialYear(d)) == print("2018/19")
    @test print(FinancialYear(d)) == 2019
    @test (FinancialYear(2019) < FinancialYear(2020)) == true
  end
  @testset "Control limits" begin
    @test controllimits(0.025, 50, .3) == 0.16613938031193454
    @test controllimits(0.975, 50, .3) == 0.4201358867168803
    @test controllimits.((0.025, 0.975), 50, .3) == (0.16613938031193454, 0.4201358867168803)
  end
end


text:
#############################################################################
# fujishige_wolfe.jl
# Implements the the fujishige_wolfe algorithm, to find the point with
# minimum norm in the base polyhedron.
#############################################################################

export fujishige_wolfe

function fujishige_wolfe(p::AssocPoly, w::AbstractArray, Tol::Float64 = 1e-3)
  # step 1 initialization
  n = length(w)
  @assert length(p.V) == n
  V = sort(p.V)         # the sorted base set
  x = zeros(n)          # the candidate
  y = zeros(n + 1)
  for i = 1:n
    y[i + 1] = evaluate(p.F, V[1:i])[1]
    x[i] = y[i + 1] - y[i]
  end
  S = zeros(n, 1)       # the set of selected extreme points that constitute the affine set
  S[:, 1] = x
  # major cycle
  major_cycle(p, w, x, S, Tol)
end

function fujishige_wolfe(p::PosPoly, w::AbstractArray, Tol::Float64 = 1e-3)
  w₁ = copy(w)
  if any(w₁ .< 0)
    negind = find(w₁ .< 0)
    w₁[negind] = zeros(length(negind))
  end
  p₁ = SubmodPoly(p.F)
  fujishige_wolfe(p₁, w₁, Tol)
end

function fujishige_wolfe(p::SymPoly, w::AbstractArray, Tol::Float64 = 1e-3)
  w₁ = sign.(w) .* w
  p₁ = SubmodPoly(p.F)
  x = fujishige_wolfe(p₁, w₁, Tol)
  return sign.(w) .* x
end

# minimum-norm-point with warmstart
function fujishige_wolfe(p::AssocPoly, w::AbstractArray, x₀::AbstractArray, Tol::Float64 = 1e-3)
  n = length(w)
  @assert length(p.V) == n
  S = zeros(n, 1)
  S[:, 1] = x₀
  major_cycle(p, w, x₀, S, Tol)
end

function fujishige_wolfe(p::PosPoly, w::AbstractArray, x₀::AbstractArray, Tol::Float64 = 1e-3)
  w₁ = copy(w)
  if any(w₁ .< 0)
    negind = find(w₁ .< 0)
    w₁[negind] = zeros(length(negind))
  end
  p₁ = SubmodPoly(p.F)
  fujishige_wolfe(p₁, w₁, x₀, Tol)
end

function fujishige_wolfe(p::SymPoly, w::AbstractArray, x₀::AbstractArray, Tol::Float64 = 1e-3)
  w₁ = sign.(w) .* w
  p₁ = SubmodPoly(p.F)
  x = fujishige_wolfe(p₁, w₁, x₀, Tol)
  return sign.(w) .* x
end

function major_cycle(p::AssocPoly, w::AbstractArray, x::AbstractArray, S::AbstractMatrix, Tol::Float64 = 1e-3, nTol::Float64 = 1e-3)
  n = length(w)
  S1 = minimum_check(p, w, x, S, Tol)
	if S1 == S
    return x
  else
    S = S1
    colnum = size(S, 2)
    # step 3 find the point closes to w in the affine hull of points in S
    if colnum > size(S, 1)              # the affine hull is in face the whole space
      y0 = copy(w)
      push!(y0, 1)
      S1 = [S; ones(1, size(S, 2))]
      μ = S1\y0
      if all(μ .>= -nTol)                # μ is in the affine hull of points in set
        return w
      else                              # delete redundant points and compute the new x
        redind0 = find(μ .< -nTol)
        remind0 = deleteat!(collect(1:colnum), redind0)
        S = S[:, remind0]
        colnum = length(remind0)
        y = aff_proj(w, S)
        major_cycle(p, w, y, S, Tol)
      end
    else
      y = aff_proj(w, S)                 # the point found)
      y0 = copy(y)
      push!(y0, 1)
      S₁ = [S; ones(1, size(S, 2))]
      μ = S₁\y0                         # the coefficients of y in terms of points in S
      while any(μ .< -nTol)
        # step 4 find the closest point to y on the segment between x and y while staying in the convex hull of points in S, and set it to be x
        x₀ = copy(x)
        push!(x₀, 1)
        λ = S₁\x₀                           # the coefficients of x in terms of points in S
        ratios = λ ./ μ
        negratios = ratios[find(μ .< -nTol)]
        α = maximum(negratios)
        redind = find((ratios .<= α + nTol) .& (ratios .<= α + nTol))    # the index of the redundant points
        remind = deleteat!(collect(1:colnum), redind)
        S = S[:, remind]                           # the new set of points
        x = y                                      # the new candidate
        y = aff_proj(w, S)                          # the new optimal solution given S
        y0 = copy(y)
        push!(y0, 1)
        S₁ = [S; ones(1, size(S, 2))]
        μ = S₁\y0
      end
      major_cycle(p, w, y, S, Tol)
    end
  end
end

function minimum_check(p::AssocPoly, w::AbstractArray, x::AbstractArray, S::AbstractMatrix, Tol::Float64 = 1e-3)
  a = x - w
  z = fenchel(p, a)
  minprod = dot(z - w, a)
  curprod = dot(a, a)
  if curprod - minprod <= Tol
    return S
  else
    return [S z]
  end
end

function minimum_check(p::SubmodPoly, w::AbstractArray, x::AbstractArray, S::AbstractMatrix, Tol::Float64 = 1e-3)
  a = x - w
  pos = false
  if any(a .> 0)
    pos = true
    posind = find(a .> 0)
  end
  z = greedy(p.F, -a)
  if pos                                                 # x - w has positive elements
    s = z[posind] - w[posind]
    if any(s .> 0)                                       # some elements of z should be adjusted
      posind1 = find(s .> 0)
      z[posind[posind1]] = w[posind[posind1]]
    end
  end
  minprod = dot(z - w, a)
  curprod = dot(a, a)
  if curprod - minprod <= Tol
    return S
  else
    S1 = [S z]
    return [S z]
  end
end

function fujishige_wolfe(f::SubmodFunc, w::AbstractArray, Tol::Float64 = 1e-3)
  p = BasePoly(f)
  return fujishige_wolfe(p, w, Tol)
end


text:
using TestProject

println("In precompile execution file.")


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

function GitStatus(repo::GitRepo; status_opts=StatusOptions())
    stat_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_status_list_new, :libgit2), Cint,
                  (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{StatusOptions}),
                  stat_ptr_ptr, repo.ptr, Ref(status_opts))
    return GitStatus(repo, stat_ptr_ptr[])
end

function Base.length(status::GitStatus)
    return Int(ccall((:git_status_list_entrycount, :libgit2), Csize_t,
                      (Ptr{Ptr{Void}}, ), status.ptr))
end

function Base.getindex(status::GitStatus, i::Integer)
    1 <= i <= length(status) || throw(BoundsError())
    entry_ptr = ccall((:git_status_byindex, :libgit2),
                      Ptr{StatusEntry},
                      (Ptr{Void}, Csize_t),
                      status.ptr, i-1)
    entry_ptr == C_NULL && throw(Error.GitError(Error.ERROR))
    return unsafe_load(entry_ptr)
end

function status(repo::GitRepo, path::String)
    status_ptr = Ref{Cuint}(0)
    ret =  ccall((:git_status_file, :libgit2), Cint,
                  (Ref{Cuint}, Ptr{Void}, Cstring),
                  status_ptr, repo.ptr, path)
    (ret == Cint(Error.ENOTFOUND) || ret == Cint(Error.EAMBIGUOUS)) && return Nullable{Cuint}()
    return Nullable(status_ptr[])
end


text:
#############################
##  These numerical experiments, included for the sake of completion,
##  showcase some numerical solutions to LINEAR Volterra integral equations
#############################

using ApproxFun, MultivariateOrthogonalPolynomials, BandedMatrices, BlockBandedMatrices, SpecialFunctions, Plots
using SparseVolterraExamples

#####################################################
## Solving Equation 18 in https://doi.org/10.1137/19M1267441
## This is a linear Volterra integral equation of FIRST kind.
#####################################################
####
## Set up g(x) and kernel K(x,y)
g1(x)= exp(-x)/4+1/4*exp(x)*(-1+2*x)
K(x,y) = exp(y-x)
####
## For linear problems the solver is wrapped into convenient functions. N is polynomial order of approximation.
N = 20
u = triVolterraEQ1FullKernelSolver(g1,K,N,true)
####
## Plot the computed and analytic solution
plot(u,grid=false,ylabel="u(x)",xlabel="x",label="sparse method")
plot!(x->x*exp(x),0,1,grid=false,ylabel="u(x)",xlabel="x",label="analytic")
####
## Plot error
plot(triVolterraEQ1FullKernelSolver(g1,K,N,true)-Fun(x->x*exp(x),Jacobi(0,1,0..1)),ylabel="error",xlabel="x",grid=false,label=false)

#####################################################
## The following is a simple linear second kind Volterra integral equation, not appearing in either paper.
## It uses the monomial coefficient vector input kernel implementation and limits of integration to 1-x instead and has analytically known solution sin(-2*pi*x).
#####################################################
####
## Set up g(x) and kernel K(x,y)
g2 = Fun(x->(x-1)/(2*pi)*cos(-2*pi*x)+(4*pi^2+1)/(4*pi^2)*sin(2*pi*(1-x)), Jacobi(0,1,0..1))
K2 = [0.,1.,0.,0.] # this is equivalent to K(x,y)=y
####
## Compute and plot solutions
N = 20
u = triVolterraEQ2Solver(g2,K2,N,false) # note the false indicating that we are not flipping the integration bounds, i.e. we are integrating from 0 to 1-x
plot(u,grid=false,ylabel="u(x)",xlabel="x",label="sparse method")
plot!(x->sin(-2*pi*x),0,1,grid=false,ylabel="u(x)",xlabel="x",label="analytic")
####
## Plot errors
plot(triVolterraEQ2Solver(g2,K2,N,false)-Fun(x->sin(-2*pi*x),Jacobi(0,1,0..1)),grid=false,ylabel="error",xlabel="x",label=false)
####
## The above is equivalent to the following much simpler input, which uses the full kernel input as a function
u2 = triVolterraEQ2FullKernelSolver(g2,(x,y)->y,N,false)
####
## Verify the two input methods give the same answer
plot(x->abs(u2(x)-u(x)),0,1,grid=false,ylabel="difference",xlabel="x",label=false)

max_stars_repo_path: examples/Example 0 - Some linear Volterra integral equation examples.jl


text:
# This file was generated, do not modify it. # hide
fp = fitted_params(mach_uni)
@show round.(fp.coefs, sigdigits=3)
@show round(fp.intercept, sigdigits=3)


text:
module SwitchTimeOpt

  # Import Necessary Modules
  using MathOptInterface
  using DiffEqBase
  using OrdinaryDiffEq
  using Ipopt
  using LinearAlgebra
  using SparseArrays
  using Plots

  export stoproblem,
         setwarmstart!,setx0!,
         solve!,
         gettau,
         getdelta,
         gettaucomplete, getdeltacomplete,
         getobjval, getstat, getsoltime,
         getnobjeval, getngradeval, getnhesseval,
         simulatelinearized, simulateinput,
         simulate, plotSolution

  # Include Source Files
  include("types.jl")
  include("interface.jl")
  include("nlp.jl")
  include("simulation.jl")

end