text:
# # 1D and 2D Hawkes processes with exponential kernel

#=
In this example, I simulate first a 1D self-exciting Hawkes process and then a 2D one.
The interaction kernel is an exponentially decaying function, defined as:

```math
g(t) = H(t) \,  \frac{1}{\tau} \, \exp\left(- \frac{t}{\tau} \right)
```

where $H(t)$ is the Heaviside function: zero for $t<0$, one for $t\geq 0$.

Note that Kernels are always normalized so that their integral
between $-\infty$ and $+\infty$ is 1.
=#

# ## Initialization
using LinearAlgebra,Statistics,StatsBase,Distributions
using Plots,NamedColors ; theme(:default) ; gr()
using FFTW

using ProgressMeter
using Random
Random.seed!(0)

using HawkesSimulator; const global H = HawkesSimulator

"""
    onedmat(x::R) where R
Generates a 1-by-1 Matrix{R} that contains `x` as only element
"""
function onedmat(x::R) where R
  return cat(x;dims=2)
end;

## #src
# First I define the kernel, and the self-interaction weight.
#
# `myw` is a scaling factor (the weight of the
# autaptic connection). The baseline rate is given by `myinput`

mytau = 0.5  # kernel time constant
myw = onedmat(0.85) # weight: this needs to be a matrix
myinput = [0.7,] # this needs to be a vector
mykernel = H.KernelExp(mytau);

# This is the plot of the (self) interaction kernel
# (before  the scaling by `myw`)

theplot = let ts = range(-1.0,5;length=150),
  y = map(t->H.interaction_kernel(t,mykernel) , ts )
  plot(ts , y ; linewidth=3,leg=false,xlabel="time (s)",
     ylabel="interaction kernel")
end;
plot(theplot)

# Now I build the network, using the simplified constructor
# for one-population networks.
# The object `PopulationState` initializes a population of
# Hawkes neurons of the specified size (1, in this case)
# and the kernel `mykernel` neurons in the same population
# all have the same kernel.
# (this might be an inconvenience if one wants qualitatively different
# self-connections... but it might be fixed in future versions)

popstate = H.PopulationState(mykernel,1)
ntw = H.RecurrentNetwork(popstate,myw,myinput);

## #src
# ## Simulation
# The length of the simulation is measured by the total number of spikes
# here called `num_spikes`.
# The function `flush_trains!(...)` is used to store older spikes as
# history and let them be ignored by the kernels.
# The time parameters should be regulated based on the kernel shape.

function run_simulation!(network,num_spikes,
    t_flush_trigger=300.0,t_flush=100.0)
  t_now = 0.0
  H.reset!(network) # clear spike trains etc
  @showprogress "Running Hawkes process..." for _ in 1:num_spikes
    t_now = H.dynamics_step_singlepopulation!(t_now,network)
    H.flush_trains!(network,t_flush_trigger;Tflush=t_flush)
  end
  H.flush_trains!(network) # flush everything into history
  return t_now
end

n_spikes = 100_000
Tmax = run_simulation!(ntw,n_spikes)
ratenum = H.numerical_rates(popstate;Tend=Tmax)[1]
@info "Simulation completed, mean rate $(round(ratenum;digits=2)) Hz"

##
# ## Visualize raster plot of the events

# the raster plot shows some correlation between the neural activities
# neuron one (lower row) excites neuron two (upper row)

function rasterplot(tlims = (1100.,1120.) )
  _train = popstate.trains_history[1]
  plt=plot()
  train = filter(t-> tlims[1]< t < tlims[2],_train)
  nspk = length(train)
  scatter!(plt,train,fill(0.1,nspk),markersize=30,
      markercolor=:black,markershape=:vline,leg=false)
  return plot!(plt,ylims=(0.0,0.2),xlabel="time (s)")
end

rasterplot()

# event times are always stored in `pops.trains_history`

##  #src
# ## Plot the instantaneous rate

#=
This is the probability of a spike given the past activity up until that
moment. It is usually denoted by $\lambda^*(t)$.
=#

function get_insta_rate(t,popstate)
  _train = popstate.trains_history[1]
  myinput[1] + H.interaction(t,_train,myw[1],popstate.unittype)
end
function plot_instarate(popstate,tlims=(1100,1120))
  tplot = range(tlims...;length=100)
  _train = popstate.trains_history[1]
  tspk = filter(t-> tlims[1]<=t<=tlims[2],_train) # add the exact spiketimes for cleaner plot
  tplot = sort(vcat(tplot,tspk,tspk .- 1E-4))
  plt=plot(xlabel="time (s)",ylabel="instantaneous rate")
  plot!(plt,t->get_insta_rate(t,popstate),tplot;linewidth=2,color=:black)
  scatter!(t->get_insta_rate(t,popstate),tspk;leg=false)
  avg_rate = H.numerical_rates(popstate;Tend=Tmax)[1]
  ylim2 = ylims()[2]
  plot!(plt,tplot, fill(avg_rate,length(tplot)),
     color=:red,linestyle=:dash,ylims=(0,ylim2))
end
plot_instarate(popstate)

# ## Plot the total event counts

# Total number of events as a function of time. It grows linearly, and the
# steepness is pretty much the rate.

function plot_counts(tlims=(0.,1000.))
  avg_rate = H.numerical_rates(popstate;Tend=Tmax)[1]
  tplot = range(tlims...;length=100)
  _train = popstate.trains_history[1]
  nevents(tnow::Real) = count(t-> t <= tnow,_train)
  plt=plot(xlabel="time (s)",ylabel="number of events",leg=false)
  plot!(plt,tplot,nevents.(tplot),color=:black,linewidth=2)
  plot!(plt,tplot , tplot .* avg_rate,color=:red,linestyle=:dash)
end
plot_counts()

# ## Rate
#=
Now I compare the numerical rate with the analytic solution.

The analytic rate corresponds to the stationary solution
of a linear dynamical system (assumung all stationary rates are above zero).
```math
\frac{\mathrm d \mathbf r}{\mathrm d t} = - \mathbf r + W\;r + \mathbf h \quad;
\qquad  r_\infty = (I-W)^{-1} \; \mathbf h
```
=#

rate_analytic = inv(I-myw)*myinput
rate_analytic = rate_analytic[1] # 1-D , just a scalar

@info "Mean rate -  numerical $(round(ratenum;digits=2)), analytic  $(round(rate_analytic;digits=2))"

##
# ## Covariance density

# TODO : write definition
#
# I compute the covariance density it numerically.
# The time inteval `mydt` should not be too small.

mytrain = popstate.trains_history[1]
mydt = 0.1
myτmax = 25.0
mytaus = H.get_times(mydt,myτmax)
ntaus = length(mytaus)
cov_num = H.covariance_self_numerical(mytrain,mydt,myτmax);

##
# Now I compute the covariance density analytically, at higher resolution,
# and I compare the two.

