text:
# test for unify.jl
using Test
#include("utils.jl") # for isinvar()
#include("unify.jl")

#==
idea:
unify([:x,:y], "P(x,y)", "P(x,f(x))")
unify([:x,:y], :(P(x,y)), :(P(x,f(x))))

:(P(x,y)) == Meta.Meta.parse("P(x,y)")
youcan readline(from console) and Meta.parse it to get the Expr
==#

@testset "unify0 find a pair unifiable" begin
 @test unify0([], :P, :P) ==()
 @test unify0([:x], :x, :x) ==()
 @test unify0([:x], :x, :a) ==(:x, :a)
 @test unify0([:x], :a, :x) ==(:x, :a)

 @test_throws ICMP unify0([:x], :P, :Q)
 @test_throws ICMP unify0([], :P, :Q)

 @test unify0([:x], :(P(x)), :(P(a))) == (:x, :a)
 @test unify0([:x], :(P(x)), :(P(f(a)))) == (:x,:(f(a)))
 @test unify0([:x], :(P(f(a))), :(P(f(x)))) == (:x,:a)
 @test unify0([:x], :(P(f(a))), :(P(x))) == (:x,:(f(a)))

 @test_throws ICMP unify0([:x],:(P(g(a))),:(P(f(x))))
 @test_throws ICMP unify0([:x],:(P(f(a))),:(P(f(b))))
 @test_throws ICMP unify0([:x],:(P(a)), :(P(f(b))))
end

@testset "unify1" begin
 @test unify1([], :(f()), :(f()), []) == []
 @test unify1([:x], :x, :x, [:x]) == [:x]
 @test unify1([:x], :x, :a, [:x]) == [:a]
 @test unify1([:x], :x, :x, [:x]) == [:x]

 @test unify1([:x], :(f(x)), :(f(x)), [:x]) == [:x]
 @test unify1([:x], :(f(x)), :(f(a)), [:x]) == [:a]
 @test unify1([:x], :(f(x)), :(f(x)), [:b]) == [:b]

 @test_throws ICMP unify1([:x], :(f(x)), :(f(a)), [:b])

 @test unify1([:x,:y], :x , :(f(a)), [:x,:y]) == [:(f(a)), :y]
 @test unify1([:x,:y], :x , :(f(a)), [:x,:a]) == [:(f(a)), :a]

 @test unify1([:x,:y,:z], :a , :z, [:y,:z,:z]) == [:a,:a,:a]
 @test unify1([:x,:y,:z], :z , :a, [:y,:z,:z]) == [:a,:a,:a]

 @test unify1([:x,:y,:z], :z , :(f(a)), [:y,:z,:z]) == [:(f(a)),:(f(a)),:(f(a))]

 @test unify1([:x,:y], :x , :(f(a)), [:x,:x]) == [:(f(a)), :(f(a))]

 @test unify1([:x,:y,:z],:x,:(f(a)),[:y,:z,:z]) == [:(f(a)),:(f(a)),:(f(a))]
 @test unify1([:x,:y,:z],:z,:(f(a)),[:y,:z,:z]) == [:(f(a)),:(f(a)),:(f(a))]
 @test unify1([:x,:y,:z,:w,:u],:z,:(f(a)),[:y,:z,:w,:u,:u]) == [:(f(a)),:(f(a)),:(f(a)),:(f(a)),:(f(a))]

 @test unify1([:x,:y], :x , :(f(a)), [:y,:y]) == [:(f(a)), :(f(a))]

 @test unify1([:x,:y], :(f(x)), :(f(x)), [:x,:a]) == [:x,:a]
 @test unify1([:x,:y], :(f(x)), :(f(b)), [:x,:a]) == [:b,:a]

 @test unify1([:x,:y], :(f(b)), :(f(x)), [:x,:a]) == [:b,:a]
 @test unify1([:x,:y], :(f(g(b))), :(f(x)), [:x,:a]) == [:(g(b)),:a]

 @test unify1([:x,:y], :(f(x)), :(f(y)), [:b,:y]) == [:b,:b]

 @test unify1([:x,:y], :(f(g(b,y))), :(f(x)), [:x,:(h(a))]) == [:(g(b,y)),:(h(a))]
 @test unify1([:x,:y], :x, :p, [:y,:y]) == [:p,:p]

 @test unify1([:x,:y], :x, :y, [:p,:y]) == [:p,:p]

end

@testset "unify" begin
 @test_throws ICMP unify([],:(P()),:(Q()))

 @test unify([:x],:(),:())== [:x]
 @test unify([:x],:(P()),:(P())) == [:x]
 @test_throws ICMP unify([:x],:(P()),:(Q()))

 @test unify([:x],:(P(x)), :(P(x))) == [:x]
 @test unify([:y],:(P(x)), :(P(x))) == [:y]

 @test unify([:x],:(P(x)), :(P(y))) == [:y]
 @test unify([:y],:(P(x)), :(P(y))) == [:x]

 @testset  "with number" begin
  @test unify([], :(f(2)), :(f(2))) == []
  @test unify([:x], :(f(x)), :(f(2))) == [2]
  @test unify([:x], :(f(3)), :(f(x))) == [3]
  @test unify([:x,:y], :(f(3,y)), :(f(x,23))) == [3,23]
 end
# in this case, x and y are constants
 @test_throws ICMP unify([:z],:(P(x)), :(P(y)))

 @test unify([],:(P(a)), :(P(a))) == []
 @test_throws ICMP unify([],:(P(a)), :(P(b)))

 (σ=unify([:x,:y],:(P(x)),:(P(y))); @test σ==[:x,:x] || σ==[:y,:y])

 @test unify([:x,:y],:(P(x)),:(P(a))) == [:a,:y]
 @test unify([:x,:y],:(P(f(x))),:(P(f(a)))) == [:a,:y]

 @test unify([:x,:y],:(P(x)),:(P(f(a)))) == [:(f(a)),:y]
 (σ = unify([:x,:y],:(P(x,f(a,x))),:(P(y,f(a,y)))); @test σ == [:y,:y] || σ == [:y,:y])

 @test unify([:x,:y],:(P(x,f(x))),:(P(a,f(y)))) == [:a,:a]
 @test unify([:x,:y],:(P(x,y)),:(P(y,f(a)))) == [:(f(a)),:(f(a))]

 @test unify([:x,:y], :(f(g(b,y),h(a))), :(f(x,y))) == [:(g(b,h(a))),:(h(a))]

 @testset "fixed point of subst" begin
   @test fp_subst([:x,:y],[:x,:y]) == [:x,:y]
   @test fp_subst([:x,:y],[:a,:x]) == [:a,:a]

   @test fp_subst([:x,:y],[:x,:x]) == [:x,:x]

   @test fp_subst([:x,:y,:z],[:a,:(f(x,z)),:b]) == [:a,:(f(a,b)),:b]
   @test fp_subst([:x,:y,:z,:w],[:x,:(f(g(z,w),h(w))),:(k(w)),:(m(a))]) == [:x,:(f(g(k(m(a)),m(a)),h(m(a)))),:(k(m(a))),:(m(a))]

   @test fp_subst([:x,:y,:z,:u],[:(f(y)),:(g(z)),:(h(u)),:(k(a))]) == [:(f(g(h(k(a))))),:(g(h(k(a)))),:(h(k(a))),:(k(a))]

   @test unify([:x,:y,:z,:u],:(P(x,y,z,u)),:(P(f(y),g(z),h(u),k(a)))) == [:(f(g(h(k(a))))),:(g(h(k(a)))),:(h(k(a))),:(k(a))]

   @testset "loop check for fixed point of subst" begin
     @test_throws Loop fp_subst([:x,:y],[:(f(y)),:x])
     @test_throws Loop fp_subst([:x,:y,:z],[:a,:(f(x,z)),:(g(y))]) == [:a,:(f(a,b)),:b]
     @test_throws Loop fp_subst([:x,:y,:z,:w],[:a,:(f(x,z)),:(g(w)),:(h(z))]) == [:a,:(f(a,b)),:b]
   end
 end

 @testset "unify intervention" begin
  @test unify([:x,:y],:(P(f(x),x)),:(P(y,a))) == [:a,:(f(a))]

  @test unify([:x,:y,:w],:(P(f(w,x),w,x)),:(P(y,h(b),g(a)))) == [:(g(a)),:(f(h(b),g(a))),:(h(b))]

  @test unify([:x,:y,:w,:u,:z,:n],:(P(f(w,z),w,g(n),z)),:(P(y,h(b),u,h(u)))) == [:x,:(f(h(b),h(g(n)))),:(h(b)),:(g(n)),:(h(g(n))),:n]

  @test unify([:x,:y,:z,:w,:u,:n,:v],:((P(f(x),y,h(y,z),z,k(n)))),:((P(w,g(w),u,q(v),v)))) == [:x,:(g(f(x))),:(q(k(n))),:(f(x)),:(h(g(f(x)),q(k(n)))),:n,:(k(n))]

  @test unify([:x,:y,:z,:w,:u,:n,:v],:((P(f(y,z),y,h(z),z,k(n)))),:((P(w,g(u,v),u,q(v),v)))) == [:x,:(g(h(q(k(n))),k(n))),:(q(k(n))),:(f(g(h(q(k(n))),k(n)),q(k(n)))),:(h(q(k(n)))),:n,:(k(n))]

 end

 @test unify([:x,:y], :(f(a=x,b=x)), :(f(a=p,b=y))) == [:p, :p]

 @testset "unify loop test fail" begin
  @test_throws Loop unify([:x,:y],:(P(x)),:(P(f(x))))
  @test_throws Loop unify([:x,:y],:(P(x,y)),:(P(y,f(y))))
 end

end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
# # Investment Science - Module 1
#
# \tableofcontents
#
# ## Problem 4
# The IRR is generally calculated using an iterative procedure.
# Suppose that we define $f(\lambda) = -a_{1} + a_{2}\lambda^{2} + \ldots + a_{n}\lambda^{n}$,
# where all $a_{i}$'s are positive and $n > 1$.
# Here is an iterative technique that generates a sequence
# $\lambda_{0}, \lambda_{1}, \lambda_{2}, \ldots \lambda_{k}, \ldots$
# of estimates that converges to the root $\overline{\lambda} > 0$, solving $f(\overline{\lambda}) = 0$.
# Start with any $\lambda_{0} > 0$ close to the solution.
# Assuming $\lambda_{k}$ has been calculated, evaluate
#
# $$
# f'(\lambda_{k}) = a_{1} + 2a_{2}\lambda_{k} + 3a_{3}\lambda_{k}^{2} + \ldots + na_{n}\lambda_{k}^{n-1}
# $$
#
# and define
#
# $$
# \lambda_{k+1} = \lambda_{k} - \frac{f(\lambda_{k})}{f'(\lambda_{k})}
# $$
#
# This is Newton's method.
# It is based on approximating the function $f$ by a line tangent to its graph at $\lambda_{k}$.
# Try the procedure on $f(\lambda) = -1 + \lambda + \lambda^{2}$.
# Start with $\lambda_{0} = 1$ and compute four additional estimates.
#
# ### Solution
#
# We know
#
# $$
# \begin{aligned}
# \lambda_{0} &= 1 \\
# \lambda_{k+1} &= \lambda_{k} - \frac{f(\lambda_{k})}{f'(\lambda_{k})} \\
# f(\lambda) &= -1 + \lambda + \lambda^{2} \\
# f'(\lambda) &= 1 + 2\lambda \\
# \end{aligned}
# $$
#
# Using all these equations we can compute four additional estimates.

## define function f
f(λ) = -1 + λ + λ^2
## define derivative of f
ḟ(λ) = 1 + 2*λ
## define function to calculate the next λ
next(λ) = λ - f(λ)/ḟ(λ)
## create a list of λs starting with 1
λs = [1.0]
## calculate next four
for n ∈ 1:4
    push!(λs, next(λs[n]))
end
## display results
@show λs

# ## Problem 5
#
# Suppose that you have the opportunity to plant trees that alter can be sold for lumber. This project requires an initial outlay of money in order to purchase and plant the seedlings. No other cash flow occurs until trees are harvested. However, you have a choice as to when to harvest. If you harvest after 1 year, you get your return quickly; but if you wait, the trees will have additional growth and the revenue generated from the sale of the trees will be greater. Assume that the case flow streams associated with these alternatives are:
#
# - Wait $1$ year: $(-1, 2)$
# - Wait $2$ year: $(-1, ,0, 3)$
# - Wait $3$ year: $(-1, 0, 0, 4)$
# - $\vdots$
# - Wait $n$ year: $(-1, 0, 0, \ldots, n+1)$
#
# The prevailing interest rate is 10%. When is the best time to cut the trees?
#
# ### Solution
#
# Let’s first define a function for calculating the Net Present Value, say on a cash flow stream that looks like $(x_{0}, \ldots, x_{n})$ where $r$ is the prevailing interest rate and $k$ is the number of years waited to harvest:
#
# $$
# NPV = -x_{0} + \frac{x_{k}}{(1+r)^{k}}
# $$
#
# We can create a function to calculate the **_NPV_** given a _cash flow stream_ and _interest rate_. Let's create a table of the first several values to determine the optimal time to harvest.

## define function to calculate net present value
function net_present_value(investment, revenue, year, rate)
    -investment + revenue/(1+rate)^year
end
## define empty array for values
values = []
## populate some values
for i ∈ 1:20
    append!(values, net_present_value(1, i+1, i, 0.10))
end
## display some values
@show values

# To make reading easier, we can create a table or `dataframe`.

using DataFrames
df = DataFrame(Dict(:year => 1:20, :revenue => values))
show(df)

# A quick glance reveals the maximum revenue.

maximum(values)
#+
using PyPlot;
figure(figsize=(8, 6));
plot(df.year, df.revenue);

xlabel("Year", fontsize=14);
ylabel("Revenue", fontsize=14);

plt.savefig(joinpath(@OUTPUT, "invsci-01-time-to-harvest.svg")); # hide

# \figalt{Time to Harvest}{invsci-01-time-to-harvest.svg}
#
# We can see that the graph levels off at the top, and upon inspection of the table we confirm that **years 9 and 10** will yield optimal revenue.

ext: jl

max_forks_repo_name: jvaverka/XranklinTest


