text:
using SubmodularMaximization
using POMDPs
using MCTS
using PyPlot
using Base.Iterators
using Printf
using Statistics

close()

steps = 100
horizon = 5
iterations = 1000

show_coverage = false

print("Enter number of partitions/planning rounds: ")
num_partitions = parse(Int64, readline())

print("Enter number of robots: ")
num_robots = parse(Int64, readline())

num_targets = default_num_targets(num_robots)

grid = Grid(num_robots=num_robots)
sensor = CoverageSensor()

configs = MultiRobotTargetCoverageConfigs(grid=grid,
                                          sensor=sensor,
                                          horizon=horizon)

robot_states = map(x->random_state(grid), 1:num_robots)
target_states = map(x->random_state(grid), 1:num_targets)

solver = x->solve_n_partitions(num_partitions, x, threaded=true)

plot_state_space(grid)
xlim([0, grid.width+1])
ylim([0, grid.height+1])

for ii = 2:steps
  println("Step ", ii)

  prior_robot_states = Array(robot_states)

  @time solution = iterate_target_coverage!(robot_states=robot_states,
                                            target_states=target_states,
                                            configs=configs,
                                            solver=solver)

  global robot_states = solution.robot_states
  global target_states = solution.target_states
  trajectories = solution.trajectories

  #= TODO: Implement and refactor to compute weight matrix

  # Compute weights
  problem = MultiRobotTargetCoverageProblem(robot_states, target_states,
                                            configs)
  @time weight = total_weight(compute_weight_matrix(problem; threaded=true))

  @printf("Weight per robot: %0.3f\n", weight / num_robots )
  =#

  #
  # Plot results
  #

  plots=[]
  # Plot covered area
  for robot in robot_states
    circle = Circle([robot...], sensor.range_limit)
    append!(plots, SubmodularMaximization.plot_filled_circle(circle, color=:blue, alpha=0.1))
  end

  # Plot robots on top
  for robot in robot_states
    append!(plots, plot_quadrotor(robot, color=:blue))
  end

  for trajectory in trajectories
    append!(plots, plot_states(trajectory, color=:blue, linestyle=":"))
  end

  # Plot motion
  for (prior, current) in zip(prior_robot_states, robot_states)
    append!(plots, plot_states([prior, current], color=:blue, linestyle="-"))
  end

  for target in target_states
    covered = is_covered(sensor, target, robot_states)

    # Dim covered targets
    alpha = covered ? 0.3 : 1.0
    append!(plots, plot_target(target, alpha=alpha))
  end

  #
  # Print a summary of the problem state
  #
  println("Objective: ", solution.objective,
          " (", solution.objective / num_robots, " per robot)")
  n = num_covered(sensor, target_states, robot_states)
  println("Num. covered: ", n, " of ", num_targets, " (", n / num_targets, ")")

  line = readline()

  if line == "q"
    break
  end

  if ii < steps
    foreach(x->x.remove(), plots)
  end
end

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
@testset "Ops" begin
    @testset "add!" begin
        d = Dag()

        rn1 = node!(d, Int32)
        rn2 = node!(d, Int64)

        # complete form using combine! inetrnally
        parents = [
            ("left", rn1, latest(0), true),
            ("right", rn2, latest(0), false),
        ]
        an = add!(d, parents)
        @test Tick.eltype(an) == Int64  # promoted

        # short form assuming init value and always triggering
        an2 = add!(d, rn1, rn2)
        @test Tick.eltype(an2) == Int64  # promoted

        res = []
        onfire!(d, an.nid, v -> push!(res, v))

        res2 = []
        onfire!(d, an2.nid, v -> push!(res2, v))

        eve = EvalDfs()
        fire(eve, d, rn1.nid, 42)
        @test res == [42]
        @test res2 == [42]

        fire(eve, d, rn2.nid, 42)
        @test res == [42]
        @test res2 == [42, 84]

        fire(eve, d, rn1.nid, 41)
        @test res == [42, 83]
        @test res2 == [42, 84, 83]
    end

    @testset "first!" begin
        d = Dag()

        pn = node!(d, Tuple{Int,String})
        an = first!(d, pn)
        @test Tick.eltype(an) == Int

        res = []
        onfire!(d, an.nid, v -> push!(res, v))

        eve = EvalDfs()
        fire(eve, d, pn.nid, (42, "fortytwo"))
        @test res == [42]
    end

    @testset "last!" begin
        d = Dag()

        pn = node!(d, Tuple{Int,String})
        an = last!(d, pn)
        @test Tick.eltype(an) == String

        res = []
        onfire!(d, an.nid, v -> push!(res, v))

        eve = EvalDfs()
        fire(eve, d, pn.nid, (42, "fortytwo"))
        @test res == ["fortytwo"]
    end

    @testset "nth!" begin
        d = Dag()

        pn = node!(d, Tuple{Int,String,Float64})
        an = nth!(d, pn, 2)
        @test Tick.eltype(an) == String

        res = []
        onfire!(d, an.nid, v -> push!(res, v))

        eve = EvalDfs()
        fire(eve, d, pn.nid, (42, "fortytwo", 3.14))
        @test res == ["fortytwo"]
    end

    @testset "win!" begin
        d = Dag()

        pn = node!(d, Int)

        an = win!(d, pn, 3)
        @test Tick.eltype(an) == Vector{Int}

        wn = win!(d, pn, 3)
        ma = tr!(d, wn, mean)
        @test Tick.eltype(ma) == Float64
        # TODO is passing transformer here ok? shouldn't transformer just be
        #      a new node? it seems more efficient this way, but is it right?
        #      maybe only for this case we can enforce a Trait/Type on
        #      transformer

        res = []
        resf = []
        onfire!(d, an.nid, v -> push!(res, copy(v)))
        onfire!(d, ma.nid, v -> push!(resf, v))

        eve = EvalDfs()
        fire(eve, d, pn.nid, 1)
        @test res == [[1]]
        @test resf == [1.0]
        fire(eve, d, pn.nid, 2)
        @test res == [[1], [1,2]]
        @test resf ≈ [1.0, 1.5]
        fire(eve, d, pn.nid, 3)
        @test res == [[1], [1,2], [1,2,3]]
        @test resf ≈ [1.0, 1.5, 2.0]
        fire(eve, d, pn.nid, 4)
        @test res == [[1], [1,2], [1,2,3], [2,3,4]]
        @test resf ≈ [1.0, 1.5, 2.0, 3.0]
    end

    # TODO can we do something like?
    #   d |>> node!(Int) |>> win!(3) |>> tr!(Float, mean)
end

ext: jl


text:
using SincFun
using Base.Test

x = linspace(-1,1,101)

sf = sincfun(exp)

@test norm(sf[x]-exp.(x)) < 100eps(eltype(sf))

println("Time for 101 samples should be <0.002 seconds.")
@time sf[x]

@test sum(sf) ≈ 2sinh(1)

for p in [2,4,6]
    norm(sf,p) ≈ (2/p*sinh(p))^(1/p)
end

sfs = sincfun(x->sinpi(5x))
sfc = sincfun(x->cospi(5x))

@test norm(sfs[x]-5π*cumsum(sfc)[x]+sinpi(5)) < 200eps(eltype(sfs))

@test norm(roots(sfc) - linspace(-0.9,0.9,10)) < 200eps(eltype(sfc))

println("Time for 10 roots should be <0.80 seconds.")
@time roots(sfc)

h = ConformalMap(1.0,0.1.^(1:4))

@test norm(h[inv(h,x)]-x) < 10eps()

println("Time for 101 samples and inverses of a ConformalMap should be < 0.0002 seconds.")
@time h[inv(h,x)]

println("Testing sinc differentiation and indefinite integration matrices.")

k=collect(-5.0:5.0)
@test sinc.(-1,k) == [-0.020107164191308535,0.025030330116344923,-0.03309323761827199,0.04858833320985967,-0.08948987223608373,0.5,1.0894898722360837,0.9514116667901403,1.033093237618272,0.9749696698836551,1.0201071641913084]
@test sinc.(0,k) == [0.0,0.0,0.0,0.0,0.0,1.0,-0.0,-0.0,-0.0,-0.0,-0.0]
@test sinc.(1,k) == [0.2,-0.25,0.3333333333333333,-0.5,1.0,0.0,-1.0,0.5,-0.3333333333333333,0.25,-0.2]
@test sinc.(2,k) == [0.08,-0.125,0.2222222222222222,-0.5,2.0,-3.289868133696453,2.0,-0.5,0.2222222222222222,-0.125,0.08]

k=collect(-500.0:500.0)
@test all(Bool[norm(sinc.(i,k)-convert(Vector{Float64},sinc.(i,convert(Vector{BigFloat},k))),Inf) < gamma(i+1)eps() for i=2:20])


text:
# This example predicts the temporal evolution of a field U,
# which has to be represented as vectors of matrices.
# Where the field comes from does not matter, but to make the example
# runnable we load one of the test systems of TimeseriesPrediction.
#
# This example uses light cone embedding and a nonlinear Barkley model.
#
# Importantly, the results are compared with the "real" evolution of the
# system.

# ### Simulate a test system
using PyPlot
using TimeseriesPrediction

testdir = dirname(dirname(pathof(TimeseriesPrediction)))*"/test"
@assert isdir(testdir)
include(testdir*"/system_defs.jl")

Ttrain = 300
Ttest = 5
T = Ttrain + Ttest

init = [ 0.6241    0.589685  0.668221   0.194882    0.687645
         0.656243  0.702544  0.476963   0.00236098  0.636111
         0.821854  0.868514  0.242682   0.2588      0.30552
         0.580972  0.355305  0.0805268  0.501724    0.728142
         0.297559  0.708676  0.583552   0.65363     0.555639]

U, V = barkley(T; tskip=100, ssize=(50,50), init = init)
summary(U)

# ### Temporal prediction of field U
γ = 2
τ = 1
r = 1
c = 1
bc = PeriodicBoundary()

pool = U[1 : Ttrain]
test  = U[ Ttrain : T]

em = light_cone_embedding(pool, γ,τ,r,c,bc)
pcaem = PCAEmbedding(pool, em; maxoutdim=5) # PCA speeds things up!

@time pred = temporalprediction(pool, em, Ttest; progress = false)

err = [abs.(test[i]-pred[i]) for i=1:Ttest+1]

println("Maximum error: ", maximum(maximum(e) for e in err))

# ### Plot prediction

# Deduce field maximum values
vmax = max(maximum(maximum(s) for s in test),
           maximum(maximum(s) for s in pred))
vmin = min(minimum(minimum(s) for s in test),
           minimum(minimum(s) for s in pred))

# plot plot plot
for i in [1, length(err)÷2, length(err)]

    fig = figure(figsize=(10,3))
    ax1 = subplot2grid((1,3), (0,0))
    ax2 = subplot2grid((1,3), (0,1))
    ax3 = subplot2grid((1,3), (0,2))

    im1 = ax1[:imshow](pred[i], cmap="viridis", vmin = vmin, vmax = vmax)
    im2 = ax2[:imshow](test[i], cmap="viridis", vmin = vmin, vmax = vmax)
    im3 = ax3[:imshow](err[i], cmap="inferno", vmin = 0, vmax = vmax-vmin)
    for (im, ax) in zip([im1,im2,im3], [ax1,ax2,ax3])
        ax[:get_xaxis]()[:set_ticks]([])
        ax[:get_yaxis]()[:set_ticks]([])
        colorbar(im, ax = ax, fraction=0.046, pad=0.04)#, format="%.1f")
    end
    ax1[:set_title]("Prediction")
    ax2[:set_title]("Real evolution")
    ax3[:set_title]("Absolute error")
    suptitle("frame $i")
    tight_layout(w_pad=0.6, h_pad=0.00001)
    subplots_adjust(top=0.75)
end
#md savefig("barkley_tempo.png"); nothing # hide
#md # ![barkley_tempo](barkley_tempo.png)


text:
export CircularQueue

mutable struct CircularQueue{T}
    data::Vector{T}
    length::Int64
    max_size::Int64
    head::Int64

    CircularQueue{T}(n::Int64) where T = new(Vector{T}(undef, n), 0, n, 0)
end

Base.length(cl::CircularQueue) = cl.length

Base.size(cl::CircularQueue) = (cl.max_size)

function Base.show(io::IO, cl::CircularQueue{T}) where T
    length = cl.length
    max_size = cl.max_size
    head = cl.head

    print(io, "CircularQueue{$(T), mem:$(length)/$(max_size)}([")

    if length < 1
        print(io, "])")
        return
    end

    for i = head.-collect(1:length) # change to 0-based index
        i = (i + max_size) % max_size # circular
        i += 1 # change to 1-based index
        print(io, "$(cl.data[i]), ")
    end

    print(io, "])")
end

function Base.pushfirst!(cl::CircularQueue, data)
    max_size = cl.max_size

    #=
        cl.head += 1 # move to next index
        cl.head -= 1 # change to 0-based index
        Skiped those to reduce calculation.
    =#
    cl.head %= max_size # circular
    cl.head += 1 # change to 1-based index

    cl.data[cl.head] = data
    cl.length += 1
    if cl.length > max_size
        cl.length = max_size
    end
end

function Base.pop!(cl::CircularQueue)
    max_size = cl.max_size
    length = cl.length

    if length < 1 return end

    index = cl.head
    index -= 1 # change to 0-based index
    index -= (length - 1) # move from head to tail
    index = (index + max_size) % max_size # circular
    index += 1 # change to 1-based index

    cl.length -= 1

    return cl.data[index]
end


text:
using LinearAlgebra

export Scale, factor

"""
    Scale{S <: Union{Number, Val}, N, BT <: AbstractBlock{N}} <: TagBlock{BT, N}

`Scale` a block with scalar. it can be either a `Number` or a compile time `Val`.

# Example

```jldoctest
julia> 2 * X
[scale: 2] X gate

julia> im * Z
[+im] Z gate

julia> -im * Z
[-im] Z gate

julia> -Z
[-] Z gate
```
"""
struct Scale{S <: Union{Number, Val}, N, BT <: AbstractBlock{N}} <: TagBlock{BT, N}
    alpha::S
    content::BT
