text:
"""
  genGaussian(mean_in, cov_in, class_label, n)

Draw from a gaussian distribution
# Arguments
- `mean_in` : means
- `cov_in` : covariances
- `class_label` : class_label
- `n` : number of samples to draw
"""
function genGaussian(mean_in, cov_in, class_label, n)
  nv = Distributions.MvNormal(mean_in, cov_in)
  X = rand(nv, n)
  y = ones(n) .* class_label
  return nv, transpose(X), y
end

"""
  genZafarData(n = 10000, d = pi/4)

Generate synthetic data from Zafar et al., 2017 Fairness Constraints: Mechanisms for Fair Classification
# Arguments
- `n=10000` : number of samples
- `d=pi/4` : discrimination factor
# Returns
- `X` : DataFrame containing features and protected attribute z
- `y` : Binary Target variable
"""
function genZafarData(n = 10000, d = pi/4)
    mu1, sigma1 = [2., 2.],  [[5., 1.] [1., 5.]]
    mu2, sigma2 = [-2.,-2.], [[10., 1.]  [1., 3.]]
    nv1, X1, y1 = genGaussian(mu1, sigma1, 1, Int64(floor(n/2))) # positive class
    nv2, X2, y2 = genGaussian(mu2, sigma2, -1, Int64(floor(n/2))) # negative class
    perm = shuffle(1:(2*Int64(floor(n/2))))
    X = vcat(X1, X2)[perm,:]
    y = vcat(y1, y2)[perm]
    rotation_mult = [[cos(d), -sin(d)] [sin(d), cos(d)]]
    X_aux = X * rotation_mult
    """ Generate the sensitive feature here """
    x_control = [] # this array holds the sensitive feature value
    for i in 1:(2*Int64(floor(n/2)))
        x = X_aux[i,:]
        # probability for each cluster that the point belongs to it
        p1 = Distributions.pdf(nv1, x)
        p2 = Distributions.pdf(nv2, x)
        # normalize the probabilities from 0 to 1
        s = p1+p2
        p1 = p1/s
        p2 = p2/s
        if rand(1)[1] < p1 # the first cluster is the positive class
            append!(x_control, 1.0)
        else
            append!(x_control, 0.0)
        end
    end
    X = DataFrame(X)
    X.z = x_control
    coerce!(X, :z => Multiclass)
    y = categorical(y)
    return X, y
end

"""
  genSubgroupData(n=10000, setting="B00")

Generate synthetic data from Loh et al., 2019 : Subgroup identification for precision medicine: A comparative review of 13 methods
# Arguments
- `n=10000` : number of samples
- `setting="B00"` : Simulation data setting: one of "B00", ..., "B02", "B1", ... , "B8"
For "B00", ..., "B02" there is no "bias" in the data, i.e. group membership has no effect on y.
whereas for  "B1", ... , "B8", there is a direct effect of group membership z on y, usually mediated by one or more features.
# Returns
- `X` : DataFrame containing features and protected attribute z
- `y` : Binary Target variable
"""
function genSubgroupData(n = 10000, setting = "B00")
    xdists = [
      Distributions.Normal(0.,1.), # x1 is N(0,1)
      Distributions.MvNormal([0., 0.], [[1, 0.5] [0.5, 1]]), # x2 and x3 are N(0,1) with cor 0.5
      Distributions.Exponential(1.), # x4 is exp(1)
      Distributions.Bernoulli(0.5), # x5 is Ber(0.5)
      Distributions.Categorical(repeat([.1], 10)), # x6 is Multinom. w. equal
      Distributions.MvNormal([0., 0., 0., 0.], [[1, 0.5, 0.5, 0.5] [0.5, 1, 0.5, 0.5] [0.5, 0.5, 1, 0.5] [0.5, 0.5, 0.5, 1]]
      ) # x7 to 10 are N(0,1) with cor 0.5
    ]
    z = rand(Distributions.Bernoulli(0.5), n)
    X = []
    for d in xdists
      if length(X) > 0
        xn = rand(d, n)
        if size(xn)[1] != n
          xn = transpose(xn)
        end
        X = hcat(X, xn)
      else
        X = rand(d, n)
      end
    end
    y = logit_fun(X, z, setting)
    X = DataFrame(X)
    X.z = z
    coerce!(X, :z => Multiclass)
    y = categorical(y)
    return X, y
end

# Logistic Link function
function log_link(x)
  return exp(x)/(1+exp(x))
end

"""
  logit_fun(X, z, setting)

Compute y from X and z according to a setting provided in Loh et al., 2019: Subgroup identification for precision medicine: A comparative review of 13 methods
# Arguments
- `X` : matrix of features
- `z` : vector of group assignments
- `setting` : Simulation data setting: one of "B00", ..., "B02", "B1", ... , "B8"
"""
function logit_fun(X, z, setting)
  if setting == "B00"
    logit = repeat([0.], length(z))
  elseif setting == "B01"
    logit = 0.5*(X[:,1] + X[:,2])
  elseif setting == "B02"
    logit = 0.5*(X[:,1] + X[:,2]) .^ 2
  elseif setting == "B1"
    logit = 0.5*(X[:,1] + X[:,2] - X[:,5]) + 2. * z .* map(x -> ifelse(isodd(Int64(x)), 1., 0.), X[:,6])
  elseif setting == "B2"
    logit = 0.5*X[:,2] + 2. * z .* map(x -> ifelse(x > 0, 1., 0.), X[:,1])
  elseif setting == "B3"
    logit = 0.3*(X[:,1] + X[:,2]) + 2. * z .* map(x -> ifelse(x > 0, 1., 0.), X[:,1])
  elseif setting == "B4"
    logit = 0.2*(X[:,2] + X[:,3] .- 2.) + 2. * z .* X[:,4]
  elseif setting == "B5"
    logit = 0.2*(X[:,1] + X[:,2] .- 3) + 2. * z .* map(x -> ifelse(x < 0, 1., 0.), X[:,1]) .* map(x -> ifelse(isodd(Int64(x)), 1., 0.), X[:,6])
  elseif setting == "B6"
    logit = 0.5*(X[:,2] .- 1.) + 2. * z .* map(x -> ifelse(abs(x) < 0.8, 1., 0.), X[:,1])
  elseif setting == "B7"
    logit = 0.2*(X[:,2] + X[:,2] .^ 2 .- 6.)  + 2. * z .* map(x -> ifelse(x < 0, 1., 0.), X[:,1])
  elseif setting == "B8"
    logit = 0.5 * X[:,2] + 2. * z .* X[:,5]
  end
  yprob = log_link.(logit)
  # Compute a random y with probability yprob
  return ifelse.(yprob .> rand(1)[1], 1, 0)
end


text:
export RigidBodyMotion, Kinematics, d_dt, rigidbodyvelocity, assign_velocity!, assign_velocity
export Oscillation, OscillationX, OscillationY, OscillationXY, RotationalOscillation,
        PitchHeave, Pitchup, EldredgeRamp, ColoniusRamp

using DocStringExtensions
import ForwardDiff
import Base: +, *, -, >>, <<, show

"""
An abstract type for types that takes in time and returns `(c, ċ, c̈, α, α̇, α̈)`.
"""
abstract type Kinematics end

"""
    RigidBodyMotion

A type to store the body's current kinematics

# Fields

- `c`: current centroid position (relative to initial position)
- `ċ`: current centroid velocity
- `c̈`: current centroid acceleration
- `α`: current angle (relative to initial angle)
- `α̇`: current angular velocity
- `α̈`: current angular acceleration
- `kin`: a [`Kinematics`](@ref) structure

The first six fields are meant as a cache of the current kinematics
while the `kin` field can be used to find the plate kinematics at any time.
"""
mutable struct RigidBodyMotion
    c::ComplexF64
    ċ::ComplexF64
    c̈::ComplexF64
    α::Float64
    α̇::Float64
    α̈::Float64

    kin::Kinematics
end

RigidBodyMotion(ċ, α̇) = RigidBodyMotion(0.0im, complex(ċ), 0.0im, 0.0, float(α̇),
                                          0.0, Constant(ċ, α̇))
RigidBodyMotion(kin::Kinematics) = RigidBodyMotion(kin(0)..., kin)
(m::RigidBodyMotion)(t) = m.kin(t)

function (m::RigidBodyMotion)(t,x̃::Tuple{Real,Real})
  # This expects coordinates in body's own coordinate system
  #
  z̃ = ComplexF64(x̃[1],x̃[2])
  m.c, m.ċ, m.c̈, m.α, m.α̇, m.α̈ = m.kin(t)
  z = exp(im*m.α)*z̃
  return m.c + z, m.ċ + im*m.α̇*z, m.c̈ + (im*m.α̈-m.α̇^2)*z
end

"""
    rigidbodyvelocity(m::RigidBodyMotion,t::Real)

Return the velocity components (as a vector) of a `RigidBodyMotion`
at the given time `t`.
"""
function rigidbodyvelocity(motion::RigidBodyMotion,t::Real)
  _,ċ,_,_,α̇,_ = motion(t)
  return [real(ċ),imag(ċ),α̇]
end

"""
    assign_velocity!(u::AbstractVector{Float64},v::AbstractVector{Float64},
                     x::AbstractVector{Float64},y::AbstractVector{Float64},
                     xc::Real,yc::Real,α::Real,
                     motion,t::Real)

Assign the components of rigid body velocity `u` and `v` (in inertial coordinate system)
at positions described by coordinates `x`, `y` (also in inertial coordinate system) at time `t`,
based on supplied motion `motion` for the body.
"""
function assign_velocity!(u::AbstractVector{Float64},v::AbstractVector{Float64},
                          x::AbstractVector{Float64},y::AbstractVector{Float64},
                          xc::Real,yc::Real,α::Real,m::RigidBodyMotion,t::Real)

  length(u) == length(v) == length(x) == length(y) || error("Inconsistent lengths of vectors")

  _,ċ,_,_,α̇,_ = m(t)
  uc = ċ .+ im*α̇*((x .- xc) .+ im*(y .- yc))
  u .= real.(uc)
  v .= imag.(uc)
  #=
  for i = 1:length(x)
      Δz = (x[i]-xc)+im*(y[i]-yc)
      ċi = ċ + im*α̇*Δz
      u[i] = real(ċi)
      v[i] = imag(ċi)
  end
  =#
  return u, v
end

"""
    assign_velocity(x::AbstractVector{Float64},y::AbstractVector{Float64},
                    xc::Real,yc::Real,α::Real,motion,t::Real)

Return the components of rigid body velocities (in inertial components) at positions
described by coordinates `x`, `y` (also in inertial coordinate system) at time `t`,
based on supplied motion `motion` for the body.
"""
assign_velocity(x::AbstractVector{Float64},y::AbstractVector{Float64},
                xc::Real,yc::Real,α::Real,m::RigidBodyMotion,t::Real) =
                assign_velocity!(similar(x),similar(y),x,y,xc,yc,α,m,t)

function show(io::IO, m::RigidBodyMotion)
    println(io, "Rigid Body Motion:")
    println(io, "  ċ = $(round(m.ċ, digits=2))")
    println(io, "  c̈ = $(round(m.c̈, digits=2))")
    println(io, "  α̇ = $(round(m.α̇, digits=2))")
    println(io, "  α̈ = $(round(m.α̈, digits=2))")
    print(io, "  $(m.kin)")
end

ext: jl


text:
module HighOrderFunctionKoans

#=
Links:
  - https://docs.julialang.org/en/v1/manual/functions/#man-anonymous-functions-1
  - https://docs.julialang.org/en/v1/base/math/#Base.round
  - https://docs.julialang.org/en/v1/base/collections/#Base.map
  - https://docs.julialang.org/en/v1/base/collections/#Base.filter
=#

function round_list(list)
    map(round, list)
end

function apply_fn(f, x)
    f(x, 5)
end

function add_n(list, n)
    f = x -> x + n
    map(f, list)
end

function remove_multiples(list, n)
    f = x -> x % n != 0
    filter(f, list)
end

end


text:
# TODO: assert a naming convention -- ??
"""
Read in power-system parameters from a Matpower, PTI, or JSON file and do some
data checks.
"""
function parsestandardfiles(file::String; kwargs...)

    # function `parse_file` is in pm_io/common.jl
    data = parse_file(file)

    # Check for at least one bus in input file
    if (length(data["bus"]) < 1)
        @error "There are no buses in this file"
    end

    # in pm2ps_parser.jl
    data = pm2ps_dict(data; kwargs...)

    return data

end

function parsestandardfiles(file::String, ts_folder::String; kwargs...)

    # TODO: assert a naming convention
    data = parsestandardfiles(file; kwargs...)

    ts_data = read_data_files(ts_folder; kwargs...)

    # assign_ts_data is in forecast_parser.jl
    data = assign_ts_data(data, ts_data)

    return data
end


text:
using Base.Test, ODBC, DataStreams, DataFrames, WeakRefStrings, CSV, SQLite, Feather

@show ODBC.listdrivers()
@show ODBC.listdsns()

@show ODBC.API.odbc_dm

@show run(`odbcinst -q -d`)

# conn_string = "Driver={MySQL ODBC Driver};SERVER=127.0.0.1;Port=3306;Database=mysql;USER=root;PASSWORD=;Option=3"

db = "ODBC-MySQL"
Username = "root"
Password = ""

run(`uname -a`)
run(`mysqlshow -uroot`)

dsn = ODBC.DSN(db, Username, Password)

# Check some basic queries
dbs = ODBC.query(dsn, "show databases")
ODBC.query(dsn, "use mysql")
data = ODBC.query(dsn, "select table_name from information_schema.tables")