text:
# # A vector of length 2 SVectors, denoting the face ID each edge is connected to. If the edge
# # is a boundary edge, face ID 0 is returned
# function edge_face_connectivity(mesh::UnstructuredMesh{Dim,Ord,T,U}) where {Dim,Ord,T,U}
#     # Each edge should only border 2 faces if it is an interior edge, and 1 face if it is
#     # a boundary edge.
#     # Loop through each face in the face_edge_connectivity vector and mark each edge with
#     # the faces that it borders.
#     if length(mesh.edges) === 0
#         @error "Does not have edges!"
#     end
#     if length(mesh.face_edge_connectivity) === 0
#         @error "Does not have face/edge connectivity!"
#     end
#     edge_face = [MVector{2, U}(0, 0) for _ in eachindex(mesh.edges)]
#     for (iface, edges) in enumerate(mesh.face_edge_connectivity)
#         for iedge in edges
#             # Add the face id in the first non-zero position of the edge_face conn. vec.
#             if edge_face[iedge][1] === U(0)
#                 edge_face[iedge][1] = iface
#             elseif edge_face[iedge][2] === U(0)
#                 edge_face[iedge][2] = iface
#             else
#                 @error "Edge $iedge seems to have 3 faces associated with it!"
#             end
#         end
#     end
#     return [SVector(sort!(two_faces).data) for two_faces in edge_face]
# end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_mono_m10(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_mono_m10!(A_copy)
end

@inline function exp_mono_m10!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=11
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    memslots9 = similar(A,T)
    memslots10 = similar(A,T)
    memslots11 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    value_one=ValueOne()
    # Computation order: B2 B3 B4 B5 B6 B7 B8 B9 B10 B11 T2k13
    # Computing B2 with operation: mult
    mul!(memslots2,memslots1,memslots1)
    # Computing B3 with operation: mult
    mul!(memslots3,memslots2,memslots1)
    # Computing B4 with operation: mult
    mul!(memslots4,memslots3,memslots1)
    # Computing B5 with operation: mult
    mul!(memslots5,memslots4,memslots1)
    # Computing B6 with operation: mult
    mul!(memslots6,memslots5,memslots1)
    # Computing B7 with operation: mult
    mul!(memslots7,memslots6,memslots1)
    # Computing B8 with operation: mult
    mul!(memslots8,memslots7,memslots1)
    # Computing B9 with operation: mult
    mul!(memslots9,memslots8,memslots1)
    # Computing B10 with operation: mult
    mul!(memslots10,memslots9,memslots1)
    # Computing B11 with operation: mult
    mul!(memslots11,memslots10,memslots1)
    # Computing T2k13 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7+x*B8+x*B9+x*B10+x*B11
    coeff1=1.0
    coeff2=1.0
    coeff3=0.5
    coeff4=0.16666666666666666
    coeff5=0.041666666666666664
    coeff6=0.008333333333333333
    coeff7=0.001388888888888889
    coeff8=0.0001984126984126984
    coeff9=2.48015873015873e-5
    coeff10=2.7557319223985893e-6
    coeff11=2.755731922398589e-7
    coeff12=2.505210838544172e-8
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots3 .+ coeff5.*memslots4 .+ coeff6.*memslots5 .+ coeff7.*memslots6 .+ coeff8.*memslots7 .+ coeff9.*memslots8 .+ coeff10.*memslots9 .+ coeff11.*memslots10 .+ coeff12.*memslots11
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 2
    # Deallocating B3 in slot 3
    # Deallocating B4 in slot 4
    # Deallocating B5 in slot 5
    # Deallocating B6 in slot 6
    # Deallocating B7 in slot 7
    # Deallocating B8 in slot 8
    # Deallocating B9 in slot 9
    # Deallocating B10 in slot 10
    # Deallocating B11 in slot 11
    return memslots1 # Returning T2k13
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


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

# RUN: julia --startup-file=no %s %t && llvm-link -S %t/* -o %t/module.ll
# RUN: cat %t/module.ll | FileCheck %s

## Notes:
# This script uses the `emit` function (defined llvmpasses.jl) to emit either
# optimized or unoptimized LLVM IR. Each function is emitted individually and
# `llvm-link` is used to create a single module that can be passed to opt.
# The order in which files are emitted and linked is important since `lit` will
# process the test cases in order.

include(joinpath("..", "testhelpers", "llvmpasses.jl"))

# CHECK-LABEL: @julia_simple_noinline
@noinline function simple_noinline(A, B)
    return A + B
end

# CHECK: attributes #{{[0-9]+}} = {{{[^}]*}} noinline {{[^}]*}}}
emit(simple_noinline, Float64, Float64)

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
struct ProductNumbers <: IncrementalInference.FunctorPairwise
  z::Distributions.Normal
end
getSample(s::ProductNumbers, N::Int=1) = (reshape(rand(s.z,N),1,:), )

function (s::ProductNumbers)(res::Array{Float64},
      userdata::FactorMetadata,
      idx::Int,
      meas::Tuple{Array{Float64,2}},
      X::Array{Float64,2},
      Y::Array{Float64,2},
      XY::Array{Float64,2}  )
  #
  res[1] = XY[1,idx] - X[1,idx]*Y[1,idx] + meas[1][1,idx]
  nothing
end

struct AreEqual <: IncrementalInference.FunctorPairwise
  z::Distributions.Normal
end
function getSample(s::AreEqual, N::Int=1)
  return (reshape(rand(s.z,N),1,:), )
end

function (s::AreEqual)(res::Array{Float64},
                       userdata::FactorMetadata,
                       idx::Int,
                       meas::Tuple{Array{Float64,2}},
                       X::Array{Float64,2},
                       Y::Array{Float64,2}  )
  #
  res[1] = X[1,idx]-Y[1,idx] + meas[1][1,idx]
  nothing
end

struct Square <: IncrementalInference.FunctorPairwise
  z::Distributions.Normal
end
getSample(s::Square, N::Int=1) = (reshape(rand(s.z,N),1,:), )

function (s::Square)(res::Array{Float64},
      userdata::FactorMetadata,
      idx::Int,
      meas::Tuple{Array{Float64,2}},
      X::Array{Float64,2},
      XX::Array{Float64,2}  )
  #
  res[1] = XX[1,idx] - X[1,idx]*X[1,idx] + meas[1][1,idx]
  nothing
end

mutable struct NumbersPrior <: IncrementalInference.FunctorSingleton
  z::BallTreeDensity
end
getSample(s::NumbersPrior, N::Int=1) = (reshape(rand(s.z,N),1,:), )


text:
using ModernGL, GeometryTypes, GLAbstraction, GLFW, Images, FileIO, Reactive, LinearAlgebra

const GLA = GLAbstraction

kitten = load(GLAbstraction.dir("tutorials", "images", "kitten.png"))
puppy  = load(GLAbstraction.dir("tutorials", "images", "puppy.png"))

# Create the window. This sets all the hints and makes the context current.
window = GLFW.Window(name="Depth and Stencils 1", resolution=(800,600),
                     windowhints=[(GLFW.DEPTH_BITS,32), (GLFW.STENCIL_BITS, 8)])

GLFW.MakeContextCurrent(window)
GLA.set_context!(window)

# The cube. This could be more efficiently represented using indexes,
# but the tutorial doesn't do it that way.
vertex_positions = Vec3f0[
    # The cube
    (-0.5f0, -0.5f0, -0.5f0),
    ( 0.5f0, -0.5f0, -0.5f0),
    ( 0.5f0,  0.5f0, -0.5f0),
    ( 0.5f0,  0.5f0, -0.5f0),
    (-0.5f0,  0.5f0, -0.5f0),
    (-0.5f0, -0.5f0, -0.5f0),

    (-0.5f0, -0.5f0,  0.5f0),
    ( 0.5f0, -0.5f0,  0.5f0),
    ( 0.5f0,  0.5f0,  0.5f0),
    ( 0.5f0,  0.5f0,  0.5f0),
    (-0.5f0,  0.5f0,  0.5f0),
    (-0.5f0, -0.5f0,  0.5f0),

    (-0.5f0,  0.5f0,  0.5f0),
    (-0.5f0,  0.5f0, -0.5f0),
    (-0.5f0, -0.5f0, -0.5f0),
    (-0.5f0, -0.5f0, -0.5f0),
    (-0.5f0, -0.5f0,  0.5f0),
    (-0.5f0,  0.5f0,  0.5f0),

    ( 0.5f0,  0.5f0,  0.5f0),
    ( 0.5f0,  0.5f0, -0.5f0),
    ( 0.5f0, -0.5f0, -0.5f0),
    ( 0.5f0, -0.5f0, -0.5f0),
    ( 0.5f0, -0.5f0,  0.5f0),
    ( 0.5f0,  0.5f0,  0.5f0),

    (-0.5f0, -0.5f0, -0.5f0),
    ( 0.5f0, -0.5f0, -0.5f0),
    ( 0.5f0, -0.5f0,  0.5f0),
    ( 0.5f0, -0.5f0,  0.5f0),
    (-0.5f0, -0.5f0,  0.5f0),
    (-0.5f0, -0.5f0, -0.5f0),

    (-0.5f0,  0.5f0, -0.5f0),
    ( 0.5f0,  0.5f0, -0.5f0),
    ( 0.5f0,  0.5f0,  0.5f0),
    ( 0.5f0,  0.5f0,  0.5f0),
    (-0.5f0,  0.5f0,  0.5f0),
    (-0.5f0,  0.5f0, -0.5f0)]

floor_positions = Vec3f0[
    # The floor
    (-1.0f0, -1.0f0, -0.5f0),
    ( 1.0f0, -1.0f0, -0.5f0),
    ( 1.0f0,  1.0f0, -0.5f0),
    ( 1.0f0,  1.0f0, -0.5f0),
    (-1.0f0,  1.0f0, -0.5f0),
    (-1.0f0, -1.0f0, -0.5f0)
]

vertex_texcoords = Vec2f0[
                          # The cube
                          (0.0f0, 0.0f0),
                          (1.0f0, 0.0f0),
                          (1.0f0, 1.0f0),
                          (1.0f0, 1.0f0),
                          (0.0f0, 1.0f0),
                          (0.0f0, 0.0f0),

                          (0.0f0, 0.0f0),
                          (1.0f0, 0.0f0),
                          (1.0f0, 1.0f0),
                          (1.0f0, 1.0f0),
                          (0.0f0, 1.0f0),
                          (0.0f0, 0.0f0),

                          (1.0f0, 0.0f0),
                          (1.0f0, 1.0f0),
                          (0.0f0, 1.0f0),
                          (0.0f0, 1.0f0),
                          (0.0f0, 0.0f0),
                          (1.0f0, 0.0f0),

                          (1.0f0, 0.0f0),
                          (1.0f0, 1.0f0),
                          (0.0f0, 1.0f0),
                          (0.0f0, 1.0f0),
                          (0.0f0, 0.0f0),
                          (1.0f0, 0.0f0),

                          (0.0f0, 1.0f0),
                          (1.0f0, 1.0f0),
                          (1.0f0, 0.0f0),
                          (1.0f0, 0.0f0),
                          (0.0f0, 0.0f0),
                          (0.0f0, 1.0f0),

                          (0.0f0, 1.0f0),
                          (1.0f0, 1.0f0),
                          (1.0f0, 0.0f0),
                          (1.0f0, 0.0f0),
                          (0.0f0, 0.0f0),
                          (0.0f0, 1.0f0)]

floor_texcoords = Vec2f0[
                          # The floor
                          (0.0f0, 0.0f0),
                          (1.0f0, 0.0f0),
                          (1.0f0, 1.0f0),
                          (1.0f0, 1.0f0),
                          (0.0f0, 1.0f0),
                          (0.0f0, 0.0f0)]

vertex_colors = fill(Vec3{Float32}(1,1,1), 36)
floor_colors = fill(Vec3{Float32}(0,0,0),6)

vertex_shader = GLA.vert"""
#version 150

in vec3 position;
in vec3 color;
in vec2 texcoord;

out vec3 Color;
out vec2 Texcoord;

uniform vec3 overrideColor;
uniform mat4 model;
uniform mat4 view;
uniform mat4 proj;

void main()
{
    Color = overrideColor * color;
    Texcoord = texcoord;
    gl_Position = proj * view * model * vec4(position, 1.0);
}
"""
fragment_shader = load(joinpath(dirname(@__FILE__), "shaders", "puppykitten_color.frag"))
fragment_shader = load(joinpath("tutorials", "shaders", "puppykitten_color.frag"))
prog = GLA.Program(GLA.Shader(vertex_shader), fragment_shader)

model1 = Mat4(diagm(0=>ones(Float32, 4)))
model2 = Mat4{Float32}(
        1, 0, 0, 0,
        0, 1, 0, 0,
        0, 0, -1, 0,
        0,0,-1,1.0f0
    )

function lookat(eyePos, lookAt, up)
    z  = normalize(eyePos-lookAt)
    x  = normalize(cross(up, z))
    y  = normalize(cross(z, x))
    T0 = 0.0f0
    return Mat4{Float32}(
        x[1], y[1], z[1], T0,
        x[2], y[2], z[2], T0,
        x[3], y[3], z[3], T0,
        -dot(x,eyePos),-dot(y,eyePos),-dot(z,eyePos),1.0f0
    )
end

function perspectiveprojection(fovy, aspect, znear::T, zfar::T) where T
    h = T(tan(deg2rad(fovy)) * znear)
    w = T(h * aspect)
    bottom = -h
    top = h
    left = -w
    right = w

    (right == left || bottom == top || znear == zfar) && return Mat4{T}(I)
    T0, T1, T2 = zero(T), one(T), T(2)
    return Mat4{T}(
        T2 * znear / (right - left), T0, T0, T0,
        T0, T2 * znear / (top - bottom), T0, T0,
        (right + left) / (right - left), (top + bottom) / (top - bottom), -(zfar + znear) / (zfar - znear), -T1,
        T0, T0, -(T2 * znear * zfar) / (zfar - znear), T0
    )
end

view = lookat(Vec3((2.5, 2.5, 2)), Vec3((0, 0, 0)), Vec3((0, 0, 1)))
proj = perspectiveprojection(45f0, 600f0/600f0, 1f0, 10f0)