end

content(x::Scale) = x.content
factor(x::Scale{<:Number}) = x.alpha
factor(x::Scale{Val{X}}) where X = X

Base.copy(x::Scale) = Scale(x.alpha, copy(x.content))
Base.adjoint(x::Scale{<:Number}) = Scale(adjoint(x.alpha), adjoint(content(x)))
Base.adjoint(x::Scale{Val{X}}) where X = Scale(Val(adjoint(X)), adjoint(content(x)))

Base.:(==)(x::Scale, y::Scale) = (factor(x) == factor(y)) && (content(x) == content(y))

chsubblocks(x::Scale, blk::AbstractBlock) = Scale(x.alpha, blk)
cache_key(x::Scale) = (factor(x), cache_key(content(x)))

mat(::Type{T}, x::Scale) where T = T(x.alpha) * mat(T, content(x))
mat(::Type{T}, x::Scale{Val{S}}) where {T, S} = T(S) * mat(T, content(x))

function apply!(r::ArrayReg{B}, x::Scale{S, N}) where {S, B, N}
    apply!(r, content(x))
    r.state .*= factor(x)
    return r
end


text:
using ..Cos
import Base: get
using ..Cos: CosDocImpl

mutable struct PDDocImpl <: PDDoc
    cosDoc::CosDocImpl
    catalog::CosIndirectObject{CosDict}
    pages::IDDN{CosDict}
    structTreeRoot::IDDN{CosDict}
    isTagged::Symbol #Valid values :tagged, :none and :suspect
    fonts::Dict{CosObject, PDFont}
    xobjs::Dict{CosObject, PDXObject}
    pager2n::Dict{CosIndirectObjectRef, Int}
    pagen2r::Dict{Int, CosIndirectObjectRef}
    function PDDocImpl(fp::AbstractString; access::Function)
        cosDoc = cosDocOpen(fp, access=access)
        catalog = cosDocGetRoot(cosDoc)
        new(cosDoc,catalog,CosNull,CosNull,:none,
            Dict{CosObject, PDFont}(), Dict{CosObject, PDXObject}(),
            Dict{CosIndirectObjectRef, Int}(),
            Dict{Int, CosIndirectObjectRef}())
    end
end

pdDocGetPage(doc::PDDocImpl, num::Int) =
    pdDocGetPage(doc, pd_doc_get_page(doc, num))

pdDocGetPage(doc::PDDoc, cosref::CosIndirectObjectRef) =
    create_pdpage(doc, cosDocGetObject(doc.cosDoc, cosref))

#=
```
    show(io::IO, doc::PDDoc)
```
Prints the PDDoc. The intent is to print lesser information from the structure.
=#
function Base.show(io::IO, doc::PDDoc)
    print(io, "\nPDDoc ==>\n")
    print(io, doc.cosDoc)
    print(io, '\n')
    print(io, "Catalog:")
    print(io, doc.catalog)
    print(io, "isTagged: $(doc.isTagged)\n")
end

#=
Recursively reads the page object and populates the indirect objects
Ensures indirect objects are read and updated in the xref Dictionary.
=#
function populate_doc_pages(doc::PDDocImpl, dict::CosIndirectObject{CosDict},
                            ncurr::Int)
    if (cn"Pages" == cosDocGetObject(doc.cosDoc, dict, cn"Type"))
        kids = cosDocGetObject(doc.cosDoc, dict, cn"Kids")
        arr = get(kids)
        len = length(arr)
        for i=1:len
            ref = splice!(arr, 1)
            obj = cosDocGetObject(doc.cosDoc, ref)
            if obj !== CosNull
                push!(arr, obj)
                ncurr = populate_doc_pages(doc, obj, ncurr)
            end
        end
    else
        ncurr += 1
        ref = CosIndirectObjectRef(dict)
        doc.pager2n[ref]   = ncurr
        doc.pagen2r[ncurr] = ref
    end
    parent = cosDocGetObject(doc.cosDoc, dict, cn"Parent")
    if (parent === CosNull)
        obj = cosDocGetObject(doc.cosDoc, parent)
        set!(dict, cn"Parent", obj)
    end
    return ncurr
end

populate_doc_pages(doc::PDDocImpl, dict::CosNullType, ncurr::Int) = nothing

pd_doc_get_pagenum(doc, pagenum::Int) = pagenum
pd_doc_get_pagenum(doc, pageref::CosIndirectObjectRef) = doc.pager2n[pageref]
pd_doc_get_page(doc, pagenum::Int) = doc.pagen2r[pagenum]

@inline function update_page_tree(doc::PDDocImpl)
    pages = cosDocGetObject(doc.cosDoc, doc.catalog, cn"Pages")
    populate_doc_pages(doc, pages, 0)
    doc.pages = pages
    return nothing
end

# The structure tree is not fully loaded but the object linkages are established for future
# correlations during text extraction.

@inline function update_structure_tree!(doc::PDDocImpl)
    catalog = pdDocGetCatalog(doc)
    marking = cosDocGetObject(doc.cosDoc, catalog, cn"MarkInfo")

    if (marking !== CosNull)
        tagged  = cosDocGetObject(doc.cosDoc, marking, cn"Marked")
        suspect = cosDocGetObject(doc.cosDoc, marking, cn"Suspect")
        doc.isTagged = (suspect === CosTrue) ? (:suspect) :
                       (tagged  === CosTrue) ? (:tagged)  : (:none)
    end

    doc.structTreeRoot = cosDocGetObject(doc.cosDoc, catalog, cn"StructTreeRoot")
    return nothing
end

function get_pd_font!(doc::PDDocImpl, cosfont::IDD{CosDict})
    font = get(doc.fonts, cosfont, nothing)
    font !== nothing && return font
    font = doc.fonts[cosfont] = PDFont(doc, cosfont)
    return font
end

function get_pd_xobject!(doc::PDDocImpl, cosxobj::CosObject)
    xobj = get(doc.xobjs, cosxobj, nothing)
    xobj !== nothing && return xobj
    xobj = doc.xobjs[cosxobj] = createPDXObject(doc, cosxobj)
    return xobj
end

function gather_sig_props(doc, fld, inherit)
    inhdown = Dict{Symbol, Any}()
    if cosDocGetObject(doc.cosDoc, fld, cn"FT") === cn"Sig"
        tobj = cosDocGetObject(doc.cosDoc, fld, cn"T")
        # Better to collect the page reference
        p_ref = cosDocGetObject(doc.cosDoc, fld, cn"P")
        page = p_ref === CosNull ? get(inherit, :P, CosNull) : p_ref
        tparent = get(inherit, :FQT, "")
        t = tobj !== CosNull ? CDTextString(tobj) : ""
        fqt = !isempty(tparent) && !isempty(t) ? tparent*"."*t :
            !isempty(tparent) ? tparent :
            !isempty(t) ? t : ""
        !isempty(fqt)    && (inhdown[:FQT] = fqt)
        page !== CosNull && (inhdown[:P]   = page)
    end
    return inhdown
end

function pd_get_signature_fields!(doc::PDDocImpl, fldroot::IDD{CosArray},
                                  inherit, sigflds)
    fldarr = get(fldroot)
    for fldobj in fldarr
        fld  = cosDocGetObject(doc.cosDoc, fldobj)
        inhdown = gather_sig_props(doc, fld, inherit)
        cosDocGetObject(doc.cosDoc, fld, cn"FT") === cn"Sig" &&
            push!(sigflds, (fld, inhdown))
        kids = cosDocGetObject(doc.cosDoc, fld, cn"Kids")
        kids !== CosNull &&
            pd_get_signature_fields!(doc, kids, inhdown, sigflds)
    end
end


text:
Pkg.update()
Pkg.add("JuMP")
Pkg.add("Ipopt")
Pkg.add("JSON")
Pkg.add("ArgParse")

max_stars_repo_name: exKAZUu-Research/SmartMotivator

max_issues_repo_name: exKAZUu-Research/SmartMotivator

max_forks_repo_name: exKAZUu-Research/WillingQuiz


text:
#
# QWTwPlot module
# qwt - based 2D plotting
#

#__precompile__(true)
module QWTWPlot

using Printf
using Libdl
using Sockets # for callbacks
import JSON3  # for settings file
import Qt_jll
import qwtw_jll
import marble_jll

function __init__()
	# this is not OK  qStart() # start in normal mode
	global cbTaskH # callback task handler
	cbTaskH = @task udpDataReader();
end

# DLLs and function handles below:
qwtwLibHandle = 0
qwtwFigureH = 0
qwtwServiceH = 0
qwtwClipGroupH = 0
qwtwRemoveLineH = 0
qwtwDebugMode = false
qwtwChangeLineH = 0
#qwtwFigure3DH = 0 # not supported because of license restrictions..    maybe later...
qwtwMapViewH = 0
qwtwsetimpstatusH = 0
qwtwCLearH = 0
qwtwPlotH = 0
qwtwPlot2H = 0
#qwtwPlot3DH = 0 # not supported because of license restrictions..    maybe later...
qwtwXlabelH = 0
qwtwYlabelH = 0
qwywTitleH = 0
qwtwVersionH = 0
qwtwMWShowH = 0
qwtEnableBroadcastH = 0
qwtDisableBroadcastH = 0

qwtMglH = 0
qwtMglLine = 0
qwtMglMesh = 0
#qwtSetCBH = 0

qwtStartH = 0
qwtStartDebugH = 0
qwtStopH = 0
started = false

old_path = ""
old_qtPath = ""
oldLdLibPath = ""

udpPort = 0 	# UDP port number
cfg = Dict()	# settings info
pleaseStopUdp = false
cbLock = ReentrantLock()
errrorBreak = false	# do 'error()' when parameters are not OK

"""
	information for the callback function
	about the mouse click
"""
struct QCBInfo
	type::Int32		# callback type ('1' for simple mouse click.. something else in case 'external UDP message info')
	plotID::Int32	# ID of the plot window
	lineID::Int32	# ID of the closest line
	index::Int32	# closest point index
	xx::Int32		# x window coord
	yy::Int32		# y window coord

	# closest point info
	x::Float64	# X coord
	y::Float64	# Y coord
	z::Float64  # Z coord (probably zero, when 'type' == 1)
	time::Float64 # time info
	label::String # label of the selected line
end
function QCBInfo()
	return QCBInfo(0, 0, 0, 0, 0, 0, 0.0, 0.0, 0.0, 0.0, "empty")
end
export QCBInfo

cbFunction  = function(info::QCBInfo) # 'picker' callback user function
	# doing nothing by default
end

"""
information about the 'clip' callback (when the user is pressing 'clip' button )
"""
struct QClipCallbackInfo
	t1::Float64 	# time 1 (left)
	t2::Float64		# time 2 (right)
	clipGroup::Int32 # clip group (every plot has its group)
	havePos::Bool       # true if x y z are valid
	x1::Float64		# point corresponding to time1, if any
	y1::Float64
	z1::Float64
	x2::Float64 	# point corresponding to time2, if any
	y2::Float64
	z2::Float64
