text:
using DimensionalData
import DimensionalData: dims, refdims, data, name, metadata, rebuild, parent, formatdims, AbstractDimArray
import Base.Array

abstract type AbstractFeatureArray{T,N,D,A} <: AbstractDimArray{T,N,D,A} end
export AbstractFeatureArray

AbstractFeatureVector = AbstractFeatureArray{T, 1} where {T}
AbstractFeatureMatrix = AbstractFeatureArray{T, 2} where {T}
export AbstractFeatureVector, AbstractFeatureMatrix

"""
    F = FeatureArray(data::AbstractArray, features::Union{Tuple{Symbol},Vector{Symbol}}, [timeseries::Union{Vector, Tuple}], args...)

Construct a `FeatureArray`, which annotates the array `data` with names of `features` along rows and, optionally, `timeseries` along columns.
Since `FeatureArray <: AbstractFeatureArray <: AbstractDimArray`, further arguments to the `FeatureArray` constructor are passed to the `DimArray` constructor.
To access feature names, use `getnames(F)`.

# Examples
```julia-repl
data = rand(Int, 2, 10) # Some feature matrix with 2 features and 10 timeseries
F = FeatureArray(data, [:sum, :length])
```
"""
struct FeatureArray{T,N,D<:Tuple,R<:Tuple,A<:AbstractArray{T,N},Na,Me} <: AbstractFeatureArray{T,N,D,A}
    data::A
    dims::D
    refdims::R
    name::Na
    metadata::Me
end
export FeatureArray

function FeatureArray(data::A, dims::D, refdims::R=(), name::Na=NoName()) where {D,R,A,Na}
    if typeof(dims[1]) <: Dim{:feature, Vector{Symbol}, A, B} where {A, B}
        FeatureArray(data, formatdims(data, dims), refdims, name, NoMetadata())
    else
        @error "Incorrect dimensions for FeatureArray"
    end
end

function FeatureArray(data::AbstractArray, features::Union{Tuple{Symbol}, Vector{Symbol}}, args...)
    FeatureArray(data, (Dim{:feature}(features), fill(AnonDim(), ndims(data)-1)...), args...)
end
function FeatureArray(data::AbstractArray, features::Union{Tuple{Symbol}, Vector{Symbol}}, timeseries::Union{Vector, Tuple}, args...)
    if data isa AbstractVector
        FeatureArray(reshape(data, :, 1), (Dim{:feature}(features), Dim{:timeseries}(timeseries)), args...)
    else
        FeatureArray(data, (Dim{:feature}(features), Dim{:timeseries}(timeseries)), args...)
    end
end

dims(A::AbstractFeatureArray) = A.dims
refdims(A::AbstractFeatureArray) = A.refdims
data(A::AbstractFeatureArray) = A.data
name(A::AbstractFeatureArray) = A.name
metadata(A::AbstractFeatureArray) = A.metadata
parent(A::AbstractFeatureArray) = data(A)
Base.Array(A::AbstractFeatureArray) = Array(parent(A))

@inline function rebuild(A::AbstractFeatureArray, data::AbstractArray, dims::Tuple, refdims::Tuple, name, metadata)
    FeatureArray(data, dims, refdims, name, metadata)
end

getindex(A::AbstractFeatureArray, 𝑓::AbstractFeature, I...) = getindex(A, getname(𝑓), I...)

getindex(A::AbstractFeatureArray, 𝒇::AbstractFeatureSet, I...) = getindex(A, getnames(𝒇), I...)

"""
    FeatureArray{T, 2} where {T}

An alias to construct a `FeatureArray` for a flat set of timeseries.

# Examples
```julia-repl
data = rand(Int, 2, 3) # Some feature matrix with 2 features and 3 timeseries
F = FeatureMatrix(data, [:sum, :length], [1, 2, 3])
```
"""
FeatureMatrix = FeatureArray{T, 2} where {T}
featureMatrix = FeatureMatrix
export FeatureMatrix, featureMatrix

"""
    FeatureArray{T, 1} where {T}

An alias to construct a `FeatureArray` for a single time series.

# Examples
```julia-repl
data = randn(2) # Feature values for 1 time series
𝐟 = FeatureVector(data, [:sum, :length])
```
"""
FeatureVector = FeatureArray{T, 1} where {T}
featureVector = FeatureVector
export FeatureVector, featureVector

FeatureArray(X::AbstractArray, 𝒇::AbstractFeatureSet) = FeatureArray(X::AbstractArray, getnames(𝒇))

(FeatureArray{T,N} where {T})(x::AbstractArray{S,N}, args...) where {S,N} = FeatureArray(x, args...)

"""
    getnames(𝒇::FeatureArray)
Get the names of features represented in the feature vector or array 𝒇 as a vector of symbols.
"""
featureDims(A::AbstractDimArray) = dims(A, :feature).val
getnames(A::AbstractFeatureArray) = featureDims(A)
export getnames

timeseriesDims(A::AbstractDimArray) = dims(A, :timeseries).val


text:
#!/bin/bash
# -*- mode: julia -*-
#=
JULIA="${JULIA:-julia}"
JULIA_CMD="${JULIA_CMD:-$JULIA --color=yes --startup-file=no}"
exec $JULIA_CMD -e 'include(popfirst!(ARGS))' "${BASH_SOURCE[0]}" "$@"
=#
using JuliaCLI
JuliaCLI.serve()

ext: jl


text:
"""
"""
abstract type FESpace <: GridapType end

# Minimal FE interface (used by FEOperator)

"""
"""
function num_free_dofs(f::FESpace)
  @abstractmethod
end

"""
"""
function FEFunction(fe::FESpace, free_values)
  @abstractmethod
end

function CellData.CellField(fe::FESpace,cell_vals)
  @abstractmethod
end

function EvaluationFunction(fe::FESpace, free_values)
  FEFunction(fe,free_values)
end

"""
"""
function zero_free_values(fs::FESpace)
  @abstractmethod
end

"""
"""
function Base.zero(f::FESpace)
  free_values = zero_free_values(f)
  FEFunction(f,free_values)
end

# Extended FEInterface used by FEOperatorFromTerms and Assemblers

"""
"""
function constraint_style(::Type{<:FESpace})
  @abstractmethod
end

constraint_style(f::T) where T<:FESpace = constraint_style(T)

"""
"""
function has_constraints(::Type{T}) where T <:FESpace
  v = constraint_style(T)
  get_val_parameter(v)
end

has_constraints(f::T) where T<:FESpace = has_constraints(T)

"""
"""
function get_cell_dofs(f::FESpace)
  @abstractmethod
end

function get_cell_dofs(f::FESpace,cellids::AbstractArray)
  reindex(get_cell_dofs(f),cellids)
end

function get_cell_dofs(f::FESpace,cellids::SkeletonPair)
  ids = reindex(get_cell_dofs(f),cellids)
  axs = reindex(get_cell_axes_with_constraints(f),cellids)
  merge_cell_dofs_at_skeleton(ids.left,ids.right,axs.left,axs.right)
end

"""
"""
function get_cell_basis(f::FESpace)
  @abstractmethod
end

"""
"""
function CellData.get_cell_axes(f::FESpace)
  @abstractmethod
end

"""
"""
function get_cell_axes_with_constraints(f::FESpace)
  _get_cell_axes_with_constraints(f,constraint_style(f))
end

function _get_cell_axes_with_constraints(f,::Val{false})
  get_cell_axes(f)
end

function _get_cell_axes_with_constraints(f,::Val{true})
  @abstractmethod
end

function get_cell_constraints(f::FESpace)
  _get_cell_constraints(f,constraint_style(f))
end

function _get_cell_constraints(f,::Val{false})
  identity_constraints(get_cell_axes(f))
end

function _get_cell_constraints(f,::Val{true})
  @abstractmethod
end

function get_cell_constraints(f::FESpace,cellids::AbstractArray)
  reindex(get_cell_constraints(f),cellids)
end

function get_cell_constraints(f::FESpace,cellids::SkeletonPair)
  constr = reindex(get_cell_constraints(f),cellids)
  axsrows = reindex(get_cell_axes_with_constraints(f),cellids)
  axscols = reindex(get_cell_axes(f),cellids)
  merge_cell_constraints_at_skeleton(
    constr.left,constr.right,
    axsrows.left,axsrows.right,
    axscols.left,axscols.right)
end

function get_cell_isconstrained(f::FESpace)
  _get_cell_isconstrained(f,constraint_style(f))
end

function _get_cell_isconstrained(f,::Val{false})
  Fill(false,length(get_cell_dofs(f)))
end

function _get_cell_isconstrained(f,::Val{true})
  @abstractmethod
end

function get_cell_isconstrained(f::FESpace,cellids::AbstractArray)
  reindex(get_cell_isconstrained(f),cellids)
end

function get_cell_isconstrained(f::FESpace,cellids::SkeletonPair)
  isconstr = reindex(get_cell_isconstrained(f),cellids)
  apply((l,r)-> l||r,isconstr.left,isconstr.right)
end

function CellData.attach_constraints_rows(f::FESpace,cellarr,cellids)
  _attach_constraints_rows(f,cellarr,cellids,constraint_style(f))
end

function _attach_constraints_rows(f::FESpace,cellarr,cellids,::Val{false})
  cellarr
end

function _attach_constraints_rows(f::FESpace,cellarr,cellids,::Val{true})
  cellconstr = get_cell_constraints(f,cellids)
  cellmask = get_cell_isconstrained(f,cellids)
  attach_constraints_rows(cellarr,cellconstr,cellmask)
end

function CellData.attach_constraints_cols(f::FESpace,cellarr,cellids)
  _attach_constraints_cols(f,cellarr,cellids,constraint_style(f))
end

function _attach_constraints_cols(f::FESpace,cellarr,cellids,::Val{false})
  cellarr
end

function _attach_constraints_cols(f::FESpace,cellarr,cellids,::Val{true})
  cellconstr = get_cell_constraints(f,cellids)
  cellmask = get_cell_isconstrained(f,cellids)
  attach_constraints_cols(cellarr,cellconstr,cellmask)
end

"""
"""
function test_fe_space(f::FESpace)
  free_values = zero_free_values(f)
  @test length(free_values) == num_free_dofs(f)
  fe_function = FEFunction(f,free_values)
  test_fe_function(fe_function)
  fe_basis = get_cell_basis(f)
  @test isa(has_constraints(f),Bool)
  @test isa(has_constraints(typeof(f)),Bool)
  @test length(get_cell_dofs(f)) == length(fe_basis)
  @test length(get_cell_axes(f)) == length(fe_basis)
  @test length(get_cell_axes_with_constraints(f)) == length(fe_basis)
  @test length(get_cell_constraints(f)) == length(fe_basis)
  @test length(get_cell_isconstrained(f)) == length(fe_basis)
  @test CellField(f,get_cell_dofs(f)) != nothing
end

function test_fe_space(f::FESpace,matvecdata,matdata,vecdata)
  test_fe_space(f)

  cellmat, cellidsrows, cellidscols = matdata
  cm = attach_constraints_cols(f,cellmat,cellidscols)
  if ! has_constraints(f)
    @test cm === cellmat
  end
  cm = attach_constraints_rows(f,cellmat,cellidsrows)
  if ! has_constraints(f)
    @test cm === cellmat
  end

  cellvec, cellidsrows = vecdata
  cv = attach_constraints_rows(f,cellvec,cellidsrows)
  if ! has_constraints(f)
    @test cv === cellvec
  end

  cellmatvec, cellidsrows, cellidscols = matvecdata
  cmv = attach_constraints_cols(f,cellmatvec,cellidscols)
  if ! has_constraints(f)
    @test cmv === cellmatvec
  end
  cmv = attach_constraints_rows(f,cellmatvec,cellidsrows)
  if ! has_constraints(f)
    @test cmv === cellmatvec
  end

end


text:
using eFEMpart
include("Laplace_example_param.jl")

function laplace_example(p::Lap_Para)
  # define mesh through eFEMpart code
  mesh = squareMesh([p.xmin,p.xmax,p.ymin,p.ymax],p.Nintervals,p.basisorder)

  # define variable name and operator type
  OperatorType = :Poisson2D

  # Define auxiliary information
  dNodes = Dirichlet(:left,:right)
  nNodes = Neumann(:top,:bottom)

  # Define dirichlet, neumann, and forcing functions
  dBCf = Dirichlet((x,y) -> (x==p.xmin ? p.lBCval : p.rBCval))
  nBCf = Neumann((x,y) -> 0.0)
  ff   = Forcing((x,y) -> 0.0)

  # required technicality for grouping auxiliary information
  Nodes = [dNodes,nNodes]; bcfun = [dBCf,nBCf,ff]

  # define problem variable
  prob = Problem(mesh,Nodes,bcfun,OperatorType)

  # solve problem
  sol = solve(prob,mesh)

  # save solution
  vtkname = Path(p.solFolder,p.solFile)
  sd = ScalarData(sol.u)
  sn = ScalarNames(p.solName)
  vtksave(mesh,sd,sn,vtkname)