# setup a test database
ODBC.execute!(dsn, "drop database if exists testdb")
ODBC.execute!(dsn, "create database testdb")
ODBC.execute!(dsn, "use testdb")
ODBC.execute!(dsn, "drop table if exists test1")
ODBC.execute!(dsn, "create table test1
                    (test_bigint bigint,
                     test_bit bit,
                     test_decimal decimal,
                     test_int int,
                     test_numeric numeric,
                     test_smallint smallint,
                     test_mediumint mediumint,
                     test_tiny_int tinyint,
                     test_float float,
                     test_real double,
                     test_date date,
                     test_datetime datetime,
                     test_timestamp timestamp,
                     test_time time,
                     test_year year,
                     test_char char(1),
                     test_varchar varchar(16),
                     test_binary binary(2),
                     test_varbinary varbinary(16),
                     test_tinyblob tinyblob,
                     test_blob blob,
                     test_mediumblob mediumblob,
                     test_longblob longblob,
                     test_tinytext tinytext,
                     test_text text,
                     test_mediumtext mediumtext,
                     test_longtext longtext
                    )")
data = ODBC.query(dsn, "select * from information_schema.columns where table_name = 'test1'")
ODBC.execute!(dsn, "insert test1 VALUES
                    (1, -- bigint
                     1, -- bit
                     1.0, -- decimal
                     1, -- int
                     1.0, -- numeric
                     1, -- smallint
                     1, -- mediumint
                     1, -- tinyint
                     1.2, -- float
                     1.2, -- double
                     '2016-01-01', -- date
                     '2016-01-01 01:01:01', -- datetime
                     '2016-01-01 01:01:01', -- timestamp
                     '01:01:01', -- time
                     2016, -- year
                     'A', -- char(1)
                     'hey there sailor', -- varchar
                     cast('12' as binary(2)), -- binary
                     NULL, -- varbinary
                     'hey there abraham', -- tinyblob
                     'hey there bill', -- blob
                     'hey there charlie', -- mediumblob
                     'hey there dan', -- longblob
                     'hey there ephraim', -- tinytext
                     'hey there frank', -- text
                     'hey there george', -- mediumtext
                     'hey there hank' -- longtext
                    )")
source = ODBC.Source(dsn, "select * from test1")
data = Data.stream!(source, DataFrame)
@test size(data) == (1,27)
@test Data.types(data) == map(x->Nullable{x},
[Int64,
 Int8,
 is_windows() ? Float64 : DecFP.Dec64,
 Int32,
 is_windows() ? Float64 : DecFP.Dec64,
 Int16,
 Int32,
 Int8,
 Float32,
 Float64,
 ODBC.API.SQLDate,
 ODBC.API.SQLTimestamp,
 ODBC.API.SQLTimestamp,
 ODBC.API.SQLTime,
 Int16,
 WeakRefString{UInt16},
 WeakRefString{UInt16},
 Array{UInt8,1},
 Array{UInt8,1},
 Array{UInt8,1},
 Array{UInt8,1},
 Array{UInt8,1},
 Array{UInt8,1},
 String,
 String,
 String,
 String])
@test data.columns[1][1] === Nullable(Int64(1))
@test data.columns[2][1] === Nullable(Int8(1))
@test data.columns[3][1] === Nullable(is_windows() ? 1.0 : ODBC.DecFP.Dec64(1))
@test data.columns[4][1] === Nullable(Int32(1))
@test data.columns[5][1] === Nullable(is_windows() ? 1.0 : ODBC.DecFP.Dec64(1))
@test data.columns[6][1] === Nullable(Int16(1))
@test data.columns[7][1] === Nullable(Int32(1))
@test data.columns[8][1] === Nullable(Int8(1))
@test data.columns[9][1] === Nullable(Float32(1.2))
@test data.columns[10][1] === Nullable(Float64(1.2))
@test data.columns[11][1] === Nullable(ODBC.API.SQLDate(2016,1,1))
@test data.columns[12][1] === Nullable(ODBC.API.SQLTimestamp(2016,1,1,1,1,1,0))
@test data.columns[13][1] === Nullable(ODBC.API.SQLTimestamp(2016,1,1,1,1,1,0))
@test data.columns[14][1] === Nullable(ODBC.API.SQLTime(1,1,1))
@test data.columns[15][1] === Nullable(Int16(2016))
@test string(get(data.columns[16][1])) == "A"
@test string(get(data.columns[17][1])) == "hey there sailor"
@test get(data.columns[18][1]) == UInt8[0x31,0x32]
@test isnull(data.columns[19][1])
@test get(data.columns[20][1]) == UInt8[0x68,0x65,0x79,0x20,0x74,0x68,0x65,0x72,0x65,0x20,0x61,0x62,0x72,0x61,0x68,0x61,0x6d]
@test get(data.columns[21][1]) == UInt8[0x68,0x65,0x79,0x20,0x74,0x68,0x65,0x72,0x65,0x20,0x62,0x69,0x6c,0x6c]
@test get(data.columns[22][1]) == UInt8[0x68,0x65,0x79,0x20,0x74,0x68,0x65,0x72,0x65,0x20,0x63,0x68,0x61,0x72,0x6c,0x69,0x65]
@test get(data.columns[23][1]) == UInt8[0x68,0x65,0x79,0x20,0x74,0x68,0x65,0x72,0x65,0x20,0x64,0x61,0x6e]
@test string(get(data.columns[24][1])) == "hey there ephraim"
@test string(get(data.columns[25][1])) == "hey there frank"
@test string(get(data.columns[26][1])) == "hey there george"
@test string(get(data.columns[27][1])) == "hey there hank"

ODBC.execute!(dsn, "insert test1 VALUES
                    (1, -- bigint
                     1, -- bit
                     1.0, -- decimal
                     1, -- int
                     1.0, -- numeric
                     1, -- smallint
                     1, -- mediumint
                     1, -- tinyint
                     1.2, -- float
                     1.2, -- double
                     '2016-01-01', -- date
                     '2016-01-01 01:01:01', -- datetime
                     '2016-01-01 01:01:01', -- timestamp
                     '01:01:01', -- time
                     2016, -- year
                     'A', -- char(1)
                     'hey there sailor', -- varchar
                     cast('12' as binary(2)), -- binary
                     NULL, -- varbinary
                     'hey there abraham', -- tinyblob
                     'hey there bill', -- blob
                     'hey there charlie', -- mediumblob
                     'hey there dan', -- longblob
                     'hey there ephraim', -- tinytext
                     'hey there frank', -- text
                     'hey there george', -- mediumtext
                     'hey there hank' -- longtext
                    )")
data = ODBC.query(dsn, "select * from test1")
@test size(data) == (2,27)

ODBC.execute!(dsn, "drop table if exists test2")
ODBC.execute!(dsn, """
CREATE TABLE test2
(
    ID INT NOT NULL PRIMARY KEY,
    first_name VARCHAR(25),
    last_name VARCHAR(25),
    Salary DECIMAL,
    `hourly rate` real,
    hireDate DATE,
    `last clockin` DATETIME
);""")
randoms = joinpath(dirname(@__FILE__), "randoms.csv")
ODBC.execute!(dsn, "load data local infile '$randoms' into table test2
                    fields terminated by ',' lines terminated by '\n'
                    (id,first_name,last_name,salary,`hourly rate`,hiredate,`last clockin`)")

data = ODBC.query(dsn, "select count(*) from test2")
@test size(data) == (1,1)
@test data.columns[1][1] === Nullable(70000)

df = ODBC.query(dsn, "select * from test2")
@test size(df) == (70000,7)
@test df.columns[1].values == [1:70000...]
@test df.columns[end][1] === Nullable(ODBC.API.SQLTimestamp(2002,1,17,21,32,0,0))

# test exporting test1 to CSV
source = ODBC.Source(dsn, "select * from test1")
csv = CSV.Sink("test1.csv")
Data.stream!(source, csv)
open("test1.csv") do f
    @test readline(f) == "\"test_bigint\",\"test_bit\",\"test_decimal\",\"test_int\",\"test_numeric\",\"test_smallint\",\"test_mediumint\",\"test_tiny_int\",\"test_float\",\"test_real\",\"test_date\",\"test_datetime\",\"test_timestamp\",\"test_time\",\"test_year\",\"test_char\",\"test_varchar\",\"test_binary\",\"test_varbinary\",\"test_tinyblob\",\"test_blob\",\"test_mediumblob\",\"test_longblob\",\"test_tinytext\",\"test_text\",\"test_mediumtext\",\"test_longtext\"\n"
    @test readline(f) == "1,1,+1E+0,1,+1E+0,1,1,1,1.2,1.2,2016-01-01,2016-01-01T01:01:01,2016-01-01T01:01:01,01:01:01,2016,\"A\",\"hey there sailor\",UInt8[0x31,0x32],\"\",UInt8[0x68,0x65,0x79,0x20,0x74,0x68,0x65,0x72,0x65,0x20,0x61,0x62,0x72,0x61,0x68,0x61,0x6d],UInt8[0x68,0x65,0x79,0x20,0x74,0x68,0x65,0x72,0x65,0x20,0x62,0x69,0x6c,0x6c],UInt8[0x68,0x65,0x79,0x20,0x74,0x68,0x65,0x72,0x65,0x20,0x63,0x68,0x61,0x72,0x6c,0x69,0x65],UInt8[0x68,0x65,0x79,0x20,0x74,0x68,0x65,0x72,0x65,0x20,0x64,0x61,0x6e],\"hey there ephraim\",\"hey there frank\",\"hey there george\",\"hey there hank\"\n"
end
rm("test1.csv")

# test exporting test2 to CSV
source = ODBC.Source(dsn, "select * from test2")
csv = CSV.Sink("test2.csv")
Data.stream!(source, csv)
open("test2.csv") do f
    @test readline(f) == "\"ID\",\"first_name\",\"last_name\",\"Salary\",\"hourly rate\",\"hireDate\",\"last clockin\"\n"
    @test readline(f) == "1,\"Lawrence\",\"Powell\",+87217E+0,26.47,2002-04-09,2002-01-17T21:32:00\n"
end

# test exporting test1 to SQLite
db = SQLite.DB()
source = ODBC.Source(dsn, "select * from test1")
SQLite.load(db, "test1", source)

data = SQLite.query(db, "select * from test1")
@test size(data) == (2,27)
@test data.columns[1][1] === Nullable(1)
@test data.columns[3][1] === Nullable(1.0)
@test data.columns[11][1] === Nullable(ODBC.API.SQLDate(2016,1,1))

# test exporting test2 to SQLite
source = ODBC.Source(dsn, "select * from test2")
SQLite.load(db, "test2", source)

data = SQLite.query(db, "select * from test2")
@test size(data) == (70000,7)
@test data.columns[1].values == [1:70000...]
@test data.columns[end][1] === Nullable(ODBC.API.SQLTimestamp(2002,1,17,21,32,0,0))

# test exporting test2 to Feather
file = "test2.feather"
Feather.write(file, ODBC.Source, dsn, "select ID, first_name, last_name, salary + 0E0 as salary, `hourly rate` from test2")
df = Feather.read(file)
@test size(df, 1) == size(data, 1)
@test df.columns[1] == [1:70000...]

ODBC.execute!(dsn, "create table test3 as select * from test2 limit 0")
ODBC.execute!(dsn, "delete from test3")
ODBC.query(dsn, "select * from test2 limit 100", ODBC.Sink, dsn, "test3")

df = ODBC.query(dsn, "select * from test3")
@test size(df) == (100,7)
@test df.columns[1].values == [1:100...]

ODBC.load(dsn, "test3", df; append=true)
df = ODBC.query(dsn, "select * from test3")
@test size(df) == (200,7)
@test df.columns[1].values == append!([1:100...], [1:100...])

stmt = ODBC.prepare(dsn, "insert into test3 values(?,?,?,?,?,?,?)")
ODBC.execute!(stmt, [101, "Steve", "McQueen", 1.0, 100.0, Date(2016,1,1), DateTime(2016,1,1)])

ODBC.execute!(stmt, [102, "Dean", "Martin", 1.5, 10.1, Date(2016,1,2), DateTime(2016,1,2)])

df = ODBC.query(dsn, "select * from test3")
@test size(df) == (202,7)
@test get(df.columns[1][end-1]) == 101
@test get(df.columns[1][end]) == 102
@test get(df.columns[2][end-1]) == "Steve"
@test get(df.columns[2][end]) == "Dean"
@test get(df.columns[3][end-1]) == "McQueen"
@test get(df.columns[3][end]) == "Martin"
@test get(df.columns[4][end-1]) == (is_windows() ? 1.0 : DecFP.Dec64(1))
@test get(df.columns[4][end]) == (is_windows() ? 2.0 : DecFP.Dec64(2))
@test get(df.columns[5][end-1]) == 100.0
@test get(df.columns[5][end]) == 10.1
@test get(df.columns[6][end-1]) == ODBC.API.SQLDate(2016,1,1)
@test get(df.columns[6][end]) == ODBC.API.SQLDate(2016,1,2)
@test get(df.columns[7][end-1]) == ODBC.API.SQLTimestamp(2016,1,1,0,0,0,0)
@test get(df.columns[7][end]) == ODBC.API.SQLTimestamp(2016,1,2,0,0,0,0)

ODBC.execute!(dsn, "drop table if exists test1")
ODBC.Source(dsn, "drop table if exists test2")
ODBC.Source(dsn, "drop table if exists test3")
rm("test2.csv")

# more robust Feather testing
mysqltypes = Dict(Int32=>"INT", Float64=>"DOUBLE", CategoricalArrays.NominalValue{String,Int32}=>"VARCHAR(255)", CategoricalArrays.OrdinalValue{String,Int32}=>"VARCHAR(255)", WeakRefString{UInt8}=>"VARCHAR(255)")
columndef(name, typ) = "`$name` $(mysqltypes[typ <: Nullable ? eltype(typ) : typ]),"

testnull{T1, T2}(v1::T1, v2::T2) = v1 == v2
testnull{T1, T2}(v1::Nullable{T1}, v2::Nullable{T2}) = (isnull(v1) && isnull(v2)) || (!isnull(v1) && !isnull(v2) && get(v1) == get(v2))
testnull{T1, T2}(v1::T1, v2::Nullable{T2}) = !isnull(v2) && get(v2) == v1
testnull{T1, T2}(v1::Nullable{T1}, v2::T2) = !isnull(v1) && get(v1) == v2

feather_dir = joinpath(dirname(@__FILE__), "../../Feather/test/data")
# feather_dir = "/Users/jacobquinn/.julia/v0.5/Feather/test/data"
# files = readdir(feather_dir)
for file in readdir(feather_dir)
    source = Feather.Source(joinpath(feather_dir,file))
    nms = Data.header(source)
    typs = Data.types(source)
    columndefs = ""
    for i = 1:length(nms)
        columndefs *= columndef(nms[i], typs[i])
    end
    createtable = "CREATE TABLE feather_test ($(chop(columndefs)))"
    ODBC.execute!(dsn, createtable)
    ODBC.load(dsn, "feather_test", source)
    df = ODBC.query(dsn, "select * from feather_test")
    df2 = Feather.read(source)
    # tests
    for (c1,c2) in zip(df.columns,df2.columns)
        for i = 1:length(c1)
            @test testnull(c1[i], c2[i])
        end
    end
    ODBC.execute!(dsn, "drop table feather_test")
end

ODBC.disconnect!(dsn)


text:
# MIT License

# Copyright (c) 2021 Dr. William A. Pisani

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# This program will read in the bonds from the designated file and then write out the
# topology (angles, dihedrals, impropers) to a file.
# Future improvements
# Enable reading from a standard LAMMPS data file
# Enable writing bonds, angles, dihedrals, and impropers to a LAMMPS data file
# Checking if correct number of arguments are supplied

# Include necessary libraries and functions
using StatsBase
include("angles.jl")
include("dihedrals.jl")
include("impropers.jl")
include("write_topology.jl")

bond_dir = pwd()
example_file = ARGS[1]

function parsebonds(directory::String,file::String)
    open(joinpath(directory,file),"r") do f
        lines = split(read(f,String),"\n")
        bonds = Array{Int64}[]
        for line in lines
            line = split(line," ")
            if length(line) == 4 # If bond type is present
                atom1 = parse(Int64,line[3])
                atom2 = parse(Int64,line[4])
            elseif length(line) == 3 # If bond type is not present
                atom1 = parse(Int64,line[2])
                atom2 = parse(Int64,line[3])
            end
            push!(bonds,[atom1,atom2])
        end
        return bonds
    end
end

bonds = parsebonds(bond_dir,example_file)

# Get angles from bonds
angle_array = angles(bonds)

# Get dihedrals
dihedral_array = dihedrals(bonds,angle_array)

# Get impropers
improper_array = impropers(bonds)

outfile = join([splitext(example_file)[1],"_topology_test.dat"])
write_topology(bond_dir,outfile,bonds,angle_array,dihedral_array,improper_array)

max_stars_repo_name: wapisani/CompMatSciTools

max_issues_repo_name: wapisani/CompMatSciTools

max_forks_repo_name: wapisani/CompMatSciTools


text:
println("Starting")

using MetaGen
using JSON
import YAML
using Pipe: @pipe
using Random

config_path = ARGS[1]
config = YAML.load_file(config_path)
output_dir = "results"
mkdir(output_dir)

include("useful_functions.jl")
dict = []
#for i = 0:config["batches_upto"]
	#to_add =  @pipe "$(config["input_file_dir"])$(i)_data_labelled.json" |> open |> read |> String |> JSON.parse
to_add =  @pipe "$(config["input_file"])" |> open |> read |> String |> JSON.parse
append!(dict, to_add)
#end
#dict = @pipe "../../scratch_work_07_16_21/0_data_labelled.json" |> open |> read |> String |> JSON.parse
#dict = @pipe "../../scratch_work_07_16_21/0_data_labelled.json" |> open |> read |> String |> JSON.parse

#Random.seed!(15)
#try to make objects_observed::Array{Array{Array{Array{Detection2D}}}} of observed objects.
#outer array is for scenes, then frames, the receptive fields, then last is an array of detections

################################################################################
num_videos = config["num_videos"]
num_frames = config["num_frames"]
threshold = config["threshold"]

categories_subset = config["categories_subset"]

params = Video_Params(n_possible_objects = length(categories_subset))
top_n = 5

receptive_fields = make_receptive_fields(params)
objects_observed, camera_trajectories = make_observations_office(dict, receptive_fields, categories_subset, num_videos, num_frames, threshold, top_n)

num_particles = config["num_particles"]
mcmc_steps_outer = config["mcmc_steps_outer"]
mcmc_steps_inner = config["mcmc_steps_inner"]
shuffle_type = config["shuffle_type"]
################################################################################
#Online MetaGen

num_videos_train = convert(Int64, num_videos/2)
half = convert(Int64, num_videos_train/2)

if shuffle_type==0
	order = collect(1:num_videos_train)
elseif shuffle_type==1
	order = vcat(reverse(collect((half+1):num_videos_train)), reverse(collect(1:half)))
elseif shuffle_type==2
	order = vcat(collect((half+1):num_videos_train), collect(1:half))
elseif shuffle_type==3
	order = vcat(reverse(collect(1:half)), reverse(collect((half+1):num_videos_train)))
end

training_objects_observed = objects_observed[order, :]
training_camera_trajectories = camera_trajectories[order, :]

#Set up the output files
online_V_file = open(output_dir * "/online_V.csv", "w")
file_header_V(online_V_file, params)
online_ws_file = open(output_dir * "/online_ws.csv", "w")
file_header_ws(online_ws_file, params, num_particles)

println("start pf for online")

traces, inferred_world_states, avg_v = unfold_particle_filter(nothing,
	num_particles, mcmc_steps_outer, mcmc_steps_inner, training_objects_observed,
	training_camera_trajectories, params, online_V_file, online_ws_file, order)
close(online_V_file)
close(online_ws_file)

println("avg_v ", avg_v)
println("done with pf for online")

################################################################################
#Retrospective MetaGen

#training set and test set
order = vcat(order, (num_videos_train+1):num_videos)
input_objects_observed = vcat(objects_observed[order, :])
input_camera_trajectories = vcat(camera_trajectories[order, :])

#Set up the output file
retro_V_file = open(output_dir * "/retro_V.csv", "w")
file_header_V(retro_V_file, params)
retro_ws_file = open(output_dir * "/retro_ws.csv", "w")
file_header_ws(retro_ws_file, params, num_particles)

println("start retrospective")

traces, inferred_world_states, avg_v = unfold_particle_filter(avg_v, num_particles, mcmc_steps_outer, mcmc_steps_inner,
	input_objects_observed, input_camera_trajectories, params, retro_V_file, retro_ws_file, order)
close(retro_V_file)
close(retro_ws_file)

println("done with pf for retrospective")

#=

################################################################################
#run Lesioned MetaGen

#Set up the output file
lesioned_outfile = output_dir * "/lesioned_output.csv"
lesioned_file = open(lesioned_outfile, "w")
file_header(lesioned_file)

v = zeros(length(params.possible_objects), 2)
v[:,1] .= 1.0
v[:,2] .= 0.5
unfold_particle_filter(v, num_particles, mcmc_steps_outer, mcmc_steps_inner,
	objects_observed, camera_trajectories, params, lesioned_file)
close(lesioned_file)

println("done with pf for lesioned metagen")

=#

################################################################################
#for writing an output file for a demo using Retro MetaGen.

#undor re-ordering of inferred_world_states
inferred_world_states = inferred_world_states[order]

###### add to dictionary
out = write_to_dict(dict, camera_trajectories, inferred_world_states, num_videos, num_frames)

#open("../../scratch_work_07_16_21/output_tiny_set_detections.json","w") do f
open(output_dir * "/output.json","w") do f
	JSON.print(f,out)
end

println("finished writing json")

max_stars_repo_name: zhangir-azerbayev/MetaGen

max_issues_repo_name: zhangir-azerbayev/MetaGen

max_forks_repo_name: zhangir-azerbayev/MetaGen


text:
include("../src/symbolics_module.jl")

module stest
using ..symbolics

symbolics.@variable v 1.0
symbolics.@variable z 0.2
symbolics.@variable x 0.5

sigma0 = 10

ex = v - z * sigma0 * abs(v)

str = symbolics.to_string(ex)

der = symbolics.derivative(ex, x)
str_der = symbolics.to_string(der)

sder = symbolics.simplify(der)
str_sder = symbolics.to_string(sder)

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
module DotaBot

using HTTP
using Random
using JSON
using CSV
using DataFrames
using Cambrian
import CartesianGeneticProgramming
import NEAT

include("util.jl")
include("server.jl")
include("action.jl")

end # module


text:
# %% CCM demonstration
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, Random
using Neighborhood, Statistics, LinearAlgebra

function ccm(x, y, d, τ, w = τ)
    Mx = embed(x, d, τ); theiler = Theiler(w); tree = KDTree(Mx)
    idxs = bulkisearch(tree, Mx, NeighborNumber(d+1), theiler)
    ỹ = copy(y)
    for i in 1:length(Mx)
        J = idxs[i]
        xᵢ = Mx[i]; n1 = norm(xᵢ - Mx[J[1]])
        w = [exp(-norm(xᵢ - Mx[j])/n1) for j in J]
        w ./= sum(w)
        ỹ[i] = sum(w[k]*y[j] for (k, j) in enumerate(J))
    end
    return cor(y, ỹ)
end

ds = Systems.lorenz()
tr = trajectory(ds, 1000; Δt = 0.05, Ttr = 1000)
x, y, z = columns(tr)
τx = estimate_delay(x, "mi_min")
τy = estimate_delay(y, "mi_min")

Ns = 50:50:3000
ρx = [ccm(x[1:N], y[1:N], 3, 5) for N in Ns]
ρy = [ccm(y[1:N], x[1:N], 3, 5) for N in Ns]
ρz = [ccm(z[1:N], y[1:N], 3, 5) for N in Ns]

# %%

fig = figure(figsize = (figx/2, figy))
ax = gca()
ax.plot(Ns, ρx; label = "\$x \\to y\$")
ax.plot(Ns, ρy; label = "\$y \\to x\$", color = "C2", ls = "dashed")
ax.set_ylim(0.68, 1)
ax.set_xticks(0:1000:3000)
ax.legend(loc = "lower right")
ax.set_ylabel("CCM measure")
ax.set_xlabel("\$N_s\$"; labelpad = -15)

fig.tight_layout(pad=0.3)
wsave(plotsdir("7", "ccm_result"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
using Literate

import Plots

function update_date(content)
    content = replace(content, "DATETIME" => Date(now()))
    return content
end

Literate.markdown(pwd() * "\\julia_sandbox\\literate_test.jl", pwd() * "\\julia_sandbox\\"; preprocess = update_date)
Literate.notebook(pwd() * "\\julia_sandbox\\literate_test.jl", pwd() * "\\julia_sandbox\\"; preprocess = update_date)
Literate.notebook(pwd() * "\\julia_sandbox\\literate_test_2.jl", pwd() * "\\julia_sandbox\\"; preprocess = update_date)


text:
using LightGraphs
using Plots
using DelimitedFiles

include("utils\\categorical_tree.jl")
include("SI\\contact_process.jl")

A = readdlm("test3reg.dat", ',', Int, '\n')

g = SimpleGraph(A)
N = nv(g)
β = 1.0
X = falses(N)
X[1] = true

ts, Xs = contact_process_gillespie(g, X, β)

plot(ts, Xs)

function test()
    t = zeros(1000)
    for i in 1:1000
        ts, Xs = contact_process_gillespie(g, X, β, tmax=9.0)
        k = 1
        while Xs[k]/N < 0.5
            k += 1
        end
        t[i] = ts[k]
    end
    return t
end

t = test()

using Statistics
mean(t)

function test2()
    tmean, Xmean, Mmean = contact_process_montecarlo(g, X, β, tmax = 10.0, nbins=1000, nsims = 1000)
    k = 1
    while Xmean[k]/N < 0.5
        k += 1
    end
    return tmean[k]
end

test2()

tmean, Xmean, Mmean = contact_process_montecarlo(g, X, β, tmax = 10.0, nbins=1000, nsims = 1000)

plot(tmean, Xmean/1000,
    label="Average (Contact Process)",
    legend=:bottomright,
    xlabel = "time",
    ylabel = "Fraction of infected individuals",
    )

max_stars_repo_name: csimal/metaplex-networks

max_issues_repo_name: csimal/metaplex-networks

max_forks_repo_name: csimal/metaplex-networks


text:
# This file was generated, do not modify it. # hide
Xs = source(X)
ys = source(y, kind=:target)

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
module acou_annulus_examples
using FinEtools
using FinEtools.AlgoAcoustModule

function acou_annulus_Q4_example_algo()
    println("""
    Annular region, pressure BC + rigid wall.
    This version uses the FinEtools algorithm module.
    Version: 08/21/2017
    """)

    t0 = time()

    rho = 1001*phun("kg/m^3");# mass density
    c  = 1500.0*phun("m/s");# sound speed
    bulk =  c^2*rho;
    omega =  2000*phun("rev/s");      # frequency of the piston
    rin =  1.0*phun("m");#internal radius

    rex =  2.0*phun("m"); #external radius
    nr = 20; nc = 120;
    Angle = 2*pi;
    thickness =  1.0*phun("m/s");
    tolerance = min(rin/nr,  rin/nc/2/pi)/10000;

    fens, fes = Q4annulus(rin, rex, nr, nc, Angle)
    fens, fes = mergenodes(fens,  fes,  tolerance);
    edge_fes = meshboundary(fes);

    # The pressure boundary condition
    l1 = selectelem(fens, edge_fes, box=[-1.1*rex -0.9*rex -0.5*rex 0.5*rex]);
    ebc1 = FDataDict("node_list"=>connectednodes(subset(edge_fes, l1)),
    "pressure"=>x -> cos(2*pi*x[2]/rin)+1im*sin(2*pi*x[2]/rin)) # entering the domain

    material = MatAcoustFluid(bulk, rho)
    femm = FEMMAcoust(IntegData(fes,  GaussRule(2, 2)),  material)
    region1 = FDataDict("femm"=>femm)

    # Make model data
    modeldata = FDataDict("fens"=>fens,
    "omega"=>omega,
    "regions"=>[region1], "essential_bcs"=>[ebc1]);

    # Call the solver
    modeldata = FinEtools.AlgoAcoustModule.steadystate(modeldata)
    geom=modeldata["geom"]
    P=modeldata["P"]
    println("Minimum/maximum pressure, real= $(minimum(real(P.values)))/$(maximum(real(P.values))))")
    println("Minimum/maximum pressure, imag= $(minimum(imag(P.values)))/$(maximum(imag(P.values))))")

    println("Total time elapsed = ",time() - t0,"s")

    # Postprocessing
    File = "acou_annulusmod.vtk"
    vtkexportmesh(File, fes.conn, geom.values, FinEtools.MeshExportModule.Q4; scalars=[("Pre", real(P.values)), ("Pim", imag(P.values))])
    @async run(`"paraview.exe" $File`)

end # acou_annulus_Q4_example_algo

function allrun()
    println("#####################################################")
    println("# acou_annulus_Q4_example_algo ")
    acou_annulus_Q4_example_algo()
    return true
end # function allrun

end # module acou_annulus_examples


text:
export getSensTMatVec

function getSensTMatVec(v::Vector,m::Vector,pFor::EikonalInvParam)

   # extract pointers
    Q     = pFor.Sources
    P     = pFor.Receivers
    nsrc = size(Q,2); nrec = size(P,2)

	n = size(Q,1);

    # allocate space for result
    JTv  = zeros(n)
	JTvi = zeros(n)
	pEiks = pFor.eikonalParams;

	if length(v)!=nsrc*nrec
		warn("wrong size of vector v in EikonalInv::getSensTMatVec.jl");
	end

    for i=1:nsrc
		t = P*v[(i-1)*nrec + 1 : i*nrec];
		if pFor.useFilesForFields
			file = matopen(getFieldsFileName());
			pEiks[i].T1 		= read(file,string("T1_",i));
			pEiks[i].ordering = read(file,string("ordering_",i));
			pEiks[i].OP 		= read(file,string("OP_",i));
			close(file);
		end
		getSensTMatVecEik(t,JTvi,pEiks[i]);
		if pFor.useFilesForFields
			FactoredEikonalFastMarching.clear!(pEiks[i]);
		end
        JTv += JTvi;
    end
    return JTv;
end


text:
"""
    abstract type Triangulation{Dc,Dp}

Abstract type representing an arbitrary tiling, tessellation,
or triangulation of a domain of parametric dimension `Dc` and
physical dimension `Dp`.

We define a triangulation from two basic ingredients:

- the cell-wise nodal coordinates of the cells in the triangulation, plus
- an interpolation of this cell-wise coordinates into the cells interior.

Note that this type represents general triangulations (not necessarily conforming),
which is the minimum geometrical information needed to perform cell-wise numerical integration.

The `Triangulation` interface is defined by overloading these methods:

- [`get_cell_coordinates(trian::Triangulation)`](@ref)
- [`get_reffes(trian::Triangulation)`](@ref)
- [`get_cell_type(trian::Triangulation)`](@ref)

Optional interface:

For triangulations living in a space of co-dimension 1, the following method can be defined:

- [`get_normal_vector(trian::Triangulation)`]

In some cases, concrete implementations want to override the default implementation of the following methods:

- [`restrict(f::AbstractArray, trian::Triangulation)`]
- [`get_cell_id(f::AbstractArray, trian::Triangulation)`]

The (mandatory) `Triangulation` interface can be tested with

- [`test_triangulation`](@ref)

"""
abstract type Triangulation{Dc,Dp} <: GridapType end

"""
    get_cell_coordinates(trian::Triangulation) -> AbstractArray{Vector{<:Point{Dp}}}
"""
function get_cell_coordinates(trian::Triangulation)
  @abstractmethod
end

"""
    get_reffes(trian::Triangulation) -> Vector{LagrangianRefFE}
"""
function get_reffes(trian::Triangulation)
  @abstractmethod
end

"""
    get_cell_type(trian::Triangulation) -> AbstractVector{<:Integer}
"""
function get_cell_type(trian::Triangulation)
  @abstractmethod
end

"""
    get_normal_vector(trian::Triangulation)
"""
function get_normal_vector(trian::Triangulation)
  Dp = num_point_dims(trian)
  Dc = num_cell_dims(trian)
  if Dp == Dc + 1
    @abstractmethod
  else
    @unreachable "get_normal_vector does not make sense for this triangulation"
  end
end

"""
    test_triangulation(trian::Triangulation)
"""
function test_triangulation(trian::Triangulation{Dc,Dp}) where {Dc,Dp}
  @test num_cell_dims(trian) == Dc
  @test num_point_dims(trian) == Dp
  @test num_cell_dims(typeof(trian)) == Dc
  @test num_point_dims(typeof(trian)) == Dp
  cell_coords = get_cell_coordinates(trian)
  @test isa(cell_coords,AbstractArray{<:AbstractVector{<:Point}})
  reffes = get_reffes(trian)
  @test isa(reffes,AbstractVector{<:LagrangianRefFE{Dc}})
  cell_types = get_cell_type(trian)
  @test isa(cell_types,AbstractArray{<:Integer})
  ncells = num_cells(trian)
  @test ncells == length(cell_coords)
  @test ncells == length(cell_types)
end

# Some API

"""
    num_cells(trian::Triangulation) -> Int
"""
num_cells(trian::Triangulation) = length(get_cell_type(trian))

"""
    num_cell_dims(::Triangulation) -> Int
    num_cell_dims(::Type{<:Triangulation}) -> Int
"""
num_cell_dims(::Triangulation{Dc,Dp}) where {Dc,Dp} = Dc
num_cell_dims(::Type{<:Triangulation{Dc,Dp}}) where {Dc,Dp} = Dc

"""
    num_point_dims(::Triangulation) -> Int
    num_point_dims(::Type{<:Triangulation}) -> Int
"""
num_point_dims(::Triangulation{Dc,Dp}) where {Dc,Dp} = Dp
num_point_dims(::Type{<:Triangulation{Dc,Dp}}) where {Dc,Dp} = Dp

"""
    num_dims(::Triangulation) -> Int
    num_dims(::Type{<:Triangulation}) -> Int

Equivalent to `num_cell_dims`.
"""
num_dims(g::Triangulation{Dc}) where Dc = Dc
num_dims(::Type{<:Triangulation{Dc}}) where Dc = Dc

"""
    is_affine(trian::Triangulation) -> Bool
"""
function is_affine(trian::Triangulation)
  reffes = get_reffes(trian)
  all(map(is_affine,reffes))
end

"""
    is_first_order(trian::Triangulation) -> Bool
"""
function is_first_order(trian::Triangulation)
  reffes = get_reffes(trian)
  all(map(is_first_order,reffes))
end

"""
    reindex(a::AbstractArray, trian::Triangulation)
"""
function reindex(a::AbstractArray,trian::Triangulation)
  reindex(a,get_cell_id(trian))
end

"""
    get_cell_id(trian::Triangulation)
"""
function get_cell_id(trian::Triangulation)
  identity_vector(num_cells(trian))
end

"""
    restrict(f::AbstractArray, trian::Triangulation)
"""
function restrict(f::AbstractArray,trian::Triangulation)
  f
end

"""
    get_cell_reffes(trian::Triangulation) -> Vector{<:LagrangianRefFE}

It is not desirable to iterate over the resulting array
for large number of cells if the underlying reference FEs
are of different Julia type.
"""
function get_cell_reffes(trian::Triangulation)
  type_to_reffe = get_reffes(trian)
  cell_to_type = get_cell_type(trian)
  _get_cell_data(type_to_reffe,cell_to_type)
end

"""
    get_cell_shapefuns(trian::Triangulation) -> Vector{<:Field}
"""
function get_cell_shapefuns(trian::Triangulation)
  type_to_reffes = get_reffes(trian)
  cell_to_type = get_cell_type(trian)
  type_to_shapefuns = map(get_shapefuns, type_to_reffes)
  _get_cell_data(type_to_shapefuns,cell_to_type)
end

"""
    get_cell_map(trian::Triangulation) -> Vector{<:Field}
"""
function get_cell_map(trian::Triangulation)
  cell_to_coords = get_cell_coordinates(trian)
  cell_to_shapefuns = get_cell_shapefuns(trian)
  lincomb(cell_to_shapefuns, cell_to_coords)
end

"""
"""
function CellField(object,trian::Triangulation)
  cm = get_cell_map(trian)
  convert_to_cell_field(object,cm)
end

"""
    get_physical_coordinate(trian::Triangulation)

In contrast to get_cell_map, the returned object:
- is a [`CellField`](@ref)
- its gradient is the identity tensor
"""
function get_physical_coordinate(trian::Triangulation)
  CellField(_phys_coord,trian)
end

_phys_coord(x) = x

_phys_coord_grad(x) = one(typeof(outer(x,x)))

gradient(::typeof(_phys_coord)) = _phys_coord_grad

# Helpers for Triangulation

function _get_cell_data(type_to_data, cell_to_type)
  CompressedArray(type_to_data,cell_to_type)
end

function _get_cell_data(type_to_data, cell_to_type::Fill)
  ncells = length(cell_to_type)
  @assert length(type_to_data) == 1 "Only one reference element expected"
  @assert cell_to_type.value == 1 "Only one type of reference element expected"
  data = first(type_to_data)
  Fill(data,ncells)
end

"""
    restrict(cf::CellField,trian::Triangulation)
"""
function restrict(cf::CellField,trian::Triangulation)
  _cf = to_ref_space(cf)
  a = get_array(_cf)
  r = restrict(a,trian)
  axs = reindex(get_cell_axes(cf),trian)
  _restrict_cell_field(r,axs,MetaSizeStyle(cf),trian)
end

function _restrict_cell_field(r::AbstractArray,axs::AbstractArray,msize_style::Val,trian)
  cm = get_cell_map(trian)
  GenericCellField(r,cm,Val(true),axs,msize_style)
end

function _restrict_cell_field(r::SkeletonPair,axs::SkeletonPair,msize_style::Val,trian)
  cm = get_cell_map(trian)
  la = r.left
  ra = r.right
  l = GenericCellField(la,cm,Val(true),axs.left,msize_style)
  r = GenericCellField(ra,cm,Val(true),axs.right,msize_style)
  merge_cell_fields_at_skeleton(l,r)
end

"""
    CellQuadrature(trian::Triangulation, degree::Integer)
"""
function CellQuadrature(trian::Triangulation, degree::Integer)
  polytopes = map(get_polytope,get_reffes(trian))
  cell_type = get_cell_type(trian)
  CellQuadrature(degree,polytopes,cell_type)
end

"""
    integrate(cell_field,trian::Triangulation,quad::CellQuadrature)

The `cell_field` is aligned with the cells in `trian`
"""
function integrate(cell_field,trian::Triangulation,quad::CellQuadrature)
  cell_map = get_cell_map(trian)
  integrate(cell_field,cell_map,quad)
end

function CellField(value::Number,trian::Triangulation,quad::CellQuadrature)
  CellField(value,get_cell_map(trian),quad)
end

"""
"""
function cell_measure(trian_Ω1::Triangulation,n_bgcells::Integer)
  trian_cut_1 = trian_Ω1
  quad_cut_1 = CellQuadrature(trian_cut_1,0)
  subcell1_to_dV = integrate(1,trian_cut_1,quad_cut_1)
  subcell1_to_bgcell = get_cell_id(trian_cut_1)
  bgcell_to_dV = zeros(n_bgcells)
  _meas_K_fill!(bgcell_to_dV,subcell1_to_dV,subcell1_to_bgcell)
  bgcell_to_dV
end

function cell_measure(trian_Ω1::Triangulation,trian::Triangulation)
  cell_measure(trian_Ω1,num_cells(trian))
end

function _meas_K_fill!(bgcell_to_dV,subcell1_to_dV,subcell1_to_bgcell)
  for (subcell1, bgcell) in enumerate(subcell1_to_bgcell)
    bgcell_to_dV[bgcell] += subcell1_to_dV[subcell1]
  end
end


text:
# Algebra generator for CliffordAlgebras.jl

import Base.show
import Combinatorics.levicivita
import PrettyTables.pretty_table

"""
    basevectorproduct(Npos::Integer, Nneg::Integer, Nzero::Integer, base::Tuple, kl::Integer, kr::Integer)

Finds the product of the kl-th base vector with the kr-th base vector. Returns a base vector index and a scale factor.
"""
function basevectorproduct(Npos::Integer, Nneg::Integer, base::Tuple, kl::Integer, kr::Integer)
    K = length(base)
    @assert kr <= K && kl <= K
    selector_left = base[kl]
    selector_right = base[kr]
    selector_prod = collect((selector_left..., selector_right...))
    permutation = sortperm(selector_prod)
    coeff = levicivita(permutation)
    v = []
    for b in unique(selector_prod[permutation])
        c = count(isequal(b), selector_prod)
        @assert c in (1, 2)
        if c == 1
            push!(v, b)
        else
            coeff *= (b <= Npos) ? 1 : (b <= Npos + Nneg) ? -1 : 0
        end
    end
    v = Tuple(v)
    vl = length(v)
    bi = findfirst(b -> length(b) == vl && length(intersect(b, v)) == vl, base)
    coeff *= levicivita(findpermutation(v, base[bi]))
    (bi, coeff)
end

"""
    multiplicationstable(Npos::Integer, Nneg::Integer, Nzero::Integer, base::Tuple)

Generates a multiplication table for the Clifford algebra with the specified signature and base vectors.
"""
function multiplicationtable(Npos::Integer, Nneg::Integer, Nzero::Integer, base::Tuple)
    N = Npos + Nneg + Nzero
    K = 2^N
    @assert length(base) == K
    M = Matrix{Tuple}(undef, K, K)
    for kl = 1:K, kr = 1:K
        M[kl, kr] = basevectorproduct(Npos, Nneg, base, kl, kr)
    end
    ntuple(row -> ntuple(col -> M[row, col], K), K)
end

"""
    CliffordAlgebra(Npos::Integer, Nneg::Integer, Nzero::Integer, S::NTuple(N,Symbol))

Singleton instance of the type CliffordAlgebra that describes a geometric algebra with the signature (Npos,Nneg,Nzero), base symbols S.
The base symbols are in order of the signature.
"""
struct CliffordAlgebra{Np,Nn,Nz,S,BT}
    function CliffordAlgebra(
        Npos::Integer,
        Nneg::Integer,
        Nzero::Integer,
        BaseSymbols::NTuple{N,Symbol},
    ) where {N}
        Npos = Int(Npos)
        Nneg = Int(Nneg)
        Nzero = Int(Nzero)
        @assert Npos >= 0 && Nneg >= 0 && Nzero >=0 "Algebra signature must be non-negative."
        @assert Npos + Nneg + Nzero == N "Base symbol count must match signature."
        BT = adaptbasefordual(enumeratebase(Int(N)))
        new{Npos,Nneg,Nzero,BaseSymbols,BT}()
    end
end

"""
    CliffordAlgebra(N::Integer)

Generates a geometric algebra with signature (N,0,0).
"""
function CliffordAlgebra(N::Integer)
    @assert N>=0 "Algbra signature must be non-negative."
    CliffordAlgebra(N, 0, 0, ntuple(i -> Symbol(:e, i), N))
end

"""
    CliffordAlgebra(Npos::Integer, Nneg::Integer)

Generates a geometric algebra with signature (Npos,Nneg,0).
"""
function CliffordAlgebra(Npos::Integer, Nneg::Integer)
    @assert Npos >= 0 && Nneg >= 0 "Algbra signature must be non-negative."
    CliffordAlgebra(Npos, Nneg, 0, ntuple(i -> Symbol(:e, i), Npos + Nneg))
end

"""
    CliffordAlgebra(Npos::Integer, Nneg::Integer, Nzero::Integer)

Generates a geometric algebra with signature (Npos,Nneg,Nzero).
"""
function CliffordAlgebra(Npos::Integer, Nneg::Integer, Nzero::Integer)
    @assert Npos >= 0 && Nneg >= 0 && Nzero >= 0 "Algbra signature must be non-negative."
    CliffordAlgebra(Npos, Nneg, Nzero, ntuple(i -> Symbol(:e, i), Npos + Nneg + Nzero))
end

"""
    CliffordAlgebra(a::Symbol)

Generates a predefined algebra from a identifier. Known algebras are
    - :Hyperbolic or :Hyper
    - :Complex or :ℂ
    - :Dual or :Grassmann
    - :Grassmann2D or :G2
    - :Grassmann3D or :G3
    - :Quaternions or :ℍ
    - :Cl2 and :Cl3
    - :Spacetime
    - :PGA2D or :Projective2D or :Plane2D
    - :PGA3D or :Projective3D or :Plane3D
    - :CGA2D or :Conformal2D
    - :CGA3D or :Conformal3D
    - :DCGA3D or :DoubleConformal3D
    - :TCGA3D or :TripleConformal3D
    - :DCGSTA or :DoubleConformalSpacetime
    - :QCGA or :QuadricConformal
"""
function CliffordAlgebra(a::Symbol)
    if a in (:Hyperbolic, :Hyper)
        return CliffordAlgebra(1, 0, 0, (:j,))
    elseif a in (:Complex, :ℂ)
        return CliffordAlgebra(0, 1, 0, (:i,))
    elseif a in (:Dual, :Grassmann)
        return CliffordAlgebra(0, 0, 1, (:ε,))
    elseif a in (:Grassmann2D, :G2)
        return CliffordAlgebra(0, 0, 2, (:ε₁, :ε₂))
    elseif a in (:Grassmann3D, :G3)
        return CliffordAlgebra(0, 0, 3, (:ε₁, :ε₂, :ε₃))
    elseif a in (:Quaternions, :ℍ)
        return CliffordAlgebra(0, 2, 0, (:i, :j))
    elseif a in (:Cl2,)
        return CliffordAlgebra(2)
    elseif a in (:Cl3,)
        return CliffordAlgebra(3)
    elseif a in (:Spacetime, :STA)
        return CliffordAlgebra(1, 3, 0, (:t, :x, :y, :z))
    elseif a in (:PGA2D, :Projective2D, :Plane2D)
        return CliffordAlgebra(2, 0, 1, (:e1, :e2, :e0))
    elseif a in (:PGA3D, :Projective3D, :Plane3D)
        return CliffordAlgebra(3, 0, 1, (:e1, :e2, :e3, :e0))
    elseif a in (:CGA2D, :Conformal2D)
        return CliffordAlgebra(3, 1, 0, (:e1, :e2, :e₊, :e₋))
    elseif a in (:CGA3D, :Conformal3D)
        return CliffordAlgebra(4, 1, 0, (:e1, :e2, :e3, :e₊, :e₋))
    elseif a in (:DCGA3D, :DoubleConformal3D)
        return CliffordAlgebra(6, 2)
    elseif a in (:TCGA3D, :TripleConformal3D)
        return CliffordAlgebra(9, 3)
    elseif a in (:DCGSTA, :DoubleConformalSpacetime)
        return CliffordAlgebra(4, 8, 0, (:t₁, :t₂, :e₊₁, :e₊₂, :x₁, :x₂, :y₁, :y₂, :z₁, :z₂, :e₋₁, :e₋₂))
    elseif a in (:QCGA, :QuadricConformal)
        return CliffordAlgebra(9, 6)
    else
        throw(ArgumentError("Unknown algebra."))
    end
end

"""
    basesymbols(::CliffordAlgebra)
    basesymbols(::Type{<:CliffordAlgebra})

Returns the 1-vector space basis symbols of the algebra.
"""
basesymbols(::Type{<:CliffordAlgebra{Np,Nn,Nz,S}}) where {Np,Nn,Nz,S} = S
basesymbols(ca::CliffordAlgebra) = basesymbols(typeof(ca))

"""
    basetable(::CliffordAlgebra)
    basetable(::Type{<:CliffordAlgebra})

Returns the internal basis table of the algebra.
"""
basetable(::Type{<:CliffordAlgebra{Np,Nn,Nz,S,BT}}) where {Np,Nn,Nz,S,BT} = BT
basetable(ca::CliffordAlgebra) = basetable(typeof(ca))

"""
    multtable(::CliffordAlgebra)
    multtable(::Type{<:CliffordAlgebra})

Returns the internal multuplication table of the geometric product of the algebra.
"""
multtable(::Type{<:CliffordAlgebra{Np,Nn,Nz,S,BT}}) where {Np,Nn,Nz,S,BT} = multiplicationtable(Np, Nn, Nz, BT)
multtable(ca::CliffordAlgebra) = multtable(typeof(ca))

"""
    baseproduct(::CliffordAlgebra, nleft::Integer, nright::Integer)
    baseproduct(::Type{<:CliffordAlgebra}, nleft::Integer, nright::Integer)

Returns a tuple holding the basis index and the scalar coefficient of the geometric product of the two basis vectors with indices nleft and nright.
"""
baseproduct(::Type{CliffordAlgebra{Np,Nn,Nz,S,BT}}, nleft::Integer, nright::Integer) where {Np,Nn,Nz,S,BT} = basevectorproduct(Np, Nn, BT, nleft, nright)
baseproduct(ca::CliffordAlgebra, nleft::Integer, nright::Integer) = baseproduct(typeof(ca), nleft, nright)

"""
    order(::CliffordAlgebra)
    order(::Type{<:CliffordAlgebra})

Returns the order of the algebra. The order is the sum of the signature.
"""
order(::Type{<:CliffordAlgebra{Np,Nn,Nz}}) where {Np,Nn,Nz} = Np + Nn + Nz
order(ca::CliffordAlgebra) = order(typeof(ca))

"""
    signature(::CliffordAlgebra)
    signature(::Type{<:CliffordAlgebra})

Returns the signature of the algebra.
"""
signature(::Type{<:CliffordAlgebra{Np,Nn,Nz}}) where {Np,Nn,Nz} = (Np, Nn, Nz)
signature(ca::CliffordAlgebra) = signature(typeof(ca))

"""
    basesignature(::CliffordAlgebra, n::Integer)
    basesignature(::Type{<:CliffordAlgebra}, n::Integer)

Returns the signature value of the n-th basis 1-vector. The return value is +1, -1 or 0.
"""
function basesignature(CA::Type{<:CliffordAlgebra{Np,Nn,Nz}}, n::Integer) where {Np,Nn,Nz}
    @assert 1 <= n <= order(CA)
    if n <= Np
        return +1
    elseif n <= Np + Nn
        return -1
    else
        return 0
    end
end

basesignature(ca::CliffordAlgebra, n::Integer) = basesignature(typeof(ca), n)

"""
    dimension(::CliffordAlgebra)
    dimension(::Type{<:CliffordAlgebra})

Returns the dimension of the algebra, i.e. the number of coefficients in a general multivector.
"""
dimension(CA::Type{<:CliffordAlgebra}) = 2^order(CA)
dimension(ca::CliffordAlgebra) = dimension(typeof(ca))

"""
    character(::CliffordAlgebra)
    character(::Type{<:CliffordAlgebra})

Returns the square of the pseudoscalar of the algebra.
"""
character(CA::Type{<:CliffordAlgebra}) = baseproduct(CA, dimension(CA), dimension(CA))[2]
character(ca::CliffordAlgebra) = character(typeof(ca))

"""
    basegrade(::CliffordAlgebra, n::Integer)
    basegrade(::Type{<:CliffordAlgebra}, n::Integer)

Returns the grade of the n-th basis multivector of the algebra.
"""
basegrade(CA::Type{<:CliffordAlgebra}, n::Integer) = length(basetable(CA)[n])
basegrade(ca::CliffordAlgebra, n::Integer) = basegrade(typeof(ca), n)

"""
    basereverse(::CliffordAlgebra, n::Integer)
    basereverse(::Type{<:CliffordAlgebra}, n::Integer)

Returns the sign change of the n-th basis multivector under reversal.
"""
basereverse(CA::Type{<:CliffordAlgebra}, n::Integer) = (-1)^(basegrade(CA, n) ÷ 2)
basereverse(ca::CliffordAlgebra, n::Integer) = basereverse(typeof(ca), n)

"""
    basesymbol(::CliffordAlgebra, n::Integer)
    basesymbol(::Type{<:CliffordAlgebra}, n::Integer)

Returns the symbol used for the n-th basis multivector of the algebra.
"""
function basesymbol(CA::Type{<:CliffordAlgebra}, n::Integer)
    @assert 1 <= n <= dimension(CA)
    s = basetable(CA)[n]
    if isempty(s)
        Symbol(:𝟏)
    else
        Symbol(map(i -> basesymbols(CA)[i], s)...)
    end
end

basesymbol(ca::CliffordAlgebra, n::Integer) = basesymbol(typeof(ca), n)

"""
    caleytable(io::IO, ca::CliffordAlgebra)
    caleytable(io::IO, CA::Type{<:CliffordAlgebra})

Generates a Cayley table view of the algebra.
"""
cayleytable(io::IO, ca::CliffordAlgebra) = cayleytable(io, typeof(ca))

function cayleytable(io::IO, CA::Type{<:CliffordAlgebra})
    dim = dimension(CA)
    bs = map(s -> s == Symbol() ? Symbol("1") : s, ntuple(k -> basesymbol(CA, k), dim))
    mt = multtable(CA)
    table = Matrix{String}(undef, dim, dim)
    for col = 1:dim
        for row = 1:dim
            baseidx, coeff = mt[row][col]
            @assert coeff in (-1, 0, 1)
            if baseidx == 1
                table[row, col] = coeff == -1 ? "-1" : coeff == 1 ? "+1" : "0"
            else
                if iszero(coeff)
                    table[row, col] = "0"
                else
                    table[row, col] = string(coeff > 0 ? :+ : :-, bs[baseidx])
                end
            end
        end
    end
    lines = vcat(0, cumsum([binomial(order(CA), i) for i = 0:order(CA)]))
    pretty_table(
        io,
        table;
        show_row_number = false,
        crop = :none,
        noheader = true,
        alignment = :c,
        vlines = lines,
        hlines = lines,
    )
end

"""
    signaturetable(io::IO, ca::CliffordAlgebra)
    signaturetable(io::IO, CA::Type{<:CliffordAlgebra})

Prints the 1-vector basis symbols and their squares.
"""
signaturetable(io::IO, ca::CliffordAlgebra) = signaturetable(io, typeof(ca))

function signaturetable(io::IO, CA::Type{<:CliffordAlgebra})
    bs = ntuple( k -> basesymbol(CA, k+1), order(CA))
    table = Matrix{String}(undef, order(CA), 2)
    for (n,s) = enumerate(bs)
        sig = basesignature(CA,n)
        table[n,1] = string(s)
        table[n,2] = sig == +1 ? "+1" : sig == -1 ? "-1" : "0"
    end
    pretty_table(
        io,
        table;
        show_row_number = false,
        crop = :none,
        noheader = true,
        alignment = :c
    )
end

function show(io::IO, ca::CliffordAlgebra)
    (Np, Nn, Nz) = signature(ca)
    print(io, "Cl(", Np, ",", Nn, ",", Nz, ")")
end

max_stars_repo_name: ATell-SoundTheory/CliffordAlgebras

max_issues_repo_name: ATell-SoundTheory/CliffordAlgebras


text:
export JLMatch, @syntax_pattern, codegen_match

"""
    JLMatch <: JLExpr

`JLMatch` describes a Julia pattern match expression defined by
[`MLStyle`](https://github.com/thautwarm/MLStyle.jl). It allows
one to construct such expression by simply assign each code block
to the corresponding pattern expression.

!!! tip
    `JLMatch` is not available in ExproniconLite since it depends
    on MLStyle's pattern matching functionality.

# Example

One can construct a `MLStyle` pattern matching expression
easily by assigning the corresponding pattern and its result
to the `map` field.

```julia
julia> jl = JLMatch(:x)
#= line 0 =#
nothing

julia> jl = JLMatch(:x)
#= line 0 =#
nothing

julia> jl.map[1] = true
true

julia> jl.map[2] = :(sin(x))
:(sin(x))

julia> jl
#= line 0 =#
@match x begin
    1 => true
    2 => sin(x)
    _ =>     nothing
end
```

to generate the corresponding Julia `Expr` object, one can call [`codegen_ast`](@ref).

```julia
julia> codegen_ast(jl)
:(let
      true
      var"##return#263" = nothing
      var"##265" = x
      if var"##265" isa Int64
          #= line 0 =#
          if var"##265" === 1
              var"##return#263" = let
                      true
                  end
              #= unused:1 =# @goto var"####final#264#266"
          end
          #= line 0 =#
          if var"##265" === 2
              var"##return#263" = let
                      sin(x)
                  end
              #= unused:1 =# @goto var"####final#264#266"
          end
      end
      #= line 0 =#
      begin
          var"##return#263" = let
                  nothing
              end
          #= unused:1 =# @goto var"####final#264#266"
      end
      (error)("matching non-exhaustive, at #= line 0 =#")
      #= unused:1 =# @label var"####final#264#266"
      var"##return#263"
  end)
```
"""
struct JLMatch <: JLExpr
    item::Any
    patterns::Vector{Any}
    actions::Vector{Any}
    fallthrough::Any
    mod::Module
    line::LineNumberNode
end

"""
    JLMatch(item)

Generate an empty `JLMatch` object with given item expression.
`item` can be a `Symbol` or an `Expr`.
"""
JLMatch(item) = JLMatch(item, [], [], nothing, Main, LineNumberNode(0))

"""
    JLMatch(;kw...)

Create a `JLMatch` object from keyword arguments.

# Kwargs

- `item`: item to match
- `map`: the pattern=>result map, should be an `OrderedDict`.
- `fallthrough`: the result of fallthrough pattern `_`.
- `mod`: module to evaluate the expression.
- `line`: line number `LineNumberNode`.
"""
JLMatch(;item, patterns=[], actions=[], fallthrough=nothing, mod=Main, line=LineNumberNode(0)) =
    JLMatch(item, patterns, actions, fallthrough, mod, line)

Base.length(jl::JLMatch) = length(jl.patterns)

function Base.getindex(jl::JLMatch, pattern)
    idx = findfirst(jl.patterns) do p
        compare_expr(p, pattern)
    end
    idx === nothing && error("cannot find pattern: $pattern")
    return jl.actions[idx]
end

function Base.setindex!(jl::JLMatch, action, pattern)
    idx = findfirst(jl.patterns) do p
        compare_expr(p, pattern)
    end

    if idx === nothing
        push!(jl.patterns, pattern)
        push!(jl.actions, action)
    else
        jl.actions[idx] = action
    end
    return action
end

function Base.iterate(jl::JLMatch, st=1)
    st > length(jl) && return
    return jl.patterns[st] => jl.actions[st], st + 1
end

function codegen_ast(def::JLMatch)
    isempty(def.patterns) && return def.fallthrough
    body = Expr(:block)
    for (pattern, code) in def
        push!(body.args, :($pattern => $code))
    end
    push!(body.args, :(_ => $(def.fallthrough)))
    return init_cfg(gen_match(def.item, body, def.line, def.mod))
end

function print_expr(io::IO, def::JLMatch, ps::PrintState, theme::Color)
    print_expr(io, def.line, ps, theme)
    println(io)
    isempty(def.patterns) && return print_expr(io, def.fallthrough, ps, theme)

    within_line(io, ps) do
        print_kw(io, "@match ", ps, theme)
        print_expr(io, def.item, ps, theme)
        print_kw(io, " begin ", ps, theme)
    end
    println(io, ps)
    within_indent(ps) do
        for (k, (pattern, action)) in enumerate(def)
            within_line(io, ps) do
                print_expr(io, pattern, ps, theme)
                print_kw(io, " => ", ps, theme)
                print_expr(io, action, ps, theme)
            end
            println(io, ps)
        end
        within_line(io, ps) do
            print(io, "_")
            print_kw(io, " => ", ps, theme)
            print_expr(io, def.fallthrough, ps, theme)
        end
    end
    println(io, ps)
    print_end(io, ps, theme)
end

"""
    codegen_match(f, x[, line::LineNumberNode=LineNumberNode(0), mod::Module=Main])

Generate a zero dependency match expression using MLStyle code generator,
the syntax is identical to MLStyle.

!!! tip
    `codegen_match` is not available in ExproniconLite since it depends
    on MLStyle's pattern matching functionality.

# Example

```julia
codegen_match(:x) do
    quote
        1 => true
        2 => false
        _ => nothing
    end
end
```

This code generates the following corresponding MLStyle expression

```julia
@match x begin
    1 => true
    2 => false
    _ => nothing
end
```
"""
function codegen_match(f, x, line::LineNumberNode=LineNumberNode(0), mod::Module=Main)
    return init_cfg(gen_match(x, f(), line, mod))
end

"""
    @syntax_pattern <syntax type> <syntax checker>

Bundle a syntax type as MLStyle pattern.

!!! tip
    this is not available in ExproniconLite since it depends
    on MLStyle's pattern matching functionality.

# Example

```julia
struct MyFunction
    ex :: Expr
    test_field :: Any
end

is_xxx(ex::Expr) = Meta.isexpr(ex, :function)
is_xxx(_) = false

julia> MyFunction(ex::Expr) = MyFunction(ex, :aaa)
julia> @syntax_pattern(MyFunction, is_xxx)

julia> @match :(function f() end) begin
        MyFunction(;test_field) => test_field
    end
:aaa
```
"""
macro syntax_pattern(type, checker)
    ex = quote
        $__source__
        function $MLStyle.pattern_uncall(
            t::Type{$type},
            self::Function,
            type_params::Any,
            type_args::Any,
            args::Any,
        )
            $__source__
            type_infer(_...) = Union{Expr, t}  # type of input
            transform(expr) = :($($checker)($expr) ? $($type)($expr) : nothing)
            $MLStyle.AbstractPatterns.P_fast_view(
                type_infer,
                transform,
                $MLStyle.Record._compile_record_pattern(t, self, type_params, type_args, args)
            )
        end
    end
    return esc(ex)
end

@syntax_pattern JLFunction  is_function
@syntax_pattern JLStruct    is_struct_not_kw_struct
@syntax_pattern JLKwStruct  is_struct
@syntax_pattern JLIfElse    is_ifelse
@syntax_pattern JLFor       is_for


text:
module Testing

using Test, Printf

using ACEbase.FIO: read_dict, write_dict, save_dict, load_dict
using LinearAlgebra: norm
using StaticArrays

export print_tf, test_fio, h0, h1, h2, h3, println_slim

function h0(str)
   dashes = "≡"^(length(str)+4)
   printstyled(dashes, color=:magenta); println()
   printstyled("  "*str*"  ", bold=true, color=:magenta); println()
   printstyled(dashes, color=:magenta); println()
end

function h1(str)
   dashes = "="^(length(str)+2)
   printstyled(dashes, color=:magenta); println()
   printstyled(" " * str * " ", bold=true, color=:magenta); println()
   printstyled(dashes, color=:magenta); println()
end

function h2(str)
   dashes = "-"^length(str)
   printstyled(dashes, color=:magenta); println()
   printstyled(str, bold=true, color=:magenta); println()
   printstyled(dashes, color=:magenta); println()
end

h3(str) = (printstyled(str, bold=true, color=:magenta); println())

print_tf(::Test.Pass) = printstyled("+", bold=true, color=:green)
print_tf(::Test.Fail) = printstyled("-", bold=true, color=:red)
print_tf(::Tuple{Test.Error,Bool}) = printstyled("x", bold=true, color=:magenta)

println_slim(::Test.Pass) = printstyled("Test Passed\n", bold=true, color=:green)
println_slim(::Test.Fail) = printstyled("Test Failed\n", bold=true, color=:red)

"""
`test_fio(obj): `  performs two tests:

- encodes `obj` as a Dict using `write_dict`, then decodes it using
`read_dict` and tests whether the two objects are equivalent using `==`
- writes `Dict` to file then reads it and decodes it and test the result is
again equivalent to `obj`

The two results are returned as Booleans.
"""
function test_fio(obj; warntype = true)
   D = write_dict(obj)
   test1 = (obj == read_dict(D))
   if !test1
      @warn("test_fio (1) fails - read_dict(write_dict(obj)) != obj")
   end
   tmpf = tempname() * ".json"
   save_dict(tmpf, D)
   D2 = load_dict(tmpf)
   # if D != D2
   #    @warn("load_dict(save_dict(D)) != D")
   # end
   obj2 = read_dict(load_dict(tmpf))
   if warntype && (typeof(obj) != typeof(obj2))
      @warn(
      """test_fio: the loaded object does not have the same type
             original : $(typeof(obj))
         deserialised : $(typeof(obj2))
      """)
   end
   test2 = (obj == obj2)
   if !test2
      @warn("test_fio (2) fails - obj2 != obj")
   end
   return test1, test2
end

_Vec(X::AbstractVector{<: StaticVector{3}}) =
      collect(reinterpret(Float64, X))

_svecs(x::AbstractVector{T}) where {T} =
      collect(reinterpret(SVector{3, T}, x))

fdtest(F, dF, X::AbstractVector{<: StaticVector{3}}; kwargs...) =
      fdtest( x -> F(_svecs(x)),
              x -> _Vec(dF(_svecs(x))),
              _Vec(X); kwargs... )

fdtest(F, dF, X::Number; kwargs...) =
      fdtest( x -> F(x[1]),
              x -> [dF(x[1])],
              [X]; kwargs... )

"""
first-order finite-difference test for scalar F
```julia
fdtest(F, dF, x; h0 = 1.0, verbose=true)
```
"""
function fdtest(F, dF, x::AbstractVector; h0 = 1.0, verbose=true)
   errors = Float64[]
   E = F(x)
   dE = dF(x)
   # loop through finite-difference step-lengths
   verbose && @printf("---------|----------- \n")
   verbose && @printf("    h    | error \n")
   verbose && @printf("---------|----------- \n")
   for p = 2:11
      h = 0.1^p
      dEh = copy(dE)
      for n = 1:length(dE)
         x[n] += h
         dEh[n] = (F(x) - E) / h
         x[n] -= h
      end
      push!(errors, norm(dE - dEh, Inf))
      verbose && @printf(" %1.1e | %4.2e  \n", h, errors[end])
   end
   verbose && @printf("---------|----------- \n")
   if minimum(errors) <= 1e-3 * maximum(errors)
      verbose && println("passed")
      return true
   else
      @warn("""It seems the finite-difference test has failed, which indicates
      that there is an inconsistency between the function and gradient
      evaluation. Please double-check this manually / visually. (It is
      also possible that the function being tested is poorly scaled.)""")
      return false
   end
end

dirfdtest(F, dF, x, u; kwargs...) =
      fdtest(t -> F(x + t * u),
             t -> dF(x + t * u) .* Ref(u),
             0.0; kwargs...)

end


text:
# E. Dynamic Parallelism

## error handling

# TODO: generalize for all uses of cudadevrt

struct CuDeviceError <: Exception
    code::cudaError
end

Base.convert(::Type{cudaError_t}, err::CuDeviceError) = err.code

name(err::CuDeviceError) = cudaGetErrorName(err)

description(err::CuDeviceError) = cudaGetErrorString(err)

@noinline function throw_device_cuerror(err::CuDeviceError)
    # the exception won't be rendered on the host, so print some details here already
    @cuprintln("ERROR: a CUDA error was thrown during kernel execution: $(description(err)) (code $(Int(err.code)), $(name(err)))")
    throw(err)
end

macro check_status(ex)
    quote
        res = $(esc(ex))
        if res != cudaSuccess
            throw_device_cuerror(CuDeviceError(res))
        end
    end
end

## streams

export CuDeviceStream

struct CuDeviceStream
    handle::cudaStream_t

    function CuDeviceStream(flags=cudaStreamNonBlocking)
        handle_ref = Ref{cudaStream_t}()
        @check_status cudaStreamCreateWithFlags(handle_ref, flags)
        return new(handle_ref[])
    end
end

Base.unsafe_convert(::Type{cudaStream_t}, s::CuDeviceStream) = s.handle

function unsafe_destroy!(s::CuDeviceStream)
    @check_status cudaStreamDestroy(s)
    return
end

## execution

# device-side counterpart of launch
@inline function device_launch(f, blocks, threads, shmem, stream, args...)
    blockdim = CuDim3(blocks)
    threaddim = CuDim3(threads)

    buf = parameter_buffer(f, blockdim, threaddim, shmem, args...)
    @check_status cudaLaunchDeviceV2(buf, stream)

    return
end

@generated function parameter_buffer(f, blocks, threads, shmem, args...)
    # allocate a buffer
    ex = quote
        Base.@_inline_meta
        buf = cudaGetParameterBufferV2(f, blocks, threads, shmem)
        ptr = Base.unsafe_convert(Ptr{UInt32}, buf)
    end

    # store the parameters
    #
    # D.3.2.2. Parameter Buffer Layout
    # > Each individual parameter placed in the parameter buffer is required to be aligned.
    # > That is, each parameter must be placed at the n-th byte in the parameter buffer,
    # > where n is the smallest multiple of the parameter size that is greater than the
    # > offset of the last byte taken by the preceding parameter. The maximum size of the
    # > parameter buffer is 4KB.
    #
    # NOTE: the above seems wrong, and we should use the parameter alignment, not its size.
    last_offset = 0
    for i in 1:length(args)
        T = args[i]
        align = Base.datatype_alignment(T)
        offset = Base.cld(last_offset, align) * align
        push!(ex.args, :(
            Base.pointerset(convert(Ptr{$T}, ptr+$offset), args[$i], 1, $align)
        ))
        last_offset = offset + sizeof(T)
    end

    push!(ex.args, :(return buf))

    return ex
end

## synchronization

"""
    device_synchronize()

Wait for the device to finish. This is the device side version,
and should not be called from the host.

`device_synchronize` acts as a synchronization point for
child grids in the context of dynamic parallelism.
"""
device_synchronize() = @check_status cudaDeviceSynchronize()


text:
using LinearAlgebra, SparseArrays
import SuiteSparse

function Successive_solves!( x, Kc, b, nsolves )
    for i=1:nsolves
        x  .= Kc\b
    end
    return
end

function ThreePointStencil( N, dx )
    # Assemble matrix (not the part of interest)
    cC =  2.0/dx^2 * ones(Float64, N); cC[1]   = 3.0/dx^2; cC[end]   = 3.0/dx^2
    cW = -1.0/dx^2 * ones(Float64, N); cW[1]   = 0.0
    cE = -1.0/dx^2 * ones(Float64, N); cE[end] = 0.0
    iC = 1:N
    iW = ones(Int64, N); iW[2:end-0] .= iC[1:end-1]
    iE = ones(Int64, N); iE[1:end-1] .= iC[2:end-0]
    I  = [iC[:]; iC[:]; iC[:]]
    J  = [iC[:]; iW[:]; iE[:]]
    V  = [cC[:]; cW[:]; cE[:]]
    K  = sparse( I, J, V )
    droptol!(K, 1e-9)
    return K
end

# Input
N  = 1000000
dx = 1

# Sparse matrix
K = ThreePointStencil( N, dx )

# Allocate arrays
b = rand(N)           # Right-hand-side array
x = zeros(Float64, N) # Solution array

# Factor
Kc = cholesky(K)

for i=1:3
@time @views x  .= Kc\b
@time @views x  .= SuiteSparse.CHOLMOD.solve(0, Kc, SuiteSparse.CHOLMOD.Dense(b) )[:]
end
xcm = SuiteSparse.CHOLMOD.zeros(N,1)
bcm = SuiteSparse.CHOLMOD.Dense(b)

for i=1:3
    @time @views x  .= Kc\b
    xcm  = SuiteSparse.CHOLMOD.solve(0, Kc, bcm )
end

# SuiteSparse.CHOLMOD.spsolve(Kc,SparseArrays.CHOLMOD.Sparse(b))
# ldiv!(Kc, b)
@time x  .= lu(K)\b

# # One solve
# @time x  .= Kc\b  # why does this allocate memory although the factorisation done and the vectors (x, b) are allocated?

# # Successive solves
# nsolves = 10
# @time Successive_solves!( x, Kc, b, nsolves ) # then, memory allocation depends on the number of solves...


text:
"""
For an initial guess from the index of a Hopf bifurcation point located in ContResult.specialpoint, returns a point which will be refined using `newtonHopf`.
"""
function HopfPoint(br::AbstractBranchResult, index::Int)
	@assert br.specialpoint[index].type == :hopf "The provided index does not refer to a Hopf point"
	specialpoint = br.specialpoint[index]			# Hopf point
	eigRes = br.eig									# eigenvector at the Hopf point
	p = specialpoint.param							# parameter value at the Hopf point
	ω = imag(eigRes[specialpoint.idx].eigenvals[specialpoint.ind_ev])	# frequency at the Hopf point
	return BorderedArray(specialpoint.x, [p, ω] )
end
####################################################################################################
# constructors
HopfProblemMinimallyAugmented(F, J, Ja, d2F, lens::Lens, a, b, linsolve::AbstractLinearSolver, linbdsolve = MatrixBLS()) = HopfProblemMinimallyAugmented(F, J, Ja, d2F, lens, a, b, 0*a, linsolve, linsolve, linbdsolve, linbdsolve, false)

HopfProblemMinimallyAugmented(F, J, Ja, lens::Lens, a, b, linsolve::AbstractLinearSolver,  linbdsolver = MatrixBLS()) = HopfProblemMinimallyAugmented(F, J, Ja, nothing, lens, a, b, linsolve)
####################################################################################################
getVec(x::BorderedArray, ::HopfProblemMinimallyAugmented) = x.u
getVec(x::AbstractVector, ::HopfProblemMinimallyAugmented) = @view x[1:end-2]
getP(x::BorderedArray, ::HopfProblemMinimallyAugmented) = x.p
getP(x::AbstractVector, ::HopfProblemMinimallyAugmented) = @view x[end-1:end]

# this function encodes the functional
function (hp::HopfProblemMinimallyAugmented)(x, p::T, ω::T, params) where T
	# These are the equations of the minimally augmented (MA) formulation of the Hopf bifurcation point
	# input:
	# - x guess for the point at which the jacobian has a purely imaginary eigenvalue
	# - p guess for the parameter for which the jacobian has a purely imaginary eigenvalue
	# The jacobian of the MA problem is solved with a BLS method
	a = hp.a
	b = hp.b
	# update parameter
	par = set(params, hp.lens, p)
	# ┌         ┐┌  ┐ ┌ ┐
	# │ J-iω  a ││v │=│0│
	# │  b    0 ││σ1│ │1│
	# └         ┘└  ┘ └ ┘
	# In the notations of Govaerts 2000, a = w, b = v
	# Thus, b should be a null vector of J - iω
	#       a should be a null vector of J'+ iω
	# we solve (J - iω)⋅v + a σ1 = 0 with <b, v> = n
	n = T(1)
	# note that the shift argument only affect J in this call:
	σ1 = hp.linbdsolver(hp.J(x, par), a, b, T(0), hp.zero, n; shift = Complex{T}(0, -ω))[2]

	# we solve (J+iω)'w + b σ2 = 0 with <a, w> = n
	# we find sigma2 = conj(sigma1)
	# w, σ2, _ = fp.linbdsolver(fp.Jadjoint(x, p) - Complex(0, ω) * I, b, a, 0., zeros(N), n)

	# the constraint is σ = <w, Jv> / n
	# σ = -dot(w, apply(fp.J(x, p) + Complex(0, ω) * I, v)) / n
	# we should have σ = σ1

	return hp.F(x, par), real(σ1), imag(σ1)
end

# this function encodes the functional
function (hopfpb::HopfProblemMinimallyAugmented)(x::BorderedArray, params)
	res = hopfpb(x.u, x.p[1], x.p[2], params)
	return BorderedArray(res[1], [res[2], res[3]])
end

# Struct to invert the jacobian of the Hopf MA problem.
struct HopfLinearSolverMinAug <: AbstractLinearSolver; end

"""
This function solves the linear problem associated with a linearization of the minimally augmented formulation of the Hopf bifurcation point. The keyword `debugArray` is used to debug the routine by returning several key quantities.
"""
function hopfMALinearSolver(x, p::T, ω::T, pb::HopfProblemMinimallyAugmented, par,
	 						duu, dup, duω;
							debugArray = nothing) where T
	################################################################################################
	# debugArray is used as a temp to be filled with values used for debugging. If debugArray = nothing, then no debugging mode is entered. If it is AbstractVector, then it is used
	################################################################################################
	# N = length(du) - 2
	# The jacobian should be passed as a tuple as Jac_hopf_MA(u0, pb::HopfProblemMinimallyAugmented) = (return (u0, pb, d2F::Bool))
	# The Jacobian J of the vector field is expressed at (x, p)
	# the jacobian expression of the hopf problem Jhopf is
	#           ┌             ┐
	#  Jhopf =  │  J  dpF   0 │
	#           │ σx   σp  σω │
	#           └             ┘
	########## Resolution of the bordered linear system ########
	# J * dX	  + dpF * dp		   = du => dX = x1 - dp * x2
	# The second equation
	#	<σx, dX> +  σp * dp + σω * dω = du[end-1:end]
	# thus becomes
	#   (σp - <σx, x2>) * dp + σω * dω = du[end-1:end] - <σx, x1>
	# This 2 x 2 system is then solved to get (dp, dω)
	############### Extraction of function names #################
	Fhandle = pb.F
	J = pb.J

	d2F = pb.d2F
	a = pb.a
	b = pb.b

	# parameter axis
	lens = pb.lens
	# update parameter
	par0 = set(par, lens, p)

	# we define the following jacobian. It is used at least 3 times below. This avoid doing 3 times the possibly costly building of J(x, p)
	J_at_xp = J(x, par0)

	# we do the following to avoid computing J_at_xp twice in case pb.Jadjoint is not provided
	JAd_at_xp = hasAdjoint(pb) ? pb.Jᵗ(x, par0) : transpose(J_at_xp)

	# normalization
	n = T(1)

	# we solve (J-iω)v + a σ1 = 0 with <b, v> = n
	v, σ1, _, itv = pb.linbdsolver(J_at_xp, a, b, T(0), pb.zero, n; shift = Complex{T}(0, -ω))

	# we solve (J+iω)'w + b σ1 = 0 with <a, w> = n
	w, σ2, _, itw = pb.linbdsolverAdjoint(JAd_at_xp, b, a, T(0), pb.zero, n; shift = Complex{T}(0, ω))

	δ = T(1e-9)
	ϵ1, ϵ2, ϵ3 = T(δ), T(δ), T(δ)
	################### computation of σx σp ####################
	################### and inversion of Jhopf ####################
	dpF   = (Fhandle(x, set(par, lens, p + ϵ1))	 - Fhandle(x, set(par, lens, p - ϵ1))) / T(2ϵ1)
	dJvdp = (apply(J(x, set(par, lens, p + ϵ3)), v) - apply(J(x, set(par, lens, p - ϵ3)), v)) / T(2ϵ3)
	σp = -dot(w, dJvdp) / n

	# case of sigma_omega
	# σω = dot(w, Complex{T}(0, 1) * v) / n
	σω = Complex{T}(0, 1) * dot(w, v) / n

	# we solve J⋅x1 = duu and J⋅x2 = dpF
	x1, x2, _, (it1, it2) = pb.linsolver(J_at_xp, duu, dpF)

	# the case of ∂_xσ is a bit more involved
	# we first need to compute the value of ∂_xσ written σx
	# σx = zeros(Complex{T}, length(x))
	σx = similar(x, Complex{T})

	if hasHessian(pb) == false
		cw = conj(w)
		vr = real(v); vi = imag(v)
		u1r = applyJacobian(pb, x + ϵ2 * vr, par0, cw, true)
		u1i = applyJacobian(pb, x + ϵ2 * vi, par0, cw, true)
		u2 = apply(JAd_at_xp,  cw)
		σxv2r = @. -(u1r-u2) / ϵ2
		σxv2i = @. -(u1i-u2) / ϵ2
		σx = @. σxv2r + Complex{T}(0, 1) * σxv2i

		σxx1 = dot(σx, x1)
		σxx2 = dot(σx, x2)
	else
		d2Fv = d2F(x, par0, v, x1)
		σxx1 = -conj(dot(w, d2Fv) / n)
		d2Fv = d2F(x, par0, v, x2)
		σxx2 = -conj(dot(w, d2Fv) / n)
	end
	# we need to be carefull here because the dot produces conjugates. Hence the + dot(σx, x2) and + imag(dot(σx, x1) and not the opposite
	dp, dω = [real(σp - σxx2) real(σω);
			  imag(σp + σxx2) imag(σω) ] \
			  [dup - real(σxx1), duω + imag(σxx1)]

	if debugArray isa AbstractVector
		debugArray .= vcat(σp, σω, σx)
	end
	return x1 .- dp .* x2, dp, dω, true, it1 + it2 + sum(itv) + sum(itw)
end

function (hopfl::HopfLinearSolverMinAug)(Jhopf, du::BorderedArray{vectype, T}; debugArray = nothing, kwargs...)  where {vectype, T}
	# kwargs is used by AbstractLinearSolver
	out = hopfMALinearSolver((Jhopf.x).u,
				(Jhopf.x).p[1],
				(Jhopf.x).p[2],
				Jhopf.hopfpb,
				Jhopf.params,
				du.u, du.p[1], du.p[2];
				debugArray = debugArray)
	# this type annotation enforces type stability
	BorderedArray{vectype, T}(out[1], [out[2], out[3]]), out[4], out[5]
end

################################################################################################### Newton / Continuation functions
"""
$(SIGNATURES)

This function turns an initial guess for a Hopf point into a solution to the Hopf problem based on a Minimally Augmented formulation. The arguments are as follows
- `F   = (x, p) -> F(x, p)` where `p` is a set of parameters.
- `dF  = (x, p) -> d_xF(x, p)` associated jacobian
- `hopfpointguess` initial guess (x_0, p_0) for the Hopf point. It should a `BorderedArray` as returned by the function `HopfPoint`.
- `par` parameters used for the vector field
- `lens` parameter axis used to locate the Hopf point.
- `eigenvec` guess for the  iω eigenvector
- `eigenvec_ad` guess for the -iω eigenvector
- `options::NewtonPar` options for the Newton-Krylov algorithm, see [`NewtonPar`](@ref).

# Optional arguments:
- `Jᵗ = (x, p) -> transpose(d_xF(x, p))` jacobian adjoint, it should be implemented in an efficient manner. For matrix-free methods, `transpose` is not readily available and the user must provide a dedicated method. In the case of sparse based jacobian, `Jᵗ` should not be passed as it is computed internally more efficiently, i.e. it avoids recomputing the jacobian as it would be if you pass `Jᵗ = (x, p) -> transpose(dF(x, p))`
- `d2F = (x, p, v1, v2) ->  d2F(x, p, v1, v2)` a bilinear operator representing the hessian of `F`. It has to provide an expression for `d2F(x, p)[v1, v2]`.
- `normN = norm`
- `bdlinsolver` bordered linear solver for the constraint equation
- `kwargs` keywords arguments to be passed to the regular Newton-Krylov solver

# Simplified call:
Simplified call to refine an initial guess for a Hopf point. More precisely, the call is as follows

	newtonHopf(F, J, br::AbstractBranchResult, ind_hopf::Int; Jᵗ = nothing, d2F = nothing, normN = norm, options = br.contparams.newtonOptions, kwargs...)

where the optional argument `Jᵗ` is the jacobian transpose and the Hessian is `d2F`. The parameters / options are as usual except that you have to pass the branch `br` from the result of a call to `continuation` with detection of bifurcations enabled and `index` is the index of bifurcation point in `br` you want to refine. You can pass newton parameters different from the ones stored in `br` by using the argument `options`.

!!! tip "Jacobian tranpose"
    The adjoint of the jacobian `J` is computed internally when `Jᵗ = nothing` by using `transpose(J)` which works fine when `J` is an `AbstractArray`. In this case, do not pass the jacobian adjoint like `Jᵗ = (x, p) -> transpose(d_xF(x, p))` otherwise the jacobian will be computed twice!

!!! tip "ODE problems"
    For ODE problems, it is more efficient to pass the Bordered Linear Solver using the option `bdlinsolver = MatrixBLS()`
"""
function newtonHopf(F, J,
			hopfpointguess::BorderedArray,
			par, lens::Lens,
			eigenvec, eigenvec_ad,
			options::NewtonPar;
			Jᵗ = nothing,
			d2F = nothing,
			normN = norm,
			bdlinsolver::AbstractBorderedLinearSolver = MatrixBLS(),
			kwargs...)
	hopfproblem = HopfProblemMinimallyAugmented(
		F, J, Jᵗ, d2F, lens,
		_copy(eigenvec_ad),	# this is pb.a ≈ null space of (J - iω I)^*
		_copy(eigenvec), 	# this is pb.b ≈ null space of  J - iω I
		options.linsolver,
		# do not change linear solver if user provides it
		@set bdlinsolver.solver = (isnothing(bdlinsolver.solver) ? options.linsolver : bdlinsolver.solver))

	# Jacobian for the Hopf problem
	Jac_hopf_MA = (x, param) -> (x = x, params = param, hopfpb = hopfproblem)

	# options for the Newton Solver
	opt_hopf = @set options.linsolver = HopfLinearSolverMinAug()

	# solve the hopf equations
	return newton(hopfproblem, Jac_hopf_MA, hopfpointguess, par, opt_hopf, normN = normN, kwargs...)..., hopfproblem
end

function newtonHopf(F, J,
			br::AbstractBranchResult, ind_hopf::Int;
			Jᵗ = nothing,
			d2F = nothing,
			normN = norm,
			options = br.contparams.newtonOptions,
			verbose = true,
			nev = br.contparams.nev,
			startWithEigen = false,
			kwargs...)
	hopfpointguess = HopfPoint(br, ind_hopf)
	ω = hopfpointguess.p[2]
	bifpt = br.specialpoint[ind_hopf]
	options.verbose && println("--> Newton Hopf, the eigenvalue considered here is ", br.eig[bifpt.idx].eigenvals[bifpt.ind_ev])
	@assert bifpt.idx == bifpt.step + 1 "Error, the bifurcation index does not refer to the correct step"
	ζ = geteigenvector(options.eigsolver, br.eig[bifpt.idx].eigenvec, bifpt.ind_ev)
	ζ ./= normN(ζ)
	ζad = LinearAlgebra.conj.(ζ)

	if startWithEigen
		# computation of adjoint eigenvalue. Recall that b should be a null vector of J-iω
		λ = Complex(0, ω)
		p = bifpt.param
		parbif = setParam(br, p)

		# jacobian at bifurcation point
		L = J(bifpt.x, parbif)

		# computation of adjoint eigenvector
		_Jt = isnothing(Jᵗ) ? adjoint(L) : Jᵗ(bifpt.x, parbif)
		ζstar, λstar = getAdjointBasis(_Jt, conj(λ), options.eigsolver; nev = nev, verbose = false)
		ζad .= ζstar ./ dot(ζstar, ζ)
	end

	# solve the hopf equations
	return newtonHopf(F, J, hopfpointguess, br.params, br.lens, ζ, ζad, options; Jᵗ = Jᵗ, d2F = d2F, normN = normN, kwargs...)
end

"""
$(SIGNATURES)

codim 2 continuation of Hopf points. This function turns an initial guess for a Hopf point into a curve of Hopf points based on a Minimally Augmented formulation. The arguments are as follows
- `F = (x, p) ->	F(x, p)` where `p` is a set of parameters
- `J = (x, p) -> d_xF(x, p)` associated jacobian
- `hopfpointguess` initial guess (x_0, p1_0) for the Hopf point. It should be a `Vector` or a `BorderedArray`
- `par` set of parameters
- `lens1` parameter axis for parameter 1
- `lens2` parameter axis for parameter 2
- `eigenvec` guess for the iω eigenvector at p1_0
- `eigenvec_ad` guess for the -iω eigenvector at p1_0
- `options_cont` keywords arguments to be passed to the regular [`continuation`](@ref)

# Optional arguments:

- `Jᵗ = (x, p) -> adjoint(d_xF(x, p))` associated jacobian adjoint
- `d2F = (x, p, v1, v2) -> d2F(x, p, v1, v2)` this is the hessian of `F` computed at `(x, p)` and evaluated at `(v1, v2)`. This helps solving the linear problem associated to the Hopf minimally augmented formulation.
- `d3F = (x, p, v1, v2, v3) -> d3F(x, p, v1, v2, v3)` this is the third derivative of `F` computed at `(x, p)` and evaluated at `(v1, v2, v3)`. This is used to detect **Bautin** bifurcation.
- `bdlinsolver` bordered linear solver for the constraint equation
- `updateMinAugEveryStep` update vectors `a,b` in Minimally Formulation every `updateMinAugEveryStep` steps
- `computeEigenElements = false` whether to compute eigenelements. If `options_cont.detecttEvent>0`, it allows the detection of ZH, HH points.
- `kwargs` keywords arguments to be passed to the regular [`continuation`](@ref)

# Simplified call:

	continuationHopf(F, J, br::AbstractBranchResult, ind_hopf::Int, lens2::Lens, options_cont::ContinuationPar ;  kwargs...)

where the parameters are as above except that you have to pass the branch `br` from the result of a call to `continuation` with detection of bifurcations enabled and `index` is the index of Hopf point in `br` that you want to refine.

!!! tip "ODE problems"
    For ODE problems, it is more efficient to pass the Bordered Linear Solver using the option `bdlinsolver = MatrixBLS()`

!!! tip "Jacobian transpose"
    The adjoint of the jacobian `J` is computed internally when `Jᵗ = nothing` by using `transpose(J)` which works fine when `J` is an `AbstractArray`. In this case, do not pass the jacobian adjoint like `Jᵗ = (x, p) -> transpose(d_xF(x, p))` otherwise the jacobian would be computed twice!

!!! tip "Detection of Bogdanov-Takens and Bautin bifurcations"
    In order to trigger the detection, pass `detectEvent = 1,2` in `options_cont`. Note that you need to provide `d3F`.
"""
function continuationHopf(F, J,
				hopfpointguess::BorderedArray{vectype, Tb}, par,
				lens1::Lens, lens2::Lens,
				eigenvec, eigenvec_ad,
				options_cont::ContinuationPar ;
				updateMinAugEveryStep = 0,
				normC = norm,
				Jᵗ = nothing,
				d2F = nothing,
				d3F = nothing,
				bdlinsolver::AbstractBorderedLinearSolver = MatrixBLS(),
				computeEigenElements = false,
				kwargs...) where {Tb, vectype}
	@assert lens1 != lens2 "Please choose 2 different parameters. You only passed $lens1"

	# we need d2F and d3F to detect Bautin bifurcation
	if options_cont.detectEvent > 0
		@assert ~isnothing(d2F) && ~isnothing(d3F) "Please provide d2F and d3F"
	end

	# options for the Newton Solver inheritated from the ones the user provided
	options_newton = options_cont.newtonOptions

	hopfPb = HopfProblemMinimallyAugmented(
		F, J, Jᵗ, d2F,
		lens1,
		_copy(eigenvec_ad),	# this is pb.a ≈ null space of (J - iω I)^*
		_copy(eigenvec), 	# this is pb.b ≈ null space of  J - iω I
		options_newton.linsolver,
		# do not change linear solver if user provides it
		@set bdlinsolver.solver = (isnothing(bdlinsolver.solver) ? options_newton.linsolver : bdlinsolver.solver))

	# Jacobian for the Hopf problem
	Jac_hopf_MA = (x, param) -> (x = x, params = param, hopfpb = hopfPb)

	opt_hopf_cont = @set options_cont.newtonOptions.linsolver = HopfLinearSolverMinAug()

	# this functions allows to tackle the case where the two parameters have the same name
	lenses = getLensSymbol(lens1, lens2)

	# current lyapunov coefficient
	l1 = Complex{eltype(Tb)}(0, 0)
	BT = one(eltype(Tb))
	GH = one(eltype(Tb))

	# this function is used as a Finalizer
	# it is called to update the Minimally Augmented problem
	# by updating the vectors a, b
	function updateMinAugHopf(z, tau, step, contResult; kUP...)
		x = z.u.u		# hopf point
		p1 = z.u.p[1]	# first parameter
		ω = z.u.p[2]	# Hopf frequency
		# we first check that the continuation step was successful
		# if not, we do not update the problem with bad information!
		success = get(kUP, :state, nothing).isconverged
		(~modCounter(step, updateMinAugEveryStep) || success == false) && return true
		p2 = z.p		# second parameter
		newpar = set(par, lens1, p1)
		newpar = set(newpar, lens2, p2)

		a = hopfPb.a
		b = hopfPb.b

		# expression of the jacobian
		J_at_xp = hopfPb.J(x, newpar)

		# compute new b
		T = typeof(p1)
		n = T(1)
		newb = hopfPb.linbdsolver(J_at_xp, a, b, T(0), hopfPb.zero, n; shift = Complex(0, -ω))[1]

		# compute new a
		JAd_at_xp = hasAdjoint(hopfPb) ? hopfPb.Jᵗ(x, newpar) : transpose(J_at_xp)
		newa = hopfPb.linbdsolver(JAd_at_xp, b, a, T(0), hopfPb.zero, n; shift = Complex(0, ω))[1]

		hopfPb.a .= newa ./ normC(newa)
		# do not normalize with dot(newb, hopfPb.a), it prevents BT detection
		hopfPb.b .= newb ./ normC(newb)

		# we stop continuation at Bogdanov-Takens points
		isbt = isnothing(contResult) ? true : isnothing(findfirst(x->x.type == :bt, contResult.specialpoint))

		# if the frequency is null, this is not a Hopf point, we halt the process
		threshBT = 100options_newton.tol
		if abs(ω) < threshBT
			@warn "[Codim 2 Hopf - Finalizer] The Hopf curve seems to be close to a BT point: ω ≈ $ω. Stopping computations at ($p1, $p2). If the BT point is not detected, trying lowering Newton tolerance or dsmax."
		end

		# call the user-passed finalizer
		finaliseUser = get(kwargs, :finaliseSolution, nothing)
		resFinal = isnothing(finaliseUser) ? true : finaliseUser(z, tau, step, contResult; prob = hopfPb, kUP...)
		return abs(ω) >= threshBT && isbt && resFinal
	end

	function testGH_BT(iter, state)
		z = getx(state)
		x = z.u					# hopf point
		p1 = z.p[1]				# first parameter
		ω  = z.p[2]				# Hopf frequency
		p2 = getp(state)		# second parameter
		newpar = set(par, lens1, p1)
		newpar = set(newpar, lens2, p2)

		a = hopfPb.a
		b = hopfPb.b

		# expression of the jacobian
		J_at_xp = hopfPb.J(x, newpar)

		# compute new b
		T = typeof(p1)
		n = T(1)
		ζ = hopfPb.linbdsolver(J_at_xp, a, b, T(0), hopfPb.zero, n; shift = Complex(0, -ω))[1]
		ζ ./= normC(ζ)

		# compute new a
		JAd_at_xp = hasAdjoint(hopfPb) ? hopfPb.Jᵗ(x, newpar) : transpose(J_at_xp)
		ζstar = hopfPb.linbdsolver(JAd_at_xp, b, a, T(0), hopfPb.zero, n; shift = Complex(0, ω))[1]
		# test function for Bogdanov-Takens
		BT = ω#real( dot(ζstar ./ normC(ζstar), ζ) )
		ζstar ./= dot(ζ, ζstar)

		hp = Hopf(x, p1, ω, newpar, lens1, ζ, ζstar, (a=Complex{T}(0, 0), b = Complex{T}(0,0)), :hopf)
		hopfNormalForm(F, J, d2F, d3F, hp, options_newton.linsolver, verbose=false)

		# lyapunov coefficient
		l1 = hp.nf.b
		# test for Bautin bifurcation.
		# If GH is too large, we take the previous value to avoid spurious detection
		# GH will be large close to BR points
		GH = abs(real(hp.nf.b))<1e10 ? real(hp.nf.b) : state.eventValue[2][1]
		return GH, BT
	end

	# the following allows to append information specific to the codim 2 continuation to the user data
	_printsol = get(kwargs, :recordFromSolution, nothing)
	_printsol2 = isnothing(_printsol) ?
		(u, p; kw...) -> (zip(lenses, (u.p[1], p))..., ω = u.p[2], l1=l1, BT = BT, GH = GH) :
		(u, p; kw...) -> begin
			(namedprintsol(_printsol(u, p;kw...))..., zip(lenses, (u.p[1], p))..., ω = u.p[2], l1 = l1, BT = BT, GH = GH)
		end

	# eigen solver
	eigsolver = HopfEig(getsolver(opt_hopf_cont.newtonOptions.eigsolver))

	# event for detecting codim 2 points
	if computeEigenElements
		event = PairOfEvents(ContinuousEvent(2, testGH_BT, computeEigenElements, ("gh", "bt"), 0), BifDetectEvent)
	else
		event = ContinuousEvent(2, testGH_BT, false, ("gh", "bt"), 0)
	end

	# solve the hopf equations
	branch, u, tau = continuation(
		hopfPb, Jac_hopf_MA,
		hopfpointguess, par, lens2,
		(@set opt_hopf_cont.newtonOptions.eigsolver = eigsolver);
		kwargs...,
		linearAlgo = BorderingBLS(solver = opt_hopf_cont.newtonOptions.linsolver, checkPrecision = false),
		normC = normC,
		recordFromSolution = _printsol2,
		finaliseSolution = updateMinAugHopf,
		event = event
	)
	@assert ~isnothing(branch) "Empty branch!"
	return correctBifurcation(setproperties(branch; type = :HopfCodim2, functional = hopfPb)), u, tau
end

function continuationHopf(F, J,
						br::AbstractBranchResult, ind_hopf::Int64,
						lens2::Lens, options_cont::ContinuationPar ;
						startWithEigen = false,
						normC = norm,
						Jᵗ = nothing,
						d2F = nothing,
						d3F = nothing,
						kwargs...)
	hopfpointguess = HopfPoint(br, ind_hopf)
	ω = hopfpointguess.p[2]
	bifpt = br.specialpoint[ind_hopf]

	@assert ~isnothing(br.eig) "The branch contains no eigen elements. This is strange because a Hopf point was detected. Please open an issue on the website."

	@assert ~isnothing(br.eig[1].eigenvec) "The branch contains no eigenvectors for the Hopf point. Please provide one."

	ζ = geteigenvector(options_cont.newtonOptions.eigsolver, br.eig[bifpt.idx].eigenvec, bifpt.ind_ev)
	ζ ./= normC(ζ)
	ζad = conj.(ζ)

	p = bifpt.param
	parbif = setParam(br, p)

	if startWithEigen
		# computation of adjoint eigenvalue
		λ = Complex(0, ω)
		# jacobian at bifurcation point
		L = J(bifpt.x, parbif)
		_Jt = isnothing(Jᵗ) ? adjoint(L) : Jᵗ(bifpt.x, parbif)
		ζstar, λstar = getAdjointBasis(_Jt, conj(λ), br.contparams.newtonOptions.eigsolver; nev = br.contparams.nev, verbose = false)
		ζad .= ζstar ./ dot(ζstar, ζ)
	end

	return continuationHopf(F, J, hopfpointguess, parbif, br.lens, lens2, ζ, ζad, options_cont ; Jᵗ = Jᵗ, d2F = d2F, d3F = d3F, normC = normC, kwargs...)
end

# structure to compute the eigenvalues along the Hopf branch
struct HopfEig{S} <: AbstractCodim2EigenSolver
	eigsolver::S
end

geteigenvector(eig::HopfEig, vectors, i::Int) = geteigenvector(eig.eigsolver, vectors, i)

function (eig::HopfEig)(Jma, nev; kwargs...)
	n = min(nev, length(Jma.x.u))

	x = Jma.x.u		# hopf point
	p1, ω = Jma.x.p	# first parameter
	newpar = set(Jma.params, getLens(Jma.hopfpb), p1)

	J = Jma.hopfpb.J(x, newpar)

	eigenelts = eig.eigsolver(J, n; kwargs...)
	# @assert 1==0 "Mettre les 2 params a jour"
	return eigenelts
end


text:
module GameZero

export
  # Actor
  Actor,
  angle,
  collide,
  distance,

  # Game
  Game,
  game,
  rungame,

  # Screen and drawing
  clear,
  draw,

  # Audio
  play_music,
  play_sound,

  # Timer
  scheduler,
  schedule_once,
  schedule_interval,
  schedule_unique,
  unschedule,

  # Keyboard modules
  Keys,
  MouseButtons,
  KeyMods,

  # Geometric structs
  Line,
  Rect,
  Triangle,
  Circle

using Random
import Base: +

using Colors
using SimpleDirectMediaLayer
const SDL2 = SimpleDirectMediaLayer

include("keyboard/Keys_MODULE.jl")
include("keyboard/Keymods_MODULE.jl")
include("keyboard/MouseButtons_MODULE.jl")
include("keyboard/Keyboard.jl")

include("timer/WallTimer.jl")
include("timer/Scheduled.jl")

include("window.jl")
include("event.jl")
include("resources.jl")
include("utility.jl")

include("Screen.jl")
include("Geom.jl")
include("Actor.jl")

include("Game.jl")

end