end
function QClipCallbackInfo()
	return QClipCallbackInfo(0.0, 0.0, 0, false, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
end

export QClipCallbackInfo

clipCallbackFunction = function(info::QClipCallbackInfo) # 'clip' callback user function
	# doing nothing by default
end

lastUdpPacket = zeros(UInt8, 80)
export lastUdpPacket # just for testing/debugging

"""
	onPickUdp(x::Vector{UInt8})

process picker UDP message
x is the message
"""
function onPickUdp(x::Vector{UInt8})
	global cbFunction
	global qwtwDebugMode
	global qwtwServiceH
	ii = QCBInfo()

	if qwtwDebugMode
		@printf "onPickUdp starting!\n"
	end

	try
		time = reinterpret(Float64, x[5:12])[1]
		iks =  reinterpret(Float64, x[13:20])[1]
		igrek =  reinterpret(Float64, x[21:28])[1]
		zet =  reinterpret(Float64, x[29:36])[1]

		index =  reinterpret(Int32, x[37:40])[1]

		xx = reinterpret(Int32, x[41:44])[1]
		yy =  reinterpret(Int32, x[45:48])[1]
		plotID =  reinterpret(Int32, x[49:52])[1]
		lineID = reinterpret(Int32, x[53:56])[1]
		type = reinterpret(Int32, x[57:60])[1]

		# and the 'label'
		eos = 84
		for k= 61:84
			if x[k] == 0
				eos = k
				break
			end
		end
		if eos == 61
			label = ""
		else
			label = String(x[61:eos-1])
		end
		if qwtwDebugMode
			@printf "\t type = %d; index = %d; label = %s\n" type index label
		end

		ii = QCBInfo(type, plotID, lineID, index, xx, yy, iks, igrek, zet, time, label)
	catch ex
		if qwtwDebugMode
			@printf "cannot process picker message\n "
			bt = backtrace()
			msg = sprint(showerror, ex, bt)
			@printf "%s\n" msg
		end
	end

	try
		stest1 = ccall(qwtwServiceH, Int32, (Int32,), 1);
		cbFunction(ii)  # call the user function
		stest2 = ccall(qwtwServiceH, Int32, (Int32,), 2);
	catch ex
		@printf "onPickUdp: callback function failed\n "
		bt = backtrace()
		msg = sprint(showerror, ex, bt)
		@printf "%s\n" msg
	end

	return
end

"""
	onClipUdp(x::Vector{UInt8})

process clip UDP message

"""
function onClipUdp(x::Vector{UInt8})
	global clipCallbackFunction
	global qwtwDebugMode
	global qwtwServiceH
	ii = QClipCallbackInfo()

	if qwtwDebugMode
		@printf "onClipUdp\n"
	end

	try
		time1 = reinterpret(Float64, x[5:12])[1]
		time2 = reinterpret(Float64, x[13:20])[1]
		clipGroup = reinterpret(Int32, x[21:24])[1]
		havePos = reinterpret(Int32, x[25:28])[1]

		x1 = reinterpret(Float64, x[29:36])[1]
		y1 = reinterpret(Float64, x[37:44])[1]
		z1 = reinterpret(Float64, x[45:52])[1]
		x2 = reinterpret(Float64, x[53:60])[1]
		y2 = reinterpret(Float64, x[61:68])[1]
		z2 = reinterpret(Float64, x[69:76])[1]

		ii = QClipCallbackInfo(time1, time2, clipGroup, havePos, x1, y1, z1, x2, y2, z2)
	catch ex
		if qwtwDebugMode
			@printf "onClipUdp: cannot process 'clip' message\n "
			bt = backtrace()
			msg = sprint(showerror, ex, bt)
			@printf "%s\n" msg
		end
	end

	if qwtwDebugMode
		@printf "\ttime1 = %f, time2 = %f\n" ii.t1  ii.t2
	end

	try
		stest1 = ccall(qwtwServiceH, Int32, (Int32,), 1);
		clipCallbackFunction(ii) # call the user function
		stest2 = ccall(qwtwServiceH, Int32, (Int32,), 2);
	catch ex
		@printf "onClipUdp: callback function failed\n "
		bt = backtrace()
		msg = sprint(showerror, ex, bt)
		@printf "%s\n" msg
	end

	return
end

"""
UDP reader task. this is a very simple UDP client.
"""
function udpDataReader()
	global udpPort
	global pleaseStopUdp
	global qwtwDebugMode
	global cbLock
	global lastUdpPacket
	global qwtwServiceH

	if udpPort == 0
		@printf " 1 cannot start udpDataReader \n"
		return
	end
	if qwtwServiceH == 0
		@printf " 2 cannot start udpDataReader \n"
		return
	end
	type = 0; plotID = 0; lineID = 0; index = 0; xx = 0; yy = 0; iks = 0.0; igrek = 0.0; zet = 0.0; time = 0.0; label = "hello";
	ii = QCBInfo(type, plotID, lineID, index, xx, yy, iks, igrek, zet, time, label)

	time1 = 0.0; time2 = 0.0
	x1 = 0.0; y1 = 0.0; z1 = 0.0;  x2 = 0.0; y2 = 0.0; z2 = 0.0; clipGroup = 0; havePos = false;
	clipInfo = QClipCallbackInfo(time1, time2, clipGroup, havePos, x1, y1, z1, x2, y2, z2)
	x = zeros(UInt8, 88)
	nx = 0; eos = 84; k = 1;

	if qwtwDebugMode
		@printf "starting udp data reader on port %d\n" udpPort
	end
	sock=UDPSocket()
	counter = 0
	bind(sock, ip"127.0.0.1", udpPort)
	lock(cbLock)
	while pleaseStopUdp == false
		#@printf "udpDataReader waiting for the data from UDP.. \n"
		x = recv(sock)
		#x=reinterpret(UInt8, packet)    # ??
		lastUdpPacket = x
		if pleaseStopUdp
			break
		end
		nx = length(x)
		if nx < 80
			if qwtwDebugMode
				@printf "udpDataReader() got %d bytes" nx
			end
			continue
		end

		if qwtwDebugMode
			@printf "got %d bytes  \n" length(x)
			print(typeof(x))
			@printf "  %s\n" x[1:4]
			#print(x)
		end

		if x[1] == 0x50 # 'P'
			onPickUdp(x)
		elseif x[1] ==  0x43 #  'C'
			onClipUdp(x)
		end

		counter +=1
    end  ## this loop will never terminate
    close(sock)
	unlock(cbLock)

	if qwtwDebugMode
		@printf "stopping udp data reader on port %d\n" udpPort
	end
end

#cbTaskH = Task(udpDataReader)

"""
function saveEnv()
	saves some of the env vars
"""
function saveEnv()
	global old_path, old_qtPath, oldLdLibPath
	try old_path = ENV["PATH"]; catch end
	try old_qtPath = ENV["QT_PLUGIN_PATH"]; catch end
	try oldLdLibPath = ENV["LD_LIBRARY_PATH"]; catch end
	return
end

"""
function restoreEnv()
	restores env vars back
"""
function restoreEnv()
	global old_path, old_qtPath, oldLdLibPath
	ENV["PATH"] = old_path
	ENV["QT_PLUGIN_PATH"] = old_qtPath
	ENV["LD_LIBRARY_PATH"] = oldLdLibPath
	return
end

"""
function printEnv()
	simply print some of the env variables
"""
function printEnv()
	try
		@printf "\n\tPATH = %s \n\n" String(ENV["PATH"])
		@printf "\n\tQT_PLUGIN_PATH = %s \n\n" String(ENV["QT_PLUGIN_PATH"])
		try
			ldp = String(ENV["LD_LIBRARY_PATH"])
			@printf "\n\tLD_LIBRARY_PATH = %s \n\n" ldp
		catch
			@printf "no LD_LIBRARY_PATH\n\n"
		end
	catch
		@printf "printEnv: not everything was printed \n"
	end
end

"""
function addEnvItem(item, var)
	add something to the ENV from the beginning
	item -  what to add
	var - where to add
"""
function addEnvItem(item, var::String; debug = false)
	dv = ":"
	@static if Sys.iswindows()
		dv = ";"
	end
	useDV = true
	try # check: do we have it?
		test = ENV[var]  # yes
	catch ex # no such variable
		useDV = false
	end

	item2Add = "" # what we will actually add
	if typeof(item) == String # simple case
		item2Add = item
		if debug
			@printf "adding [%s] to [%s] \n" item2Add  var
		end
	elseif typeof(item) == Base.RefValue{String} # this can happen also
		item2Add = item[]
		if debug
			@printf "adding [%s] to [%s] \n" item2Add  var
		end
	else  # strange case
		@printf "WARNING: trying to add following item to %s : " var
		print(item)
		try
			item2Add = String(item)
		catch ex
			item2Add = item
		end
	end

	try
		if useDV
			ENV[var] = item2Add * dv * ENV[var]
		else
			ENV[var] = item2Add
			if debug
				@printf " new ENV %s created \n " var
			end
		end
	catch ex
		@printf "ERROR while adding item to [%s]: " var
		print(ex)
		print(item)
	end
end

#=
function qCallback(fID, lineID, index, fx, fy, x, y, z, t, legend) :: Cvoid
	@printf "qCallback  t = %f\n" t
end

function qcbTest1() :: Cvoid
#	@printf "qcbTest1 \n"
end

function qcbTest2(index) :: Cvoid
	@printf "qcbTest2  %d \n"  index
end
=#

"""
	qstart(;debug = false, qwtw_test = false, libraryName = "libqwtw")::Int32
starts the qwtw "C" library.

Without it, nothing will work. Call it before any other functions.\\
if debug, will try to enable debug print out. A lot of info. Just for debugging.\\
qwtw_test	-> if true, will try to not modify env variables. \\
"""
function qstart(;debug = false, qwtw_test = false, libraryName = "libqwtw")::Int32
	global qwtwDebugMode
	qwtwDebugMode = debug
	qwtw_libName = "nolib"
	if debug
		@printf "starting qwtw; current path: %s\n\n" ENV["PATH"]
	end

	# this could be still useful:
	#if debug qwtw_libName *= "d"; end

	global qwtwLibHandle, qwtwFigureH, qwtwMapViewH,  qwtwsetimpstatusH, qwtwCLearH, qwtwPlotH, qwtwServiceH
	global qwtwPlot2H, qwtwXlabelH, qwtwYlabelH, qwywTitleH, qwtwVersionH, qwtwMWShowH, qwtwRemoveLineH
	global qwtEnableBroadcastH, qwtDisableBroadcastH, qwtwChangeLineH, qwtwClipGroupH
	#global qwtwPlot3DH, qwtwFigure3DH
	global qwtStartH, qwtStopH, started
	global old_path, old_qtPath, oldLdLibPath
	global qwtMglH, qwtMglLine, qwtMglMesh
	#global qwtSetCBH
	global udpPort
	global cfg
	global pleaseStopUdp

	if started
		@printf "qwtw already started\n"
		return 0
	end

	settingsFileName = joinpath(homedir(), ".qwtw", "settings.json")
	try
		settings = read(settingsFileName)
		cfg = JSON3.read(String(settings))
		udpPort = parse(Int32, cfg["udp_client_port"])
	catch
		udpPort = 0
		@printf "error while reading settings file %s\n" settingsFileName
	end
	if udpPort != 0
		pleaseStopUdp = false
		#udpDataReader();
	end

	@static if Sys.iswindows() #  this part will handle OS differences
		qwtw_libName = libraryName * ".dll"
	else # hopefully this may be Linux
		if debug
			@printf "\t non-Windows detected\n"
		end
		qwtw_libName = libraryName * ".so"
	end

	saveEnv()

	marbleDataPath = joinpath(marble_jll.artifact_dir, "data")
	marblePluginPath = joinpath(marble_jll.artifact_dir, "plugins")

	if qwtw_test
		@printf "qwtw test; loading %s .. \n" qwtw_libName
		@printf "\nPATH = %s\n\n" ENV["PATH"]
	else
		@static if Sys.iswindows()
			ENV["QT_PLUGIN_PATH"]=Qt_jll.artifact_dir * "\\plugins"
			addEnvItem(qwtw_jll.PATH, "PATH", debug = debug)
			addEnvItem(qwtw_jll.LIBPATH, "PATH", debug = debug)

			#ENV["PATH"]= boost_jll.artifact_dir * "\\bin;" *  ENV["PATH"]
			#ENV["PATH"]= qwt_jll.artifact_dir * "\\bin;" *  ENV["PATH"]
			#ENV["PATH"]= Qt_jll.artifact_dir * "\\bin;" *  ENV["PATH"]

			#ENV["PATH"]= CompilerSupportLibraries_jll.artifact_dir * "\\bin;" *  ENV["PATH"]
			#new_env["PATH"]= FreeType2_jll.artifact_dir * "\\bin;" *  ENV["PATH"]
		else
			ENV["QT_PLUGIN_PATH"]=string(Qt_jll.artifact_dir) * "/plugins"
			addEnvItem(qwtw_jll.PATH, "PATH", debug = debug)
			addEnvItem(qwtw_jll.LIBPATH, "LD_LIBRARY_PATH", debug = debug)

			#ENV["PATH"]= boost_jll.artifact_dir * "/bin;" *  ENV["PATH"]

			#ENV["LD_LIBRARY_PATH"] = string(Qt_jll.LIBPATH) * ":" * ENV["LD_LIBRARY_PATH"] # do not sure why this is not already there
			#ENV["LD_LIBRARY_PATH"] = string(CompilerSupportLibraries_jll.LIBPATH) * ":" * ENV["LD_LIBRARY_PATH"] # do not sure why this is not already there
			#ENV["LD_LIBRARY_PATH"] = string(FreeType2_jll.LIBPATH) * ":" * ENV["LD_LIBRARY_PATH"] # do not sure why this is not already there
		end
	end

	if debug
		@printf "\nloading %s \n" qwtw_libName
		@printf "corrected ENV:\n"
		@static if Sys.iswindows()

		else
			try
				@printf "LD_LIBRARY_PATH: %s \n\n" String(ENV["LD_LIBRARY_PATH"])
			catch
				@printf "NO LD_LIBRARY_PATH \n"
			end
		end
		@printf "PATH: %s \n\n" String(ENV["PATH"])
		@printf "\n"
	end

	try
		qwtwLibHandle = Libdl.dlopen(qwtw_libName)
	catch ex
		if debug
			printEnv();
		end
		restoreEnv()
		@printf "Sorry, dlopen for %s failed; something is wrong\n" qwtw_libName
		throw(ex)
	end

	if debug
		@printf "\nlibrary %s opened from %s \n" qwtw_libName  Libdl.dlpath(qwtwLibHandle)
	end
	qwtwFigureH = Libdl.dlsym(qwtwLibHandle, "qwtfigure")
	#qwtwFigure3DH = Libdl.dlsym(qwtwLibHandle, "qwtfigure3d")

	try
		qwtwMapViewH = Libdl.dlsym(qwtwLibHandle, "qwtmap")
	catch
		qwtwMapViewH = 0
		@printf "WARNING: topview functions disabled (looks like no [marble] support)\n"
	end

	qwtwsetimpstatusH = Libdl.dlsym(qwtwLibHandle, "qwtsetimpstatus")
	qwtwCLearH = Libdl.dlsym(qwtwLibHandle, "qwtclear")
	qwtwPlotH = Libdl.dlsym(qwtwLibHandle, "qwtplot")
	#qwtwPlot3DH = Libdl.dlsym(qwtwLibHandle, "qwtplot3d")
	qwtwPlot2H = Libdl.dlsym(qwtwLibHandle, "qwtplot2")
	qwtwXlabelH = Libdl.dlsym(qwtwLibHandle, "qwtxlabel")
	qwtwYlabelH = Libdl.dlsym(qwtwLibHandle, "qwtylabel")
	qwywTitleH = Libdl.dlsym(qwtwLibHandle, "qwttitle")
	qwtwVersionH = Libdl.dlsym(qwtwLibHandle, "qwtversion")
	qwtwMWShowH = Libdl.dlsym(qwtwLibHandle, "qwtshowmw")
	qwtStartH = Libdl.dlsym(qwtwLibHandle, "qtstart")
	qwtStartDebugH = Libdl.dlsym(qwtwLibHandle, "qtstart_debug")
	qwtStopH = Libdl.dlsym(qwtwLibHandle, "qwtclose")
	qwtwRemoveLineH = Libdl.dlsym(qwtwLibHandle, "qwtremove")
	qwtwChangeLineH = Libdl.dlsym(qwtwLibHandle, "qwtchange")

	qwtwServiceH = Libdl.dlsym(qwtwLibHandle, "qwtservice")

	try
		qwtwClipGroupH = Libdl.dlsym(qwtwLibHandle, "qwtclipgroup")
	catch
		@printf "WARNING: clip groups disabled (wrong qwtw version?) \n"
	end

	try
		qwtEnableBroadcastH = Libdl.dlsym(qwtwLibHandle, "qwtEnableCoordBroadcast")
		qwtDisableBroadcastH = Libdl.dlsym(qwtwLibHandle, "qwtDisableCoordBroadcast")
	catch
		@printf "WARNING: UDP broacast disabled \n"
	end

	try
		qwtMglH = Libdl.dlsym(qwtwLibHandle, "qwtmgl")
		qwtMglLine = Libdl.dlsym(qwtwLibHandle, "qwtmgl_line")
		qwtMglMesh = Libdl.dlsym(qwtwLibHandle, "qwtmgl_mesh")
	catch
		@printf "WARNING: 3D features disabled \n"
	end
#=
	try
		qwtSetCBH = Libdl.dlsym(qwtwLibHandle, "setcallback")

		qcbTest1H = Libdl.dlsym(qwtwLibHandle, "setcallback_t1")
		qcbTest2H = Libdl.dlsym(qwtwLibHandle, "setcallback_t2")

		cbTest1 = @cfunction(qcbTest1, Cvoid, ())
		ccall(qcbTest1H, Cvoid, (Ptr{Cvoid},), cbTest1);

		#if qwtSetCBH != 0
		#	cb_c = @cfunction(qCallback, Cvoid, (Cint, Cint, Cint, Cint, Cint,
		#		Cdouble, Cdouble, Cdouble, Cdouble, Cstring))
#
		#	ccall(qwtSetCBH, Cvoid, (Ptr{Cvoid},), cb_c);
		#end
	catch
		@printf "WARNING: setcallback features disabled \n"
	end
=#
	# hangs tight!!!!
	if debug
		@printf "starting qstart in debug mode \n"
		@printf "\tmarbleDataPath %s \n" marbleDataPath
		@printf "\tmarblePluginPath %s \n" marblePluginPath
	end
	test = 1
	try
		if debug
			test = ccall(qwtStartDebugH, Int32, (Ptr{UInt8}, Ptr{UInt8}, Int32),
				marbleDataPath, marblePluginPath, 10);
		else
			test = ccall(qwtStartH, Int32, (Ptr{UInt8}, Ptr{UInt8}),
				marbleDataPath, marblePluginPath); # very important to call this in the very beginning
		end
	catch ex
		if debug
			printEnv();
		end
		restoreEnv()
		@printf "Sorry, cannot start qwtw. Something is wrong\n"
		throw(ex)
	end
	#end
	#@printf "qtstart = %d \n" test

	#version = qversion();
	#println(version);
	if test == 0
		started = true
		if debug
			version = qversion();
			println(version);
		end
	else
		@printf "\nERROR: library was not started (return code %d) \n" test
		if debug
			printEnv();
		end
	end

	restoreEnv()

	return test
end

"""
	qstop()
close everything and detach from qwtw library.

It maybe useful for debugging. What it does actually, it sends a command to the QT process to exit.\\
Have to call `qstart()` before, though.
"""
function qstop()
	global qwtwLibHandle, qwtStopH
	global started
	global pleaseStopUdp

	pleaseStopUdp = true

	if qwtwLibHandle != 0
		if (qwtStopH != 0)
			ccall(qwtStopH,  Cvoid, ())
			#@printf "qwtw stopped\n"
		else
			@printf "error: was not started correctly\n"
		end
		Libdl.dlclose(qwtwLibHandle)
	else
		@printf "not started (was qstart() called?)\n"
		return
	end
	qwtwLibHandle = 0
	started = false
	return
end

"""
	qErrorBreak(b::Bool)

if b == true, then do 'error()' in case of bad input parameters
"""
function qErrorBreak(b::Bool)
	global errrorBreak
	errrorBreak = b
	return
end
export qErrorBreak

"""
	qversion() :: String

useful for debugging.\\
return version info (as a string);
"""
function qversion() :: String
	global qwtwVersionH, qwtwLibHandle
	if qwtwLibHandle == 0
		return "not started yet"
	end
	v = zeros(UInt8, 1024)
	#@printf "qwtwc version: "
	bs = ccall(qwtwVersionH, Int32, (Ptr{UInt8}, Int32), v, 1024);
	#return bytestring(pointer(v))
	cmd = unsafe_string(pointer(v), bs);
	#@printf " bs = %d .... " bs
	return cmd
end

# looks like not used anymore
function traceit( msg )
     # global g_bTraceOn
      if ( true )
         bt = backtrace() ;
         s = sprint(io->Base.show_backtrace(io, bt))
         println( "debug: $s: $msg" )
      end
end

"""
	qfigure(n)

create a new plot window, OR make plot (with this ID `n`) 'active'.

looks like `n` have to be an integer number (this plot ID, or zero if you do not care).\\
after this function, this plot is the "active plot". \\
If `n == 0` then another new plot will be created.
xAxisType and yAxisType could be :aLinear (default) or :aLog (logarithmic)

Returns ID of the created plot
"""
function qfigure(n::Integer = 0; xAxisType=:aLinear, yAxisType=:aLinear)::Int32
	global qwtwFigureH, started
	if !started
		@printf "not started (was qstart() called?)\n"
		return 0
	end
	flags::UInt32 = 0
	if xAxisType == :aLog
		flags |= 1
	end
	if yAxisType == :aLog
		flags |= 2
	end
	test = ccall(qwtwFigureH, Int32, (Int32, UInt32), n, flags);
	return test
end;

"""
	qclipgrp(gr)

set current 'clip group'.

gr is the ID of this new group.
all (new) subsequent plots/figires will belong to this group. if press 'clip' button,
then if will work only for plots from the same group.

useful in case some of the plots have another time range.
"""
function qclipgrp(gr)
	global qwtwClipGroupH, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return
	end
	if qwtwClipGroupH == 0
		@printf "clip groups are not supported, sorry\n"
		return
	end

	ccall(qwtwClipGroupH, Cvoid, (Int32,), Int32(gr));
end
export qclipgrp

"""
	qremove(id::Int32)

remove a line from a plot.
"""
function qremove(id::Int32)
	global qwtwRemoveLineH, started
	if !started
		@printf "not started (was qstart() called?)\n"
		return
	end
	ccall(qwtwRemoveLineH, Cvoid, (Int32,), id);
	return
end
export qremove

"""
	qfmap(n)

create a new  plot window to draw on a map (with specific window ID)\\
'n' is Int32
"""
function qfmap(n)::Int32
	global qwtwMapViewH, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return 0
	end
	if qwtwMapViewH == 0
		@printf "map vew not supported, sorry\n"
		return 0
	end
	test = 0
	try
		test = ccall(qwtwMapViewH, Int32, (Int32,), n);
	catch ex
		@printf "ERROR in QWTWPlot::qfmap()\n"
		println(ex)
	end
	return test
end

"""
	qfmap()

create a new  plot window to draw on a map.
"""
function qfmap()::Int32
	return qfmap(0)
end;

"""
	qmgl(n = 0)

create a new  plot window to draw 3D lines / surfaces.

currently works only for Linux (?)
"""
function qmgl(n = 0)::Int32
	global qwtMglH, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return 0
	end
	if qwtMglH == 0
		@printf "3D not supported, sorry\n"
		return 0
	end

	test = ccall(qwtMglH, Int32, (Int32,), n);
	return test
end;

export qmgl

"""
	qmgline(x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64}, style::String = "-sb"; name = "")

draw a 3D line.

x, y, and z are the vectors with point coordinates.\\
style - how to draw a line. Explanation is here: http://mathgl.sourceforge.net/doc_en/Line-styles.html\\
name - a legend for this line
"""
function qmgline(x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64}, style::String = "-sb"; name = "")
	global qwtMglLine, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return
	end
	if qwtMglLine == 0
		@printf "3D not supported, sorry\n"
		return
	end
	if length(x) < 1
		@printf "QWTWPlot::qmgline empty X value\n"
		return
	end
	n = length(x)
	@assert length(x) == length(y)
	@assert length(z) == length(x)

	try
	ccall(qwtMglLine, Cvoid, (Int32, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{UInt8}, Ptr{UInt8}),
		n, x, y, z, name, style);

		sleep(0.025)
	catch
		@printf "ERROR in QWTWPlot::qmgline\n"
	end
	return