function four_high_res(dt::Real,Tmax::Real) # higher time resolution, longer time
  k1,k2 = 2 , 0.01
  myτmax = Tmax * k1
  dt *= k2
  mytaus = H.get_times(dt,myτmax)
  nkeep = div(length(mytaus),k1)
  myfreq = H.get_frequencies_centerzero(dt,myτmax)
  gfou = myw[1,1] .* H.interaction_kernel_fourier.(myfreq,Ref(popstate)) |> ifftshift
  ffou = let r=rate_analytic
    covf(g) = r/((1-g)*(1-g'))
    map(covf,gfou)
  end
  retf = real.(ifft(ffou)) ./ dt
  retf[1] *= dt  # first element is rate
  return mytaus[1:nkeep],retf[1:nkeep]
end

(taush,covfou)=four_high_res(mydt,myτmax)

function doplot()
  plt = plot(xlabel="time delay (s)",ylabel="Covariance density")
  plot!(plt,mytaus[2:end], cov_num[2:end] ; linewidth=3, label="simulation" )
  plot!(plt,taush[2:end],covfou[2:end]; label="analytic",linewidth=3,linestyle=:dash)
  return plt
end
doplot()

# 1D system completed !
##
# ## Same results, but in a 2D system

myτ = 1/2.33
mywmat = [ 0.31   0.3
           0.9  0.15 ]
myin = [1.0,0.1]
p1 = H.KernelExp(myτ)
ps1 = H.PopulationState(p1,2)
ntw = H.RecurrentNetwork(ps1,mywmat,myin)

# ## Start the simulation
# the function `run_simulation!(...)` has been defined above
# Note that `n_spikes` is the total number of spikes
# among **all** units in the system.

n_spikes = 500_000

Tmax = run_simulation!(ntw,n_spikes,100.0,10.0);

# ## Check the rates
# The analytic rate is from  Eq between 6 and  7 in Hawkes 1971

num_rates = H.numerical_rates(ps1)
myspikes_both = ps1.trains_history

ratefou = let G0 =  mywmat .* H.interaction_kernel_fourier(0,p1)
  inv(I-G0)*myin |> real
end
rate_analytic = inv(I-mywmat)*myin

@info "Total duration $(round(Tmax;digits=1)) s"
@info "Rates are $(round.(num_rates;digits=2))"
@info "Analytic rates are $(round.(rate_analytic;digits=2)) Hz"

##    #src
# ## Covariance density
# there are 4 combinations, therefore I will compare 4 lines.

mydt = 0.1
myτmax = 15.0
mytaus = H.get_times(mydt,myτmax)
ntaus = length(mytaus)
cov_num = H.covariance_density_numerical(myspikes_both,mydt,myτmax)

function doplot()
  plt=plot(xlabel="time delay (s)",ylabel="Covariance density")
  for i in 1:2, j in 1:2
    plot!(plt,mytaus[2:end-1],cov_num[i,j,2:end-1], linewidth = 3, label="cov $i-$j")
  end
  return plt
end

doplot()

# The analytic solution is eq 12 from Hawkes 1971

function four_high_res(dt::Real,Tmax::Real)
  k1 = 2
  k2 = 0.005
  myτmax,mydt = Tmax * k1, dt*k2
  mytaus = H.get_times(mydt,myτmax)
  nkeep = div(length(mytaus),k1)
  myfreq = H.get_frequencies_centerzero(mydt,myτmax)
  G_omega = map(mywmat) do w
    ifftshift( w .* H.interaction_kernel_fourier.(myfreq,Ref(p1)))
  end
  D = Diagonal(ratefou)
  M = Array{ComplexF64}(undef,2,2,length(myfreq))
  Mt = similar(M,Float64)
  for i in eachindex(myfreq)
    G = getindex.(G_omega,i)
    M[:,:,i] = (I-G)\D/(I-G')
  end
  for i in 1:2,j in 1:2
    Mt[i,j,:] = real.(ifft(M[i,j,:]))
    Mt[i,j,2:end] ./= mydt # diagonal of t=0 contains the rate
  end
  return mytaus[1:nkeep],Mt[:,:,1:nkeep]
end

taush,Cfou=four_high_res(mydt,myτmax)

function oneplot(i,j)
  plt=plot(xlabel="time delay (s)",ylabel="Covariance density",title="cov $i - $j")
  plot!(plt,mytaus[2:end],cov_num[i,j,2:end] ; linewidth = 3, label="simulation")
  plot!(plt,taush[2:end],Cfou[i,j,2:end]; linestyle=:dash, linewidth=3, label="analytic")
end

# 1
oneplot(1,1)

# 2
oneplot(1,2)

# 3
oneplot(2,1)

# 4
oneplot(2,2)

# **THE END**

# Literate.markdown("examples/exp_1and2D.jl","docs/src";documenter=true,repo_root_url="https://github.com/dylanfesta/HawkesSimulator.jl/blob/master") #src


text:
#=
This file contains the definition of the generalized dimension and related concepts
and also uses linear_regions.jl
=#
export generalized_dim

"""
    generalized_dim(dataset [, sizes]; q = 1, base = MathConstants.e) -> Δ_q
Return the `q` order generalized dimension of the `dataset`, by calculating
the [`genentropy`](@ref) for each `ε ∈ sizes`.

The case of `q = 0` is often called "capacity" or "box-counting" dimension, while
`q = 1` is the "information" dimension.

## Description
The returned dimension is approximated by the
(inverse) power law exponent of the scaling of the [`genentropy`](@ref) ``H_q``
versus the box size `ε`, where `ε ∈ sizes`:

```math
H_q \\sim -\\Delta_q\\log(\\varepsilon)
```

Calling this function performs a lot of automated steps:

  1. A vector of box sizes is decided by calling `sizes = estimate_boxsizes(dataset)`,
     if `sizes` is not given.
  2. For each element of `sizes` the appropriate entropy is
     calculated, through `H = genentropy.(Ref(dataset), sizes; q, base)`.
     Let `x = -log.(sizes)`.
  3. The curve `H(x)` is decomposed into linear regions,
     using [`linear_regions`](@ref)`(x, h)`.
  4. The biggest linear region is chosen, and a fit for the slope of that
     region is performed using the function [`linear_region`](@ref),
     which does a simple linear regression fit using [`linreg`](@ref).
     This slope is the return value of `generalized_dim`.

By doing these steps one by one yourself, you can adjust the keyword arguments
given to each of these function calls, refining the accuracy of the result.
"""
function generalized_dim(data::AbstractDataset, sizes = estimate_boxsizes(data);
        base = Base.MathConstants.e, q = 1.0
    )
    dd = [genentropy(data, ε; q, base) for ε ∈ sizes]
    return linear_region(-log.(base, sizes), dd)[2]
end

ext: jl


text:
# Julia wrapper for header: ogr_core.h
# Automatically generated using Clang.jl

"""
    OGRMalloc(size_t size) -> void *
"""
function ogrmalloc()
    aftercare(ccall((:OGRMalloc, libgdal), Ptr{Cvoid}, ()))
end

"""
    OGRCalloc(size_t count,
              size_t size) -> void *
"""
function ogrcalloc()
    aftercare(ccall((:OGRCalloc, libgdal), Ptr{Cvoid}, ()))
end

"""
    OGRRealloc(void * pOld,
               size_t size) -> void *
"""
function ogrrealloc(arg1, size_t)
    aftercare(ccall((:OGRRealloc, libgdal), Ptr{Cvoid}, (Ptr{Cvoid}, Cint), arg1, size_t))
end

function ogrstrdup(arg1)
    aftercare(ccall((:OGRStrdup, libgdal), Cstring, (Cstring,), arg1), false)
end

"""
    OGRFree(void * pMemory) -> void
"""
function ogrfree(arg1)
    aftercare(ccall((:OGRFree, libgdal), Cvoid, (Ptr{Cvoid},), arg1))
end

"""
    OGRGeometryTypeToName(OGRwkbGeometryType eType) -> const char *

Fetch a human readable name corresponding to an OGRwkbGeometryType value.

### Parameters
* **eType**: the geometry type.

### Returns
internal human readable string, or NULL on failure.
"""
function ogrgeometrytypetoname(eType)
    aftercare(ccall((:OGRGeometryTypeToName, libgdal), Cstring, (OGRwkbGeometryType,), eType), false)
end

"""
    OGRMergeGeometryTypes(OGRwkbGeometryType eMain,
                          OGRwkbGeometryType eExtra) -> OGRwkbGeometryType

Find common geometry type.

### Parameters
* **eMain**: the first input geometry type.
* **eExtra**: the second input geometry type.

### Returns
the merged geometry type.
"""
function ogrmergegeometrytypes(eMain, eExtra)
    aftercare(ccall((:OGRMergeGeometryTypes, libgdal), OGRwkbGeometryType, (OGRwkbGeometryType, OGRwkbGeometryType), eMain, eExtra))
end

"""
    OGRMergeGeometryTypesEx(OGRwkbGeometryType eMain,
                            OGRwkbGeometryType eExtra,
                            int bAllowPromotingToCurves) -> OGRwkbGeometryType

Find common geometry type.

### Parameters
* **eMain**: the first input geometry type.
* **eExtra**: the second input geometry type.
* **bAllowPromotingToCurves**: determine if promotion to curve type must be done.

### Returns
the merged geometry type.
"""
function ogrmergegeometrytypesex(eMain, eExtra, bAllowPromotingToCurves)
    aftercare(ccall((:OGRMergeGeometryTypesEx, libgdal), OGRwkbGeometryType, (OGRwkbGeometryType, OGRwkbGeometryType, Cint), eMain, eExtra, bAllowPromotingToCurves))
end

"""
    OGR_GT_Flatten(OGRwkbGeometryType eType) -> OGRwkbGeometryType

Returns the 2D geometry type corresponding to the passed geometry type.

### Parameters
* **eType**: Input geometry type

### Returns
2D geometry type corresponding to the passed geometry type.
"""
function ogr_gt_flatten(eType)
    aftercare(ccall((:OGR_GT_Flatten, libgdal), OGRwkbGeometryType, (OGRwkbGeometryType,), eType))
end

"""
    OGR_GT_SetZ(OGRwkbGeometryType eType) -> OGRwkbGeometryType

Returns the 3D geometry type corresponding to the passed geometry type.

### Parameters
* **eType**: Input geometry type

### Returns
3D geometry type corresponding to the passed geometry type.
"""
function ogr_gt_setz(eType)
    aftercare(ccall((:OGR_GT_SetZ, libgdal), OGRwkbGeometryType, (OGRwkbGeometryType,), eType))
end

"""
    OGR_GT_SetM(OGRwkbGeometryType eType) -> OGRwkbGeometryType

Returns the measured geometry type corresponding to the passed geometry type.

### Parameters
* **eType**: Input geometry type

### Returns
measured geometry type corresponding to the passed geometry type.
"""
function ogr_gt_setm(eType)
    aftercare(ccall((:OGR_GT_SetM, libgdal), OGRwkbGeometryType, (OGRwkbGeometryType,), eType))
end

"""
    OGR_GT_SetModifier(OGRwkbGeometryType eType,
                       int bHasZ,
                       int bHasM) -> OGRwkbGeometryType

Returns a XY, XYZ, XYM or XYZM geometry type depending on parameter.

### Parameters
* **eType**: Input geometry type
* **bHasZ**: TRUE if the output geometry type must be 3D.
* **bHasM**: TRUE if the output geometry type must be measured.

### Returns
Output geometry type.
"""
function ogr_gt_setmodifier(eType, bSetZ, bSetM)
    aftercare(ccall((:OGR_GT_SetModifier, libgdal), OGRwkbGeometryType, (OGRwkbGeometryType, Cint, Cint), eType, bSetZ, bSetM))
end

"""
    OGR_GT_HasZ(OGRwkbGeometryType eType) -> int

Return if the geometry type is a 3D geometry type.

### Parameters
* **eType**: Input geometry type

### Returns
TRUE if the geometry type is a 3D geometry type.
"""
function ogr_gt_hasz(eType)
    aftercare(ccall((:OGR_GT_HasZ, libgdal), Cint, (OGRwkbGeometryType,), eType))
end

"""
    OGR_GT_HasM(OGRwkbGeometryType eType) -> int

Return if the geometry type is a measured type.

### Parameters
* **eType**: Input geometry type

### Returns
TRUE if the geometry type is a measured type.
"""
function ogr_gt_hasm(eType)
    aftercare(ccall((:OGR_GT_HasM, libgdal), Cint, (OGRwkbGeometryType,), eType))
end

"""
    OGR_GT_IsSubClassOf(OGRwkbGeometryType eType,
                        OGRwkbGeometryType eSuperType) -> int

Returns if a type is a subclass of another one.

### Parameters
* **eType**: Type.
* **eSuperType**: Super type

### Returns
TRUE if eType is a subclass of eSuperType.
"""
function ogr_gt_issubclassof(eType, eSuperType)
    aftercare(ccall((:OGR_GT_IsSubClassOf, libgdal), Cint, (OGRwkbGeometryType, OGRwkbGeometryType), eType, eSuperType))
end

"""
    OGR_GT_IsCurve(OGRwkbGeometryType eGeomType) -> int

Return if a geometry type is an instance of Curve.

### Parameters
* **eGeomType**: the geometry type

### Returns
TRUE if the geometry type is an instance of Curve
"""
function ogr_gt_iscurve(arg1)
    aftercare(ccall((:OGR_GT_IsCurve, libgdal), Cint, (OGRwkbGeometryType,), arg1))
end

"""
    OGR_GT_IsSurface(OGRwkbGeometryType eGeomType) -> int

Return if a geometry type is an instance of Surface.

### Parameters
* **eGeomType**: the geometry type

### Returns
TRUE if the geometry type is an instance of Surface
"""
function ogr_gt_issurface(arg1)
    aftercare(ccall((:OGR_GT_IsSurface, libgdal), Cint, (OGRwkbGeometryType,), arg1))
end

"""
    OGR_GT_IsNonLinear(OGRwkbGeometryType eGeomType) -> int

Return if a geometry type is a non-linear geometry type.

### Parameters
* **eGeomType**: the geometry type

### Returns
TRUE if the geometry type is a non-linear geometry type.
"""
function ogr_gt_isnonlinear(arg1)
    aftercare(ccall((:OGR_GT_IsNonLinear, libgdal), Cint, (OGRwkbGeometryType,), arg1))
end

"""
    OGR_GT_GetCollection(OGRwkbGeometryType eType) -> OGRwkbGeometryType

Returns the collection type that can contain the passed geometry type.

### Parameters
* **eType**: Input geometry type

### Returns
the collection type that can contain the passed geometry type or wkbUnknown
"""
function ogr_gt_getcollection(eType)
    aftercare(ccall((:OGR_GT_GetCollection, libgdal), OGRwkbGeometryType, (OGRwkbGeometryType,), eType))
end

"""
    OGR_GT_GetCurve(OGRwkbGeometryType eType) -> OGRwkbGeometryType

Returns the curve geometry type that can contain the passed geometry type.

### Parameters
* **eType**: Input geometry type

### Returns
the curve type that can contain the passed geometry type
"""
function ogr_gt_getcurve(eType)
    aftercare(ccall((:OGR_GT_GetCurve, libgdal), OGRwkbGeometryType, (OGRwkbGeometryType,), eType))
end

"""
    OGR_GT_GetLinear(OGRwkbGeometryType eType) -> OGRwkbGeometryType

Returns the non-curve geometry type that can contain the passed geometry type.

### Parameters
* **eType**: Input geometry type

### Returns
the non-curve type that can contain the passed geometry type
"""
function ogr_gt_getlinear(eType)
    aftercare(ccall((:OGR_GT_GetLinear, libgdal), OGRwkbGeometryType, (OGRwkbGeometryType,), eType))
end

"""
    OGRParseDate(const char * pszInput,
                 OGRField * psField,
                 int nOptions) -> int

Parse date string.

### Parameters
* **pszInput**: the input date string.
* **psField**: the OGRField that will be updated with the parsed result.
* **nOptions**: parsing options, for now always 0.

### Returns
TRUE if apparently successful or FALSE on failure.
"""
function ogrparsedate(pszInput, psOutput, nOptions)
    aftercare(ccall((:OGRParseDate, libgdal), Cint, (Cstring, Ptr{OGRField}, Cint), pszInput, psOutput, nOptions))
end

"""
    GDALVersionInfo(const char * pszRequest) -> const char *

Get runtime version information.

### Parameters
* **pszRequest**: the type of version info desired, as listed above.

### Returns
an internal string containing the requested information.
"""
function gdalversioninfo(arg1)
    aftercare(ccall((:GDALVersionInfo, libgdal), Cstring, (Cstring,), arg1), false)
end

"""
    GDALCheckVersion(int nVersionMajor,
                     int nVersionMinor,
                     const char * pszCallingComponentName) -> int

Return TRUE if GDAL library version at runtime matches nVersionMajor.nVersionMinor.

### Parameters
* **nVersionMajor**: Major version to be tested against
* **nVersionMinor**: Minor version to be tested against
* **pszCallingComponentName**: If not NULL, in case of version mismatch, the method will issue a failure mentioning the name of the calling component.

### Returns
TRUE if GDAL library version at runtime matches nVersionMajor.nVersionMinor, FALSE otherwise.
"""
function gdalcheckversion(nVersionMajor, nVersionMinor, pszCallingComponentName)
    aftercare(ccall((:GDALCheckVersion, libgdal), Cint, (Cint, Cint, Cstring), nVersionMajor, nVersionMinor, pszCallingComponentName))
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "zfp"
version = v"0.5.5"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/LLNL/zfp/releases/download/0.5.5/zfp-0.5.5.tar.gz", "fdf7b948bab1f4e5dccfe2c2048fd98c24e417ad8fb8a51ed3463d04147393c5")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/zfp-*
mkdir build && cd build/
if [[ "${target}" == *-freebsd* ]]; then
    # Give a hint to the linker to explicitly use `-lm`
    export LDFLAGS="-lm"
fi
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release ..
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libzfp", :libzfp)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# Goal: Make sure some combinations that raise errors do so
using Test
using OrdinaryDiffEq
using LinearAlgebra
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems;
importodeproblems();
import DiffEqProblemLibrary.ODEProblemLibrary:
    prob_ode_linear, prob_ode_2Dlinear, prob_ode_lotkavoltera, prob_ode_fitzhughnagumo

import ProbNumDiffEq: remake_prob_with_jac

@testset "Fixed-timestep requires dt" begin
    prob = prob_ode_lotkavoltera
    @test_throws ErrorException solve(prob, EK0(), adaptive=false)
    @test solve(prob, EK0(), adaptive=false, dt=0.05) isa ProbNumDiffEq.ProbODESolution
end

@testset "`dense=true` requires `smooth=true`" begin
    prob = prob_ode_lotkavoltera
    @test_throws ErrorException solve(prob, EK0(smooth=false))
end

@testset "`save_everystep=false` requires `smooth=false`" begin
    prob = prob_ode_lotkavoltera
    @test_throws ErrorException solve(prob, EK0(smooth=true), save_everystep=false)
end

@testset "`dense=false` warns if `smooth=true`" begin
    prob = prob_ode_lotkavoltera
    @test_logs (
        :warn,
        "If you set dense=false for efficiency, you might also want to set smooth=false.",
    ) solve(prob, EK0(smooth=true), dense=false)
end


text:
#=
2D advection on unstructured grid using FV
=#

### If the Finch package has already been added, use this line #########
using Finch # Note: to add the package, first do: ]add "https://github.com/paralab/Finch.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Finch)
#     include("../Finch.jl");
#     using .Finch
# end
##########################################################################

init_finch("FVadvection2d");

useLog("FVadvection2dlog")

# Configuration setup
domain(2)
solverType(FV)
timeStepper(RK4)

# Mesh
# This rectangle covers [0, 0.1]x[0, 0.3]
# Uncomment the desired mesh.
mesh("utriangle.msh")  # Using triangles
#mesh("uquad.msh")     # Using quads

add_boundary_ID(2, (x,y) -> (x >= 0.1));
add_boundary_ID(3, (x,y) -> (y <= 0));
add_boundary_ID(4, (x,y) -> (y >= 0.3));

# Variables and BCs
u = variable("u", location=CELL)
boundary(u, 1, FLUX, 0) # x=0
boundary(u, 2, NO_BC) # x=1
boundary(u, 3, FLUX, "sin(10*pi*x)*sin(2*pi*t)^2") # y=0
boundary(u, 4, NO_BC) # y=1

# Time interval and initial condition
T = 0.7;
timeInterval(T)
initial(u, "0")

# The flux and source terms of the conservation equation
# F and S in the following equation:
# Dt(int(u dx)) = int(S dx) - int(F.n ds)
coefficient("a", [0.1, 1], type=VECTOR) # advection velocity
# The "upwind" function applies upwinding to the term (a.n)*u with flow velocity a.
# The optional third parameter is for tuning. Default upwind = 0, central = 1. Choose something between these.
flux(u, "upwind(a,u)")
# Note that there is no source() for this problem

#exportCode("fvad2dcode") # uncomment to export generated code to a file
#importCode("fvad2dcode") # uncomment to import code from a file

solve(u)

finalize_finch()

##### Uncomment below to plot

# xy = Finch.fv_info.cellCenters

# using Plots
# pyplot();
# display(plot(xy[1,:], xy[2,:], u.values[:], st=:surface))

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
const N_GRADS_4D_EXPONENT = 9
const N_GRADS_4D = 1 << N_GRADS_4D_EXPONENT

const NORMALIZER_4D = 0.11127401889945551

struct Grad4
    x::Float32
    y::Float32
    z::Float32
    w::Float32
end

const grad4 = Grad4[
 Grad4(-0.6740059517812944,   -0.3239847771997537,   -0.3239847771997537,    0.5794684678643381),
 Grad4(-0.7504883828755602,   -0.4004672082940195,    0.15296486218853164,   0.5029860367700724),
 Grad4(-0.7504883828755602,    0.15296486218853164,  -0.4004672082940195,    0.5029860367700724),
 Grad4(-0.8828161875373585,    0.08164729285680945,   0.08164729285680945,   0.4553054119602712),
 Grad4(-0.4553054119602712,   -0.08164729285680945,  -0.08164729285680945,   0.8828161875373585),
 Grad4(-0.5029860367700724,   -0.15296486218853164,   0.4004672082940195,    0.7504883828755602),
 Grad4(-0.5029860367700724,    0.4004672082940195,   -0.15296486218853164,   0.7504883828755602),
 Grad4(-0.5794684678643381,    0.3239847771997537,    0.3239847771997537,    0.6740059517812944),
 Grad4(-0.6740059517812944,   -0.3239847771997537,    0.5794684678643381,   -0.3239847771997537),
 Grad4(-0.7504883828755602,   -0.4004672082940195,    0.5029860367700724,    0.15296486218853164),
 Grad4(-0.7504883828755602,    0.15296486218853164,   0.5029860367700724,   -0.4004672082940195),
 Grad4(-0.8828161875373585,    0.08164729285680945,   0.4553054119602712,    0.08164729285680945),
 Grad4(-0.4553054119602712,   -0.08164729285680945,   0.8828161875373585,   -0.08164729285680945),
 Grad4(-0.5029860367700724,   -0.15296486218853164,   0.7504883828755602,    0.4004672082940195),
 Grad4(-0.5029860367700724,    0.4004672082940195,    0.7504883828755602,   -0.15296486218853164),
 Grad4(-0.5794684678643381,    0.3239847771997537,    0.6740059517812944,    0.3239847771997537),
 Grad4(-0.6740059517812944,    0.5794684678643381,   -0.3239847771997537,   -0.3239847771997537),
 Grad4(-0.7504883828755602,    0.5029860367700724,   -0.4004672082940195,    0.15296486218853164),
 Grad4(-0.7504883828755602,    0.5029860367700724,    0.15296486218853164,  -0.4004672082940195),
 Grad4(-0.8828161875373585,    0.4553054119602712,    0.08164729285680945,   0.08164729285680945),
 Grad4(-0.4553054119602712,    0.8828161875373585,   -0.08164729285680945,  -0.08164729285680945),
 Grad4(-0.5029860367700724,    0.7504883828755602,   -0.15296486218853164,   0.4004672082940195),
 Grad4(-0.5029860367700724,    0.7504883828755602,    0.4004672082940195,   -0.15296486218853164),
 Grad4(-0.5794684678643381,    0.6740059517812944,    0.3239847771997537,    0.3239847771997537),
 Grad4( 0.5794684678643381,   -0.6740059517812944,   -0.3239847771997537,   -0.3239847771997537),
 Grad4( 0.5029860367700724,   -0.7504883828755602,   -0.4004672082940195,    0.15296486218853164),
 Grad4( 0.5029860367700724,   -0.7504883828755602,    0.15296486218853164,  -0.4004672082940195),
 Grad4( 0.4553054119602712,   -0.8828161875373585,    0.08164729285680945,   0.08164729285680945),
 Grad4( 0.8828161875373585,   -0.4553054119602712,   -0.08164729285680945,  -0.08164729285680945),
 Grad4( 0.7504883828755602,   -0.5029860367700724,   -0.15296486218853164,   0.4004672082940195),
 Grad4( 0.7504883828755602,   -0.5029860367700724,    0.4004672082940195,   -0.15296486218853164),
 Grad4( 0.6740059517812944,   -0.5794684678643381,    0.3239847771997537,    0.3239847771997537),
 #-----------------------------------------------------------------------------------------------#
 Grad4(-0.753341017856078,    -0.37968289875261624,  -0.37968289875261624,  -0.37968289875261624),
 Grad4(-0.7821684431180708,   -0.4321472685365301,   -0.4321472685365301,    0.12128480194602098),
 Grad4(-0.7821684431180708,   -0.4321472685365301,    0.12128480194602098,  -0.4321472685365301),
 Grad4(-0.7821684431180708,    0.12128480194602098,  -0.4321472685365301,   -0.4321472685365301),
 Grad4(-0.8586508742123365,   -0.508629699630796,     0.044802370851755174,  0.044802370851755174),
 Grad4(-0.8586508742123365,    0.044802370851755174, -0.508629699630796,     0.044802370851755174),
 Grad4(-0.8586508742123365,    0.044802370851755174,  0.044802370851755174, -0.508629699630796),
 Grad4(-0.9982828964265062,   -0.03381941603233842,  -0.03381941603233842,  -0.03381941603233842),
 Grad4(-0.37968289875261624,  -0.753341017856078,    -0.37968289875261624,  -0.37968289875261624),
 Grad4(-0.4321472685365301,   -0.7821684431180708,   -0.4321472685365301,    0.12128480194602098),
 Grad4(-0.4321472685365301,   -0.7821684431180708,    0.12128480194602098,  -0.4321472685365301),
 Grad4( 0.12128480194602098,  -0.7821684431180708,   -0.4321472685365301,   -0.4321472685365301),
 Grad4(-0.508629699630796,    -0.8586508742123365,    0.044802370851755174,  0.044802370851755174),
 Grad4( 0.044802370851755174, -0.8586508742123365,   -0.508629699630796,     0.044802370851755174),
 Grad4( 0.044802370851755174, -0.8586508742123365,    0.044802370851755174, -0.508629699630796),
 Grad4(-0.03381941603233842,  -0.9982828964265062,   -0.03381941603233842,  -0.03381941603233842),
 Grad4(-0.37968289875261624,  -0.37968289875261624,  -0.753341017856078,    -0.37968289875261624),
 Grad4(-0.4321472685365301,   -0.4321472685365301,   -0.7821684431180708,    0.12128480194602098),
 Grad4(-0.4321472685365301,    0.12128480194602098,  -0.7821684431180708,   -0.4321472685365301),
 Grad4( 0.12128480194602098,  -0.4321472685365301,   -0.7821684431180708,   -0.4321472685365301),
 Grad4(-0.508629699630796,     0.044802370851755174, -0.8586508742123365,    0.044802370851755174),
 Grad4( 0.044802370851755174, -0.508629699630796,    -0.8586508742123365,    0.044802370851755174),
 Grad4( 0.044802370851755174,  0.044802370851755174, -0.8586508742123365,   -0.508629699630796),
 Grad4(-0.03381941603233842,  -0.03381941603233842,  -0.9982828964265062,   -0.03381941603233842),
 Grad4(-0.37968289875261624,  -0.37968289875261624,  -0.37968289875261624,  -0.753341017856078),
 Grad4(-0.4321472685365301,   -0.4321472685365301,    0.12128480194602098,  -0.7821684431180708),
 Grad4(-0.4321472685365301,    0.12128480194602098,  -0.4321472685365301,   -0.7821684431180708),
 Grad4( 0.12128480194602098,  -0.4321472685365301,   -0.4321472685365301,   -0.7821684431180708),
 Grad4(-0.508629699630796,     0.044802370851755174,  0.044802370851755174, -0.8586508742123365),
 Grad4( 0.044802370851755174, -0.508629699630796,     0.044802370851755174, -0.8586508742123365),
 Grad4( 0.044802370851755174,  0.044802370851755174, -0.508629699630796,    -0.8586508742123365),
 Grad4(-0.03381941603233842,  -0.03381941603233842,  -0.03381941603233842,  -0.9982828964265062),
 Grad4(-0.3239847771997537,   -0.6740059517812944,   -0.3239847771997537,    0.5794684678643381),
 Grad4(-0.4004672082940195,   -0.7504883828755602,    0.15296486218853164,   0.5029860367700724),
 Grad4( 0.15296486218853164,  -0.7504883828755602,   -0.4004672082940195,    0.5029860367700724),
 Grad4( 0.08164729285680945,  -0.8828161875373585,    0.08164729285680945,   0.4553054119602712),
 Grad4(-0.08164729285680945,  -0.4553054119602712,   -0.08164729285680945,   0.8828161875373585),
 Grad4(-0.15296486218853164,  -0.5029860367700724,    0.4004672082940195,    0.7504883828755602),
 Grad4( 0.4004672082940195,   -0.5029860367700724,   -0.15296486218853164,   0.7504883828755602),
 Grad4( 0.3239847771997537,   -0.5794684678643381,    0.3239847771997537,    0.6740059517812944),
 Grad4(-0.3239847771997537,   -0.3239847771997537,   -0.6740059517812944,    0.5794684678643381),
 Grad4(-0.4004672082940195,    0.15296486218853164,  -0.7504883828755602,    0.5029860367700724),
 Grad4( 0.15296486218853164,  -0.4004672082940195,   -0.7504883828755602,    0.5029860367700724),
 Grad4( 0.08164729285680945,   0.08164729285680945,  -0.8828161875373585,    0.4553054119602712),
 Grad4(-0.08164729285680945,  -0.08164729285680945,  -0.4553054119602712,    0.8828161875373585),
 Grad4(-0.15296486218853164,   0.4004672082940195,   -0.5029860367700724,    0.7504883828755602),
 Grad4( 0.4004672082940195,   -0.15296486218853164,  -0.5029860367700724,    0.7504883828755602),
 Grad4( 0.3239847771997537,    0.3239847771997537,   -0.5794684678643381,    0.6740059517812944),
 Grad4(-0.3239847771997537,   -0.6740059517812944,    0.5794684678643381,   -0.3239847771997537),
 Grad4(-0.4004672082940195,   -0.7504883828755602,    0.5029860367700724,    0.15296486218853164),
 Grad4( 0.15296486218853164,  -0.7504883828755602,    0.5029860367700724,   -0.4004672082940195),
 Grad4( 0.08164729285680945,  -0.8828161875373585,    0.4553054119602712,    0.08164729285680945),
 Grad4(-0.08164729285680945,  -0.4553054119602712,    0.8828161875373585,   -0.08164729285680945),
 Grad4(-0.15296486218853164,  -0.5029860367700724,    0.7504883828755602,    0.4004672082940195),
 Grad4( 0.4004672082940195,   -0.5029860367700724,    0.7504883828755602,   -0.15296486218853164),
 Grad4( 0.3239847771997537,   -0.5794684678643381,    0.6740059517812944,    0.3239847771997537),
 Grad4(-0.3239847771997537,   -0.3239847771997537,    0.5794684678643381,   -0.6740059517812944),
 Grad4(-0.4004672082940195,    0.15296486218853164,   0.5029860367700724,   -0.7504883828755602),
 Grad4( 0.15296486218853164,  -0.4004672082940195,    0.5029860367700724,   -0.7504883828755602),
 Grad4( 0.08164729285680945,   0.08164729285680945,   0.4553054119602712,   -0.8828161875373585),
 Grad4(-0.08164729285680945,  -0.08164729285680945,   0.8828161875373585,   -0.4553054119602712),
 Grad4(-0.15296486218853164,   0.4004672082940195,    0.7504883828755602,   -0.5029860367700724),
 Grad4( 0.4004672082940195,   -0.15296486218853164,   0.7504883828755602,   -0.5029860367700724),
 Grad4( 0.3239847771997537,    0.3239847771997537,    0.6740059517812944,   -0.5794684678643381),
 Grad4(-0.3239847771997537,    0.5794684678643381,   -0.6740059517812944,   -0.3239847771997537),
 Grad4(-0.4004672082940195,    0.5029860367700724,   -0.7504883828755602,    0.15296486218853164),
 Grad4( 0.15296486218853164,   0.5029860367700724,   -0.7504883828755602,   -0.4004672082940195),
 Grad4( 0.08164729285680945,   0.4553054119602712,   -0.8828161875373585,    0.08164729285680945),
 Grad4(-0.08164729285680945,   0.8828161875373585,   -0.4553054119602712,   -0.08164729285680945),
 Grad4(-0.15296486218853164,   0.7504883828755602,   -0.5029860367700724,    0.4004672082940195),
 Grad4( 0.4004672082940195,    0.7504883828755602,   -0.5029860367700724,   -0.15296486218853164),
 Grad4( 0.3239847771997537,    0.6740059517812944,   -0.5794684678643381,    0.3239847771997537),
 Grad4(-0.3239847771997537,    0.5794684678643381,   -0.3239847771997537,   -0.6740059517812944),
 Grad4(-0.4004672082940195,    0.5029860367700724,    0.15296486218853164,  -0.7504883828755602),
 Grad4( 0.15296486218853164,   0.5029860367700724,   -0.4004672082940195,   -0.7504883828755602),
 Grad4( 0.08164729285680945,   0.4553054119602712,    0.08164729285680945,  -0.8828161875373585),
 Grad4(-0.08164729285680945,   0.8828161875373585,   -0.08164729285680945,  -0.4553054119602712),
 Grad4(-0.15296486218853164,   0.7504883828755602,    0.4004672082940195,   -0.5029860367700724),
 Grad4( 0.4004672082940195,    0.7504883828755602,   -0.15296486218853164,  -0.5029860367700724),
 Grad4( 0.3239847771997537,    0.6740059517812944,    0.3239847771997537,   -0.5794684678643381),
 Grad4( 0.5794684678643381,   -0.3239847771997537,   -0.6740059517812944,   -0.3239847771997537),
 Grad4( 0.5029860367700724,   -0.4004672082940195,   -0.7504883828755602,    0.15296486218853164),
 Grad4( 0.5029860367700724,    0.15296486218853164,  -0.7504883828755602,   -0.4004672082940195),
 Grad4( 0.4553054119602712,    0.08164729285680945,  -0.8828161875373585,    0.08164729285680945),
 Grad4( 0.8828161875373585,   -0.08164729285680945,  -0.4553054119602712,   -0.08164729285680945),
 Grad4( 0.7504883828755602,   -0.15296486218853164,  -0.5029860367700724,    0.4004672082940195),
 Grad4( 0.7504883828755602,    0.4004672082940195,   -0.5029860367700724,   -0.15296486218853164),
 Grad4( 0.6740059517812944,    0.3239847771997537,   -0.5794684678643381,    0.3239847771997537),
 Grad4( 0.5794684678643381,   -0.3239847771997537,   -0.3239847771997537,   -0.6740059517812944),
 Grad4( 0.5029860367700724,   -0.4004672082940195,    0.15296486218853164,  -0.7504883828755602),
 Grad4( 0.5029860367700724,    0.15296486218853164,  -0.4004672082940195,   -0.7504883828755602),
 Grad4( 0.4553054119602712,    0.08164729285680945,   0.08164729285680945,  -0.8828161875373585),
 Grad4( 0.8828161875373585,   -0.08164729285680945,  -0.08164729285680945,  -0.4553054119602712),
 Grad4( 0.7504883828755602,   -0.15296486218853164,   0.4004672082940195,   -0.5029860367700724),
 Grad4( 0.7504883828755602,    0.4004672082940195,   -0.15296486218853164,  -0.5029860367700724),
 Grad4( 0.6740059517812944,    0.3239847771997537,    0.3239847771997537,   -0.5794684678643381),
 Grad4( 0.03381941603233842,   0.03381941603233842,   0.03381941603233842,   0.9982828964265062),
 Grad4(-0.044802370851755174, -0.044802370851755174,  0.508629699630796,     0.8586508742123365),
 Grad4(-0.044802370851755174,  0.508629699630796,    -0.044802370851755174,  0.8586508742123365),
 Grad4(-0.12128480194602098,   0.4321472685365301,    0.4321472685365301,    0.7821684431180708),
 Grad4( 0.508629699630796,    -0.044802370851755174, -0.044802370851755174,  0.8586508742123365),
 Grad4( 0.4321472685365301,   -0.12128480194602098,   0.4321472685365301,    0.7821684431180708),
 Grad4( 0.4321472685365301,    0.4321472685365301,   -0.12128480194602098,   0.7821684431180708),
 Grad4( 0.37968289875261624,   0.37968289875261624,   0.37968289875261624,   0.753341017856078),
 Grad4( 0.03381941603233842,   0.03381941603233842,   0.9982828964265062,    0.03381941603233842),
 Grad4(-0.044802370851755174,  0.044802370851755174,  0.8586508742123365,    0.508629699630796),
 Grad4(-0.044802370851755174,  0.508629699630796,     0.8586508742123365,   -0.044802370851755174),
 Grad4(-0.12128480194602098,   0.4321472685365301,    0.7821684431180708,    0.4321472685365301),
 Grad4( 0.508629699630796,    -0.044802370851755174,  0.8586508742123365,   -0.044802370851755174),
 Grad4( 0.4321472685365301,   -0.12128480194602098,   0.7821684431180708,    0.4321472685365301),
 Grad4( 0.4321472685365301,    0.4321472685365301,    0.7821684431180708,   -0.12128480194602098),
 Grad4( 0.37968289875261624,   0.37968289875261624,   0.753341017856078,     0.37968289875261624),
 Grad4( 0.03381941603233842,   0.9982828964265062,    0.03381941603233842,   0.03381941603233842),
 Grad4(-0.044802370851755174,  0.8586508742123365,   -0.044802370851755174,  0.508629699630796),
 Grad4(-0.044802370851755174,  0.8586508742123365,    0.508629699630796,    -0.044802370851755174),
 Grad4(-0.12128480194602098,   0.7821684431180708,    0.4321472685365301,    0.4321472685365301),
 Grad4( 0.508629699630796,     0.8586508742123365,   -0.044802370851755174, -0.044802370851755174),
 Grad4( 0.4321472685365301,    0.7821684431180708,   -0.12128480194602098,   0.4321472685365301),
 Grad4( 0.4321472685365301,    0.7821684431180708,    0.4321472685365301,   -0.12128480194602098),
 Grad4( 0.37968289875261624,   0.753341017856078,     0.37968289875261624,   0.37968289875261624),
 Grad4( 0.9982828964265062,    0.03381941603233842,   0.03381941603233842,   0.03381941603233842),
 Grad4( 0.8586508742123365,   -0.044802370851755174, -0.044802370851755174,  0.508629699630796),
 Grad4( 0.8586508742123365,   -0.044802370851755174,  0.508629699630796,    -0.044802370851755174),
 Grad4( 0.7821684431180708,   -0.12128480194602098,   0.4321472685365301,    0.4321472685365301),
 Grad4( 0.8586508742123365,    0.508629699630796,    -0.044802370851755174, -0.044802370851755174),
 Grad4( 0.7821684431180708,    0.4321472685365301,   -0.12128480194602098,   0.4321472685365301),
 Grad4( 0.7821684431180708,    0.4321472685365301,    0.4321472685365301,   -0.12128480194602098),
 Grad4( 0.753341017856078,     0.37968289875261624,   0.37968289875261624,   0.37968289875261624),
]

const GRADIENTS_4D = Vector{Grad4}(undef, N_GRADS_4D)

function grad(a, seed::Int64, xsvp::Int64, ysvp::Int64, zsvp::Int64, wsvp::Int64,
              dx::Float32, dy::Float32, dz::Float32, dw::Float32)
    hash = xor(seed, xsvp, ysvp, zsvp, wsvp) * HASH_MULTIPLIER
    g = GRADIENTS_4D[(xor(hash, hash >> (64 - N_GRADS_4D_EXPONENT)) & (N_GRADS_4D - 1)) + 1]
    Float32(p4(a) * (g.x * dx + g.y * dy + g.z * dz + g.w * dw) / NORMALIZER_4D)
end

function fill_4d()
    j = 0
    len = length(grad4)
    for i = 1:N_GRADS_4D
        GRADIENTS_4D[i] = grad4[j += 1]
        j == len && (j = 0)
    end
end

fill_4d()


text:
# Libraries
using DataFrames, Dates, FileIO, JLD;
using Contour, DecisionTree, Random, StableRNGs, Statistics, MessyTimeSeries, MessyTimeSeriesOptim;
using PGFPlotsX, LaTeXStrings;

function learn_grid_structure(x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64}, rng::AbstractRNG=StableRNG(1), n_trees::Int64=1000)

    # Initial settings
    n = length(x);
    trees = Vector{DecisionTreeRegressor}(undef, n_trees);

    # Bagging
    for i=1:n_trees
        trees[i] = DecisionTreeRegressor(min_samples_leaf=25, rng=rng);
        current_selection = rand(rng, 1:n, n);
        current_x = x[current_selection];
        current_y = y[current_selection];
        current_z = z[current_selection];
        DecisionTree.fit!(trees[i], [current_x current_y], current_z);
    end

    return trees;
end

function interpolate_z(x::Float64, y::Float64, trees::Vector{DecisionTreeRegressor})

    # Initial settings
    output = 0.0;
    n_trees = length(trees);
    covariates = [x y];

    # Bagging prediction
    for tree in trees
        output += DecisionTree.predict(tree, covariates)[1];
    end
    output /= n_trees;

    return output;
end

# Plot: grid
push!(PGFPlotsX.CUSTOM_PREAMBLE,
        raw"\usepgfplotslibrary{colorbrewer}",
        raw"\usepgfplotslibrary{colormaps}",
        raw"\usepgfplotslibrary{patchplots}",
        raw"\pgfplotsset
            {
                tick label style = {font = {\fontsize{12 pt}{12 pt}\selectfont}},
                label style = {font = {\fontsize{12 pt}{12 pt}\selectfont}},
                legend style = {font = {\fontsize{12 pt}{12 pt}\selectfont}},
                title style = {font = {\fontsize{12 pt}{12 pt}\selectfont}},
            }"
)

# GroupPlot axs per row
gpr = Array{Any,1}(undef, 2);
for i=1:2
    gpr[i] = Array{Any,1}(undef, 5);
end

# Models and methods
axs_titles = ["In-sample", "Pseudo out-of-sample", "Block jackknife (c/T=0.10)", "Block jackknife (c/T=0.20)", "Artificial jackknife"];

# BJK specific options
t0 = 52;
adj_bjk = true;

# Interpolated grid
x_grid = range(0.001, stop=1, length=50);
y_grid = range(0.001, stop=20, length=50);

# Loop over each model
for (j, model_prefix) in enumerate(["VAR", "VMA"])

    if model_prefix == "VAR"
        ztick_distance = 5
        zmin_array = [0; 49; 50; 50; 59];
        zmax_array = [70; 70; 65; 63; 68];

    else
        ztick_distance = 10
        zmin_array = [0; 10; 30; 30; 50];
        zmax_array = [120; 100; 70; 60; 125];
    end

    for i in 1:5

        #=
        ----------------------------------------------------------------------------------------------------
        Load input data and interpolate grid
        ----------------------------------------------------------------------------------------------------
        =#

        # Raw output

        if i <= 3
            err_type = i;
        else
            err_type = i-1;
        end

        if err_type != 3
            f = load("./$(model_prefix)_output/err_type_$(err_type).jld");
        else
            f = load("./$(model_prefix)_output/3_$(10+(i-3)*10)pct/err_type_$(err_type).jld");
        end

        candidates = f["candidates"];
        errors = f["errors"];

        if err_type == 3
            if i==3
                bjk_sample = block_jackknife(f["selection_sample"], 0.1);
            else
                bjk_sample = block_jackknife(f["selection_sample"], 0.2);
            end

            T_validation = size(bjk_sample, 2)-t0;
            missings_after_t0 = [sum(sum(ismissing.(bjk_sample[:, t0+1:end, i]), dims=1) .== size(bjk_sample, 1)) for i=1:size(bjk_sample, 3)];
            errors .*= mean(T_validation ./ (T_validation .- missings_after_t0));
        end

        # Explored grid
        x = candidates[3,:];
        y = candidates[2,:] .* (candidates[4,:].^(candidates[1,:].-1));
        z = errors;

        # Learn structure of the grid and interpolate z
        trees = learn_grid_structure(x, y, z);
        z_grid = [interpolate_z(xx, yy, trees) for xx in x_grid, yy in y_grid];

        # Reference points for zmax and zmin in the chart
        z_min_chart = zmin_array[i];
        z_max_chart = zmax_array[i];

        # Meta source: minimum meta_source* is equal to 0
        min_z = min(minimum(z), minimum(z_grid));
        meta_source = z_grid .- min_z;
        meta_source_scatter = z .- min_z;

        # Meta source: maximum meta_source* is equal to 1
        max_meta = max(maximum(meta_source), maximum(meta_source_scatter));
        meta_source ./= max_meta;
        meta_source_scatter ./= max_meta;

        #=
        ----------------------------------------------------------------------------------------------------
        Current subplot settings
        ----------------------------------------------------------------------------------------------------
        =#

        x_dir_settings = "normal";
        view_angle = (45, 15);

        subtitle = ifelse(j==1, axs_titles[i], "");
        axis_on_top = "false";
        colormap_settings = ifelse((j==1) & (i==5), "true", "false");
        z_label_name = ifelse(i==2, L"Loss", "");

        #=
        ----------------------------------------------------------------------------------------------------
        Generate subplot
        ----------------------------------------------------------------------------------------------------
        =#

        subplot_p1 = @pgf Plot3({surf, point_meta="explicit"}, Coordinates(x_grid, y_grid, z_grid, meta=meta_source));
        subplot_p2 = @pgf Plot3({only_marks, scatter, mark_size="1pt", point_meta="explicit"}, Coordinates(x, y, z_min_chart .* ones(length(x)), meta=meta_source_scatter));

        gpr[j][i] = @pgf Axis(
        {
            "axis on top"=axis_on_top,
            height = "10cm",
            width = "10cm",
            grid = "major",
            grid_style = "{darkgray, very thin}",
            tick_style = "{darkgray, very thin}",
            "colormap"="{turbo}{rgb=(0.18995,0.07176,0.23217), rgb=(0.19483,0.08339,0.26149), rgb=(0.19956,0.09498,0.29024), rgb=(0.20415,0.10652,0.31844), rgb=(0.20860,0.11802,0.34607), rgb=(0.21291,0.12947,0.37314), rgb=(0.21708,0.14087,0.39964), rgb=(0.22111,0.15223,0.42558), rgb=(0.22500,0.16354,0.45096), rgb=(0.22875,0.17481,0.47578), rgb=(0.23236,0.18603,0.50004), rgb=(0.23582,0.19720,0.52373), rgb=(0.23915,0.20833,0.54686), rgb=(0.24234,0.21941,0.56942), rgb=(0.24539,0.23044,0.59142), rgb=(0.24830,0.24143,0.61286), rgb=(0.25107,0.25237,0.63374), rgb=(0.25369,0.26327,0.65406), rgb=(0.25618,0.27412,0.67381), rgb=(0.25853,0.28492,0.69300), rgb=(0.26074,0.29568,0.71162), rgb=(0.26280,0.30639,0.72968), rgb=(0.26473,0.31706,0.74718), rgb=(0.26652,0.32768,0.76412), rgb=(0.26816,0.33825,0.78050), rgb=(0.26967,0.34878,0.79631), rgb=(0.27103,0.35926,0.81156), rgb=(0.27226,0.36970,0.82624), rgb=(0.27334,0.38008,0.84037), rgb=(0.27429,0.39043,0.85393), rgb=(0.27509,0.40072,0.86692), rgb=(0.27576,0.41097,0.87936), rgb=(0.27628,0.42118,0.89123), rgb=(0.27667,0.43134,0.90254), rgb=(0.27691,0.44145,0.91328), rgb=(0.27701,0.45152,0.92347), rgb=(0.27698,0.46153,0.93309), rgb=(0.27680,0.47151,0.94214), rgb=(0.27648,0.48144,0.95064), rgb=(0.27603,0.49132,0.95857), rgb=(0.27543,0.50115,0.96594), rgb=(0.27469,0.51094,0.97275), rgb=(0.27381,0.52069,0.97899), rgb=(0.27273,0.53040,0.98461), rgb=(0.27106,0.54015,0.98930), rgb=(0.26878,0.54995,0.99303), rgb=(0.26592,0.55979,0.99583), rgb=(0.26252,0.56967,0.99773), rgb=(0.25862,0.57958,0.99876), rgb=(0.25425,0.58950,0.99896), rgb=(0.24946,0.59943,0.99835), rgb=(0.24427,0.60937,0.99697), rgb=(0.23874,0.61931,0.99485), rgb=(0.23288,0.62923,0.99202), rgb=(0.22676,0.63913,0.98851), rgb=(0.22039,0.64901,0.98436), rgb=(0.21382,0.65886,0.97959), rgb=(0.20708,0.66866,0.97423), rgb=(0.20021,0.67842,0.96833), rgb=(0.19326,0.68812,0.96190), rgb=(0.18625,0.69775,0.95498), rgb=(0.17923,0.70732,0.94761), rgb=(0.17223,0.71680,0.93981), rgb=(0.16529,0.72620,0.93161), rgb=(0.15844,0.73551,0.92305), rgb=(0.15173,0.74472,0.91416), rgb=(0.14519,0.75381,0.90496), rgb=(0.13886,0.76279,0.89550), rgb=(0.13278,0.77165,0.88580), rgb=(0.12698,0.78037,0.87590), rgb=(0.12151,0.78896,0.86581), rgb=(0.11639,0.79740,0.85559), rgb=(0.11167,0.80569,0.84525), rgb=(0.10738,0.81381,0.83484), rgb=(0.10357,0.82177,0.82437), rgb=(0.10026,0.82955,0.81389), rgb=(0.09750,0.83714,0.80342), rgb=(0.09532,0.84455,0.79299), rgb=(0.09377,0.85175,0.78264), rgb=(0.09287,0.85875,0.77240), rgb=(0.09267,0.86554,0.76230), rgb=(0.09320,0.87211,0.75237), rgb=(0.09451,0.87844,0.74265), rgb=(0.09662,0.88454,0.73316), rgb=(0.09958,0.89040,0.72393), rgb=(0.10342,0.89600,0.71500), rgb=(0.10815,0.90142,0.70599), rgb=(0.11374,0.90673,0.69651), rgb=(0.12014,0.91193,0.68660), rgb=(0.12733,0.91701,0.67627), rgb=(0.13526,0.92197,0.66556), rgb=(0.14391,0.92680,0.65448), rgb=(0.15323,0.93151,0.64308), rgb=(0.16319,0.93609,0.63137), rgb=(0.17377,0.94053,0.61938), rgb=(0.18491,0.94484,0.60713), rgb=(0.19659,0.94901,0.59466), rgb=(0.20877,0.95304,0.58199), rgb=(0.22142,0.95692,0.56914), rgb=(0.23449,0.96065,0.55614), rgb=(0.24797,0.96423,0.54303), rgb=(0.26180,0.96765,0.52981), rgb=(0.27597,0.97092,0.51653), rgb=(0.29042,0.97403,0.50321), rgb=(0.30513,0.97697,0.48987), rgb=(0.32006,0.97974,0.47654), rgb=(0.33517,0.98234,0.46325), rgb=(0.35043,0.98477,0.45002), rgb=(0.36581,0.98702,0.43688), rgb=(0.38127,0.98909,0.42386), rgb=(0.39678,0.99098,0.41098), rgb=(0.41229,0.99268,0.39826), rgb=(0.42778,0.99419,0.38575), rgb=(0.44321,0.99551,0.37345), rgb=(0.45854,0.99663,0.36140), rgb=(0.47375,0.99755,0.34963), rgb=(0.48879,0.99828,0.33816), rgb=(0.50362,0.99879,0.32701), rgb=(0.51822,0.99910,0.31622), rgb=(0.53255,0.99919,0.30581), rgb=(0.54658,0.99907,0.29581), rgb=(0.56026,0.99873,0.28623), rgb=(0.57357,0.99817,0.27712), rgb=(0.58646,0.99739,0.26849), rgb=(0.59891,0.99638,0.26038), rgb=(0.61088,0.99514,0.25280), rgb=(0.62233,0.99366,0.24579), rgb=(0.63323,0.99195,0.23937), rgb=(0.64362,0.98999,0.23356), rgb=(0.65394,0.98775,0.22835), rgb=(0.66428,0.98524,0.22370), rgb=(0.67462,0.98246,0.21960), rgb=(0.68494,0.97941,0.21602), rgb=(0.69525,0.97610,0.21294), rgb=(0.70553,0.97255,0.21032), rgb=(0.71577,0.96875,0.20815), rgb=(0.72596,0.96470,0.20640), rgb=(0.73610,0.96043,0.20504), rgb=(0.74617,0.95593,0.20406), rgb=(0.75617,0.95121,0.20343), rgb=(0.76608,0.94627,0.20311), rgb=(0.77591,0.94113,0.20310), rgb=(0.78563,0.93579,0.20336), rgb=(0.79524,0.93025,0.20386), rgb=(0.80473,0.92452,0.20459), rgb=(0.81410,0.91861,0.20552), rgb=(0.82333,0.91253,0.20663), rgb=(0.83241,0.90627,0.20788), rgb=(0.84133,0.89986,0.20926), rgb=(0.85010,0.89328,0.21074), rgb=(0.85868,0.88655,0.21230), rgb=(0.86709,0.87968,0.21391), rgb=(0.87530,0.87267,0.21555), rgb=(0.88331,0.86553,0.21719), rgb=(0.89112,0.85826,0.21880), rgb=(0.89870,0.85087,0.22038), rgb=(0.90605,0.84337,0.22188), rgb=(0.91317,0.83576,0.22328), rgb=(0.92004,0.82806,0.22456), rgb=(0.92666,0.82025,0.22570), rgb=(0.93301,0.81236,0.22667), rgb=(0.93909,0.80439,0.22744), rgb=(0.94489,0.79634,0.22800), rgb=(0.95039,0.78823,0.22831), rgb=(0.95560,0.78005,0.22836), rgb=(0.96049,0.77181,0.22811), rgb=(0.96507,0.76352,0.22754), rgb=(0.96931,0.75519,0.22663), rgb=(0.97323,0.74682,0.22536), rgb=(0.97679,0.73842,0.22369), rgb=(0.98000,0.73000,0.22161), rgb=(0.98289,0.72140,0.21918), rgb=(0.98549,0.71250,0.21650), rgb=(0.98781,0.70330,0.21358), rgb=(0.98986,0.69382,0.21043), rgb=(0.99163,0.68408,0.20706), rgb=(0.99314,0.67408,0.20348), rgb=(0.99438,0.66386,0.19971), rgb=(0.99535,0.65341,0.19577), rgb=(0.99607,0.64277,0.19165), rgb=(0.99654,0.63193,0.18738), rgb=(0.99675,0.62093,0.18297), rgb=(0.99672,0.60977,0.17842), rgb=(0.99644,0.59846,0.17376), rgb=(0.99593,0.58703,0.16899), rgb=(0.99517,0.57549,0.16412), rgb=(0.99419,0.56386,0.15918), rgb=(0.99297,0.55214,0.15417), rgb=(0.99153,0.54036,0.14910), rgb=(0.98987,0.52854,0.14398), rgb=(0.98799,0.51667,0.13883), rgb=(0.98590,0.50479,0.13367), rgb=(0.98360,0.49291,0.12849), rgb=(0.98108,0.48104,0.12332), rgb=(0.97837,0.46920,0.11817), rgb=(0.97545,0.45740,0.11305), rgb=(0.97234,0.44565,0.10797), rgb=(0.96904,0.43399,0.10294), rgb=(0.96555,0.42241,0.09798), rgb=(0.96187,0.41093,0.09310), rgb=(0.95801,0.39958,0.08831), rgb=(0.95398,0.38836,0.08362), rgb=(0.94977,0.37729,0.07905), rgb=(0.94538,0.36638,0.07461), rgb=(0.94084,0.35566,0.07031), rgb=(0.93612,0.34513,0.06616), rgb=(0.93125,0.33482,0.06218), rgb=(0.92623,0.32473,0.05837), rgb=(0.92105,0.31489,0.05475), rgb=(0.91572,0.30530,0.05134), rgb=(0.91024,0.29599,0.04814), rgb=(0.90463,0.28696,0.04516), rgb=(0.89888,0.27824,0.04243), rgb=(0.89298,0.26981,0.03993), rgb=(0.88691,0.26152,0.03753), rgb=(0.88066,0.25334,0.03521), rgb=(0.87422,0.24526,0.03297), rgb=(0.86760,0.23730,0.03082), rgb=(0.86079,0.22945,0.02875), rgb=(0.85380,0.22170,0.02677), rgb=(0.84662,0.21407,0.02487), rgb=(0.83926,0.20654,0.02305), rgb=(0.83172,0.19912,0.02131), rgb=(0.82399,0.19182,0.01966), rgb=(0.81608,0.18462,0.01809), rgb=(0.80799,0.17753,0.01660), rgb=(0.79971,0.17055,0.01520), rgb=(0.79125,0.16368,0.01387), rgb=(0.78260,0.15693,0.01264), rgb=(0.77377,0.15028,0.01148), rgb=(0.76476,0.14374,0.01041), rgb=(0.75556,0.13731,0.00942), rgb=(0.74617,0.13098,0.00851), rgb=(0.73661,0.12477,0.00769), rgb=(0.72686,0.11867,0.00695), rgb=(0.71692,0.11268,0.00629), rgb=(0.70680,0.10680,0.00571), rgb=(0.69650,0.10102,0.00522), rgb=(0.68602,0.09536,0.00481), rgb=(0.67535,0.08980,0.00449), rgb=(0.66449,0.08436,0.00424), rgb=(0.65345,0.07902,0.00408), rgb=(0.64223,0.07380,0.00401), rgb=(0.63082,0.06868,0.00401), rgb=(0.61923,0.06367,0.00410), rgb=(0.60746,0.05878,0.00427), rgb=(0.59550,0.05399,0.00453), rgb=(0.58336,0.04931,0.00486), rgb=(0.57103,0.04474,0.00529), rgb=(0.55852,0.04028,0.00579), rgb=(0.54583,0.03593,0.00638), rgb=(0.53295,0.03169,0.00705), rgb=(0.51989,0.02756,0.00780), rgb=(0.50664,0.02354,0.00863), rgb=(0.49321,0.01963,0.00955), rgb=(0.47960,0.01583,0.01055)}",
            "colorbar"=colormap_settings,
            "colorbar style"=raw"
            {
                height = 2*\pgfkeysvalueof{/pgfplots/parent axis height} + \pgfkeysvalueof{/pgfplots/group/vertical sep},
                ytick={0,0.1,...,1},
                yticklabel style=
                {
                    text width=2.5em,
                    align=right,
                    /pgf/number format/.cd,
                        fixed,
                        fixed zerofill
                }
            }",
            "colorbar/width" = "0.8cm",
            point_meta_min = "0", #"35",
            point_meta_max = "1", #"72",
            title = subtitle,
            view = view_angle,
            xlabel = L"\alpha",
            ylabel = L"\lambda\,\beta^3",
            zlabel = z_label_name,
            zmin = "$(z_min_chart)",
            zmax = "$(z_max_chart)",
            x_dir=x_dir_settings,
            xtick_distance = "0.25",              # set the distance betwen each tick
            ytick_distance = "5",                 # set the distance betwen each tick
            ztick_distance = "$(ztick_distance)", # set the distance betwen each tick
        },
            subplot_p1,
            subplot_p2,
        );
    end
end

fig = @pgf TikzPicture(GroupPlot(
    { group_style = {group_size="4 by 2", vertical_sep="50pt", horizontal_sep="50pt"},
      height="150pt",
      width="225pt"
    },
    gpr[1][2:end]..., gpr[2][2:end]...,)
)

pgfsave("./img/appendix/surface_$(adj_bjk).pdf", fig);
#fig


text:
function affordable_real(
    irreducible_characters,
    multiplicities=fill(1, length(irreducible_characters)),
)
    irr_real = similar(irreducible_characters, 0)
    mls_real = similar(multiplicities, 0)
    for (i, χ) in pairs(irreducible_characters)
        ι = Characters.frobenius_schur(χ)
        if abs(ι) == 1 # real or quaternionic
            @debug "real/quaternionic:" χ
            push!(irr_real, χ)
            push!(mls_real, multiplicities[i])
        else # complex one...
            cχ = conj(χ)
            k = findfirst(==(cχ), irreducible_characters)
            @assert k !== nothing
            @debug "complex" χ conj(χ)=irreducible_characters[k]
            if k > i # ... we haven't already observed a conjugate of
                @assert multiplicities[i] == multiplicities[k]
                push!(irr_real, χ + cχ)
                push!(mls_real, multiplicities[i])
            end
        end
    end

    return irr_real, mls_real
end

"""
    symmetry_adapted_basis([T::Type,] G::AbstractPermutationGroup[, S=Rational{Int};
        semisimple=false])
Compute a basis for the linear space `ℝⁿ` (where `n = degree(G)`) which is
invariant under the symmetry of `G`.

The basis is returned as a vector of `DirectSummand{T}`s (blocks) corresponding
to the irreducible characters of `G`. The blocks are orthogonal to **each other**,
however vectors within a single block may *not* be orthogonal.
If `T<:LinearAlgebra.BlasFloat` BLAS routines will be used to orthogonalize
vectors within each `DirectSummand`.

Arguments:
* `S` controls the types of `Cyclotomic`s used in the computation of
character table. Exact type are preferred. For larger groups `G` `Rational{BigInt}`
might be necessary.
* `T` controls the type of coefficients of the returned basis.
* `semisimple`: if set to `false` (the default) an effort to find minimal
projection system is made, so that each block defines a projection to a single
simple summand within each (isotypical) block. Otherwise an isomorphism to the
semisimple decomposition is computed.

!!! Note:
Each returned block (a `DirectSummand`) is invariant under the action of `G`,
which means that the action may still e.g. permute (row) vectors , but only
*within* each block.

When `semisimple=true` the blocks constitute an isomorphism. Otherwise blocks
may represent only a projection onto the commutant of the appropriate matrix
algebra (which has in general lower dimension). This happens precisely when
`issimple` on a block returns `true` and `SymbolicWedderburn.degree(ds) == 1`.
"""
function symmetry_adapted_basis(
    G::PermutationGroups.AbstractPermutationGroup,
    S::Type = Rational{Int};
    semisimple::Bool = false,
)
    tbl = CharacterTable(S, G)
    return symmetry_adapted_basis(eltype(tbl), tbl, semisimple = semisimple)
end

symmetry_adapted_basis(
    T::Type,
    G::PermutationGroups.AbstractPermutationGroup,
    S::Type = Rational{Int};
    semisimple::Bool = false,
) = symmetry_adapted_basis(T, CharacterTable(S, G), semisimple = semisimple)

function symmetry_adapted_basis(
    T::Type,
    tbl::CharacterTable;
    semisimple::Bool = false,
)
    irr, multips = _constituents_decomposition(
        action_character(conjugacy_classes(tbl), tbl),
        tbl,
    )

    if T <: Real
        irr, multips = affordable_real(irr, multips)
    end

    if semisimple || all(isone ∘ degree, irr)
        return _symmetry_adapted_basis(T, irr, multips)
    else
        RG = _group_algebra(parent(tbl))
        return _symmetry_adapted_basis(T, irr, multips, RG)
    end
end

"""
    symmetry_adapted_basis([T::Type,] G::Group, action, basis[, S=Rational{Int}];
        semisimple=false)
Compute a decomposition of basis into (semi)simple subspaces which are invariant under
the action of `G`.

It is assumed that `G` acts on a subset of basis and the action needs to be
extended to the whole `basis`. If `G` is a permutation group already acting on
the whole `basis`, a call to `symmetry_adapted_basis(G)` is preferred.
* For inducing the action `basis` needs to be indexable and iterable
(e.g. in the form of an `AbstractVector`).
"""
function symmetry_adapted_basis(
    G::Group,
    action::Action,
    basis,
    S::Type = Rational{Int};
    semisimple=false,
)
    tbl = CharacterTable(S, G)
    ehom = CachedExtensionHomomorphism(parent(tbl), action, basis, precompute=true)
    return symmetry_adapted_basis(eltype(tbl), tbl, ehom, semisimple=semisimple)
end

function symmetry_adapted_basis(
    T::Type,
    G::Group,
    action::Action,
    basis,
    S::Type = Rational{Int};
    semisimple=false,
)
    tbl = CharacterTable(S, G)
    ehom = CachedExtensionHomomorphism(parent(tbl), action, basis, precompute=true)
    return symmetry_adapted_basis(T, tbl, ehom, semisimple=semisimple)
end

function symmetry_adapted_basis(
    T::Type,
    tbl::CharacterTable,
    ehom::InducedActionHomomorphism;
    semisimple=false,
)
    ψ = action_character(ehom, tbl)

    irr, multips = _constituents_decomposition(ψ, tbl)
    if T <: Real
        irr, multips = affordable_real(irr, multips)
        @debug "Decomposition into real character spaces:
        degrees:        $(join([lpad(d, 6) for d in degrees], ""))
        multiplicities: $(join([lpad(m, 6) for m in multiplicities], ""))"

    end

    if semisimple || all(isone ∘ degree, irr)
        return _symmetry_adapted_basis(T, irr, multips, ehom)
    else
        RG = _group_algebra(parent(tbl))
        return _symmetry_adapted_basis(T, irr, multips, RG, ehom)
    end
end

function _constituents_decomposition(ψ::Character, tbl::CharacterTable)
    irr = irreducible_characters(tbl)
    degrees = degree.(irr)
    multiplicities = constituents(ψ)

    @debug "Decomposition into character spaces:
    degrees:        $(join([lpad(d, 6) for d in degrees], ""))
    multiplicities: $(join([lpad(m, 6) for m in multiplicities], ""))"

    @assert dot(multiplicities, degrees) == degree(ψ)
    "Something went wrong: characters do not constitute a complete basis for action:
    $(dot(multiplicities, degrees)) ≠ $(degree(ψ))"

    present_irreps = [i for (i, m) in enumerate(multiplicities) if m ≠ 0]
    return irr[present_irreps], multiplicities[present_irreps]
end

function _group_algebra(G::Group)
    @assert isfinite(G)
    b = StarAlgebras.Basis{UInt16}(vec(collect(G)))
    RG = if order(Int, G) <= (typemax(UInt16)>>2)
        StarAlgebra(G, b, (length(b), length(b)), precompute=true)
        # cache is about ~ 1Gb
    else
        StarAlgebra(G, b)
    end
    return RG
end

function _symmetry_adapted_basis(
    T::Type,
    irr::AbstractVector{<:Character},
    multiplicities::AbstractVector{<:Integer},
    hom=nothing
)
    res = map(zip(irr, multiplicities)) do (µ, m)
        @spawn_compat begin
            µT = eltype(µ) == T ? µ : Character{T}(µ)
            image = hom === nothing ? image_basis(µT) : image_basis(hom, µT)
            simple = size(image, 1) == m
            deg = degree(µ)
            @assert size(image, 1) == (simple ? m : m*deg) "incompatible projection dimension: $(size(image, 1)) ≠ $(simple ? m : m*deg)"
            if deg == 1
                @assert simple "Central projection associated to character is not simple unless its degree == 1"
            end
            DirectSummand(image, m, deg, simple)
        end
    end
    return fetch.(res)
end

function _symmetry_adapted_basis(
    T::Type,
    irr::AbstractVector{<:Character},
    multiplicities::AbstractVector{<:Integer},
    RG::StarAlgebra{<:Group},
    hom=nothing,
)
    mps, simples = minimal_projection_system(irr, RG)
    degrees = degree.(irr)
    res = map(zip(mps, multiplicities, degrees, simples)) do (µ, m, deg, simple)
        @spawn_compat begin
            µT = eltype(µ) == T ? µ : AlgebraElement{T}(µ)
            image = hom === nothing ? image_basis(µT) : image_basis(hom, µT)
            @assert size(image, 1) == (simple ? m : m*deg) "incompatible projection dimension: $(size(image, 1)) ≠ $(simple ? m : m*deg)"
            DirectSummand(image, m, deg, simple)
        end
    end
    direct_summands = fetch.(res)

    for (χ, ds) in zip(irr, direct_summands)
        if issimple(ds) && (d = size(ds, 1)) != (e = multiplicity(ds)*sum(constituents(χ).>0))
            throw("The dimension of the projection doesn't match with simple summand multiplicity: $d ≠ $e")
        end
    end

    return direct_summands
end


text:
module Meshes
using DocStringExtensions
export EquispacedRectangleMesh,
    rectangular_mesh,
    cube_panel_mesh,
    sphere_mesh,
    equispaced_rectangular_mesh,
    TensorProductMesh,
    Mesh2D,
    AbstractWarp,
    AbstractSphereWarp,
    EquiangularSphereWarp,
    EquidistantSphereWarp

import ..Domains:
    Domains, IntervalDomain, RectangleDomain, CubePanelDomain, SphereDomain
import IntervalSets: ClosedInterval
import ..Geometry: Geometry

"""
    AbstractMesh

A `Mesh` is an object which represents how we discretize a domain into elements.

It should be lightweight (i.e. exists on all MPI ranks), e.g for meshes stored
in a file, it would contain the filename.
"""
abstract type AbstractMesh{FT} end

Base.eltype(::AbstractMesh{FT}) where {FT} = FT

domain(mesh::AbstractMesh) = getfield(mesh, :domain)
coordinate_type(mesh::AbstractMesh) = Domains.coordinate_type(domain(mesh))

abstract type AbstractWarp end
abstract type AbstractSphereWarp <: AbstractWarp end
struct NoWarp <: AbstractWarp end
struct EquiangularSphereWarp <: AbstractSphereWarp end
struct EquidistantSphereWarp <: AbstractSphereWarp end

struct IntervalMesh{FT, I <: IntervalDomain, V <: AbstractVector, B} <:
       AbstractMesh{FT}
    domain::I
    faces::V
    boundaries::B
end

IntervalMesh{FT}(domain::I, faces::V, boundaries::B) where {FT, I, V, B} =
    IntervalMesh{FT, I, V, B}(domain, faces, boundaries)

abstract type Stretching end
struct Uniform <: Stretching end

function IntervalMesh(
    domain::IntervalDomain{CT},
    ::Uniform;
    nelems,
) where {CT <: Geometry.Abstract1DPoint{FT}} where {FT}
    faces = range(domain.coord_min, domain.coord_max; length = nelems + 1)
    boundaries = NamedTuple{domain.boundary_tags}((5, 6))
    IntervalMesh{FT}(domain, faces, boundaries)
end

# 3.1.2 in the design docs
"""
    ExponentialStretching(H)

Apply exponential stretching to the  domain. `H` is the scale height (a typical atmospheric scale height `H ≈ 7.5e3`km).
"""
struct ExponentialStretching{FT} <: Stretching
    H::FT
end

function IntervalMesh(
    domain::IntervalDomain{CT},
    stretch::ExponentialStretching;
    nelems,
) where {CT <: Geometry.Abstract1DPoint{FT}} where {FT}
    cmin = Geometry.component(domain.coord_min, 1)
    cmax = Geometry.component(domain.coord_max, 1)
    R = cmax - cmin
    h = stretch.H / R
    η(ζ) = -h * log1p(-(1 - exp(-1 / h)) * ζ)
    faces =
        [CT(cmin + R * η(ζ)) for ζ in range(FT(0), FT(1); length = nelems + 1)]
    boundaries = NamedTuple{domain.boundary_tags}((5, 6))
    IntervalMesh{FT, typeof(domain), typeof(faces), typeof(boundaries)}(
        domain,
        faces,
        boundaries,
    )
end

IntervalMesh(domain::IntervalDomain; nelems) =
    IntervalMesh(domain, Uniform(); nelems)

function Base.show(io::IO, mesh::IntervalMesh)
    nelements = length(mesh.faces) - 1
    print(io, nelements, " IntervalMesh of ")
    print(io, mesh.domain)
end

struct EquispacedLineMesh{FT, ID <: IntervalDomain, R} <: AbstractMesh{FT}
    domain::ID
    n1::Int64 # number of elements in x1 direction
    n2::Int64 # always 1
    range1::R
    range2::R # always 1:1
end

function EquispacedLineMesh(domain::IntervalDomain, n1)
    FT = eltype(Domains.coordinate_type(domain))
    cmin = Geometry.component(domain.coord_min, 1)
    cmax = Geometry.component(domain.coord_max, 1)
    range1 = range(cmin, cmax; length = n1 + 1)
    range2 = range(zero(FT), zero(FT), length = zero(n1))
    return EquispacedLineMesh(domain, n1, zero(n1), range1, range2)
end

function Base.show(io::IO, mesh::EquispacedLineMesh)
    print(io, "(", mesh.n1, " × ", " ) EquispacedLineMesh of ")
    print(io, mesh.domain)
end

"""
    EquispacedRectangleMesh(domain::RectangleDomain, n1::Integer, n2::Integer)

A regular `AbstractMesh` of `domain` with `n1` elements in dimension 1, and `n2`
in dimension 2.
"""
struct EquispacedRectangleMesh{FT, RD <: RectangleDomain, R} <: AbstractMesh{FT}
    domain::RD
    n1::Int64 # number of elements in x1 direction
    n2::Int64 # number of elements in x2 direction
    range1::R
    range2::R
end

function EquispacedRectangleMesh(
    domain::RectangleDomain,
    n1::Integer,
    n2::Integer,
)
    FT = eltype(Domains.coordinate_type(domain))
    x1min = Geometry.component(domain.x1x2min, 1)
    x2min = Geometry.component(domain.x1x2min, 2)
    x1max = Geometry.component(domain.x1x2max, 1)
    x2max = Geometry.component(domain.x1x2max, 2)
    range1 = range(x1min, x1max; length = n1 + 1)
    range2 = range(x2min, x2max; length = n2 + 1)
    EquispacedRectangleMesh{FT, typeof(domain), typeof(range1)}(
        domain,
        n1,
        n2,
        range1,
        range2,
    )
end

function Base.show(io::IO, mesh::EquispacedRectangleMesh)
    print(io, mesh.n1, "×", mesh.n2, " EquispacedRectangleMesh of ")
    print(io, mesh.domain)
end

#struct EquiangularCubedSphereMesh{FT} <: AbstractMesh{FT}
#    domain::SphereDomain{FT}
#    n::Int64
#end

"""
    Mesh2D{I,IA2D,FT,FTA2D} <: AbstractMesh{FT}

Conformal mesh for a 2D manifold. The manifold can be
embedded in a higher dimensional space.

                        Quadrilateral

                v3            f4           v4
                  o------------------------o
                  |                        |		  face    vertices
                  |                        |
                  |                        |		   f1 =>  v1 v3
               f1 |                        | f2        f2 =>  v2 v4
                  |                        |		   f3 =>  v1 v2
                  |                        |           f4 =>  v3 v4
                  |                        |
                  |                        |
                  o------------------------o
                 v1           f3           v2

z-order numbering convention for 2D quadtrees

Reference:

p4est: SCALABLE ALGORITHMS FOR PARALLEL ADAPTIVE
MESH REFINEMENT ON FORESTS OF OCTREES∗
CARSTEN BURSTEDDE†, LUCAS C. WILCOX‡ , AND OMAR GHATTAS§
SIAM J. Sci. Comput. Vol. 33, No. 3, pp. 1103-1133

https://p4est.github.io/papers/BursteddeWilcoxGhattas11.pdf

# Fields
$(DocStringExtensions.FIELDS)
"""
struct Mesh2D{D, W, I, IA1D, IA2D, FT, FTA2D, SNT, NB} <: AbstractMesh{FT}
    "domain"
    domain::D
    "warping type"
    warp_type::W
    "# of unique vertices in the mesh"
    nverts::I
    "# of unique faces in the mesh"
    nfaces::I
    "# of elements in the mesh"
    nelems::I
    "# of zones in the mesh"
    nbndry::I
    "x₁, x₂, ... coordinates of vertices `(nverts, dim)`, dim can be greater than 2 for 2D manifolds embedded in higher dimensional space"
    coordinates::FTA2D
    "unique vertices `(n_uniquevertices)`"
    unique_verts::IA1D
    "connectivity information for unique vertices"
    uverts_conn::IA1D
    "offset information for uverts_conn `(n_unique_verts + 1)`"
    uverts_offset::IA1D
    "face vertices numbers `(nfaces, 2)`"
    face_verts::IA2D
    "boundary elems for each face `(nfaces, 5)` -> [elem1, localface1, elem2, localface2, relative orientation]"
    face_neighbors::IA2D
    "face numbers on each boundary `(nfaces)`"
    face_boundary::IA1D
    "boundary tags"
    boundary_tags::IA1D
    "boundary tag names"
    boundary_tag_names::SNT
    "face boundary offset for each boundary for face_boundary array"
    face_boundary_offset::IA1D
    "vertices numbers for each elem `(nelems, 4)`"
    elem_verts::IA2D
    "face numbers for each elem `(nelems, 4)`"
    elem_faces::IA2D
end

Mesh2D(
    domain,
    warp_type,
    nverts,
    nfaces,
    nelems,
    nbndry,
    coordinates,
    unique_verts,
    uverts_conn,
    uverts_offset,
    face_verts,
    face_neighbors,
    face_boundary,
    boundary_tags,
    boundary_tag_names,
    face_boundary_offset,
    elem_verts,
    elem_faces,
) = Mesh2D{
    typeof(domain),
    typeof(warp_type),
    eltype(nverts),
    typeof(face_boundary),
    typeof(face_verts),
    eltype(coordinates),
    typeof(coordinates),
    typeof(boundary_tag_names),
    length(boundary_tag_names),
}(
    domain,
    warp_type,
    nverts,
    nfaces,
    nelems,
    nbndry,
    coordinates,
    unique_verts,
    uverts_conn,
    uverts_offset,
    face_verts,
    face_neighbors,
    face_boundary,
    boundary_tags,
    boundary_tag_names,
    face_boundary_offset,
    elem_verts,
    elem_faces,
)

include("box_mesh.jl")
include("warp_cube_to_sphere.jl")
include("sphere_mesh.jl")

# implementations
include("tensorproductmesh.jl")

end # module


text:
# coding: utf-8

# In[1]:

#shaker sort, or cocktail sort, or cocktail shaker sort
#in this script, we will call it cocktail shaker sort
#with so many names, it is merely a variation of bubble sort
#go the link below to check out the details of bubble sort
# https://github.com/je-suis-tm/search-and-sort/blob/master/bubble,%20selection%20and%20insertion%20sort.py
#bubble sort repeatedly takes comparison from one side to the other
#usually left to right, takes n*(n-1) steps, assuming n is the size
#cocktail shaker sort has a similar mechanism
#except it compares its elements from left to right
#then right to left, later left to right
#the process looks like a cocktail shaker (not really)
#that is how the sort gets its name
#in theory, for a partially sorted list
#cocktail shaker sort takes less time than a bubble sort

# In[2]:

function cocktail_shaker_sort(arr)

    swap=true

    #we use left and right to define each round's start and end point
    left=1
    right=length(arr)

    while (left<right) || !(swap)

        swap=false

        #as usual, we start from left to right
        for i in left+1:right

            if arr[i]<arr[i-1]

                arr[i],arr[i-1]=arr[i-1],arr[i]
                swap=true

            end

        end

        right-=1

        #swap is the key to increase the efficiency
        #once there is no swap happened in the loop
        #we can conclude the target list is fully sorted
        if !(swap)
            return arr
        end

        swap=false

        #then right to left
        for i in right:-1:left+1

            if arr[i]<arr[i-1]

                arr[i],arr[i-1]=arr[i-1],arr[i]
                swap=true

            end

        end

        left+=1

    end

    return arr

end

# In[3]:

for _ in 1:100

    test_arr=rand(1000)

    if !(cocktail_shaker_sort(test_arr)==sort(test_arr))
        printstyled("Erreur",color=:red)
    end

end

max_stars_repo_path: shaker sort.jl

max_stars_repo_name: je-suis-tm/search-and-sort

max_issues_repo_path: shaker sort.jl

max_issues_repo_name: HafizSabia/search-and-sort

max_forks_repo_path: shaker sort.jl

max_forks_repo_name: HafizSabia/search-and-sort


text:
using EcoNet
using Test

@testset "EcoNet.jl" begin
    # Write your tests here.
end


text:
"""
Abstract type representing an amount of data, i.e. a certain number of bytes,
with a unit prefix (also "metric prefix"). Examples include the SI prefixes,
like KB, MB, and GB, but also the binary prefixes (ISO/IEC 80000), like KiB, MiB, and GiB.

See https://en.wikipedia.org/wiki/Binary_prefix for more information.
"""
abstract type UnitPrefixedBytes <: Number end

"Bytes"
struct B <: UnitPrefixedBytes
    value::Float64
end
"Kibibytes, i.e. `2^10 = 1024` bytes"
struct KiB <: UnitPrefixedBytes
    value::Float64
end
"Mebibytes, i.e. `2^20 = 1024^2` bytes"
struct MiB <: UnitPrefixedBytes
    value::Float64
end
"Gibibytes, i.e. `2^30 = 1024^3` bytes"
struct GiB <: UnitPrefixedBytes
    value::Float64
end
"Tebibytes, i.e. `2^40 = 1024^4` bytes"
struct TiB <: UnitPrefixedBytes
    value::Float64
end
"Kilobytes, i.e. `10^3 = 1000` bytes"
struct KB <: UnitPrefixedBytes
    value::Float64
end
"Megabytes, i.e. `10^6 = 1000^2` bytes"
struct MB <: UnitPrefixedBytes
    value::Float64
end
"Gigabytes, i.e. `10^9 = 1000^3` bytes"
struct GB <: UnitPrefixedBytes
    value::Float64
end
"Terabytes, i.e. `10^12 = 1000^4` bytes"
struct TB <: UnitPrefixedBytes
    value::Float64
end

# accessor
value(x::UnitPrefixedBytes) = x.value

# pretty printing
function Base.show(io::IO, x::T) where {T<:UnitPrefixedBytes}
    v = value(x)
    v_rounded = round(v; digits=2)
    if abs(v - v_rounded) > 0
        print(io, "~", v_rounded, " ", T)
    else
        print(io, v, " ", T)
    end
end

# zero
Base.zero(::Type{T}) where {T<:UnitPrefixedBytes} = T(0.0)
Base.zero(::T) where {T<:UnitPrefixedBytes} = T(0.0)

# one
Base.one(::Type{T}) where {T<:UnitPrefixedBytes} = T(1.0)
Base.one(::T) where {T<:UnitPrefixedBytes} = T(1.0)

# conversion
Base.convert(::Type{B}, x::B) = B(value(x)) # to avoid infinite recursion / stackoverflow
Base.convert(::Type{B}, x::KB) = B(value(x) * 1e3)
Base.convert(::Type{B}, x::MB) = B(value(x) * 1e6)
Base.convert(::Type{B}, x::GB) = B(value(x) * 1e9)
Base.convert(::Type{B}, x::TB) = B(value(x) * 1e12)
Base.convert(::Type{B}, x::KiB) = B(value(x) * 2^10)
Base.convert(::Type{B}, x::MiB) = B(value(x) * 2^20)
Base.convert(::Type{B}, x::GiB) = B(value(x) * 2^30)
Base.convert(::Type{B}, x::TiB) = B(value(x) * 2^40)

Base.convert(::Type{KB}, x::B) = KB(value(x) * 1e-3)
Base.convert(::Type{MB}, x::B) = MB(value(x) * 1e-6)
Base.convert(::Type{GB}, x::B) = GB(value(x) * 1e-9)
Base.convert(::Type{TB}, x::B) = TB(value(x) * 1e-12)
Base.convert(::Type{KiB}, x::B) = KiB(value(x) * 2^(-10))
Base.convert(::Type{MiB}, x::B) = MiB(value(x) * 2^(-20))
Base.convert(::Type{GiB}, x::B) = GiB(value(x) * 2^(-30))
Base.convert(::Type{TiB}, x::B) = TiB(value(x) * 2^(-40))

function Base.convert(::Type{T}, x::UnitPrefixedBytes) where {T<:UnitPrefixedBytes}
    return convert(T, convert(B, x))
end

function Base.isapprox(x::UnitPrefixedBytes, y::UnitPrefixedBytes; kwargs...)
    return isapprox(bytes(x), bytes(y); kwargs...)
end
# Base.:(==)(x::UnitPrefixedBytes, y::UnitPrefixedBytes) = bytes(x) == bytes(y) # unnecessary
# Base.:(===)(x::UnitPrefixedBytes, y::UnitPrefixedBytes) = (typeof(x) == typeof(y)) && (bytes(x) == bytes(y)) # unnecessary

# promotion
Base.promote_rule(::Type{B}, ::Type{KB}) = KB
Base.promote_rule(::Type{B}, ::Type{MB}) = MB
Base.promote_rule(::Type{B}, ::Type{GB}) = GB
Base.promote_rule(::Type{B}, ::Type{TB}) = TB
Base.promote_rule(::Type{B}, ::Type{KiB}) = KiB
Base.promote_rule(::Type{B}, ::Type{MiB}) = MiB
Base.promote_rule(::Type{B}, ::Type{GiB}) = GiB
Base.promote_rule(::Type{B}, ::Type{TiB}) = TiB

Base.promote_rule(::Type{KB}, ::Type{MB}) = MB
Base.promote_rule(::Type{KB}, ::Type{GB}) = GB
Base.promote_rule(::Type{KB}, ::Type{TB}) = TB
Base.promote_rule(::Type{KB}, ::Type{KiB}) = KiB
Base.promote_rule(::Type{KB}, ::Type{MiB}) = MiB
Base.promote_rule(::Type{KB}, ::Type{GiB}) = GiB
Base.promote_rule(::Type{KB}, ::Type{TiB}) = TiB

Base.promote_rule(::Type{MB}, ::Type{GB}) = GB
Base.promote_rule(::Type{MB}, ::Type{TB}) = TB
Base.promote_rule(::Type{MB}, ::Type{KiB}) = MB
Base.promote_rule(::Type{MB}, ::Type{MiB}) = MiB
Base.promote_rule(::Type{MB}, ::Type{GiB}) = GiB
Base.promote_rule(::Type{MB}, ::Type{TiB}) = TiB

Base.promote_rule(::Type{GB}, ::Type{TB}) = TB
Base.promote_rule(::Type{GB}, ::Type{KiB}) = GB
Base.promote_rule(::Type{GB}, ::Type{MiB}) = GB
Base.promote_rule(::Type{GB}, ::Type{GiB}) = GiB
Base.promote_rule(::Type{GB}, ::Type{TiB}) = TiB

Base.promote_rule(::Type{TB}, ::Type{KiB}) = TB
Base.promote_rule(::Type{TB}, ::Type{MiB}) = TB
Base.promote_rule(::Type{TB}, ::Type{GiB}) = TB
Base.promote_rule(::Type{TB}, ::Type{TiB}) = TiB

Base.promote_rule(::Type{KiB}, ::Type{MiB}) = MiB
Base.promote_rule(::Type{KiB}, ::Type{GiB}) = GiB
Base.promote_rule(::Type{KiB}, ::Type{TiB}) = TiB

Base.promote_rule(::Type{MiB}, ::Type{GiB}) = GiB
Base.promote_rule(::Type{MiB}, ::Type{TiB}) = TiB

Base.promote_rule(::Type{GiB}, ::Type{TiB}) = TiB

# basic arithmetics
Base.:+(x::T, y::T) where {T<:UnitPrefixedBytes} = T(value(x) + value(y))
Base.:-(x::T, y::T) where {T<:UnitPrefixedBytes} = T(value(x) - value(y))

Base.:*(x::Number, y::T) where {T<:UnitPrefixedBytes} = T(x * value(y))
Base.:*(x::T, y::Number) where {T<:UnitPrefixedBytes} = y * x # symmetry / commutativity

Base.:/(x::T, y::Number) where {T<:UnitPrefixedBytes} = T(value(x) / y)

# utility

"""
    bytes(x::UnitPrefixedBytes)
Return the number of bytes (without prefix) as `Float64`.
"""
bytes(x::UnitPrefixedBytes) = value(convert(B, x))

"""
    bytes(x::Number)
Returns an appropriate `UnitPrefixedBytes` object, representing the number of bytes.

**Note:** This function is type unstable by construction!

See [`simplify`](@ref) for what "appropriate" means here.
"""
bytes(x::Number) = simplify(B(x))

"""
Toggle between
* Base 10, SI prefixes, i.e. factors of 1000
* Base 2, ISO/IEC prefixes, i.e. factors of 1024

**Example:**
```julia
julia> change_base(KB(13))
~12.7 KiB

julia> change_base(KiB(13))
~13.31 KB
```
"""
change_base(x::T) where {T<:UnitPrefixedBytes} = convert(_base_partner(T), x)
_base_partner(::Type{B}) = B
_base_partner(::Type{KB}) = KiB
_base_partner(::Type{MB}) = MiB
_base_partner(::Type{GB}) = GiB
_base_partner(::Type{TB}) = TiB
_base_partner(::Type{KiB}) = KB
_base_partner(::Type{MiB}) = MB
_base_partner(::Type{GiB}) = GB
_base_partner(::Type{TiB}) = TB

"""
    simplify(x::UnitPrefixedBytes[; base])
Given a `UnitPrefixedBytes` number `x`, finds a more appropriate
`UnitPrefixedBytes` that represents the same number of bytes
but with a smaller value.

The optional keyword argument `base` can be used to switch between
base 2, i.e. ISO/IEC prefixes (default), and base 10. Allowed values are
`2`, `10`, `:SI`, `:ISO`, and `:IEC`.

**Note:** This function is type unstable by construction!

**Example:**
```julia
julia> simplify(B(40_000_000))
~38.15 MiB

julia> simplify(B(40_000_000); base=10)
40.0 MB
```
"""
function simplify(x::B; base=2)
    if base == 2 || base == :SI
        hierarchy = (B, KiB, MiB, GiB, TiB)
        factor = 1024
    elseif base == 10 || base == :IEC || base == :ISO
        hierarchy = (B, KB, MB, GB, TB)
        factor = 1000
    else
        throw(ArgumentError("Only base 2 or base 10 supported."))
    end
    v = value(x)
    for T in hierarchy
        v /= factor
        v < 1 && return convert(T, x)
    end
    return convert(last(hierarchy), x)
end
simplify(x::UnitPrefixedBytes; kwargs...) = simplify(convert(B, x); kwargs...)


text:
"""
Tests for the function `estimate_h2`.
Tests all fail because `ld_score_regression` is only half-implemented.
"""

include("../src/LDScore.jl")

using Main
using Main.LDScore
using Test

approx(x, y, eps) = abs(x - y) <= eps

function vector_approx(x, y, eps)
    if size(x) != size(y)
        return false
    end
    for (i, j) in enumerate(x)
        if !approx(x[i], y[i], eps) return false end
    end
    return true
end

LDScore.args["ref-ld"] = "test/test_ldscore/oneld_onefile"
LDScore.args["w-ld"] = "test/test_ldscore/w"
LDScore.args["h²"] = "test/test_sumstats/1"
LDScore.args["out"] = "test/test_out/1"

# Test case from original ldsc repo
x = LDScore.estimate_h2()

# TODO finish implementing `ld_score_regression`


text:
export Fourier,Taylor,Hardy,CosSpace,SinSpace,Laurent

for T in (:CosSpace,:SinSpace)
    @eval begin
        doc"""
            `CosSpace` is the basis `[1,cos θ,cos 2θ,...]`
            `SinSpace` is the basis `[sin θ,sin 2θ,...]`
        """
        immutable $T{D<:Domain} <: RealUnivariateSpace{D}
            domain::D
        end
        $T()=$T(PeriodicInterval())
        spacescompatible(a::$T,b::$T)=domainscompatible(a,b)
        hasfasttransform(::$T)=true
        canonicalspace(S::$T)=Fourier(domain(S))
        setdomain(S::$T,d::Domain)=$T(d)
    end
end
# s == true means analytic inside, taylor series
# s == false means anlytic outside and decaying at infinity

doc"""
    `Hardy{true}` is the basis `[1,z,z^2,...]`
    `Hardy{false}` is the basis `[1/z,1/z^2,...]`
"""
immutable Hardy{s,D<:Domain} <: UnivariateSpace{ComplexBasis,D}
    domain::D
    Hardy(d)=new(d)
    Hardy()=new(D())
end

Base.promote_rule{T<:Number,S<:Union{Hardy{true},CosSpace},V}(::Type{Fun{S,V}},::Type{T}) =
    Fun{S,promote_type(V,T)}
Base.promote_rule{T<:Number,S<:Union{Hardy{true},CosSpace}}(::Type{Fun{S}},::Type{T}) =
    Fun{S,T}

@compat (H::Type{Hardy{s}}){s}(d::Domain) = Hardy{s,typeof(d)}(d)
@compat (H::Type{Hardy{s}}){s}() = Hardy{s}(Circle())

canonicalspace(S::Hardy) = S
setdomain{s}(S::Hardy{s},d::Domain) = Hardy{s}(d)

spacescompatible{s}(a::Hardy{s},b::Hardy{s}) = domainscompatible(a,b)
hasfasttransform(::Hardy) = true

# The <: Domain is crucial for matching Basecall overrides
typealias Taylor{D<:Domain} Hardy{true,D}

plan_transform(::Taylor,x::Vector)=plan_fft(x)
plan_itransform(::Taylor,x::Vector)=plan_ifft(x)
transform(::Taylor,vals::Vector,plan)=alternatesign!(plan*vals/length(vals))
itransform(::Taylor,cfs::Vector,plan)=plan*alternatesign!(cfs)*length(cfs)

plan_transform(::Hardy{false},x::Vector)=plan_fft(x)
plan_itransform(::Hardy{false},x::Vector)=plan_ifft(x)
transform(::Hardy{false},vals::Vector,plan)=-alternatesign!(flipdim(plan*vals,1)/length(vals))
itransform(::Hardy{false},cfs::Vector,plan)=plan*flipdim(alternatesign!(-cfs),1)*length(cfs)

evaluate{D<:Domain}(f::AbstractVector,S::Taylor{D},z) = horner(f,fromcanonical(Circle(),tocanonical(S,z)))
function evaluate{D<:Circle}(f::AbstractVector,S::Taylor{D},z)
    z=mappoint(S,𝕌,z)
    d=domain(S)
    horner(f,z)
end

function evaluate{D<:Domain}(f::AbstractVector,S::Hardy{false,D},z)
    z=mappoint(S,𝕌,z)
    z=1./z
    z.*horner(f,z)
end
function evaluate{D<:Circle}(f::AbstractVector,S::Hardy{false,D},z)
    z=mappoint(S,𝕌,z)
    z=1./z
    z.*horner(f,z)
end

##TODO: fast routine

function horner(c::AbstractVector,kr::Range{Int64},x)
    T = promote_type(eltype(c),eltype(x))
    if isempty(c)
        return zero(x)
    end

    ret = zero(T)
    @inbounds for k in reverse(kr)
        ret = muladd(x,ret,c[k])
    end

    ret
end

function horner(c::AbstractVector,kr::Range{Int64},x::AbstractVector)
    n,T = length(x),promote_type(eltype(c),eltype(x))
    if isempty(c)
        return zero(x)
    end

    ret = zeros(T,n)
    @inbounds for k in reverse(kr)
        ck = c[k]
        @simd for i = 1:n
            ret[i] = muladd(x[i],ret[i],ck)
        end
    end

    ret
end

horner(c::AbstractVector,x) = horner(c,1:length(c),x)
horner(c::AbstractVector,x::AbstractArray) = horner(c,1:length(c),x)
horner(c::AbstractVector,kr::Range{Int64},x::AbstractArray) = reshape(horner(c,kr,vec(x)),size(x))

## Cos and Sin space

points(sp::CosSpace,n)=points(domain(sp),2n-2)[1:n]
plan_transform(::CosSpace,x::Vector)=plan_chebyshevtransform(x;kind=2)
plan_itransform(::CosSpace,x::Vector)=plan_ichebyshevtransform(x;kind=2)
transform(::CosSpace,vals,plan)=chebyshevtransform(vals,plan;kind=2)
itransform(::CosSpace,cfs,plan)=ichebyshevtransform(cfs,plan;kind=2)
evaluate(f::Vector,S::CosSpace,t)=clenshaw(Chebyshev(),f,cos(tocanonical(S,t)))

points(sp::SinSpace,n)=points(domain(sp),2n+2)[n+3:2n+2]
plan_transform{T<:FFTW.fftwNumber}(::SinSpace,x::Vector{T})=FFTW.plan_r2r(x,FFTW.RODFT00)
plan_itransform{T<:FFTW.fftwNumber}(::SinSpace,x::Vector{T})=FFTW.plan_r2r(x,FFTW.RODFT00)

plan_transform{D}(::SinSpace{D},x::Vector)=error("transform for Fourier only implemented for fftwNumbers")
plan_itransform{D}(::SinSpace{D},x::Vector)=error("transform for Fourier only implemented for fftwNumbers")

transform(::SinSpace,vals,plan)=plan*vals/(length(vals)+1)
itransform(::SinSpace,cfs,plan)=plan*cfs/2
evaluate(f::AbstractVector,S::SinSpace,t)=sineshaw(f,tocanonical(S,t))

## Laurent space

typealias Laurent{DD} SumSpace{Tuple{Hardy{true,DD},Hardy{false,DD}},ComplexBasis,DD,1}

plan_transform{DD}(::Laurent{DD},x::Vector)=plan_svfft(x)
plan_itransform{DD}(::Laurent{DD},x::Vector)=plan_isvfft(x)
transform{DD}(::Laurent{DD},vals,plan)=svfft(vals,plan)
itransform{DD}(::Laurent{DD},cfs,plan)=isvfft(cfs,plan)

function evaluate{DD}(f::AbstractVector,S::Laurent{DD},z)
    z = mappoint(domain(S),Circle(),z)
    invz = 1./z
    horner(f,1:2:length(f),z) + horner(f,2:2:length(f),invz).*invz
end

function Base.conj{DD}(f::Fun{Laurent{DD}})
    cfs=Array(eltype(f),iseven(ncoefficients(f))?ncoefficients(f)+1:ncoefficients(f))
    cfs[1]=conj(f.coefficients[1])
    for k=2:2:ncoefficients(f)-1
        cfs[k]=conj(f.coefficients[k+1])
    end
    for k=3:2:ncoefficients(f)
        cfs[k]=conj(f.coefficients[k-1])
    end
    Fun(cfs,space(f))
end

## Fourier space

typealias Fourier{DD} SumSpace{Tuple{CosSpace{DD},SinSpace{DD}},RealBasis,DD,1}

for TYP in (:Laurent,:Fourier)
    @eval begin
        @compat (::Type{$TYP})(d::Domain) = $TYP{typeof(d)}(d)
        @compat (::Type{$TYP})() = $TYP(PeriodicInterval())
        @compat (::Type{$TYP}){T<:Number}(d::Vector{T}) = $TYP(PeriodicDomain(d))

        hasfasttransform{D}(::$TYP{D}) = true
    end
end

for T in (:CosSpace,:SinSpace)
    @eval begin
        # override default as canonicalspace must be implemented
        maxspace{D}(::$T,::Fourier{D}) = NoSpace()
        maxspace{D}(::Fourier{D},::$T) = NoSpace()
    end
end

points{D}(sp::Fourier{D},n)=points(domain(sp),n)
plan_transform{T<:FFTW.fftwNumber,D}(::Fourier{D},x::Vector{T}) = FFTW.plan_r2r(x, FFTW.R2HC)
plan_itransform{T<:FFTW.fftwNumber,D}(::Fourier{D},x::Vector{T}) = FFTW.plan_r2r(x, FFTW.HC2R)

plan_transform{D}(::Fourier{D},x::Vector)=error("transform for Fourier only implemented for fftwNumbers")
plan_itransform{D}(::Fourier{D},x::Vector)=error("transform for Fourier only implemented for fftwNumbers")

function transform{T<:Number,D}(::Fourier{D},vals::Vector{T},plan)
    n=length(vals)
    cfs=2plan*vals/n
    cfs[1]/=2
    if iseven(n)
        cfs[div(n,2)+1]/=2
    end

    fouriermodalt!(cfs)

    ret=Array(T,n)
    if iseven(n)
        ret[1:2:end]=cfs[1:div(n,2)]
        ret[2:2:end]=cfs[end:-1:div(n,2)+1]
    else
        ret[1:2:end]=cfs[1:div(n+1,2)]
        ret[2:2:end]=cfs[end:-1:div(n+3,2)]
    end
    ret
end

function itransform{T<:Number,D}(::Fourier{D},a::Vector{T},plan)
    n=length(a)
    cfs=[a[1:2:end];
            flipdim(a[2:2:end],1)]
    fouriermodalt!(cfs)
    if iseven(n)
        cfs[div(n,2)+1]*=2
    end
    cfs[1]*=2
    plan*cfs/2
end

function fouriermodalt!(cfs)
    n=length(cfs)
    if iseven(n)
        for k=2:2:div(n,2)+1
            cfs[k]*=-1
        end
    else
        for k=2:2:div(n+1,2)
            cfs[k]*=-1
        end
    end

    if mod(n,4)==0
        for k=div(n,2)+3:2:n
            cfs[k]*=-1
        end
    elseif mod(n,4)==2
        for k=div(n,2)+2:2:n
            cfs[k]*=-1
        end
    elseif mod(n,4)==1
        for k=div(n+3,2):2:n
            cfs[k]*=-1
        end
    else #mod(n,4)==3
        for k=div(n+5,2):2:n
            cfs[k]*=-1
        end
    end
    cfs
end

canonicalspace{DD<:PeriodicInterval}(S::Laurent{DD})=Fourier(domain(S))
canonicalspace{DD<:Circle}(S::Fourier{DD})=Laurent(domain(S))
canonicalspace{DD<:PeriodicLine}(S::Laurent{DD})=S

for TYP in (:CosSpace,:Taylor)
    @eval union_rule(A::ConstantSpace,B::$TYP)=B
end

## Ones and zeros

for sp in (:Fourier,:CosSpace,:Laurent,:Taylor)
    @eval begin
        Base.ones{T<:Number,D}(::Type{T},S::$sp{D})=Fun(ones(T,1),S)
        Base.ones{D}(S::$sp{D})=Fun(ones(1),S)
    end
end

function identity_fun{DD<:Circle}(S::Taylor{DD})
    d=domain(S)
    if d.orientation
        Fun([d.center,d.radius],S)
    else
        error("Cannot create identity on $S")
    end
end

reverseorientation{D}(f::Fun{Fourier{D}})=Fun(alternatesign!(copy(f.coefficients)),Fourier(reverse(domain(f))))
function reverseorientation{D}(f::Fun{Laurent{D}})
    # exp(im*k*x) -> exp(-im*k*x), or equivalentaly z -> 1/z
    n=ncoefficients(f)
    ret=Array(eltype(f),iseven(n)?n+1:n)  # since z -> 1/z we get one more coefficient
    ret[1]=f.coefficients[1]
    for k=2:2:length(ret)-1
        ret[k+1]=f.coefficients[k]
    end
    for k=2:2:n-1
        ret[k]=f.coefficients[k+1]
    end
    iseven(n) && (ret[n] = 0)

    Fun(ret,Laurent(reverse(domain(f))))
end

include("calculus.jl")
include("specialfunctions.jl")
include("FourierOperators.jl")
include("LaurentOperators.jl")
include("LaurentDirichlet.jl")


text:
############################################
#
#    dolfin.cpp.function
# https://fenicsproject.org/docs/dolfin/latest/python/_autogenerated/dolfin.cpp.function.html
#
#    dolfin.function
# https://fenicsproject.org/docs/dolfin/2017.2.0/python/programmers-reference/functions/index.html
#
############################################

@pyclass dolfin Expression
@pyclass ufl Coefficient Expression
@pyclass dolfin Constant Expression
@pyclass dolfin MeshCoordinates Expression
@pyclass dolfin FacetArea Expression

@pyclass dolfin FunctionSpace
@pyclass dolfin Function FeObject FeFunction
@pyclass dolfin MultiMeshFunctionSpace
@pyclass dolfin MultiMeshFunction

@pyfunc dolfin interpolate
@pyfunc dolfin TensorFunctionSpace
@pyfunc dolfin VectorFunctionSpace
@pyfunc dolfin TrialFunction
@pyfunc dolfin TrialFunctions
@pyfunc dolfin TestFunction
@pyfunc dolfin TestFunctions

OpType = Union{Expression, Variable, FeFunction}

export Expression, Coefficient, Constant, FeFunction, FunctionSpace, MultiMeshFunction,
       FacetArea, MeshCoordinates, interpolate, VectorFunctionSpace, MultiMeshFunctionSpace,
       TensorFunctionSpace, TrialFunction, TrialFunctions, TestFunction, TestFunctions

############################################
#
#    ufi Expression
#
# https://fenicsproject.org/docs/ufl/1.6.0/ufl.html
#
############################################

# no need to export
# not exported in `dolfin` or `ufl`
@pyclass dolfin Equation Expression

# export
@pyclass ufl Argument Expression
@pyclass dolfin Form Expression
@pyclass ufl Measure
@pyclass ufl Identity Expression

array(vec::GenericVector) = vec.pyobject.gather_on_zero()
array(vec::GenericMatrix) = vec.pyobject.array()
array(form::Expression) =  array(assemble(form)).gather_on_zero()
function array(solution::FeFunction)
    generic_vector = solution.pyobject.vector()
    instantiated_vector = dolfin.Vector(generic_vector)
    instantiated_vector.gather_on_zero()
end

len(U::OpType) = length(U.pyobject)

export Argument, Form, Measure, Identity, array, len


text:
using Test, FinancialMonteCarlo;
@show "Black Scholes Model"
S0=100.0;
K=100.0;
r=0.02;
T=1.0;
d=0.01;
D=90.0;

underlying_name="ENJ"
underlying_name_2="ABPL"

underlying_=underlying_name*"_"*underlying_name_2*"_TESL"

Nsim=10000;
Nstep=30;
sigma=0.2;
mc=MonteCarloConfiguration(Nsim,Nstep);
mc1=MonteCarloConfiguration(Nsim,Nstep,FinancialMonteCarlo.AntitheticMC());
toll=0.8

rfCurve=ZeroRate(r);

FwdData=Forward(T)
EUData=EuropeanOptionND(T,K)
AMData=AmericanOption(T,K)
BarrierData=BarrierOptionDownOut(T,K,D)
AsianFloatingStrikeData=AsianFloatingStrikeOption(T)
AsianFixedStrikeData=AsianFixedStrikeOption(T,K)
Model_enj=BlackScholesProcess(sigma,Underlying(S0,d));
Model_abpl=BlackScholesProcess(sigma,Underlying(S0,d));
Model_tesl=BlackScholesProcess(sigma,Underlying(S0,d));
rho=[1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0];
Model=GaussianCopulaNVariateProcess(rho,Model_enj,Model_abpl,Model_tesl)

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

p=EuropeanOption(1.0,70.0)+Forward(1.0)
p2=3.0*p;
port_=Array{typeof(p)}(undef,7)
port_[1]=p
port_[2]=p2

pricer(Model,rfCurve,mc,port_)