text:
using Base.Test, ODBC, DataStreams, Missings, WeakRefStrings, DataFrames, DecFP

@show ODBC.drivers()
@show ODBC.dsns()

@show ODBC.API.odbc_dm

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

run(`uname -a`)

dsn = if haskey(ENV, "ODBC_DSN")
    ODBC.DSN(ENV["ODBC_DSN"])
else
    props = open(f->Dict(split(ln, ':') for ln in eachline(f)) , "/Users/jacobquinn/.tug/application.properties")
    RIG, server, port = map(strip, split(props["environment.subdomain"], '.'))[1], strip(props["routing.services.mysql.host"]), strip(props["routing.services.mysql.port"])
    ODBC.DSN("Driver={MySQL ODBC Driver}; user=domo; password=popchart; server=$server; port=$port;")
end

# dsn = ODBC.DSN("Driver=MySQL;uid=root")
# dsn = ODBC.DSN("MySQL-test", "root", "")

# 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
println("testing all mysql types...")
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 = ODBC.query(source)
@test size(Data.schema(data)) == (1,27)
@test Data.types(Data.schema(data)) == (
 Union{Int64, Missing},
 Union{Int8, Missing},
 Union{DecFP.Dec64, Missing},
 Union{Int32, Missing},
 Union{DecFP.Dec64, Missing},
 Union{Int16, Missing},
 Union{Int32, Missing},
 Union{Int8, Missing},
 Union{Float32, Missing},
 Union{Float64, Missing},
 Union{ODBC.API.SQLDate, Missing},
 Union{ODBC.API.SQLTimestamp, Missing},
 Union{ODBC.API.SQLTimestamp, Missing},
 Union{ODBC.API.SQLTime, Missing},
 Union{Int16, Missing},
 Union{WeakRefString{UInt8}, Missing},
 Union{WeakRefString{UInt8}, Missing},
 Union{Array{UInt8,1}, Missing},
 Union{Array{UInt8,1}, Missing},
 Union{Array{UInt8,1}, Missing},
 Union{Array{UInt8,1}, Missing},
 Union{Array{UInt8,1}, Missing},
 Union{Array{UInt8,1}, Missing},
 Union{String, Missing},
 Union{String, Missing},
 Union{String, Missing},
 Union{String, Missing})
