text:
using Crazyflie
using ProgressMeter

# temporary
using PyCall
using DelimitedFiles

figure8 = [
    [1.050000, 0.000000, -0.000000, 0.000000, -0.000000, 0.830443, -0.276140, -0.384219, 0.180493, -0.000000, 0.000000, -0.000000, 0.000000, -1.356107, 0.688430, 0.587426, -0.329106, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000],  # noqa
    [0.710000, 0.396058, 0.918033, 0.128965, -0.773546, 0.339704, 0.034310, -0.026417, -0.030049, -0.445604, -0.684403, 0.888433, 1.493630, -1.361618, -0.139316, 0.158875, 0.095799, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000],  # noqa
    [0.620000, 0.922409, 0.405715, -0.582968, -0.092188, -0.114670, 0.101046, 0.075834, -0.037926, -0.291165, 0.967514, 0.421451, -1.086348, 0.545211, 0.030109, -0.050046, -0.068177, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000],  # noqa
    [0.700000, 0.923174, -0.431533, -0.682975, 0.177173, 0.319468, -0.043852, -0.111269, 0.023166, 0.289869, 0.724722, -0.512011, -0.209623, -0.218710, 0.108797, 0.128756, -0.055461, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000],  # noqa
    [0.560000, 0.405364, -0.834716, 0.158939, 0.288175, -0.373738, -0.054995, 0.036090, 0.078627, 0.450742, -0.385534, -0.954089, 0.128288, 0.442620, 0.055630, -0.060142, -0.076163, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000],  # noqa
    [0.560000, 0.001062, -0.646270, -0.012560, -0.324065, 0.125327, 0.119738, 0.034567, -0.063130, 0.001593, -1.031457, 0.015159, 0.820816, -0.152665, -0.130729, -0.045679, 0.080444, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000],  # noqa
    [0.700000, -0.402804, -0.820508, -0.132914, 0.236278, 0.235164, -0.053551, -0.088687, 0.031253, -0.449354, -0.411507, 0.902946, 0.185335, -0.239125, -0.041696, 0.016857, 0.016709, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000],  # noqa
    [0.620000, -0.921641, -0.464596, 0.661875, 0.286582, -0.228921, -0.051987, 0.004669, 0.038463, -0.292459, 0.777682, 0.565788, -0.432472, -0.060568, -0.082048, -0.009439, 0.041158, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000],  # noqa
    [0.710000, -0.923935, 0.447832, 0.627381, -0.259808, -0.042325, -0.032258, 0.001420, 0.005294, 0.288570, 0.873350, -0.515586, -0.730207, -0.026023, 0.288755, 0.215678, -0.148061, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000],  # noqa
    [1.053185, -0.398611, 0.850510, -0.144007, -0.485368, -0.079781, 0.176330, 0.234482, -0.153567, 0.447039, -0.532729, -0.855023, 0.878509, 0.775168, -0.391051, -0.713519, 0.391628, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000],  # noqa
]