end
export qmgline

"""
	qmglmesh(data::Array{Float64, 2}, xMin = -10.0, xMax = 10.0, yMin = -10.0, yMax = 10.0,  style::String = ""; name = "", type = 0)

draw a 3D mesh/surface.   Currently works only for Linux (is this true?). \\

data - vertical (z) coordinates of the surface. should be Array{Float64, 2}\\
xMin, xMax, yMin, yMax is the definition of the range where to draw a surface. All the
coordinates from 'data' are evenly distributed inside this range.\\
style how to draw the surface. Some hints could be taken from here: http://mathgl.sourceforge.net/doc_en/Color-scheme.html \\
type  -   0 or 1;  mesh/grid  or surface\\
name - not used yet (maybe will add later)

"""
function qmglmesh(data::Array{Float64, 2}, xMin = -10.0, xMax = 10.0, yMin = -10.0, yMax = 10.0,  style::String = ""; name = "", type = 0)
	global qwtMglMesh, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return
	end
	if qwtMglMesh == 0
		@printf "3D not supported, sorry\n"
		return
	end
	if length(data) < 1
		@printf "QWTWPlot::qmgline empty X value\n"
		return
	end
	d = size(data)
	if (d[1] < 1) || (d[2] < 1)
		@printf "QWTWPlot::qmglmesh empty data\n"
		return
	end
	xSize = d[1]
	ySize = d[2]

	try
	ccall(qwtMglMesh, Cvoid, (Int32, Int32,
		Float64, Float64, Float64, Float64,
		Ptr{Float64},
		Ptr{UInt8}, Ptr{UInt8}, Int32),
		xSize, ySize,
		xMin, xMax, yMin, yMax, data,  name, style, type);

		sleep(0.025)
	catch
		@printf "ERROR in QWTWPlot::qwtMglMesh\n"
	end
	return
end

export qmglmesh

# create a new  window to draw a 3D points (QT engine)
#function qf3d(n)
#	global qwtwFigure3DH
#	ccall(qwtwFigure3DH, Cvoid, (Int32,), n);
#end;

"""
	qimportant(i)

set up an importance status for next lines.

looks like `0` means 'not important', `1` means "important.\\
'not important' will not participate in 'clipping':\\
	'not important' lines may be not completely inside the window.
"""
function qimportant(i)
	global qwtwsetimpstatusH, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return
	end

	ccall(qwtwsetimpstatusH, Cvoid, (Int32,), i);
	return
end

"""
	qclear()

close all the plots
"""
function qclear()
	global qwtwCLearH, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return
	end

	ccall(qwtwCLearH, Cvoid, ());
	return
end
#=
function qsetcallback(cb)
	global qwtSetCBH, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return
	end
	if qwtSetCBH == 0
		@printf "qwtSetCBH == 0 \n"
		return
	end

	#ccall(qwtSetCBH, Cvoid, (cb,));
	return
end
export qsetcallback
=#

"""
	qsmw()

open/show "main control window".
"""
function qsmw()
	global qwtwMWShowH, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return
	end

	ccall(qwtwMWShowH, Cvoid, ());
	return
end

"""
qplot(x::Vector{Float64}, y::Vector{Float64}, name::String, style::String, lineWidth, symSize):: Int32

plot normal lines.

#Parameters:

x and y:   the points.\\

name:      	name for this line.\\

style: 	 	how to draw a line.\\

lineWidth:	is a line width.\\

symSize:	size of the symbols, if they are used in 'style' spec.\\

what does 'style' parameter means? It's a string which has 1 or 2 or 3 symbols. \\

Look at two places for the detail explanation:\\
* example code  https://github.com/ig-or/QWTWPlot.jl/blob/master/src/qwexample.jl
* https://github.com/ig-or/QWTWPlot.jl/blob/master/docs/line-styles.md

this function returns ID of the line created. hopefully you can use it in some other functions\\
ID supposed to be >=0 is all is OK, or <0 in case of error
"""
function qplot(x::Vector{Float64}, y::Vector{Float64}, name::String, style::String="-b", lineWidth=1, symSize=1) :: Int32
	global qwtwPlotH, qwtwLibHandle, started, qwtwDebugMode
	global errrorBreak
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return -88
	end

	n = length(x)
	n1 = length(y)
	if  (n < 1) || (n != n1)
		mm = @sprintf "QWTWPlot qplot: bad input array length x[%d], y[%d]; name = %s" n n1 name
		if errrorBreak
			error(mm)
		else
			println(mm)
		end
		return -89
	end

	ww::Int32 = lineWidth;
	s::Int32 = symSize
	test::Int32 = -50

	if qwtwDebugMode
		@printf "qplot: (%s) (%s) n=%d \n " name style n
	end
	try
		test = ccall(qwtwPlotH, Int32, (Ptr{Float64}, Ptr{Float64}, Int32, Ptr{UInt8}, Ptr{UInt8}, Int32, Int32),
			x, y, n, name, style, ww, s);
		sleep(0.025)
	catch
		@printf "qplot: error #2  n = %d;  name = %s style = %s\n" n name style
		traceit("error #2")
		return -42
	end
	return test
end