tex_kitten = GLA.Texture(collect(kitten'))
tex_puppy  = GLA.Texture(collect(puppy'))

## Now render the distinct objects. Rather than always using std_renderobject,
## here we control the settings manually.
# The cube
buffers_cube = GLA.generate_buffers(prog, GLA.GEOMETRY_DIVISOR, position = vertex_positions,
                                                           color = vertex_colors,
                                                           texcoord=vertex_texcoords)

vao_cube = GLA.VertexArray(buffers_cube)
# The floor. This is drawn without writing to the depth buffer, but we
# write stencil values.
buffers_floor = GLA.generate_buffers(prog, GLA.GEOMETRY_DIVISOR, position = floor_positions,
                                                           color = floor_colors,
                                                           texcoord=floor_texcoords)
vao_floor = GLA.VertexArray(buffers_floor)

glEnable(GL_DEPTH_TEST)
glDepthFunc(GL_LEQUAL)
glClearColor(1,1,1,1) # make the background white, so we can see the floor
glClearStencil(0)     # clear the stencil buffer with 0

while !GLFW.WindowShouldClose(window)
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    GLA.bind(prog)

    GLA.gluniform(prog, :texKitten, 0, tex_kitten) #first texture sampler
    GLA.gluniform(prog, :texPuppy, 1, tex_puppy) # second texture sampler

    # Render Cube
    GLA.gluniform(prog, :model, model1)
    GLA.gluniform(prog, :view, view)
    GLA.gluniform(prog, :proj, proj)
    GLA.gluniform(prog, :overrideColor, Vec3{Float32}((1,1,1)))

    glDisable(GL_STENCIL_TEST)
    GLA.bind(vao_cube)
    GLA.draw(vao_cube)

    # Render Floor
    GLA.gluniform(prog, :model, model1)
    glDepthMask(GL_FALSE)                  # don't write to depth buffer
    glEnable(GL_STENCIL_TEST)              # use stencils
    glStencilMask(0xff)                    # do write to stencil buffer
    glStencilFunc(GL_ALWAYS, 1, 0xff)      # all pass
    glStencilOp(GL_KEEP, GL_KEEP, GL_REPLACE)  # replace stencil value
    glClear(GL_STENCIL_BUFFER_BIT)         # start with empty buffer
    GLA.bind(vao_floor)
    GLA.draw(vao_floor)

    # Render reflection
    GLA.gluniform(prog, :model, model2)
    GLA.gluniform(prog, :overrideColor, Vec3{Float32}((0.3,0.3,0.3)))

    glStencilFunc(GL_EQUAL, 1, 0xff)
    glStencilMask(0x00)

    GLA.bind(vao_cube)
    GLA.draw(vao_cube)

    GLFW.SwapBuffers(window)
    GLFW.PollEvents()
    if GLFW.GetKey(window, GLFW.KEY_ESCAPE) == GLFW.PRESS
        GLFW.SetWindowShouldClose(window, true)
    end
end
GLFW.DestroyWindow(window)  # needed if you're running this from the REPL


text:
################################################################################
#
#  Abstract types
#
################################################################################

# abstract spaces
abstract type AbsSpace{S} end

# abstract lattices
abstract type AbsLat{S} end

################################################################################
#
#  Quadratic spaces
#
################################################################################

@attributes mutable struct QuadSpace{S, T} <: AbsSpace{S}
  K::S
  gram::T

  function QuadSpace(K::S, G::T) where {S, T}
    # I also need to check if the gram matrix is Hermitian
    if dense_matrix_type(elem_type(S)) === T
      z = new{S, T}(K, G)
    else
      try
        Gc = change_base_ring(K, G)
        if typeof(Gc) !== dense_matrix_type(elem_type(S))
          error("Cannot convert entries of the matrix to the number field")
        end
        @assert base_ring(Gc) === K
        z = new{S, dense_matrix_type(elem_type(S))}(K, Gc)
        return z
      catch e
        rethrow(e)
      end
    end
  end
end

################################################################################
#
#  Hermitian spaces
#
################################################################################

@attributes mutable struct HermSpace{S, T, U, W} <: AbsSpace{S}
  E::S
  K::T
  gram::U
  involution::W

  function HermSpace(E::S, gram::U) where {S, U}
    # I also need to check if the gram matrix is Hermitian
    if dense_matrix_type(elem_type(S)) === U
      gramc = gram
    else
      try
        gramc = change_base_ring(E, gram)
        if typeof(gramc) !== dense_matrix_type(elem_type(S))
          error("Cannot convert entries of the matrix to the number field")
        end
      catch e
        rethrow(e)
      end
    end

    @assert degree(E) == 2
    A = automorphisms(E)
    a = gen(E)
    if A[1](a) == a
      involution = A[2]
    else
      involution = A[1]
    end

    K = base_field(E)

    z = new{S, typeof(K), dense_matrix_type(elem_type(S)), typeof(involution)}(E, K, gramc, involution)
    return z
  end
end


text:
function IntersectedDomain{T}(
    domain::IntersectedDomain{T},
    parent::ArbField,
) where {T<:AbstractPlanarDomain}
    return IntersectedDomain(
        T(domain.exterior, parent),
        [typeof(d)(d, parent) for d in domain.interiors],
    )
end

IntersectedDomain(exterior::AbstractPlanarDomain, interior::AbstractPlanarDomain) =
    IntersectedDomain(exterior, [interior])

function Base.show(io::IO, domain::IntersectedDomain)
    println(io, "Intersected domain")
    if !get(io, :compact, false)
        println(io, "Exterior: $(domain.exterior)")
        println(io, "Interiors:")
        recur_io = IOContext(io, :compact => true)
        for d in domain.interiors
            println(recur_io, d)
        end
    end
end

function Base.getproperty(domain::IntersectedDomain, name::Symbol)
    if name == :parent
        return domain.exterior.parent
    else
        return getfield(domain, name)
    end
end

exterior_boundaries(domain::IntersectedDomain) = boundaries(domain.exterior)

function interior_boundaries(domain::IntersectedDomain)
    boundaries_exterior = boundaries(domain.exterior)
    boundaries_interior = boundaries.(domain.interiors)
    return (1:sum(length.(boundaries_interior))) .+ maximum(boundaries_exterior)
end

boundaries(domain::IntersectedDomain) =
    1:(length(exterior_boundaries(domain))+length(interior_boundaries(domain)))

function get_domain_and_boundary(domain::IntersectedDomain, i::Integer)
    boundaries_exterior = boundaries(domain.exterior)
    if i ∈ boundaries_exterior
        return domain.exterior, i
    end

    i_local = i - maximum(boundaries_exterior)
    for d in domain.interiors
        bs = boundaries(d)
        if i_local ∈ bs
            return (d, i_local)
        end
        i_local -= maximum(bs)
    end

    throw(ArgumentError("attempt to get boundary $i from a $(typeof(domain))"))
end

function vertex(domain::IntersectedDomain, i::Integer)
    d, j = get_domain_and_boundary(domain, i)
    return vertex(d, j)
end

vertices(domain::IntersectedDomain) =
    vcat(collect(vertices(domain.exterior)), collect.(vertices.(domain.interiors))...)

center(domain::IntersectedDomain) = center(domain.exterior)

area(domain::IntersectedDomain) = area(domain.exterior) - sum(area.(domain.interiors))

Base.in(xy, domain::IntersectedDomain) =
    xy ∈ domain.exterior && !any(d -> xy ∈ d, domain.interiors)

boundary_parameterization(t, domain::IntersectedDomain, i::Integer) =
    boundary_parameterization(t, get_domain_and_boundary(domain, i)...)

boundary_points(
    domain::IntersectedDomain,
    i::Integer,
    n::Integer;
    distribution = :chebyshev,
) = (boundary_points(get_domain_and_boundary(domain, i)..., n; distribution)[1], fill(i, n))

# TODO: Exclude points landing in the interior domain
function interior_points(domain::IntersectedDomain, n::Integer; rng = MersenneTwister(42))
    # We generate points in the interior of domain.exterior and remove
    # those which are also in domain.interior. The issue is that we
    # then don't know how many interior points we need to generate.
    gotten = 0
    points = SVector{2,arb}[]
    while gotten < n
        new_points = Iterators.take(
            filter(
                xy -> !any(d -> xy ∈ d, domain.interiors),
                interior_points(domain.exterior, 2(n - gotten), rng = rng),
            ),
            n - gotten,
        )

        gotten += length(new_points)
        append!(points, new_points)
    end

    return points
end


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function test_epitrochoid(fname)
    Drawing(1400, 1400, fname)
    origin()
    background("grey20")
    sethue("antiquewhite")
    setopacity(0.8)
    tiles = Tiler(1200, 1200, 6, 6)
    for (pos, n) in tiles
        gsave()
        setline(.5)
        translate(pos)

        R = 20; r = n; d = 20
        sethue("antiquewhite")
        poly(epitrochoid(R, n, d, :stroke, stepby=pi/300, vertices=true), :stroke, close=true)
        text(string("R = $(R), r = $(r), d = $(d)"), halign=:center, O.x, O.y+tiles.tileheight/2)

        # test perimeterlength
        ptslist = epitrochoid(R, n, d, :stroke, stepby=pi/300, vertices=true)
        polyperimeter(ptslist)

        R = 20; r = 12; d = 5
        sethue("orange")
        (R == r) && (r += 1)
        p = epitrochoid(R, r, d, :none, stepby=pi/305, vertices=true)
        polydistances(p)
        text(string("R = $(R), r = $(r), d = $(d)"), halign=:center, O.x, O.y+tiles.tileheight/2 + 10)
        poly(p, :stroke, close=true)

        # testing portions of polys
        sethue("magenta")
        portion = rand()
        poly(polyportion(p, portion), :stroke, close=false)
        text(string("portion ", round(portion, digits=2)), halign=:center, O.x, O.y+tiles.tileheight/2 + 30)

        sethue("cyan")
        setline(0.3)
        poly(offsetpoly(p, 5), :stroke, close=true)
        text("offset 15", halign=:center, O.x, O.y+tiles.tileheight/2 + 20)

        grestore()
    end
    @test finish() == true
end

Random.seed!(42)

fname = "epitrochoid-test.pdf"
test_epitrochoid(fname)
println("...finished test: output in $(fname)")


text:
export MolOrbital, getMolOrbitals, Molecule, nnRepulsions

using LinearAlgebra: norm

"""

    MolOrbital{N} <: AbstractMolOrbital

`Struct` of molecular orbitals.

≡≡≡ Field(s) ≡≡≡

`symmetry::String`: The symmetry of the orbital. The default value is "A" for being
antisymmetric.

`energy::Float64`: Molecular energy.

`spin::String`: Spin function of the orbital. Available values: "Alpha", "Beta".

`occupancy::Real`: Occupation number.

`orbitalCoeffs::NTuple{N, Float64}`: coefficients of the basis functions to form the
molecular orbital.

≡≡≡ Initialization Method(s) ≡≡≡

    MolOrbital(energy::Float64, occupancy::Real, orbitalCoeffs::Array{Float64, 1},
               spin::String="Alpha", symmetry::String="A") -> MolOrbital{N}

"""
struct MolOrbital{N} <: AbstractMolOrbital
    symmetry::String
    energy::Float64
    spin::String
    occupancy::Real
    orbitalCoeffs::NTuple{N, Float64}

    function MolOrbital(energy::Float64, occupancy::Real, orbitalCoeffs::Vector{Float64},
                        spin::String="Alpha", symmetry::String="A")
        spin != "Alpha" && spin != "Beta" && error("Keyword argument \"spin\" can only"*
                                                   " be \"Alpha\" or \"Beta\"")
        new{length(orbitalCoeffs)}(symmetry, energy, spin, occupancy, orbitalCoeffs|>Tuple)
    end
end

"""

    getMolOrbitals(ens::Array{Float64, 1}, occus::Array{<:Real, 1}, C::Matrix{Float64},
                   spins::Array{String, 1},
                   symms::Array{String, 1}=repeat(["A"], length(occus))) ->
    Tuple{Vararg{getMolOrbitals}}

A function that returns the molecular orbitals.
"""
function getMolOrbitals(ens::Vector{Float64}, occus::Vector{<:Real}, C::Matrix{Float64},
                        spins::Vector{String},
                        symms::Vector{String}=repeat(["A"], length(occus)))
    MolOrbital.(ens, occus, [ C[:,i] for i = 1:size(C, 2) ], spins, symms) |> Tuple
end

"""

    Molecule{Nc, Ne, Nb} <:MolecularHartreeFockCoefficient{Nc, Ne}

Container for the information of a molecule.

≡≡≡ Field(s) ≡≡≡

`nuc::Tuple{Vararg{String}}`: Nuclei of the molecule.

`nucCoords::Tuple{Vararg{NTuple{3, Real}}}`: The coordinates of the nuclei.

`Ne::Int`: The number of electrons.

`orbital::Tuple{Vararg{MolOrbital}}`: Molecular orbitals.

`basis::Tuple{Vararg{FloatingGTBasisFuncs}}`: The basis set for the molecular orbitals.

`E0HF::Float64`: Hartree-Fock energy of the electronic Hamiltonian from the basis set.

`EnnR::Float64`: The nuclear-nuclear repulsion energy.

≡≡≡ Initialization Method(s) ≡≡≡

    Molecule(basis::Array{FloatingGTBasisFuncs, 1}, nuc::Array{String, 1},
             nucCoords::Array{<:AbstractArray, 1}, Ne::Int, E0HF::Float64,
             Emos::Array{Float64, 1}, occus::Array{<:Real, 1}, C::Array{Float64, 2},
             spins::Array{String, 1},
             symms::Array{String, 1}=repeat(["A"], length(occus))) ->
    Molecule{<:Any, Ne, <:Any}

`Emos` are the energies of corresponding molecular energies. `occus` are the occupation
numbers of the orbitals. `C` is the coefficient matrix, which does not need to be a square
matrix since the number of rows is the size of the (spatial) basis set whereas the number
of the columns represents the number of molecular orbitals. `spin` specifies the spin
functions of the orbitals, entries of which can be set to "Alpha" or "Beta". `symms` are
symmetries of the orbitals where the default entry value is "A" for being antisymmetric.

    Molecule(basis::Array{<:FloatingGTBasisFuncs, 1}, nuc::Array{String, 1},
             nucCoords::Array{<:AbstractArray, 1}, HFfVars::HFfinalVars) ->
    Molecule

Construct a `Molecule` from a basis set, nuclei information, and the result from the
corresponding Hartree-Fock SCF procedure, specifically a `HFfinalVars` `struct`.
"""
struct Molecule{Nc, Ne, Nb} <:MolecularHartreeFockCoefficient{Nc, Ne}
    nuc::Tuple{Vararg{String}}
    nucCoords::Tuple{Vararg{NTuple{3, Real}}}
    Ne::Int
    orbital::Tuple{Vararg{MolOrbital}}
    basis::Tuple{Vararg{FloatingGTBasisFuncs}}
    E0HF::Float64
    EnnR::Float64

    function Molecule(basis::Vector{<:FloatingGTBasisFuncs}, nuc::Vector{String},
                      nucCoords::Vector{<:AbstractArray}, Ne::Int, E0HF::Float64,
                      Emos::Vector{Float64}, occus::Vector{<:Real}, C::Matrix{Float64},
                      spins::Vector{String},
                      symms::Vector{String}=repeat(["A"], length(occus)))
        @assert length(nuc) == length(nucCoords)
        Nb = basisSize(basis) |> sum
        coeff = spins |> unique |> length
        @assert (coeff*Nb .==
                 length.([Emos, occus, C[1,:], spins, symms]) .==
                 coeff*length(C[:,1])) |> prod
        EnnR = nnRepulsions(nuc, nucCoords)
        iSorted = mapPermute(basis, sortBasisFuncs)
        basis = basis[iSorted]
        offset = pushfirst!(accumulate(+, collect(basisSize(basis)) .- 1), 0)
        iSortedExtended = Union{Int, Vector{Int}}[]
        append!(iSortedExtended, iSorted)
        for (i, idx) in zip(iSorted, 1:length(iSorted))
            if basis[i] isa BasisFuncs
                iSortedExtended[idx] = collect(i : i+length(basis[i])-1) .+ offset[i]
            else
                iSortedExtended[idx] += offset[i]
            end
        end
        C = C[vcat(iSortedExtended...), :]
        new{getCharge(nuc), Ne, Nb}(nuc |> Tuple,
                                    nucCoords .|> Tuple |> Tuple,
                                    Ne,
                                    getMolOrbitals(Emos, occus, C, spins, symms),
                                    deepcopy(basis) |> Tuple,
                                    E0HF,
                                    EnnR)
    end
end

function Molecule(basis::Vector{<:FloatingGTBasisFuncs}, nuc::Vector{String},
                  nucCoords::Vector{<:AbstractArray}, HFfVars::HFfinalVars)
    t = typeof(HFfVars)
    len = t.parameters[3]
    if t <: HFfinalVars{:UHF}
        Emos = HFfVars.Emo |> flatten
        occus = HFfVars.occu |> flatten
        C = hcat(HFfVars.C...)
        spins = vcat(fill("Alpha", len), fill("Beta", len))
    else
        Emos = HFfVars.Emo
        occus = HFfVars.occu
        C = HFfVars.C
        spins = fill("Alpha", len)
    end
    Molecule(basis, nuc, nucCoords, t.parameters[2], HFfVars.E0HF, Emos|>collect,
             occus|>collect, C, spins, fill("A", length(spins)))
end

"""

    nnRepulsions(nuc::Array{String, 1}, nucCoords::Array{<:AbstractArray, 1}) -> Float64

Calculate the nuclear-nuclear repulsion energy given the nuclei and their coordinates of a
molecule.
"""
function nnRepulsions(nuc::Vector{String}, nucCoords::Vector{<:AbstractArray})
    E = 0.0
    len = length(nuc)
    for i = 1:len, j=i+1:len
        E += getCharge(nuc[i]) * getCharge(nuc[j]) / norm(nucCoords[i] - nucCoords[j])
    end
    E
end


text:
module DocHosting

"""
    greet()

Output "Hello World!" onto the console.
"""
greet() = print("Hello World!")

export greet

end # module


text:
# This file contains the function generator tools to drive other tools of DsSimulator.

import UUIDs: uuid4

"""
    @def_source ex

where `ex` is the expression to define to define a new AbstractSource component type. The usage is as follows:
```julia
@def_source struct MySource{T1,T2,T3,...,TN,OP, RO} <: AbstractSource
    param1::T1 = param1_default     # optional field
    param2::T2 = param2_default     # optional field
    param3::T3 = param3_default     # optional field
        ⋮
    paramN::TN = paramN_default     # optional field
    output::OP = output_default     # mandatory field
    readout::RO = readout_function  # mandatory field
end
```
Here, `MySource` has `N` parameters, an `output` port and a `readout` function.

!!! warning
    `output` and `readout` are mandatory fields to define a new source. The rest of the fields are the parameters of the source.

!!! warning
    `readout` must be a single-argument function, i.e. a fucntion of time `t`.

!!! warning
    New source must be a subtype of `AbstractSource` to function properly.

# Example
```julia
julia> @def_source struct MySource{OP, RO} <: AbstractSource
       a::Int = 1
       b::Float64 = 2.
       output::OP = Outport()
       readout::RO = t -> (a + b) * sin(t)
       end

julia> gen = MySource();

julia> gen.a
1

julia> gen.output
1-element Outport{Outpin{Float64}}:
 Outpin(eltype:Float64, isbound:false)
```
"""
macro def_source(ex)
    ex.args[2].head == :(<:) && ex.args[2].args[2] == :AbstractSource ||
        error("Invalid usage. The type should be a subtype of AbstractSource.\n$ex")
    foreach(nex -> appendex!(ex, nex), [
        :( trigger::$TRIGGER_TYPE_SYMBOL = Inpin() ),
        :( handshake::$HANDSHAKE_TYPE_SYMBOL = Outpin{Bool}() ),
        :( callbacks::$CALLBACKS_TYPE_SYMBOL = nothing ),
        :( name::Symbol = Symbol() ),
        :( id::$ID_TYPE_SYMBOL = Causal.uuid4() )
        ])
    quote
        Base.@kwdef $ex
    end |> esc
end

##### Define Sources library
"""
    FunctionGenerator(; readout, output=Outport())

Constructs a generic function generator with `readout` function and `output` port.

# Example
```jldoctest
julia> gen = FunctionGenerator(readout = t -> [t, 2t], output = Outport(2));

julia> gen.readout(1.)
2-element Array{Float64,1}:
 1.0
 2.0
```
"""
@def_source struct FunctionGenerator{RO, OP} <: AbstractSource
    readout::RO
    output::OP = Outport(1)
end

 @doc raw"""
    SinewaveGenerator(;amplitude=1., frequency=1., phase=0., delay=0., offset=0.)

Constructs a `SinewaveGenerator` with output of the form
```math
    x(t) = A sin(2 \pi f  (t - \tau) + \phi) + B
```
where ``A`` is `amplitude`, ``f`` is `frequency`, ``\tau`` is `delay` and ``\phi`` is `phase` and ``B`` is `offset`.
"""
@def_source struct SinewaveGenerator{RO, OP} <: AbstractSource
    amplitude::Float64 = 1.
    frequency::Float64 = 1.
    phase::Float64 = 0.
    delay::Float64 = 0.
    offset::Float64 = 0.
    output::OP = Outport()
    readout::RO = (t, amplitude=amplitude, frequency=frequency, delay=delay, offset=offset) ->
        amplitude * sin(2 * pi * frequency * (t - delay) + phase) + offset
end

@doc raw"""
    DampedSinewaveGenerator(;amplitude=1., decay=-0.5, frequency=1., phase=0., delay=0., offset=0.)

Constructs a `DampedSinewaveGenerator` which generates outputs of the form
```math
    x(t) = A e^{\alpha t} sin(2 \pi f (t - \tau) + \phi) + B
```
where ``A`` is `amplitude`, ``\alpha`` is `decay`, ``f`` is `frequency`, ``\phi`` is `phase`, ``\tau`` is `delay` and ``B`` is `offset`.
"""
@def_source struct DampedSinewaveGenerator{RO, OP} <: AbstractSource
    amplitude::Float64 = 1.
    decay::Float64 = 0.5
    frequency::Float64 = 1.
    phase::Float64 = 0.
    delay::Float64 = 0.
    offset::Float64 = 0.
    output::OP = Outport()
    readout::RO = (t, amplitude=amplitude, decay=decay, frequency=frequency, phase=phase, delay=delay, offset=offset) ->
        amplitude * exp(decay * t) * sin(2 * pi * frequency * (t - delay)) + offset
end

@doc raw"""
    SquarewaveGenerator(;level1=1., level2=0., period=1., duty=0.5, delay=0.)

Constructs a `SquarewaveGenerator` with output of the form
```math
    x(t) = \left\{\begin{array}{lr}
	A_1 + B, &  kT + \tau \leq t \leq (k + \alpha) T + \tau \\
	A_2 + B,  &  (k + \alpha) T + \tau \leq t \leq (k + 1) T + \tau
	\end{array} \right. \quad k \in Z
```
where ``A_1``, ``A_2`` is `level1` and `level2`, ``T`` is `period`, ``\tau`` is `delay` ``\alpha`` is `duty`.
"""
@def_source struct SquarewaveGenerator{OP, RO} <: AbstractSource
    high::Float64 = 1.
    low::Float64 = 0.
    period::Float64 = 1.
    duty::Float64 = 0.5
    delay::Float64 = 0.
    output::OP = Outport()
    readout::RO = (t, high=high, low=low, period=period, duty=duty, delay=delay) ->
        t <= delay ? low : ( ((t - delay) % period <= duty * period) ? high : low )
end

@doc raw"""
    TriangularwaveGenerator(;amplitude=1, period=1, duty=0.5, delay=0, offset=0)

Constructs a `TriangularwaveGenerator` with output of the form
```math
    x(t) = \left\{\begin{array}{lr}
	\dfrac{A t}{\alpha T} + B, &  kT + \tau \leq t \leq (k + \alpha) T + \tau \\[0.25cm]
	\dfrac{A (T - t)}{T (1 - \alpha)} + B,  &  (k + \alpha) T + \tau \leq t \leq (k + 1) T + \tau
	\end{array} \right. \quad k \in Z
```
where ``A`` is `amplitude`, ``T`` is `period`, ``\tau`` is `delay` ``\alpha`` is `duty`.
"""
@def_source struct TriangularwaveGenerator{OP, RO} <: AbstractSource
    amplitude::Float64 =  1.
    period::Float64 = 1.
    duty::Float64 = 0.5
    delay::Float64 = 0.
    offset::Float64 = 0.
    output::OP = Outport()
    readout::RO = (t, amplitude=amplitude, period=period, duty=duty, delay=delay, offset=offset) -> begin
        if t <= delay
            return offset
        else
            t = (t - delay) % period
            if t <= duty * period
                amplitude / (duty * period) * t + offset
            else
                (amplitude * (period - t)) / (period * (1 - duty)) + offset
            end
        end
    end
end

@doc raw"""
    ConstantGenerator(;amplitude=1.)

Constructs a `ConstantGenerator` with output of the form
```math
    x(t) = A
```
where ``A`` is `amplitude.
"""
@def_source struct ConstantGenerator{OP, RO} <: AbstractSource
    amplitude::Float64 = 1.
    output::OP = Outport()
    readout::RO = (t, amplitude=amplitude) -> amplitude
end

@doc raw"""
    RampGenerator(;scale=1, delay=0.)

Constructs a `RampGenerator` with output of the form
```math
    x(t) = \alpha (t - \tau)
```
where ``\alpha`` is the `scale` and ``\tau`` is `delay`.
"""
@def_source struct RampGenerator{OP, RO} <: AbstractSource
    scale::Float64 = 1.
    delay::Float64 = 0.
    offset::Float64 = 0.
    output::OP = Outport()
    readout::RO = (t, scale=scale, delay=delay, offset=offset) ->  scale * (t - delay) + offset
end

@doc raw"""
    StepGenerator(;amplitude=1, delay=0, offset=0)

Constructs a `StepGenerator` with output of the form
```math
    x(t) = \left\{\begin{array}{lr}
	B, &  t \leq \tau  \\
	A + B,  &  t > \tau
	\end{array} \right.
```
where ``A`` is `amplitude`, ``B`` is the `offset` and ``\tau`` is the `delay`.
"""
@def_source struct StepGenerator{OP, RO} <: AbstractSource
    amplitude::Float64 = 1.
    delay::Float64 = 0.
    offset::Float64 = 0.
    output::OP = Outport()
    readout::RO = (t, amplitude=amplitude, delay=delay, offset=offset) ->
        t - delay >= 0 ? amplitude + offset : offset
end

@doc raw"""
    ExponentialGenerator(;scale=1, decay=-1, delay=0.)

Constructs an `ExponentialGenerator` with output of the form
```math
    x(t) = A e^{\alpha (t - \tau)}
```
where ``A`` is `scale`, ``\alpha`` is `decay` and ``\tau`` is `delay`.
"""
@def_source struct ExponentialGenerator{OP, RO} <: AbstractSource
    scale::Float64 = 1.
    decay::Float64 = -1.
    delay::Float64 = 0.
    offset::Float64 = 0.
    output::OP = Outport()
    readout::RO = (t, scale=scale, decay=decay, delay=delay, offset=offset) -> scale * exp(decay * (t - delay)) + offset
end

@doc raw"""
    DampedExponentialGenerator(;scale=1, decay=-1, delay=0.)

Constructs an `DampedExponentialGenerator` with outpsuts of the form
```math
    x(t) = A (t - \tau) e^{\alpha (t - \tau)}
```
where ``A`` is `scale`, ``\alpha`` is `decay`, ``\tau`` is `delay`.
"""
@def_source struct DampedExponentialGenerator{OP, RO} <: AbstractSource
    scale::Float64 = 1.
    decay::Float64 = -1.
    delay::Float64 = 0.
    offset::Float64 = 0.
    output::OP = Outport()
    readout::RO = (t, scale=scale, decay=decay, delay=delay, offset=offset) ->
        scale * (t - delay) * exp(decay * (t - delay)) + offset
end

##### Pretty-Printing of generators.

show(io::IO, gen::FunctionGenerator) = print(io,
    "FunctionGenerator(readout:$(gen.readout),  output:$(gen.output))")
show(io::IO, gen::SinewaveGenerator) = print(io,
    "SinewaveGenerator(amp:$(gen.amplitude), freq:$(gen.frequency), phase:$(gen.phase), ",
    "offset:$(gen.offset), delay:$(gen.delay))")
show(io::IO, gen::DampedSinewaveGenerator) = print(io,
    "DampedSinewaveGenerator(amp:$(gen.amplitude), decay:$(gen.delay), freq:$(gen.frequency), ",
    "phase:$(gen.phase), delay:$(gen.delay), offset:$(gen.offset))")
show(io::IO, gen::SquarewaveGenerator) = print(io,
    "SquarewaveGenerator(high:$(gen.high), low:$(gen.low), period:$(gen.period), duty:$(gen.duty), ",
    "delay:$(gen.delay))")
show(io::IO, gen::TriangularwaveGenerator) = print(io,
    "TriangularwaveGenerator(amp:$(gen.amplitude), period:$(gen.period), duty:$(gen.duty), ",
    "delay:$(gen.delay), offset:$(gen.offset))")
show(io::IO, gen::ConstantGenerator) = print(io,
    "ConstantGenerator(amp:$(gen.amplitude))")
show(io::IO, gen::RampGenerator) = print(io, "RampGenerator(scale:$(gen.scale), delay:$(gen.delay))")
show(io::IO, gen::StepGenerator) = print(io,
    "StepGenerator(amp:$(gen.amplitude), delay:$(gen.delay), offset:$(gen.offset))")
show(io::IO, gen::ExponentialGenerator) = print(io,
    "ExponentialGenerator(scale:$(gen.scale), decay:$(gen.decay), delay:$(gen.delay))")
show(io::IO, gen::DampedExponentialGenerator) = print(io,
    "DampedExponentialGenerator(scale:$(gen.scale), decay:$(gen.decay), delay:$(gen.delay))")


text:
module Multigrid
using jInv.Mesh;
using jInv.LinearSolvers;
using KrylovMethods

# check if MUMPS can be used
const minMUMPSversion = VersionNumber(0,0,1)
hasMUMPS=false
vMUMPS = VersionNumber(0,0,0)
try
	vMUMPS = Pkg.installed("MUMPS")
	hasMUMPS = vMUMPS >= minMUMPSversion
	if hasMUMPS
		using MUMPS;
	end
catch
end

# check if ParSPMatVec is available
hasParSpMatVec = false
const minVerParSpMatVec = VersionNumber(0,0,1)
	vParSpMatVec = VersionNumber(0,0,0)
try
 	using ParSpMatVec
 	hasParSpMatVec = true
 	if myid()==1
 		vParSpMatVec = Pkg.installed("ParSpMatVec")
 		if vParSpMatVec < minVerParSpMatVec;
 			warn("ParSpMatVec is outdated! Please update to version $(minVerParSpMatVec) or higher.")
 		end
 	end
catch
end
if hasParSpMatVec
	println("Using ParSpMatVec")
	using ParSpMatVec
end

export MGparam;
export getMGparam, MGsetup, clear!
export BlockFGMRES, ArrayTypes,hierarchyExists,copySolver,destroyCoarsestLU
export multilevelOperatorConstructor, getMultilevelOperatorConstructor

SparseCSCTypes = Union{SparseMatrixCSC{Complex128,Int64},SparseMatrixCSC{Float64,Int64},SparseMatrixCSC{Complex64,Int32},SparseMatrixCSC{Float32,Int32}}
ArrayTypes = Union{Array{Complex128},Array{Complex64},Array{Float64},Array{Float32}}

include("SpMatMul.jl");
include("FGMRES.jl");
include("BlockFGMRES.jl");

"""
type Multigrid.multilevelOperatorConstructor

Fields:

	param - parameters for the PDE
	getOperator(mesh,param) - a function that creates the PDE from
	restrictParams(mesh_fine,mesh_coarse,param_fine,level) - returns param_coarse for level+1.

"""

type multilevelOperatorConstructor
	param
	getOperator		::Function
	restrictParams  ::Function
end

function getMultilevelOperatorConstructor(param,getOperator::Function,restrictParams)
if restrictParams == []
	## PDE with no params
	restrictParams2(mesh_fine,mesh_coarse,param_fine,level) = [];
	getOperator2(Mesh,param) = getOperator(Mesh);
	return multilevelOperatorConstructor([],getOperator2,restrictParams2);
else
	return multilevelOperatorConstructor(param,getOperator,restrictParams);
end
end
"""
type Multigrid.CYCLEmem

Fields:

	b::ArrayTypes - memory for the right-hand-side
	r::ArrayTypes - memory for the residual
	x::ArrayTypes - memory for the iterated solution
"""
type CYCLEmem
	b 					::ArrayTypes
	r					::ArrayTypes
	x					::ArrayTypes
end

"""
type Multigrid.MGparam

Fields:

	levels::Int64            - Maximum number of multigrid levels
	numCores::Int64          - Number of OMP cores to work with. Some operations (setup) are not parallelized.
	maxOuterIter::Int64      - Maximum outer iterations.
	relativeTol::Float64	 - Relative L2/Frobenius norm stopping criterion.
	relaxType:: String	 - Relax type. Can be "Jac", "Jac-GMRES" or "SPAI".
	relaxParam::Float64	     - Relax damping parameter.
	relaxPre::Function 	     - pre and post relaxation numbers
	relaxPost::Function	     - Can be 'V', 'F', 'W', 'K' (Krylov cycles are done with FGMRES).
	Ps::Array{SparseCSCTypes} - all matrices here are transposed/conjugated so that parallel multiplication is efficient
	Rs::Array{SparseCSCTypes} - all matrices here are transposed/conjugated so that parallel multiplication is efficient
	As::Array{SparseCSCTypes} - all matrices here are transposed/conjugated so that parallel multiplication is efficient
	relaxPrecs -  an array of relaxation preconditioners for all levels.
	memCycle::Array{CYCLEmem} - Space for x,b and r for each level.
	memRelax::Union{Array{FGMRESmem},Array{BlockFGMRESmem}}  - This is used just in case of GMRES relaxation.
	memKcycle::Union{Array{FGMRESmem},Array{BlockFGMRESmem}} - Memory for the Krylov-cycle FGMRES. First field is ignored.
	coarseSolveType::String - Can be "MUMPS" or "NoMUMPS" for Julia backslash.
	LU - Factorization of coarsest level.
	doTranspose::Int64
	strongConnParam::Float64  - (for SA-AMG only) A threshold for determining a strong connection should >0.25, and <0.85.
	FilteringParam::Float64	  - (for SA-AMG only) A threshold for prolongation filtering >0.0, and <0.2.
	Meshes::Array{RegularMesh} 		 - Array of Regular Meshes for geometric multigrid.
	transferOperatorType:: String	- (for geometric MG only) May be "FullWeighting", "SystemsFacesLinear"
"""
type MGparam
	levels				:: Int64
	numCores			:: Int64
	maxOuterIter		:: Int64
	relativeTol			:: Float64
	relaxType		    :: String
	relaxParam			:: Union{Array{Float64},Float64}
	relaxPre			:: Function
	relaxPost			:: Function
	cycleType			:: Char
	Ps					:: Array{SparseCSCTypes}
	Rs					:: Array{SparseCSCTypes}
	As					:: Array{SparseCSCTypes}
	relaxPrecs
	memCycle			:: Array{CYCLEmem}
	memRelax			:: Union{Array{FGMRESmem},Array{BlockFGMRESmem}}
	memKcycle			:: Union{Array{FGMRESmem},Array{BlockFGMRESmem}}
	coarseSolveType		:: String
	LU
	doTranspose			:: Int64
	strongConnParam		:: Float64
	FilteringParam		:: Float64
	Meshes				:: Array{RegularMesh}
	transferOperatorType:: String
end

include("MGsetup.jl");
include("SA-AMG.jl");
include("MGcycle.jl");
include("SolveFuncs.jl");
include("SAAMGWrapper.jl");
include("Systems.jl");
include("Vanka.jl");
#include("MulticolorDomainDecomposition/MulticolorDomainDecomposition.jl");

"""
function Multigrid.copySolver(MG::MGparam)

copies the solver parameters without the setup and allocated memory.
"""
function copySolver(MG::MGparam)
return getMGparam(MG.levels,MG.numCores,MG.maxOuterIter,MG.relativeTol,MG.relaxType,MG.relaxParam,
					MG.relaxPre,MG.relaxPost,MG.cycleType,MG.coarseSolveType,MG.strongConnParam,MG.FilteringParam,MG.transferOperatorType);
end

function getMGparam(levels::Int64,numCores::Int64,maxIter::Int64,relativeTol:: Float64,relaxType::String,relaxParam,
					relaxPre::Function,relaxPost::Function,cycleType::Char='V',coarseSolveType::String="NoMUMPS",strongConnParam::Float64=0.4,FilteringParam::Float64 = 0.0,transferOperatorType = "FullWeighting")
return MGparam(levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,[],[],[],[],Array{CYCLEmem}(0),
				Array{FGMRESmem}(0),Array{FGMRESmem}(0),coarseSolveType,[],0,strongConnParam,FilteringParam,Array{RegularMesh}(0),transferOperatorType);
end

function getMGparam(levels::Int64=3,numCores::Int64=8,maxIter::Int64=20,relativeTol::Float64=1e-6,relaxType::String="SPAI",relaxParam=1.0,
					relaxPre::Int64=2,relaxPost::Int64=2,cycleType::Char='V',coarseSolveType::String="NoMUMPS",strongConnParam::Float64=0.4,FilteringParam::Float64 = 0.0,transferOperatorType = "FullWeighting")
relaxPreFun(x) = relaxPre;
relaxPostFun(x) = relaxPost;
return getMGparam(levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPreFun,relaxPostFun,cycleType,coarseSolveType,strongConnParam,FilteringParam,transferOperatorType);
end

function getCYCLEmem(n::Int64,m::Int64,T::Type,withB::Bool=true)
b = zeros(T,0);
if m==1
	if withB
		b = zeros(T,n);
	end
	return CYCLEmem(b,zeros(T,n),zeros(T,n));
else
	if withB
		b = zeros(T,n,m);
	end
	return CYCLEmem(b,zeros(T,n,m),zeros(T,n,m));
end
end

import jInv.Utils.clear!
function clear!(param::MGparam)
param.Ps = Array{SparseMatrixCSC}(0);
param.Rs = Array{SparseMatrixCSC}(0);
param.As = Array{SparseMatrixCSC}(0);
param.relaxPrecs = [];
param.memCycle = Array{CYCLEmem}(0);
param.memRelax = Array{FGMRESmem}(0);
param.memKcycle = Array{FGMRESmem}(0);
param.Meshes = Array{RegularMesh}(0);
destroyCoarsestLU(param);
end

function destroyCoarsestLU(param::MGparam)
if param.LU==[]
	return;
end
if param.coarseSolveType=="MUMPS"
	if isa(param.LU,MUMPSfactorization)
		destroyMUMPS(param.LU);
		param.LU = [];
	end
else
	param.LU = [];
end
return;
end

function hierarchyExists(param::MGparam)
return length(param.As) > 0;
end

end # Module


text:
"""
    BernLik <: Likelihood

Bernoulli likelihood
```math
p(y = k | f) = θᵏ (1 - θ)^{1-k}
```
for ``k ∈ \\{0,1\\}``, where ``θ = Φ(f)`` and ``f`` is the latent Gaussian process.
"""
struct BernLik <: Likelihood end

function log_dens(bernoulli::BernLik, f::VecF64, y::Vector{Bool})
    return Float64[yi ? log(Φ(fi)) : log(1.0 - Φ(fi)) for (fi,yi) in zip(f,y)]
end

function dlog_dens_df(bernoulli::BernLik, f::VecF64, y::Vector{Bool})
    return Float64[yi ? φ(fi)/Φ(fi) : -φ(fi)/(1.0 - Φ(fi)) for (fi,yi) in zip(f,y)]
end

#mean and variance under the likelihood
mean_lik(bernoulli::BernLik, f::VecF64) = Float64[Φ(fi) for fi in f]
var_lik(bernoulli::BernLik, f::VecF64) = Float64[Φ(fi)*(1-Φ(fi)) for fi in f]

get_params(bernoulli::BernLik) = []
num_params(bernoulli::BernLik) = 0

#Computes the predictive mean and variance
function predict_obs(bernoulli::BernLik, fmean::VecF64, fvar::VecF64)
    p = Float64[Φ(fm./sqrt(1+fv)) for (fm,fv) in zip(fmean,fvar)]
    return p, p-p.*p
end


text:
"""
    read_dump_prim(dump_name::String)

This will read data in `dump_name` file directly, returning a Dict variable, without spliting them into different array.
"""
function read_dump_prim(dump_name::String)
    io = open(dump_name, "r")
    str = read(io, String)

    # Reading all digital data
    reg = r"\b[\d\se\+ -\.]+\n"
    res = findall(reg, str)
    res = [str[n] for n in res]
    num_steps = Int64(length(res) / 4)

    # Reading Simulation Info
    step_vec = parse(Int64, res[1][2:end])
    num_atoms = parse(Int64, res[2][2:end])
    box_info = [str2float(x) for x in split(res[3], "\n")[2:end-1]]
    box_info = Array(hcat(box_info...)')
    box_temp = copy_array(box_info)
    atom_info = [str2float(x) for x in split(res[4], "\n")[2:end-1]]
    atom_info = Array(hcat(atom_info...)')
    atom_temp = copy_array(atom_info)

    box_info = Array{Float64}(undef, (num_steps, size(box_info, 1), size(box_info, 2)))
    box_info[1, :, :] .= box_temp
    atom_info = Array{Float64}(undef, (num_steps, size(atom_info, 1), size(atom_info, 2)))
    atom_info[1, :, :] .= sortslices(atom_temp, dims=1)

    # Rearanging Data
    for step = 2:num_steps
        step_vec = vcat(step_vec, parse(Int64, res[(step-1)*4+1][2:end]))
        box_info[step, :, :] = Array(hcat([str2float(x) for x in split(res[(step-1)*4+3], "\n")[2:end-1]]...)')
        atom_temp = Array(hcat([str2float(x) for x in split(res[(step-1)*4+4], "\n")[2:end-1]]...)')
        atom_info[step, :, :] = sortslices(atom_temp, dims=1)
    end

    # Output
    res = Dict(
        "step_vec" => step_vec,
        "num_atoms" => num_atoms,
        "num_steps" => num_steps,
        "box_info" => box_info,
        "atom_info" => atom_info,
    )
end

"""
    split_info!(data::Dict; id=false, atom_type=false, mol=false, element=false, charge=false, mass=false, coord=false, coord_scl=false, vel=false, force=false, acc=false)

This will split data readed from dump file, adding variable named same as `kwg` before to `data`.
The value of `kwg` can be `false` or column range of property of interest.

# Example

```julia-repl
dump_name = "xxxx" # name of dump file
data = read_dump_prim(dump_name)
split_info!(data, id=1, atom_type=2, coord=3:5)
```
Code above will add three elements, named as "id", "atom_type", and "coord", to data which contains variable in column 1, 2, 3:5 of `data["atom_info"]`
"""
function split_info!(data::Dict; id=false, atom_type=false, mol=false, element=false, charge=false, mass=false, coord=false, coord_scl=false, vel=false, force=false, acc=false)
    atom_info = data["atom_info"]
    if id != false
        data["id"] = convert.(Int64, atom_info[1, :, id])
        println("\"id\"\t\thas been added to `data` ")
    end
    if atom_type != false
        data["atom_type"] = convert.(Int64, atom_info[1, :, atom_type])
        println("\"atom_type\"\thas been added to `data` ")
    end
    if mol != false
        data["mol"] = convert.(Int64, atom_info[1, :, mol])
        println("\"mol\"\thas been added to `data`")
    end
    if element != false
        data["element"] = convert.(Int64, atom_info[1, :, element])
        println("\"element\"\thas been added to `data`")
    end
    if charge != false
        data["charge"] = atom_info[1, :, charge]
        println("\"charge\"\thas been added to `data`")
    end
    if mass != false
        data["mass"] = atom_info[1, :, mass]
        println("\"mass\"\t\thas been added to `data`")
    end
    if coord != false
        data["coord"] = atom_info[:, :, coord]
        println("\"coord\"\t\thas been added to `data`")
    end
    if coord_scl != false
        data["coord_scl"] = atom_info[:, :, coord_scl]
        println("\"coord_scl\"\thas been added to `data` ")
    end
    if vel != false
        data["vel"] = atom_info[:, :, vel]
        println("\"vel\"\t\thas been added to `data`")
    end
    if acc != false
        data["acc"] = atom_info[:, :, acc]
        println("\"acc\"\thas been added to `data`")
    end
    if force != false
        data["force"] = atom_info[:, :, force]
        println("\"force\"\t\thas been added to `data`")
    end
end

"""
    read_dump(dump_name::String; id=false, atom_type=false, mol=false, element=false, charge=false, mass=false, coord=false, coord_scl=false, vel=false, force=false, acc=false)

This will read data in `dump_name` file and return a Dict which contains elements
- `step_vec`
- `num_atoms`
- `num_steps`
- `box_info`
- `atom_info`

# `kwg`

Each `kwg` can be `false` or column range of property of interest in `atom_info`. If a non-`false` `kwg` is activated, a element with the same name as `kwg` will be add to Dict that will return.

# Example

```julia-repl
dump_name = "xxxx"  # name of dump file
data = read_dump(dump_name, id=1, atom_type=2, coord=3:5)
```
Code above will return a Dict contains elements:
- `step_vec`
- `num_atoms`
- `num_steps`
- `box_info`
- `atom_info`
- `id`
- `atom_type`
- `coord`

Where `id`, `atom_type`, and `coord` are values of `atom_info` in column 1, 2, 3:5 respectively.
"""
function read_dump(dump_name::String; id=false, atom_type=false, mol=false, element=false, charge=false, mass=false, coord=false, coord_scl=false, vel=false, force=false, acc=false)
    data = read_dump_prim(dump_name)
    split_info!(data, id=id, atom_type=atom_type, mol=mol, element=element, charge=charge, mass=mass,
        coord=coord, coord_scl=coord_scl, vel=vel, acc=acc, force=force)
    return data
end

"""
function split_atom(data::Dict, atom_type::Int64)
This will return a new Dict contain elements below:
- "step_vec"
- "num_atoms"
- "num_steps"
- "box_info"
- "atom_info"
# Notice
- 1) "atom_type" should be included in `data`, which is created by `read_dump`
- 2) `split_info!` should be called to get details of "atom_info"
"""
function split_atom(data::Dict, atom_type::Int64)
    try
        data["atom_type"]
    catch
        error("\"atom_type\" has not been included in `data` !")
    end
    pos = findall(x->x==atom_type, data["atom_type"])
    res = Dict()
    # Output
    res = Dict(
        "step_vec" => data["step_vec"],
        "num_atoms" => length(pos),
        "num_steps" => data["num_steps"],
        "box_info" => data["box_info"],
        "atom_info" => data["atom_info"][:, pos, :],
    )
    @printf("\nInformation of atom with \"atom_type=%d\" has been created!\n", atom_type)
    return res
end

"""
    function momentum!(data::Dict)

This will add a component called "momentum" into `data` while `vel` and `mass` is contained in `data`
"""
function momentum!(data::Dict)
    try
        data["mass"]
    catch
        error("Info of \"mass\" is not contained in `data`")
    end
    try
        data["vel"]
    catch
        error("Info of \"vel\" is not contained in `data`")
    end

    momentum = zeros(size(data["vel"]))
    mass_diag = diag(data["mass"])
    for dim = 1 : 3
        momentum[:, :, dim] .= data["vel"][:, :, dim] * mass_diag
    end
    data["momentum"] = momentum
    println("\"momentum\"\thas been added to `data`")
end

"""
    function mass!(data::Dict, mass_vec, id_vec=false)

This will add a component called "mass" into `data`. Elements in `mass_vec` are the mass of atoms with id correspond to `id_vec`. If `id_vec` is set to default, it means `id_vec = [1, 2, 3, ... length(mass_vec)]'
"""
function mass!(data::Dict, mass_vec, id_vec=false)
    if id_vec == false
        id_vec = [i for i = 1:length(mass_vec)]
    end
    id = copy_array(data["atom_type"][1, :])
    mass = zeros(size(id))
    for i = 1:length(mass_vec)
        mass[findall(x->x==id_vec[i], id)] .= mass_vec[i]
    end
    data["mass"] = mass
    println("\"mass\"\thas been added to `data`")
end

"""
    function time_step!(data::Dict, time_len, converter, dump_len)

This will add a component, which is an array contain all information needed to calulate the time in MD_simulation, called "time_step" into `data`.

# Parameters:
- `data`: Dict variable create by `read_dump()`.
- `time_len`: length of time step of simulation.
- `converter`: unit converter that convert time_len to second.
- `dump_len`: length of dump step.
"""
function time_step!(data::Dict, time_len, converter)
    data["time_step"] = [time_len, converter, data["num_steps"]]
    println("\"time_step\"\thas been added to `data`")
end

"""
    function coord_scl!(data::Dict)

This will add a component, which is the scaled coordinate within range [0, 1], called "coord_scl" into `data`.
"""
function coord_scl!(data::Dict)
    try
        data["coord"]
    catch
        error("\"coord\" has not been included in `data` !")
    end
    box_inv = genr_box_inv(data)
    coord_scl = zeros(size(data["coord"]))

    for step = 1:data["num_steps"]
        coord_scl[step, :, :] = data["coord"][step, :, :] * box_inv
    end

    data["coord_scl"] = coord_scl
    println("\"coord_scl\"\thas been added to `data`")
end


text:
int_rules_4_3_4_2 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.3.4.2*(a+b*tan)^m*(c+d*tan)^n*(A+B*tan+C*tan^2) =#
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~(m') * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~(n') * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x) + ~(C') * tan(~(e') + ~(f') * ~x) ^ 2), ~x) => (1 / (~b) ^ 2) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ (~m + 1) * (~c + ~d * tan(~e + ~f * ~x)) ^ ~n * ((~b * ~B - ~a * ~C) + ~b * ~C * tan(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~C, ~m, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && EqQ((~A * (~b) ^ 2 - ~a * ~b * ~B) + (~a) ^ 2 * ~C, 0))
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~(m') * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~(n') * (~(A') + ~(C') * tan(~(e') + ~(f') * ~x) ^ 2), ~x) => (-(~C) / (~b) ^ 2) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ (~m + 1) * (~c + ~d * tan(~e + ~f * ~x)) ^ ~n * (~a - ~b * tan(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~C, ~m, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && EqQ(~A * (~b) ^ 2 + (~a) ^ 2 * ~C, 0))
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~(m') * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~(n') * (~A + ~(C') * tan(~(e') + ~(f') * ~x) ^ 2), ~x) => (~A / ~f) * Subst(Antiderivative((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~n, ~x), ~x, tan(~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~C, ~m, ~n], ~x) && EqQ(~A, ~C)
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x) + ~(C') * tan(~(e') + ~(f') * ~x) ^ 2), ~x) => (-((~b * ~c - ~a * ~d)) * (((~c) ^ 2 * ~C - ~B * ~c * ~d) + ~A * (~d) ^ 2) * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1)) / ((~d) ^ 2 * ~f * (~n + 1) * ((~c) ^ 2 + (~d) ^ 2)) + (1 / (~d * ((~c) ^ 2 + (~d) ^ 2))) * Antiderivative((~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1) * Simp(~a * ~d * ((~A * ~c - ~c * ~C) + ~B * ~d) + ~b * (((~c) ^ 2 * ~C - ~B * ~c * ~d) + ~A * (~d) ^ 2) + ~d * ((((~A * ~b * ~c + ~a * ~B * ~c) - ~b * ~c * ~C) - ~a * ~A * ~d) + ~b * ~B * ~d + ~a * ~C * ~d) * tan(~e + ~f * ~x) + ~b * ~C * ((~c) ^ 2 + (~d) ^ 2) * tan(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~C], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~c) ^ 2 + (~d) ^ 2, 0) && LtQ(~n, -1)))
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(C') * tan(~(e') + ~(f') * ~x) ^ 2), ~x) => (-((~b * ~c - ~a * ~d)) * ((~c) ^ 2 * ~C + ~A * (~d) ^ 2) * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1)) / ((~d) ^ 2 * ~f * (~n + 1) * ((~c) ^ 2 + (~d) ^ 2)) + (1 / (~d * ((~c) ^ 2 + (~d) ^ 2))) * Antiderivative((~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1) * Simp(~a * ~d * (~A * ~c - ~c * ~C) + ~b * ((~c) ^ 2 * ~C + ~A * (~d) ^ 2) + ~d * (((~A * ~b * ~c - ~b * ~c * ~C) - ~a * ~A * ~d) + ~a * ~C * ~d) * tan(~e + ~f * ~x) + ~b * ~C * ((~c) ^ 2 + (~d) ^ 2) * tan(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~C], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~c) ^ 2 + (~d) ^ 2, 0) && LtQ(~n, -1)))
    @apply_utils Antiderivative((~a + ~(b') * tan(~(e') + ~(f') * ~x)) * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~(n') * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x) + ~(C') * tan(~(e') + ~(f') * ~x) ^ 2), ~x) => (~b * ~C * tan(~e + ~f * ~x) * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~d * ~f * (~n + 2)) - (1 / (~d * (~n + 2))) * Antiderivative((~c + ~d * tan(~e + ~f * ~x)) ^ ~n * Simp(((~b * ~c * ~C - ~a * ~A * ~d * (~n + 2)) - ((~A * ~b + ~a * ~B) - ~b * ~C) * ~d * (~n + 2) * tan(~e + ~f * ~x)) - (~a * ~C * ~d * (~n + 2) - ~b * (~c * ~C - ~B * ~d * (~n + 2))) * tan(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~C, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~c) ^ 2 + (~d) ^ 2, 0) && Not(LtQ(~n, -1))))
    @apply_utils Antiderivative((~a + ~(b') * tan(~(e') + ~(f') * ~x)) * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~(n') * (~(A') + ~(C') * tan(~(e') + ~(f') * ~x) ^ 2), ~x) => (~b * ~C * tan(~e + ~f * ~x) * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~d * ~f * (~n + 2)) - (1 / (~d * (~n + 2))) * Antiderivative((~c + ~d * tan(~e + ~f * ~x)) ^ ~n * Simp(((~b * ~c * ~C - ~a * ~A * ~d * (~n + 2)) - (~A * ~b - ~b * ~C) * ~d * (~n + 2) * tan(~e + ~f * ~x)) - (~a * ~C * ~d * (~n + 2) - ~b * ~c * ~C) * tan(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~C, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~c) ^ 2 + (~d) ^ 2, 0) && Not(LtQ(~n, -1))))
    @apply_utils Antiderivative((~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~(n') * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x) + ~(C') * tan(~(e') + ~(f') * ~x) ^ 2), ~x) => (((~a * ~A + ~b * ~B) - ~a * ~C) * (~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1)) / (2 * ~f * ~m * (~b * ~c - ~a * ~d)) + (1 / (2 * ~a * ~m * (~b * ~c - ~a * ~d))) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ (~m + 1) * (~c + ~d * tan(~e + ~f * ~x)) ^ ~n * Simp(~b * (~c * (~A + ~C) * ~m - ~B * ~d * (~n + 1)) + ~a * ((~B * ~c * ~m + ~C * ~d * (~n + 1)) - ~A * ~d * (2 * ~m + ~n + 1)) + (~b * ~C * ~d * ((~m - ~n) - 1) + ~A * ~b * ~d * (~m + ~n + 1) + ~a * (2 * ~c * ~C * ~m - ~B * ~d * (~m + ~n + 1))) * tan(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~C, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && (LtQ(~m, 0) || EqQ(~m + ~n + 1, 0))))
    @apply_utils Antiderivative((~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~(n') * (~(A') + ~(C') * tan(~(e') + ~(f') * ~x) ^ 2), ~x) => (~a * (~A - ~C) * (~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1)) / (2 * ~f * ~m * (~b * ~c - ~a * ~d)) + (1 / (2 * ~a * ~m * (~b * ~c - ~a * ~d))) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ (~m + 1) * (~c + ~d * tan(~e + ~f * ~x)) ^ ~n * Simp(~b * ~c * (~A + ~C) * ~m + ~a * (~C * ~d * (~n + 1) - ~A * ~d * (2 * ~m + ~n + 1)) + (~b * ~C * ~d * ((~m - ~n) - 1) + ~A * ~b * ~d * (~m + ~n + 1) + 2 * ~a * ~c * ~C * ~m) * tan(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~C, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && (LtQ(~m, 0) || EqQ(~m + ~n + 1, 0))))
    @apply_utils Antiderivative((~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~(m') * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x) + ~(C') * tan(~(e') + ~(f') * ~x) ^ 2), ~x) => ((((~c) ^ 2 * ~C - ~B * ~c * ~d) + ~A * (~d) ^ 2) * (~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~d * ~f * (~n + 1) * ((~c) ^ 2 + (~d) ^ 2)) - (1 / (~a * ~d * (~n + 1) * ((~c) ^ 2 + (~d) ^ 2))) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1) * Simp((~b * (((~c) ^ 2 * ~C - ~B * ~c * ~d) + ~A * (~d) ^ 2) * ~m - ~a * ~d * (~n + 1) * ((~A * ~c - ~c * ~C) + ~B * ~d)) - ~a * (~d * (~B * ~c - ~A * ~d) * (~m + ~n + 1) - ~C * ((~c) ^ 2 * ~m - (~d) ^ 2 * (~n + 1))) * tan(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~C, ~m], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && (Not(LtQ(~m, 0)) && (LtQ(~n, -1) && NeQ((~c) ^ 2 + (~d) ^ 2, 0)))))
    @apply_utils Antiderivative((~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~(m') * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(C') * tan(~(e') + ~(f') * ~x) ^ 2), ~x) => (((~c) ^ 2 * ~C + ~A * (~d) ^ 2) * (~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~d * ~f * (~n + 1) * ((~c) ^ 2 + (~d) ^ 2)) - (1 / (~a * ~d * (~n + 1) * ((~c) ^ 2 + (~d) ^ 2))) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1) * Simp((~b * ((~c) ^ 2 * ~C + ~A * (~d) ^ 2) * ~m - ~a * ~d * (~n + 1) * (~A * ~c - ~c * ~C)) - ~a * (-(~A) * (~d) ^ 2 * (~m + ~n + 1) - ~C * ((~c) ^ 2 * ~m - (~d) ^ 2 * (~n + 1))) * tan(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~C, ~m], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && (Not(LtQ(~m, 0)) && (LtQ(~n, -1) && NeQ((~c) ^ 2 + (~d) ^ 2, 0)))))
    @apply_utils Antiderivative((~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~(m') * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~(n') * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x) + ~(C') * tan(~(e') + ~(f') * ~x) ^ 2), ~x) => (~C * (~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~d * ~f * (~m + ~n + 1)) + (1 / (~b * ~d * (~m + ~n + 1))) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ ~n * Simp((~A * ~b * ~d * (~m + ~n + 1) + ~C * (~a * ~c * ~m - ~b * ~d * (~n + 1))) - (~C * ~m * (~b * ~c - ~a * ~d) - ~b * ~B * ~d * (~m + ~n + 1)) * tan(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~C, ~m, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && (Not(LtQ(~m, 0)) && NeQ(~m + ~n + 1, 0))))
    @apply_utils Antiderivative((~a + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~(m') * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~(n') * (~(A') + ~(C') * tan(~(e') + ~(f') * ~x) ^ 2), ~x) => (~C * (~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~d * ~f * (~m + ~n + 1)) + (1 / (~b * ~d * (~m + ~n + 1))) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ ~n * Simp((~A * ~b * ~d * (~m + ~n + 1) + ~C * (~a * ~c * ~m - ~b * ~d * (~n + 1))) - ~C * ~m * (~b * ~c - ~a * ~d) * tan(~e + ~f * ~x), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~C, ~m, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ((~a) ^ 2 + (~b) ^ 2, 0) && (Not(LtQ(~m, 0)) && NeQ(~m + ~n + 1, 0))))
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x) + ~(C') * tan(~(e') + ~(f') * ~x) ^ 2), ~x) => ((~A * (~d) ^ 2 + ~c * (~c * ~C - ~B * ~d)) * (~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~d * ~f * (~n + 1) * ((~c) ^ 2 + (~d) ^ 2)) - (1 / (~d * (~n + 1) * ((~c) ^ 2 + (~d) ^ 2))) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ (~m - 1) * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1) * Simp(((~A * ~d * (~b * ~d * ~m - ~a * ~c * (~n + 1)) + (~c * ~C - ~B * ~d) * (~b * ~c * ~m + ~a * ~d * (~n + 1))) - ~d * (~n + 1) * ((~A - ~C) * (~b * ~c - ~a * ~d) + ~B * (~a * ~c + ~b * ~d)) * tan(~e + ~f * ~x)) - ~b * (~d * (~B * ~c - ~A * ~d) * (~m + ~n + 1) - ~C * ((~c) ^ 2 * ~m - (~d) ^ 2 * (~n + 1))) * tan(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~C], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && (NeQ((~c) ^ 2 + (~d) ^ 2, 0) && (GtQ(~m, 0) && LtQ(~n, -1)))))
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(C') * tan(~(e') + ~(f') * ~x) ^ 2), ~x) => ((~A * (~d) ^ 2 + (~c) ^ 2 * ~C) * (~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~d * ~f * (~n + 1) * ((~c) ^ 2 + (~d) ^ 2)) - (1 / (~d * (~n + 1) * ((~c) ^ 2 + (~d) ^ 2))) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ (~m - 1) * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1) * Simp(((~A * ~d * (~b * ~d * ~m - ~a * ~c * (~n + 1)) + ~c * ~C * (~b * ~c * ~m + ~a * ~d * (~n + 1))) - ~d * (~n + 1) * ((~A - ~C) * (~b * ~c - ~a * ~d)) * tan(~e + ~f * ~x)) + ~b * (~A * (~d) ^ 2 * (~m + ~n + 1) + ~C * ((~c) ^ 2 * ~m - (~d) ^ 2 * (~n + 1))) * tan(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~C], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && (NeQ((~c) ^ 2 + (~d) ^ 2, 0) && (GtQ(~m, 0) && LtQ(~n, -1)))))
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~(m') * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x) + ~(C') * tan(~(e') + ~(f') * ~x) ^ 2), ~x) => (~C * (~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~d * ~f * (~m + ~n + 1)) + (1 / (~d * (~m + ~n + 1))) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ (~m - 1) * (~c + ~d * tan(~e + ~f * ~x)) ^ ~n * Simp(((~a * ~A * ~d * (~m + ~n + 1) - ~C * (~b * ~c * ~m + ~a * ~d * (~n + 1))) + ~d * ((~A * ~b + ~a * ~B) - ~b * ~C) * (~m + ~n + 1) * tan(~e + ~f * ~x)) - (~C * ~m * (~b * ~c - ~a * ~d) - ~b * ~B * ~d * (~m + ~n + 1)) * tan(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~C, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && (NeQ((~c) ^ 2 + (~d) ^ 2, 0) && (GtQ(~m, 0) && Not(IGtQ(~n, 0) && (Not(IntegerQ(~m)) || EqQ(~c, 0) && NeQ(~a, 0)))))))
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~(m') * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(C') * tan(~(e') + ~(f') * ~x) ^ 2), ~x) => (~C * (~a + ~b * tan(~e + ~f * ~x)) ^ ~m * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~d * ~f * (~m + ~n + 1)) + (1 / (~d * (~m + ~n + 1))) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ (~m - 1) * (~c + ~d * tan(~e + ~f * ~x)) ^ ~n * Simp(((~a * ~A * ~d * (~m + ~n + 1) - ~C * (~b * ~c * ~m + ~a * ~d * (~n + 1))) + ~d * (~A * ~b - ~b * ~C) * (~m + ~n + 1) * tan(~e + ~f * ~x)) - ~C * ~m * (~b * ~c - ~a * ~d) * tan(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~C, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && (NeQ((~c) ^ 2 + (~d) ^ 2, 0) && (GtQ(~m, 0) && Not(IGtQ(~n, 0) && (Not(IntegerQ(~m)) || EqQ(~c, 0) && NeQ(~a, 0)))))))
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x) + ~(C') * tan(~(e') + ~(f') * ~x) ^ 2), ~x) => ((~A * (~b) ^ 2 - ~a * (~b * ~B - ~a * ~C)) * (~a + ~b * tan(~e + ~f * ~x)) ^ (~m + 1) * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~f * (~m + 1) * (~b * ~c - ~a * ~d) * ((~a) ^ 2 + (~b) ^ 2)) + (1 / ((~m + 1) * (~b * ~c - ~a * ~d) * ((~a) ^ 2 + (~b) ^ 2))) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ (~m + 1) * (~c + ~d * tan(~e + ~f * ~x)) ^ ~n * Simp(((~A * (~a * (~b * ~c - ~a * ~d) * (~m + 1) - (~b) ^ 2 * ~d * (~m + ~n + 2)) + (~b * ~B - ~a * ~C) * (~b * ~c * (~m + 1) + ~a * ~d * (~n + 1))) - (~m + 1) * (~b * ~c - ~a * ~d) * ((~A * ~b - ~a * ~B) - ~b * ~C) * tan(~e + ~f * ~x)) - ~d * (~A * (~b) ^ 2 - ~a * (~b * ~B - ~a * ~C)) * (~m + ~n + 2) * tan(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~C, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && (NeQ((~c) ^ 2 + (~d) ^ 2, 0) && (LtQ(~m, -1) && Not(ILtQ(~n, -1) && (Not(IntegerQ(~m)) || EqQ(~c, 0) && NeQ(~a, 0)))))))
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(C') * tan(~(e') + ~(f') * ~x) ^ 2), ~x) => ((~A * (~b) ^ 2 + (~a) ^ 2 * ~C) * (~a + ~b * tan(~e + ~f * ~x)) ^ (~m + 1) * (~c + ~d * tan(~e + ~f * ~x)) ^ (~n + 1)) / (~f * (~m + 1) * (~b * ~c - ~a * ~d) * ((~a) ^ 2 + (~b) ^ 2)) + (1 / ((~m + 1) * (~b * ~c - ~a * ~d) * ((~a) ^ 2 + (~b) ^ 2))) * Antiderivative((~a + ~b * tan(~e + ~f * ~x)) ^ (~m + 1) * (~c + ~d * tan(~e + ~f * ~x)) ^ ~n * Simp(((~A * (~a * (~b * ~c - ~a * ~d) * (~m + 1) - (~b) ^ 2 * ~d * (~m + ~n + 2)) - ~a * ~C * (~b * ~c * (~m + 1) + ~a * ~d * (~n + 1))) - (~m + 1) * (~b * ~c - ~a * ~d) * (~A * ~b - ~b * ~C) * tan(~e + ~f * ~x)) - ~d * (~A * (~b) ^ 2 + (~a) ^ 2 * ~C) * (~m + ~n + 2) * tan(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~C, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && (NeQ((~c) ^ 2 + (~d) ^ 2, 0) && (LtQ(~m, -1) && Not(ILtQ(~n, -1) && (Not(IntegerQ(~m)) || EqQ(~c, 0) && NeQ(~a, 0)))))))
    @apply_utils Antiderivative((~(A') + ~(B') * tan(~(e') + ~(f') * ~x) + ~(C') * tan(~(e') + ~(f') * ~x) ^ 2) / ((~a + ~(b') * tan(~(e') + ~(f') * ~x)) * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x))), ~x) => (((~a * ((~A * ~c - ~c * ~C) + ~B * ~d) + ~b * ((~B * ~c - ~A * ~d) + ~C * ~d)) * ~x) / (((~a) ^ 2 + (~b) ^ 2) * ((~c) ^ 2 + (~d) ^ 2)) + (((~A * (~b) ^ 2 - ~a * ~b * ~B) + (~a) ^ 2 * ~C) / ((~b * ~c - ~a * ~d) * ((~a) ^ 2 + (~b) ^ 2))) * Antiderivative((~b - ~a * tan(~e + ~f * ~x)) / (~a + ~b * tan(~e + ~f * ~x)), ~x)) - ((((~c) ^ 2 * ~C - ~B * ~c * ~d) + ~A * (~d) ^ 2) / ((~b * ~c - ~a * ~d) * ((~c) ^ 2 + (~d) ^ 2))) * Antiderivative((~d - ~c * tan(~e + ~f * ~x)) / (~c + ~d * tan(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~C], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && NeQ((~c) ^ 2 + (~d) ^ 2, 0)))
    @apply_utils Antiderivative((~(A') + ~(C') * tan(~(e') + ~(f') * ~x) ^ 2) / ((~a + ~(b') * tan(~(e') + ~(f') * ~x)) * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x))), ~x) => (((~a * (~A * ~c - ~c * ~C) - ~b * (~A * ~d - ~C * ~d)) * ~x) / (((~a) ^ 2 + (~b) ^ 2) * ((~c) ^ 2 + (~d) ^ 2)) + ((~A * (~b) ^ 2 + (~a) ^ 2 * ~C) / ((~b * ~c - ~a * ~d) * ((~a) ^ 2 + (~b) ^ 2))) * Antiderivative((~b - ~a * tan(~e + ~f * ~x)) / (~a + ~b * tan(~e + ~f * ~x)), ~x)) - (((~c) ^ 2 * ~C + ~A * (~d) ^ 2) / ((~b * ~c - ~a * ~d) * ((~c) ^ 2 + (~d) ^ 2))) * Antiderivative((~d - ~c * tan(~e + ~f * ~x)) / (~c + ~d * tan(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~C], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && NeQ((~c) ^ 2 + (~d) ^ 2, 0)))
    @apply_utils Antiderivative(((~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x) + ~(C') * tan(~(e') + ~(f') * ~x) ^ 2)) / (~(a') + ~(b') * tan(~(e') + ~(f') * ~x)), ~x) => (1 / ((~a) ^ 2 + (~b) ^ 2)) * Antiderivative((~c + ~d * tan(~e + ~f * ~x)) ^ ~n * Simp(~b * ~B + ~a * (~A - ~C) + (~a * ~B - ~b * (~A - ~C)) * tan(~e + ~f * ~x), ~x), ~x) + (((~A * (~b) ^ 2 - ~a * ~b * ~B) + (~a) ^ 2 * ~C) / ((~a) ^ 2 + (~b) ^ 2)) * Antiderivative(((~c + ~d * tan(~e + ~f * ~x)) ^ ~n * (1 + tan(~e + ~f * ~x) ^ 2)) / (~a + ~b * tan(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~C, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && (NeQ((~c) ^ 2 + (~d) ^ 2, 0) && (Not(GtQ(~n, 0)) && Not(LeQ(~n, -1))))))
    @apply_utils Antiderivative(((~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(C') * tan(~(e') + ~(f') * ~x) ^ 2)) / (~(a') + ~(b') * tan(~(e') + ~(f') * ~x)), ~x) => (1 / ((~a) ^ 2 + (~b) ^ 2)) * Antiderivative((~c + ~d * tan(~e + ~f * ~x)) ^ ~n * Simp(~a * (~A - ~C) - (~A * ~b - ~b * ~C) * tan(~e + ~f * ~x), ~x), ~x) + ((~A * (~b) ^ 2 + (~a) ^ 2 * ~C) / ((~a) ^ 2 + (~b) ^ 2)) * Antiderivative(((~c + ~d * tan(~e + ~f * ~x)) ^ ~n * (1 + tan(~e + ~f * ~x) ^ 2)) / (~a + ~b * tan(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~C, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && (NeQ((~c) ^ 2 + (~d) ^ 2, 0) && (Not(GtQ(~n, 0)) && Not(LeQ(~n, -1))))))
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(B') * tan(~(e') + ~(f') * ~x) + ~(C') * tan(~(e') + ~(f') * ~x) ^ 2), ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], (ff / ~f) * Subst(Antiderivative(((~a + ~b * ff * ~x) ^ ~m * (~c + ~d * ff * ~x) ^ ~n * (~A + ~B * ff * ~x + ~C * ff ^ 2 * (~x) ^ 2)) / (1 + ff ^ 2 * (~x) ^ 2), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~B, ~C, ~m, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && NeQ((~c) ^ 2 + (~d) ^ 2, 0)))
    @apply_utils Antiderivative((~(a') + ~(b') * tan(~(e') + ~(f') * ~x)) ^ ~m * (~(c') + ~(d') * tan(~(e') + ~(f') * ~x)) ^ ~n * (~(A') + ~(C') * tan(~(e') + ~(f') * ~x) ^ 2), ~x) => With([ff = FreeFactors(tan(~e + ~f * ~x), ~x)], (ff / ~f) * Subst(Antiderivative(((~a + ~b * ff * ~x) ^ ~m * (~c + ~d * ff * ~x) ^ ~n * (~A + ~C * ff ^ 2 * (~x) ^ 2)) / (1 + ff ^ 2 * (~x) ^ 2), ~x), ~x, tan(~e + ~f * ~x) / ff)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~A, ~C, ~m, ~n], ~x) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ((~a) ^ 2 + (~b) ^ 2, 0) && NeQ((~c) ^ 2 + (~d) ^ 2, 0)))
end

max_stars_repo_path: src/rules/4 Trig functions/4.3 Tangent/4.3.4.2 (a+b tan)^m (c+d tan)^n (A+B tan+C tan^2).jl


text:
push!(LOAD_PATH,"../src/")
using Documenter, ARCHModels
DocMeta.setdocmeta!(ARCHModels, :DocTestSetup, :(using ARCHModels; using Random; Random.seed!(1)); recursive=true)
makedocs(modules=[ARCHModels],
        sitename="ARCHModels.jl",
        format = Documenter.HTML(assets=["assets/invenia.css"]),
        doctest=true,
        strict=true,
        pages = ["Home" => "index.md",
                 "introduction.md",
                 "Type Hierarchy" => Any[
                        "univariatetypehierarchy.md",
                        "multivariatetypehierarchy.md"
                        ],
                 "usage.md",
                 "reference.md"
                 ]
        )

deploydocs(repo="github.com/s-broda/ARCHModels.jl.git")


text:
# For the ambiguity test, we need to check the dependencies first
using Interpolations, OffsetArrays
using Test, Random
using OffsetArrays: IdentityUnitRange

if !isdefined(@__MODULE__, :ambs)
    const ambs = detect_ambiguities(Base, Interpolations, OffsetArrays)
end

using CenterIndexedArrays
using CenterIndexedArrays: SymRange

@testset "Ambiguities" begin
    ambscia = detect_ambiguities(Base, Interpolations, OffsetArrays, CenterIndexedArrays)
    VERSION >= v"1.1" && @test isempty(setdiff(ambscia, ambs))
end

@testset "SymRange" begin
    r = SymRange(3)
    @test first(r) == -3
    @test last(r)  ==  3
    @test axes(r) == (r,)
    @test axes(r, 1) == r
    @test size(r) == (7,)

    items = []
    for i in r
        push!(items, i)
    end
    @test items == -3:3

    @test r[-3] == -3
    @test r[3] == 3
    @test_throws BoundsError r[-4]
    @test_throws BoundsError r[4]
    s = SymRange(2)
    @test r[s] == s
    @test r[1:2] == 1:2
    @test intersect(r, s) == s

    @test convert(SymRange, -2:2) === SymRange(2)
    @test_throws ErrorException convert(SymRange, -2:1)
    # The following seems a bit wrong
    @test @inferred(promote(1:2, SymRange(2))) === (1:2, -2:2)

    io = IOBuffer()
    print(io, r)
    @test String(take!(io)) == "SymRange(3)"
end

@testset "Uninitialized" begin
    @test isa(CenterIndexedArray{Float32,2}(undef, 3, 5), CenterIndexedArray)
    @test isa(CenterIndexedArray{Float32,2}(undef, (3, 5)), CenterIndexedArray)
    @test isa(CenterIndexedArray{Float32}(undef, 3, 5), CenterIndexedArray)
    @test isa(CenterIndexedArray{Float32}(undef, (3, 5)), CenterIndexedArray)
    @test_throws ErrorException CenterIndexedArray{Float32}(undef, 4, 5)
end

@testset "Construction & traits" begin
    dat = rand(3,5)
    A = CenterIndexedArray(dat)
    @test size(A) == size(dat)
    @test axes(A) === (SymRange(1), SymRange(2))
    @test length(A) == length(dat)
    @test ndims(A) == 2
    @test eltype(A) == eltype(dat)
end

@testset "Indexing & iteration" begin
    dat = rand(3,5)
    A = CenterIndexedArray(dat)
    @test A[0,0] == dat[2,3]
    k = 0
    for j = -2:2, i = -1:1
        @test @inferred(A[i,j]) == dat[k+=1]
    end
    @test_throws BoundsError A[3,5]
    @test @inferred(A[:,:]) == A
    @test @inferred(A[:,SymRange(1)]) == CenterIndexedArray(dat[:,2:4])
    @test @inferred(A[SymRange(1),:]) == A
    @test @inferred(A[:,-2:0]) == OffsetArray(dat[:,1:3], -1:1, 1:3)  # axes-of-the-axes
    @test @inferred(A[:,IdentityUnitRange(-2:0)]) == OffsetArray(dat[:,1:3], -1:1, -2:0)
    k = 0
    for j = -2:2, i = -1:1
        A[i,j] = (k+=1)
    end
    @test dat == reshape(1:15, 3, 5)
    @test_throws BoundsError A[3,5] = 15

    rand!(dat)
    iall = (-1:1).*ones(Int, 5)'
    jall = ones(Int,3).*(-2:2)'
    k = 0
    for I in eachindex(A)
        k += 1
        @test I[1] == iall[k]
        @test I[2] == jall[k]
    end

    # Iteration
    for (a, d) in zip(A, dat)
        @test a == d
    end
end

@testset "Display" begin
    A = CenterIndexedArray(reshape(1:15, 3, 5))
    io = IOBuffer()
    show(io, MIME("text/plain"), A)
    str = String(take!(io))
    @test endswith(str, "CenterIndexedArray(reshape(::UnitRange{$Int}, 3, 5)) with eltype $Int with indices SymRange(1)×SymRange(2):\n 1  4  7  10  13\n 2  5  8  11  14\n 3  6  9  12  15")
end

@testset "Operations" begin
    dat = rand(3,5)
    A = CenterIndexedArray(dat)
    # Standard julia operations
    B = copy(A)

    @test B.data == dat
    @test B == A
    @test isequal(B, A)

    @test vec(A) == vec(dat)

    @test minimum(A) == minimum(dat)
    @test maximum(A) == maximum(dat)
    @test minimum(A,dims=1) == CenterIndexedArray(minimum(dat,dims=1))
    @test maximum(A,dims=2) == CenterIndexedArray(maximum(dat,dims=2))

    amin, iamin = findmin(A)
    dmin, idmin = findmin(dat)
    @test amin == dmin
    @test A[iamin] == amin
    @test amin == dat[idmin]

    amax, iamax = findmax(A)
    dmax, idmax = findmax(dat)
    @test amax == dmax
    @test A[iamax] == amax
    @test amax == dat[idmax]

    fill!(A, 2)
    @test all(x->x==2, A)

    ii, jj = begin
        II = findall(!iszero, A)
        (getindex.(II, 1), getindex.(II, 2))
    end
    iall = (-1:1).*ones(Int, 5)'
    jall = ones(Int,3).*(-2:2)'
    @test vec(ii) == vec(iall)
    @test vec(jj) == vec(jall)

    rand!(dat)

    # @test cat(1, A, dat) == cat(1, dat, dat)
    # @test cat(2, A, dat) == cat(2, dat, dat)

    @test permutedims(A, (2,1)) == CenterIndexedArray(permutedims(dat, (2,1)))
    # @test ipermutedims(A, (2,1)) == CenterIndexedArray(ipermutedims(dat, (2,1)))

    @test cumsum(A, dims=1) == CenterIndexedArray(cumsum(dat, dims=1))
    @test cumsum(A, dims=2) == CenterIndexedArray(cumsum(dat, dims=2))

    @test mapslices(v->sort(v), A, dims=1) == CenterIndexedArray(mapslices(v->sort(v), dat, dims=1))
    @test mapslices(v->sort(v), A, dims=2) == CenterIndexedArray(mapslices(v->sort(v), dat, dims=2))

    @test reverse(A, dims=1) == CenterIndexedArray(reverse(dat, dims=1))
    @test reverse(A, dims=2) == CenterIndexedArray(reverse(dat, dims=2))

    @test A .+ 1 == CenterIndexedArray(dat .+ 1)
    @test 2*A == CenterIndexedArray(2*dat)
    @test A+A == CenterIndexedArray(dat+dat)
    @test isa(A .+ 1, CenterIndexedArray)
end


text:
# # Thermal Generation Dispatch Example

#md # [![](https://img.shields.io/badge/GitHub-100000?style=for-the-badge&logo=github&logoColor=white)](@__REPO_ROOT_URL__/examples/Thermal_Generation_Dispatch_Example.jl)

# This example illustrates the sensitivity analysis of thermal generation dispatch problem.

# This problem can be described as the choice of thermal generation `g` given a demand `d`, a price for thermal generation `c` and a penalty price `c_{ϕ}` for any demand not attended ϕ.

# ```math
# \begin{split}
# \begin{array} {ll}
# \mbox{minimize} & \sum_{i=1}^{N} c_{i} g_{i} + c_{\phi} \phi \\
# \mbox{s.t.} & g_{i} \ge 0 \quad i=1..N  \\
#             & g_{i} \le G_{i} \quad i=1..N  \\
#             & \sum_{i=1}^{N} g_{i} + \phi = d\\
# \end{array}
# \end{split}
# ```
# where
# - `G_{i}` is the maximum possible generation for a thermal generator `i`

# ## Define and solve the Thermal Dispatch Problem

# First, import the libraries.

using Test
using JuMP
import DiffOpt
import LinearAlgebra: dot
import HiGHS
import MathOptInterface
import Plots
const MOI = MathOptInterface

# Define the model that will be construct given a set of parameters.

function generate_model(d::Float64; g_sup::Vector{Float64}, c_g::Vector{Float64}, c_ϕ::Float64)
    ## Creation of the Model and Parameters
    model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer))
    set_silent(model)
    I = length(g_sup)

    ## Variables
    @variable(model, g[i in  1:I] >= 0.0)
    @variable(model, ϕ >= 0.0)

    ## Constraints
    @constraint(model, limit_constraints_sup[i in 1:I], g[i] <= g_sup[i])
    @constraint(model, demand_constraint, sum(g) + ϕ == d)

    ## Objectives
    @objective(model, Min, dot(c_g, g) + c_ϕ * ϕ)

    ## Solve the model
    optimize!(model)

    ## Return the solved model
    return model
end

# Define the functions that will get the primal values `g` and `\phi` and sensitivity analysis of the demand `dg/dd` and `d\phi/dd` from a optimized model.

function diff_forward(model::Model, ϵ::Float64 = 1.0)
    ## Initialization of parameters and references to simplify the notation
    vect_ref = [model[:g]; model[:ϕ]]
    I = length(model[:g])

    ## Get the primal solution of the model
    vect = MOI.get.(model, MOI.VariablePrimal(), vect_ref)

    ## Pass the perturbation to the DiffOpt Framework and set the context to Forward
    constraint_equation = convert(MOI.ScalarAffineFunction{Float64}, ϵ)
    MOI.set(model, DiffOpt.ForwardConstraintFunction(), model[:demand_constraint], constraint_equation)
    DiffOpt.forward_differentiate!(model)

    ## Get the derivative of the model
    dvect = MOI.get.(model, DiffOpt.ForwardVariablePrimal(), vect_ref)

    ## Return the values as a vector
    return [vect; dvect]
end

function diff_reverse(model::Model, ϵ::Float64 = 1.0)
    ## Initialization of parameters and references to simplify the notation
    vect_ref = [model[:g]; model[:ϕ]]
    I = length(model[:g])

    ## Get the primal solution of the model
    vect = MOI.get.(model, MOI.VariablePrimal(), vect_ref)

    ## Set variables needed for the DiffOpt Backward Framework
    dvect = Array{Float64, 1}(undef, I + 1)
    perturbation = zeros(I + 1)

    ## Loop for each primal variable
    for i in 1:I+1
        ## Set the perturbation in the Primal Variables and set the context to Backward
        perturbation[i] = ϵ
        MOI.set.(model, DiffOpt.ReverseVariablePrimal(), vect_ref, perturbation)
        DiffOpt.reverse_differentiate!(model)

        ## Get the value of the derivative of the model
        dvect[i] = JuMP.constant(MOI.get(model, DiffOpt.ReverseConstraintFunction(), model[:demand_constraint]))
        perturbation[i] = 0.0
    end

    ## Return the values as a vector
    return [vect;dvect]
end

# Initialize of Parameters

g_sup = [10.0, 20.0, 30.0]
I = length(g_sup)
d = 0.0:0.1:80
d_size = length(d)
c_g = [1.0, 3.0, 5.0]
c_ϕ = 10.0
;

# Generate models for each demand `d`
models = generate_model.(d; g_sup = g_sup, c_g = c_g, c_ϕ = c_ϕ);

# Get the results of models with the DiffOpt Forward and Backward context

result_forward = diff_forward.(models)
optimize!.(models)
result_reverse = diff_reverse.(models);

# Organization of results to plot
# Initialize data_results that will contain every result
data_results = Array{Float64,3}(undef, 2, d_size, 2*(I+1));

# Populate the data_results array
for k in 1:d_size
    data_results[1,k,:] = result_forward[k]
    data_results[2,k,:] = result_reverse[k]
end

# ## Results with Plot graphs
# ### Results for the forward context
# Result Primal Values:
Plots.plot(d,data_results[1,:,1:I+1],
    title="Generation by Demand",label=["Thermal Generation 1" "Thermal Generation 2" "Thermal Generation 3" "Generation Deficit"],
    xlabel="Demand [unit]",ylabel= "Generation [unit]"
)

# Result Sensitivity Analysis:
Plots.plot(d,data_results[1,:,I+2:2*(I+1)],
    title="Sensitivity of Generation by Demand",label=["T. Gen. 1 Sensitivity" "T. Gen. 2 Sensitivity" "T. Gen. 3 Sensitivity" "Gen. Deficit Sensitivity"],
    xlabel="Demand [unit]",ylabel= "Sensitivity [-]"
)

# ### Results for the reverse context
# Result Primal Values:
Plots.plot(d,data_results[2,:,1:I+1],
    title="Generation by Demand",label=["Thermal Generation 1" "Thermal Generation 2" "Thermal Generation 3" "Generation Deficit"],
    xlabel="Demand [unit]",ylabel= "Generation [unit]"
)

# Result Sensitivity Analysis:
Plots.plot(d,data_results[2,:,I+2:2*(I+1)],
    title="Sensitivity of Generation by Demand",label=["T. Gen. 1 Sensitivity" "T. Gen. 2 Sensitivity" "T. Gen. 3 Sensitivity" "Gen. Deficit Sensitivity"],
    xlabel="Demand [unit]",ylabel= "Sensitivity [-]"
)