end

function main()
  p = Lap_Para()
  laplace_example(p)
end

main()

max_stars_repo_name: pseastham/eFEMpart

max_issues_repo_name: pseastham/eFEMpart


text:
# This is a lightly-modified version of code in DocumenterTools.jl (thanks!)

"""
    pubkey, privkey = genkeys(; keycomment="Documenter")

Generate a public/private key pair. `pubkey` can be used as a [`deploykey`](@ref) and `privkey` as a [`secret`](@ref).
"""
function genkeys(; keycomment="Documenter")
    # Error checking. Do the required programs exist?
    if Sys.iswindows()
        success(`where where`)      || error("'where' not found.")
        success(`where ssh-keygen`) || error("'ssh-keygen' not found.")
    else
        success(`which which`)      || error("'which' not found.")
        success(`which ssh-keygen`) || error("'ssh-keygen' not found.")
    end

    directory = pwd()
    filename  = "github-private-key"

    isfile(filename) && error("temporary file '$(filename)' already exists in working directory")
    isfile("$(filename).pub") && error("temporary file '$(filename).pub' already exists in working directory")

    # Generate the ssh key pair.
    success(`ssh-keygen -N "" -C $keycomment -f $filename`) || error("failed to generate a SSH key pair.")

    pubkey = chomp(read("$filename.pub", String))
    privkey = base64encode(read(filename, String))
    rm("$filename.pub")
    rm(filename)

    return pubkey, privkey
end


text:
module Roassal
using Infiltrator

# ------------------------------------
# Graphic
using Gtk, Graphics
# ------------------------------------
# Modeling
export Shape
export pos, extent, compute_encompassing_rectangle
export translate_to!, extent!

export RBox, get_color, set_color!

export RLine

export RColor, RColor_BLUE, RColor_GREEN, RColor_RED

export RCanvas
export number_of_shapes, add!, remove_shape!, rshow
export rendererVisitor
export get_shape_at_position
export offset_from_canvas_to_screen, offsetFromScreenToCanvas
export get_shapes

export Callback
export numberOfCallbacks, add_callback!, trigger_callback

export highlightable
export popup
# ------------------------------------
"""
Shape definitions
"""
abstract type Shape end
abstract type BoundedShape <: Shape end
# ------------------------------------
"""
Model a box
"""
mutable struct RBox <: BoundedShape
    color
    x
    y
    width
    height
    callbacks
    canvas
    outgoing_edges
    incoming_edges
end
RBox(;color=RColor(),
        x=0, y=0,
        width=10, height=10) = RBox(color, x, y, width, height, [], nothing, [], [])

pos(s::BoundedShape) = (s.x, s.y)
extent(s::BoundedShape) = (s.width, s.height)

function set_color!(s::Shape, color)
    s.color = color
end

function get_color(s::Shape)
    return s.color
end

function extent!(s::BoundedShape, width, height)
    s.width = width
    s.height = height
end

function translate_to!(s::BoundedShape, p::Tuple{Number, Number})
    translate_to!(s, p[1], p[2])
end

function translate_to!(s::BoundedShape, x::Number, y::Number)
    s.x = x
    s.y = y
end

# Return (x, y, w, h)
function compute_encompassing_rectangle(s::BoundedShape)
    x = s.x - (s.width / 2)
    y = s.y - (s.height / 2)
    return (x, y, s.width, s.height)
end

# ------------------------------------

mutable struct RLine <: Shape
    from
    to
    color
    canvas
end

function RLine(from::Shape, to::Shape; color=RColor_BLUE)
    a_line = RLine(from, to, color, nothing)
    push!(from.outgoing_edges, a_line)
    push!(to.incoming_edges, a_line)
    return a_line
end

# ------------------------------------
"""
Color
"""
mutable struct RColor
    r
    g
    b
end
RColor() = RColor(0.8, 0.8, 0.8)

RColor_BLUE = RColor(0, 0, 1)
RColor_GREEN = RColor(0, 1, 0)
RColor_RED = RColor(1, 0, 0)

# ------------------------------------
"""
Canvas
"""
mutable struct RCanvas
    shapes::Array{Shape}
    callbacks
    shapeBeingPointed
end
RCanvas() = RCanvas([], [], RBox())

number_of_shapes(c::RCanvas) = length(c.shapes)

function add!(c::RCanvas, s::Shape)
    push!(c.shapes, s)
    s.canvas = c
end

get_shapes(c::RCanvas) = c.shapes

function redraw(canvas::RCanvas, c::GtkCanvas)
    @guarded draw(c) do widget
        h = height(c)
        w = width(c)
        ctx = getgc(c)
        rectangle(ctx, 0, 0, w, h)
        set_source_rgb(ctx, 0.2, 0.2, 0.2)
        fill(ctx)

        rendererVisitor(canvas, c)
    end
end

function rshow(canvas::RCanvas)
    c = @GtkCanvas()
    win = GtkWindow(c, "Roassal")
    redraw(canvas, c)

    signal_connect(win, "key-press-event") do widget, event
        println("You pressed key ", event.keyval)
    end

    c.mouse.motion = @guarded (widget, event) -> begin
        offset = offsetFromScreenToCanvas(c)
        shape_or_canvas_under_mouse = get_shape_at_position(canvas, event.x + offset[1], event.y + offset[2])
        #println("($(event.x), $(event.y)) -> $(typeof(shape_or_canvas_under_mouse))")
        trigger_callback(shape_or_canvas_under_mouse, :mouseMove, event)

        # Emit the enter / leave event
        if (canvas.shapeBeingPointed !== shape_or_canvas_under_mouse)
            #println("CHANGE SHAPE/CANVAS\n")
            trigger_callback(canvas.shapeBeingPointed, :mouseLeave, event)
            canvas.shapeBeingPointed = shape_or_canvas_under_mouse
            trigger_callback(shape_or_canvas_under_mouse, :mouseEnter, event)
        end

        #Probably trigger_callback should indicates whether there has been some trigger.
        redraw(canvas, c)
        reveal(widget)
        #println("refresh!!")
    end

    c.mouse.button1press = @guarded (widget, event) -> begin
        # ctx = getgc(widget)
        # set_source_rgb(ctx, 0, 1, 0)
        # arc(ctx, event.x, event.y, 5, 0, 2pi)
        # stroke(ctx)
        #reveal(widget)
        offset = offsetFromScreenToCanvas(c)
        shape_or_canvas_under_mouse = get_shape_at_position(canvas, event.x + offset[1], event.y + offset[2])
        trigger_callback(shape_or_canvas_under_mouse, :mouseClick, event)

        #Probably trigger_callback should indicates whether there has been some trigger.
        redraw(canvas, c)
        reveal(widget)
        #println("refresh!!")
    end

    show(c)
end

function get_shape_at_position(canvas::RCanvas, x::Number, y::Number)
    for shape in canvas.shapes
        c = compute_encompassing_rectangle(shape)
        if c[1] <= x && c[2] <= y && (c[3] + c[1]) > x && (c[4] + c[2]) > y
            return shape
        end
    end
    return canvas
end

function remove_shape!(canvas::RCanvas, shape::Shape)
    deleteat!(canvas.shapes, findall(s -> s == shape, canvas.shapes))
end
# ------------------------------------
"""
Rendering using a visitor
"""
function rendererVisitor(canvas::RCanvas, gtk::GtkCanvas=GtkCanvas())
    for shape in canvas.shapes
        rendererVisitor(shape, gtk)
    end
end

function rendererVisitor(box::RBox, gtk::GtkCanvas=GtkCanvas())
    ctx = getgc(gtk)
    encompassingRectangle = compute_encompassing_rectangle(box)
    _offsetFromCameraToScreen = offset_from_canvas_to_screen(gtk)
    rectangle(ctx,
                encompassingRectangle[1] + _offsetFromCameraToScreen[1],
                encompassingRectangle[2] + _offsetFromCameraToScreen[2],
                encompassingRectangle[3],
                encompassingRectangle[4])
    color = box.color
    set_source_rgb(ctx, color.r, color.g, color.b)
    fill(ctx)
end

function rendererVisitor(line::RLine, gtk::GtkCanvas=GtkCanvas())
    ctx = getgc(gtk)

    color = line.color
    set_source_rgb(ctx, color.r, color.g, color.b)

    _offsetFromCameraToScreen = offset_from_canvas_to_screen(gtk)
    from_position = pos(line.from) .+ _offsetFromCameraToScreen
    to_position = pos(line.to) .+ _offsetFromCameraToScreen
    move_to(ctx, from_position...)
    line_to(ctx, to_position...)
    set_line_width(ctx, 2.0)
    stroke(ctx)

    println("DEBUG: $color $from_position $to_position")

#=     fill(ctx)
    paint(cr)
 =##=
    rectangle(ctx,
                encompassingRectangle[1] + _offsetFromCameraToScreen[1],
                encompassingRectangle[2] + _offsetFromCameraToScreen[2],
                encompassingRectangle[3],
                encompassingRectangle[4])
    color = box.color
    set_source_rgb(ctx, color.r, color.g, color.b)
    fill(ctx) =#
end

function offset_from_canvas_to_screen(gtk::GtkCanvas)
    return (width(gtk) / 2, height(gtk) / 2)
end

function offsetFromScreenToCanvas(gtk::GtkCanvas)
    return offset_from_canvas_to_screen(gtk) .* -1
end

# ------------------------------------
"""
Callbacks
"""
mutable struct Callback
    name::Symbol
    f
end

function add_callback!(shape_or_canvas_under_mouse, callback::Callback)
    push!(shape_or_canvas_under_mouse.callbacks, callback)
end

function add_callback!(anArray::Array{Shape}, callback::Callback)
    for aShape in anArray
        add_callback!(aShape, callback::Callback)
    end
end

function numberOfCallbacks(shape_or_canvas_under_mouse)
    return length(shape_or_canvas_under_mouse.callbacks)
end

function trigger_callback(shape_or_canvas_under_mouse, name::Symbol, event=nothing)
    # Can be better written
    for c in shape_or_canvas_under_mouse.callbacks
        if (c.name == name)
            #c.f(event, shape_or_canvas_under_mouse)
            c.f()
        end
    end
end

# ------------------------------------
# Interactions
function highlightable(canvas::RCanvas)
    foreach(shape -> highlightable(shape), get_shapes(canvas))
end

function highlightable(shapes::Vector{Shape})
    foreach(shape -> highlightable(shape), shapes)
end
#=
function highlightable(shape::Shape)
    oldColor = nothing
    function recordOld()
        oldColor = get_color(s)
        set_color!(shape, RColor_BLUE)
        print("Recorded!")
        println(oldColor)
    end
    function giveOldColor()
        println("mouse leave!")
        set_color!(shape, oldColor)
    end

    #add_callback!(shape, Callback(:mouseEnter, (event, s) -> recordOld(s)))
    #add_callback!(shape, Callback(:mouseLeave, (event, s) -> giveOldColor(s)))
    add_callback!(shape, Callback(:mouseEnter, recordOld))
    add_callback!(shape, Callback(:mouseLeave, giveOldColor))
    return shape
end =#

highlighted_shapes = Dict{Shape,RColor}()
function highlightable(shape::Shape)
    function recordColor()
        global highlighted_shapes[shape] = get_color(shape)
        #print("Recorded! $shape \n")
        set_color!(shape, RColor_BLUE)
    end
    function restoreColor()
        #println("mouse leave! $(haskey(highlighted_shapes, shape))\n")
        if(haskey(highlighted_shapes, shape))
            set_color!(shape, highlighted_shapes[shape])
            delete!(highlighted_shapes, shape)
        end
    end

    add_callback!(shape, Callback(:mouseEnter, recordColor))
    add_callback!(shape, Callback(:mouseLeave, restoreColor))
    return shape
end

function reset_highlight()
    highlighted_shapes = Dict{Shape,RColor}()
end

# ------------------------------------
last_popup = nothing
function popup(shape::Shape)
    function addPopup()
        removePopup()
        global last_popup = RBox(width=40, height=15)
        add!(shape.canvas, last_popup)
        translate_to!(last_popup, pos(shape) .- (20, 20))
    end
    function removePopup()
        if(!isnothing(last_popup))
            remove_shape!(shape.canvas, last_popup)
            global last_popup = nothing
        end
    end

    add_callback!(shape, Callback(:mouseEnter, addPopup))
    add_callback!(shape, Callback(:mouseLeave, removePopup))
    return shape
end
# ------------------------------------

end