"""
qplot1(x::Vector{Float64}, y::Vector{Float64}, name::String, style::String, symSize):: Int32

plot lines with line width == 1.

```

x and y:   the points
name:      	name for this line
style: 	 	how to draw a line
symSize:	size of the symbols

```

what does 'style' parameter means? It's a string which has 1 or 2 or 3 symbols. \\
Look at two places for the detail explanation:

* example code  https://github.com/ig-or/QWTWPlot.jl/blob/master/src/qwexample.jl
* https://github.com/ig-or/QWTWPlot.jl/blob/master/docs/line-styles.md

this function returns ID of the line created. hopefully you can use it in some other functions\\
ID supposed to be >=0 is all is OK, or <0 in case of error

"""
function qplot1(x::Vector{Float64}, y::Vector{Float64}, name::String, style::String, symSize) :: Int32
	global qwtwPlotH, qwtwLibHandle, started, qwtwDebugMode
	global errrorBreak
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return -88
	end

	n = length(x)
	n1 = length(y)
	if (n < 1) || (n1 < 1) || (n != n1)
		mm = @sprintf "QWTWPlot qplot1: bad input array length (x:%d; y:%d) name=%s style=%s symSize=%d\n" n n1 name style symSize
		if errrorBreak
			error(mm)
		else
			println(mm)
		end
		return -89
	end

	ww::Int32 = symSize;
	test::Int32 = -50
	if qwtwDebugMode
		@printf "qplot1: (%s) (%s) n=%d \n " name style n
	end

	try
		test = ccall(qwtwPlotH, Int32, (Ptr{Float64}, Ptr{Float64}, Int32, Ptr{UInt8}, Ptr{UInt8}, Int32, Int32),
			x, y, n, name, style, 1, ww);
	catch ex
		@printf "qplot1 ERROR\n"
		throw(ex)
		return -42
	end
	sleep(0.025)
	return test
end

"""
qchange(id::Int32, x::Vector{Float64}, y::Vector{Float64}, t::Vector{Float64} = []) :: Int32

change existing line.
"""
function qchange(id::Int32, x::Vector{Float64}, y::Vector{Float64}, t::Vector{Float64} = Vector{Float64}([])) :: Int32
	global qwtwChangeLineH, started
	global errrorBreak
	if !started
		@printf "not started (was qstart() called?)\n"
		return
	end
	n = length(x)
	n1 = length(y)
	n2 = length(t)

	if (n < 1) || (n != n1)
		mm = @sprintf "QWTWPlot qchange: bad input array length (x:%d; y:%d)\n" n n1
		if errrorBreak
			error(mm)
		else
			println(mm)
		end
		return -89
	end

	zz = Ptr{Float64}(C_NULL)
	tt = Ptr{Float64}(C_NULL)
	if n2 > 0
		tt = t
		if n2 != n
			mm = @sprintf "qchange: wrong _t_ array length (t: %d; n:%d)" n2  n
			if errrorBreak
				error(mm)
			else
				println(mm)
			end
			return 84
		end
	end
	test = -50
	try
		test = ccall(qwtwChangeLineH, Int32, (Int32, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Int32),
			id, x, y, zz, tt, n);
	catch ex
		@printf "qchange ERROR\n"
		throw(ex)
		return -42
	end
	return test
end
export qchange

# draw symbols in 3D space
# currently style and 'w' are not used
#function qplot3d(x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64},
#	 		name::String, style::String, w, time::Vector{Float64})
#	global qwtwPlotH
#	@assert (length(x) == length(y))
#	n = length(x)
#	ww::Int32 = w;
#	ccall(qwtwPlot3DH, Cvoid, (Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
#			Int32, Ptr{UInt8}, Ptr{UInt8}, Int32, Int32,  Ptr{Float64}),
#		x, y, z, n, name, style, 1, ww, time);
#	sleep(0.025)
#
#end;
#

"""
	qEnableCoordBroadcast(x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64}, time::Vector{Float64})

Enable UDP client.\\

In case you'd like to connect to some external software from this library,
it has UDP client and server inside. If/when you will move a marker, it will
broadcast 'marker information' via UDP. \\
Other way should also work: if some other application will broadcast UDP with marker info,
all the markers in this library supposed to be updated.

parameters: x, y,  and z are the vectors describing the line in 3D space;
and 'time' is an additional time information, for every point of the line.

It will send out marker info to UDP port 49561.\\
Incoming (to this library) UDP  port number is 49562, and IP address is "127.0.0.1".

Incoming message format:  'CRDS"<x><y><z>"FFFF', where
<x>, <y>, and <z> are  Float64 point coordinates, 8 bytes each,
so this message contains 4 + 3*8 + 4 bytes\\
outgoing message format is approximately the same: "EEEE"<x><y><z>"FFFF"

If UDP client is enabled, server will also start. For all this to work, some firewall rules have to be added probably,
line 'allow 49561 and 49562 ports for the local host'.
"""
function qEnableCoordBroadcast(x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64},
	 		time::Vector{Float64})
	global qwtEnableBroadcastH, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return
	end

	@assert (length(x) == length(y))
	@assert (length(x) == length(z))
	@assert (length(x) == length(time))

	n = length(time)

	if n < 1
		@printf "QWTWPlot qEnableCoordBroadcast: bad input array length (%d) \n" n
	end

	ccall(qwtEnableBroadcastH, Cvoid, (Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
			 Ptr{Float64}, Int32),
		x, y, z , time, n);
	sleep(0.025)
	return
end;

"""
	qDisableCoordBroadcast()
disable all this UDP -related stuff.\\

This function will stop UDP server and client.
"""
function qDisableCoordBroadcast()
	global qwtDisableBroadcastH, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return
	end

	ccall(qwtDisableBroadcastH, Cvoid, ());
	sleep(0.025);
	return
end;

"""
	qplot2(x::Array{Float64}, y::Array{Float64}, time::Array{Float64}, name::String, style::String, lineWidth=1, symSize=1):: Int32

plot with additional parameter (time?) info.

## Parameters:

```

'x', 'y': the line
'name' name of this line
'style'  the line style
'lineWidth' line width
'symSize' symbol size;
'time' time info, for every point

```

## Example

```

julia-repl
julia> time=collect(0.0:0.01:10.0)
x = sin.(time)
y = cos.(time)
qfigure()
qplot(time, x + y, "function 1", "-b", 3)
qfigure()
qplot2(x, y, time, "function 2", "-m", 3)

```
Now use marker on both plots and see that it moves on both plots.

this function returns ID of the line created. hopefully you can use it in some other functions\\
ID supposed to be >=0 is all is OK, or <0 in case of error

"""
function qplot2(x::Array{Float64}, y::Array{Float64}, time::Array{Float64}, name::String, style::String, lineWidth=1, symSize=1):: Int32
	global qwtwPlot2H, qwtwLibHandle, started, qwtwDebugMode
	global errrorBreak
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return -88
	end
	n = length(x)
	n1 = length(y)
	n2 = length(time)
	if (n < 1) || (n != n1) || (n != n2)
		mm = @sprintf "QWTWPlot qplot2 error: bad array length; x:%d; y:%d; time:%d; name = %s" n n1 n2 name
		if errrorBreak
			error(mm)
		else
			println(mm)
		end
		return -76
	end

	ww::Int32 = lineWidth;
	s::Int32 = symSize;
	test::Int32 = -50
	if qwtwDebugMode
		@printf "qplot2: (%s) (%s) n=%d \n " name style n
	end
	try
		test = ccall(qwtwPlot2H, Int32, (Ptr{Float64}, Ptr{Float64}, Int32, Ptr{UInt8}, Ptr{UInt8}, Int32, Int32, Ptr{Float64}),
			x, y, n, name, style, ww, s, time);
	catch ex
		@printf "qplot1 ERROR\n"
		throw(ex)
		return -42
	end

	sleep(0.025)
	return test
end

"""
	qxlabel(s::String)

put a label on the horizontal axis on the bottom.
"""
function qxlabel(s::String)
	global qwtwXlabelH, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return
	end

	ccall(qwtwXlabelH, Cvoid, (Ptr{UInt8},), s);
	return
end

"""
	qylabel(s::String)

put a label on the left vertical axis
"""
function qylabel(s::String)
	global qwtwYlabelH, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return
	end

	ccall(qwtwYlabelH, Cvoid, (Ptr{UInt8},), s);
	return
end

"""
	qtitle(s::String)

put a title on the current plot.
"""
function qtitle(s::String)
	global qwywTitleH, qwtwLibHandle, started
	if (qwtwLibHandle == 0) || (!started)
		@printf "not started (was qstart() called?)\n"
		return
	end

	ccall(qwywTitleH, Cvoid, (Ptr{UInt8},), s);
	return
end

"""
	qStarted()

does it started or not yet.
"""
function qStarted()
	global started
	return started
end

export qfigure, qfmap, qplot, qplot1, qplot2, qxlabel,  qylabel, qtitle
export qimportant, qclear, qstart, qstop, qversion, qsmw
export traceit
export  qEnableCoordBroadcast, qDisableCoordBroadcast
export qStarted
#export qplot3d, qf3d,

"""
	stopUdpThread()

stops UDP thread.
"""
function stopUdpThread()
	global pleaseStopUdp
	global udpPort
	global cbTaskH
	global qwtwDebugMode
	global cbLock

	if qwtwDebugMode
		@printf "stopUdpThread: istaskdone1(cbTaskH): %d\n" istaskdone(cbTaskH)
	end
	#if istaskstarted(cbTaskH) # stop the task
	if islocked(cbLock)
		if qwtwDebugMode
			@printf "stopping the task..\n"
		end
		pleaseStopUdp = true
		sock=UDPSocket()
		for i = 1:25
			send(sock, ip"127.0.0.1", udpPort, "stop\n")
		end

		if qwtwDebugMode
			@printf "waiting for the task to finish..\n"
		end
		#while !istaskdone(cbTaskH)
		#end

		#lock(cbLock)
		#unlock(cbLock)

		wait(cbTaskH)
		if qwtwDebugMode
			@printf "finished ..\n"
			@printf "stopUdpThread: istaskdone2(cbTaskH): %d\n" istaskdone(cbTaskH)
		end

		cbTaskH = @task udpDataReader();
	else

	end
end

function startUdpThread()
	global pleaseStopUdp, cbTaskH
	pleaseStopUdp = false
	schedule(cbTaskH)
	return
end

"""
	qsetCallback(cb)

set up another 'picker' callback
cb is the callback function, which will take one single parameter of type QCBInfo
see cbtest.jl for example how to use this.
"""
function qsetCallback(cb)
	global cbFunction
	global qwtwDebugMode
	#global cbLock

	stopUdpThread()

	#@printf "qsetCallback  locking.. \n"
	#lock(cbLock)
	#@printf "qsetCallback locked ! \n"
	cbFunction = cb
	if qwtwDebugMode
		@printf "setting up cbFunction \n"
	end
	#unlock(cbLock)
	#@printf "qsetCallback unlocked ! \n"
	startUdpThread()
	return
end
export qsetCallback

"""
	qsetClipCallback(cb)

set up a 'clip' callback
cb is the callback function, which will take one single parameter of type QClipCallbackInfo
see cbtest.jl for example how to use this.
"""
function qsetClipCallback(cb)
	global clipCallbackFunction
	global qwtwDebugMode

	stopUdpThread()
	clipCallbackFunction = cb
	if qwtwDebugMode
		@printf "setting up clipCallbackFunction! \n"
	end
	startUdpThread()
	return
end
export qsetClipCallback

"""
	setDebugMode(m::Bool)

for testing/defugging
"""
function setDebugMode(m::Bool)
	global qwtwDebugMode
	qwtwDebugMode = m
	return
end
export setDebugMode

end # module


text:
module JuliaMBD

import DifferentialEquations
import Plots

export toJulia, @xmlmodel, @modelcontents

include("TypeDef.jl")

include("Port.jl")
include("Line.jl")

include("predefined_blocks/AddBlock.jl")
include("predefined_blocks/SubBlock.jl")
include("predefined_blocks/GainBlock.jl")
include("predefined_blocks/IntegratorBlock.jl")
include("predefined_blocks/ModBlock.jl")
include("predefined_blocks/DivisionBlock.jl")
include("predefined_blocks/TrigonometricFunctionBlock.jl")
include("predefined_blocks/MathFunctionBlock.jl")
include("predefined_blocks/ProductBlock.jl")

include("predefined_blocks/ConstantBlock.jl")
include("predefined_blocks/StepBlock.jl")
include("predefined_blocks/RampBlock.jl")
include("predefined_blocks/PulseGeneratorBlock.jl")
include("predefined_blocks/SaturationBlock.jl")
include("predefined_blocks/QuantizerBlock.jl")

include("predefined_blocks/InBlock.jl")
include("predefined_blocks/OutBlock.jl")

include("SystemBlock.jl")
include("SimMacro.jl")

include("xmltojulia/xmltojulia.jl")

import .xmlToJulia

function toJulia(s)
    xmlToJulia.toJulia(s)
end

macro xmlmodel(fn)
    expr = quote
        f = open($fn, "r")
        data = read(f, String)
        close(f)
        expr = Meta.parse(toJulia(data))
        eval(expr)
    end
esc(expr)
end

macro modelcontents(fn)
    expr2 = quote
        f = open($fn, "r")
        data = read(f, String)
        close(f)
        println(toJulia(data))
    end
end

end