heart = [
[0.608998,0,0,0,0,8.42094,-20.1848,17.4023,-5.13733,0,0,0,0,-5.25294,19.5805,-26.0513,12.0229,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0.779122,0.195632,0.637789,0.0169511,-0.469071,1.13898,-5.83124,8.3787,-3.63444,-0.0377946,-0.0794284,-0.0646026,0.0244159,-3.78373,9.99798,-9.63602,3.27284,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0.757723,0.467339,-0.036739,-0.110484,0.0639573,0.108074,-1.39392,2.11362,-0.935232,-0.236257,-0.347344,-0.032209,-0.0141377,-1.22068,4.0166,-4.3725,1.62188,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0.957797,0.357255,-0.274841,-0.111489,-0.0115323,-3.42962,8.38992,-7.1466,2.09672,-0.518182,-0.328702,0.0913417,0.0150219,0.733687,-1.05708,0.479298,-0.0521867,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0.95006,-0.108941,-0.502428,-0.0278039,-0.0132172,-2.03991,5.51081,-5.00167,1.53592,-0.639182,0.0542293,0.051927,0.014188,3.04728,-7.84289,6.93303,-2.09398,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0.998431,-0.624157,-0.461845,-0.0138318,-0.00967601,-2.04884,5.5933,-4.99039,1.48711,-0.481169,0.114548,0.029441,0.000554473,3.31264,-7.6003,6.14741,-1.71994,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0.812602,-1.06789,-0.279874,-0.0292519,-0.0503488,-3.7103,10.0379,-9.64748,3.24897,-0.197235,0.269061,-0.00363612,0.0151035,0.916496,-3.79814,4.50747,-1.72174,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0.419508,-1.42038,-0.468108,0.0856745,0.221818,3.92471,-0.277235,-23.883,26.7762,-0.0240585,0.0921,-0.0459625,0.0111955,-3.86908,19.5572,-35.3542,22.4094,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0.64054,-1.5363,-0.00325354,0.514808,0.0803983,-2.16099,9.25223,-14.0528,6.98868,0.000127304,0.0443542,0.0111982,0.0261656,1.57639,-5.36123,6.53792,-2.78901,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0.893144,-1.33358,0.558221,0.00558749,-0.182813,1.17326,-3.86353,4.22095,-1.49454,0.05546,0.116935,0.0458194,0.000332267,2.56277,-6.40182,5.67683,-1.7497,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[1.0137,-0.944994,0.264934,0.0564772,0.0488614,3.09438,-7.21848,5.83434,-1.62251,0.277414,0.275656,0.00454676,0.0105108,1.06699,-2.79915,2.41189,-0.698664,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0.941904,-0.480716,0.458841,0.0299512,0.00400315,2.13691,-5.30337,4.62042,-1.38829,0.551501,0.16792,-0.0455683,-0.000158787,-0.00812595,-0.31599,0.457477,-0.175203,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0.902676,0.0448938,0.563892,0.0356819,-0.0204628,3.81762,-12.3689,12.5856,-4.21939,0.632664,-0.0105303,-0.0563382,-0.0149921,-4.47262,11.1707,-9.8176,2.99179,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0.696637,0.437821,0.0271313,-0.0955769,0.0683252,-3.88281,13.7328,-16.7036,6.91373,0.441299,-0.197977,0.0795611,-0.0410307,-17.0929,56.5153,-65.8397,26.5171,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0.790226,0.413439,0.0024842,-0.0197612,-0.0573184,-9.94595,26.5283,-25.7537,8.8682,0.160997,-0.308277,0.0562815,0.00236697,1.27222,-2.73425,2.12724,-0.57832,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0.418106,0.106269,-0.544119,0.234554,1.02684,3.93774,-3.17014,-30.7186,40.0474,0.0139615,-0.0875489,0.0363838,0.0507818,5.39905,-25.4338,43.0935,-25.8168,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
]

py"""
import time
from cflib.crazyflie.mem import MemoryElement
from cflib.crazyflie.mem import Poly4D

class Uploader:
    def __init__(self):
        self._is_done = False

    def upload(self, trajectory_mem):
        print('Uploading data')
        trajectory_mem.write_data(self._upload_done)

        while not self._is_done:
            time.sleep(0.2)

    def _upload_done(self, mem, addr):
        print('Data uploaded')
        self._is_done = True

def upload_trajectory(cf, trajectory_id, trajectory):
    trajectory_mem = cf.mem.get_mems(MemoryElement.TYPE_TRAJ)[0]

    total_duration = 0
    for row in trajectory:
        duration = row[0]
        x = Poly4D.Poly(row[1:9])
        y = Poly4D.Poly(row[9:17])
        z = Poly4D.Poly(row[17:25])
        yaw = Poly4D.Poly(row[25:33])
        trajectory_mem.poly4Ds.append(Poly4D(duration, x, y, z, yaw))
        total_duration += duration

    Uploader().upload(trajectory_mem)
    cf.high_level_commander.define_trajectory(trajectory_id, 0,
                                              len(trajectory_mem.poly4Ds))
    return total_duration
"""

function reset_estimator(cf)
    cf.param.set_value("kalman.resetEstimation", "1")
    sleep(0.1)
    cf.param.set_value("kalman.resetEstimation", "0")
    lconfig = logger.LogConfig(name="Kalman Variance", period_in_ms=500)
    lconfig.add_variable("kalman.varPX", "float")
    lconfig.add_variable("kalman.varPY", "float")
    lconfig.add_variable("kalman.varPZ", "float")
    varx = ones(10)*1000
    vary = ones(10)*1000
    varz = ones(10)*1000
    record(cf, lconfig) do logs
        for entry in logs
            data = entry[2]
            push!(varx, data["kalman.varPX"])
            push!(vary, data["kalman.varPY"])
            push!(varz, data["kalman.varPZ"])
            popfirst!(varx); popfirst!(vary); popfirst!(varz)
            threshx = (maximum(varx) - minimum(varx)) < 0.001
            threshy = (maximum(vary) - minimum(vary)) < 0.001
            threshz = (maximum(varz) - minimum(varz)) < 0.001
            if threshx && threshy && threshz
                break
            end
        end
    end
end

function enable_commander(cf)
    cf.param.set_value("commander.enHighLevel", "1")
    cf.high_level_commander
end

# because commander.land is not reliable
function land(commander, duration)
    print("Landing\n")
    for i = 1:10
        commander.land(0.0, duration)
        sleep(0.1)
    end
    sleep(duration)
    commander.stop()
end

# because commander.takeoff is not reliable
function takeoff(commander, height, duration)
    print("Taking Off\n")
    for i = 1:5
        commander.takeoff(height, duration)
        sleep(0.1)
    end
    sleep(duration)
end

# because commander.start_trajectory is not reliable
function start_trajectory(commander, id, rel)
    print("Starting Trajectory\n")
    for i = 1:5
        commander.start_trajectory(id, relative=rel)
        sleep(0.1)
    end
    sleep(1.0)
end

function traj(uri="radio://0/80/2M")
    play(uri) do cf
        commander = enable_commander(cf)
        duration = py"upload_trajectory($cf, $1, $heart)"
        reset_estimator(cf)
        takeoff(commander, 1.0, 2.0) # height, duration
        start_trajectory(commander, 1, true)
        sleep(duration)
        land(commander, 5.0) # duration
    end
    return nothing
end


text:
## header to provide surface and context
using Cairo
c = CairoRGBSurface(256,256);
cr = CairoContext(c);

save(cr);
set_source_rgb(cr,0.8,0.8,0.8);    # light gray
rectangle(cr,0.0,0.0,256.0,256.0); # background
fill(cr);
restore(cr);

save(cr);
## original example, following here
set_line_width(cr, 40.96);
move_to(cr, 76.8, 84.48);
rel_line_to(cr, 51.2, -51.2);
rel_line_to(cr, 51.2, 51.2);
set_line_join(cr, Cairo.CAIRO_LINE_JOIN_MITER); # default
stroke(cr);

move_to(cr, 76.8, 161.28);
rel_line_to(cr, 51.2, -51.2);
rel_line_to(cr, 51.2, 51.2);
set_line_join(cr, Cairo.CAIRO_LINE_JOIN_BEVEL);
stroke(cr);

move_to(cr, 76.8, 238.08);
rel_line_to(cr, 51.2, -51.2);
rel_line_to(cr, 51.2, 51.2);
set_line_join(cr, Cairo.CAIRO_LINE_JOIN_ROUND);
stroke(cr);

## mark picture with current date
restore(cr);
move_to(cr,0.0,12.0);
set_source_rgb(cr, 0,0,0);
show_text(cr,Libc.strftime(time()));
write_to_png(c,"sample_set_line_join.png");

ext: jl


text:
module prettytbl
#=  adapted from https://github.com/djsegal/Fussy.jl/blob/
 2cde77e2422dce590e770d2937a444c1eabfea7a/src/modules/HTMLElements.jl  =#
  using ArgCheck

  function _table(cur_matrix::Matrix; header_row=[], header_col=[], title=[],
      header=[], dp::Int=Inf)

    if isfinite(dp) && dp > 0
        cur_matrix = round.(cur_matrix*10^dp)/10^dp
    end

    cur_table = ""
    !isempty(header_col) && @argcheck size(cur_matrix, 1) == length(header_col)
    !isempty(header_row) && @argcheck size(cur_matrix, 2) == length(header_row)

    if !isempty(header)
      cur_table *= "

$header

"
    end

    cur_table *= "
" #"
"
    if !isempty(title)
        cur_table *= "
$title
"
    end

    if !isempty(header_row)
      cur_table *= "

"
      !isempty(header_col) && (cur_table *= "")

      for cur_header in header_row
        cur_table *= "$cur_header	"
      end

      cur_table *= "

"
    end

    cur_table *= "
"

    for ii in 1:size(cur_matrix, 1)

      cur_table *= "
"
      if !isempty(header_col)
        cur_table *= ""
        cur_table *= header_col[ii]
      end

      for jj in 1:size(cur_matrix, 2)
        cur_element = cur_matrix[ii, jj]

        if isa(cur_element, Dict)
          cur_table *= ""
          cur_value = cur_element["value"]
        else
          cur_table *= ""
          cur_value = cur_element
        end

        cur_table *= string(cur_value)

        cur_table *= "	"
      end

      cur_table *= "
"

    end

    cur_table *= "
"

    cur_table *= "
"

    return cur_table

  end

table(cur_matrix::Matrix; header_row=[], header_col=[], title=[], dp=Inf) =
    HTML(_table(cur_matrix; header_row=header_row, header_col=header_col,
    title=title, header=[], dp=dp))

function table(cur_matrix1::Matrix, cur_matrix2::Matrix; header_row=[], header_col=[],
    title=[], header=[], dp=Inf)
    title = isempty(title) ? ["", ""] : (title isa String ? repeat([title], 2) : title)
    t1 = _table(cur_matrix1; header_row=header_row, header_col=header_col, title=title[1],
            header=[], dp=dp)
    t2 = _table(cur_matrix2; header_row=header_row, title=title[2], header=[], dp=dp)

    head = !isempty(header) ? "

$header

" : ""
    HTML(head * "\n" * t1 * "\n" * t2)
end

function table(cur_matrix1::Array, cur_matrix2::Array; header_row=[], header_col=[], title=[],
    header=[], dp=Inf)
    if cur_matrix1[1] isa Real
        cur_matrix1 = reshape(cur_matrix1, 1, :)
    elseif cur_matrix1[1] isa Array
        cur_matrix1 = reduce(vcat, [reshape(x, 1, :) for x in cur_matrix1])
    end
    table(cur_matrix1; header_row=header_row, header_col=header_col, title=title,
        header=header, dp=dp)
end

function table(cur_matrix1::Array, cur_matrix2::Array; header_row=[], header_col=[],
    title=[], header=[], dp=Inf)
    if cur_matrix1[1] isa Real
        cur_matrix1 = reshape(cur_matrix1, 1, :)
    elseif cur_matrix1[1] isa Array
        cur_matrix1 = reduce(vcat, [reshape(x, 1, :) for x in cur_matrix1])
    end

    if cur_matrix2[1] isa Real
        cur_matrix2 = reshape(cur_matrix2, 1, :)
    elseif cur_matrix2[1] isa Array
        cur_matrix2 = reduce(vcat, [reshape(x, 1, :) for x in cur_matrix2])
    end
    table(cur_matrix1, cur_matrix2; header_row=header_row, header_col=header_col,
        title=title, header=header, dp=dp)
end
end

max_stars_repo_name: ornithos/pytorch-mtds-mocap

max_issues_repo_name: ornithos/pytorch-mtds-mocap

max_forks_repo_name: ornithos/pytorch-mtds-mocap


text:
abstract type Parameter end

mutable struct BoolParameter <: Parameter
    value::Bool
    name::String
end

mutable struct NumberParameter{T <: Number} <: Parameter
    min::T
    max::T
    value::T
    name::String

    function NumberParameter{T}(min::T, max::T, value::T, name::String) where T <: Number
        if min > max
            error("invalid bounds: min > max")
        elseif value < min || value > max
            error("value out of bounds.")
        end
        new(min, max, value, name)
    end

    function NumberParameter{T}(min::T, max::T, name::String) where T <: Number
        if min > max
            error("invalid bounds: min > max")
        end
        new(min, max, rand_in(min, max), name)
    end

    function NumberParameter{T}(value::T, name::String) where T <: Number
        new(value, value, value, name)
    end
end

mutable struct EnumParameter{P <: Parameter, T <: AbstractArray} <: Parameter
    values::T
    value::Int
    name::String
    current::P

    function EnumParameter{P, T}(values::T, value::Int, name::String) where {P <: Parameter, T <: AbstractArray}
        if value > length(values) || value < 1
            error("value is out of bounds.")
        end
        new{P, T}(values, value, name, values[value])
    end

    function EnumParameter{P, T}(values::T, name::String) where {P <: Parameter, T <: AbstractArray}
        value = rand(1:length(values))
        new{P, T}(values, value, name, values[value])
    end
end

function EnumParameter(values::T, value::Int, name::String) where T <: AbstractArray
    EnumParameter{Parameter, T}(values, value, name)
end

function EnumParameter(values::T, name::String) where T <: AbstractArray
    EnumParameter{Parameter, T}(values, name)
end

mutable struct PermutationParameter{T <: AbstractArray} <: Parameter
    value::T
    size::Int
    name::String

    function PermutationParameter{T}(value::T, name::String) where T <: AbstractArray
        new(value, length(value), name)
    end
end

function PermutationParameter(value::T, name::String) where T <: AbstractArray
    PermutationParameter{T}(value, name)
end

mutable struct StringParameter <: Parameter
    value::String
    name::String
end

function StringParameter(value::String)
    StringParameter(value, value)
end

const IntegerParameter = NumberParameter{Integer}
const FloatParameter   = NumberParameter{AbstractFloat}


text:
module Pade
    using SymPy
    using Main.Wynn

    # data structure to store the pade table
    # R[l/m] indexed by ints l & m, stored within dict.
    struct PadeTable{T}
        series::T
        terms::Vector{T}
        ptable::Dict{Tuple{Int,Int},T}
    end

    # filter function, when epsilon j-value is even
    flt(tuple::Tuple{Int,Int}) = iseven(tuple[2]) && (tuple[2] >= 0) &&
    (tuple[1]+tuple[2]/2 >= 0)

    function PadeTable(terms::Vector{T}; simplified::Bool = true) where T<:Union{Real,Sym}
        # compute epsilon algorithm coefficients
        series = sum(terms)
        etable = EpsilonTable(terms, simplified = simplified).etable

        # filter away auxiliary indexes
        idxs = Iterators.filter(flt, keys(etable))
        PadeTable(series, terms, Dict((Int(i+j/2),Int(j/2)) => etable[i,j] for (i,j) in idxs))
    end

    function pade(series::Vector{T}, l::Int, m::Int) where T<:Union{Real,Sym}
        if length(series) >= l+m+1
            return PadeTable(series[1:l+m+1]).ptable[l,m]
        else
            @error "A Pade approximation of order [$l/$m] requires at least $(l+m+1) series terms."
        end
    end

    export PadeTable, pade
end


text:
"""
````
function simplexgrid(coord::Array{Tc,2},
                     cellnodes::Array{Ti,2},
                     cellregions::Array{Ti,1},
                     bfacenodes::Array{Ti,2},
                     bfaceregions::Array{Ti,1}
                     ) where {Tc,Ti}
````

    Create simplex grid from five arrays.
"""
function simplexgrid(coord::Array{Tc,2},
                     cellnodes::Array{Ti,2},
                     cellregions::Array{Ti,1},
                     bfacenodes::Array{Ti,2},
                     bfaceregions::Array{Ti,1}
                     ) where {Tc,Ti}
    @assert size(coord,2)>0
    dim=size(coord,1)
    if dim==1
        eltype=Edge1D
        btype=Vertex0D
        csys=Cartesian1D
    elseif dim==2
        eltype=Triangle2D
        btype=Edge1D
        csys=Cartesian2D
    elseif dim==3
        eltype=Tetrahedron3D
        btype=Triangle2D
        csys=Cartesian3D
    end

    grid=ExtendableGrid{Tc,Ti}()
    grid[Coordinates]=coord
    grid[CellNodes]=cellnodes
    grid[CellRegions]=cellregions
    grid[CellGeometries]=VectorOfConstants(eltype,length(cellregions))
    grid[BFaceNodes]=bfacenodes
    grid[BFaceRegions]=bfaceregions
    grid[BFaceGeometries]=VectorOfConstants(btype,length(bfaceregions))
    grid[CoordinateSystem]=csys
    return grid
end

"""
````
function simplexgrid(coord::Array{Tc,2},
                     cellnodes::Array{Ti,2},
                     cellregions::Array{Ti,1},
                     bfacenodes::Array{Ti,2},
                     bfaceregions::Array{Ti,1}
                     bedgenodes::Array{Ti,2},
                     bedgeregions::Array{Ti,1}
                     ) where {Tc,Ti}
````

    Create simplex grid from coordinates, cell-nodes-adjancency, cell-region-numbers,
    boundary-face-nodes adjacency, boundary-face-region-numbers, boundary-edge-nodes, and
    boundary-edge-region-numbers arrays.
"""
function simplexgrid(coord::Array{Tc,2},
                     cellnodes::Array{Ti,2},
                     cellregions::Array{Ti,1},
                     bfacenodes::Array{Ti,2},
                     bfaceregions::Array{Ti,1},
                     bedgenodes::Array{Ti,2},
                     bedgeregions::Array{Ti,1}
                     ) where {Tc,Ti}
    grid = simplexgrid(coord, cellnodes, cellregions, bfacenodes, bfaceregions)
    grid[BEdgeNodes]   = bedgenodes
    grid[BEdgeRegions] = bedgeregions
    return grid
end

simplexgrid(C,CN,CR,BFN,BFR)=simplexgrid(collect(C),collect(CN),collect(CR),collect(BFN), collect(BFR))

##########################################################
abstract type XCoordinates <: AbstractGridFloatArray1D end
abstract type YCoordinates <: AbstractGridFloatArray1D end
abstract type ZCoordinates <: AbstractGridFloatArray1D end

"""
$(TYPEDSIGNATURES)

Constructor for 1D grid.

Construct 1D grid from an array of node cordinates.
It creates two boundary regions with index 1 at the left end and
index 2 at the right end.

Primal grid holding unknowns: marked by `o`, dual
grid marking control volumes: marked by `|`.

```@raw html
 o-----o-----o-----o-----o-----o-----o-----o-----o
 |--|-----|-----|-----|-----|-----|-----|-----|--|
```
"""
function simplexgrid(_X::AbstractVector)
    X=collect_or_assign(_X)
    #    is_monotone(X) || error("X not monotone")
    coord=reshape(X,1,length(X))
    cellnodes=zeros(Int32,2,length(X)-1)
    cellregions=zeros(Int32,length(X)-1)
    for i=1:length(X)-1
        cellnodes[1,i]=i
        cellnodes[2,i]=i+1
        cellregions[i]=1
    end
    bfacenodes=Array{Int32}(undef,1,2)
    bfaceregions=zeros(Int32,2)
    bfacenodes[1,1]=1
    bfacenodes[1,2]=length(X)
    bfaceregions[1]=1
    bfaceregions[2]=2
    grid=simplexgrid(coord,
                     cellnodes,
                     cellregions,
                     bfacenodes,
                     bfaceregions)
    grid[XCoordinates]=X
    grid
end

##########################################################
"""
$(TYPEDSIGNATURES)

Constructor for 2D grid
from coordinate arrays.
Boundary region numbers count counterclockwise:

| location  |  number |
| --------- | ------- |
| south     |       1 |
| east      |       2 |
| north     |       3 |
| west      |       4 |

"""
function  simplexgrid(_X::AbstractVector, _Y::AbstractVector)
    X=collect_or_assign(_X)
    Y=collect_or_assign(_Y)
    is_monotone(X) || error("X not monotone")
    is_monotone(Y) || error("Y not monotone")

    Tc=promote_type(eltype(X),eltype(Y))

    function leq(x, x1, x2)
        if (x>x1)
            return false
        end
        if (x>x2)
            return false
        end
        return true
    end

    function geq(x, x1, x2)
        if (x<x1)
            return false
        end
        if (x<x2)
            return false
        end
        return true
    end

    nx=length(X)
    ny=length(Y)

    hmin=X[2]-X[1]
    for i=1:nx-1
        h=X[i+1]-X[i]
        if h <hmin
            hmin=h
        end
    end
    for i=1:ny-1
        h=Y[i+1]-Y[i]
        if h <hmin
            hmin=h
        end
    end

    @assert(hmin>0.0)
    eps=1.0e-5*hmin

    x1=X[1]+eps
    xn=X[nx]-eps
    y1=Y[1]+eps
    yn=Y[ny]-eps

    function  check_insert_bface(ibface,coord,n1,n2)
        if (geq(x1,coord[1,n1],coord[1,n2]))
            ibface=ibface+1
            bfacenodes[1,ibface]=n1
            bfacenodes[2,ibface]=n2
	    bfaceregions[ibface]=4
        elseif (leq(xn,coord[1,n1],coord[1,n2]))
            ibface=ibface+1
            bfacenodes[1,ibface]=n1
            bfacenodes[2,ibface]=n2
	    bfaceregions[ibface]=2
        elseif (geq(y1,coord[2,n1],coord[2,n2]))
            ibface=ibface+1
            bfacenodes[1,ibface]=n1
            bfacenodes[2,ibface]=n2
	    bfaceregions[ibface]=1
        elseif (leq(yn,coord[2,n1],coord[2,n2]))
            ibface=ibface+1
            bfacenodes[1,ibface]=n1
            bfacenodes[2,ibface]=n2
	    bfaceregions[ibface]=3
        end
        ibface
    end

    num_nodes=nx*ny
    num_cells=2*(nx-1)*(ny-1)
    num_bfacenodes=2*(nx-1)+2*(ny-1)

    coord=zeros(Tc,2,num_nodes)
    cellnodes=zeros(Int32,3,num_cells)
    cellregions=zeros(Int32,num_cells)
    bfacenodes=Array{Int32}(undef,2,num_bfacenodes)
    bfaceregions=zeros(Int32,num_bfacenodes)

    ipoint=0
    for iy=1:ny
        for ix=1:nx
            ipoint=ipoint+1
            coord[1,ipoint]=X[ix]
            coord[2,ipoint]=Y[iy]
        end
    end
    @assert(ipoint==num_nodes)

    icell=0
    for iy=1:ny-1
        for ix=1:nx-1
	    ip=ix+(iy-1)*nx
	    p00 = ip
	    p10 = ip+1
	    p01 = ip  +nx
	    p11 = ip+1+nx

            icell=icell+1
            cellnodes[1,icell]=p00
            cellnodes[2,icell]=p10
            cellnodes[3,icell]=p11
            cellregions[icell]=1

            icell=icell+1
            cellnodes[1,icell]=p11
            cellnodes[2,icell]=p01
            cellnodes[3,icell]=p00
            cellregions[icell]=1
        end
    end
    @assert(icell==num_cells)

    #lazy way to  create boundary grid

    ibface=0
    for icell=1:num_cells
        n1=cellnodes[1,icell]
	n2=cellnodes[2,icell]
	n3=cellnodes[3,icell]
        ibface=check_insert_bface(ibface,coord,n1,n2)
	ibface=check_insert_bface(ibface,coord,n1,n3)
	ibface=check_insert_bface(ibface,coord,n2,n3)
    end
    @assert(ibface==num_bfacenodes)

    grid=simplexgrid(coord,
                     cellnodes,
                     cellregions,
                     bfacenodes,
                     bfaceregions)

    grid[XCoordinates]=X
    grid[YCoordinates]=Y
    grid
end

##########################################################
"""
$(TYPEDSIGNATURES)

Constructor for 3D grid
from coordinate arrays.
Boundary region numbers:

| location  |  number |
| --------- | ------- |
| south     |       1 |
| east      |       2 |
| north     |       3 |
| west      |       4 |
| bottom    |       5 |
| top       |       6 |

"""
function  simplexgrid(_X::AbstractVector,_Y::AbstractVector,_Z::AbstractVector)
    X=collect_or_assign(_X)
    Y=collect_or_assign(_Y)
    Z=collect_or_assign(_Z)

    is_monotone(X) || error("X not monotone")
    is_monotone(Y) || error("Y not monotone")
    is_monotone(Z) || error("Z not monotone")

    Tc=promote_type(eltype(X),eltype(Y),eltype(Z))

    leq(x, x1, x2, x3)=x≤x1 && x≤x2 && x≤x3
    geq(x, x1, x2, x3)=x≥x1 && x≥x2 && x≥x3

    nx=length(X)
    ny=length(Y)
    nz=length(Z)

    hmin=X[2]-X[1]
    for i=1:nx-1
        h=X[i+1]-X[i]
        if h <hmin
            hmin=h
        end
    end
    for i=1:ny-1
        h=Y[i+1]-Y[i]
        if h <hmin
            hmin=h
        end
    end

    for i=1:nz-1
        h=Z[i+1]-Z[i]
        if h <hmin
            hmin=h
        end
    end

    @assert(hmin>0.0)
    eps=1.0e-5*hmin

    x1=X[1]+eps
    xn=X[end]-eps

    y1=Y[1]+eps
    yn=Y[end]-eps

    z1=Z[1]+eps
    zn=Z[end]-eps

    function  check_insert_bface(ibface,coord,bfacenodes,n1,n2,n3)
        if geq(x1,coord[1,n1],coord[1,n2],coord[1,n3])
            ibface=ibface+1
            bfacenodes[1,ibface]=n1
            bfacenodes[2,ibface]=n2
            bfacenodes[3,ibface]=n3
	    bfaceregions[ibface]=4
        elseif leq(xn,coord[1,n1],coord[1,n2],coord[1,n3])
            ibface=ibface+1
            bfacenodes[1,ibface]=n1
            bfacenodes[2,ibface]=n2
            bfacenodes[3,ibface]=n3
	    bfaceregions[ibface]=2
        elseif geq(y1,coord[2,n1],coord[2,n2],coord[2,n3])
            ibface=ibface+1
            bfacenodes[1,ibface]=n1
            bfacenodes[2,ibface]=n2
            bfacenodes[3,ibface]=n3
	    bfaceregions[ibface]=1
        elseif leq(yn,coord[2,n1],coord[2,n2],coord[2,n3])
            ibface=ibface+1
            bfacenodes[1,ibface]=n1
            bfacenodes[2,ibface]=n2
            bfacenodes[3,ibface]=n3
	    bfaceregions[ibface]=3
        elseif geq(z1,coord[3,n1],coord[3,n2],coord[3,n3])
            ibface=ibface+1
            bfacenodes[1,ibface]=n1
            bfacenodes[2,ibface]=n2
            bfacenodes[3,ibface]=n3
	    bfaceregions[ibface]=5
        elseif leq(zn,coord[3,n1],coord[3,n2],coord[3,n3])
            ibface=ibface+1
            bfacenodes[1,ibface]=n1
            bfacenodes[2,ibface]=n2
            bfacenodes[3,ibface]=n3
	    bfaceregions[ibface]=6
        end
        ibface
    end

    num_nodes=nx*ny*nz
    num_cells=6*(nx-1)*(ny-1)*(nz-1)
    num_bfacenodes=4*(nx-1)*(ny-1)+4*(nx-1)*(nz-1)+4*(ny-1)*(nz-1)

    Ti=Int32
    coord=zeros(Tc,3,num_nodes)
    cellnodes=zeros(Ti,4,num_cells)
    cellregions=zeros(Ti,num_cells)
    bfacenodes=zeros(Ti,3,num_bfacenodes)
    bfaceregions=zeros(Ti,num_bfacenodes)

    ipoint=0
    for iz=1:nz
        for iy=1:ny
            for ix=1:nx
                ipoint=ipoint+1
                coord[1,ipoint]=X[ix]
                coord[2,ipoint]=Y[iy]
                coord[3,ipoint]=Z[iz]
            end
        end
    end

    @assert(ipoint==num_nodes)

    icell=0
    nxy=nx*ny
    for iz=1:nz-1
        for iy=1:ny-1
            for ix=1:nx-1

	        ip=ix+(iy-1)*nx+(iz-1)*nxy;

	        p000 = ip ;
	        p100 = ip+1;
	        p010 = ip  +nx;
	        p110 = ip+1+nx;
	        p001 = ip     +nxy;
	        p101 = ip+1   +nxy;
	        p011 = ip  +nx+nxy;
	        p111 = ip+1+nx+nxy;

                icell=icell+1; cellregions[icell]=1; @. cellnodes[:,icell]=(p000,p100,p110,p111)
                icell=icell+1; cellregions[icell]=1; @. cellnodes[:,icell]=(p000,p100,p101,p111)
                icell=icell+1; cellregions[icell]=1; @. cellnodes[:,icell]=(p000,p010,p011,p111)
                icell=icell+1; cellregions[icell]=1; @. cellnodes[:,icell]=(p000,p010,p110,p111)
                icell=icell+1; cellregions[icell]=1; @. cellnodes[:,icell]=(p000,p001,p101,p111)
                icell=icell+1; cellregions[icell]=1; @. cellnodes[:,icell]=(p000,p001,p011,p111)
            end
        end
    end
    @assert(icell==num_cells)
    #lazy way to  create boundary grid
    # lazy but easy... bad for partitioning !!!

    ibface=0
    for icell=1:num_cells
        n1=cellnodes[1,icell]
	n2=cellnodes[2,icell]
	n3=cellnodes[3,icell]
	n4=cellnodes[4,icell]
        ibface=check_insert_bface(ibface,coord,bfacenodes,n1,n2,n3)
	ibface=check_insert_bface(ibface,coord,bfacenodes,n1,n2,n4)
	ibface=check_insert_bface(ibface,coord,bfacenodes,n1,n3,n4)
	ibface=check_insert_bface(ibface,coord,bfacenodes,n2,n3,n4)
    end
    @assert(ibface==num_bfacenodes)

    grid=simplexgrid(coord,
                     cellnodes,
                     cellregions,
                     bfacenodes,
                     bfaceregions)

    grid[XCoordinates]=X
    grid[YCoordinates]=Y
    grid[ZCoordinates]=Z
    grid
end

"""
$(TYPEDSIGNATURES)
Create tensor product of 2D  grid and 1D coordinate arrays.

Cellregions and outer facet regions are taken over from 2D grid.
Top an bottom facet region given in ibot and itop, by default these
are detected from the largest cell region number of the 2D grid.
"""
function simplexgrid(grid2::ExtendableGrid, coordZ; ibot=0, itop=0)

    coord2=grid2[Coordinates]
    nnodes2=size(coord2,2)
    nnodesZ=length(coordZ)
    nnodes3=nnodes2*nnodesZ
    coord3=zeros(3,nnodes3)

    cells2=grid2[CellNodes]
    ncells2=size(cells2,2)
    ncells3=ncells2*3*(nnodesZ-1)
    cells3=zeros(Int, 4,ncells3)

    cellregions2=grid2[CellRegions]
    cellregions3=zeros(Int, ncells3)

    bfaces2=grid2[BFaceNodes]
    nbfaces2=size(bfaces2,2)
    nbfaces3=2*ncells2+2*nbfaces2*(nnodesZ-1)
    bfaces3=zeros(Int,3,nbfaces3)

    bfaceregions2=grid2[BFaceRegions]
    bfaceregions3=zeros(Int,nbfaces3)

    # handle defaults of top and bottom bc numbers
    nbfr=maximum(bfaceregions2)
    if ibot==0
        nbfr=nbfr+1
        ibot=nbfr
    end

    if itop==0
        nbfr=nbfr+1
        itop=nbfr
    end

    # 3D coordinates
    i3=1
    for iZ=1:nnodesZ
        for i2=1:nnodes2
            coord3[1,i3]=coord2[1,i2]
            coord3[2,i3]=coord2[2,i2]
            coord3[3,i3]=coordZ[iZ]
            i3+=1
        end
    end

    # 3D cells
    i3=1
    ishift=0
    for iZ=1:nnodesZ-1
        for i2=1:ncells2
            n1 = cells2[1,i2]+ishift
            n2 = cells2[2,i2]+ishift
            n3 = cells2[3,i2]+ishift

            n1>n2 ? (n1,n2) = (n2,n1) : nothing
            n2>n3 ? (n2,n3) = (n3,n2) : nothing
            n1>n2 ? (n1,n2) = (n2,n1) : nothing

            cells3[:,i3  ].=(n1,n2,        n3,        n3+nnodes2)
            cells3[:,i3+1].=(n1,n2,        n2+nnodes2,n3+nnodes2)
            cells3[:,i3+2].=(n1,n1+nnodes2,n2+nnodes2,n3+nnodes2)

            cellregions3[i3]   = cellregions2[i2]
            cellregions3[i3+1] = cellregions2[i2]
            cellregions3[i3+2] = cellregions2[i2]

            i3+=3
        end
        ishift+=nnodes2
    end

    #xy parallel  facets
    i3=1
    for iZ ∈ [1,nnodesZ]
        ishift=(iZ-1)*nnodes2
        for i2=1:ncells2
            @views bfaces3[:,i3].=cells2[:,i2].+ishift
            bfaceregions3[i3] = iZ==1 ? ibot : itop
            i3+=1
        end
    end

    # facets vertical to the x-y plane
    ishift=0
    for iZ=1:nnodesZ-1
        for i2=1:nbfaces2
            n1=bfaces2[1,i2]+ishift
            n2=bfaces2[2,i2]+ishift
            n1>n2 ? (n1,n2)=(n2,n1) : nothing
            @views bfaces3[:,i3  ].=(n1,n2,        n2+nnodes2)
            @views bfaces3[:,i3+1].=(n1,n1+nnodes2,n2+nnodes2)
            bfaceregions3[i3]   = bfaceregions2[i2]
            bfaceregions3[i3+1] = bfaceregions2[i2]
            i3+=2
        end
        ishift+=nnodes2
    end
    simplexgrid(coord3,cells3,cellregions3,bfaces3,bfaceregions3)
end

######################################################
"""
$(TYPEDSIGNATURES)

Read grid from file. Currently for pdelib sg format only
"""
function simplexgrid(file::String;format="")
    (fbase,fext)=splitext(file)
    if format==""
        format=fext[2:end]
    end
    @assert format=="sg"

    tks=TokenStream(file)
    expecttoken(tks,"SimplexGrid")
    version=parse(Float64,gettoken(tks))
    version20=false;

    if (version==2.0)
        version20=true;
    elseif (version==2.1)
        version20=false;
    else
        error("Read grid: wrong format version: $(version)")
    end

    dim::Int32=0
    coord=Array{Float64,2}(undef,0,0)
    cells=Array{Int32,2}(undef,0,0)
    regions=Array{Int32,1}(undef,0)
    faces=Array{Int32,2}(undef,0,0)
    bregions=Array{Int32,1}(undef,0)
    while(true)
        if (trytoken(tks,"DIMENSION"))
            dim=parse(Int32,gettoken(tks));
        elseif (trytoken(tks,"NODES"))
            nnodes=parse(Int32,gettoken(tks));
            embdim=parse(Int32,gettoken(tks));
            if(dim!=embdim)
                error("Dimension error (DIMENSION $(dim)) in section NODES")
            end
            coord=Array{Float64,2}(undef,dim,nnodes)
            for inode=1:nnodes
                for idim=1:embdim
                    coord[idim,inode]=parse(Float64,gettoken(tks))
                end
            end
        elseif (trytoken(tks,"CELLS"))
            ncells=parse(Int32,gettoken(tks));
            cells=Array{Int32,2}(undef,dim+1,ncells)
            regions=Array{Int32,1}(undef,ncells)
            for icell=1:ncells
                for inode=1:dim+1
                    cells[inode,icell]=parse(Int32,gettoken(tks));
                end
                regions[icell]=parse(Int32,gettoken(tks));
	        if version20
		    for j=1:dim+1
		        gettoken(tks);  # skip file format garbage
                    end
                end
            end
        elseif (trytoken(tks,"FACES"))
            nfaces=parse(Int32,gettoken(tks));
            faces=Array{Int32,2}(undef,dim,nfaces)
            bregions=Array{Int32,1}(undef,nfaces)
            for iface=1:nfaces
                for inode=1:dim
                    faces[inode,iface]=parse(Int32,gettoken(tks));
                end
                bregions[iface]=parse(Int32,gettoken(tks));
	        if (version20)
		    for j=1:dim+2
		        gettoken(tks); #skip file format garbage
                    end
                end
            end
        else
            expecttoken(tks,"END")
            break
        end
    end
    simplexgrid(coord,cells,regions,faces,bregions);
end

"""
$(TYPEDSIGNATURES)

Write grid to file. Currently for pdelib sg format only
"""
function Base.write(fname::String, g::ExtendableGrid; format="")
    (fbase,fext)=splitext(fname)
    if format==""
        format=fext[2:end]
    end
    @assert format=="sg"

    dim_g=dim_grid(g)
    dim_s=dim_space(g)
    nn=num_nodes(g)
    nc=num_cells(g)
    nbf=num_bfaces(g)
    coord=g[Coordinates]
    cellnodes=g[CellNodes]
    bfacenodes=g[BFaceNodes]
    cellregions=g[CellRegions]
    bfaceregions=g[BFaceRegions]

    # TODO: replace @sprintf by someting non-allocating
    open(fname, "w") do file
        write(file,@sprintf("SimplexGrid"))
        write(file,@sprintf(" "))
        write(file,@sprintf("2.1\n"))
        write(file,@sprintf("#created by ExtendableGrids.jl (c) J.Fuhrmann et al\n"))
        write(file,@sprintf("#mailto:{fuhrmann|streckenbach}@wias-berlin.de\n"))
        write(file,@sprintf("#%s\n",Dates.format(Dates.now(),"yyyy-mm-ddTHH-mm-SS")))

        write(file,@sprintf("DIMENSION\n%d\n",dim_g))
        write(file,@sprintf("NODES\n%d %d\n",nn,dim_s))

        for inode=1:nn
            for idim=1:dim_s
	        write(file,@sprintf("%.20e ",coord[idim,inode]))
                write(file,@sprintf("\n"))
            end
        end

        write(file,@sprintf("CELLS\n%d\n",nc))
        for icell=1:nc
            for inode=1:dim_g+1
	        write(file,@sprintf("%d ", cellnodes[inode,icell]))
            end
            write(file,@sprintf("%d\n",cellregions[icell]))
        end

        write(file,@sprintf("FACES\n%d\n",nbf))
        for ibface=1:nbf
            for inode=1:dim_g
	        write(file,@sprintf("%d ", bfacenodes[inode,ibface]))
            end
            write(file,@sprintf("%d\n",bfaceregions[ibface]))
        end
        write(file,@sprintf("END\n"))
        flush(file)
        flush(file)
    end
    nothing
end

"""
$(TYPEDSIGNATURES)

Merge two grids along their common boundary facets.

- g1: First grid to be merged
- g2: Second grid to be merged
- breg:  Interior boundary region number of merged facets. If zero (default value), no extra facets are generated.
- tol:  Distance below which two points are seen as identical. Default: 1.0e-10

"""
function glue(g1,g2; breg=0, tol=1.0e-10)

    dim=dim_space(g1)

    bfn1=g1[BFaceNodes]
    bfn2=g2[BFaceNodes]

    nbf1=size(bfn1,2)
    nbf2=size(bfn2,2)

    coord1=g1[Coordinates]
    coord2=g2[Coordinates]

    nn1=size(coord1,2)
    nn2=size(coord2,2)

    # numbers of faces in grid1 which match nodes in grid2
    matching_faces=zeros(Int,nbf2)
    n_matching_faces=0

    # numbers of nodes in grid1 which match nodes in grid2
    matching_nodes=zeros(Int,nn2)
    n_matching_nodes=0

    # Add matching pair to index of matching pairs
    function add_match!(match_list, nmatch, i1,i2)
        @assert match_list[i2]==0 || match_list[i2]==i1
        if match_list[i2]==0
            match_list[i2]=i1;
            nmatch+=1
        end
        nmatch
    end

    # Add two point indices to list of matching points
    add_matching_points(ip1,ip2) = n_matching_nodes = add_match!(matching_nodes,n_matching_nodes,ip1,ip2)

    # Add two face indices to list of matching faces
    add_matching_faces(if1,if2) = n_matching_faces = add_match!(matching_faces,n_matching_faces,if1,if2)

    ## Check if two points match (distance < tol)
    points_match(ip1, ip2) = @views norm(coord1[:,ip1]-coord2[:,ip2])<tol

    # Check if points in faces match, if so, add them to point matching index
    # Check if faces match, if so, add them to face matching index
    function match_faces(if1,if2)
        nmatch=0
        for ip1=1:dim
            for ip2=1:dim
                xip1=bfn1[ip1,if1]
                xip2=bfn2[ip2,if2]
                if points_match(xip1,xip2)
                    add_matching_points(xip1,xip2)
                    nmatch+=1
                end
            end
        end

      if nmatch==dim
          add_matching_faces(if1,if2);
      end
    end

    # Run over all pairs of boundary faces and try to match them
    for ibf1=1:nbf1
        for ibf2=1:nbf2
            match_faces(ibf1,ibf2)
        end
    end
    @info "glue: matches found: nodes $(n_matching_nodes) bfaces: $(n_matching_faces)"

    creg1=g1[CellRegions]
    creg2=g2[CellRegions]

    breg1=g1[BFaceRegions]
    breg2=g2[BFaceRegions]

    cn1=g1[CellNodes]
    cn2=g2[CellNodes]

    nc1=size(cn1,2)
    nc2=size(cn2,2)

    # transposed list of matching faces
    mf1=zeros(Int,nbf1)
    for if2=1:nbf2
        if matching_faces[if2]!=0
            mf1[matching_faces[if2]]=if2
        end
    end

    mfac = breg == 0 ? 2 : 1

    coordx=zeros(dim,nn1+nn2-n_matching_nodes)
    cnx=zeros(Int, dim+1,nc1+nc2)
    cregx=zeros(Int,nc1+nc2)
    bfnx=zeros(Int,dim,nbf1+nbf2-mfac*n_matching_faces)
    bregx=zeros(Int,nbf1+nbf2-mfac*n_matching_faces)

    #   copy all data from g1 into new fields
    for ix=1:nn1
        @views coordx[:,ix].=coord1[:,ix]
    end

    for ix=1:nc1
        cregx[ix]=creg1[ix]
        @views cnx[:,ix].=cn1[:,ix]
    end

    ibfx=1
    for ibf1=1:nbf1
        if mf1[ibf1]!=0 && breg == 0
            continue
        end
        if mf1[ibf1]!=0
            bregx[ibfx]=breg
        else
            bregx[ibfx]=breg1[ibf1]
        end
        @views bfnx[:,ibfx].=bfn1[:,ibf1]
        ibfx+=1
    end

    # re-calculate matching nodes with new global numbers
    # add missing coordinates
    ix=nn1+1
    for in2=1:nn2
        if matching_nodes[in2]!=0
            continue
        end
        matching_nodes[in2]=ix
        @views coordx[:,ix].=coord2[:,in2]
        ix=ix+1
    end

    # copy missing data from g2 into arrays.
    ix=nc1+1
    for ic2=1:nc2
        cregx[ix]=creg2[ic2]
        for id=1:dim+1
            cnx[id,ix]=matching_nodes[cn2[id,ic2]]
        end
        ix=ix+1
    end

    for ibf2=1:nbf2
        if matching_faces[ibf2]!=0
            continue
        end
        bregx[ibfx]=breg2[ibf2]
        for id=1:dim
            bfnx[id,ibfx]=matching_nodes[bfn2[id,ibf2]]
        end
        ibfx+=1
    end
    @assert ibfx == nbf1+nbf2-mfac*n_matching_faces+1
    simplexgrid(coordx,cnx,cregx,bfnx,bregx)
end


text:
# This code and all turorial code can be found on https://github.com/forio/julia-tutorials

data = readcsv("gasoline.csv")

# Get the columns that correspond with our X values
# and create a new matrix to hold them
# syntax: data[row, column]
# to get a whole row do data[row, :]
# to get a whole column do data[:, column]
# to get a range of columns do data[:, column_a:column_b]

x = data[:, 2:4]
y = data[:, 6]

# Call linreg
coefs = linreg(x, y)


text:
function jump_mimo(solver, seed, n; verbose = false, test = false)

    # n = 3
    m = 10n
    s, H, y, L = mimo_data(seed, m, n)

    nvars = sympackedlen(n + 1)

    model = Model(with_optimizer(solver))
    @variable(model, X[1:n+1, 1:n+1], PSD)
    for j in 1:(n+1), i in j:(n+1)
        @constraint(model, X[i, j] <=  1.0)
        @constraint(model, X[i, j] >= -1.0)
    end
    @objective(model, Min, sum(L[i, j] * X[i, j] for j in 1:n+1, i in 1:n+1))
    @constraint(model, ctr[i in 1:n+1], X[i, i] == 1.0)

    teste = @time optimize!(model)

    XX = value.(X)

    if test
        for j in 1:n+1, i in 1:n+1
            @test 1.01 > abs(XX[i,j]) > 0.99
        end
    end

    verbose && mimo_eval(s,H,y,L,XX)

    objval = objective_value(model)
    stime = MOI.get(model, MOI.SolveTime())

    # @show tp = typeof(model.moi_backend.optimizer.model.optimizer)
    # @show fieldnames(tp)
    rank = -1
    try
        @show rank = model.moi_backend.optimizer.model.optimizer.sol.final_rank
    catch
    end
    status = 0
    if JuMP.termination_status(model) == MOI.OPTIMAL
        status = 1
    end
    return (objval, stime, rank, status)
    # return (objval, stime)
end


text:
using MPI, Test

include("../testhelpers.jl")

@testset "Ocean" begin
    runmpi(joinpath(@__DIR__, "HydrostaticBoussinesq/test_ocean_gyre_short.jl"))
    runmpi(joinpath(@__DIR__, "SplitExplicit/test_spindown_short.jl"))
    include(joinpath("OceanProblems", "test_initial_value_problem.jl"))
    include(joinpath(
        "HydrostaticBoussinesqModel",
        "test_hydrostatic_boussinesq_model.jl",
    ))
end


text:
# ------------------------------------------------------------------------
# categorical

MMI.categorical(::FI, a...; kw...) = categorical(a...; kw...)

# ------------------------------------------------------------------------
# matrix

MMI.matrix(::FI, ::Val{:table}, X; kw...) = Tables.matrix(X; kw...)

# ------------------------------------------------------------------------
# int

MMI.int(::FI, x) = throw(
    DomainError(x, "Can only convert categorical elements to integers. "))

MMI.int(::FI, x::Missing)       = missing
MMI.int(::FI, x::AbstractArray) = int.(x)

MMI.int(::FI, x::CategoricalElement) = CategoricalArrays.order(x.pool)[x.level]

# ------------------------------------------------------------------------
# classes

MMI.classes(::FI, p::CategoricalPool) =
    [p[i] for i in invperm(CategoricalArrays.order(p))]

# older method that avoids inverting a permutation but has dict lookup:
# classes(p::CategoricalPool) = [p.valindex[p.invindex[v]] for v in p.levels]
MMI.classes(::FI, x::CategoricalElement) = classes(x.pool)

# ------------------------------------------------------------------------
# decoder

struct CategoricalDecoder{T,R}
    pool::CategoricalPool{T,R}
    invorder::Vector{Int}
end

MMI.decoder(::FI, x::CategoricalElement) =
    CategoricalDecoder(x.pool, sortperm(x.pool.order))

(d::CategoricalDecoder{T,R})(i::Integer) where {T,R} =
    CategoricalValue{T,R}(d.invorder[i], d.pool)
(d::CategoricalDecoder{String,R})(i::Integer) where R =
    CategoricalString{R}(d.invorder[i], d.pool)
(d::CategoricalDecoder)(a::AbstractArray{<:Integer}) = d.(a)

# ------------------------------------------------------------------------
# table

function MMI.table(::FI, cols::NamedTuple; prototype=NamedTuple())
    Tables.istable(prototype) || error("`prototype` is not a table. ")
    if !Tables.istable(cols)
        tuple_of_vectors = tuple((collect(v) for v in values(cols))...)
        names = keys(cols)
        cols = NamedTuple{names}(tuple_of_vectors)
    end
    return Tables.materializer(prototype)(cols)
end

function MMI.table(::FI, A::AbstractMatrix; names=nothing, prototype=nothing)
    if names === nothing
        _names = [Symbol(:x, j) for j in 1:size(A, 2)]
    else
        _names = collect(names)
    end
    matrix_table = Tables.table(A, header=_names)
    prototype === nothing && return matrix_table
    return Tables.materializer(prototype)(matrix_table)
end

# ------------------------------------------------------------------------
# nrows, selectrows, selectcols

function MMI.nrows(::FI, ::Val{:table}, X)
    Tables.columnaccess(X) || return length(collect(X))
    # if has columnaccess
    cols = Tables.columntable(X)
    !isempty(cols) || return 0
    return length(cols[1])
end

MMI.selectrows(::FI, ::Val{:table}, X, ::Colon) = X
MMI.selectcols(::FI, ::Val{:table}, X, ::Colon) = X

function MMI.selectrows(::FI, ::Val{:table}, X, r)
    r = r isa Integer ? (r:r) : r
    # next uncommented line is a hack; see
    # https://github.com/alan-turing-institute/MLJBase.jl/issues/151
    isdataframe(X) && return X[r, :]
    cols = Tables.columntable(X)
    new_cols = NamedTuple{keys(cols)}(tuple((c[r] for c in values(cols))...))
    return Tables.materializer(X)(new_cols)
end

function MMI.selectcols(::FI, ::Val{:table}, X, c::Union{Symbol,Integer})
    cols = Tables.columntable(X) # named tuple of vectors
    return cols[c]
end

function MMI.selectcols(::FI, ::Val{:table}, X, c::AbstractArray)
    cols = Tables.columntable(X) # named tuple of vectors
    newcols = project(cols, c)
    return Tables.materializer(X)(newcols)
end

# -------------------------------
# utils for select*

# project named tuple onto a tuple with only specified `labels` or indices:
project(t::NamedTuple, labels::AbstractArray{Symbol}) =
    NamedTuple{tuple(labels...)}(t)
project(t::NamedTuple, label::Colon) = t
project(t::NamedTuple, label::Symbol) = project(t, [label,])
project(t::NamedTuple, indices::AbstractArray{<:Integer}) =
    NamedTuple{tuple(keys(t)[indices]...)}(tuple([t[i] for i in indices]...))
project(t::NamedTuple, i::Integer) = project(t, [i,])

# utils for selectrows
typename(X)    = split(string(supertype(typeof(X)).name), '.')[end]
isdataframe(X) = typename(X) == "AbstractDataFrame"


text:
# Test script to test solver for an lp

using COSMO, LinearAlgebra, SparseArrays, Test

# Linear program example
# min c'x
# s.t. Ax <= b
#      x >= 1,  x2 >= 5, x1+x3 >= 4
# with data matrices
c = [1; 2; 3; 4.]
A = Matrix(1.0I, 4, 4)
b = [10; 10; 10; 10]
n = 4
# -------------------
# create constraints A * x + b in set
# -------------------
# Ax <= b
c1 = COSMO.Constraint(-A, b, COSMO.Nonnegatives)
# x >= 1
c2 = COSMO.Constraint(Matrix(1.0I, n, n), -ones(n), COSMO.Nonnegatives)
# x2 >= 5
c3 = COSMO.Constraint(1, -5, COSMO.Nonnegatives, n, 2:2)
# x1 + x3 >= 4
c4 = COSMO.Constraint([1 0 1 0], -4, COSMO.Nonnegatives)

# -------------------
# define cost function
# -------------------
P = spzeros(4, 4)
q = c

# -------------------
# assemble solver model
# -------------------
settings = COSMO.Settings(max_iter=2500, verbose=true, eps_abs = 1e-4, eps_rel = 1e-5)
model = COSMO.Model()
assemble!(model, P, q, [c1; c2; c3; c4], settings = settings)
res = COSMO.optimize!(model);

@testset "Linear Problem" begin
  @test isapprox(res.x[1:4], [3; 5; 1; 1], atol=1e-2, norm = (x -> norm(x, Inf)))
  @test isapprox(res.obj_val, 20.0, atol=1e-2)
end
nothing


text:
@testset "Double64(x)" begin
  @test Hi(Double64(1.0)) === 1.0
  @test Lo(Double64(1.0)) === 0.0
  @test Hio(Double64(1.0)) === (1.0, 0.0)
  @test Double64(1.0e7, 1.0e-7) === Double64(1.00000000000001e7, -5.828380584716842e-10)
  @test HiLo(Double64(1.0e7, 1.0e-7)) === (1.00000000000001e7, -5.828380584716842e-10)
  @test HiLo(Double64((1.0e7, 1.0e-7))) === (1.0e7, 1.0e-7)
  @test Double64(BigFloat(pi)) === Double64((3.141592653589793, 1.2246467991473532e-16))
end

@testset "ComplexD64(x)" begin
end

@testset "T(::Double64)" begin
end

@testset "T(::ComplexD64)" begin
end


text:
sudoku_start = [
    9  3  0  0  0  0  0  4  0
    0  0  0  0  4  2  0  9  0
    8  0  0  1  9  6  7  0  0
    0  0  0  4  7  0  0  0  0
    0  2  0  0  0  0  0  6  0
    0  0  0  0  2  3  0  0  0
    0  0  8  5  3  1  0  0  2
    0  9  0  2  8  0  0  0  0
    0  7  0  0  0  0  0  5  3
]

qap_weights = [
    0  1  2  2  3  4  4  5  3  5  6  7
    1  0  1  1  2  3  3  4  2  4  5  6
    2  1  0  2  1  2  2  3  1  3  4  5
    2  1  2  0  1  2  2  3  3  3  4  5
    3  2  1  1  0  1  1  2  2  2  3  4
    4  3  2  2  1  0  2  3  3  1  2  3
    4  3  2  2  1  2  0  1  3  1  2  3
    5  4  3  3  2  3  1  0  4  2  1  2
    3  2  1  3  2  3  3  4  0  4  5  6
    5  4  3  3  2  1  1  2  4  0  1  2
    6  5  4  4  3  2  2  1  5  1  0  1
    7  6  5  5  4  3  3  2  6  2  1  0
]

qap_distances = [
    0  3  4  6  8  5  6  6  5  1  4  6
    3  0  6  3  7  9  9  2  2  7  4  7
    4  6  0  2  6  4  4  4  2  6  3  6
    6  3  2  0  5  5  3  3  9  4  3  6
    8  7  6  5  0  4  3  4  5  7  6  7
    5  9  4  5  4  0  8  5  5  5  7  5
    6  9  4  3  3  8  0  6  8  4  6  7
    6  2  4  3  4  5  6  0  1  5  5  3
    5  2  2  9  5  5  8  1  0  4  5  2
    1  7  6  4  7  5  4  5  4  0  7  7
    4  4  3  3  6  7  6  5  5  7  0  9
    6  7  6  6  7  5  7  3  2  7  9  0
]

atoms_compounds = [
    1 0 0
    2 0 0
    2 0 1
    0 1 0
    0 2 0
    1 1 0
    0 1 1
    0 0 1
    0 0 2
    1 0 1
]

elements_weights = [2 1 1]

standard_free_energy = [
    -6.0890
    -17.164
    -34.054
    -5.9140
    -24.721
    -14.986
    -24.100
    -10.708
    -26.662
    -22.179
]

compounds_names = "x₁⋅H + x₂⋅H₂ + x₃⋅H₂O + x₄⋅N + x₅⋅N₂ + x₆⋅NH + x₇⋅NO + x₈⋅O + x₉⋅O₂ + x₁₀⋅OH"

mixture_name =  "½⋅N₂H₄ + ½⋅O₂"

equation = compounds_names * " = " * mixture_name

processing_times = [3, 2, 1]
due_times = [5, 6, 8]


text:
function claytonCopula(θ::Int64, samples::Int64, dim::Int64)
	res = zeros(dim, samples)
	U_dist = Uniform(0,1)
	γ_dist = Gamma(1/θ,1)
	for i = 1:samples
		V = rand(U_dist, dim)
		γ = rand(γ_dist, 1)
		Ψ(t) = (t+1)^(-1/θ)
		U_γ = Ψ.(-log.(V)./γ)
        for j = 1:dim
            res[j,i] = U_γ[j]
        end
	end
    return res
end


text:
abstract Optimizer
immutable Options{TCallback <: Union{Void, Function}}
    x_tol::Float64
    f_tol::Float64
    g_tol::Float64
    iterations::Int
    store_trace::Bool
    show_trace::Bool
    extended_trace::Bool
    autodiff::Bool
    show_every::Int
    callback::TCallback
    time_limit::Float64
end

function Options(;
        x_tol::Real = 1e-32,
        f_tol::Real = 1e-32,
        g_tol::Real = 1e-8,
        iterations::Integer = 1_000,
        store_trace::Bool = false,
        show_trace::Bool = false,
        extended_trace::Bool = false,
        autodiff::Bool = false,
        show_every::Integer = 1,
        callback = nothing,
        time_limit = NaN)
    show_every = show_every > 0 ? show_every: 1
    if extended_trace && callback == nothing
        show_trace = true
    end
    Options{typeof(callback)}(
        Float64(x_tol), Float64(f_tol), Float64(g_tol), Int(iterations),
        store_trace, show_trace, extended_trace, autodiff, Int(show_every),
        callback, time_limit)
end

function print_header(options::Options)
    if options.show_trace
        @printf "Iter     Function value   Gradient norm \n"
    end
end

function print_header(method::Optimizer)
        @printf "Iter     Function value   Gradient norm \n"
end

immutable OptimizationState{T <: Optimizer}
    iteration::Int
    value::Float64
    g_norm::Float64
    metadata::Dict
end

typealias OptimizationTrace{T} Vector{OptimizationState{T}}

abstract OptimizationResults

type MultivariateOptimizationResults{T,N,M} <: OptimizationResults
    method::String
    initial_x::Array{T,N}
    minimizer::Array{T,N}
    minimum::Float64
    iterations::Int
    iteration_converged::Bool
    x_converged::Bool
    x_tol::Float64
    f_converged::Bool
    f_tol::Float64
    g_converged::Bool
    g_tol::Float64
    trace::OptimizationTrace{M}
    f_calls::Int
    g_calls::Int
    h_calls::Int
end

type UnivariateOptimizationResults{T,M} <: OptimizationResults
    method::String
    initial_lower::T
    initial_upper::T
    minimizer::T
    minimum::Float64
    iterations::Int
    iteration_converged::Bool
    converged::Bool
    rel_tol::Float64
    abs_tol::Float64
    trace::OptimizationTrace{M}
    f_calls::Int
end

immutable NonDifferentiableFunction
    f::Function
end

immutable DifferentiableFunction
    f::Function
    g!::Function
    fg!::Function
end

immutable TwiceDifferentiableFunction
    f::Function
    g!::Function
    fg!::Function
    h!::Function
end

function Base.show(io::IO, t::OptimizationState)
    @printf io "%6d   %14e   %14e\n" t.iteration t.value t.g_norm
    if !isempty(t.metadata)
        for (key, value) in t.metadata
            @printf io " * %s: %s\n" key value
        end
    end
    return
end

function Base.show(io::IO, tr::OptimizationTrace)
    @printf io "Iter     Function value   Gradient norm \n"
    @printf io "------   --------------   --------------\n"
    for state in tr
        show(io, state)
    end
    return
end

function Base.show(io::IO, r::MultivariateOptimizationResults)
    @printf io "Results of Optimization Algorithm\n"
    @printf io " * Algorithm: %s\n" method(r)
    if length(join(initial_state(r), ",")) < 40
        @printf io " * Starting Point: [%s]\n" join(initial_state(r), ",")
    else
        @printf io " * Starting Point: [%s, ...]\n" join(initial_state(r)[1:2], ",")
    end
    if length(join(minimizer(r), ",")) < 40
        @printf io " * Minimizer: [%s]\n" join(minimizer(r), ",")
    else
        @printf io " * Minimizer: [%s, ...]\n" join(minimizer(r)[1:2], ",")
    end
    @printf io " * Minimum: %e\n" minimum(r)
    @printf io " * Iterations: %d\n" iterations(r)
    @printf io " * Convergence: %s\n" converged(r)
    if r.method == "Nelder-Mead"
        @printf io "   *  √(Σ(yᵢ-ȳ)²)/n < %.1e: %s\n" g_tol(r) g_converged(r)
    else
        @printf io "   * |x - x'| < %.1e: %s\n" x_tol(r) x_converged(r)
        @printf io "   * |f(x) - f(x')| / |f(x)| < %.1e: %s\n" f_tol(r) f_converged(r)
        @printf io "   * |g(x)| < %.1e: %s\n" g_tol(r) g_converged(r)
    end
    @printf io "   * Reached Maximum Number of Iterations: %s\n" iteration_limit_reached(r)
    @printf io " * Objective Function Calls: %d" f_calls(r)
    if !(r.method in ("Nelder-Mead", "Simulated Annealing"))
        @printf io "\n * Gradient Calls: %d" g_calls(r)
    end
    return
end

function Base.show(io::IO, r::UnivariateOptimizationResults)
    @printf io "Results of Optimization Algorithm\n"
    @printf io " * Algorithm: %s\n" method(r)
    @printf io " * Search Interval: [%f, %f]\n" lower_bound(r) upper_bound(r)
    @printf io " * Minimizer: %e\n" minimizer(r)
    @printf io " * Minimum: %e\n" minimum(r)
    @printf io " * Iterations: %d\n" iterations(r)
    @printf io " * Convergence: max(|x - x_upper|, |x - x_lower|) <= 2*(%.1e*|x|+%.1e): %s\n" rel_tol(r) abs_tol(r) converged(r)
    @printf io " * Objective Function Calls: %d" f_calls(r)
    return
end

function Base.append!(a::MultivariateOptimizationResults, b::MultivariateOptimizationResults)
    a.iterations += iterations(b)
    a.minimizer = minimizer(b)
    a.minimum = minimum(b)
    a.iteration_converged = iteration_limit_reached(b)
    a.x_converged = x_converged(b)
    a.f_converged = f_converged(b)
    a.g_converged = g_converged(b)
    append!(a.trace, b.trace)
    a.f_calls += f_calls(b)
    a.g_calls += g_calls(b)
end

# TODO: Expose ability to do forward and backward differencing
function DifferentiableFunction(f::Function)
    function g!(x::Array, storage::Array)
        Calculus.finite_difference!(f, x, storage, :central)
        return
    end
    function fg!(x::Array, storage::Array)
        g!(x, storage)
        return f(x)
    end
    return DifferentiableFunction(f, g!, fg!)
end

function DifferentiableFunction(f::Function, g!::Function)
    function fg!(x::Array, storage::Array)
        g!(x, storage)
        return f(x)
    end
    return DifferentiableFunction(f, g!, fg!)
end

function TwiceDifferentiableFunction(f::Function)
    function g!(x::Vector, storage::Vector)
        Calculus.finite_difference!(f, x, storage, :central)
        return
    end
    function fg!(x::Vector, storage::Vector)
        g!(x, storage)
        return f(x)
    end
    function h!(x::Vector, storage::Matrix)
        Calculus.finite_difference_hessian!(f, x, storage)
        return
    end
    return TwiceDifferentiableFunction(f, g!, fg!, h!)
end

function TwiceDifferentiableFunction(f::Function,
                                     g!::Function,
                                     h!::Function)
    function fg!(x::Vector, storage::Vector)
        g!(x, storage)
        return f(x)
    end
    return TwiceDifferentiableFunction(f, g!, fg!, h!)
end


text:
# This file has an abstract interface for a Tree structure
# and some traversal iterators.
# I wrote it for another project, but use it for the Finite-Difference Taylor models here.

module Trees

abstract type Node end

# methods that can/should be customized
children(n :: N ) where{N<:Node}= isdefined( n, :children ) ? n.children : N[];
parent( n :: Node ) = isdefined( n, :parent ) ? n.parent : nothing;

@doc """
	depth(node)

Return depth of `node` in a tree.
Root node is assigned depth 0, its children have depth 1 and so forth.

This default implementation goas up the tree starting at `node` as long
as there is a parent node. \n
It can be overwritten for custom Node/Tree types, e.g. when depth is
stored explicitly.
"""
function depth( node :: Node )
	depth = 0;
	iter_node = parent(node);
	while !isnothing( iter_node )
		depth += 1;
		iter_node = parent(iter_node)
	end
	return depth
end

# DFS Iterator
@doc """
	PreOrderDFS{T}
		start_node :: T
		traverse_subtree_filter :: Union{Nothing, F where F<:Function};
	end

An iterator to perform a depth first traversal of the subtree that has
`start_node` as its root.
PreOrder means that parent nodes are visited before their respective children.
The field `traverse_subtree_filter` defaults to `nothing`.
If a function is provided it should accept an argument of type `T` (a parent node)
and return `true` if its children should be pushed to the stack for visiting.
"""
struct PreOrderDFS{T}
	start_node :: T
	traverse_subtree_filter :: Union{Nothing, F where F<:Function};
end
# default constructor without a filter:
PreOrderDFS( sn :: T ) where T = PreOrderDFS{T}( sn, nothing );

Base.IteratorSize( ::PreOrderDFS{T} ) where T = Base.SizeUnknown();

function _children_to_stack( node :: T ) where T
    node_children = children( node )
    if isnothing( node_children )
        return T[]
    else
        return reverse(collect(node_children))
    end
end

function Base.iterate( iter :: PreOrderDFS{T} ) where T
	# in first iteration, collect all shildren as candidates to visit
	# `reverse` to imitate stack behavior
	init_stack = _children_to_stack( iter.start_node )
	(iter.start_node, init_stack);
end

function Base.iterate( iter :: PreOrderDFS{T}, stack :: AbstractVector ) where T
	if isempty( stack )
		return nothing
	else
		return_node = pop!(stack)
		if (
				isnothing( iter.traverse_subtree_filter ) ||
				iter.traverse_subtree_filter(return_node)
			)
			push!( stack, _children_to_stack(return_node)... );
		end
		return (return_node, stack );
	end
end

@doc "An abstract super type for iterators that are based on the
PreOrderDFS iterator."
abstract type PreOrderDFSSubIterator end;
Base.IteratorSize( ::PreOrderDFSSubIterator ) = Base.SizeUnknown();

# `node_filter` should be implemented by Iterators subtyping PreOrderDFSSubIterator
function node_filter( ::PreOrderDFSSubIterator ) end;

# Iteration procedure for every `PreOrderDFSSubIterator`
# (including the special cases like `Leaves` and `DepthFilterNodes` etc.)
function Base.iterate( iter::PreOrderDFSSubIterator, stack = nothing )
	if isnothing(stack)
		# first iteration …
		# initilize using the base iterator `iter.dfs_iterator`
		node, stack = iterate( iter.dfs_iterator );
	else
		if isempty(stack)
			# no more nodes left to visit; return
			return nothing
		else
			# let `iter.dfs_iterator` handle iteration
			# this takes care of "subtree filtering"
			node, stack = iterate( iter.dfs_iterator, stack);
		end
	end
	next = (node, stack);
	# as long as there is a non-nothing *next* …
	while !isnothing( next )
		node, stack = next;
		if node_filter( iter, node )
			# … return the node if it is deemed suitable …
			return node, stack
		end
		# … and progress iteration
		next = iterate( iter.dfs_iterator, stack )
	end
	return nothing
end

@doc """
	PreOrderDFSFiltered{T} <: PreOrderDFSSubIterator

An iterator that is based on the `PreOrderDFS` iterator and allows
for filtering out nodes that are returned.
The most basic constructor is
```
iter = PreOrderDFSFiltered(
	start_node :: T,
	func :: F where F <: Function )
```
The function `func` should take one argument of type `T` (a node)
and return true if it should be returned by `iter`.
"""
struct PreOrderDFSFiltered{T} <: PreOrderDFSSubIterator
	start_node :: T
	dfs_iterator :: Union{Nothing, PreOrderDFS{T}}
	node_filter :: F where F<:Function
end

# outer constructor without `dfs_iterator` provided
function PreOrderDFSFiltered( sn::T, fn :: F where F<:Function ) where T
	PreOrderDFSFiltered{T}( sn, PreOrderDFS(sn), fn )
end

# set default behavior for `node_filter` for type `PreOrderDFSFiltered`
node_filter( iter::PreOrderDFSFiltered, node ) = iter.node_filter( node )

@doc """
	Leaves{T} <: PreOrderDFSSubIterator

Initialize as per
```iter = Leaves( start_node )```
`iter` is an iterator that returns the leaves of a subtree
starting at `start_node`.
It is based on a PreOrderDFS traversal.
"""
struct Leaves{T} <: PreOrderDFSSubIterator
	start_node :: T
	dfs_iterator :: Union{Nothing, PreOrderDFS{T}}
end
Leaves( sn :: T ) where T = Leaves{T}( sn, PreOrderDFS( sn ) );

# determine a `node` a leave it it has no children
node_filter( :: Leaves, node ) = let childs = children(node);
    return isnothing(childs) || isnothing( iterate( childs ) )
end

@doc """
	DepthFilterNodes{T} <: PreOrderDFSSubIterator

An iterator based on `PreOrderDFS` that returns only nodes
with a certain depth(s). Traversal starts at `start_node` but depth
is measured by `depth( any_node )` and hence in relation to the
root of the tree (which possibly does not equal `start_node`) by default.\n
A subtree is only visited for traversal if the depth of its respective
parent is below the maximum value of desired return depths.

Initialize via
```
dfs_iterator = DepthFilterNodes( start_node, desired_depths :: Vector{Int} );
```
"""
struct DepthFilterNodes{T} <: PreOrderDFSSubIterator
	start_node :: T
	dfs_iterator :: Union{Nothing, PreOrderDFS{T}}
	desired_depths :: Union{Nothing, Vector{Int}}
end
function DepthFilterNodes( sn :: T ) where T
	return PreOrderDFS(sn)
end
function DepthFilterNodes( sn :: T, desired_depths :: Vector{Int64} ) where T
	max_depth = maximum( desired_depths )

	# setup a filter function for the base iterator so that
	# only substrees are visited where the respective parent
	# has sufficiently small depth
	children_filter = function( parent_node :: T )
		return depth(parent_node) < max_depth
	end
	base_iterator = PreOrderDFS( sn, children_filter );

	DepthFilterNodes{T}( sn, base_iterator, desired_depths );
end

function node_filter(iter::DepthFilterNodes, node )
	depth(node) in iter.desired_depths
end

#=
begin
	struct MyTree <: Node
		data :: Int64
		parent :: Union{Nothing, MyTree}
		children :: Vector{MyTree}
	end

	children( n::MyTree ) = n.children;
	parent( n :: MyTree ) = n.parent;

	function add_child!( parent :: MyTree, child_data :: Int64 )
		child_node = MyTree( child_data, parent, [] );
		push!( parent.children, child_node )
		return child_node
	end

	root_node = MyTree( 0, nothing, [] )

	child_l1_1 = add_child!( root_node, 1 )
	child_l1_2 = add_child!( root_node, 2 )
	child_l2_1 = add_child!( child_l1_1, 3 )
	child_l2_2 = add_child!( child_l1_1, 4 )

	@show [n.data for n in DepthFilterNodes(root_node, collect(0:2))];
end
=#

end


text:
using ModelingToolkit, LinearAlgebra

ModelingToolkit.@variables s[1:5], σ[1:5], V[1:5, 1:5]

F = permutedims(s-σ)*V*(s-σ)

F = simplify(F)

nobs = 2
nlat = 1

### compute K

Sig = ones(Bool, nobs, nobs)
Sig = LowerTriangular(Sig)

ind = findall(Sig)

Es = []

for i in ind
    E = zeros(Bool, nobs,nobs)
    E[i] = true
    E = vec(E)
    push!(Es, E)
end

L = hcat(Es...)

K = L*inv(transpose(L)*L)

K = convert(BitArray, K)

# compute L's

Lat = ones(Bool, nlat, nlat)
Lat = LowerTriangular(Lat)

ind = findall(Lat)

Es = []
for i in ind
    E = zeros(Bool, nlat, nlat)
    E[i] = true
    E = vec(E)
    push!(Es, E)
end

L_lat = hcat(Es...)

Obse = zeros(Bool, nobs,nobs)
Obse[diagind(Obse)] .= true

ind = findall(Obse)

Es = []
for i in ind
    E = zeros(Bool, nobs, nobs)
    E[i] = true
    E = vec(E)
    push!(Es, E)
end

L_obs = hcat(Es...)

# Define Model

ModelingToolkit.@variables λ[1:2, 1:2], ε[1:6], ϕ[1:3]

Λ =[1 0
    λ[1,1] 0
    λ[1,2] 0
    0 1
    0 λ[2,1]
    0 λ[2,2]]

Θ = [ε[1] 0 0 0 0 0
    0 ε[2] 0 0 0 0
    0 0 ε[3] 0 0 0
    0 0 0 ε[4] 0 0
    0 0 0 0 ε[5] 0
    0 0 0 0 0 ε[6]]

Φ = [ϕ[1] ϕ[3]
    ϕ[3] ϕ[2]]

transpose(K)*(kron(Λ,Λ)*L_lat*L_obs)

G = transpose(K)*(kron(Λ,Λ)*L_lat)

G = simplify.(G)

G

### in RAM notation

F = zeros(2, 3)
F[diagind(F)] .= 1.0

kron(F, F)

ModelingToolkit.@variables λ[1:2], s[1:3]

A = [0 0 λ[1]
    0 0 λ[2]
    0 0 0]

A = I + A

A = simplify.(A)

kronprod = simplify.(kron(F, F)*kron(A, A))

kronprod

simplify.(kron(F*A, F*A))

Obse = zeros(Bool, nobs+nlat, nobs+nlat)
Obse[diagind(Obse)] .= true

ind = findall(Obse)

Es = []
for i in ind
    E = zeros(Bool, nobs+nlat, nobs+nlat)
    E[i] = true
    E = vec(E)
    push!(Es, E)
end

L_obs = hcat(Es...)

G = transpose(K)*kronprod*L_obs

G = simplify.(G)

simplify.(G*s)

### in RAM notation

F = zeros(2, 3)
F[diagind(F)] .= 1.0

kron(F, F)

ModelingToolkit.@variables λ[1:2], s[1:3], A[1:3, 1:3]

A = [0 0 λ[1]
    0 0 λ[2]
    0 0 0]

A = I + A

A = simplify.(A)

kronprod = simplify.(kron(F, F)*kron(A, A))

kronprod

simplify.(kron(F*A, F*A))

Obse = zeros(Bool, nobs+nlat, nobs+nlat)
Obse[diagind(Obse)] .= true

ind = findall(Obse)

Es = []
for i in ind
    E = zeros(Bool, nobs+nlat, nobs+nlat)
    E[i] = true
    E = vec(E)
    push!(Es, E)
end

L_obs = hcat(Es...)

G = transpose(K)*kronprod*L_obs

G = simplify.(G)

simplify.(G*s)

### in RAM notation

F = zeros(2, 3)
F[diagind(F)] .= 1.0

kron(F, F)

ModelingToolkit.@variables λ[1:2], s[1:3]

A = [0 0 λ[1]
    0 0 λ[2]
    0 0 0]

A = I + A

A = simplify.(A)

kronprod = simplify.(kron(F, F)*kron(A, A))

kronprod

simplify.(kron(F*A, F*A))

Obse = zeros(Bool, nobs+nlat, nobs+nlat)
Obse[diagind(Obse)] .= true

ind = findall(Obse)

Es = []
for i in ind
    E = zeros(Bool, nobs+nlat, nobs+nlat)
    E[i] = true
    E = vec(E)
    push!(Es, E)
end

L_obs = hcat(Es...)

G = transpose(K)*kronprod*L_obs

G = simplify.(G)

simplify.(G*s)

### in RAM notation
using SparseArrays

nobs = 9
nlat = 6

F = zeros(9, 15)
F[diagind(F)] .= 1.0
F = sparse(F)
#kron(F, F)

ModelingToolkit.@variables λ[1:12], ϕ[1:6, 1:6], σ[1:9]

A =[0 0 0 0 0 0 0 0 0 1.0 0 0 1 0 0
    0 0 0 0 0 0 0 0 0 0 1 0 1 0 0
    0 0 0 0 0 0 0 0 0 0 0 1 1 0 0
    0 0 0 0 0 0 0 0 0 1 0 0 0 1 0
    0 0 0 0 0 0 0 0 0 0 1 0 0 1 0
    0 0 0 0 0 0 0 0 0 0 0 1 0 1 0
    0 0 0 0 0 0 0 0 0 1 0 0 0 0 1
    0 0 0 0 0 0 0 0 0 0 1 0 0 0 1
    0 0 0 0 0 0 0 0 0 0 0 1 0 0 1
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
]

A =[0 0 0 0 0 0 0 0 0 1 0 0 1 0 0
    0 0 0 0 0 0 0 0 0 0 1 0 0.91 0 0
    0 0 0 0 0 0 0 0 0 0 0 1 0.92 0 0
    0 0 0 0 0 0 0 0 0 0.93 0 0 0 1 0
    0 0 0 0 0 0 0 0 0 0 0.96 0 0 0.76 0
    0 0 0 0 0 0 0 0 0 0 0 0.72 0 0.72 0
    0 0 0 0 0 0 0 0 0 0.87 0 0 0 0 1
    0 0 0 0 0 0 0 0 0 0 0.76 0 0 0 0.842
    0 0 0 0 0 0 0 0 0 0 0 0.97 0 0 0.767
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
]

A =[0 0 0 0 0 0 0 0 0 1 0 0 1 0 0
    0 0 0 0 0 0 0 0 0 0 1 0 λ[7] 0 0
    0 0 0 0 0 0 0 0 0 0 0 1 λ[8] 0 0
    0 0 0 0 0 0 0 0 0 λ[1] 0 0 0 1 0
    0 0 0 0 0 0 0 0 0 0 λ[2] 0 0 λ[9] 0
    0 0 0 0 0 0 0 0 0 0 0 λ[3] 0 λ[10] 0
    0 0 0 0 0 0 0 0 0 λ[4] 0 0 0 0 1
    0 0 0 0 0 0 0 0 0 0 λ[5] 0 0 0 λ[11]
    0 0 0 0 0 0 0 0 0 0 0 λ[6] 0 0 λ[12]
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
]

A = sparse(A)

S = zeros(Num, nobs+nlat, nobs+nlat)
S[diagind(S)[1:9]] .= σ
S[13:15, 13:15] .= ϕ[4:6, 4:6]
S[10:12, 10:12] .= ϕ[1:3, 1:3]

θ = sparse(LowerTriangular(S)).nzval
θ

rank(A)

A = I + A

#A = simplify.(A)

A = inv(I-A)

kronprod = simplify.(kron(F, F)*kron(A, A))

kronprod = kron(F, F)*kron(A, A)

simplify.(kron(F*A, F*A))

Sig = ones(Bool, nobs, nobs)
Sig = LowerTriangular(Sig)
ind = findall(Sig)
Es = []
for i in ind
    E = zeros(Bool, nobs,nobs)
    E[i] = true
    E = vec(E)
    push!(Es, E)
end
L = hcat(Es...)
K = L*inv(transpose(L)*L)
K = sparse(K)

Obse = zeros(Bool, nobs+nlat, nobs+nlat)
Obse[diagind(Obse)] .= true
Obse[13:15, 13:15] .= true
Obse[10:12, 10:12] .= true
Obse = LowerTriangular(Obse)
ind = findall(Obse)

Obse

Es = []
for i in ind
    E = zeros(Bool, nobs+nlat, nobs+nlat)
    E[i] = true
    E = vec(E)
    push!(Es, E)
end

L_obs = hcat(Es...)
L_obs = sparse(L_obs)

G = transpose(K)*kronprod*L_obs

G = simplify.(Matrix(G))

G = sparse(G)

reshape(simplify.(L_obs*θ), (15, 15))[10:15, 10:15]

simplify.(G*θ)

V = rand(45,45)
V = V*V'

test = transpose(G)*V*G

using BenchmarkTools

rank(test)

test2 = V[1:27, 1:27]

@benchmark pinv(test)

@benchmark inv(test2)

@benchmark pinv(test2)

pinv(test2) == inv(test2)


text:
module EventQueue # max heap used for the event queue

using ..Diagram

abstract type Event end

mutable struct SiteEvent <: Event
   region::Diagram.Region
end

function key(ev::SiteEvent)
   return ev.region.generator[2] # y coordinate
end

function coordinates(ev::SiteEvent)
   return ev.region.generator
end

mutable struct CircleEvent <: Event
   coordinates::Tuple{Real, Real}
   center::Tuple{Real, Real}
   disappearingArc#::Union{BeachLine.Arc, Nothing} --- I'm not type checking because I'm getting a weird Julia-specific error I don't wanna solve now
                                                    # Cannot `convert` an object of type BeachLine.Arc to an object of type BeachLine.Arc
   removed::Bool
end

function CircleEvent(center::Tuple{Real, Real}, r, arc)
   return CircleEvent((center[1], center[2] - r), center, arc, false)
end

function key(ev::CircleEvent)
   return ev.coordinates[2] # y coordinate
end

function coordinates(ev::CircleEvent)
   return ev.coordinates
end

function key_gt(ev1, ev2)
    if coordinates(ev1)[2] > coordinates(ev2)[2]
        return true
    elseif coordinates(ev1)[2] == coordinates(ev2)[2]
        if coordinates(ev1)[1] < coordinates(ev2)[1]
            return true
        end
    end

    return false
end

mutable struct Heap
   data::Array{Event}
   pos::Int # first empty position in the array
end

function Heap(n::Int)
   if n == 0
      n = 1
   end

   return Heap(Array{Event}(undef, n), 1)
end

function Heap()
   return Heap(1)
end

function leftChild(i)
   return 2i
end

function rightChild(i)
   return 2i+1
end

function parent(i)
   return Int(floor(i/2))
end

function largestChild(h::Heap, i::Int)
   l = leftChild(i)
   r = rightChild(i)

   if h.pos <= l
      return -1 # no children
   elseif h.pos  == r
      return l # no right child
   elseif !key_gt(h.data[r], h.data[l])
   #elseif key(h.data[l]) >= key(h.data[r])
      return l
   else
      return r
   end
end

function push!(h::Heap, ev::Event)
   n = size(h.data)[1]

   if h.pos >= n # full
      resize!(h.data, 2n + 1) # space for another complete level
   end

   i = h.pos
   h.data[i] = ev

   p = parent(i)
   #while i > 1 && key(h.data[i]) > key(h.data[p])
   while i > 1 && key_gt(h.data[i], h.data[p])
      # bubble up
      swap = h.data[i]
      h.data[i] = h.data[p]
      h.data[p] = swap

      i = p
      p = parent(i)
   end

   h.pos += 1
end

function pop(h::Heap)
   resolved = false

   while !resolved
      if h.pos == 1
         return nothing # empty
      end

      h.pos -= 1

      # swap last element with root
      last = h.pos
      swap = h.data[1]
      h.data[1] = h.data[last]
      h.data[last] = swap

      i = 1
      c = largestChild(h, i)
      #while c != -1 && key(h.data[c]) > key(h.data[i])
      while c != -1 && key_gt(h.data[c], h.data[i])
         # trickle down
         swap = h.data[i]
         h.data[i] = h.data[c]
         h.data[c] = swap

         i = c
         c = largestChild(h, i)
      end

      if isa(h.data[h.pos], SiteEvent) || !(h.data[h.pos].removed)
         resolved = true
      end
   end

   return h.data[h.pos]
end

function remove(h::Heap, event::EventQueue.CircleEvent)
   event.removed = true
end

#function heapPropertyIsOk(h::Heap) # for debugging
#   ok = true
#
#   for i in 2:h.pos-1
#      if key(h.data[i]) > key(h.data[parent(i)])
#         ok = false
#      end
#   end
#
#   return ok
#end

end # module

max_stars_repo_name: arthursantana/coveringWithCircles

max_issues_repo_name: arthursantana/coveringWithCircles

max_forks_repo_name: arthursantana/coveringWithCircles


text:
"""
`f` from page `c.jl`.

Links:

- [`ccall`](@ref)
- [`while`](@ref)
- [`@time`](@ref)
- [`T`](@ref)
- [`f`](@ref)
- [`f(::Any)`](@ref)
- [`f(::Any, ::Any)`](@ref)
- [`f(::Any, ::Any, ::Any)`](@ref)

"""
f(x, y, z) = x + y + z


text:
module BaseModule
export Position, LQuantumNumber, MQuantumNumber, MQuantumNumbers, distance, origin, floatregex, @T_str, doublefactorial
import Base.*, Base.+, Base./, Base.-, Base.isless, Base.convert

immutable Position
	x::Float64
	y::Float64
	z::Float64
end

const origin = Position(0.,0.,0.)

"""
True String: Taking the string as is without interpolation or escape interpretation.
Necessary for regex interpolation, e.g.: Regex(T"\wo"*"rd") == r"\word" != Regex("\wo"*"rd") == r"word"
"""
macro T_str(s)
  s
end
const floatregex = T"[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?"

*(c::Real,p::Position) = Position(c*p.x,c*p.y,c*p.z)
*(p::Position,c::Real) = Position(c*p.x,c*p.y,c*p.z)
/(p::Position,c::Real) = Position(p.x/c,p.y/c,p.z/c)
+(p::Position,q::Position) = Position(p.x+q.x,p.y+q.y,p.z+q.z)
-(p::Position,q::Position) = Position(p.x-q.x,p.y-q.y,p.z-q.z)

convert(::Type{Vector{Float64}},p::Position) = [p.x,p.y,p.z]

function distance(p::Position,q::Position)
  pq = p - q
  return sqrt((pq.x)^2+(pq.y)^2+(pq.z)^2)
end

"""
MQuantumNumber describes a single orientation of an atomic orbital in cartesian space by denoting the exponent of x, y and z of the radial part
E.g.:   a d_xy orbital can be described by MQuantumNumber(1,1,0)
"""
immutable MQuantumNumber
  x::Int
  y::Int
  z::Int
end

lqnExponent = Dict{ASCIIString,Int}(
  "S" => 	0,
  "P" => 	1,
  "D" => 	2,
  "F" =>	3,
  "G" =>	4,
  "H" =>	5,
  "I" =>	6,
  "K" =>	7,
)

exponentLqn = Dict{Int,ASCIIString}(
  0 => 	"S",
  1 => 	"P",
  2 => 	"D",
  3 =>	"F",
  4 =>	"G",
  5 =>	"H",
  6 =>	"I",
  7 =>	"K",
)

"""
The LQuantumNumber desribes whether a function is spherical (s, lqn=0), polarized (p, lqn=1), etc.
As is generally known the total MQuantumNumber in pure coordinates lies between -lqn and +lqn.
This means that the sum of the three cartesian components can at most b. LQuantumNumber.
"""
immutable LQuantumNumber
	symbol::ASCIIString
	exponent::Int

	LQuantumNumber(sym::AbstractString) = new(sym,lqnExponent[sym])
	LQuantumNumber(exp::Int) = new(exponentLqn[exp],exp)
end

"""
MQuantumNumber**s** returns all MQuantumNumber objects corresponding to the LQuantumNumber object given to the constructor as an iterable collection.
Order is for example xx,xy,xz,yy,yz,zz for LQuantumNumber = 2
"""
immutable MQuantumNumbers
  lqn::LQuantumNumber
  mqnarray::Array{MQuantumNumber,1}
  count::Int

  function MQuantumNumbers(lqn::LQuantumNumber)
    maxmqn = lqn.exponent
    mqnarray=MQuantumNumber[]
		for x in maxmqn:-1:0
			for y in (maxmqn-x):-1:0
				z = maxmqn-x-y
				append!(mqnarray,[MQuantumNumber(x,y,z)])
			end
		end
		count = ((maxmqn+1)^2+(maxmqn+1))/2
		new(lqn,mqnarray,count)
	end
end
Base.start(mqns::MQuantumNumbers) = 1
Base.next(mqns::MQuantumNumbers,state) = (mqns.mqnarray[state],state+1)
Base.done(mqns::MQuantumNumbers,state) = state > mqns.count
Base.length(mqns::MQuantumNumbers) = mqns.count
#Base.eltype(::Type{MQuantumNumbers}) = Int

isless(lqn1::LQuantumNumber, lqn2::LQuantumNumber) = lqn1.exponent<lqn2.exponent

function evaluateFunction(x::Position, f::Function)
  return f(x)
end

doublefactorial(n::Int) = prod(n:-2:1)

end # module


text:
using Documenter, Aswing

makedocs(
    sitename = "Aswing",
    modules = [Aswing],
    format=Documenter.HTML(prettyurls = get(ENV, "CI", nothing)=="true"),
    pages   = ["Home" => "index.md",
        "Installation" => "install.md",
        "Example Usage" => "example.md",
        "Library" => "library.md"]
)

# deploydocs(
#     root = ".",
#     target = "./build",
#     repo   = "github.com/byuflowlab/Aswing.jl",
#     branch = "gh-pages",
#     devbranch = "master"
# )


text:
module twisting_circle_examples
# Twisting of a circle into shape with two turns
# G. Rebel; Finite Rotation Shell Theory including Drill Rotations and its
# Finite Element Implementation; PhD Thesis, Delft University Press (1998).

using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CrossSectionModule: CrossSectionRectangle
using FinEtoolsFlexStructures.MeshFrameMemberModule: frame_member, merge_members
using FinEtoolsFlexStructures.FEMMCorotBeamModule
using FinEtoolsFlexStructures.FEMMCorotBeamModule: FEMMCorotBeam
stiffness = FEMMCorotBeamModule.stiffness
geostiffness = FEMMCorotBeamModule.geostiffness
mass = FEMMCorotBeamModule.mass
distribloads_global = FEMMCorotBeamModule.distribloads_global
restoringforce = FEMMCorotBeamModule.restoringforce
gyroscopic = FEMMCorotBeamModule.gyroscopic
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using LinearAlgebra: dot
using Arpack
using LinearAlgebra
using SparseArrays
using VisualStructures: plot_points, plot_nodes, plot_midline, render, plot_space_box, plot_solid, space_aspectratio
using PlotlyJS
using JSON

num3dig(a) = round(a * 1000) / 1000

function twisting_circle()
    # Parameters:
    E = 200000.0;
    nu = 0.3;# Poisson ratio
    rho = 7.8e-9;
    b = 0.6; h = 6.0;  # cross-sectional dimensions
    radius = 120 # radius of the circle
    utol = 1e-3;
    maxit = 30
    rotmag = 2*pi
    # Select the number of elements around the circumference.
    nel = 64;

    # Cross-sectional properties
    cs = CrossSectionRectangle(s -> h, s -> b, s -> [0.0, 0.0, 1.0])

    tolerance = radius/nel/1000;
    fens, fes = frame_member([0 0 0; 2*pi 0 0], nel, cs)
    for i in 1:count(fens)
        a = fens.xyz[i, 1]
        fens.xyz[i, :] .= (radius+radius*cos(a), radius*sin(a), 0)
    end
    fens, fes = mergenodes(fens, fes, tolerance, [1, nel+1])

    # Material properties
    material = MatDeforElastIso(DeforModelRed3D, 0.0, E, nu, 0.0)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Locate some nodes
    box = fill(0.0, 6)
    # The pinned end is
    initbox!(box, [0.0,0.0,0.0])
    clampedn = selectnode(fens; box = box, inflate = tolerance)
    # The force is applied at the elbow
    initbox!(box, [2*radius,0,0])
    torquen = selectnode(fens; box = box, inflate = tolerance)

    # Apply EBC's
    for i in [1, 2, 3, 4, 5, 6]
        setebc!(dchi, clampedn, true, i)
    end
    for i in [2, 3, 4, 5, 6]
        setebc!(dchi, torquen, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Additional fields
    u1 = deepcopy(u0)
    Rfield1 = deepcopy(Rfield0)
    rhs = gathersysvec(dchi);
    TMPv = deepcopy(rhs)
    utol = 1e-13*dchi.nfreedofs;

    femm = FEMMCorotBeam(IntegDomain(fes, GaussRule(1, 2)), material)

    tbox = plot_space_box([[0 -radius -radius]; [2.5*radius radius radius]])
    tshape0 = plot_solid(fens, fes; x = geom0.values, u = 0.0.*dchi.values[:, 1:3], R = Rfield0.values, facecolor = "rgb(125, 155, 125)", opacity = 0.3);
    plots = cat(tbox,  tshape0; dims = 1)
    pl = render(plots; width = 800, height = 800)
    sleep(0.5)
    # savefig(pl, "twisting_circle-$(0).png", width = 300, height = 300)

    load_parameter_delta = 0.002
    load_parameters = 0.0:load_parameter_delta:1;
    step = 0
    for load_parameter in load_parameters
        dchi.values[:] .= 0.0
        dchi.fixed_values[torquen[1], 4] = rotmag * load_parameter_delta
        applyebc!(dchi) # Apply boundary conditions
        u1.values[:] = u0.values[:]; # guess
        Rfield1.values[:] = Rfield0.values[:]; # guess
        update_rotation_field!(Rfield1, dchi)

        println("Load: $load_parameter")
        iter = 1;
        while true
            Fr = restoringforce(femm, geom0, u1, Rfield1, dchi);       # Internal forces
            rhs = Fr;
            K = stiffness(femm, geom0, u1, Rfield1, dchi) + geostiffness(femm, geom0, u1, Rfield1, dchi);
            dchi = scattersysvec!(dchi, (K)\rhs); # Disp. incr
            dchi.values[torquen[1], 4] = 0.0
            u1.values[:] += (dchi.values[:,1:3])[:];   # increment displacement
            update_rotation_field!(Rfield1, dchi)
            print("$iter: ||du||=$(maximum(abs.(dchi.values[:])))\n")
            if maximum(abs.(dchi.values[:])) < utol # convergence check
                break;
            end
            if (iter > maxit)# bailout for failed convergence
                error("Possible failed convergence");
            end
            iter += 1;
        end
        u0.values[:] = u1.values[:];       # update the displacement
        Rfield0.values[:] = Rfield1.values[:]; # update the rotations

        step = step + 1
        if mod(step, 5) == 1
            @show step
            tshape1 = plot_solid(fens, fes; x = geom0.values, u = u1.values, R = Rfield1.values, facecolor = "rgb(125, 15, 15)");
            plots = cat(tbox,  tshape0,  tshape1; dims = 1)
            pl.plot.layout[:title] = "Angle=$(num3dig(rotmag * load_parameter / pi))*pi"
            react!(pl, plots, pl.plot.layout)
            sleep(0.25)
            # savefig(pl, "twisting_circle-$(step).png", width = 300, height = 300)
        end

    end

    return true
end # twisting_circle

function twisting_circle_frames()
    # Parameters:
    E = 200000.0;
    nu = 0.3;# Poisson ratio
    rho = 7.8e-9;
    b = 0.6; h = 6.0;  # cross-sectional dimensions
    radius = 120 # radius of the circle
    utol = 1e-3;
    maxit = 30
    rotmag = 2*pi
    # Select the number of elements around the circumference.
    nel = 64;

    # Cross-sectional properties
    cs = CrossSectionRectangle(s -> h, s -> b, s -> [0.0, 0.0, 1.0])

    tolerance = radius/nel/1000;
    fens, fes = frame_member([0 0 0; 2*pi 0 0], nel, cs)
    for i in 1:count(fens)
        a = fens.xyz[i, 1]
        fens.xyz[i, :] .= (radius+radius*cos(a), radius*sin(a), 0)
    end
    fens, fes = mergenodes(fens, fes, tolerance, [1, nel+1])

    # Material properties
    material = MatDeforElastIso(DeforModelRed3D, 0.0, E, nu, 0.0)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Locate some nodes
    box = fill(0.0, 6)
    # The pinned end is
    initbox!(box, [0.0,0.0,0.0])
    clampedn = selectnode(fens; box = box, inflate = tolerance)
    # The force is applied at the elbow
    initbox!(box, [2*radius,0,0])
    torquen = selectnode(fens; box = box, inflate = tolerance)

    # Apply EBC's
    for i in [1, 2, 3, 4, 5, 6]
        setebc!(dchi, clampedn, true, i)
    end
    for i in [2, 3, 4, 5, 6]
        setebc!(dchi, torquen, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Additional fields
    u1 = deepcopy(u0)
    Rfield1 = deepcopy(Rfield0)
    rhs = gathersysvec(dchi);
    TMPv = deepcopy(rhs)
    utol = 1e-13*dchi.nfreedofs;

    femm = FEMMCorotBeam(IntegDomain(fes, GaussRule(1, 2)), material)

    tbox = plot_space_box([[0 -radius -radius]; [2.5*radius radius radius]])
    tshape0 = plot_solid(fens, fes; x = geom0.values, u = 0.0.*dchi.values[:, 1:3], R = Rfield0.values, facecolor = "rgb(125, 155, 125)", opacity = 0.3);
    plots = cat(tbox,  tshape0; dims = 1)
    pl = render(plots; width = 800, height = 800)
    sleep(0.5)
    savefig(pl, "twisting_circle-$(0).png", width = 300, height = 300)

    load_parameter_delta = 0.002
    load_parameters = 0.0:load_parameter_delta:1;
    step = 0
    for load_parameter in load_parameters
        dchi.values[:] .= 0.0
        dchi.fixed_values[torquen[1], 4] = rotmag * load_parameter_delta
        applyebc!(dchi) # Apply boundary conditions
        u1.values[:] = u0.values[:]; # guess
        Rfield1.values[:] = Rfield0.values[:]; # guess
        update_rotation_field!(Rfield1, dchi)

        println("Load: $load_parameter")
        iter = 1;
        while true
            Fr = restoringforce(femm, geom0, u1, Rfield1, dchi);       # Internal forces
            rhs = Fr;
            K = stiffness(femm, geom0, u1, Rfield1, dchi) + geostiffness(femm, geom0, u1, Rfield1, dchi);
            dchi = scattersysvec!(dchi, (K)\rhs); # Disp. incr
            dchi.values[torquen[1], 4] = 0.0
            u1.values[:] += (dchi.values[:,1:3])[:];   # increment displacement
            update_rotation_field!(Rfield1, dchi)
            print("$iter: ||du||=$(maximum(abs.(dchi.values[:])))\n")
            if maximum(abs.(dchi.values[:])) < utol # convergence check
                break;
            end
            if (iter > maxit)# bailout for failed convergence
                error("Possible failed convergence");
            end
            iter += 1;
        end
        u0.values[:] = u1.values[:];       # update the displacement
        Rfield0.values[:] = Rfield1.values[:]; # update the rotations

        step = step + 1
        if mod(step, 5) == 1
            @show step
            tshape1 = plot_solid(fens, fes; x = geom0.values, u = u1.values, R = Rfield1.values, facecolor = "rgb(125, 15, 15)");
            plots = cat(tbox,  tshape0,  tshape1; dims = 1)
            pl.plot.layout[:title] = "Angle=$(num3dig(rotmag * load_parameter / pi))*pi"
            react!(pl, plots, pl.plot.layout)
            sleep(0.5)
            savefig(pl, "twisting_circle-$(step).png", width = 300, height = 300)
        end

    end

    return true
end # twisting_circle_frames

function twisting_circle_frames_alt()
    # Parameters:
    E = 200000.0;
    nu = 0.3;# Poisson ratio
    rho = 7.8e-9;
    b = 0.6; h = 6.0;  # cross-sectional dimensions
    radius = 120 # radius of the circle
    utol = 1e-3;
    maxit = 30
    rotmag = 2*pi
    # Select the number of elements around the circumference.
    nel = 64;

    # Cross-sectional properties
    cs = CrossSectionRectangle(s -> h, s -> b, s -> [0.0, 0.0, 1.0])

    tolerance = radius/nel/1000;
    fens, fes = frame_member([0 0 0; 2*pi 0 0], nel, cs)
    for i in 1:count(fens)
        a = fens.xyz[i, 1]
        fens.xyz[i, :] .= (radius+radius*cos(a), radius*sin(a), 0)
    end
    fens, fes = mergenodes(fens, fes, tolerance, [1, nel+1])

    # Material properties
    material = MatDeforElastIso(DeforModelRed3D, 0.0, E, nu, 0.0)

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Locate some nodes
    box = fill(0.0, 6)
    # The pinned end is
    initbox!(box, [0.0,0.0,0.0])
    clampedn = selectnode(fens; box = box, inflate = tolerance)
    # The force is applied at the elbow
    initbox!(box, [2*radius,0,0])
    torquen = selectnode(fens; box = box, inflate = tolerance)

    # Apply EBC's
    for i in [1, 2, 3, 4, 5, 6]
        setebc!(dchi, clampedn, true, i)
    end
    for i in [2, 3, 4, 5, 6]
        setebc!(dchi, torquen, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Additional fields
    u1 = deepcopy(u0)
    Rfield1 = deepcopy(Rfield0)
    rhs = gathersysvec(dchi);
    TMPv = deepcopy(rhs)
    utol = 1e-13*dchi.nfreedofs;

    femm = FEMMCorotBeam(IntegDomain(fes, GaussRule(1, 2)), material)

    tbox = plot_space_box([[0 -radius -radius]; [2.05*radius radius radius]])
    tshape0 = plot_solid(fens, fes; x = geom0.values, u = 0.0.*dchi.values[:, 1:3], R = Rfield0.values, facecolor = "rgb(125, 155, 125)", opacity = 0.3);
    plots = cat(tbox,  tshape0; dims = 1)

    layout = Layout(width="500", height="500", autosize=true, title="title",
        scene=attr(
            xaxis=attr(title="X", gridcolor="rgb(255, 255, 255)",
              zerolinecolor="rgb(255, 255, 255)",
              showbackground=true,
              backgroundcolor="rgb(255, 255, 255)"),
            yaxis=attr(title="Y", gridcolor="rgb(255, 255, 255)",
             zerolinecolor="rgb(255, 255, 255)",
             showbackground=true,
             backgroundcolor="rgb(255, 255, 255)"),
            zaxis=attr(title="Z", gridcolor="rgb(255, 255, 255)",
             zerolinecolor="rgb(255, 255, 255)",
             showbackground=true,
             backgroundcolor="rgb(255, 255, 255)"),
            camera = attr(eye = attr(x = 0.98, y = 0.9, z = 0.1),
                center = attr(x = 0, y = 0, z = 0),
                up = attr(x = 0, y = 0, z = 1)),
            aspectmode = "manual"),
        showlegend=false,
        )
    pl = plot(plots, layout)
    display(pl)
    sleep(0.5)
    # savejson(pl, "plot.json")
    # return pl

    load_parameter_delta = 0.004
    load_parameters = 0.0:load_parameter_delta:1;
    step = 0
    for load_parameter in load_parameters
        dchi.values[:] .= 0.0
        dchi.fixed_values[torquen[1], 4] = rotmag * load_parameter_delta
        applyebc!(dchi) # Apply boundary conditions
        u1.values[:] = u0.values[:]; # guess
        Rfield1.values[:] = Rfield0.values[:]; # guess
        update_rotation_field!(Rfield1, dchi)

        println("Load: $load_parameter")
        iter = 1;
        while true
            Fr = restoringforce(femm, geom0, u1, Rfield1, dchi);       # Internal forces
            rhs = Fr;
            K = stiffness(femm, geom0, u1, Rfield1, dchi) + geostiffness(femm, geom0, u1, Rfield1, dchi);
            dchi = scattersysvec!(dchi, (K)\rhs); # Disp. incr
            dchi.values[torquen[1], 4] = 0.0
            u1.values[:] += (dchi.values[:,1:3])[:];   # increment displacement
            update_rotation_field!(Rfield1, dchi)
            print("$iter: ||du||=$(maximum(abs.(dchi.values[:])))\n")
            if maximum(abs.(dchi.values[:])) < utol # convergence check
                break;
            end
            if (iter > maxit)# bailout for failed convergence
                error("Possible failed convergence");
            end
            iter += 1;
        end
        u0.values[:] = u1.values[:];       # update the displacement
        Rfield0.values[:] = Rfield1.values[:]; # update the rotations

        step = step + 1
        if mod(step, 5) == 1
            @show step
            tshape1 = plot_solid(fens, fes; x = geom0.values, u = u1.values, R = Rfield1.values, facecolor = "rgb(125, 15, 15)");
            plots = cat(tbox,  tshape0,  tshape1; dims = 1)
            pl.plot.layout[:title] = "Angle=$(num3dig(rotmag * load_parameter / pi))*pi"
            react!(pl, plots, pl.plot.layout)
            sleep(0.5)
            savefig(pl, "twisting_circle-$(step).png", width = 300, height = 300)
        end

    end

    return true
end # twisting_circle_frames_alt

function allrun()
    println("#####################################################")
    println("# twisting_circle ")
    twisting_circle()
    println("#####################################################")
    println("# twisting_circle_frames ")
    twisting_circle_frames()
    println("#####################################################")
    println("# twisting_circle_frames_alt ")
    twisting_circle_frames_alt()
    return true
end # function allrun

@info "All examples may be executed with "
println("using .$(@__MODULE__); $(@__MODULE__).allrun()")

end # module
nothing


text:
function main()
    # Define beam
    λ = LinRange(200,5000,1000) # wavelength range [nm]
    θ = LinRange(0,90,90)
    beam = PlaneWave(λ,θ)
    # Define layers
    l0 = LayerTMMO(RIdb.air(beam.λ))
    l1 = LayerTMMO(RIdb.dummy(beam.λ,1.45,0.0); type=:OT, d=1/4.)
    l2 = LayerTMMO(RIdb.dummy(beam.λ,3.45,0.0); type=:OT, d=1/4.)
    layers = [l0, l1, l2, l1, l2, l1, l2, l1, l2, l0]
    # Reference wavelenth
    λ0 = 700.0
    # call main script
    sol = tmm_optics(beam, layers; λ0=λ0, emfflag=true, h=10, pbgflag=true)
    return sol
end

sol = main()