text:
"""
    gauss_jordan(A::AbstractMatrix{T}) where T<:Number

Gaussian elimination, also known as row reduction, is an algorithm for solving systems of linear equations.
It consists of a sequence of operations performed on the corresponding matrix of coefficients.
This method can also be used to compute the rank of a matrix, the determinant of a square matrix, and the inverse of an invertible matrix.
https://en.wikipedia.org/wiki/Gaussian_elimination

# Examples/Tests
```julia
julia> M1 = [1 2 3; 4 5 6];
julia> M2 = [1 2 3; 4 8 12];

julia> @test gauss_jordan(M1) == [1 0 -1; 0 1 2]        # Test Passed
julia> @test_throws AssertionError gauss_jordan(M2)     # Test Passed - Thrown: AssertionError
```

# Contributed By:- [AugustoCL](https://github.com/AugustoCL)
"""
function gauss_jordan(A::AbstractMatrix{T}) where {T<:Number}

    # check if matrix is singular
    m, n = size(A)
    if m == n
        @assert det(A) ≠ 0.0 "Must insert a non-singular matrix"
    else
        @assert det(A[:, 1:end-1]) ≠ 0.0 "Must insert a non-singular matrix or a system matrix [A b]"
    end

    # execute the gauss-jordan elimination
    for i ∈ axes(A, 1)

        if A[i, i] == 0.0
            for n ∈ (i + 1):m                           # iterate in lines below to check if could be swap
                if A[n, i] ≠ 0.0                        # condition to swap row
                    L = copy(A[i, :])                   # copy line to swap
                    A[i, :] = A[n, :]                   # swap occur
                    A[n, :] = L
                    break
                end
            end
        end

        @. A[i, :] = A[i, :] ./ A[i, i]                  # divide pivot line by pivot element

        for j ∈ axes(A, 1)                               # iterate each line for each pivot column, except pivot line
            if j ≠ i                                     # jump pivot line
                @. A[j, :] = A[j, :] - A[i, :]*A[j, i]   # apply gauss jordan in each line
            end
        end
    end

    return A
end

# using multiple dispatch to avoid InexactError with Integers
gauss_jordan(A::AbstractMatrix{T}) where {T<:Integer} = gauss_jordan(convert(Matrix{Float64}, A))

max_stars_repo_name: frankschmitt/Julia

max_issues_repo_name: frankschmitt/Julia

max_forks_repo_name: frankschmitt/Julia


text:
function CatValMul(A::Assoc,B::Assoc)
    if isa(Val(A)[1],String) && isa(Val(B)[1],String)
        AB = sortedintersect(A.col,B.row)
        A = A[:,AB]
        B = B[AB,:]
        rrr,ccc,vvv = findnz(Adj(A*B))
        ABVal = Array(Union{AbstractString,Number},length(rrr))
        for i in 1:length(rrr)
            r = rrr[i]
            c = ccc[i]
            ABIntersect = sortedintersect(Col(A[r,:]),Row(B[:,c]))
            ABValList = Array{Union{AbstractString,Number},1}()
            print(ABIntersect)
            for x in ABIntersect
                print(x)
                push!(ABValList,Val(A[r,x*","])[1])
                push!(ABValList,Val(B[x*",",c])[1])
            end
            if length(ABValList) > 0
                val = join(ABValList,";")*";"
                ABVal[i] = val
            end
        end
        return Assoc(Row(A)[rrr],Col(B)[ccc],ABVal)
    else
        return A*B
    end

end


text:
include("runtests-prerequisite.jl")

module Test_util_get_ids_cols_names

    using PostgresORM

    mutable struct MySuperMovie <: PostgresORM.IEntity

      title::Union{Missing,String}
      year::Union{Missing,Int32}

      MySuperMovie(args::NamedTuple) = MySuperMovie(;args...)
      MySuperMovie(;title = missing,
               year = missing,) = (
                x = new(missing, missing);
                x.title = title;
                x.year = year;
                return x
                )
    end

    mutable struct MyNewStruct <: PostgresORM.IEntity

      actor_id::Union{Missing,Int32}
      first_name::Union{Missing,String}
      favorite_movie::Union{Missing, MySuperMovie}

      MyNewStruct(args::NamedTuple) = MyNewStruct(;args...)
      MyNewStruct(;
            actor_id = missing,
            first_name = missing,
            favorite_movie = missing,
     ) = (
                x = new(missing, missing, missing);
                x.actor_id = actor_id;
                x.first_name = first_name;
                x.favorite_movie = favorite_movie;
                return x
                )

    end
end # module

# ############################################### #
# TEST 1: the ORM has its function 'get_id_props()' #
# ############################################### #
module MyNewStructORM
    data_type = Main.Test_util_get_ids_cols_names.MyNewStruct
    Main.PostgresORM.get_orm(x::Main.Test_util_get_ids_cols_names.MyNewStruct) = return(MyNewStructORM)
    get_table_name() = "public.my_new_struct"
    const columns_selection_and_mapping =
      Dict(
         :actor_id => "actor_id",
         :first_name => "first_name",
         :favorite_movie => ["favorite_movie_id"],
         :favorite_movies => ["favorite_movie1_id","favorite_movie2_id"]
      )
    get_id_props() = return [:actor_id,:first_name,:favorite_movie, :favorite_movies]

    # A dictionnary of mapping between fields symbols and overriding types
    #   Left hanside is the field symbol ; right hand side is the type override
    const types_override = Dict()

    const track_changes = false
end

@testset "Test `util_get_ids_cols_names(o::IEntity)` with a ORM with a function `get_id_props()``" begin
    myNewStruct = Test_util_get_ids_cols_names.MyNewStruct(actor_id = 3,first_name="toto")
    PostgresORM.get_orm(myNewStruct)
    @test PostgresORM.Controller.util_get_ids_cols_names(myNewStruct) ==
      ["actor_id", "first_name", "favorite_movie_id", "favorite_movie1_id", "favorite_movie2_id"]
end


text:
function my_lin_interp(grid, vals)
    function func(x::Real)
        a=searchsortedfirst(grid, x)
        b=searchsortedlast(grid, x)
        if a==length(grid)+1||b==0
            return 0
        elseif a==1
            return vals[1]
        else
            return (vals[a]-vals[a-1])/(grid[a]-grid[a-1])*(x-grid[a-1])+vals[a-1]
        end
    end

    function func{T<:Real}(x::AbstractVector{T})
        n = length(x)
        out = Array(Float64, n)
        for i in 1:n
            out[i] = func(x[i])
        end
        return out
    end

    return func
end


text:
"""
    cd!(rbm, data)

Trains the RBM on data using contrastive divergence.
"""
function cd!(rbm::RBM, data::AbstractArray;
    batchsize = 1,
    epochs = 1,
    optim = ADAM(), # optimizer algorithm
    wts = nothing, # weighted data points; named wts to avoid conflicts with RBM nomenclature
    steps::Int = 1, # Monte Carlo steps to update fantasy particles
    stats = suffstats(rbm, data; wts),
    callback = empty_callback
)
    @assert size(data) == (size(visible(rbm))..., size(data)[end])
    @assert isnothing(wts) || _nobs(data) == _nobs(wts)
    for epoch in 1:epochs
        batches = minibatches(data, wts; batchsize = batchsize)
        for (batch_idx, (vd, wd)) in enumerate(batches)
            vm = sample_v_from_v(rbm, vd; steps = steps)
            ∂ = ∂contrastive_divergence(rbm, vd, vm; wd = wd, wm = wd, stats)
            update!(∂, rbm, optim)
            update!(rbm, ∂)
            callback(; rbm, optim, epoch, batch_idx, vm, vd, wd)
        end
    end
    return rbm
end

"""
    contrastive_divergence(rbm, vd, vm; wd = 1, wm = 1)

Contrastive divergence loss.
`vd` is a data sample, and `vm` are samples from the model.
"""
function contrastive_divergence(
    rbm::RBM, vd::AbstractArray, vm::AbstractArray; wd = nothing, wm = nothing
)
    Fd = mean_free_energy(rbm, vd; wts = wd)
    Fm = mean_free_energy(rbm, vm; wts = wm)
    return Fd - Fm
end

function mean_free_energy(rbm::RBM, v::AbstractArray; wts = nothing)::Number
    # TODO: Allow passing sufficient statistics here, so that AD can exploit them too
    F = free_energy(rbm, v)
    if ndims(visible(rbm)) == ndims(v)
        wts::Nothing
        return F
    else
        return wmean(F; wts)
    end
end

function ∂contrastive_divergence(
    rbm::RBM, vd::AbstractArray, vm::AbstractArray;
    wd = nothing, wm = nothing,
    stats = suffstats(rbm, vd; wts = wd)
)
    ∂d = ∂free_energy(rbm, vd; wts = wd, stats)
    ∂m = ∂free_energy(rbm, vm; wts = wm)
    return subtract_gradients(∂d, ∂m)
end


text:
include("spharm.jl")

function ylmr2( lmax::Int, Gvec::Array{Float64,2} )

  Ng = size(Gvec)[2]
  lmax2 = (lmax+1)^2

  SMALL = 1e-9

  cost = zeros(Float64, Ng)
  sint = zeros(Float64, Ng)
  phi  = zeros(Float64, Ng)

  for ig = 1:Ng
    gmod = sqrt( Gvec[1,ig]^2 + Gvec[2,ig]^2 + Gvec[3,ig]^2 )
    if gmod < SMALL
       cost[ig] = 0.0
    else
       cost[ig] = Gvec[3,ig]/gmod
    end
    #  beware the arc tan, it is defined modulo pi
    if Gvec[1,ig] > SMALL
      phi[ig] = atan( Gvec[2,ig]/Gvec[1,ig] )
    elseif Gvec[1,ig] < -SMALL
      phi[ig] = atan( Gvec[2,ig]/Gvec[1,ig] ) + pi
    else
      phi[ig] = if Gvec[2,ig] >= 0 pi/2 else -pi/2 end
    end
    sint[ig] = sqrt( max( 0.0, 1-cost[ig]^2 ) )
  end

  ylm = zeros( Float64, Ng, lmax2 )

  ilm = 0
  for l = 0:lmax
    # m = 0
    ilm = ilm + 1
    ylm[:,ilm] = spharm( l, 0, Gvec, cost, sint, phi )
    # m > 0 and m < 0
    for m = 1:l
      ilm = ilm + 1
      ylm[:,ilm] = spharm( l, m, Gvec, cost, sint, phi )*(-1.0)^m
      #
      ilm = ilm + 1
      ylm[:,ilm] = spharm( l, -m, Gvec, cost, sint, phi )*(-1.0)^m
    end
  end

  return ylm

end


text:
include("bgzip_file.jl")
include("readbigvcf.jl")
include("vcf2dataframe.jl")
include("get_pedfile_trios.jl")
include("str_allele_length.jl")

"""
    (t, dt, p, nd, mutCon) = STR_iht(vcffile.vcf.gz, pedigree, STR_location, iid; showhaps=false)

Read a small **phased** vcf data array and return dataframes containing 1)
transmitting chromosomes and the STR allele differences and 2) a dataframe
of all phased individuals.

Individuals in the pedigree must match those in the vcf file.  The pedigree
file may contain multiple pedigrees.

STR location format is 1:258132133

Options:
    showhaps        show haplotypes for offspring with de novo STR [false]

"""
function STR_iht(vcf, pedfile, str, iid, mutSTR, othpar; showhaps=false)

    vcf_data_array = readbigvcf(vcf)

    hap_df = vcf2dataframe(vcf_data_array)

    (trios, peds) = get_pedfile_trios(pedfile)

    ftdf = DataFrame(ped=Any[], chr=Any[], pos=Any[], parent_chr=Any[],
    parent_sex=Any[], parent_str=Any[], parent_bp=Any[], child_chr=Any[],
    child_str=Any[], child_bp=Any[], str_diff=Any[], str_bp_diff=Any[] )

    fftdf = DataFrame(ped=Any[], chr=Any[], pos=Any[], parent_chr=Any[],
    parent_sex=Any[], parent_str=Any[], parent_bp=Any[], child_chr=Any[],
    child_str=Any[], child_bp=Any[], str_diff=Any[], str_bp_diff=Any[] )

    strchr, strpos = split(str, ':')
    stridx = argmax((hap_df[:, :POS] .== strpos))    #idx of str site

    str_sizes = str_allele_length(vcf, str)