text:
"""
Estimate a linear model with high dimensional categorical variables / instrumental variables

### Arguments
* `df`: a Table
* `FormulaTerm`: A formula created using [`@formula`](@ref)
* `CovarianceEstimator`: A method to compute the variance-covariance matrix

### Keyword arguments
* `contrasts::Dict = Dict()` An optional Dict of contrast codings for each categorical variable in the `formula`.  Any unspecified variables will have `DummyCoding`.
* `weights::Union{Nothing, Symbol}` A symbol to refer to a columns for weights
* `save::Symbol`: Should residuals and eventual estimated fixed effects saved in a dataframe? Default to `:none` Use `save = :residuals` to only save residuals, `save = :fe` to only save fixed effects, `save = :all` for both. Once saved, they can then be accessed using `residuals()` or `fe()`. The returned DataFrame is automatically aligned with the original DataFrame.
* `method::Symbol`: A symbol for the method. Default is :cpu. Alternatively,  :gpu requires `CuArrays`. In this case, use the option `double_precision = false` to use `Float32`.
* `nthreads::Integer` Number of threads to use in the estimation. If `method = :cpu`, defaults to `Threads.nthreads()`. If `method = :gpu`, defaults to 256.
* `double_precision::Bool`: Should the demeaning operation use Float64 rather than Float32? Default to true.
* `tol::Real` Tolerance. Default to 1e-6.
* `maxiter::Integer = 10000`: Maximum number of iterations
* `drop_singletons::Bool = true`: Should singletons be dropped?
* `progress_bar::Bool = true`: Should the regression show a progressbar
* `first_stage::Bool = true`: Should the first-stage F-stat and p-value be computed?
* `dof_add::Integer = 0`:
* `subset::Union{Nothing, AbstractVector} = nothing`: select specific rows

### Details
Models with instruments variables are estimated using 2SLS. `reg` tests for weak instruments by computing the Kleibergen-Paap rk Wald F statistic, a generalization of the Cragg-Donald Wald F statistic for non i.i.d. errors. The statistic is similar to the one returned by the Stata command `ivreg2`.

### Examples
```julia
using RDatasets, FixedEffectModels
df = dataset("plm", "Cigar")
reg(df, @formula(Sales ~ NDI + fe(State) + fe(State)&Year))
reg(df, @formula(Sales ~ NDI + fe(State)*Year))
reg(df, @formula(Sales ~ (Price ~ Pimin)))
reg(df, @formula(Sales ~ NDI), Vcov.robust())
reg(df, @formula(Sales ~ NDI), Vcov.cluster(:State))
reg(df, @formula(Sales ~ NDI), Vcov.cluster(:State , :Year))
df.YearC = categorical(df.Year)
reg(df, @formula(Sales ~ YearC), contrasts = Dict(:YearC => DummyCoding(base = 80)))
```
"""