@test data[1][1] === Int64(1)
@test data[2][1] === Int8(1)
@test data[3][1] === DecFP.Dec64(1)
@test data[4][1] === Int32(1)
@test data[5][1] === DecFP.Dec64(1)
@test data[6][1] === Int16(1)
@test data[7][1] === Int32(1)
@test data[8][1] === Int8(1)
@test data[9][1] === Float32(1.2)
@test data[10][1] === Float64(1.2)
@test data[11][1] === ODBC.API.SQLDate(2016,1,1)
@test data[12][1] === ODBC.API.SQLTimestamp(2016,1,1,1,1,1,0)
@test data[13][1] === ODBC.API.SQLTimestamp(2016,1,1,1,1,1,0)
@test data[14][1] === ODBC.API.SQLTime(1,1,1)
@test data[15][1] === Int16(2016)
@test string(data[16][1]) == "A"
@test string(data[17][1]) == "hey there sailor"
@test data[18][1] == UInt8[0x31,0x32]
@test ismissing(data[19][1])
@test data[20][1] == UInt8[0x68,0x65,0x79,0x20,0x74,0x68,0x65,0x72,0x65,0x20,0x61,0x62,0x72,0x61,0x68,0x61,0x6d]
@test data[21][1] == UInt8[0x68,0x65,0x79,0x20,0x74,0x68,0x65,0x72,0x65,0x20,0x62,0x69,0x6c,0x6c]
@test data[22][1] == UInt8[0x68,0x65,0x79,0x20,0x74,0x68,0x65,0x72,0x65,0x20,0x63,0x68,0x61,0x72,0x6c,0x69,0x65]
@test data[23][1] == UInt8[0x68,0x65,0x79,0x20,0x74,0x68,0x65,0x72,0x65,0x20,0x64,0x61,0x6e]
@test string(data[24][1]) == "hey there ephraim"
@test string(data[25][1]) == "hey there frank"
@test string(data[26][1]) == "hey there george"
@test string(data[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.schema(data)) == (2,27)
println("passed. testing large query...")

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")
# randoms = joinpath(Pkg.dir("ODBC"), "test/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.schema(data)) == (1,1)
@test data[1][1] === 70000

df = ODBC.query(dsn, "select * from test2")
@test size(Data.schema(df)) == (70000,7)
@test df[1] == collect(1:70000)
@test df[end][1] === ODBC.API.SQLTimestamp(2002,1,17,21,32,0,0)
println("passed. testing prepared statement...")

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

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(Data.schema(df)) == (2,7)
@test df[1][end-1] == 101
@test df[1][end] == 102
@test df[2][end-1] == "Steve"
@test df[2][end] == "Dean"
@test df[3][end-1] == "McQueen"
@test df[3][end] == "Martin"
@test df[4][end-1] == DecFP.Dec64(1)
@test df[4][end] == DecFP.Dec64(2)
@test df[5][end-1] == 100.0
@test df[5][end] == 10.1
@test df[6][end-1] == ODBC.API.SQLDate(2016,1,1)
@test df[6][end] == ODBC.API.SQLDate(2016,1,2)
@test df[7][end-1] == ODBC.API.SQLTimestamp(2016,1,1,0,0,0,0)
@test df[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")

ODBC.disconnect!(dsn)
println("mysql passed.")

# PostgreSQL
# dsn = ODBC.DSN("PgSQL-test", "postgres", "")
# dbs = ODBC.query(dsn, "SELECT datname FROM pg_database WHERE datistemplate = false;")
# data = ODBC.query(dsn, "SELECT table_schema,table_name FROM information_schema.tables ORDER BY table_schema,table_name;")
# ODBC.execute!(dsn, "drop table if exists test1")
# ODBC.execute!(dsn, "create table test1
#                     (test_bigint bigint,
#                      test_decimal decimal,
#                      test_int integer,
#                      test_numeric numeric,
#                      test_smallint smallint,
#                      test_float real,
#                      test_real double precision,
#                      test_money money,
#                      test_date date,
#                      test_timestamp timestamp,
#                      test_time time,
#                      test_char char(1),
#                      test_varchar varchar(16),
#                      test_bytea bytea,
#                      test_boolean boolean,
#                      test_text text,
#                      test_array integer[]
#                     )")
# data = ODBC.query(dsn, "select * from information_schema where table_name = 'test1'")
# showall(data)
# ODBC.execute!(dsn, "insert into test1 VALUES
#                     (1, -- bigint,
#                      1.2, -- decimal,
#                      2, -- integer,
#                      1.4, -- numeric,
#                      3, -- smallint,
#                      1.6, -- real,
#                      1.8, -- double precision,
#                      2.0, -- money,
#                      '2016-01-01', -- date,
#                      '2016-01-01 01:01:01', -- timestamp,
#                      '01:01:01', -- time,
#                      'A', -- char(1),
#                      'hey there sailor', -- varchar(16),
#                      NULL, -- bytea,
#                      TRUE, -- boolean,
#                      'hey there abraham', -- text
#                      ARRAY[1, 2, 3] -- integer array
#                     )")
# source = ODBC.Source(dsn, "select * from test1")
# data = ODBC.query(source)
# @test size(Data.schema(data)) == (1,17)
# @test Data.types(data, Data.Field) == (
#   Union{Int64, Missing}
#  ,Union{DecFP.Dec64, Missing}
#  ,Union{Int32, Missing}
#  ,Union{DecFP.Dec64, Missing}
#  ,Union{Int16, Missing}
#  ,Union{Float32, Missing}
#  ,Union{Float64, Missing}
#  ,Union{Float64, Missing}
#  ,Union{ODBC.API.SQLDate, Missing}
#  ,Union{ODBC.API.SQLTimestamp, Missing}
#  ,Union{ODBC.API.SQLTime, Missing}
#  ,Union{WeakRefStrings.WeakRefString{UInt8}, Missing}
#  ,Union{WeakRefStrings.WeakRefString{UInt8}, Missing}
#  ,Union{Array{UInt8,1}, Missing}
#  ,Union{WeakRefStrings.WeakRefString{UInt8}, Missing}
#  ,Union{String, Missing}
#  ,Union{WeakRefStrings.WeakRefString{UInt8}, Missing})
# @test data.test_array[1] == "{1,2,3}"
# showall(data)

# println("postgresql passed.")


text:
using ELF2
using Test

function execute(cmd::Cmd)
  out = Pipe()
  err = Pipe()

  process = run(pipeline(ignorestatus(cmd), stdout=out, stderr=err))
  close(out.in)
  close(err.in)

  (
    stdout = String(read(out)),
    stderr = String(read(err)),
    code = process.exitcode
  )
end

function check_command(command::String)
    which = `which $(command)`
    println(which)
    try
        run(which)
    catch
        println("Cannot find $(command) skip this test")
        return false
    end
    return true
end

function test_command_itself(command::String)
    if !check_command(command)
        return
    end

    path, _, _ = execute(`which $(command)`)
    path = chop(path, tail=1)
    f = open(path, "r")
    elf = ELF2.read_elf(f)
    show(elf)
end

function test_cc(command::String)
    if !check_command(command)
        return
    end

    cmd = `$(command) main.c -o main.$(command).out`
    run(cmd)
    f = open("main.$(command).out", "r")
    elf = ELF2.read_elf(f)
    show(elf)
end

function test_cxx(command::String)
    if !check_command(command)
        return
    end

    cmd = `$(command) main.cc -o main.$(command).out`
    run(cmd)
    f = open("main.$(command).out", "r")
    elf = ELF2.read_elf(f)
    show(elf)
end

function test_magic()
    f = open("libmax.so", "r")

    elf = ELF2.read_elf(f)
    @test elf.ehdr.e_ident1 == ELF2.ELFMAG0
    @test elf.ehdr.e_ident2 == ELF2.ELFMAG1
    @test elf.ehdr.e_ident3 == ELF2.ELFMAG2
    @test elf.ehdr.e_ident4 == ELF2.ELFMAG3
    @test elf.ehdr.e_data == ELF2.ELFDATA2LSB
    @test elf.ehdr.e_fversion == ELF2.EV_CURRENT
    @test elf.ehdr.e_type == ELF2.ET_DYN

    show(elf)

    close(f)
end

function test_demangle()
    @test ELF2.demangle("_Z41__static_initialization_and_destruction_0ii") == "__static_initialization_and_destruction_0"
    @test ELF2.demangle("_ZNSolsEPFRSoS_E") == "std::basic_ostream<char, std::char_traits<char> >::operator<<(std::basic_ostream<char, std::char_traits<char> >& (*)(std::basic_ostream<char, std::char_traits<char> >&))"
    @test ELF2.demangle("_ZNSt8ios_base4InitC1Ev") == "std::ios_base::Init::Init()"
    @test ELF2.demangle("_ZNSt8ios_base4InitD1Ev") == "std::ios_base::Init::~Init()"
    @test ELF2.demangle("_ZSt4cout") == "_ZSt4cout"
    @test ELF2.demangle("_ZSt4endlIcSt11char_traitsIcEERSt13basic_ostreamIT_T0_ES6_") == "_ZSt4endlIcSt11char_traitsIcEERSt13basic_ostreamIT_T0_ES6_"
    @test ELF2.demangle("_ZStL8__ioinit") == "_ZStL8__ioinit"
    @test ELF2.demangle("_ZStlsISt11char_traitsIcEERSt13basic_ostreamIcT_ES5_PKc") == "_ZStlsISt11char_traitsIcEERSt13basic_ostreamIcT_ES5_PKc"
end

@testset "ELF2.jl" begin
    # test_magic()
    # for cc in ["gcc", "aarch64-linux-gnu-gcc", "clang"]
    #     test_cc(cc)
    # end
    # for cxx in ["g++", "aarch64-linux-gnu-g++", "clang++"]
    #     test_cxx(cxx)
    # end
    # test_command_itself("gcc")
    # test_command_itself("g++")
    # test_command_itself("clang")
    # test_command_itself("clang++")
    test_demangle()
end

ext: jl


text:
## changingnames.jl  methods for NamedArray that change some of the names

## (c) 2013 David A. van Leeuwen

## This code is licensed under the MIT license
## See the file LICENSE.md in this distribution

## sum etc: keep names in one dimension
# import Base.sum, Base.prod, Base.maximum, Base.minimum, Base.mean, Base.std
for f = (:sum, :prod, :maximum, :minimum, :mean, :std, :var)
    eval(Expr(:import, :Base, f))
    @eval function ($f)(a::NamedArray, d::Dims)
        s = ($f)(a.array, d)
        dicts = [issubset(i,d) ? OrderedDict(string($f,"(",a.dimnames[i],")") => 1) : a.dicts[i] for i=1:ndims(a)]
        NamedArray(s, tuple(dicts...), a.dimnames)
    end
    @eval ($f)(a::NamedArray, d::Int) = ($f)(a, (d,))
end

function fan{T,N}(f::Function, fname::AbstractString, a::NamedArray{T,N}, dim::Int)
    dimnames = Array(Any, N)
    for i=1:N
        if i==dim
            dimnames[i] = string(fname, "(", a.dimnames[dim], ")")
        else
            dimnames[i] = a.dimnames[i]
        end
    end
     NamedArray(f(a.array,dim), a.dicts, tuple(dimnames...))
end

## rename a dimension
for f in (:cumprod, :cumsum, :cumsum_kbn, :cummin, :cummax)
    eval(Expr(:import, :Base, f))
    @eval ($f)(a::NamedArray, d::Integer=1) = fan($f, string($f), a, d)
end


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.

"""
    write_topology(directory::String,file::String,bonds::Array{Array{Int64,N} where N,1},angles::Array{Tuple{Int64,Int64,Int64},1},dihedrals::Array{NTuple{4,Int64},1},impropers::Array{NTuple{4,Int64},1})

Writes out the topology to a user-specified text file
"""
function write_topology(directory::String,file::String,bonds::Array{Array{Int64,N} where N,1},angles::Array{Tuple{Int64,Int64,Int64},1},dihedrals::Array{NTuple{4,Int64},1},impropers::Array{NTuple{4,Int64},1})
    outfile = open(joinpath(directory,file),"w")
    write(outfile,"# Please note that the types of the topology have not been determined and have not been included here.\n\n")
    nbonds = length(bonds)
    nangles = length(angles)
    ndihedrals = length(dihedrals)
    nimpropers = length(impropers)
    write(outfile,"$nbonds bonds\n$nangles angles\n$ndihedrals dihedrals\n$nimpropers impropers\n\n")
    count = 1
    write(outfile,"Bonds \n\n")
    for bond in bonds
        atom1 = bond[1]
        atom2 = bond[2]
        write(outfile,"$count $atom1 $atom2\n")
        count += 1
    end
    count = 1
    write(outfile,"\nAngles\n\n")
    for angle in angles
        atom1 = angle[1]
        atom2 = angle[2]
        atom3 = angle[3]
        write(outfile,"$count $atom1 $atom2 $atom3\n")
        count += 1
    end
    count = 1
    write(outfile,"\nDihedrals\n\n")
    for dihedral in dihedrals
        atom1 = dihedral[1]
        atom2 = dihedral[2]
        atom3 = dihedral[3]
        atom4 = dihedral[4]
        write(outfile,"$count $atom1 $atom2 $atom3 $atom4\n")
        count += 1
    end
    count = 1
    write(outfile,"\nImpropers\n\n")
    for improper in impropers
        atom1 = improper[1]
        atom2 = improper[2]
        atom3 = improper[3]
        atom4 = improper[4]
        write(outfile,"$count $atom1 $atom2 $atom3 $atom4\n")
        count += 1
    end
    close(outfile)
    println("Finished writing topology to $file in $directory")
end

max_stars_repo_name: wapisani/CompMatSciTools

max_issues_repo_name: wapisani/CompMatSciTools

max_forks_repo_name: wapisani/CompMatSciTools


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    VariableJoiner()

Join spatial data with shared domain to produce
spatial data with all variables included.
"""
struct VariableJoiner <: AbstractJoiner
end

function Base.join(sdata₁::AbstractData, sdata₂::AbstractData, joiner::VariableJoiner)
  @assert npoints(sdata₁) == npoints(sdata₂) "cannot join different number of points"

  # retrieve variable names and types
  vars₁ = variables(sdata₁)
  vars₂ = variables(sdata₂)

  # variable names as vectors
  vnames₁ = collect(keys(vars₁))
  vnames₂ = collect(keys(vars₂))

  # find common variable names
  allvars = vnames₁ ∪ vnames₂
  comvars = vnames₁ ∩ vnames₂

  # find common names with same type
  stvars, dtvars = [], []
  for v in comvars
    if nonmissingtype(vars₁[v]) == nonmissingtype(vars₂[v])
      push!(stvars, v)
    else
      push!(dtvars, v)
    end
  end

  # warn if same name has two different types
  isempty(dtvars) || @warn "variables $dtvars have multiple types"

  # create single dictionary with all variables
  pairs = []
  for v in vnames₁
    if v ∈ comvars
      push!(pairs, Symbol(v, 1) => sdata₁[v])
    else
      push!(pairs, v => sdata₁[v])
    end
  end
  for v in vnames₂
    if v ∈ comvars
      push!(pairs, Symbol(v, 2) => sdata₂[v])
    else
      push!(pairs, v => sdata₂[v])
    end
  end

  dict = OrderedDict([(k, v) for (k, v) in pairs])

  georef(dict, domain(sdata₁))
end


text:
VTKMultiblockData(dataset::AbstractVTKSimpleData) = VTKMultiblockData(AbstractStaticVTKData[dataset])

function VTKRectilinearData(dataset::VTKUniformRectilinearData)
    _extents = extents(dataset)
    _dim = dim(dataset)
    point_coords = ntuple(j -> [dataset.origin[j] + (i-1)*dataset.spacing[j] for i in 1:_extents[j]], Val(_dim))
    return VTKRectilinearData(point_coords, dataset.point_data, dataset.cell_data)
end
function Base.convert(::Type{VTKRectilinearData}, dataset::VTKUniformRectilinearData)
    return VTKRectilinearData(dataset)
end

function VTKStructuredData(dataset::VTKUniformRectilinearData)
    return VTKStructuredData(VTKRectilinearData(dataset))
end
function VTKStructuredData(dataset::VTKRectilinearData)
    if dim(dataset) == 2
        xcoord = dataset.point_coords[1]
        ycoord = dataset.point_coords[2]

        point_coords = zeros(dim(dataset), extents(dataset)...)
        for j in 1:extents(dataset, 2), i in 1:extents(dataset, 1)
            point_coords[:, i, j] = [xcoord[i], ycoord[j]]
        end
    elseif dim(dataset) == 3
        xcoord = dataset.point_coords[1]
        ycoord = dataset.point_coords[2]
        zcoord = dataset.point_coords[3]

        point_coords = zeros(dim(dataset), extents(dataset)...)
        for k in 1:extents(dataset, 3), j in 1:extents(dataset, 2), i in 1:extents(dataset, 1)
            point_coords[:, i, j, k] = [xcoord[i], ycoord[j], zcoord[k]]
        end
    else
        throw("Invalid dataset dimensions.")
    end
    return VTKStructuredData(point_coords, dataset.point_data, dataset.cell_data)
end
function Base.convert(::Type{VTKStructuredData}, dataset::VTKUniformRectilinearData)
    return VTKStructuredData(dataset)
end
function Base.convert(::Type{VTKStructuredData}, dataset::VTKRectilinearData)
    return VTKStructuredData(dataset)
end

function VTKUnstructuredData(dataset::VTKPolyData)
    return VTKUnstructuredData(dataset.point_coords, dataset.cell_types,
        dataset.cell_connectivity, dataset.point_data, dataset.cell_data)
end

function VTKUnstructuredData(dataset::T) where {T<:AbstractVTKStructuredData}
    point_coords = zeros(dim(dataset), num_of_points(dataset))
    if T <: VTKStructuredData
        point_coords = unstructured_point_coords_from_structured(dataset)
    elseif T <: VTKRectilinearData
        point_coords = unstructured_point_coords_from_rectilinear(dataset)
    elseif T <: VTKUniformRectilinearData
        point_coords = unstructured_point_coords_from_uniform_rectilinear(dataset)
    end
    point_data = unstructured_point_data_from_structured(dataset)
    cell_types, _cell_connectivity, cell_data = unstructured_cell_info_from_structured(dataset)
    return VTKUnstructuredData(point_coords, cell_types, _cell_connectivity, point_data, cell_data)
end

function VTKUnstructuredData(data_blocks::VTKMultiblockData)
    combined_block = VTKUnstructuredData(deepcopy(data_blocks[1]))
    for block in data_blocks[2:end]
        append!(combined_block, VTKUnstructuredData(block))
    end
    return combined_block
end

function Base.convert(::Type{VTKUnstructuredData}, dataset::VTKPolyData)
    return VTKUntructuredData(dataset)
end
function Base.convert(::Type{VTKUnstructuredData}, dataset::AbstractVTKStructuredData)
    return VTKUntructuredData(dataset)
end
function Base.convert(::Type{VTKUnstructuredData}, dataset::VTKMultiblockData)
    return VTKUntructuredData(dataset)
end

function VTKPolyData(dataset::AbstractStaticVTKData)
    return decompose(VTKUnstructuredData(dataset), "Faces", true)
end
function unstructured_point_data_from_structured(dataset::AbstractVTKStructuredData)
    _num_of_points = num_of_points(dataset)
    _dim = dim(dataset)
    _extents = extents(dataset)

    point_data = empty(dataset.point_data)
    for m in keys(dataset.point_data)
        _var_dim = var_dim(dataset, m, "Point")
        if _var_dim == 1
            point_data[m] = reshape(dataset.point_data[m], (_num_of_points,))
        else
            point_data[m] = reshape(dataset.point_data[m], (_var_dim, _num_of_points))
        end
    end
    point_data
end

function unstructured_cell_info_from_structured(dataset::T) where {T<:AbstractVTKStructuredData}
    _num_of_cells = num_of_cells(dataset)
    _dim = dim(dataset)
    pextents = extents(dataset)
    cextents = cell_extents(dataset)

    local_cell_connectivity(cell_ind) = cell_connectivity(T, pextents, cell_ind)
    if _dim == 2
        if T <: VTKUniformRectilinearData
            cell_types = fill(8, _num_of_cells) #VTK_PIXEL
        else
            cell_types = fill(9, _num_of_cells) #VTK_QUAD
        end
        _cell_connectivity = reshape(map(local_cell_connectivity, ((i,j) for j in 1:cextents[2], i in 1:cextents[1])), (_num_of_cells,))
    elseif _dim == 3
        if T <: VTKUniformRectilinearData
            cell_types = fill(11, _num_of_cells) #VTK_VOXEL
        else
            cell_types = fill(12, _num_of_cells) #VTK_HEXAHEDRON
        end
        _cell_connectivity = reshape(map(local_cell_connectivity, ((i,j,k) for k in 1:cextents[3], j in 1:cextents[2], i in 1:cextents[1])), (_num_of_cells,))
    end

    cell_data = empty(dataset.cell_data)
    for m in keys(dataset.cell_data)
        _var_dim = var_dim(dataset, m, "Cell")
        if _var_dim == 1
            cell_data[m] = reshape(dataset.cell_data[m], (num_of_cells(dataset),))
        else
            cell_data[m] = reshape(dataset.cell_data[m], (_var_dim, num_of_cells(dataset)))
        end
    end
    return cell_types, _cell_connectivity, cell_data
end

function unstructured_point_coords_from_structured(dataset::VTKStructuredData)
    _extents = extents(dataset)
    _num_of_points = num_of_points(dataset)
    _dim = dim(dataset)

    return reshape(dataset.point_coords, (_dim, _num_of_points))
end

function unstructured_point_coords_from_rectilinear(dataset::VTKRectilinearData)
    _num_of_points = num_of_points(dataset)
    _dim = dim(dataset)

    _point_coords = zeros(_dim, _num_of_points)
    p = 1
    if _dim == 2
        for y in dataset.point_coords[2], x in dataset.point_coords[1]
            _point_coords[1,p] = x
            _point_coords[2,p] = y
            p += 1
        end
    elseif _dim == 3
        for z in dataset.point_coords[3], y in dataset.point_coords[2], x in dataset.point_coords[1]
            _point_coords[1,p] = x
            _point_coords[2,p] = y
            _point_coords[3,p] = z
            p += 1
        end
    else
        throw("Invalid dimension.")
    end
    return _point_coords
end

function unstructured_point_coords_from_uniform_rectilinear(dataset::VTKUniformRectilinearData)
    _extents = extents(dataset)
    _num_of_points = num_of_points(dataset)
    _dim = dim(dataset)
    _spacing = dataset.spacing
    _origin = dataset.origin

    _point_coords = zeros(_dim, _num_of_points)
    p = 1
    if _dim == 2
        x, y = _origin
        for j in 1:_extents[2]
            x = _origin[1]
            for i in 1:_extents[1]
                _point_coords[1,p] = x
                _point_coords[2,p] = y
                p += 1
                x += _spacing[1]
            end
            y += _spacing[2]
        end
    elseif _dim == 3
        x, y, z = _origin
        for k in 1:_extents[3]
            y = _origin[2]
            for j in 1:_extents[2]
                x = _origin[1]
                for i in 1:_extents[1]
                    _point_coords[1,p] = x
                    _point_coords[2,p] = y
                    _point_coords[3,p] = z
                    p += 1
                    x += _spacing[1]
                end
                y += _spacing[2]
            end
            z += _spacing[3]
        end
    else
        throw("Invalid dimension.")
    end

    return _point_coords
end


text:
"""
Load the columns of the JDF as vectors
"""

load_columns(jdf::JDFFile; args...) = load_columns(path(jdf); args...)

function load_columns(indir; cols = Symbol[], verbose = false)
    # starting from DataFrames.jl 0.21 the colnames are strings
    cols = string.(cols)

    if VERSION < v"1.3.0"
        return sload_columns(indir, cols = cols, verbose = verbose)
    end

    if verbose
        println("loading $indir in parallel")
    end

    metadatas = open(joinpath(indir, "metadata.jls")) do io
        deserialize(io)
    end

    if length(cols) == 0
        cols = string.(metadatas.names)
    else
        scmn = setdiff(cols, string.(metadatas.names))
        if length(scmn) > 0
            throw("columns $(reduce((x,y) -> string(x) * ", " * string(y), scmn)) are not available, please ensure you have spelt them correctly")
        end
    end

    # get the maximum number of bytes needs to read
    # bytes_needed = maximum(get_bytes.(metadatas.metadatas))

    results = Vector{Any}(undef, length(cols))

    # rate limit channel
    c1 = Channel{Bool}(Threads.nthreads())
    atexit(() -> close(c1))

    i = 1
    for (name, metadata) in zip(metadatas.names, metadatas.metadatas)
        name_str = string(name)
        if name_str in cols
            put!(c1, true)
            results[i] = @spawn begin
                io = BufferedInputStream(open(joinpath(indir, string(name)), "r"))
                new_result = column_loader(metadata.type, io, metadata)
                close(io)
                (name = name_str, task = new_result)
            end
            take!(c1)
            i += 1
        end
    end

    # run the collection of results this serially
    result_vectors = Vector{Any}(undef, length(cols))
    for (i, result) in enumerate(results)
        if verbose
            println("Extracting $(fetch(result).name)")
        end

        new_result = fetch(result).task
        colname = fetch(result).name
        if new_result == nothing
            result_vectors[i] = Vector{Missing}(missing, metadatas.rows)
        else
            result_vectors[i] = new_result
        end
    end

    return result_vectors
end

function sload_columns(indir; cols=Symbol[], verbose=false)
    metadatas = open(joinpath(indir, "metadata.jls")) do io
        deserialize(io)
    end

    if length(cols) == 0
        cols = metadatas.names
    else
        scmn = setdiff(cols, string.(metadatas.names))
        if length(scmn) > 0
            throw("columns $(reduce((x,y) -> string(x) * ", " * string(y), scmn[1:min(8, length(scmn))])) are not available, please ensure you have spelt them correctly")
        end
    end

    # get the maximum number of bytes needs to read
    # bytes_needed = maximum(get_bytes.(metadatas.metadatas))

    # rate limit channel
    #results = Vector{Any}(undef, length(metadatas.names))
    results = Vector{Any}(undef, length(cols))

    i = 1
    for (name, metadata) in zip(metadatas.names, metadatas.metadatas)
        if string(name) in cols
            if verbose
                println("Loading $name")
            end
            results[i] = begin
                io = BufferedInputStream(open(joinpath(indir, string(name)), "r"))
                new_result = column_loader(metadata.type, io, metadata)
                close(io)
                (name = name, task = new_result)
            end
            i += 1
        end
    end

    result_vectors = Vector{Any}(undef, length(results))

    # run this serially
    for (i, result) in enumerate(results)
        if verbose
            println("Extracting $(result.name)")
            println(result.task)
        end

        new_result = result.task
        colname = result.name
        if new_result === nothing
           result_vectors[i] = Vector{Missing}(missing, metadatas.rows)
        else
            result_vectors[i] = new_result
        end
    end

    return result_vectors
end


text:
#
# inventory control problem
#
using DiscreteEvents, Distributions, Random

mutable struct Station
    q::Float64           # fuel amount
    t::Vector{Float64}   # time vector
    qt::Vector{Float64}  # fuel/time vector
    cs::Int              # customers served
    cl::Int              # customers lost
    qs::Float64          # fuel sold
    ql::Float64          # lost sales
end

function customer(c::Clock, s::Station, X::Distribution)
    function fuel(s::Station, x::Float64)
        s.q -= x             # take x from tank
        push!(s.t, c.time)   # record time, amount, customer, sale
        push!(s.qt, s.q)
        s.cs += 1
        s.qs += x
    end

    x = rand(X)              # calculate demand
    if s.q ≥ x               # demand can be met
        fuel(s, x)
    elseif s.q ≥ a           # only partially in stock
        s.ql += x - s.q      # count the loss
        fuel(s, s.q)         # give em all we have
    else
        s.cl += 1            # count the lost customer
        s.ql += x            # count the lost demand
    end
end

function replenish(c::Clock, s::Station, Q::Float64)
    if s.q < a
        push!(s.t, c.time)
        push!(s.qt, s.q)
        s.q += Q
        push!(s.t, c.time)
        push!(s.qt, s.q)
    end
end

Random.seed!(123)
const λ = 0.5      # ~ every two minutes a customer
const ρ = 1/180    # ~ every 3 hours a replenishment truck
const μ = 30       # ~ mean demand per customer
const σ = 10       #   standard deviation
const a = 5        #   minimum amount
const Q = 6000.0   # replenishment amount
const M₁ = Exponential(1/λ)  # customer arrival time distribution
const M₂ = Exponential(1/ρ)  # replenishment time distribution
const X = truncated(Normal(μ, σ), a, Inf)  # demand distribution

clock = Clock()
s = Station(Q, Float64[0.0], Float64[Q], 0, 0, 0.0, 0.0)
@event replenish(clock, s, Q) every M₂
@event customer(clock, s, X) every M₁
println(@run! clock 5000)

@show fuel_sold = s.qs;
@show loss_rate = s.ql/s.qs;
@show served_customers = s.cs;
@show lost_customers = s.cl;

# using Plots
# plot(s.t, s.qt, title="Filling Station", xlabel="time [min]", ylabel="inventory [L]", legend=false)
# savefig("invctrl.png")


text:
import LinearAlgebra
import FileIO
import GeometryBasics
import Observables
import Makie
import GLMakie

using Porta

startframe = 1
FPS = 60
resolution = (3840, 2160)
segments = 30
speed = 1
factor = 0.0025
scale = 1.0
radius = 0.01
number = 360
frames = 3600
basemapradius = 1.025
basemapsegments = 2segments
basepointradius = 0.02
basepointsegments = 15
basepointtransparency = false
exportmode = ["gif", "frames", "video"]
exportmode = exportmode[2]
modelname = "latitudinaltori"

frequencies = begin
    array = []
    N = 5
    n = Int(number ÷ N)
    for i in 1:N
        θ = (i / N) * π - (π / 2)
        θ = 0.99 * θ
        points = [Geographic(1, ϕ, θ) for ϕ in range(0, stop = 2π * 7 / 10, length = n)]
        points = ComplexLine.(points)
        points = [p.z for p in points]
        append!(array, points)
    end
    convert(Array{Complex,1}, array)
end

number = length(frequencies)

# The scene object that contains other visual objects
#Makie.reasonable_resolution() = (800, 800)
scene = Makie.Scene(backgroundcolor = :white,
                               show_axis = false,
                               resolution = resolution)

fmap(b::S²) = b

"""
    getpointonpath(t)

Get a point on a path on S² with the given parameter `t`. t ∈ [0, 1].
"""
getpointonpath(t::Real) = begin
    t = t * 12π
    d = ℯ^cos(t) - 2cos(4t) - sin(t / 12)^5
    x, y = sin(t) * d, cos(t) * d
    ComplexLine(x + im * y)
end

"""
    getbutterflycurve(points)

Get butterfly curve by Temple H. Fay (1989) with the given number of `points`.
"""
function getbutterflycurve(points::Int)
    array = Array{ComplexLine,1}(undef, points)
    # 0 ≤ t ≤ 12π
    for i in 1:points
        t = (i - 1) / points * 12π
        d = ℯ^cos(t) - 2cos(4t) - sin(t / 12)^5
        x, y = sin(t) * d, cos(t) * d
        array[i] = ComplexLine(x + im * y)
    end
    array
end

quaternionicrgb(p::Quaternion, q::Quaternion, s::Real, v::Real) = begin
    hue = acos(dot(p, q)) / π
    hsvtorgb([hue * 360; s; v])
end

"""
    getpoints(z [, segments])

Get a circle of point in the base space with the given center `z` in the Riemann
sphere. The optional argument `segments` determines how many points should the
circle have.
"""
function getpoints(z::S²; segments::Int = 30)
    g = Geographic(z)
    r, ϕ, θ = vec(g)
    lspace = collect(range(float(-pi), stop = float(pi), length = segments))
    array = Array{Geographic,1}(undef, segments)
    for i in 1:segments
        ϕ′, θ′ = g.ϕ + factor * cos(lspace[i]), g.θ + factor * sin(lspace[i])
        array[i] = Geographic(r, ϕ′, θ′)
    end
    array
end

#basemapconfig = Biquaternion(ℝ³(0, -1.25, 0))
basemapconfig = Biquaternion(Quaternion(1, 0, 0, 0), ℝ³(0, -1.1, 0))
basemapcolor = Makie.RGBAf0(0.75, 0.75, 0.75, 0.25)
basemaptransparency = true
basemap = Sphere(basemapconfig,
                 scene,
                 radius = basemapradius,
                 segments = basemapsegments,
                 color = basemapcolor,
                 transparency = basemaptransparency)

curvepoints = getbutterflycurve(number)
solidcolors = Array{Makie.RGBAf0,1}(undef, number)
ghostcolors = similar(solidcolors)
basepoints = Array{Sphere,1}(undef, number)
for i in 1:number
    config = Biquaternion(ℝ³(Cartesian(curvepoints[i])) * basemapradius)
    hue = (i - 1) / number * 360
    rgb = hsvtorgb([hue, 1, 1])
    solidcolors[i] = Makie.RGBAf0(rgb..., 0.9)
    ghostcolors[i] = Makie.RGBAf0(rgb..., 0.1)
    sphere = Sphere(basemapconfig * config,
                    scene,
                    radius = basepointradius,
                    segments = basepointsegments,
                    color = solidcolors[i],
                    transparency = basepointtransparency)
    basepoints[i] = sphere
end

#bundleconfig = Biquaternion(ℝ³(0, 1.25, 0))
bundleconfig = Biquaternion(ℝ³(0, 1.1, 0))
s3rotation = Quaternion(1, 0, 0, 0)
solidfibers = Array{Fiber,1}(undef, number)
ghostfibers = Array{Fiber,1}(undef, number)
α = 40 / 180 * pi
solidtop = U1(pi - 2α)
solidbottom = U1(-pi)
ghosttop = U1(pi)
ghostbottom = solidtop
for i in 1:number
    point = curvepoints[i]
    transparency = false
    solidfiber = Fiber(scene,
                       point,
                       τmap,
                       fmap,
                       radius = radius,
                       top = solidtop,
                       bottom = solidbottom,
                       s3rotation = s3rotation,
                       config = bundleconfig,
                       segments = segments,
                       color = solidcolors[i],
                       transparency = transparency,
                       scale = scale)
    solidfibers[i] = solidfiber
    transparency = true
    ghostfiber = Fiber(scene,
                       point,
                       τmap,
                       fmap,
                       radius = radius,
                       top = ghosttop,
                       bottom = ghostbottom,
                       s3rotation = s3rotation,
                       config = bundleconfig,
                       segments = segments,
                       color = ghostcolors[i],
                       transparency = transparency,
                       scale = scale)
    ghostfibers[i] = ghostfiber
end

basespacepoints = Array{ComplexLine,1}(undef, number)
phases = Array{U1,2}(undef, number, 4)
s3rotations = Array{Quaternion,1}(undef, number)

len = 0.1
width = 0.01
transparency = false
tail, head = ℝ³(rand(3)), ℝ³(rand(3))
arrows = []
for i in 1:number
    hue = (i - 1) / number * 360
    rgb = hsvtorgb([hue, 1, 1])
    color = Makie.RGBAf0(rgb..., 0.9)
    push!(arrows, Arrow(tail, head, scene, width = width, color = color))
end

a, b = ComplexLine(1 + im), ComplexLine(0.5 - 1.5 * im)
c, d = ComplexLine(-0.5 + 3 * im), ComplexLine(-1 + 0 * im)
f(z, a, b, c, d) = (a * z + b) / (c * z + d)

"""
    animate(i)

Update the state of observables with the given frame number `i`.
"""
function animate(i::Int)
    step = (i - 1) / frames
    τ = -step * speed * 4pi
    u = ℝ³(0, 1, 0)
    q = Quaternion(τ, u)
    a′, b′, c′, d′ = [ComplexLine(Cartesian(rotate(ℝ³(Cartesian(point)), q))).z
                      for point in (a, b, c, d)]
    point = getpointonpath(step)
    point = ComplexLine(f(point.z, a′, b′, c′, d′))
    if isnan(vec(point)[1]) || isnan(vec(point)[2])
        point = Geographic(1, 0, π/2)
    end
    point = Cartesian(point)
    s3rotation = Quaternion(1, 0, 0, 0)
    for (index, item) in enumerate(frequencies)
        z = item
        z = ComplexLine(f(2π * z, a′, b′, c′, d′))
        #println("z = ", z)
        if isnan(vec(z)[1]) || isnan(vec(z)[2])
            z = ComplexLine(Geographic(1, 0, -π / 2))
        end
        r, ϕ, θ = vec(Geographic(z))
        if isnan(r)
            z = ComplexLine(Geographic(1, ϕ, θ * 0.99))
        end

        basespacepoints[index] = z
        #basespacecircles[index] = getpoints(z′, segments = basepointsegments)
        realpart = real(z) / abs(z)
        imaginarypart = imag(z) / abs(z)
        #magnitude = abs(item)
        #α = magnitude * 2π
        #s = U1(α)
        α = angle(z)
        g = Quaternion(0, ℝ³(Cartesian(z)))
        q = Quaternion(α, ℝ³(0, 0, 1))
        #h = S¹action(τmap(point), s) * g
        rgb = quaternionicrgb(g, q, 0.1 * abs(realpart) + 0.1 * rand() + 0.8,
                              0.1 * abs(imaginarypart) + 0.1 * rand() + 0.8)
        solidcolors[index] = Makie.RGBAf0(rgb..., abs(realpart))
        ghostcolors[index] = Makie.RGBAf0(rgb..., abs(imaginarypart))
        solidtop = U1(0)
        solidbottom = U1(angle(z))
        ghosttop = solidtop
        ghostbottom = U1(2π)
        phases[index, 1] = solidtop
        phases[index, 2] = solidbottom
        phases[index, 3] = ghosttop
        phases[index, 4] = ghostbottom
        #s3rotations[index] = g
    end
    for (index, item) in enumerate(basepoints)
        z = basespacepoints[index]
        p = ℝ³(Cartesian(z)) * basemapradius
        config = Biquaternion(p)
        update(item, basemapconfig * config)
        update(item, solidcolors[index])
    end
    for (index, item) in enumerate(solidfibers)
        solidtop = phases[index, 1]
        solidbottom = phases[index, 2]
        update(item, basespacepoints[index], s3rotation, radius, solidtop,
               solidbottom)
        update(item, solidcolors[index])
    end
    for (index, item) in enumerate(ghostfibers)
        ghosttop = phases[index, 3]
        ghostbottom = phases[index, 4]
        update(item, basespacepoints[index], s3rotation, radius, ghosttop,
               ghostbottom)
        update(item, ghostcolors[index])
    end

    # update tangent vectors
    for (index, item) in enumerate(arrows)
        w = basespacepoints[index]
        tail = ℝ³(Cartesian(w)) * basemapradius
        initial = ℝ³(0, 0, 1) * basemapradius
        q = getrotation(normalize(tail), normalize(initial))
        α = phases[index, 2]
        head = initial + ℝ³(1, 0, 0) * (0.01 + angle(α) / 5)
        head = rotate(head, q)
        q = Quaternion(angle(α) / 2, normalize(tail))
        head = rotate(head, q)
        tail, head = applyconfig([tail; head], basemapconfig)
        update(item, tail, head - tail)
        update(item, solidcolors[index])
    end

#=     z = Complex(0 + 0im)
    z = ComplexLine(f(π * z, a′, b′, c′, d′))
    z = ℝ³(Cartesian(z))

    z₀ = Complex(0 + im)
    z₀ = ComplexLine(f(π * z₀, a′, b′, c′, d′))
    z₀ = ℝ³(Cartesian(z₀)) =#

    #n = normalize(z₀ - z)
    #v = ℝ³(-1, 0, 1) * 2.5
    #v = ℝ³(-1, 0, 1) * 2.1
    #distance = norm(v)
    #v = normalize(z) * distance
    # update eye position
    # scene.camera.eyeposition.val
    #upvector = GeometryBasics.Vec3f0(vec(n)...)
    #eyeposition = GeometryBasics.Vec3f0(vec(v)...)
    #lookat = GeometryBasics.Vec3f0(0, 0, 0)
    #Makie.update_cam!(scene, eyeposition, lookat, upvector)
    #scene.center = false # prevent scene from recentering on display
end

n = ℝ³(0, 0, 1)
v = ℝ³(-1, 0, 1) * 2.5
# update eye position
# scene.camera.eyeposition.val
upvector = GeometryBasics.Vec3f0(vec(n)...)
eyeposition = GeometryBasics.Vec3f0(vec(v)...)
lookat = GeometryBasics.Vec3f0(0, 0, 0)
Makie.update_cam!(scene, eyeposition, lookat, upvector)
scene.center = false # prevent scene from recentering on display
if exportmode ∈ ["gif", "video"]
    outputextension = exportmode == "gif" ? "gif" : "mkv"
    Makie.record(scene, "gallery/$modelname.$outputextension",
                            framerate = FPS) do io
        for i in startframe:frames
            animate(i) # animate the scene
            Makie.recordframe!(io) # record a new frame
            step = (i - 1) / frames
            println("Completed step $(100step).\n")
        end
    end
elseif exportmode == "frames"
    directory = joinpath("gallery", modelname)
    !isdir(directory) && mkdir(directory)
    directory = joinpath("gallery", modelname)
    !isdir(directory) && mkdir(directory)
    for i in startframe:frames
        start = time()
        animate(i) # animate the scene
        elapsed = time() - start
        println("Generating frame $i took $elapsed (s).")

        start = time()

        paddingnumber = length(digits(frames))
        imageid = lpad(i, paddingnumber, "0")
        imagename = "$(modelname)_$(resolution[2])p_$(FPS)fps_$imageid.jpeg"
        filepath = joinpath(directory, imagename)
        FileIO.save(filepath,
                    scene;
                    resolution = resolution,
                    pt_per_unit = 400.0,
                    px_per_unit = 400.0)
        elapsed = time() - start
        println("Saving file $filepath took $elapsed (s).")

        step = (i - 1) / frames
        println("Completed step $(100step).\n")
        sleep(1)

        stitch() = begin
            part = i ÷ (frames / 10)
            exportdir = joinpath(directory, "export")
            !isdir(exportdir) && mkdir(exportdir)
            WxH = "$(resolution[1])x$(resolution[2])"
            #WxH = "1920x1080"
            commonpart = "$(modelname)_$(resolution[2])p_$(FPS)fps"
            inputname = "$(commonpart)_%0$(paddingnumber)d.jpeg"
            inputpath = joinpath(directory, inputname)
            outputname = "$(commonpart)_$(part).mp4"
            outputpath = joinpath(exportdir, outputname)
            command1 = `ffmpeg -y -f image2 -framerate $FPS -i $inputpath -s $WxH -pix_fmt yuvj420p $outputpath`
            run(command1)
            sleep(1)
        end

        if isapprox(i % (frames / 10), 0)
            stitch()
            sleep(1)
        end
    end
end

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
@testset "log" begin
    all_types = (RotMatrix3, RotMatrix{3}, AngleAxis, RotationVec,
                 QuatRotation, RodriguesParam, MRP,
                 RotXYZ, RotYZX, RotZXY, RotXZY, RotYXZ, RotZYX,
                 RotXYX, RotYZY, RotZXZ, RotXZX, RotYXY, RotZYZ,
                 RotX, RotY, RotZ,
                 RotXY, RotYZ, RotZX, RotXZ, RotYX, RotZY,
                 RotMatrix2, RotMatrix{2}, Angle2d)

    @testset "$(T)" for T in all_types, F in (one, rand)
        R = F(T)
        @test R ≈ exp(log(R))
        @test log(R) isa RotationGenerator
        @test exp(log(R)) isa Rotation
    end

    @testset "$(N)-dim" for N in 1:5
        M = @SMatrix rand(N,N)
        R = nearest_rotation(M)
        @test isrotationgenerator(log(R))
        @test log(R) isa RotMatrixGenerator
        @test exp(log(R)) isa RotMatrix
    end
end

@testset "exp(zero)" begin
    all_types = (RotMatrixGenerator{3}, RotationVecGenerator,
                 RotMatrixGenerator{2}, Angle2dGenerator)

    @testset "$(T)" for T in all_types
        r = zero(T)
        @test one(exp(r)) ≈ exp(r)
        @test exp(r) isa Rotation
    end
end

@testset "exp(::RotMatrixGenerator)" begin
    for N in 2:3
        r = zero(RotMatrixGenerator{N})
        @test r isa RotMatrixGenerator{N}
        @test exp(r) isa RotMatrix{N}
    end
end


text:
__precompile__()

module MiniPETSc

using Reexport

@reexport using MPI

const library = string(ENV["PETSC_DIR"], "/", ENV["PETSC_ARCH"], "/lib/libpetsc")

include("PetscTypes.jl")

function finalize()
    ccall((:PetscFinalize, library), PetscErrorCode, ())
    # Somehow doesn't work here
    # MPI.Finalize()
end

function __init__()
    args = vcat("julia", ARGS)
    nargs = length(args)
    MPI.Init()
    ccall((:PetscInitializeNoPointers, library), PetscErrorCode, (Cint, Ptr{Ptr{UInt8}}, Cstring, Cstring), nargs, args, C_NULL, C_NULL)

    # Cleanup at the end
    atexit(finalize)
end

export PetscMat
export setSize!
export setPreallocation!
export assemble!
export viewMat
export zero!
export zeroRows!

export PetscVec
export GhostedPetscVec
export setSize!
export assemble!
export viewVec
export plusEquals!
export zero!
export copy!
export serializeToZero

export PetscKSP
export setOperators
export solve!

include("Mat.jl")
include("Vec.jl")
include("KSP.jl")

end


text:
"""
    Auroral()

Aurora with linear regression.
"""
struct Auroral <: AbstractAurora end

Base.@kwdef struct FittedAuroral{A, Ts, Ms, Ss} <: AbstractFittedAurora
    βs::Ts
    βs_list::Vector{Ts}
    μs_mat::Ms
    μs::Ss
    method::A
end

function StatsBase.fit(auroral::Auroral, Zs::AbstractVector{<:ReplicatedSample})
    K = nobs(Zs[1])
    n = length(Zs)
    any(nobs.(Zs) != K) || throw("All Zs should have the same number of replicates.")
    μs_mat = zeros(n,K)

    # See https://github.com/dmbates/CopenhagenEcon/blob/master/jmd/03-LinearAlgebra.jmd
    cache_vec = zeros(K+1)
    cache_mat = zeros(K+1,n)
    cache_cov_mat = zeros(K+1, K+1)

    βs_list = Vector{typeof(cache_vec)}(undef, K)

    for j in Base.OneTo(K)
        fill!(cache_mat, 0)
        fill!(cache_cov_mat, 0)
        fill!(cache_vec, 0)

        for (i,Z) in enumerate(Zs)
            design_response_row!(view(cache_mat, :, i), Z, j)
        end

        mul!(cache_cov_mat, cache_mat, cache_mat')
        chr = cholesky!(cache_cov_mat)

        RXX = UpperTriangular(view(chr.U, 1:K, 1:K))
        βs_list[j] = ldiv!(RXX, copy(chr.U[1:K, end]))

        for (i,Z) in enumerate(Zs)
            design_response_row!(cache_vec, Z, j)
            μs_mat[i,j] = dot(βs_list[j], view(cache_vec,1:K))
        end
    end
    μs = vec(mean(μs_mat; dims=2))
    βs = vec(mean(hcat(βs_list...), dims=2))
    FittedAuroral(βs = βs, βs_list = βs_list,
                μs = μs, μs_mat = μs_mat,
                method = auroral)
end

struct CoeyCunningham <: AbstractAurora end

function StatsBase.fit(auroral::CoeyCunningham, Zs::AbstractVector{<:ReplicatedSample})
    K = nobs(Zs[1])
    n = length(Zs)
    any(nobs.(Zs) != K) || throw("All Zs should have the same number of replicates.")
    μs_mat = zeros(n,K)

    cache_means = zeros(n)
    cache_responses = zeros(n)
    cache_vec = zeros(K-1)

    βs_list = Vector{typeof(cache_vec)}(undef, K)

    for j in Base.OneTo(K)
        cache_responses .= StatsBase.response.(Zs, j)
        fill!(cache_means, 0)

        for (i,Z) in enumerate(Zs)
            cache_means[i] = mean(design_row!(cache_vec, Z, j))
        end

        β_slope = cov(cache_means, cache_responses)/var(cache_means)
        β_intercept = mean(cache_responses) - β_slope*mean(cache_means)

        βs_list[j] = [β_intercept; β_slope]

        μs_mat[:,j] = β_intercept .+ β_slope .* cache_means
    end
    μs = vec(mean(μs_mat; dims=2))
    βs = vec(mean(hcat(βs_list...), dims=2))
    FittedAuroral(βs = βs, βs_list = βs_list,
                μs = μs, μs_mat = μs_mat,
                method = auroral)
end

# Plotting

@recipe function f(aur::Aurora.FittedAuroral{Auroral}; stick_color=:purple, x_nudge = 0.0)

    βs = aur.βs
    K = length(βs) - 1

    xs = (0:K) .+ x_nudge

    xlabs = ["Int."; [L"X^{(%$k)}" for k in Base.OneTo(K)]]

	background_color_legend --> :transparent
	foreground_color_legend --> :transparent
    grid --> nothing
    label --> ""
    linewidth --> 2.5
    seriestype := :sticks
    seriescolor := stick_color

    @series begin
        seriestype := :scatter
        seriescolor := stick_color
        markerstrokewidth := 0
        markersize :=5
        label := ""
        xs, βs
    end

    xticks := (xs, xlabs)
    yguide --> L"\hat{\beta}"
    xs, βs
end


text:
# D(u)_function.jl

"""
Implementation of FHIP's susceptibility formula D(u) for the Optical Absorption of Polarons.

See FHIP 1962, equation (35c):
https://link.aps.org/doi/10.1103/PhysRev.127.1004.
"""

"""
ℜD(x::Float64, y::Float64, v::Float64, w::Float64, β::Float64)

    Calculate the real part of D(u) (equation (35c) in FHIP) for a complex argument u = x + iy. v and w are the variational Polaron parameters that minimise the free energy, for the supplied α Frohlich coupling. β is the reduced thermodynamical beta = ħω/kT with k being Boltzman's constant.
"""
function ℜD(x, y, v, w, β)
    x = BigFloat(x)
    y = BigFloat(y)
    β = BigFloat(β)
    v = BigFloat(v)
    w = BigFloat(w)

    R = (v^2 - w^2) / (w^2 * v)
    a_squared = β^2 / 4 + R * β * coth(β * v / 2)
    b = R * β / sinh(β * v / 2)

    w^2 * (a_squared - β^2 / 4 - b * cos(v * x) * cosh(v * (y - β / 2)) + x^2 + y * (β - y)) / (β * v^2)
end

"""
ℑD(x::Float64, y::Float64, v::Float64, w::Float64, β::Float64)

    Calculate the imaginary part of D(u) (equation (35c) in FHIP) for a complex argument u = x + iy. v and w are the variational Polaron parameters that minimise the free energy, for the supplied α Frohlich coupling. β is the reduced thermodynamical beta = ħω/kT with k being Boltzman's constant.
"""
function ℑD(x, y, v, w, β)
    x = BigFloat(x)
    y = BigFloat(y)
    β = BigFloat(β)
    v = BigFloat(v)
    w = BigFloat(w)

    R = (v^2 - w^2) / (w^2 * v)
    a_squared = β^2 / 4 + R * β * coth(β * v / 2)
    b = R * β / sinh(β * v / 2)

    w^2 * (b * sin(v * x) * sinh(v * (y - β / 2)) + 2 * x * (y - β / 2)) / (β * v^2)
end

# Create a complex number from the real and imaginary parts of D(u).

D(x, y, v, w, β) = ℜD(x, y, v, w, β) + 1im * ℑD(x, y, v, w, β)


text:
module robot_arm
using ReachabilityAnalysis
@taylorize function robot_arm!(dx, x, params, t)
    dx[1] = x[3]
    dx[2] = x[4]
    dx[3] = (-2x[2] * x[3] * x[4] - 2x[1] - 2x[3] + 4) / (x[2]^2 + 1)
    dx[4] = x[2] * x[3]^2 - x[2] - x[4] + 1
    return dx
end

function model(X0)
    S = @system(x' = robot_arm!(x), dim:4)
    return IVP(S, X0)
end
end # module

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
###############################################################
# Part of BDisposal.jl package
###############################################################

using LinearAlgebra

# Dispatching to either convexProblem or nonConvexProblem
function problem(gI₀,bI₀,gO₀,bO₀,gI,bI,gO,bO;
                 retToScale,prodStructure,convexAssumption=true,
                 directions=(),
                 startValues=(),crossTime=false,forceLinearModel=false)

     if convexAssumption
         return convexProblem(gI₀,bI₀,gO₀,bO₀,gI,bI,gO,bO;
                          retToScale=retToScale,prodStructure=prodStructure,
                          directions=directions,
                          startValues=startValues,crossTime=crossTime,forceLinearModel=forceLinearModel)
     else
         return nonConvexProblem(gI₀,bI₀,gO₀,bO₀,gI,bI,gO,bO;
                          retToScale=retToScale,prodStructure=prodStructure,
                          directions=directions)
     end

end

function convexProblem(inp₀,bInp₀,gO₀,bO₀,inp,bInp,gO,bO;
                       retToScale,prodStructure,
                       directions,
                       startValues,crossTime=false,forceLinearModel=false)

   nI, nGO, nBO, nDMUs, nbI = length(inp₀), length(gO₀), length(bO₀), size(inp,1), length(bInp₀)
   (dirGI,dirBI,dirGO,dirBO) = directions
   if !forceLinearModel
       (startθ,startμ,startλ)    = startValues
   end

   if prodStructure == "multiplicative"
       if !forceLinearModel
               # Redirecting the output to /deV/null as no way to suppress ipopt output
               oldstd = stdout
               redirect_stdout(open("nul", "w"))

               # Model declaration (efficiency model)
               effmodel = Model(() -> AmplNLWriter.Optimizer(Ipopt.amplexe))
               #effmodel = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
               MOI.set(effmodel, MOI.RawParameter("print_level"), 0)

               # Defining variables
               @variables effmodel begin
                   0 <= θ[z in 1:nDMUs] <= Inf, (start = startθ) # thetas (nDMUs)
                   0 <= μ[z in 1:nDMUs] <= Inf, (start = startμ) # mus (nDMUs)
               end
               if !crossTime
                  @variable(effmodel, 1 <= λ <= Inf, start = startλ)
               else
                  @variable(effmodel, 0 <= λ <= Inf, start = startλ)
               end

               # Defining constraints
               @NLconstraints effmodel begin
                   cgInp1[i in 1:nI],   # constraint on good inputs
                      λ^dirGI * inp₀[i] >= sum(θ[z]*inp[z,i] for z in 1:nDMUs)
                   cbInp1[i in 1:nbI], # constraint on bad inputs
                      λ^dirBI * bInp₀[i] <= sum(θ[z]*bInp[z,i] for z in 1:nDMUs)
                   cgInp2[i in 1:nI],   # second constraint on good inputs
                      λ^dirGI * inp₀[i] >= sum(μ[z]*inp[z,i] for z in 1:nDMUs)
                   cbInp2[i in 1:nbI], # second constraint on bad inputs
                      λ^dirBI * bInp₀[i] >= sum(μ[z]*bInp[z,i] for z in 1:nDMUs)
                   cgO1[j in 1:nGO],   # first constraint on good outputs
                      λ^dirGO * gO₀[j] <= sum(θ[z]*gO[z,j] for z in 1:nDMUs)
                   cbO1[j in 1:nBO],   # first constraint on bad outputs
                      λ^dirBO * bO₀[j] >= sum(θ[z]*bO[z,j] for z in 1:nDMUs)
                   cgO2[j in 1:nGO],   # second constraint on good outputs
                      λ^dirGO * gO₀[j] <= sum(μ[z]*gO[z,j] for z in 1:nDMUs)
                   cbO2[j in 1:nBO],   # second constraint on bad outputs
                      λ^dirBO * bO₀[j] <= sum(μ[z]*bO[z,j] for z in 1:nDMUs)
               end
               if retToScale == "variable"
                   @constraints effmodel begin
                       sumtheta,
                         sum(θ[z] for z in 1:1:nDMUs) == 1
                       summu,
                         sum(μ[z] for z in 1:1:nDMUs) == 1
                   end
               end
               # Defining objective
               @objective effmodel Max begin
                   λ
               end
               # Printing the model (for debugging)
               #print(effmodel) # The model in mathematical terms is printed

               # Solving the model
               JuMP.optimize!(effmodel)
               redirect_stdout(oldstd) # recover original stdout

               # Copy the results to the output matrix...
               status = termination_status(effmodel)

               #if (status == MOI.OPTIMAL || status == MOI.LOCALLY_SOLVED || status == MOI.TIME_LIMIT) && has_values(effmodel)
               if (status == MOI.OPTIMAL || status == MOI.LOCALLY_SOLVED) && has_values(effmodel)
                   return value(λ)
               else
                   return missing
               end
       else #forceLinearModel is then true

               # Model declaration (efficiency model)
               #effmodel = Model(() -> AmplNLWriter.Optimizer(path, ["print_level=0"]))
               #effmodel = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
               effmodel = Model(optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => GLPK.GLP_MSG_OFF)) # we choose GLPK with a verbose output
                #effmodel = Model(optimizer_with_attributes(GLPK.Optimizer))
               #effmodel = Model(optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => GLPK.GLP_MSG_OFF))

               # Defining variables
               @variables effmodel begin
                   0 <= θ[z in 1:nDMUs] <= Inf# thetas (nDMUs)
                   0 <= μ[z in 1:nDMUs] <= Inf # mus (nDMUs)
               end

               if dirGO == 1
                   if !crossTime # we  keep the original λ
                  @variable(effmodel, 1 <= λ <= Inf)
                   else
                      @variable(effmodel, 0 <= λ <= Inf)
                   end
                   # @variable(effmodel, λ)

               else
                   if !crossTime # we use as variable its reverse 1/λ
                      @variable(effmodel, 0 <= invλ <= 1)
                   else
                      @variable(effmodel, 0 <= invλ <= Inf)
                   end
                   # @variable(effmodel, invλ)
               end

               if directions == (-1,0,0,0)
                   # Defining constraints
                   @constraints effmodel begin
                       cgInp1[i in 1:nI],   # constraint on good inputs
                          invλ * inp₀[i] >= sum(θ[z]*inp[z,i] for z in 1:nDMUs)
                       cbInp1[i in 1:nbI], # constraint on bad inputs
                           bInp₀[i] <= sum(θ[z]*bInp[z,i] for z in 1:nDMUs)
                       cgInp2[i in 1:nI],   # second constraint on good inputs
                          invλ * inp₀[i] >= sum(μ[z]*inp[z,i] for z in 1:nDMUs)
                       cbInp2[i in 1:nbI], # second constraint on bad inputs
                           bInp₀[i] >= sum(μ[z]*bInp[z,i] for z in 1:nDMUs)
                       cgO1[j in 1:nGO],   # first constraint on good outputs
                           gO₀[j] <= sum(θ[z]*gO[z,j] for z in 1:nDMUs)
                       cbO1[j in 1:nBO],   # first constraint on bad outputs
                           bO₀[j] >= sum(θ[z]*bO[z,j] for z in 1:nDMUs)
                       cgO2[j in 1:nGO],   # second constraint on good outputs
                           gO₀[j] <= sum(μ[z]*gO[z,j] for z in 1:nDMUs)
                       cbO2[j in 1:nBO],   # second constraint on bad outputs
                           bO₀[j] <= sum(μ[z]*bO[z,j] for z in 1:nDMUs)
                   end
              elseif directions == (0,-1,0,0)
                  # Defining constraints
                  @constraints effmodel begin
                      cgInp1[i in 1:nI],   # constraint on good inputs
                          inp₀[i] >= sum(θ[z]*inp[z,i] for z in 1:nDMUs)
                      cbInp1[i in 1:nbI], # constraint on bad inputs
                         invλ * bInp₀[i] <= sum(θ[z]*bInp[z,i] for z in 1:nDMUs)
                      cgInp2[i in 1:nI],   # second constraint on good inputs
                          inp₀[i] >= sum(μ[z]*inp[z,i] for z in 1:nDMUs)
                      cbInp2[i in 1:nbI], # second constraint on bad inputs
                         invλ * bInp₀[i] >= sum(μ[z]*bInp[z,i] for z in 1:nDMUs)
                      cgO1[j in 1:nGO],   # first constraint on good outputs
                          gO₀[j] <= sum(θ[z]*gO[z,j] for z in 1:nDMUs)
                      cbO1[j in 1:nBO],   # first constraint on bad outputs
                          bO₀[j] >= sum(θ[z]*bO[z,j] for z in 1:nDMUs)
                      cgO2[j in 1:nGO],   # second constraint on good outputs
                          gO₀[j] <= sum(μ[z]*gO[z,j] for z in 1:nDMUs)
                      cbO2[j in 1:nBO],   # second constraint on bad outputs
                          bO₀[j] <= sum(μ[z]*bO[z,j] for z in 1:nDMUs)
                  end
              elseif directions == (0,0,1,0)
                  @constraints effmodel begin
                      cgInp1[i in 1:nI],   # constraint on good inputs
                          inp₀[i] >= sum(θ[z]*inp[z,i] for z in 1:nDMUs)
                      cbInp1[i in 1:nbI], # constraint on bad inputs
                          bInp₀[i] <= sum(θ[z]*bInp[z,i] for z in 1:nDMUs)
                      cgInp2[i in 1:nI],   # second constraint on good inputs
                          inp₀[i] >= sum(μ[z]*inp[z,i] for z in 1:nDMUs)
                      cbInp2[i in 1:nbI], # second constraint on bad inputs
                          bInp₀[i] >= sum(μ[z]*bInp[z,i] for z in 1:nDMUs)
                      cgO1[j in 1:nGO],   # first constraint on good outputs
                         λ * gO₀[j] <= sum(θ[z]*gO[z,j] for z in 1:nDMUs)
                      cbO1[j in 1:nBO],   # first constraint on bad outputs
                          bO₀[j] >= sum(θ[z]*bO[z,j] for z in 1:nDMUs)
                      cgO2[j in 1:nGO],   # second constraint on good outputs
                         λ * gO₀[j] <= sum(μ[z]*gO[z,j] for z in 1:nDMUs)
                      cbO2[j in 1:nBO],   # second constraint on bad outputs
                         bO₀[j] <= sum(μ[z]*bO[z,j] for z in 1:nDMUs)
                  end
              elseif directions == (0,0,0,-1)
                  # Defining constraints
                  @constraints effmodel begin
                      cgInp1[i in 1:nI],   # constraint on good inputs
                         inp₀[i] >= sum(θ[z]*inp[z,i] for z in 1:nDMUs)
                      cbInp1[i in 1:nbI], # constraint on bad inputs
                         bInp₀[i] <= sum(θ[z]*bInp[z,i] for z in 1:nDMUs)
                      cgInp2[i in 1:nI],   # second constraint on good inputs
                         inp₀[i] >= sum(μ[z]*inp[z,i] for z in 1:nDMUs)
                      cbInp2[i in 1:nbI], # second constraint on bad inputs
                          bInp₀[i] >= sum(μ[z]*bInp[z,i] for z in 1:nDMUs)
                      cgO1[j in 1:nGO],   # first constraint on good outputs
                          gO₀[j] <= sum(θ[z]*gO[z,j] for z in 1:nDMUs)
                      cbO1[j in 1:nBO],   # first constraint on bad outputs
                         invλ * bO₀[j] >= sum(θ[z]*bO[z,j] for z in 1:nDMUs)
                      cgO2[j in 1:nGO],   # second constraint on good outputs
                         gO₀[j] <= sum(μ[z]*gO[z,j] for z in 1:nDMUs)
                      cbO2[j in 1:nBO],   # second constraint on bad outputs
                         invλ * bO₀[j] <= sum(μ[z]*bO[z,j] for z in 1:nDMUs)
                  end

              else
                  @error "Directions not suppported when the option forceLinearModel is true"
              end
               if retToScale == "variable"
                   @constraints effmodel begin
                       sumtheta,
                         sum(θ[z] for z in 1:1:nDMUs) == 1
                       summu,
                         sum(μ[z] for z in 1:1:nDMUs) == 1
                   end
               end

               # Defining objective
               if dirGO == 1
                   @objective effmodel Max begin
                       λ
                   end
               else
                   @objective effmodel Min begin
                       invλ
                   end
               end
               # Printing the model (for debugging)
               #print(effmodel) # The model in mathematical terms is printed

               JuMP.optimize!(effmodel)

               # Copy the results to the output matrix...
               status = termination_status(effmodel)

               #if (status == MOI.OPTIMAL || status == MOI.LOCALLY_SOLVED || status == MOI.TIME_LIMIT) && has_values(effmodel)
               if (status == MOI.OPTIMAL || status == MOI.LOCALLY_SOLVED) && has_values(effmodel)
                   if dirGO == 1
                       return value(λ)
                   else
                       return 1/value(invλ)
                   end
               else
                   return missing
               end

        end   # end of forCeLinearModel == true
   else # prodStructure is then additive...
       # Model declaration (efficiency model)
       effmodel = Model(optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => GLPK.GLP_MSG_OFF)) # we choose GLPK with a verbose output
       # Defining variables
       @variables effmodel begin
           θ[z in 1:nDMUs] >= 0 # thetas (nDMUs)
           μ[z in 1:nDMUs] >= 0 # mus (nDMUs)
       end
       if !crossTime
          @variable(effmodel, 0 <= λ <= Inf)
       else
          @variable(effmodel, -Inf <= λ <= Inf)
       end
       # Defining constraints
       # additive structure
       @constraints effmodel begin
           cgInp1[i in 1:nI],   # constraint on good inputs
              inp₀[i]- λ*dirGI * inp₀[i] >= sum(θ[z]*inp[z,i] for z in 1:nDMUs)
           cbInp1[i in 1:nbI], # constraint on bad inputs
              bInp₀[i]- λ*dirBI * bInp₀[i] <= sum(θ[z]*bInp[z,i] for z in 1:nDMUs)
           cgInp2[i in 1:nI],   # second constraint on good inputs
               inp₀[i] - λ*dirGI * inp₀[i] >= sum(μ[z]*inp[z,i] for z in 1:nDMUs)
           cbInp2[i in 1:nbI], # second constraint on bad inputs
              bInp₀[i] - λ*dirBI * bInp₀[i] >= sum(μ[z]*bInp[z,i] for z in 1:nDMUs)
           cgO1[j in 1:nGO],   # first constraint on good outputs
              gO₀[j] + λ * dirGO * gO₀[j] <= sum(θ[z]*gO[z,j] for z in 1:nDMUs)
           cbO1[j in 1:nBO],   # first constraint on bad outputs
              bO₀[j] - λ * dirBO * bO₀[j] >= sum(θ[z]*bO[z,j] for z in 1:nDMUs)
           cgO2[j in 1:nGO],   # second constraint on good outputs
              gO₀[j] + λ*dirGO * gO₀[j] <= sum(μ[z]*gO[z,j] for z in 1:nDMUs)
           cbO2[j in 1:nBO],   # second constraint on bad outputs
              bO₀[j] - λ*dirBO * bO₀[j] <= sum(μ[z]*bO[z,j] for z in 1:nDMUs)
       end
       if retToScale == "variable"
           @constraints effmodel begin
               sumtheta,
                 sum(θ[z] for z in 1:1:nDMUs) == 1
               summu,
                 sum(μ[z] for z in 1:1:nDMUs) == 1
           end
       end
       # Defining objective
       @objective effmodel Max begin
           λ
       end
       # Printing the model (for debugging)
       # print(effmodel) # The model in mathematical terms is printed

       # Solving the model
       optimize!(effmodel)

       # Copy the results to the output matrix...
       status = termination_status(effmodel)

       #if (status == MOI.OPTIMAL || status == MOI.LOCALLY_SOLVED || status == MOI.TIME_LIMIT) && has_values(effmodel)
       if (status == MOI.OPTIMAL || status == MOI.LOCALLY_SOLVED) && has_values(effmodel)
           return value(λ)
       else
           return missing
       end
   end # ending addictive  prodStructure
end # ending convexProblem

function nonConvexProblem(gI₀,bI₀,gO₀,bO₀,gI,bI,gO,bO;
                       retToScale,prodStructure,
                       directions)

    if retToScale == "constant" && directions != (0,0,1,-1) @error "Not supported" end
    nGI, nGO, nBO, nDMUs, nBI = length(gI₀), length(gO₀), length(bO₀), size(gI,1), length(bI₀)
    (dirGI,dirBI,dirGO,dirBO) = directions

    if directions == (-1,0,0,0) && prodStructure == "multiplicative"
        gI_ratio  =  gI ./ gI₀' # nDMU x ngI
        # Normal dist function
        bIConstraint = all(bI₀'  .>=   bI, dims=2)
        gOConstraint = all(gO₀'  .<=   gO, dims=2)
        bOConstraint = all(bO₀'  .<=   bO, dims=2)
        globalContraint = dropdims(all(hcat(bIConstraint,gOConstraint,bOConstraint), dims=2), dims=2)
        effScore_normal = ! any(globalContraint) ? missing : minimum(maximum(gI_ratio[globalContraint,:],dims=2))
        # Bdisposal Distance function
        bIConstraint = all(bI₀'  .>=   bI, dims=2)
        gOConstraint = all(gO₀'  .<=   gO, dims=2)
        bOConstraint = all(bO₀'  .>=   bO, dims=2)
        globalContraint = dropdims(all(hcat(bIConstraint,gOConstraint,bOConstraint), dims=2), dims=2)
        effScore_bfrontier = ! any(globalContraint) ? missing : minimum(maximum(gI_ratio[globalContraint,:],dims=2))
        effscore = max(effScore_normal,effScore_bfrontier)
        return effscore
    elseif  directions == (1,0,0,0) && prodStructure == "additive"
        gI_ratio  =  1 .- (gI ./ gI₀') # nDMU x ngI
        # Normal dist function
        bIConstraint = all(bI₀'  .>=   bI, dims=2)
        gOConstraint = all(gO₀'  .<=   gO, dims=2)
        bOConstraint = all(bO₀'  .<=   bO, dims=2)
        globalContraint = dropdims(all(hcat(bIConstraint,gOConstraint,bOConstraint), dims=2), dims=2)
        effScore_normal = ! any(globalContraint) ? missing : maximum(minimum(gI_ratio[globalContraint,:],dims=2))
        # Bdisposal Distance function
        bIConstraint = all(bI₀'  .>=   bI, dims=2)
        gOConstraint = all(gO₀'  .<=   gO, dims=2)
        bOConstraint = all(bO₀'  .>=   bO, dims=2)
        globalContraint = dropdims(all(hcat(bIConstraint,gOConstraint,bOConstraint), dims=2), dims=2)
        effScore_bfrontier = ! any(globalContraint) ? missing : maximum(minimum(gI_ratio[globalContraint,:],dims=2))
        effscore = min(effScore_normal,effScore_bfrontier)
        return effscore

    elseif  directions == (0,-1,0,0) && prodStructure == "multiplicative"
        bI_ratio  =  bI ./ bI₀' # nDMU x ngI
        # Normal dist function
        gIConstraint = all(gI₀'  .>=   gI, dims=2)
        gOConstraint = all(gO₀'  .<=   gO, dims=2)
        bOConstraint = all(bO₀'  .<=   bO, dims=2)
        globalContraint = dropdims(all(hcat(gIConstraint,gOConstraint,bOConstraint), dims=2), dims=2)
        effScore_normal = ! any(globalContraint) ? missing : minimum(maximum(bI_ratio[globalContraint,:],dims=2))
        # Bdisposal Distance function
        gIConstraint = all(gI₀'  .>=   gI, dims=2)
        gOConstraint = all(gO₀'  .<=   gO, dims=2)
        bOConstraint = all(bO₀'  .>=   bO, dims=2)
        globalContraint = dropdims(all(hcat(gIConstraint,gOConstraint,bOConstraint), dims=2), dims=2)
        effScore_bfrontier = ! any(globalContraint) ? missing : minimum(maximum(bI_ratio[globalContraint,:],dims=2))
        effscore = max(effScore_normal,effScore_bfrontier)
        return effscore
    elseif  directions == (0,1,0,0) && prodStructure == "additive"
        bI_ratio  =  1 .- (bI ./ bI₀') # nDMU x ngI
        # Normal dist function
        gIConstraint = all(gI₀'  .>=   gI, dims=2)
        gOConstraint = all(gO₀'  .<=   gO, dims=2)
        bOConstraint = all(bO₀'  .<=   bO, dims=2)
        globalContraint = dropdims(all(hcat(gIConstraint,gOConstraint,bOConstraint), dims=2), dims=2)
        effScore_normal = ! any(globalContraint) ? missing : maximum(minimum(bI_ratio[globalContraint,:],dims=2))
        # Bdisposal Distance function
        gIConstraint = all(gI₀'  .>=   gI, dims=2)
        gOConstraint = all(gO₀'  .<=   gO, dims=2)
        bOConstraint = all(bO₀'  .>=   bO, dims=2)
        globalContraint = dropdims(all(hcat(gIConstraint,gOConstraint,bOConstraint), dims=2), dims=2)
        effScore_bfrontier = ! any(globalContraint) ? missing : maximum(minimum(bI_ratio[globalContraint,:],dims=2))
        effscore = min(effScore_normal,effScore_bfrontier)
        return effscore
    elseif  directions == (0,0,1,0) && prodStructure == "multiplicative"
        gO_ratio  =  gO₀' ./ gO  # nDMU x ngI
        # Normal dist function
        gIConstraint = all(gI₀'  .>=   gI, dims=2)
        bIConstraint = all(bI₀'  .>=   bI, dims=2)
        bOConstraint = all(bO₀'  .<=   bO, dims=2)
        globalContraint = dropdims(all(hcat(gIConstraint,bIConstraint,bOConstraint), dims=2), dims=2)
        effScore_normal = ! any(globalContraint) ? missing : minimum(maximum(gO_ratio[globalContraint,:],dims=2))
        # Bdisposal Distance function
        gIConstraint = all(gI₀'  .>=   gI, dims=2)
        bIConstraint = all(bI₀'  .<=   bI, dims=2)
        bOConstraint = all(bO₀'  .>=   bO, dims=2)
        globalContraint = dropdims(all(hcat(gIConstraint,bIConstraint,bOConstraint), dims=2), dims=2)
        effScore_bfrontier = ! any(globalContraint) ? missing : minimum(maximum(gO_ratio[globalContraint,:],dims=2))
        effscore = max(effScore_normal,effScore_bfrontier)
        return 1/effscore
    elseif  directions == (0,0,1,0) && prodStructure == "additive"
        gO_ratio  =  (gO₀' ./ gO ) .- 1 # nDMU x ngI
        # Normal dist function
        gIConstraint = all(gI₀'  .>=   gI, dims=2)
        bIConstraint = all(bI₀'  .>=   bI, dims=2)
        bOConstraint = all(bO₀'  .<=   bO, dims=2)
        globalContraint = dropdims(all(hcat(gIConstraint,bIConstraint,bOConstraint), dims=2), dims=2)
        effScore_normal = ! any(globalContraint) ? missing : maximum(minimum(gO_ratio[globalContraint,:],dims=2)) #todo check this
        # Bdisposal Distance function
        gIConstraint = all(gI₀'  .>=   gI, dims=2)
        bIConstraint = all(bI₀'  .<=   bI, dims=2)
        bOConstraint = all(bO₀'  .>=   bO, dims=2)
        globalContraint = dropdims(all(hcat(gIConstraint,bIConstraint,bOConstraint), dims=2), dims=2)
        effScore_bfrontier = ! any(globalContraint) ? missing : maximum(minimum(gO_ratio[globalContraint,:],dims=2)) # todo check this
        effscore = min(effScore_normal,effScore_bfrontier)
        return -effscore
    elseif  directions == (0,0,0,-1) && prodStructure == "multiplicative"
        bO_ratio  =  bO ./ bO₀' # nDMU x ngI
        # Normal dist function
        gIConstraint = all(gI₀'  .>=   gI, dims=2)
        bIConstraint = all(bI₀'  .>=   bI, dims=2)
        gOConstraint = all(gO₀'  .<=   gO, dims=2)
        globalContraint = dropdims(all(hcat(gIConstraint,bIConstraint,gOConstraint), dims=2), dims=2)
        effScore_normal = ! any(globalContraint) ? missing : minimum(minimum(bO_ratio[globalContraint,:],dims=2))
        # Bdisposal Distance function
        gIConstraint = all(gI₀'  .>=   gI, dims=2)
        bIConstraint = all(bI₀'  .<=   bI, dims=2)
        gOConstraint = all(gO₀'  .<=   gO, dims=2)
        globalContraint = dropdims(all(hcat(gIConstraint,bIConstraint,gOConstraint), dims=2), dims=2)
        effScore_bfrontier = ! any(globalContraint) ? missing : minimum(maximum(bO_ratio[globalContraint,:],dims=2))
        effscore = max(effScore_normal,effScore_bfrontier)
        return 1/effscore
    elseif  directions == (0,0,0,1) && prodStructure == "additive"
        bO_ratio  =  1 .- (bO ./ bO₀') # nDMU x ngI
        # Normal dist function
        gIConstraint = all(gI₀'  .>=   gI, dims=2)
        bIConstraint = all(bI₀'  .>=   bI, dims=2)
        gOConstraint = all(gO₀'  .<=   gO, dims=2)
        globalContraint = dropdims(all(hcat(gIConstraint,bIConstraint,gOConstraint), dims=2), dims=2)
        effScore_normal = ! any(globalContraint) ? missing : maximum(minimum(bO_ratio[globalContraint,:],dims=2)) #todo check this
        # Bdisposal Distance function
        gIConstraint = all(gI₀'  .>=   gI, dims=2)
        bIConstraint = all(bI₀'  .<=   bI, dims=2)
        gOConstraint = all(gO₀'  .<=   gO, dims=2)
        globalContraint = dropdims(all(hcat(gIConstraint,bIConstraint,gOConstraint), dims=2), dims=2)
        effScore_bfrontier = ! any(globalContraint) ? missing : minimum(minimum(bO_ratio[globalContraint,:];dims=2)) # todo check this
        effscore = min(effScore_normal,effScore_bfrontier)
        return effscore
    elseif  directions == (0,0,1,-1) && prodStructure == "multiplicative"
        gI_ratio  =  gI₀' ./ gI
        bI_ratio  =  bI₀' ./ bI
        gO_ratio  =  gO ./ gO₀'
        bO_ratio  =  bO ./ bO₀'

        if retToScale == "constant"
            normalFrontierDistances   = [minimum(gI_ratio[z,:]) * minimum(hcat(gO_ratio,bO_ratio)[z,:]) for z in 1:nDMUs]
            bFrontierDistances        = [minimum(gI_ratio[z,:]) * minimum(gO_ratio[z,:]) * (minimum(gO_ratio[z,:]) >= minimum(bO_ratio[z,:])) for z in 1:nDMUs]
        else
            normalFrontierDistances   = [minimum(gI_ratio[z,:]) * minimum(hcat(gO_ratio,bO_ratio)[z,:] * ( minimum(gI_ratio[z,:]) >= (1.0 - eps()) ) ) for z in 1:nDMUs]
            bFrontierDistances        = [minimum(gI_ratio[z,:]) * minimum(gO_ratio[z,:]) * (minimum(gO_ratio[z,:]) >= minimum(bO_ratio[z,:]) &&  minimum(gI_ratio[z,:]) >= (1.0 - eps())    ) for z in 1:nDMUs]
        end
        return min(maximum(normalFrontierDistances),maximum(bFrontierDistances))

    else
        @error "Directions not supported"
        return missing
    end

end # end function

"""
     dmuEfficiency(I₀,O₀,I,O)

Compute the efficiency score for a DMU using vanilla Data Envelope Analysis.

## Parameters:
  *  `I₀`: This DMU inputs (nI)
  *  `O₀`: This DMU outputs (n0)
  *  `I`: All DMUs inputs (nDMU x nI)
  *  `O`: All DMUs outputs (nDMU x n0)

## Returns:
- A named tuple with:
  - `computed`: a boolean to indicase wheter the underlying optimisation succeeded
  - `eff`: a boolean to indicate if the DMU is efficient or not
  - `obj`: a scalar representing the efficiency score of the DMU
  - `wI`: a vector of the optimal input weights
  - `wO`: a vector of the optimal output weights
  - `refSet`: a dictionary having as key the numeral of the reference DMUs and
  values the relative constraint duals
"""
function dmuEfficiency(I₀,O₀,I,O) # # Cooper & oth., p.23-24 and 43
    (nDMU,nI,nO) = size(I,1), size(I,2), size(O,2)

    effmodel = Model(optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => GLPK.GLP_MSG_OFF)) # we choose GLPK with a null output
    # Defining variables
    @variables effmodel begin
        wI[i in 1:nI] >= 0  # inputs weigths
        wO[o in 1:nO] >= 0  # outputs weigths
    end

    # Defining constraints
    @constraints effmodel begin
        othDMUsEffConstr[d in 1:nDMU],   #
           dot(O[d,:], wO) - dot(I[d,:], wI) <=  0.0
        iRegularisationConstr, # regularisation of inputs
            dot(I₀,wI) == 1.0

    end

    # Defining objective
    @objective effmodel Max begin
         dot(O₀, wO)
    end

    # Printing the model (for debugging)
    #print(effmodel) # The model in mathematical terms is printed

    # Solving the model
    optimize!(effmodel)

    status = termination_status(effmodel)

    if (status == MOI.OPTIMAL || status == MOI.LOCALLY_SOLVED) && has_values(effmodel)
        wI = value.(wI)
        wO = value.(wO)
        othDMUsEffConstrDuals = -dual.(othDMUsEffConstr)
        iRegConstrDual = -dual(iRegularisationConstr)
        obj = dot(O₀,wO) / dot(I₀,wI)
        refSet = Dict{Int64,Float64}()
        [refSet[i] = f for (i,f) in enumerate(othDMUsEffConstrDuals) if f != 0.0]
        eff  = (isapprox.(obj,1.0,atol=1e-15) && !any(isapprox.(wI,0.0,atol=1e-15)) && !any(isapprox.(wI,0.0,atol=1e-15))) # efficient if the objective is 1 AND no weigths are zero (otherwise it is on one periphal side of the frontier and shadowed by other points)
        return (computed=true,eff=eff, obj=obj, wI=wI, wO = wO, refSet=refSet, othDMUsEffConstrDuals=othDMUsEffConstrDuals,iRegConstrDual=iRegConstrDual)
    else
        return (computed=false,eff=missing,obj=missing,wI=missing,wO = missing,refSet=missing,othDMUsEffConstrDuals=missing,iRegConstrDual=missing)
    end
end

# Cooper & oth., p.43
"""
     dmuEfficiencyDual(I₀,O₀,I,O)

Compute the efficiency score for a DMU using vanilla Data Envelope Analysis in the dual representation using the two pass method.

## Parameters:
  *  `I₀`: This DMU inputs (nI)
  *  `O₀`: This DMU outputs (n0)
  *  `I`: All DMUs inputs (nDMU x nI)
  *  `O`: All DMUs outputs (nDMU x n0)

## Returns:
- A named tuple with:
  - `computed`: a boolean to indicase wheter the underlying optimisation succeeded
  - `eff`: a boolean to indicate if the DMU is efficient or not
  - `obj`: a scalar representing the efficiency score of the DMU
  - `λ`
  - `θ`
  - `dualsI`
  - `dualsO`
  - `s⁻`
  - `s⁺`
"""
function dmuEfficiencyDual(I₀,O₀,I,O)
    (nDMU,nI,nO) = size(I,1), size(I,2), size(O,2)

    effmodeldual = Model(optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => GLPK.GLP_MSG_OFF)) # we choose GLPK with a verbose output
    # Defining variables
    @variables effmodeldual begin
        λ[i in 1:nDMU] >= 0  # the constrain dual in the primal
        θ                    # the costs per output ???
    end

    # Defining constraints
    @constraints effmodeldual begin
        inConstr[i in 1:nI],   #
           θ*I₀[i] - dot(I[:,i], λ) >=  0.0
        oRegularisationConstr[j in 1:nO], # regularisation of outputs ??
            dot(O[:,j],λ) >= O₀[j]
    end

    # Defining objective
    @objective effmodeldual Min begin
         θ
    end

    # Printing the model (for debugging)
    #print(effmodel) # The model in mathematical terms is printed

    # Solving the model
    optimize!(effmodeldual)

    status = termination_status(effmodeldual)

    if (status == MOI.OPTIMAL || status == MOI.LOCALLY_SOLVED) && has_values(effmodeldual)
        λv = value.(λ)
        θv = value(θ)
        dualsI = dual.(inConstr)
        dualsO = dual.(oRegularisationConstr)
        obj = objective_value(effmodeldual)
        outPass2 = dmuPass2(θv,I₀,O₀,I,O)
        s⁻ = outPass2.s⁻ # input excesses
        s⁺ = outPass2.s⁺ # output shortfalls
        eff  = (isapprox.(obj,1.0,atol=1e-15) && all(isapprox.(s⁻,0.0,atol=1e-15)) && all(isapprox.(s⁺,0.0,atol=1e-15))) # efficient if the objective is 1 AND no weigths are zero (otherwise it is on one periphal side of the frontier and shadowed by other points)
        #refSet = Dict{Int64,Float64}()
        #[refSet[i] = -f for (i,f) in enumerate(duals) if f != 0.0]
        #eff  = (isapprox.(obj,1.0,atol=1e-15) && !any(isapprox.(wI,0.0,atol=1e-15)) && !any(isapprox.(wI,0.0,atol=1e-15))) # efficient if the objective is 1 AND no weigths are zero (otherwise it is on one periphal side of the frontier and shadowed by other points)
        #return (eff=eff, obj=obj, wI=wI, wO = wO, refSet=refSet)
        return (computed=true, eff=eff, obj=obj,  λ=λv, θ=θv, dualsI=dualsI, dualsO=dualsO, s⁻=s⁻, s⁺=s⁺)
    else
        #return (eff=missing,obj=missing,wI=missing,wO = missing,refSet=missing)
        return (computed=false, eff=eff, obj=missing, λ=missing, θ=missing, dualsI=missing , dualsO=missing , s⁻=missing, s⁺=missing)
    end
end

# Cooper & oth., p.44
function dmuPass2(θ,I₀,O₀,I,O)
    (nDMU,nI,nO) = size(I,1), size(I,2), size(O,2)

    effmodeldual = Model(optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => GLPK.GLP_MSG_OFF)) # we choose GLPK with a verbose output
    # Defining variables
    @variables effmodeldual begin
        λ[i in 1:nDMU] >= 0  # the constrain dual in the primal
        s⁻[i in 1:nI]  >= 0  # input excesses
        s⁺[i in 1:nO]  >= 0  # output shortfalls
    end

    # Defining constraints
    @constraints effmodeldual begin
        inConstr[i in 1:nI],   #
           θ*I₀[i] - dot(I[:,i], λ) ==  s⁻[i]
        oRegularisationConstr[j in 1:nO], # regularisation of outputs ??
            dot(O[:,j],λ) - O₀[j] == s⁺[j]
    end

    # Defining objective
    @objective effmodeldual Max begin
         sum(s⁻)+sum(s⁺)
    end

    # Printing the model (for debugging)
    #print(effmodel) # The model in mathematical terms is printed

    # Solving the model
    optimize!(effmodeldual)

    status = termination_status(effmodeldual)

    if (status == MOI.OPTIMAL || status == MOI.LOCALLY_SOLVED) && has_values(effmodeldual)
        λv = value.(λ)
        s⁻v = value.(s⁻)
        s⁺v = value.(s⁺)
        obj = objective_value(effmodeldual)
        return (computed=true, obj=obj, λ=λv, s⁻=s⁻v, s⁺=s⁺v )
    else
        #return (eff=missing,obj=missing,wI=missing,wO = missing,refSet=missing)
        return (computed=false, obj=missing, λ=missing,s⁻=missing, s⁺=missing )
    end
end


text:
module Libxkb

import xkbcommon_jll

using CEnum

const libxkbcommon = xkbcommon_jll.libxkbcommon
const libxkbcommon_x11 = first(filter( x -> occursin("libxkbcommon-x11", basename(x)), readdir(joinpath(xkbcommon_jll.artifact_dir, "lib"), join=true)))

const FILE=Ptr{Cvoid}

include("xcb_common.jl")
include("xkb_common.jl")
include("xkb_api.jl")

# export everything
foreach(names(@__MODULE__, all=true)) do s
    if startswith(string(s), r"(?:X|XKB|xkb)") && !startswith(string(s), r"(?:XCB|xcb)")
        @eval export $s
    end
end

end


text:
# ---
# title: 149. Max Points on a Line
# id: problem149
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Hash Table, Math
# link: <https://leetcode.com/problems/max-points-on-a-line/description/>
# hidden: true
# ---
#
# Given _n_ points on a 2D plane, find the maximum number of points that lie on
# the same straight line.
#
# **Example 1:**
#
#
#
#     Input: [[1,1],[2,2],[3,3]]
#     Output: 3
#     Explanation:
#     ^
#     |
#     |        o
#     |     o
#     |  o
#     +------------->
#     0  1  2  3  4
#
#
# **Example 2:**
#
#
#
#     Input: [[1,1],[3,2],[5,3],[4,1],[2,3],[1,4]]
#     Output: 4
#     Explanation:
#     ^
#     |
#     |  o
#     |     o        o
#     |        o
#     |  o        o
#     +------------------->
#     0  1  2  3  4  5  6
#
#
# **NOTE:**  input types have been changed on April 15, 2019. Please reset to
# default code definition to get new method signature.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Codecs
using Test

include("tests.jl")


text:
# ------------ #
# Original API #
# ------------ #

# chebdef.m -- DONE
function chebdef(n::Int, a::Real, b::Real)
    p = ChebParams(n, a, b)
    return (p.n, p.a, p.b, Any[p.n, p.a, p.b], p)
end

# chebnode.m -- DONE
chebnode(n::Int, a::Real, b::Real, nodetype=0) =
    nodes(ChebParams(n, a, b), nodetype)

# chebdop.m -- DONE
function chebdop(n, a, b, order=1)
    D, p = derivative_op(ChebParams(n, a, b), order)
    D, p.n-order, p.a, p.b, Any[p.n-order, p.a, p.b]
end

# chebbase.m -- DONE
function chebbase(n, a, b, x=chebnode(n, a, b, 1), order=0, nodetype=1)
    B = evalbase(ChebParams(n, a, b), x, order, nodetype)
    B, x
end

# chebbasex.m -- DONE
chebbasex(n, a, b, x) = evalbasex(ChebParams(n, a, b), x)


text:
function f_model(
    f_x::R,
    grad_f_x::AbstractArray{R},
    res::AbstractArray{R},
    gamma::R,
) where {R<:Real}
    return f_x - dot(grad_f_x, res) + (0.5 / gamma) * norm(res)^2
end


text:
module CMIM

export cmim

function jointProbability(x::AbstractArray{Bool}, y::AbstractArray{Bool})

    counts = zeros(Int, 2, 2)

    for ii in 1:length(x)
        counts[x[ii]+1, y[ii]+1] += 1
    end

    counts ./ length(x)

end

function jointProbability(x::AbstractArray{Bool}, y::AbstractArray{Bool}, z::AbstractArray{Bool})

    counts = zeros(Int, 2, 2, 2)

    for ii in 1:length(x)
        counts[x[ii]+1, y[ii]+1, z[ii]+1] += 1
    end

    counts ./ length(x)

end

function mutualInformation(x::AbstractArray{Bool}, y::AbstractArray{Bool})

    pxy = jointProbability(x, y)

    mutualInformation = 0.0

    for xx in 1:2, yy in 1:2
        if pxy[xx,yy] > 0.0
            px = sum(pxy[xx,:])
            py = sum(pxy[:,yy])
            mutualInformation += pxy[xx,yy] * log(pxy[xx,yy]/px/py)
        end
    end

    mutualInformation

end

function conditionalEntropy(x::AbstractArray{Bool}, y::AbstractArray{Bool})

    pxy = jointProbability(x, y)

    conditionalEtropy = 0.0

    for xx in 1:2, yy in 1:2
        if pxy[xx,yy] > 0.0
            py = sum(pxy[:,yy])
            conditionalEtropy -= pxy[xx,yy] * log(pxy[xx,yy]/py)
        end
    end

    conditionalEtropy

end

function conditionalEntropy(x::AbstractArray{Bool}, y::AbstractArray{Bool}, z::AbstractArray{Bool})

    pxyz = jointProbability(x, y, z)

    conditionalEntropy = 0.0

    for xx in 1:2, yy in 1:2, zz in 1:2
        if pxyz[xx,yy,zz] > 0.0
            pyz = sum(pxyz[:,yy,zz])
            conditionalEntropy -= pxyz[xx,yy,zz] * log(pxyz[xx,yy,zz]/pyz)
        end
    end

    conditionalEntropy

end

function conditionalMutualInformation(x::AbstractArray{Bool}, y::AbstractArray{Bool}, z::AbstractArray{Bool})

    conditionalEntropy(x, z) - conditionalEntropy(x, y, z)
#     conditionalEntropy(y, z) - conditionalEntropy(y, x, z)

end

function conditionalMutualInformation(cacheConditionalEntropy::Vector, X::AbstractArray{Bool,2}, x::Int, y::AbstractArray{Bool,1}, z::Int)

    if cacheConditionalEntropy[z] < 0.0
        cacheConditionalEntropy[z] = conditionalEntropy(y, X[:,z])
    end

    cacheConditionalEntropy[z] - conditionalEntropy(y, X[:,x], X[:,z])

end

function cmim(k::Int, X::AbstractArray{Bool,2}, y::AbstractArray{Bool,1})

    @assert k <= size(X, 2)
    @assert size(X, 1) == length(y)

    outputFeatures = Vector{Int}(k)
    outputScores = Vector{Float64}(k)

    selectedFeatures = falses(size(X, 2))
    lastUsedFeature = ones(Int, size(X, 2))

    cacheConditionalEntropy = -ones(size(X, 2))

    classMI = [mutualInformation(X[:,ff], y) for ff in 1:size(X, 2)]

    for kk in 1:k
        score = -1.0
        for ff in 1:size(X, 2)
            if !selectedFeatures[ff]

                while (classMI[ff] > score) && (lastUsedFeature[ff] < kk)
                    classMI[ff] = min(
                        classMI[ff],
#                         conditionalMutualInformation(X[:,ff], y, X[:,outputFeatures[lastUsedFeature[ff]]])
                        conditionalMutualInformation(cacheConditionalEntropy, X, ff, y, outputFeatures[lastUsedFeature[ff]])
                    )
                    lastUsedFeature[ff] += 1
                end

                if classMI[ff] > score
                    score = classMI[ff]
                    outputFeatures[kk] = ff
                    outputScores[kk] = classMI[ff]
                end

            end
        end
        selectedFeatures[outputFeatures[kk]] = true
    end

    collect(zip(outputFeatures, outputScores))

end

end # module


text:
# function to normalize values
function equalize(vec)
    map = Array{Int64,2}(undef, length(unique(vec)),2) .= 0
    vec = round.(Int64, vec)
    map[:,1] = unique(vec)
    for i = 1:size(map,1)
        map[i,2] = i
    end
    for i = 1:size(vec,1)
        for j = 1:size(map,1)
            if vec[i] == map[j,1]
                vec[i] = map[j,2]
            end
        end
    end
    normalize(vec)
    return vec
end


text:
#=
I used to define AL/ARview setters (and still can if needed), but the problem is that they're generically not really defined
our tensor decompositions are only unique when some conditions are satisfied
Therefore I've decided to remove them
=#
struct ALView{S}
    parent::S
end

function Base.getindex(v::ALView{S},i::Int)::eltype(S) where S
    ismissing(v.parent.ALs[i]) && v.parent.CR[i] # by getting CL[i+1], we are garantueeing that AL[i] exists
    return v.parent.ALs[i]
end

struct ARView{S}
    parent::S
end

function Base.getindex(v::ARView{S},i::Int)::eltype(S) where S
    ismissing(v.parent.ARs[i]) && v.parent.CR[i-1] # by getting CL[i], we are garantueeing that AR[i] exists
    return v.parent.ARs[i]
end

struct CRView{S}
    parent::S
end

function Base.getindex(v::CRView{S},i::Int)::bond_type(S) where S
    if ismissing(v.parent.CLs[i+1])
        if i == 0 || !ismissing(v.parent.ALs[i])
            (v.parent.CLs[i+1],temp) = rightorth(_permute_tail(v.parent.AC[i+1]),alg=LQpos())
            v.parent.ARs[i+1] = _permute_front(temp);
        else
            (v.parent.ALs[i],v.parent.CLs[i+1]) = leftorth(v.parent.AC[i],alg=QRpos())
        end
    end
    return v.parent.CLs[i+1]
end

function Base.setindex!(v::CRView,vec,i::Int)
    if ismissing(v.parent.CLs[i+1])
        if !ismissing(v.parent.ALs[i])
            (v.parent.CLs[i+1],temp) = rightorth(_permute_tail(v.parent.AC[i+1]),alg=LQpos())
            v.parent.ARs[i+1] = _permute_front(temp);
        else
            (v.parent.ALs[i],v.parent.CLs[i+1]) = leftorth(v.parent.AC[i],alg=QRpos())
        end
    end

    v.parent.CLs.=missing;
    v.parent.ACs.=missing;
    v.parent.ALs[i+1:end].=missing;
    v.parent.ARs[1:i].=missing;

    v.parent.CLs[i+1] = vec
end

struct ACView{S}
    parent::S
end

function Base.getindex(v::ACView{S},i::Int)::eltype(S) where S
    if ismissing(v.parent.ACs[i]) && !ismissing(v.parent.ARs[i])
        c = v.parent.CR[i-1];
        ar = v.parent.ARs[i];
        v.parent.ACs[i] = _permute_front(c*_permute_tail(ar))
    elseif ismissing(v.parent.ACs[i]) && !ismissing(v.parent.ALs[i])
        c = v.parent.CR[i];
        al = v.parent.ALs[i];
        v.parent.ACs[i] = al*c;
    end
    return v.parent.ACs[i]
end

function Base.setindex!(v::ACView,vec::GenericMPSTensor,i::Int)
    if ismissing(v.parent.ACs[i])
        i<length(v) && v.parent.AR[i+1]
        i>1 && v.parent.AL[i-1]
    end

    v.parent.ACs.=missing;
    v.parent.CLs.=missing;
    v.parent.ALs[i:end].=missing;
    v.parent.ARs[1:i].=missing;
    v.parent.ACs[i] = vec;
end

function Base.setindex!(v::ACView,vec::Tuple{<:GenericMPSTensor,<:GenericMPSTensor},i::Int)
    if ismissing(v.parent.ACs[i])
        i<length(v) && v.parent.AR[i+1]
        i>1 && v.parent.AL[i-1]
    end

    v.parent.ACs.=missing;
    v.parent.CLs.=missing;
    v.parent.ALs[i:end].=missing;
    v.parent.ARs[1:i].=missing;

    (a,b) = vec
    if isa(a,MPSBondTensor) #c/ar
        v.parent.CLs[i] = a;
        v.parent.ARs[i] = b;
    else #al/c
        @assert isa(b,MPSBondTensor)

        v.parent.CLs[i+1] = b;
        v.parent.ALs[i] = a;
    end

end

Base.firstindex(psi::Union{ACView,ALView,ARView}, i...) = firstindex(psi.parent.ALs, i...)
Base.lastindex(psi::Union{ACView,ALView,ARView}, i...) = lastindex(psi.parent.ALs, i...)
Base.length(psi::Union{ACView,ALView,ARView}) = length(psi.parent);
Base.size(psi::Union{ACView,ALView,ARView},args...) = size(psi.parent,args...);

Base.firstindex(psi::Union{CRView}) = 0;
Base.lastindex(psi::Union{CRView}, i...) = lastindex(psi.parent.ALs, i...)
Base.length(psi::Union{CRView}) = length(psi.parent.CLs);
Base.size(psi::Union{CRView},args...) = size(psi.parent.CLs,args...);

Base.IteratorSize(::Type{<:Union{ACView,ALView,ARView,CRView}}) = Base.HasShape{1}()
Base.IteratorEltype(::Type{<:Union{ACView,ALView,ARView,CRView}}) = Base.HasEltype()
Base.iterate(view::Union{ACView,ALView,ARView,CRView},istate = firstindex(view)) = istate > lastindex(view) ? nothing : (view[istate],istate+1)

Base.getindex(psi::Union{ACView,ALView,ARView,CRView},r::AbstractRange{Int64}) = [psi[ri] for ri in r]