for i in peds

    pid = i; header = []; tp = []; t_all=[]; tdf=DataFrame();

    t_all = get(trios, i, 0)

    for i in t_all
        line = "-------------------------\n"

        ids = split(i, ',')

        for a in eachindex(ids)
            push!(header, Symbol(ids[a] * "p"), Symbol(ids[a] * "q"))
        end

        tdf = hap_df[ :, header]  #df of each child,father,mother trio

        c1f1 = sum(tdf[!,1] .== tdf[!,3]) / size(tdf, 1) #fraction matching
        c1f2 = sum(tdf[!,1] .== tdf[!,4]) / size(tdf, 1) #across all sites
        c1m1 = sum(tdf[!,1] .== tdf[!,5]) / size(tdf, 1)
        c1m2 = sum(tdf[!,1] .== tdf[!,6]) / size(tdf, 1)
        c2f1 = sum(tdf[!,2] .== tdf[!,3]) / size(tdf, 1)
        c2f2 = sum(tdf[!,2] .== tdf[!,4]) / size(tdf, 1)
        c2m1 = sum(tdf[!,2] .== tdf[!,5]) / size(tdf, 1)
        c2m2 = sum(tdf[!,2] .== tdf[!,6]) / size(tdf, 1)

        c1altidx = findall(tdf[!,1] .!= "0")  #alt alleles matching
        c2altidx = findall(tdf[!,2] .!= "0")  #across sites, >1 allele
        f1altidx = findall(tdf[!,3] .!= "0")  #allowed
        f2altidx = findall(tdf[!,4] .!= "0")
        m1altidx = findall(tdf[!,5] .!= "0")
        m2altidx = findall(tdf[!,6] .!= "0")

        c1f1alt = length(intersect( Set(c1altidx), Set(f1altidx) ) ) #alt allele
        c1f2alt = length(intersect( Set(c1altidx), Set(f2altidx) ) ) #matched
        c1m1alt = length(intersect( Set(c1altidx), Set(m1altidx) ) ) #counts
        c1m2alt = length(intersect( Set(c1altidx), Set(m2altidx) ) )
        c2f1alt = length(intersect( Set(c2altidx), Set(f1altidx) ) )
        c2f2alt = length(intersect( Set(c2altidx), Set(f2altidx) ) )
        c2m1alt = length(intersect( Set(c2altidx), Set(m1altidx) ) )
        c2m2alt = length(intersect( Set(c2altidx), Set(m2altidx) ) )

        c1pfm =  [c1f1, c1f2, c1m1, c1m2]  #fractions for data above
        c2pfm =  [c2f1, c2f2, c2m1, c2m2]

        c1altmatch = [c1f1alt, c1f2alt, c1m1alt, c1m2alt] ./ size(tdf, 1)
        c2altmatch = [c2f1alt, c2f2alt, c2m1alt, c2m2alt] ./ size(tdf, 1)

        chriht = ["f1", "f2", "m1", "m2"]

        #pick all max indices for 1) total fraction of sites matched
        #and total alt alleles matched
        c1maxidx = [i for (i, x) in enumerate(c1pfm) if x == maximum(c1pfm)]
        c2maxidx = [i for (i, x) in enumerate(c2pfm) if x == maximum(c2pfm)]
        c1altmaxidx = [i for (i, x) in enumerate(c1altmatch) if x == maximum(c1altmatch)]
        c2altmaxidx = [i for (i, x) in enumerate(c2altmatch) if x == maximum(c2altmatch)]

        c1max = chriht[c1maxidx]   #translated idx to f1, f2, ...
        c2max = chriht[c2maxidx]

        c1altmax = chriht[c1altmaxidx]
        c2altmax = chriht[c2altmaxidx]

        if maximum(c1pfm) < 0.80     #declare phasing/rec. error for >20% mismatch
            c1max = ["pe"]; c1altmax = ["pe"]
        end

        if maximum(c2pfm) < 0.80
            c2max = ["pe"]; c2altmax = ["pe"]
        end

        c1_i = intersect(c1max, c1altmax)  #take isect of best frac & best alts
        c2_i = intersect(c2max, c2altmax)  #->most likely inheritance

        function pickone(c)  #if f1,f2 or m1,m2 equal probs then choose randomly
            if length(c) == 2 && ( c[1][1] == c[2][1] )
                nc = rand(c, 1)
                return nc
            else
                return c
            end
        end

        c1_i = pickone(c1_i)
        c2_i = pickone(c2_i)

        if length(c1_i) == 1  #logic to test other parent if one parent works
            if sum(occursin.(r"pe", c1_i)) > 0
                c1_i = ["nd"]
            elseif sum(occursin.(r"f", c1_i)) > 0
                c2_i = c2_i[occursin.(r"m", c2_i), :]
            elseif sum(occursin.(r"m", c1_i)) > 0
                c2_i = c2_i[occursin.(r"f", c2_i), :]
            end
        end

        if length(c2_i) == 1
            if sum(occursin.(r"pe", c2_i)) > 0
                c2_i = ["nd"]
            elseif sum(occursin.(r"f", c2_i)) > 0
                c1_i = c1_i[occursin.(r"m", c1_i), :]
            elseif sum(occursin.(r"m", c2_i)) > 0
                c1_i = c1_i[occursin.(r"f", c1_i), :]
            end
        end

        if ( (length(c1_i) == 0) || (length(c1_i) ≥ 2) )
            c1_i = []
            push!(c1_i, "nd")
        end

        if ( (length(c2_i) == 0) || (length(c2_i) ≥ 2) )
            c2_i = []
            push!(c2_i, "nd")
        end

        h = String.(header)

        function mapid(c, h)  #translate back to ids
            k = replace.(c, "f1" => h[3])
            l = replace.(k, "f2" => h[4])
            m = replace.(l, "m1" => h[5])
            n = replace.(m, "m2" => h[6])
            return n
        end

        c1_bestids = mapid(c1_i, h)
        c2_bestids = mapid(c2_i, h)

        #flatten arrays and push results to a df
        c1_i = join(reshape(c1_i, (1, length(c1_i) ) ), ",")
        c2_i = join(reshape(c2_i, (1, length(c2_i) ) ), ",")
        c1_bestids = join(reshape(c1_bestids, (1, length(c1_bestids) ) ), ",")
        c2_bestids = join(reshape(c2_bestids, (1, length(c2_bestids) ) ), ",")

        #SET TO TRUE TO SEE THE HAPLOTYPES FOR PARENT AND GRANDPARENTS!
        if showhaps == true
            if occursin(Regex("$iid[p|q]"), names(tdf)[1]) #print only denovo parent trio
                println("TABLE: phased chromosomes for parent with de novo STR and both grandparents.")
                println("TABLE INFO: genotypes for the (de novo) STR locus are shown as -1 (unphasable).")
                println(tdf)
                println("$iid p transmitted from: ", c1_bestids,
                "\n", "$iid q transmitted from: ", c2_bestids, "\n")
            end
        end

        c1chr = h[1]
        c2chr = h[2]

        if (c1_bestids != "nd")
            p1_allele = parse(Int64, (tdf[ stridx, "$c1_bestids" ] ) )
            c1_allele = parse(Int64, (tdf[ stridx, "$c1chr" ] ) )
            p1_length = str_sizes[p1_allele]
            c1_length = str_sizes[c1_allele]
            c1strdiff = abs(p1_allele - c1_allele)
            c1bpdiff  = p1_length - c1_length
        else
            c1strdiff = "nd"; c1_allele = "nd"; c1_length = "nd"
            c1bpdiff  = "nd"; p1_allele = "nd"; p1_length = "nd"
        end

        if (c2_bestids != "nd")
            p2_allele = parse(Int64, (tdf[ stridx, "$c2_bestids" ] ) )
            c2_allele = parse(Int64, (tdf[ stridx, "$c2chr" ] ) )
            p2_length = str_sizes[p2_allele]
            c2_length = str_sizes[c2_allele]
            c2strdiff = abs(p2_allele - c2_allele)
            c2bpdiff  = p2_length - c2_length
        else
            c2strdiff = "nd"; c2_allele = "nd"; c2_length = "nd"
            c2bpdiff  = "nd"; p2_allele = "nd"; p2_length = "nd"
        end

        if occursin(r"f", c1_i)
            p1_sex = "male"
        elseif occursin(r"m", c1_i)
            p1_sex = "female"
        else
            p1_sex = "nd"
        end

        if occursin(r"f", c2_i)
            p2_sex = "male"
        elseif occursin(r"m", c2_i)
            p2_sex = "female"
        else
            p2_sex = "nd"
        end

        push!(ftdf, [pid strchr strpos c1_bestids p1_sex p1_allele p1_length c1chr c1_allele c1_length c1strdiff c1bpdiff] )
        push!(ftdf, [pid strchr strpos c2_bestids p2_sex p2_allele p2_length c2chr c2_allele c2_length c2strdiff c2bpdiff] )

        r1 = []   #reinit for next trio
        r2 = []
        header = []
    end

    fracnd = round(sum(ftdf.str_diff .== "nd") / size(ftdf, 1), digits=2) #frac bad

    fftdf = ftdf[ .&(ftdf.str_diff .!= 0, ftdf.str_diff .!= "nd"), :]

    #Next, extract the rest of the family that carries the mutant STR allele
    #and get the consensus haplotype with the mutation and compare
    #to the grandparents later.

    mutCon = nothing
    num_mut_iht = 0

    if size(hap_df,2) > 7  #run fam but not trio here

        #remove othpar if present
        othpar1 = join([othpar, "p"], "")
        othpar2 = join([othpar, "q"], "")
        ftdf = ftdf[findall(!in([ othpar1, othpar2 ]), ftdf.child_chr), : ]

        mut_offspr = (ftdf[ ftdf.child_str .== mutSTR, [:child_chr]]) #df
        mut_ids = vec(Array(mut_offspr))  #array, 2209p ...
        mut_df = hap_df[:, Symbol.(mut_ids)]

        if sum(occursin.(Regex(othpar), mut_ids)) > 0 #other parent has de novo
            println("WARN: one or more phased chromosomes from parent $othpar
cannot be distiguished from the chromosome with the de novo STR allele.")
            println("INFO: omitting phased chromosomes containing the de novo STR
allele that match the other parent $othpar at a rate of ≥ 0.85")
            othchrid = mut_offspr[ occursin.(Regex(othpar), mut_ids), :]
            othchrid = othchrid[1,1]
            othchr = Array(mut_df[:, Symbol(othchrid)])

            p_oth = Float64[]
            not_othchr = Int64[]

            if showhaps == true
                println("TABLE: Phased offspring chromosomes with de novo allele:")
                println(mut_df)
            end

            for i in 1:size(mut_df, 2)
                p = sum(othchr .== mut_df[:,i]) / size(mut_df, 1)
                push!(p_oth, p)
            end

            for i in 1:length(p_oth)
                if p_oth[i] ≥ 0.85
                    continue
                else
                    push!(not_othchr, i)
                end
            end

            mut_df = mut_df[:, not_othchr]

        end

        if size(mut_df, 2) == 0

            println("WARN: The offsprings' consensus phased chromosome with the
de novo STR is very similar the other parent. This result can be
caused by 1) high similarity among grandparent haplotypes,
2) incorrect phasing, or 3) a false positive de novo STR.
Try different haplotype lengths.\n")

            num_mut_iht = -9
            return(fftdf, ftdf, hap_df, fracnd, mutCon, str_sizes, num_mut_iht)
        end

        av = sum(eachcol(mut_df .== "1")) / size(mut_df,2) #generate a consensus
        mutCon = Array{String,1}(undef, length(av) )       #haplotype (mutCon)

        for i in eachindex(av)
            if av[i] ≥ 0.5
                mutCon[i] = "1"
            else
                mutCon[i] = "0"
            end
        end

        num_mut_iht = size(mut_df,2)

        println("INFO: Found $num_mut_iht offspring carrying the de novo allele $mutSTR\n")

    end

    sort!(fftdf, :parent_chr)

    return (fftdf, ftdf, hap_df, fracnd, mutCon, str_sizes, num_mut_iht)

end
end

max_stars_repo_name: ScottWatkins/STRdiff

max_issues_repo_name: ScottWatkins/STRdiff

max_forks_repo_name: ScottWatkins/STRdiff


text:
# compilation cache

using Core.Compiler: retrieve_code_info, CodeInfo, MethodInstance, SSAValue, SlotNumber, ReturnNode
using Base: _methods_by_ftype

using Serialization, Scratch

const cache_lock = ReentrantLock()

# whether to cache compiled kernels on disk or not
const disk_cache = Ref(false)

@inline function check_cache(cache, job, prekey,
                             @nospecialize(compiler), @nospecialize(linker))
    key = hash(job, prekey)
    force_compilation = compile_hook[] !== nothing

    # XXX: by taking the hash, we index the compilation cache directly with the world age.
    #      that's wrong; we should perform an intersection with the entry its bounds.

    # NOTE: no use of lock(::Function)/@lock/get! to keep stack traces clean
    lock(cache_lock)
    try
        obj = get(cache, key, nothing)
        if obj === nothing || force_compilation
            asm = nothing

            # can we load from the disk cache?
            if disk_cache[] && !force_compilation
                path = joinpath(@get_scratch!("kernels"), "$key.jls")
                if isfile(path)
                    try
                        asm = deserialize(path)
                        @debug "Loading compiled kernel for $spec from $path"
                    catch ex
                        @warn "Failed to load compiled kernel at $path" exception=(ex, catch_backtrace())
                    end
                end
            end

            # compile
            if asm === nothing
                if compile_hook[] !== nothing
                    compile_hook[](job)
                end

                asm = compiler(job)

                if disk_cache[] && !isfile(path)
                    serialize(path, asm)
                end
            end

            # link (but not if we got here because of forced compilation)
            if obj === nothing
                obj = linker(job, asm)
                cache[key] = obj
            end
        end
        obj
    finally
        unlock(cache_lock)
    end
end