function reg(
    @nospecialize(df),
    @nospecialize(formula::FormulaTerm),
    @nospecialize(vcov::CovarianceEstimator = Vcov.simple());
    @nospecialize(contrasts::Dict = Dict{Symbol, Any}()),
    @nospecialize(weights::Union{Symbol, Nothing} = nothing),
    @nospecialize(save::Union{Bool, Symbol} = :none),
    @nospecialize(method::Symbol = :cpu),
    @nospecialize(nthreads::Integer = method == :cpu ? Threads.nthreads() : 256),
    @nospecialize(double_precision::Bool = true),
    @nospecialize(tol::Real = 1e-6),
    @nospecialize(maxiter::Integer = 10000),
    @nospecialize(drop_singletons::Bool = true),
    @nospecialize(progress_bar::Bool = true),
    @nospecialize(dof_add::Integer = 0),
    @nospecialize(subset::Union{Nothing, AbstractVector} = nothing),
    @nospecialize(first_stage::Bool = true))

    df = DataFrame(df; copycols = false)
    N = size(df, 1)

    ##############################################################################
    ##
    ## Parse formula
    ##
    ##############################################################################

    formula_origin = formula
    if  !omitsintercept(formula) & !hasintercept(formula)
        formula = FormulaTerm(formula.lhs, InterceptTerm{true}() + formula.rhs)
    end
    formula, formula_endo, formula_iv = parse_iv(formula)
    has_iv = formula_iv !== nothing
    has_weights = weights !== nothing

    ##############################################################################
    ##
    ## Save keyword argument
    ##
    ##############################################################################
    if save == true
        save = :all
    elseif save == false
        save = :none
    end
    if save ∉ (:all, :residuals, :fe, :none)
            throw("the save keyword argument must be a Symbol equal to :all, :none, :residuals or :fe")
    end
    save_residuals = (save == :residuals) | (save == :all)

    if method == :cpu && nthreads > Threads.nthreads()
        @warn "Keyword argument nthreads = $(nthreads) is ignored (Julia was started with only $(Threads.nthreads()) threads)."
        nthreads = Threads.nthreads()
    end
    ##############################################################################
    ##
    ## Construct new dataframe after removing missing values
    ##
    ##############################################################################

    # create a dataframe without missing values & negative weights
    vars = StatsModels.termvars(formula)
    iv_vars = Symbol[]
    endo_vars = Symbol[]
    if has_iv
        iv_vars = StatsModels.termvars(formula_iv)
        endo_vars = StatsModels.termvars(formula_endo)
    end
    # create a dataframe without missing values & negative weights
    all_vars = unique(vcat(vars, endo_vars, iv_vars))

    esample = completecases(df, all_vars)
    if has_weights
        esample .&= BitArray(!ismissing(x) && (x > 0) for x in df[!, weights])
    end
    if subset !== nothing
        if length(subset) != N
            throw("df has $(N) rows but the subset vector has $(length(subset)) elements")
        end
        esample .&= BitArray(!ismissing(x) && x for x in subset)
    end
    esample .&= Vcov.completecases(df, vcov)
    fes, ids, fekeys, formula = parse_fixedeffect(df, formula)
    has_fes = !isempty(fes)
    if has_fes
        if drop_singletons
            for fe in fes
                drop_singletons!(esample, fe)
            end
        end
    end
    save_fe = (save == :fe) | ((save == :all) & has_fes)

    nobs = sum(esample)
    (nobs > 0) || throw("sample is empty")

    if nobs == N
        esample = Colon()
    end

    # Compute weights
    if has_weights
        weights = Weights(convert(Vector{Float64}, view(df, esample, weights)))
        all(isfinite, weights) || throw("Weights are not finite")
    else
        weights = uweights(nobs)
    end

    # Compute feM, an AbstractFixedEffectSolver
    has_intercept = hasintercept(formula)
    has_fe_intercept = false
    if has_fes
        if any(fe.interaction isa UnitWeights for fe in fes)
            has_fe_intercept = true
        end
        fes = FixedEffect[fe[esample] for fe in fes]
        feM = AbstractFixedEffectSolver{double_precision ? Float64 : Float32}(fes, weights, Val{method}, nthreads)
    end
    # Compute data for std errors
    vcov_method = Vcov.materialize(view(df, esample, :), vcov)
    ##############################################################################
    ##
    ## Dataframe --> Matrix
    ##
    ##############################################################################
    exo_vars = unique(StatsModels.termvars(formula))
    subdf = Tables.columntable((; (x => disallowmissing(view(df[!, x], esample)) for x in exo_vars)...))
    s = schema(formula, subdf, contrasts)
    formula_schema = apply_schema(formula, s, FixedEffectModel, has_fe_intercept)

    # Obtain y
    # for a Vector{Float64}, convert(Vector{Float64}, y) aliases y
    y = convert(Vector{Float64}, response(formula_schema, subdf))
    all(isfinite, y) || throw("Some observations for the dependent variable are infinite")

    # Obtain X
    Xexo = convert(Matrix{Float64}, modelmatrix(formula_schema, subdf))
    all(isfinite, Xexo) || throw("Some observations for the exogeneous variables are infinite")

    response_name, coef_names = coefnames(formula_schema)
    if !(coef_names isa Vector)
        coef_names = typeof(coef_names)[coef_names]
    end

    if has_iv
        subdf = Tables.columntable((; (x => disallowmissing(view(df[!, x], esample)) for x in endo_vars)...))
        formula_endo_schema = apply_schema(formula_endo, schema(formula_endo, subdf, contrasts), StatisticalModel)
        Xendo = convert(Matrix{Float64}, modelmatrix(formula_endo_schema, subdf))
        all(isfinite, Xendo) || throw("Some observations for the endogenous variables are infinite")
        _, coefendo_names = coefnames(formula_endo_schema)
        append!(coef_names, coefendo_names)

        subdf = Tables.columntable((; (x => disallowmissing(view(df[!, x], esample)) for x in iv_vars)...))
        formula_iv_schema = apply_schema(formula_iv, schema(formula_iv, subdf, contrasts), StatisticalModel)
        Z = convert(Matrix{Float64}, modelmatrix(formula_iv_schema, subdf))
        all(isfinite, Z) || throw("Some observations for the instrumental variables are infinite")

        # modify formula to use in predict
        formula = FormulaTerm(formula.lhs, (tuple(eachterm(formula.rhs)..., (term for term in eachterm(formula_endo.rhs) if term != ConstantTerm(0))...)))
    end
    # compute tss now before potentially demeaning y
    tss_total = tss(y, has_intercept | has_fe_intercept, weights)
    # create unitilaized
    iterations, converged, r2_within = nothing, nothing, nothing
    F_kp, p_kp = nothing, nothing

    if has_fes
        # used to compute tss even without save_fe
        if save_fe
            oldy = deepcopy(y)
            if has_iv
                oldX = hcat(Xexo, Xendo)
            else
                oldX = deepcopy(Xexo)
            end
        end

        # initialize iterations and converged
        iterations = Int[]
        convergeds = Bool[]
        if has_iv
            Xall = Combination(y, Xexo, Xendo, Z)
        else
            Xall = Combination(y, Xexo)
        end

        _, iterations, convergeds = solve_residuals!(Xall, feM; maxiter = maxiter, tol = tol, progress_bar = progress_bar)

        iterations = maximum(iterations)
        converged = all(convergeds)
        if converged == false
            @info "Convergence not achieved in $(iterations) iterations; try increasing maxiter or decreasing tol."
        end

        tss_partial = tss(y, has_intercept | has_fe_intercept, weights)
    end

    if has_weights
        sqrtw = sqrt.(weights)
        y .= y .* sqrtw
        Xexo .= Xexo .* sqrtw
        if has_iv
            Xendo .= Xendo .* sqrtw
            Z .= Z .* sqrtw
        end
    end

    ##############################################################################
    ##
    ## Get Linearly Independent Components of Matrix
    ##
    ##############################################################################
    # Compute linearly independent columns + create the Xhat matrix
    if has_iv
        perm = 1:(size(Xexo, 2) + size(Xendo, 2))
        # first pass: remove colinear variables in Xendo
    	basis_endo = basis(eachcol(Xendo)...)
    	Xendo = getcols(Xendo, basis_endo)

    	# second pass: remove colinear variable in Xexo, Z, and Xendo
    	basis_all = basis(eachcol(Xexo)..., eachcol(Z)..., eachcol(Xendo)...)
        basis_Xexo = basis_all[1:size(Xexo, 2)]
        basis_Z = basis_all[(size(Xexo, 2) +1):(size(Xexo, 2) + size(Z, 2))]
        basis_endo_small = basis_all[(size(Xexo, 2) + size(Z, 2) + 1):end]
        if !all(basis_endo_small)
            # if adding Xexo and Z makes Xendo collinar, consider these variables are exogeneous
            Xexo = hcat(Xexo, getcols(Xendo, .!basis_endo_small))
            Xendo = getcols(Xendo, basis_endo_small)

            # out returns false for endo collinear with instruments
            basis_endo2 = trues(length(basis_endo))
            basis_endo2[basis_endo] = basis_endo_small

            # Change coef_names and oldX
            # TODO: I should probably also change formula in this case so that predict still works
            ans = 1:length(basis_endo)
            ans = vcat(ans[.!basis_endo2], ans[basis_endo2])
            perm = vcat(1:length(basis_Xexo), length(basis_Xexo) .+ ans)

            out = join(coefendo_names[.!basis_endo2], " ")
            @info "Endogenous vars collinear with ivs. Recategorized as exogenous: $(out)"

            # third pass
            basis_all = basis(eachcol(Xexo)..., eachcol(Z)..., eachcol(Xendo)...)
            basis_Xexo = basis_all[1:size(Xexo, 2)]
            basis_Z = basis_all[(size(Xexo, 2) +1):(size(Xexo, 2) + size(Z, 2))]
        end

    	Xexo = getcols(Xexo, basis_Xexo)
    	Z = getcols(Z, basis_Z)
        size(Z, 2) >= size(Xendo, 2) || throw("Model not identified. There must be at least as many ivs as endogeneous variables")
        basis_coef = vcat(basis_Xexo, basis_endo[basis_endo_small])

        # Build
        newZ = hcat(Xexo, Z)
        Pi = cholesky!(Symmetric(newZ' * newZ)) \ (newZ' * Xendo)
        Xhat = hcat(Xexo, newZ * Pi)
        X = hcat(Xexo, Xendo)

        # prepare residuals used for first stage F statistic
        ## partial out Xendo in place wrt (Xexo, Z)
        Xendo_res = BLAS.gemm!('N', 'N', -1.0, newZ, Pi, 1.0, Xendo)
        ## partial out Z in place wrt Xexo
        Pi2 = cholesky!(Symmetric(Xexo' * Xexo)) \ (Xexo' * Z)
        Z_res = BLAS.gemm!('N', 'N', -1.0, Xexo, Pi2, 1.0, Z)
    else
        # get linearly independent columns
        perm = 1:size(Xexo, 2)
        basis_Xexo = basis(eachcol(Xexo)...)
        Xexo = getcols(Xexo, basis_Xexo)
        Xhat = Xexo
        X = Xexo
        basis_coef = basis_Xexo
    end

    ##############################################################################
    ##
    ## Do the regression
    ##
    ##############################################################################

    crossx = cholesky!(Symmetric(Xhat' * Xhat))
    coef = crossx \ (Xhat' * y)

    ##############################################################################
    ##
    ## Optionally save objects in a new dataframe
    ##
    ##############################################################################
    residuals = y - X * coef
    residuals2 = nothing
    if save_residuals
        residuals2 = Vector{Union{Float64, Missing}}(missing, N)
        if has_weights
            residuals2[esample] .= residuals ./ sqrt.(weights)
        else
            residuals2[esample] .= residuals
        end
    end

    augmentdf = DataFrame()
    if save_fe
        oldX = getcols(oldX[:, perm], basis_coef)
        newfes, b, c = solve_coefficients!(oldy - oldX * coef, feM; tol = tol, maxiter = maxiter)
        for fekey in fekeys
            augmentdf[!, fekey] = df[:, fekey]
        end
        for j in eachindex(fes)
            augmentdf[!, ids[j]] = Vector{Union{Float64, Missing}}(missing, N)
            augmentdf[esample, ids[j]] = newfes[j]
        end
    end

    ##############################################################################
    ##
    ## Test Statistics
    ##
    ##############################################################################
    # Compute degrees of freedom
    dof_fes = 0
    if has_fes
        for fe in fes
            # adjust degree of freedom only if fe is not fully nested in a cluster variable:
            if (vcov isa Vcov.ClusterCovariance) && any(isnested(fe, v.groups) for v in values(vcov_method.clusters))
                dof_fes += 1 # if fe is nested you still lose 1 degree of freedom
            else
                #only count groups that exists
                dof_fes += nunique(fe)
            end
        end
    end
    dof_residual_ = max(1, nobs - size(X, 2) - dof_fes - dof_add)

    nclusters = nothing
    if vcov isa Vcov.ClusterCovariance
        nclusters = Vcov.nclusters(vcov_method)
    end

    # Compute standard error
    vcov_data = Vcov.VcovData(Xhat, crossx, residuals, dof_residual_)
    matrix_vcov = StatsBase.vcov(vcov_data, vcov_method)

    # Compute Fstat
    F = Fstat(coef, matrix_vcov, has_intercept)
    df_FStat_ = max(1, Vcov.df_FStat(vcov_data, vcov_method, has_intercept | has_fe_intercept))
    p = fdistccdf(max(length(coef) - (has_intercept | has_fe_intercept), 1), df_FStat_, F)
    # Compute Fstat of First Stage
    if has_iv && first_stage
        Pip = Pi[(size(Pi, 1) - size(Z_res, 2) + 1):end, :]
        try
            r_kp = Vcov.ranktest!(Xendo_res, Z_res, Pip,
                              vcov_method, size(X, 2), dof_fes)
            p_kp = chisqccdf(size(Z_res, 2) - size(Xendo_res, 2) + 1, r_kp)
            F_kp = r_kp / size(Z_res, 2)
        catch
            @info "ranktest failed ; first-stage statistics not estimated"
            p_kp, F_kp = NaN, NaN
        end
    end

    # Compute rss, tss, r2, r2 adjusted
    rss = sum(abs2, residuals)
    mss = tss_total - rss
    r2 = 1 - rss / tss_total
    adjr2 = 1 - rss / tss_total * (nobs - (has_intercept | has_fe_intercept)) / dof_residual_
    if has_fes
        r2_within = 1 - rss / tss_partial
    end

    ##############################################################################
    ##
    ## Return regression result
    ##
    ##############################################################################
    # add omitted variables
    if !all(basis_coef)
        newcoef = zeros(length(basis_coef))
        newmatrix_vcov = fill(NaN, (length(basis_coef), length(basis_coef)))
        newindex = [searchsortedfirst(cumsum(basis_coef), i) for i in 1:length(coef)]
        for i in eachindex(newindex)
            newcoef[newindex[i]] = coef[i]
            for j in eachindex(newindex)
                newmatrix_vcov[newindex[i], newindex[j]] = matrix_vcov[i, j]
            end
        end
        coef = newcoef
        matrix_vcov = Symmetric(newmatrix_vcov)
    end
    if any(perm[i] != i for i in perm)
        _invperm = invperm(perm)
        coef = coef[_invperm]
        newmatrix_vcov = zeros(size(matrix_vcov))
        for i in 1:size(newmatrix_vcov, 1)
            for j in 1:size(newmatrix_vcov, 1)
                newmatrix_vcov[i, j] = matrix_vcov[_invperm[i], _invperm[j]]
            end
        end
        matrix_vcov = Symmetric(newmatrix_vcov)
    end

    if esample == Colon()
        esample = trues(N)
    end

    return FixedEffectModel(coef, matrix_vcov, vcov, nclusters, esample, residuals2, augmentdf, fekeys, coef_names, response_name, formula_origin, formula, contrasts, nobs, dof_residual_,  rss, tss_total, r2, adjr2, F, p, iterations, converged, r2_within, F_kp, p_kp)
end


text:
# The SteinPolyharmonicSplineKernel for R^d is given by
#
# k(r) = Gamma(d/2 - m) / (2^{2m} pi^{d/2} (m-1)!) r^{2m-d}
#
# if d is odd and
#
# k(r) = (-1)^{m + d/2 - 1} / (2^{2m-1} pi^{d/2} (m-1)!(m-d/2)!) r^{2m-d} log r
#
# if d is even, where r = ||x-y||_2. This is a norm
# for the order m Beppo Levi Space on R^d.
#
# See Section 5.3 of http://arxiv.org/pdf/1204.6448.pdf.

type SteinPolyharmonicSplineKernel <: SteinKernel
    # the order of the derivative operator
    m::Int
end

# default will be 2
SteinPolyharmonicSplineKernel() = SteinPolyharmonicSplineKernel(2)

function getcons(m::Int, d::Int)
    if (d % 2 == 0)
        (-1)^(m + d/2 - 1) / (2^(2*m - 1) * pi^(d/2) *
            factorial(m - 1) * factorial(int(m - d/2)))
    else
        gamma(d/2 - m) / (2^(2*m) * pi^(d/2) * factorial(m - 1))
    end
end

function k(ker::SteinPolyharmonicSplineKernel, x::Array{Float64,1}, y::Array{Float64,1})
    m = ker.m; d = length(x)
    @assert m > (d/2)
    r = norm(x - y)
    if r == 0.0
        0.0
    end
    c = getcons(m, d)
    if (d % 2 == 0)
        c * r^(2*m - d) * log(r)
    else
        c * r^(2*m - d)
    end
end

function gradxk(ker::SteinPolyharmonicSplineKernel, x::Array{Float64,1}, y::Array{Float64,1})
    m = ker.m; d = length(x)
    @assert m > (d/2)
    r = norm(x - y)
    if r == 0.0
        0.0
    end
    c = getcons(m, d)
    l = 2*m - d
    if (d % 2 == 0)
        c * (x - y) * r^(l - 2) * (l * log(r) + 1.0)
    else
        c * l * (x - y) * r^(l - 2)
    end
end

function gradxyk(ker::SteinPolyharmonicSplineKernel, x::Array{Float64,1}, y::Array{Float64,1})
    m = ker.m; d = length(x)
    @assert m > (d/2)
    r = norm(x - y)

    c = getcons(m, d)
    l = 2*m - d
    if (d % 2 == 0)
        -c * r^(l - 2) * ((l^2 - 2*l + l*d)*log(r) + 2*l + d - 2)
    else
        -c * l * (2*m - 2) * r^(l - 2)
    end
end


text:
using Documenter
using TianJun

include("common.jl")

recursive_replace(@__DIR__, r"\.zh(?<ext>(\.md)?)$" => s"\g<ext>")

const BUILD_DIR = "build.zh"

makedocs(
    modules = [TianJun],
    format = Documenter.HTML(
        prettyurls = true,
        analytics = "UA-132847825-3",
        lang = "zh-CN",
        footer = "本站基于 [Documenter.jl](https://github.com/JuliaDocs/Documenter.jl) 和 [Julia 编程语言](https://julialang.org/) 构建，所有内容默认遵循[CC BY 4.0](https://creativecommons.org/licenses/by/4.0/)协议。英文版请访问[juntian.me](https://juntian.me)。",
        assets = [
            "assets/favicon.ico",
            "assets/custom.css",
        ]
    ),
    sitename = "田俊",
    linkcheck = haskey(ENV, "CI"),
    build = BUILD_DIR,
    pages = [
        "👋 关于" => "index.md",
        "💻 编程" => [
            "如何在Julia中计算点积?" => "programming/Dot_Product_in_Julia/index.md"
        ],
        "🙋 提问" => "AMA.md",
        "🔗 友链" => "blogroll.md",
        "🗃 存档" => [
            "20210812" => "essays/archive.md",
            hide("notebook_demo/pluto.md"),
            hide("essays/A_Draft_Design_of_Distributed_Reinforcement_Learning_in_Julia/index.md"),
            hide("essays/A_Dream/index.md"),
            hide("essays/A_Guide_to_Wrap_a_C++_Library_with_CXXWrap.jl_and_BinaryBuilder.jl_in_Julia/index.md"),
            hide("essays/A_Pratical_Guide_to_Distributed_Computing_in_Julia/index.md"),
            hide("essays/A_Short_Introduction_to_AI_for_children/index.md"),
            hide("essays/About_Recent/index.md"),
            hide("essays/About_This_Site/index.md"),
            hide("essays/All_About_Zipper/index.md"),
            hide("essays/An_Introduction_to_Flux.jl/index.md"),
            hide("essays/An_Introduction_to_Parallel_Computing_in_Julia_From_Bottom_Up/index.md"),
            hide("essays/An_Overview_of_Existing_Reinforcement_Learning_Libraries/index.md"),
            hide("essays/CUDA_in_One_Picture/index.md"),
            hide("essays/Categorical_Sampling_on_GPU_with_Julia/index.md"),
            hide("essays/Data_Science_London_Scikit-Learn/index.md"),
            hide("essays/From_Python_to_Julia/index.md"),
            hide("essays/Implementing_Readers_Writer_Lock_in_Julia/index.md"),
            hide("essays/Increase_Your_PPT_Skills/index.md"),
            hide("essays/Introduction_to_Deep_Learning_Libraries/index.md"),
            hide("essays/Learn_Some_Julia_Packages/index.md"),
            hide("essays/Lets_Play_Hanabi/index.md"),
            hide("essays/Locks_Threads_Tasks_Channels_Actors/index.md"),
            hide("essays/My_Interview_Questions/index.md"),
            hide("essays/Notes_on_Artificial_Intelligence_Foundations_of_Computational_Agents/index.md"),
            hide("essays/Notes_on_CS3110/index.md"),
            hide("essays/Notes_on_Computer_Age_Statistical_Inference/index.md"),
            hide("essays/Notes_on_Machine_Learning_A_Bayesian_and_Optimization_Perspective/index.md"),
            hide("essays/Notes_on_Statistical_Rethinking/index.md"),
            hide("essays/Oneday_Out_of_School/index.md"),
            hide("essays/Optional_in_Java_8/index.md"),
            hide("essays/Parallel_Computing/index.md"),
            hide("essays/Paraphrase_Generation/index.md"),
            hide("essays/Quantum_Computing/index.md"),
            hide("essays/Random_Thoughts_on_Chatbot/index.md"),
            hide("essays/Recursion_Memoize_Y-Combinator_in_Clojure/index.md"),
            hide("essays/Reinforcement_Learning_in_Action/index.md"),
            hide("essays/RelationExtraction/index.md"),
            hide("essays/Some_Interesting_Papers/index.md"),
            hide("essays/Some_Interesting_Quil_Examples/index.md"),
            hide("essays/Some_Thoughts_Recently/index.md"),
            hide("essays/Test/index.md"),
            hide("essays/The_Chinese_Translation_of_Bayesian_Analysis_with_Python/index.md"),
            hide("essays/The_EM_Algorithm/index.md"),
            hide("essays/The_Gambler_Ruin_Problem/index.md"),
            hide("essays/The_Past_Year_at_DidiChuxing/index.md"),
            hide("essays/The_Senior_Software_Engineer/index.md"),
            hide("essays/Thinking_in_NLP/index.md"),
            hide("essays/Throwing_Eggs_from_a_Building/index.md"),
            hide("essays/Understand_Buddy_in_Clojure/index.md"),
            hide("essays/Understanding_Variational_Autoencoder/index.md"),
            hide("essays/Write_Python_in_Lisp/index.md"),
            hide("essays/Write_a_Reinforcement_Learning_Package_in_Julia_from_Scratch/index.md"),
            hide("essays/[4]_Learn_Python_Together_Basic/index.md"),
            hide("essays/[6]_Learn_Python_Together_Data_Structure/index.md"),
        ]
    ]
)

cp(joinpath(@__DIR__, "CNAME"), joinpath(@__DIR__, BUILD_DIR, "CNAME");force=true)


text:
using Text, Stage, Ollam, DataStructures

function text(fn)
  res = ""
  for l in map(l -> chomp(l), eachline(`iconv -f latin1 -t utf8 $fn`))
    ismatch(r"^M[ -`]{60}$", l) && continue
    if res != ""
      res *= " " * l
    else
      res = l
    end
  end
  return res
end

function getinstances(dir)
  docs   = AbstractString[]
  truth  = AbstractString[]

  for t in filter(d -> d != "." && d != "..", readdir(dir))
    for d in filter(d -> d != "." && d != "..", readdir("$dir/$t"))
      push!(docs, "$dir/$t/$d")
      names = split(t, ".")
      push!(truth, names[end-1][1:1] * "." * names[end][1:min(4, end)])
    end
  end
  docs, truth
end

function tokenize(text)
  cleaned = english_tokenizer(text) # lowercase(text))
  res = filter(w -> !(lowercase(w) in english_stoplist), cleaned)
  res
end
tokenize_file(fn) = tokenize(text(fn))

# -------------------------------------------------------------------------------------------------------------------------
# download and prep data
# -------------------------------------------------------------------------------------------------------------------------
getfile("http://qwone.com/~jason/20Newsgroups/20news-bydate.tar.gz", "20news-bydate.tar.gz", expected_size = 14_464_277)
if !isdir("20ng")
  run(`tar -xzf 20news-bydate.tar.gz`)
  mkdir("20ng")
  mv("20news-bydate-train", "20ng/train")
  mv("20news-bydate-test", "20ng/test")
end
train, train_truth = getinstances("20ng/train")
test, test_truth   = getinstances("20ng/test")
@info "train: $(length(train)), test: $(length(test))"

# -------------------------------------------------------------------------------------------------------------------------
# topic training
# -------------------------------------------------------------------------------------------------------------------------
bkgmodel, fextractor, model = tc_train(train, train_truth, tokenize_file, mincount = 2, cutoff = 1e10,
                                       trainer = (fvs, truth, init_model) -> train_mira(fvs, truth, init_model, iterations = 20, k = 19, C = 0.01, average = true),
                                       iteration_method = :eager)
confmat = DefaultDict(AbstractString, DefaultDict{AbstractString, Int32}, () -> DefaultDict(AbstractString, Int32, 0))
res     = test_classification(model, lazy_map(x -> fextractor(tokenize_file(x)), test), test_truth,
                              record = (t, h) -> confmat[t][h] += 1) * 100.0
@info @sprintf("mira test set error rate: %7.3f", res)
print_confusion_matrix(confmat, width = 6)
@expect abs(res - 14.485) < 0.01

# List specific errors
# for (tr, t) in zip(test, test_truth)
#   trtext  = text(tr)
#   fv      = fextractor(tokenize(trtext))
#   scores  = score(model, fv)
#   bidx, b = best(scores)
#   if model.index_class[bidx] != t
#     @debug "ERROR: (ref: $t, hyp: $(model.index_class[bidx])) $(replace(trtext[1:70], r"\s+", " "))"
#   end
# end


text:
# VERSION >= v"0.4" && __precompile__(true)
module GeneNameGenieJ

using Neo4j, DataFrames, DataFramesMeta, Lazy;

global const DEFAULT_PORT = 7474;
global const DEFAULT_PATH = "db/data/";
global const DEFAULT_HOST = "localhost";

export
   setNeo4jConnection,
   getOfficialGeneSymbol,
   convertFromToExtended,
   convertFromTo,
   getValidDatabases,
   getValidAttributes,
   convertToCurrentMirbaseVersion,
   convertMatureMirnasToVersions,
   getValidMirnaMetadataValues,
   getCurrentMirbaseVersion,
   getEnsemblVersion;

function establishNeo4jConnection()

   # Try and connect to default host
   con = Neo4j.Connection(DEFAULT_HOST);
   global NEO4J_CONNECTION;
   NEO4J_CONNECTION = con;
   try
      graph = getgraph(con);
      println("Neo4j connection with default URL $DEFAULT_HOST successfully established.");
      true;
   catch
      warn("Was not able to connect to Neo4j database on default host $DEFAULT_HOST");
      false;
   end

   return con;
end

global NEO4J_CONNECTION = establishNeo4jConnection();

"""
    setNeo4jConnection(;host::String = DEFAULT_HOST, port::Integer = DEFAULT_PORT, path::String = DEFAULT_PATH)

Set desired global Neo4j connection. Default is `localhost:7474/db/data/`.

"""
function setNeo4jConnection(;host::T = DEFAULT_HOST, port::I = DEFAULT_PORT, path::T = DEFAULT_PATH) where
   {T<:String, I<:Integer}

   # Add leading and trailing slash if missing
   if !startswith(path, "/")
      path = "/" * path;
   end
   if !endswith(path, "/")
      path = path * "/";
   end

   con = Neo4j.Connection(host, port = port, path = path);
   # conUrl = con.host * ":" * string(con.port);
   try
      print("test connection for $(con.url)");
      graph = getgraph(con);
      println("Neo4j connection successfully established and set.");
      global NEO4J_CONNECTION;
      NEO4J_CONNECTION = con;
   catch
      warn("Was not able to connect to Neo4j database $(con.url)");
   end

end

# ids = String["P10415", "PUMA", "Bcl-2", "MYH14", "IGHV1-2", "AMPK", "ENST00000589955"];
# ids = String["Bcl-2", "MYH14", "IGHV1-2"];
# # ids = String["HOX1A", "CDK1", "ENST00000589950", "MYH14", "IGHV1-2", "PRKAA2", "HGNC:5550"];
# # dbs = String["HGNC", "GeneSymbolDB", "Uniprot/SWISSPROT", "Uniprot/SPTREMBL", "ArrayExpress", "Ens_Hs_transcript"];
# dbs = String["HGNC", "GeneSymbolDB", "Uniprot/SWISSPROT", "Uniprot/SPTREMBL"];
# # params = String["ensg:version", "ensg:region", "enst:region", "enst:regionStart", "enst:regionEnd", "aliases"];
# params = String["aliases"];

global addDbColDict = Dict{String, Symbol}(
   "ArrayExpress" => :EnsemblGeneId,
   "Ens_Hs_transcript" => :EnsemblTranscriptId,
   "Ens_Hs_translation" => :EnsemblProteinId
);

#########################################
# Include files
#########################################
include("validationFunctions.jl");
include("convertFromToExtendedParser.jl");
include("parameterFunctions.jl");
include("molecularIdentifierFunctions.jl");
include("mirnas.jl");

end # module

max_stars_repo_name: StefanHaunsberger/GeneNameGenieJ

max_issues_repo_name: StefanHaunsberger/GeneNameGenieJ

max_forks_repo_name: StefanHaunsberger/GeneNameGenieJ


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

"""
    Variable(name, [type])

A spatial variable with given `name` and machine `type`.
Default machine type is `Float64`.
"""
struct Variable
  name::Symbol
  type::DataType
end

Variable(name) = Variable(name, Float64)

"""
    name(variable)

Return the name of the spatial `variable`.
"""
name(var::Variable) = var.name

"""
    mactype(variable)

Return the machine type of the spatial `variable`.
"""
mactype(var::Variable) = var.type


text:
# # Tomography

#-

# Tomography is the process of reconstructing a density distribution from given
# integrals over sections of the distribution. In our example, we will
# work with tomography on black and white images.
# Suppose $x$ be the vector of $n$ pixel densities, with $x_j$
# denoting how white pixel $j$ is.
# Let $y$ be the vector of $m$ line integrals over the image, with $y_i$
# denoting the integral for line $i$.
# We can define a matrix $A$ to describe the geometry of the lines. Entry
# $A_{ij}$ describes how much of pixel $j$ is intersected by line $i$.
# Assuming our measurements of the line integrals are perfect, we have the relationship that
#
# $$
#   y = Ax
# $$
#
# However, anytime we have measurements, there are usually small errors that occur.
# Therefore it makes sense to try to minimize
#
# $$
#  \|y - Ax\|_2^2.
# $$
#
# This is simply an unconstrained least squares problem; something we can readily solve!

using Convex, ECOS, DelimitedFiles, SparseArrays
aux(str) = joinpath(@__DIR__, "aux_files", str) # path to auxiliary files
line_mat_x = readdlm(aux("tux_sparse_x.txt"))
summary(line_mat_x)

line_mat_y = readdlm(aux("tux_sparse_y.txt"))
summary(line_mat_y)

line_mat_val = readdlm(aux("tux_sparse_val.txt"))
summary(line_mat_val)

line_vals = readdlm(aux("tux_sparse_lines.txt"))
summary(line_vals)

# Form the sparse matrix from the data
# Image is 50 x 50
img_size = 50

# The number of pixels in the image
num_pixels = img_size * img_size

line_mat = spzeros(length(line_vals), num_pixels)

num_vals = length(line_mat_val)

for i in 1:num_vals
  x = Int(line_mat_x[i])
  y = Int(line_mat_y[i])
  line_mat[x + 1, y + 1] = line_mat_val[i]
end

pixel_colors = Variable(num_pixels)
## line_mat * pixel_colors should be close to the line_integral_values
## to reflect that, we minimize a norm
objective = sumsquares(line_mat * pixel_colors - line_vals)
problem = minimize(objective)
solve!(problem, () -> ECOS.Optimizer(verbose=0))

rows = zeros(img_size*img_size)
cols = zeros(img_size*img_size)
for i = 1:img_size
  for j = 1:img_size
    rows[(i-1)*img_size + j] = i
    cols[(i-1)*img_size + j] = img_size + 1 - j
  end
end

# Plot the image using the pixel values obtained!

using Plots
image = reshape(evaluate(pixel_colors), img_size, img_size)
heatmap(image, yflip=true, aspect_ratio=1, colorbar=nothing, color=:grays)


text:
module NativeFileDialog

using FilePathsBase
using NativeFileDialog_jll

export pick_file, pick_multi_file, pick_folder, save_file

@enum NFDResult NFD_ERROR = 0 NFD_OKAY = 1 NFD_CANCEL = 2

mutable struct NFDPathSet
    buf::Ptr{UInt8}
    indices::Ptr{Csize_t}
    count::Csize_t
end

NFDPathSet() = NFDPathSet(C_NULL, C_NULL, 0)

include("wrapper.jl")

Base.length(pathset::NFDPathSet) = pathset.count

function getindices(pathset::NFDPathSet)
    count = length(pathset)
    unsafe_load.(pathset.indices, 1:count)
end

function getpaths(pathset::NFDPathSet)
    unsafe_string.(pathset.buf .+ getindices(pathset))
end

"""
    pick_file(path=""; filterlist="")

Open an OS native window to select a file.

Returns the absolute path of the selected file or an empty string in case the
operation was cancelled.

# Examples
```jldoctest
julia> pick_file() # the selection window is open
"/home/user/file.txt" # MacOS/Linux

julia> pick_file()
"C\\:Users\\user\\Documents\\file.txt" # Windows

julia> pick_file()
"" # cancelled
```
"""
function pick_file(path = ""; filterlist = "")
    GC.@preserve path filterlist begin
        path = isempty(path) ? Cstring(C_NULL) : Cstring(pointer(path))
        filters = isempty(filterlist) ? Cstring(C_NULL) : Cstring(pointer(filterlist))
        status, outpath = OpenDialog(filters, path)
    end
    if status == NFD_ERROR
        error(unsafe_string(GetError()))
    elseif status == NFD_OKAY
        out = unsafe_string(outpath)
        Free(outpath)
    elseif status == NFD_CANCEL
        out = ""
    else
        error("Unkown status return ($status)")
    end
    return out
end

function pick_file(path::AbstractPath; filterlist = "")
    outpath = pick_file(string(path); filterlist)
    if isempty(outpath)
        return nothing
    end
    Path(outpath)
end

"""
    pick_multi_file(path=""; filterlist="")

Open an OS native window to select multiple files.

Returns an array of the absolute paths of the selected files or an empty vector
in case the operation was cancelled.

# Examples
```jldoctest
julia> pick_multi_file() # the selection window is open
["/home/user/file1.txt", "/home/user/file2.txt"] # MacOS/Linux

julia> pick_multi_file()
["C\\:Users\\user\\Documents\\file1.txt", "C\\:Users\\user\\Documents\\file2.txt"] # Windows

julia> pick_multi_file()
String[] # cancelled
```
"""
function pick_multi_file(path = ""; filterlist = "")
    GC.@preserve path filterlist begin
        path = isempty(path) ? Cstring(C_NULL) : Cstring(pointer(path))
        filters = isempty(filterlist) ? Cstring(C_NULL) : Cstring(pointer(filterlist))
        status, outpath = OpenDialogMultiple(filters, path)
    end
    if status == NFD_ERROR
        error(unsafe_string(GetError()))
    elseif status == NFD_OKAY
        out = getpaths(outpath)
        PathSetFree(outpath)
    elseif status == NFD_CANCEL
        out = String[]
    else
        error("Unkown status return ($status)")
    end
    return out
end

function pick_multi_file(path::AbstractPath; filterlist = "")
    outpathset = pick_multi_file(string(path); filterlist)
    if isempty(outpathset)
        return nothing
    end
    Path.(outpathset)
end

"""
    save_file(path=""; filterlist="")

Open an OS native window to set the location and name to save a file.

Returns the absolute path of the selected name and location for the file or an
empty string in case the operation was cancelled.

# Examples
```jldoctest
julia> save_file() # the selection window is open
"/home/user/file1.txt" # MacOS/Linux example

julia> save_file()
"C\\:Users\\user\\Documents\\file1.txt" # Windows

julia> save_file()
"" # cancelled
```
"""
function save_file(path = ""; filterlist = "")
    GC.@preserve path filterlist begin
        path = isempty(path) ? Cstring(C_NULL) : Cstring(pointer(path))
        filters = isempty(filterlist) ? Cstring(C_NULL) : Cstring(pointer(filterlist))
        status, outpath = SaveDialog(filters, path)
    end
    if status == NFD_ERROR
        error(unsafe_string(GetError()))
    elseif status == NFD_OKAY
        out = unsafe_string(outpath)
        Free(outpath)
    elseif status == NFD_CANCEL
        out = ""
    else
        error("Unkown status return ($status)")
    end
    return out
end

function save_file(path::AbstractPath; filterlist = "")
    outpath = save_file(string(path); filterlist)
    if isempty(outpath)
        return nothing
    end
    Path(outpath)
end

"""
    pick_folder(path="")

Open an OS native window to select a directory.

Returns the absolute path of the selected directory or an empty string in case
the operation was cancelled.

# Examples
```jldoctest
julia> pick_folder() # the selection window is open
"/home/user/dir" # MacOS/Linux

julia> pick_folder()
"C\\:Users\\user\\Documents\\dir" # Windows

julia> pick_folder()
"" # cancelled
```
"""
function pick_folder(path = "")
    GC.@preserve path begin
        path = isempty(path) ? Cstring(C_NULL) : Cstring(pointer(path))
        status, outpath = PickFolder(path)
    end
    if status == NFD_ERROR
        error(unsafe_string(GetError()))
    elseif status == NFD_OKAY
        out = unsafe_string(outpath)
        Free(outpath)
    elseif status == NFD_CANCEL
        out = ""
    else
        error("Unkown status return ($status)")
    end
    return out
end

function pick_folder(path::AbstractPath)
    outpath = pick_folder(string(path))
    if isempty(outpath)
        return nothing
    end
    Path(outpath)
end

end


text:
using StaticTrafficAssignment
using Test

@testset "StaticTrafficAssignment.jl" begin
    # Write your tests here.
end


text:
using Sellmeier
using Test

const testfile = joinpath(dirname(@__FILE__), "literate_org_tangled_tests.jl")
if isfile(testfile)
    include(testfile)
else
    error("Sellmeier not properly installed. Please run Pkg.build(\"Sellmeier\") then restart Julia.")
end