# generated function that crafts a custom code info to call the actual compiler.
# this gives us the flexibility to insert manual back edges for automatic recompilation.
#
# we also increment a global specialization counter and pass it along to index the cache.

const specialization_counter = Ref{UInt}(0)

@generated function cached_compilation(cache::AbstractDict,
                                       job::CompilerJob{<:Any,<:Any,FunctionSpec{f,tt}},
                                       compiler::Function, linker::Function) where {f,tt}
    # get a hold of the method and code info of the kernel function
    sig = Tuple{f, tt.parameters...}
    # XXX: instead of typemax(UInt) we should use the world-age of the fspec
    mthds = _methods_by_ftype(sig, -1, typemax(UInt))
    Base.isdispatchtuple(tt) || return(:(error("$tt is not a dispatch tuple")))
    length(mthds) == 1 || return (:(throw(MethodError(job.source.f,job.source.tt))))
    mtypes, msp, m = mthds[1]
    mi = ccall(:jl_specializations_get_linfo, Ref{MethodInstance}, (Any, Any, Any), m, mtypes, msp)
    ci = retrieve_code_info(mi)
    @assert isa(ci, CodeInfo)

    # generate a unique id to represent this specialization
    id = (specialization_counter[] += 1)
    # TODO: save the mi/ci here (or embed it in the AST to pass to the compiler)
    #       and use that to drive compilation

    # prepare a new code info
    new_ci = copy(ci)
    empty!(new_ci.code)
    empty!(new_ci.codelocs)
    resize!(new_ci.linetable, 1)                # see note below
    empty!(new_ci.ssaflags)
    new_ci.ssavaluetypes = 0
    new_ci.edges = MethodInstance[mi]
    # XXX: setting this edge does not give us proper method invalidation, see
    #      JuliaLang/julia#34962 which demonstrates we also need to "call" the kernel.
    #      invoking `code_llvm` also does the necessary codegen, as does calling the
    #      underlying C methods -- which GPUCompiler does, so everything Just Works.

    # prepare the slots
    new_ci.slotnames = Symbol[Symbol("#self#"), :cache, :job, :compiler, :linker]
    new_ci.slotflags = UInt8[0x00 for i = 1:5]
    cache = SlotNumber(2)
    job = SlotNumber(3)
    compiler = SlotNumber(4)
    linker = SlotNumber(5)

    # call the compiler
    append!(new_ci.code, [Expr(:call, check_cache, cache, job, id, compiler, linker),
                          ReturnNode(SSAValue(1))])
    append!(new_ci.codelocs, [1, 1])   # see note below
    new_ci.ssavaluetypes += 2

    # NOTE: we keep the first entry of the original linetable, and use it for location info
    #       on the call to check_cache. we can't not have a codeloc (using 0 causes
    #       corruption of the back trace), and reusing the target function's info
    #       has as advantage that we see the name of the kernel in the backtraces.

    return new_ci
end


text:
# Must add CSV, ReadStat, DataFrames and Gadfly packages needed for this project
# import Pkg; Pkg.add("Pkg")
# Pkg.add("CSV")
# Pkg.add("ReadStat")
# Pkg.add("DataFrames")
# Pkg.add("Gadfly")
# Pkg.add("Compose")
# Pkg.add("StatFiles")
# Pkg.add("VegaLite")
# Pkg.add("VegaDatasets")
#import Pkg; Pkg.add("JSON")
#import Pkg; Pkg.add("Requests")
#Pkg.clone("https://github.com/plotly/Plotly.jl")
#import Pkg; Pkg.add("Plotly")

# Don't add NodeJS but know that it is needed.
# Pkg.add("NodeJS")
# Pkg.build("NodeJS")

# Pkg.add("VegaLite")
# Pkg.add("VegaDatasets")
# Pkg.add("DataVoyager")
# Pkg.add("IJulia")
# Pkg.build("VegaLite")

# Load required libraries
# using DataFrames, CSV, Gadfly, Compose, ReadStat, StatFiles, VegaLite, VegaDatasets;
using DataFrames, CSV, Gadfly, StatFiles, ReadStat

# Read tables using the updated CSV.jl package v0.5.3
   demo17q2 = CSV.read("C:\\Users\\churley\\data\\DEMO17Q2.txt"; header = 1, delim = '$')
   drug17q2 = CSV.read("C:\\Users\\churley\\data\\DRUG17Q2.txt"; header = 1, delim = '$')
   indi17q2 = CSV.read("C:\\Users\\churley\\data\\INDI17Q2.txt"; header = 1, delim = '$')
   outc17q2 = CSV.read("C:\\Users\\churley\\data\\OUTC17Q2.txt"; header = 1, delim = '$')
   reac17q2 = CSV.read("C:\\Users\\churley\\data\\REAC17Q2.txt"; header = 1, delim = '$')
   rpsr17q2 = CSV.read("C:\\Users\\churley\\data\\RPSR17Q2.txt"; header = 1, delim = '$')
   ther17q2 = CSV.read("C:\\Users\\churley\\data\\THER17Q2.txt"; header = 1, delim = '$')
#    ther17q2 = CSV.read("C:\\Users\\churley\\data\\THER17Q2.txt"; header = 1, delim = '$', datarow=195880)

# Read SAS data using a function from the StatFiles package
# using StatFiles, DataFrames
# Can't use MedDRA because of licensing
# meddra20 = load("C:\\Users\\churley\\data\\meddra20_0.sas7bdat") |> DataFrame
# m20=read_sas7bdat("C:\\Users\\churley\\data\\meddra20_0.sas7bdat")
# meddra20=DataFrame(m20.data,m20.headers)

# There are some columns that need to be renamed.
# We can envision a future state where this is automated... for now manual.
indi17q2 = rename(indi17q2, :indi_drug_seq => :drug_seq)
ther17q2 = rename(ther17q2, :dsg_drug_seq => :drug_seq);

# For proof-of-concept speed purposes, we will use only the first 10,000 rows of each dataset.
demo17q2 = first(demo17q2, 10000)
drug17q2 = first(drug17q2, 10000)
indi17q2 = first(indi17q2, 10000)
outc17q2 = first(outc17q2, 10000)
reac17q2 = first(reac17q2, 10000)
rpsr17q2 = first(rpsr17q2, 10000)
ther17q2 = first(ther17q2, 10000);

########## Perform joins
# Pull in indications for use into drug file table
drug_ind = join(
    drug17q2[:, [:primaryid, :caseid, :drug_seq, :drugname, :route, :cum_dose_chr, :dose_amt,
        :dose_unit, :dose_form, :dose_freq, :role_cod]],
    indi17q2[:, [:primaryid, :drug_seq, :indi_pt]],
    on = [:primaryid, :drug_seq], kind = :left
)

drug_ind_ther = join(
    drug_ind, ther17q2[:, [:primaryid, :drug_seq, :start_dt, :end_dt, :dur, :dur_cod]],
    on = [:primaryid, :drug_seq], kind = :left
)

# Pull reaction data into patient demographic table
demo_reac = join(
    demo17q2[:, [:primaryid, :caseid, :event_dt, :age, :age_cod, :age_grp, :sex, :wt, :occr_country]],
    reac17q2[:, [:primaryid, :pt, :drug_rec_act]],
    on = :primaryid, kind = :left,
)

# Pull outcomes information into the previous table
demo_reac_out = join(
    demo_reac, outc17q2[:, [:primaryid, :outc_cod]], on = :primaryid, kind = :left
)

# Concatenate previous table with report source information
demo_reac_out_source = join(
    demo_reac_out, rpsr17q2[:, [:primaryid, :rpsr_cod]], on = :primaryid, kind = :left
)

# Lastly, join the (Drug/Therapy/Indication) and (Demographic/Reaction/Outcome/Report Source)
# combined sets into one set with all this information in one place.
#faers17 = join(drug_ind_ther, delete!(demo_reac_out_source, :caseid), on = :primaryid, kind = :left);
faers17 = join(drug_ind_ther, deletecols!(demo_reac_out_source, :caseid), on = :primaryid, kind = :left);
########## End joins

########## Perform transformations
# Begin by storing a vector of counts by PT, where PT is "Preferred Term", or the event.
faers17q2_1 = by(faers17, :pt, d -> DataFrame(count_pt_17q2 = nrow(d)))

# Append this vector to the combined FAERS dataset.
faers17q2_2 = join(faers17, faers17q2_1, on = :pt, kind = :left, makeunique = true)

# Create a vector of counts by combination of Drug, and PT.
faers17q2_3 = by(faers17q2_2, [:drugname, :pt],
    d -> DataFrame(count_dpt_17q2 = nrow(d)))

#
faers17q2_4 = join(faers17q2_2, faers17q2_3, on = [:drugname, :pt], kind = :left,
    makeunique = true)

faers17q2_4[:count_percentage_17q2] =
    (faers17q2_4[:count_dpt_17q2] ./ faers17q2_4[:count_pt_17q2]) .* 100;
########## End transformations

#Show dataframe faers17q2_1

withenv("LINES" => 5) do
    display(faers17q2_1)
end

#Show dataframe faers17q2_1

withenv("LINES" => 5) do
    display(faers17q2_2)
end

#Show dataframe faers17q2_1

withenv("LINES" => 5) do
    display(faers17q2_3)
end

#Show dataframe faers17q2_1

withenv("LINES" => 5) do
    display(faers17q2_4)
end

# View final dataset.  ETL is complete.
first(faers17q2_4, 10)

# Exploring with the describe() function:
drug_case_pct = faers17q2_4[:count_percentage_17q2]

describe(drug_case_pct)

# Isolate the pilot drug
pilot_data = faers17q2_4[faers17q2_4[:drugname] .== "ATACAND", :]

# Adverse events by country
# Bar chart for now.  Consider the JuliaGeo library
d1 = unique(faers17q2_4[:, [:primaryid, :pt, :occr_country]])
d2 = by(d1, :occr_country,
    d -> DataFrame(Frequency = nrow(d)))
sort!(d2, :Frequency, rev = true)

Gadfly.plot(d2, x = :occr_country, y = :Frequency, Geom.bar,
    Guide.xlabel("Country"),
    Guide.ylabel("Frequency"),
    Guide.title("Adverse Events by Country"),
    Scale.y_continuous(format=:plain),
    Theme(default_color = "red"))

#VegaLite data presentation will use data from all drugs
d1 = unique(faers17q2_4[:, [:primaryid, :sex, :pt, :occr_country, :drugname]])
d2 = by(d1, [:occr_country, :sex, :drugname],
    d -> DataFrame(Frequency = nrow(d)))
sort!(d2, :Frequency, rev = true)
rename!(d2, :sex => :gender)
allc = first(d2,60)

#Julia visualization remove after figuring out how to use vegalite
using VegaLite, VegaDatasets

allc |> @vlplot(:bar,
        x={"occr_country:n", axis={title="Country"}},
        title="Adverse Events by Country and Gender",
        y={:Frequency, axis={title="Number of AEs"}},
        color={"gender:n", scale={range=["#1f77b4","#e377c2"]}},
        height=400,
        width=600)

allc |> @vlplot(:bar,

        x={"occr_country:n", axis={title="Country"}},
        title="Adverse Events by Country and Gender - Atacand Only",
        y={:Frequency, axis={title="Number of AEs"}},
        color={"gender:n", scale={range=["#1f77b4","#e377c2"]}},
        height=400,
        width=600)

Pkg.add("PlotlyJS")

using PlotlyJS

# Adverse events by gender
d2 = unique(faers17q2_4[:, [:primaryid, :pt, :sex]])
d2 = by(d2, :sex,
    d -> DataFrame(Frequency = nrow(d)))
d2[1, :sex] = "Unknown"

palette = ["grey", "blue", "pink"]
# labels = string(d2[:, :Frequency])
# labels

Gadfly.plot(d2, x = :sex, label = 1, y = :Frequency, color = :sex, Geom.bar,
    Geom.label(position = :above),
    Guide.xlabel("Gender"),
    Guide.ylabel("Frequency"),
    Guide.title("Adverse Events by Gender"),
    Scale.y_continuous(format=:plain),
    Scale.color_discrete_manual(palette...)
    #Guide.annotation(compose(context(), text(:sex, :Frequency, labels)))
    )

describe(d2)

# Adverse events by age group
d3 = unique(faers17q2_4[:, [:primaryid, :pt, :age_grp]])
d3 = by(d3, :age_grp,
    d -> DataFrame(Frequency = nrow(d)))
d3[1, :age_grp] = "Unknown"

Gadfly.plot(d3, x = :age_grp, y = :Frequency, color = :age_grp, Geom.bar,
    Guide.xlabel("Age Group"),
    Guide.ylabel("Frequency"),
    Guide.title("Adverse Events by Age Group"),
    Scale.y_continuous(format=:plain),
    Guide.colorkey("Legend")
    )

# Adverse events by role of the drug
d4 = unique(faers17q2_4[:, [:primaryid, :pt, :role_cod]])
d4 = by(d4, :role_cod,
    d -> DataFrame(Frequency = nrow(d)))

Gadfly.plot(d4, x = :role_cod, y = :Frequency, color = :role_cod, Geom.bar,
    Guide.xlabel("Role of the Drug"),
    Guide.ylabel("Frequency"),
    Guide.title("Adverse Events by Drug Role"),
    Scale.y_continuous(format=:plain),
    Guide.colorkey("Legend")
    )

d5 = unique(faers17q2_4[:, [:primaryid, :pt]])
d5 = by(d5, :pt,
    d -> DataFrame(Frequency = nrow(d)))
names!(d5, [:Adverse_Event, :AEs_in_2017_Q2])
#rename!(d5, :pt => :Adverse_Event)
sort!(d5, :AEs_in_2017_Q2, rev = true)

#VegaLite data presentation
head(d5,10)

#Julia visualization remove after figuring out how to use vegalite
using VegaLite, VegaDatasets

data = dataset("cars")

data |> @vlplot(:point,
        x= "Cylinders:o",
        title="Number of Cylinders",
        y=:Miles_per_Gallon,
        color="Origin:n",
        height=400,
        width=300)

#data |> @vlplot(:point) |> display

#using VegaLite

#d5 |> @vlplot(:point, x=:AEs_in_2017_Q2, y=:Adverse_Event)

Pkg.add("Query")

using VegaLite, VegaDatasets

dataset("gapminder-health-income") |>
@vlplot(
    :circle,
    width=500,height=300,
    selection={
        view={typ=:interval, bind=:scales}
    },
    y={:health, scale={zero=false}},
    x={:income, scale={typ=:log}},
    size=:population,
    color={value="#000"}
)

using VegaLite, VegaDatasets

us10m = dataset("us-10m").path
unemployment = dataset("unemployment.tsv").path

@vlplot(
    :geoshape,
    width=500, height=300,
    data={
        url=us10m,
        format={
            typ=:topojson,
            feature=:counties
        }
    },
    transform=[{
        lookup=:id,
        from={
            data=unemployment,
            key=:id,
            fields=["rate"]
        }
    }],
    projection={
        typ=:albersUsa
    },
    color="rate:q"
)

using VegaLite, DataFrames

x = [j for i in -5:4, j in -5:4]
y = [i for i in -5:4, j in -5:4]
z = x.^2 .+ y.^2
data = DataFrame(x=vec(x'),y=vec(y'),z=vec(z'))

data |> @vlplot(:rect, x="x:o", y="y:o", color=:z)

max_stars_repo_name: andyjessen/phuse-scripts

max_issues_repo_name: andyjessen/phuse-scripts

max_forks_repo_name: andyjessen/phuse-scripts


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using LinearAlgebra, Random
using Distributions, PDMats, StatsBase

@testset "distribution_functions" begin
    @testset "_check_rand_compat" begin
        @test BAT._check_rand_compat(MvNormal(ones(2)), ones(2,10)) == nothing
        @test_throws DimensionMismatch BAT._check_rand_compat(MvNormal(ones(3)), ones(2,10))
    end

    @testset "BAT.issymmetric_around_origin" begin
        @test BAT.issymmetric_around_origin(Normal(0.0, 3.2)) == true
        @test BAT.issymmetric_around_origin(Normal(1.7, 3.2)) == false
        @test BAT.issymmetric_around_origin(Gamma(4.2, 2.2)) == false
        @test BAT.issymmetric_around_origin(Chisq(20.3)) == false
        @test BAT.issymmetric_around_origin(MvNormal(zeros(2), ones(2))) == true
        @test BAT.issymmetric_around_origin(MvNormal(ones(2), ones(2))) == false
        @test BAT.issymmetric_around_origin(MvTDist(1.5, zeros(2), PDMat(Matrix{Float64}(I, 2, 2))))
        @test BAT.issymmetric_around_origin(MvTDist(1.5, ones(2), PDMat(Matrix{Float64}(I, 2, 2)))) == false
    end

    @testset "BAT.cov2pdmat" begin
        @test BAT.cov2pdmat(Float64, PDMat(Matrix(1.0I, 2, 2))) == PDMat(Matrix(1.0I, 2, 2))
        @test BAT.cov2pdmat(Float64, PDiagMat([1, 1])) == BAT.cov2pdmat(Float64, Matrix(PDiagMat([1, 1])))
        @test BAT.cov2pdmat(Float64, ScalMat(2, 3)) == BAT.cov2pdmat(Float64, Matrix(ScalMat(2, 3)))
        @test BAT.cov2pdmat(Float64, cholesky(Matrix(1.0I, 2, 2))) == PDMat(cholesky(Matrix(1.0I, 2, 2)))
        @test BAT.cov2pdmat(Float16, cholesky(Matrix(1.0I, 2, 2))) == BAT.cov2pdmat(Float16, Matrix(1.0I, 2, 2))
    end
end


text:
using BenchmarkTools

#%%

function test1(x::Int)
    println("$x is an Int")
end

function test1(x::Float64)
    println("$x is a Float64")
end

function test1(x)
    println("$x is neither an Int nor a Float64")
end
#%%
test1(1)

test1(1.0)

test1("techytok")

#%%
"""
Returns x if x>0 else returns 0
"""
function test3(x)
    if x>0
        return x
    else
        return 0
    end
end
#%%
test3(2)
test3(-1)

test3(2.0)
test3(-1.0)
#%%
"""
Returns x if x>0 else returns 0
"""
function test4(x)
    if x>0
        return x
    else
        return zero(x)
    end
end
#%%
test4(-1.0)

#%%
@code_warntype test4(1.0)

#%%
function test5()
    r=0
    for i in 1:10
        r+=sin(i)
    end
    return r
end

function test6()
    r=0.0
    for i in 1:10
        r+=sin(i)
    end
    return r
end
#%%
@code_warntype test5()
@code_warntype test6()

@btime test5()
@btime test6()

#%%
function test7(x)
    result = 0
    if x>0
        result = x
    else
        result = 0
    end
    return result::typeof(x)
end
#%%
test7(-2.0)

#%%
function take_a_breath()
    sleep(22*1e-3)
    return
end

function test8()
    r=zeros(100,100)
    take_a_breath()
    for i in 1:100
        A=rand(100,100)
        r+=A
    end
    return r
end

#%%
test8()
@profiler test8()
Juno.Profile.print()

#%%
using BenchmarkTools

arr1=zeros(10000)

function put1!(arr)
    for i in 1:length(arr)
        arr[i]=1.0
    end
end

function put1_inbounds!(arr)
    @inbounds for i in 1:length(arr)
        arr[i]=1.0
    end
end

@btime put1!($arr1)

@btime put1_inbounds!($arr1)

#%%

#%%
using BenchmarkTools
using LinearAlgebra
using StaticArrays

add!(C, A, B) = (C .= A .+ B)

function simple_bench(N, T=Float64)
    A = rand(T,N,N)
    A = A'*A
    B = copy(A)
    SA = SMatrix{N,N}(A)
    MA = MMatrix{N,N}(A)
    MB = copy(MA)

    print("""
============================================
    Benchmarks for $N×$N $T matrices
============================================
""")
    ops = [
           ("Matrix multiplication              ", *,        (A, A),     (SA, SA)),
           ("Matrix multiplication (mutating)   ", mul!,     (B, A, A),  (MB, MA, MA)),
           ("Matrix addition                    ", +,        (A, A),     (SA, SA)),
           ("Matrix addition (mutating)         ", add!,     (B, A, A),  (MB, MA, MA)),
           ("Matrix determinant                 ", det,      (A,),       (SA,)),
           ("Matrix inverse                     ", inv,      (A,),       (SA,)),
           ("Matrix symmetric eigendecomposition", eigen,    (A,),       (SA,)),
           ("Matrix Cholesky decomposition      ", cholesky, (A,),       (SA,)),
           ("Matrix LU decomposition            ", lu,       (A,),       (SA,)),
           ("Matrix QR decomposition            ", qr,       (A,),       (SA,)),
    ]
    for (name, op, Aargs, SAargs) in ops
        # We load from Ref's here to avoid the compiler completely removing the
        # benchmark in some cases.
        #
        # Like any microbenchmark, the speedups you see here should only be
        # taken as roughly indicative of the speedup you may see in real code.
        if length(Aargs) == 1
            A1  = Ref(Aargs[1])
            SA1 = Ref(SAargs[1])
            speedup = @belapsed($op($A1[])) / @belapsed($op($SA1[]))
        elseif length(Aargs) == 2
            A1  = Ref(Aargs[1])
            A2  = Ref(Aargs[2])
            SA1 = Ref(SAargs[1])
            SA2 = Ref(SAargs[2])
            speedup = @belapsed($op($A1[], $A2[])) / @belapsed($op($SA1[], $SA2[]))
        elseif length(Aargs) == 3
            A1  = Ref(Aargs[1])
            A2  = Ref(Aargs[2])
            A3  = Ref(Aargs[3])
            SA1 = Ref(SAargs[1])
            SA2 = Ref(SAargs[2])
            SA3 = Ref(SAargs[3])
            speedup = @belapsed($op($A1[], $A2[], $A3[])) / @belapsed($op($SA1[], $SA2[], $SA3[]))
        else
        end
        println(name*" -> $(round(speedup, digits=1))x speedup")
    end
end
#%%
simple_bench(3)

#%%
arr1 = zeros(100, 200)

@btime for i in 1:100
    for j in 1:200
        arr1[i,j] = 1
    end
end

@btime for j in 1:200
    for i in 1:100
        arr1[i,j] = 1
    end
end

#%%
function test_positional(a, b, c)
    return a + b + c
end

function test_keyword(;a, b, c)
    return a + b + c
end

@code_llvm test_positional(1,2,3)

@code_llvm test_keyword(a=1,b=2,c=3)

@btime test_positional(1,2,3)

@btime test_keyword(a=1,b=2,c=3)

max_stars_repo_name: aurelio-amerio/techytok-examples

max_issues_repo_name: aurelio-amerio/techytok-examples

max_forks_repo_name: aurelio-amerio/techytok-examples


text:
#=
Implement an autocomplete system. That is, given a query string s and a set of all possible query strings, return all strings in the set that have s as a prefix.

For example, given the query string de and the set of strings [dog, deer, deal], return [deer, deal].

Hint: Try preprocessing the dictionary into a more efficient data structure to speed up queries.
=#

function autocomplete(string::String, dictionary::Array{String,1})
    string_length = length(string)
    possible_query = dictionary[first.(dictionary, string_length) .== string]
    return possible_query
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
function aetherlang_show(io::IO, v::Vector)
    print(io, "(arr "*join([aetherlang_repr(x) for x in v], ' ')*")")
end

max_stars_repo_name: Aetherlang/aetherlang

max_issues_repo_name: gojakuch/aetherlang

max_forks_repo_name: gojakuch/aetherlang


text:
#=
surface_extraction:
- Julia version: 1.1.0
- Author: Jirik
- Date: 2019-03-19
=#

using SparseArrays
using Distributed
_b3_size = nothing
# ch = nothing
# using LarSurf

function set_block_size(block_size)
    global _b3_size
    _b3_size = size_as_array(block_size)
end

# function set_channel(channel)
#     if channel == nothing
#         ch = RemoteChannel(()->Channel{Int}(32));
#     else
#         ch = channel
#     end
# @everywhere begin
    const chnnel_big_fids = Channel{Int}(32);

    """
    Calculate multiplication linearized volume cells with boundary matrix and
    filter it for number 1.
    """
    function grid_get_surf_Fvec_larmodel(segmentation::AbstractArray)
        segClin = grid_to_linear(segmentation, 0)

        block_size = size_as_array(size(segmentation))

        b3, larmodel = get_boundary3(block_size)
        Flin = segClin' * b3
        # Matrix(Flin)
        sparse_filter!(Flin, 1, 1, 0)
        dropzeros!(Flin)
        return Flin, larmodel
    end

    function grid_get_surf_Fvec_larmodel_old(segmentation::AbstractArray)
        segClin = grid_to_linear(segmentation, 0)

        block_size = size_as_array(size(segmentation))

        b3, larmodel = get_boundary3(block_size)
        V, topology = larmodel
        (VV, EV, FV, CV) = topology
        # println("segmentation: ", size(segmentation))
        # println("segClin: ", size(segClin), " ", typeof(segClin))
        # println("b3: ", size(b3), " ", typeof(b3))
        # println("==========")
        Flin = segClin' * b3
        # Matrix(Flin)
        sparse_filter_old!(Flin, 1, 1, 0)
        dropzeros!(Flin)
        return Flin, larmodel
    end

    """
    Construct B matrix and make multiplication with boundary3.
    B matrix composed from faces per block. The size of block is constant.
    """
    function grid_get_surface_Bchar_loc_fixed_block_size(segmentation::AbstractArray, block_size::AbstractArray{Int,1})
        data_size = size_as_array(size(segmentation))
        margin_size = 0
        block_number, blocks_number_axis = number_of_blocks_per_axis(
            data_size, block_size)

        tmp_img_size = blocks_number_axis::AbstractArray{Int, 1} .* block_size::Array{Int,1}
        # numF = grid_number_of_faces(block_size)
        numC = prod(block_size)
        # println("numC = $numC, block_number = $block_number")
        Slin = spzeros(Int8, numC, block_number)
        offsets = Array{Array,1}(undef, block_number)

        # block_size = size_as_array(size(segmentation))
        # oneS=1
        fixed_block_size=true
        for block_i=1:block_number
            block1, offset1, block_size1 = get_block(block_i,
                segmentation, block_size, margin_size, blocks_number_axis,
                fixed_block_size
            )
            oneS =  grid_to_linear(block1, 0)
            # println("Slin[$block_i, :] = ", oneS)
            # display(block1)
            Slin[:, block_i] = oneS[:, 1]
            offsets[block_i] = offset1
        end
        b3, larmodel = get_boundary3(block_size)
        # ------
        Bchar = Slin' * b3
        # Bchar is similar to Flin
        # ------
        sparse_filter!(Bchar, 1, 1, 0)
        dropzeros!(Bchar)
        return Bchar, offsets, blocks_number_axis, larmodel
        # return Slin, oneS, b3
    end

    """
    Construct B matrix and make multiplication with boundary3.
    B matrix composed from faces per block. The size of block is constant.
    """
    function grid_get_surface_Bchar_loc_fixed_block_size_parallel(segmentation::AbstractArray, block_size::AbstractArray{Int,1})
        # TODO parallelization
        data_size = size_as_array(size(segmentation))
        margin_size = 0
        block_number, blocks_number_axis = number_of_blocks_per_axis(
            data_size, block_size)

        tmp_img_size = blocks_number_axis::AbstractArray{Int, 1} .* block_size::Array{Int,1}
        # numF = grid_number_of_faces(block_size)
        numC = prod(block_size)
        # println("numC = $numC, block_number = $block_number")
        Slin = spzeros(Int8, numC, block_number)
        offsets = Array{Array,1}(undef, block_number)

        # block_size = size_as_array(size(segmentation))
        # oneS=1
        fixed_block_size=true
        @sync @distributed for block_i=1:block_number
            block1, offset1, block_size1 = get_block(block_i,
                segmentation, block_size, margin_size, blocks_number_axis,
                fixed_block_size
            )
            oneS =  grid_to_linear(block1, 0)
            # println("Slin[$block_i, :] = ", oneS)
            # display(block1)
            Slin[:, block_i] = oneS[:, 1]
            offsets[block_i] = offset1
        end
        b3, larmodel = get_boundary3(block_size)
        # ------
        Bchar = Slin' * b3
        # Bchar is similar to Flin
        # ------
        sparse_filter!(Bchar, 1, 1, 0)
        dropzeros!(Bchar)
        return Bchar, offsets, blocks_number_axis, larmodel
        # return Slin, oneS, b3
    end

    function __make_return(V, filteredFV, Flin, larmodel, return_all)
        if return_all
            return (V, [filteredFV]), Flin, larmodel
        else
            return (V, [filteredFV])
        end
    end

    """
    return
    surfaceLARmodel, Flin, fullLARmodel
    """
    function get_surface_grid(segmentation::AbstractArray; return_all::Bool=false)
        Flin, larmodel = grid_get_surf_Fvec_larmodel(segmentation)

        V = larmodel[1]
        FV = larmodel[2][3]
        inds1, inds2, vals = findnz(Flin)

        filteredFV = [FV[inds2[i]] for i=1:length(inds2)]

        # filteredFV = [FV[i] for i=1:length(Flin) if (Flin[i] == 1)]
        return __make_return(V, filteredFV, Flin, larmodel, return_all)
        # if return_all
        #     return (V, [filteredFV]), Flin, larmodel
        # else
        #     return (V, [filteredFV])
        # end
    end

    function get_surface_grid_old(segmentation::AbstractArray; return_all::Bool=false)
        Flin, larmodel = grid_get_surf_Fvec_larmodel_old(segmentation)
        V, topology = larmodel
        (VV, EV, FV, CV) = topology

        filteredFV = [FV[i] for i=1:length(Flin) if (Flin[i] == 1)]
        return __make_return(V, filteredFV, Flin, larmodel, return_all)
        # if return_all
        #     return (V, [filteredFV]), Flin, larmodel
        # else
        #     return (V, [filteredFV])
        # end
    end

    """
    Based on input segmentation and block size calculate filtered FV and full sparse FV.
    In sparse FV is number 1 where is the surface. There is also number 2 where is
    the edge between blocks.
    """
    function __grid_get_surface_Fchar_per_block(segmentation::AbstractArray, block_size::Array{Int,1})
        data_size = size_as_array(size(segmentation))
        numF = grid_number_of_faces(data_size)
        # print("size vs length vs grid_number_of_faces: ", szF, " ", lenF, " ", numF)
        # print("size vs length vs grid_number_of_faces: ", typeof(szF), " ", typeof(lenF), " ", typeof(numF))

        # block_size = [2, 3, 4]
        margin_size = 0
        block_number, blocks_number_axis = number_of_blocks_per_axis(
            data_size, block_size)

        bigFchar = spzeros(Int8, numF)
        # println("bigFchar ", size(bigFchar))
        Flin = Nothing
        # println("block number ", block_number, " block size: ", block_size, "margin size: ", margin_size,
        #     " block number axis: ", blocks_number_axis)
        for block_i=1:block_number
            block_seg, offset1, block_size1 = get_block(
                block_i,
                segmentation, block_size, margin_size, blocks_number_axis, false
            )

            # filteredFVi, Flin, (V, model) = get_surface_grid(segmentation1; return_all=true)
            # (VV, EV, FV, CV) = model
            # println("=== segmentation")
            # display(block_seg)
            Flin, larmodel = grid_get_surf_Fvec_larmodel(block_seg)

        # face from small to big
            i, j, v = findnz(Flin)
            # println("=== findnz 2")
            # display(Flin)
            # println("Flin i j v, ", length(i)," ", length(v), " bigFchar ", nnz(bigFchar))
            for fid in j
                big_fid, voxel_cart = sub_grid_face_id_to_orig_grid_face_id(data_size, block_size1, offset1, fid)
                # if it is 0 set it one
                # if it is 1 there are two faces (prev and current) so remove
                # if bigFchar[big_fid] == 0
                bigFchar[big_fid] = (bigFchar[big_fid] + 1) % 2
                # bigFchar[big_fid] += 1
                # print(".")
            end
        end
        dropzeros!(bigFchar)
        return bigFchar
    end

    """
        block_params = (block_i, segmentation, block_size, margin_size,
            blocks_number_axis, fixed_block_size)

        __grid_get_surface_get_Fids_used_in_block(block_i, block_getter...)
    """
    function __grid_get_surface_get_Fids_used_in_block(
        block_i, segmentation, block_size, margin_size, blocks_number_axis, fixed_block_size
        )
        # = block_params

        segmentation_block, offset1, block_size1 = get_block(block_i,
            segmentation, block_size, margin_size, blocks_number_axis, fixed_block_size
        )
        if fixed_block_size
            block_size1 = block_size
        end
        data_size = size_as_array(size(segmentation))
        Flin, larmodel = grid_get_surf_Fvec_larmodel(segmentation_block)

    # face from small to big
        i, j, v = findnz(Flin)
        return [
            sub_grid_face_id_to_orig_grid_face_id(data_size, block_size1, offset1, fid)[1]
            for fid in j
        ]
    end

    """
    Put Fids used in block into channel
        block_params = (block_i, segmentation, block_size, margin_size,
            blocks_number_axis, fixed_block_size)

        __grid_get_surface_get_Fids_used_in_block(block_i, block_getter...)
    """
    function __grid_get_surface_channel_Fids_used_in_block(
        channel,
        block_i, segmentation, block_size, margin_size, blocks_number_axis, fixed_block_size
        )
        # = block_params

        segmentation_block, offset1, block_size1 = get_block(block_i,
            segmentation, block_size, margin_size, blocks_number_axis, fixed_block_size
        )
        if fixed_block_size
            block_size1 = block_size
        end
        data_size = size_as_array(size(segmentation))
        Flin, larmodel = grid_get_surf_Fvec_larmodel(segmentation_block)

    # face from small to big
        i, j, v = findnz(Flin)
        for fid in j
            print("_")
            put!(channel,
                sub_grid_face_id_to_orig_grid_face_id(data_size, block_size1, offset1, fid)[1]
            )
        end
        put!(channel, -1)

    end
    function __fcn_doing_nothing(block)
        @info "info doing nothing"
        return block + 1
    end
    function __fcn_doing_print(block)
        println(block)
        return block
    end

    """
    Based on input segmentation and block size calculate filtered FV and full sparse FV.
    In sparse FV is number 1 where is the surface. There is also number 2 where is
    the edge between blocks.
    """
    function __grid_get_surface_Fchar_per_block_parallel_pmap(segmentation::AbstractArray, block_size::Array{Int,1}; fixed_block_size=false)
        data_size = size_as_array(size(segmentation))
        numF = grid_number_of_faces(data_size)

        block_number, bgetter = block_getter(segmentation, block_size; fixed_block_size=fixed_block_size)
        # data3d, block_size, u, blocks_per_axis, fixed_block_s = bgetter
        # block_number, blocks_number_axis = number_of_blocks_per_axis(
        #     data_size, block_size)

        # @everywhere function get_Fids(block_i)
        #     return __grid_get_surface_get_Fids_used_in_block(block_i, bgetter...)
        # end

        # get_Fids(block_i) = __grid_get_surface_get_Fids_used_in_block(block_i, bgetter...)
        # get_Fids = function(block_i)
        # end
        # variable1 = 2
        get_Fids = function(block_i)
            # __fcn_doing_print(bgetter)
             # return __grid_get_surface_get_Fids_used_in_block(block_i, data3d, block_size, u, blocks_per_axis, fixed_block_s)
             return __grid_get_surface_get_Fids_used_in_block(block_i, bgetter...)
            # return variable1 * __fcn_doing_nothing(block_i)

        end

        bigFchar = spzeros(Int8, numF)
        # println("bigFchar ", size(bigFchar))
        # Flin = Nothing

        faces_per_blocks = pmap(get_Fids, 1:block_number)
        for  block_i=1:block_number
            for big_fid in faces_per_blocks[block_i]

                bigFchar[big_fid] = (bigFchar[big_fid] + 1) % 2
            end
        end
        dropzeros!(bigFchar)
        return bigFchar
    end

    function __temp(ch, block_i, bgetter...)
        # This __grid_get ... function is invisible for other workers
        return __grid_get_surface_channel_Fids_used_in_block(ch, block_i, bgetter...)
    end

    """
    Based on input segmentation and block size calculate filtered FV and full sparse FV.
    In sparse FV is number 1 where is the surface. There is also number 2 where is
    the edge between blocks.
    """
    function __grid_get_surface_Fchar_per_block_parallel_channel(
        segmentation::AbstractArray, block_size::Array{Int,1},
        channel::RemoteChannel
        ;
        fixed_block_size=false
        )
        # println("begin get surface")
        data_size = size_as_array(size(segmentation))
        numF = grid_number_of_faces(data_size)

        block_number, bgetter = block_getter(segmentation, block_size; fixed_block_size=fixed_block_size)
        # c = Channel{Int64}(32)

        # @everywhere put_Fids(block_i) = __grid_get_surface_channel_Fids_used_in_block(ch, block_i, bgetter...)

        bigFchar = spzeros(Int8, numF)
        # println("bigFchar ", size(bigFchar))
        # put!(channel, 1)
        # put!(channel, 3)
        # put!(channel, -2)
        # put!(channel, 2)
        # println("put in channedl")
        Flin = Nothing
        # println("before distributed")
        @distributed for block_i=1:block_number
            # print(".")
            put!(channel, -3)
            __fcn_doing_print(block_i)
            # __grid_get_surface_channel_Fids_used_in_block(ch, block_i, bgetter...)
            # __temp(ch, block_i, bgetter...)
            # print(",")
            # put_Fids(block_i)
        end

        # println("after distributed")
        n = 0
        # println("parallel processing, expected n = $block_number")
        # while n < block_number
        while n < 3
            big_fid = take!(channel)
            # print("$big_fid,")
            if big_fid < 0
                n += 1
            else
                bigFchar[big_fid] = (bigFchar[big_fid] + 1) % 2
            end
        end
        dropzeros!(bigFchar)
        return bigFchar
    end
    """
    Compatibility with interface of all other get_Fchar functions
    """
    function __grid_get_surface_Fchar_per_block_parallel_channel(
        segmentation::AbstractArray, block_size::Array{Int,1}
        ;
        channel::RemoteChannel=nothing,
        fixed_block_size=false
        )
        if channel == nothing
            channel = RemoteChannel(()->Channel{Int}(32));
        end
        return __grid_get_surface_Fchar_per_block_parallel_channel(
            segmentation::AbstractArray, block_size::Array{Int,1},
            channel::RemoteChannel
            ;
            fixed_block_size=false
        )
    end

    """
    Based on input segmentation and block size calculate filtered FV and full sparse FV.
    In sparse FV is number 1 where is the surface. There is also number 2 where is
    the edge between blocks.
    """
    __grid_get_surface_Fchar_per_block_parallel = __grid_get_surface_Fchar_per_block_parallel_pmap

    function __grid_get_surface_Fchar_per_block_old_implementation(segmentation::AbstractArray, block_size::Array{Int,1})
        data_size = size_as_array(size(segmentation))
        numF = grid_number_of_faces(data_size)
        # print("size vs length vs grid_number_of_faces: ", szF, " ", lenF, " ", numF)
        # print("size vs length vs grid_number_of_faces: ", typeof(szF), " ", typeof(lenF), " ", typeof(numF))

        # block_size = [2, 3, 4]
        margin_size = 0
        block_number, blocks_number_axis = number_of_blocks_per_axis(
            data_size, block_size)

        # bigFchar = spzeros(Int8, lenF)
        bigFchar = spzeros(Int8, numF)
        Flin = Nothing
        @debug "block number ", block_number, " block size: ", block_size, "margin size: ", margin_size,
            " block number axis: ", blocks_number_axis
        for block_i=1:block_number
            block1, offset1, block_size1 = get_block(block_i,
                segmentation, block_size, margin_size, blocks_number_axis
            )
            segmentation1 = block1

            filteredFVi, Flin, (V, model) = get_surface_grid_old(segmentation1; return_all=true)
            (VV, EV, FV, CV) = model

        # face from small to big

            for fid=1:length(Flin)
                if (Flin[fid] == 1)

                    big_fid, voxel_cart = sub_grid_face_id_to_orig_grid_face_id(data_size, block_size1, offset1, fid)
                    bigFchar[big_fid] += 1

                end
            end
        end
        return bigFchar
    end

    """
    Calculate Flin using multiplication b3 * S.
    Output is Flin and new size of image data. It is bigger than the original one
    because of limitation of fixed block size.
    """
    function __grid_get_surface_Fchar_per_fixed_block_size(segmentation::AbstractArray, block_size::AbstractArray{Int,1})
        # B is F per bricks
        B, offsets, blocks_number_axis, larmodel1 = grid_get_surface_Bchar_loc_fixed_block_size(segmentation, block_size)
        data_size = size_as_array(size(segmentation))

        tmp_data_size = block_size .* blocks_number_axis
        numF = grid_number_of_faces(tmp_data_size)
        bigFchar = spzeros(Int8, numF)

        for nzout in zip(findnz(B)...)
            block_id, fid, val = nzout
            offset1 = offsets[block_id]
            big_fid, voxel_cart = sub_grid_face_id_to_orig_grid_face_id(tmp_data_size, block_size, offset1, fid)
            bigFchar[big_fid] = (bigFchar[big_fid] + 1) % 2

        end

        dropzeros!(bigFchar)
        return bigFchar, tmp_data_size
    end
    # =====================

    function get_surface_grid_per_block_general(
        segmentation::AbstractArray, block_size::ArrayOrTuple, face_extraction_fcn, model_extraction_fcn; return_all::Bool=false)
        Fchar = face_extraction_fcn(segmentation, block_size)
        data_size = size_as_array(size(segmentation))
        bigV, FVreduced = model_extraction_fcn(Fchar, data_size)
        if return_all
            return (bigV,[FVreduced]), Fchar, (bigV, model)
        else
            return (bigV,[FVreduced])
        end
    end

    get_surface_grid_per_block_full2(segmentation::AbstractArray, block_size::ArrayOrTuple; return_all::Bool=false
        ) = get_surface_grid_per_block_general(
            segmentation, block_size,
            __grid_get_surface_Fchar_per_block,
            grid_Fchar_to_V_FVfulltoreduced,
            )

    function get_surface_grid_per_block_full(segmentation::AbstractArray, block_size::ArrayOrTuple; return_all::Bool=false)
        Fchar = __grid_get_surface_Fchar_per_block(segmentation, block_size)

        data_size = size_as_array(size(segmentation))
        # bigV, (bigVV, bigEV, bigFV, bigCV) = Lar.cuboidGrid(data_size, true)
        # model = (bigVV, bigEV, bigFV, bigCV)
        bigV, FVreduced = grid_Fchar_to_V_FVfulltoreduced(Fchar, data_size)
        model = [FVreduced]

        # return filtered_bigFV, Fchar, (bigV, model)

        if return_all
            return (bigV,[FVreduced]), Fchar, (bigV, model)
        else
            return (bigV,[FVreduced])
        end
    end

    """
    Construction of FV is reduced. The V
    """
    function get_surface_grid_per_block_Vreduced_FVreduced(segmentation::AbstractArray, block_size::ArrayOrTupleOrNothing=nothing; return_all::Bool=false)
        block_size = get_global_block_size_if_necessary(block_size)
        Fchar = __grid_get_surface_Fchar_per_block(segmentation, block_size)
        # println("Fchar ", size(Fchar))

        data_size = size_as_array(size(segmentation))
        # bigV, (bigVV, bigEV, bigFV, bigCV) = Lar.cuboidGrid(data_size, true)
        # model = (bigVV, bigEV, bigFV, bigCV)
        bigV, FVreduced = grid_Fchar_to_Vreduced_FVreduced(Fchar, data_size)

        # return filtered_bigFV, Fchar, (bigV, model)

        if return_all
            return (bigV,[FVreduced]), Fchar, (bigV, [FVreduced])
        else
            return (bigV,[FVreduced])
        end
    end

    """
    Construction of FV is reduced. The V
    """
    function get_surface_grid_per_block_Vreduced_FVreduced_parallel(segmentation::AbstractArray, block_size::ArrayOrTupleOrNothing=nothing; return_all::Bool=false)
        block_size = get_global_block_size_if_necessary(block_size)
        Fchar = __grid_get_surface_Fchar_per_block_parallel_pmap(segmentation, block_size)
        # println("Fchar ", size(Fchar))

        data_size = size_as_array(size(segmentation))
        # bigV, (bigVV, bigEV, bigFV, bigCV) = Lar.cuboidGrid(data_size, true)
        # model = (bigVV, bigEV, bigFV, bigCV)
        bigV, FVreduced = grid_Fchar_to_Vreduced_FVreduced(Fchar, data_size)

        # return filtered_bigFV, Fchar, (bigV, model)

        if return_all
            return (bigV,[FVreduced]), Fchar, (bigV, [FVreduced])
        else
            return (bigV,[FVreduced])
        end
    end
    """
    Construction of FV is reduced. The V
    """
    function get_surface_grid_per_block_Vreduced_FVreduced_fixed_block_size(segmentation::AbstractArray, block_size::ArrayOrTupleOrNothing=nothinge; return_all::Bool=false)
        block_size = get_global_block_size_if_necessary(block_size)
        # grid_get_surface_Bchar_loc_fixed_block_size
        Fchar, new_data_size = __grid_get_surface_Fchar_per_fixed_block_size(segmentation, block_size)
        bigV, FVreduced = grid_Fchar_to_Vreduced_FVreduced(Fchar, new_data_size)

        if return_all
            return (bigV,[FVreduced]), Fchar, (bigV, [FVreduced])
        else
            return (bigV,[FVreduced])
        end
    end

    function get_surface_grid_per_block_Vreduced_FVreduced_old(segmentation::AbstractArray, block_size::ArrayOrTupleOrNothing=nothing; return_all::Bool=false)
        block_size = get_global_block_size_if_necessary(block_size)
        Fchar = __grid_get_surface_Fchar_per_block_old_implementation(segmentation, block_size)

        data_size = size_as_array(size(segmentation))
        # bigV, (bigVV, bigEV, bigFV, bigCV) = Lar.cuboidGrid(data_size, true)
        # model = (bigVV, bigEV, bigFV, bigCV)
        bigV, FVreduced = grid_Fchar_to_Vreduced_FVreduced_old(Fchar, data_size)

        # return filtered_bigFV, Fchar, (bigV, model)

        if return_all
            return (bigV,[FVreduced]), Fchar, (bigV, [FVreduced])
        else
            return (bigV,[FVreduced])
        end
    end

    function get_surface_grid_per_block_FVreduced(segmentation::AbstractArray, block_size::ArrayOrTupleOrNothing=nothing; return_all::Bool=false)
        block_size = get_global_block_size_if_necessary(block_size)
        Fchar = __grid_get_surface_Fchar_per_block(segmentation, block_size)

        data_size = size_as_array(size(segmentation))
        bigV, FVreduced = grid_Fchar_to_V_FVreduced(Fchar, data_size)

        # return filtered_bigFV, Fchar, (bigV, model)

        if return_all
            return (bigV,[FVreduced]), Fchar, (bigV, [FVreduced])
        else
            return (bigV,[FVreduced])
        end
    end

    function get_surface_grid_per_block(segmentation::AbstractArray, block_size::ArrayOrTuple; return_all::Bool=false)
        return get_surface_grid_per_block_Vreduced_FVreduced(segmentation, block_size; return_all=return_all)
    end

# end


text:
# C++

function cxxsetup()
    quote
        const cxx_filename = Ref{AbstractString}("code.cxx")
        function set_cxx_filename!(filename::AbstractString)
            return cxx_filename[] = filename
        end

        const cxx_chunks = Dict{AbstractString,Vector{String}}()
        function cxx_add_code!(cxx_code::AbstractString)
            return push!(get!(cxx_chunks, cxx_filename[], String[]), cxx_code)
        end
        function cxx_get_code()
            iobuffer = IOBuffer()
            allcode = Dict{AbstractString,String}()
            for (filename, chunks) in cxx_chunks
                for chunk in chunks
                    println(iobuffer, chunk)
                end
                allcode[filename] = String(take!(iobuffer))
            end
            return allcode
        end
        function cxx_write_code!()
            println("Generating C++ code:")
            allcode = cxx_get_code()
            for (filename, content) in allcode
                println("Generating $filename...")
                write(filename, content)
            end
            println("Done.")
            return nothing
        end
    end
end
export cxxsetup

function cxxnewfile(filename::AbstractString, cxx_stmts::AbstractString)
    cxx_code = """
               #include <ccomplex>
               #include <cstdint>

               $cxx_stmts
               """
    quote
        set_cxx_filename!($filename)
        cxx_add_code!($cxx_code)
    end
end
export cxxnewfile

function cxxcode(cxx_stmts::AbstractString)
    cxx_code = """
               $cxx_stmts
               """
    quote
        cxx_add_code!($cxx_code)
    end
end
export cxxcode

function cxxvariable(name::VarName, type::VarType, cxx_stmts::AbstractString)
    julia_code = quote
        const $(name.julia_name) = begin
            lib = CxxInterface.dlopen($(name.cxx_library))
            sym = CxxInterface.dlsym(lib, $(name.cxx_name); throw_error=false)
            if sym ≡ nothing
                libname = $(name.cxx_library)
                symname = $(name.cxx_name)
                println("Load error (library \"$libname\", symbol \"$symname\"). You can ignore this error when generating C++ code.")
                nothing
            else
                # ptr = cglobal(($(name.cxx_name), $(name.cxx_library)), $(type.julia_type))
                ptr = cglobal(sym, $(type.julia_type))
                res = unsafe_load(ptr)
                $(type.convert_to_final(:res))::$(type.final_julia_type)
            end
        end
    end
    julia_code = clean_code(julia_code)
    simple_julia_code = simplify_code(julia_code)

    cxx_code = """
        /*
        $(string(simple_julia_code))
        */
        extern "C" const $(type.cxx_type) $(name.cxx_name) = [] { $cxx_stmts }();
        """
    quote
        cxx_add_code!($cxx_code)
        $julia_code
    end
end
export cxxvariable

function cxxfunction(name::FnName, result::FnResult, arguments::AbstractVector{FnArg}, cxx_stmts::AbstractString)
    julia_code = quote
        function $(name.julia_name)($([:($(arg.julia_name)::$(arg.initial_julia_type)) for arg in arguments]...))
            res = ccall(($(name.cxx_name), $(name.cxx_library)), $(result.julia_type),
                        ($([arg.julia_type for arg in arguments if !arg.skip]...),),
                        $([arg.convert_from_initial(arg.julia_name) for arg in arguments if !arg.skip]...))
            return $(result.convert_to_final(:res))::$(result.final_julia_type)
        end
    end
    julia_code = clean_code(julia_code)
    simple_julia_code = simplify_code(julia_code)

    cxx_code = """
        /*
        $(string(simple_julia_code))
        */
        extern "C" $(result.cxx_type) $(name.cxx_name)(
            $(join(["$(arg.cxx_type) $(arg.cxx_name)" for arg in arguments if !arg.skip], ",\n    "))
        ) {
            $cxx_stmts
        }
        """
    quote
        cxx_add_code!($cxx_code)
        $julia_code
    end
end
export cxxfunction