text:
#= MIT License

Copyright (c) 2020, 2021 Uwe Fechner

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

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

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

using GeometryBasics, Rotations, GLMakie, FileIO, LinearAlgebra, Printf
AbstractPlotting.__init__()

using Revise
includet("./Utils.jl")
using .Utils

includet("./Plot2D.jl")
using .Plot2D

const SCALE = 1.2
const INITIAL_HEIGHT =  80.0*se().zoom # meter, for demo
const MAX_HEIGHT     = 200.0*se().zoom # meter, for demo
const KITE = FileIO.load(se().model)
const FLYING    = [false]
const PLAYING    = [false]
const GUI_ACTIVE = [false]
const AXIS_LABEL_SIZE = 30
const TEXT_SIZE = 16
const running = Node(false)
const starting = [0]
const zoom = [1.0]
const textnode = Node("")
const textsize = Node(TEXT_SIZE)
const textsize2 = Node(AXIS_LABEL_SIZE)
const status = Node("")
const p1 = Node(Vector{Point2f0}(undef, 6000)) # 5 min
const p2 = Node(Vector{Point2f0}(undef, 6000)) # 5 min
const pos_x = Node(0.0f0)
const y_label1 = Node("")
const y_label2 = Node("")

const points          = Vector{Point3f0}(undef, se().segments+1)
const quat            = Node(Quaternionf0(0,0,0,1))                        # orientation of the kite
const kite_pos        = Node(Point3f0(1,0,0))                              # position of the kite
const positions       = Node([Point3f0(x,0,0) for x in 1:se().segments])   # positions of the tether segments
const part_positions  = Node([Point3f0(x,0,0) for x in 1:se().segments+1]) # positions of the tether particles
const markersizes     = Node([Point3f0(1,1,1) for x in 1:se().segments])   # includes the segment length
const rotations       = Node([Point3f0(1,0,0) for x in 1:se().segments])   # unit vectors corresponding with
                                                                           #   the orientation of the segments
const energy = [0.0]

function create_coordinate_system(scene, points = 10, max_x = 15.0)
    # create origin
    mesh!(scene, Sphere(Point3f0(0, 0, 0), 0.1 * SCALE), color=RGBf0(0.7, 0.7, 0.7))

    # create x-axis in red
    points += 2
    for x in range(1, length=points)
        mesh!(scene, Sphere(Point3f0(x * max_x/points, 0, 0), 0.1 * SCALE), color=:red)
    end
    mesh!(scene, Cylinder(Point3f0(-max_x/points, 0, 0), Point3f0(points * max_x/points, 0, 0), Float32(0.05 * SCALE)), color=:red)
    for i in range(0, length=points)
        start = Point3f0((points + 0.07 * (i-0.5)) * max_x/points, 0, 0)
        stop = Point3f0((points + 0.07 * (i+0.5)) * max_x/points, 0, 0)
        mesh!(scene, Cylinder(start, stop, Float32(0.018 * (10 - i) * SCALE)), color=:red)
    end

    # create y-axis in green
    points -= 3
    for y in range(0, length = 2points + 1)
        if y - points != 0
            mesh!(scene, Sphere(Point3f0(0, (y - points) * SCALE, 0), 0.1 * SCALE), color=:green)
        end
    end
    mesh!(scene, Cylinder(Point3f0(0, -(points+1) * SCALE, 0), Point3f0(0, (points+1) * SCALE , 0), Float32(0.05 * SCALE)), color=:green)
    for i in range(0, length=10)
        start = Point3f0(0, (points+1 + 0.07 * (i-0.5)) * SCALE, 0)
        stop = Point3f0(0, (points+1 + 0.07 * (i+0.5)) * SCALE, 0)
        mesh!(scene, Cylinder(start, stop, Float32(0.018 * (10 - i) * SCALE)), color=:green)
    end

    # create z-axis in blue
    points += 1
    for z in range(2, length=points)
        mesh!(scene, Sphere(Point3f0(0, 0, (z - 1) * SCALE), 0.1 * SCALE), color=:mediumblue)
    end
    mesh!(scene, Cylinder(Point3f0(0, 0, -SCALE), Point3f0(0, 0, (points+1) * SCALE), Float32(0.05 * SCALE)), color=:mediumblue)
    for i in range(0, length=10)
        start = Point3f0(0, 0, (points+1 + 0.07 * (i-0.5)) * SCALE)
        stop = Point3f0(0, 0, (points+1 + 0.07 * (i+0.5)) * SCALE)
        mesh!(scene, Cylinder(start, stop, Float32(0.018 * (10 - i) * SCALE)), color=:dodgerblue3)
    end
end

# draw the kite power system, consisting of the tether, the kite and the state (text and numbers)
function init_system(scene)
    sphere = Sphere(Point3f0(0, 0, 0), Float32(0.07 * SCALE))
    meshscatter!(scene, part_positions, marker=sphere, markersize=1.0, color=:yellow)
    cyl = Cylinder(Point3f0(0,0,-0.5), Point3f0(0,0,0.5), Float32(0.035 * SCALE))
    meshscatter!(scene, positions, marker=cyl, rotations=rotations, markersize=markersizes, color=:yellow)
    meshscatter!(scene, kite_pos, marker=KITE, markersize = 0.25, rotations=quat, color=:blue)
    if Sys.islinux()
        lin_font="/usr/share/fonts/truetype/ttf-bitstream-vera/VeraMono.ttf"
        if isfile(lin_font)
            font=lin_font
        else
            font="/usr/share/fonts/truetype/freefont/FreeMono.ttf"
        end
    else
        font="Courier New"
    end
    if se().fixed_font != ""
        font=se().fixed_font
    end
    text!(scene, textnode, position = Point3f0(-5.2, 3.5, -1), textsize = textsize, font=font, align = (:left, :top))
end

# update the kite power system, consisting of the tether, the kite and the state (text and numbers)
function update_system(scene, state, step=0)

    # move the particles to the correct position
    for i in range(1, length=se().segments+1)
        points[i] = Point3f0(state.X[i], state.Y[i], state.Z[i])
    end
    part_positions[] = [(points[k]) for k in 1:se().segments+1]

    # move, scale and turn the cylinder correctly
    positions[] = [(points[k] + points[k+1])/2 for k in 1:se().segments]
    markersizes[] = [Point3f0(1, 1, norm(points[k+1] - points[k])) for k in 1:se().segments]
    rotations[] = [normalize(points[k+1] - points[k]) for k in 1:se().segments]

    # move and turn the kite to the new position
    q0 = UnitQuaternion(state.orient)
    quat[]     = Quaternionf0(q0.x, q0.y, q0.z, q0.w)
    kite_pos[] = points[end]

    # print state values
    power = state.force * state.v_reelout
    energy[1] += (power / se().sample_freq * 2)
    if mod(step, 2) == 0
        height = points[end][3]/se().zoom
        msg = "time:      $(@sprintf("%7.2f", state.time)) s\n" *
            "height:    $(@sprintf("%7.2f", height)) m\n" *
            "elevation: $(@sprintf("%7.2f", state.elevation/pi*180.0)) °\n" *
            "azimuth:   $(@sprintf("%7.2f", state.azimuth/pi*180.0)) °\n" *
            "v_reelout: $(@sprintf("%7.2f", state.v_reelout)) m/s   " * "p_mech:  $(@sprintf("%8.2f", state.force*state.v_reelout)) W\n" *
            "force:     $(@sprintf("%7.2f", state.force    )) N     " * "energy:  $(@sprintf("%8.2f", energy[1]/3600)) Wh\n"
        textnode[] = msg
    end
end

function reset_view(cam, scene3D)
    update_cam!(scene3D.scene, [-15.425113, -18.925116, 5.5], [-1.5, -5.0, 5.5])
end

function zoom_scene(camera, scene, zoom=1.0f0)
    @extractvalue camera (fov, near, projectiontype, lookat, eyeposition, upvector)
    dir_vector = eyeposition - lookat
    new_eyeposition = lookat + dir_vector * (2.0f0 - zoom)
    update_cam!(scene, new_eyeposition, lookat)
end

function reset_and_zoom(camera, scene3D, zoom)
    reset_view(camera, scene3D)
    if ! (zoom ≈ 1.0)
        zoom_scene(camera, scene3D.scene, zoom)
    end
end

function main(gl_wait=true)
    scene, layout = layoutscene(resolution = (840+800, 900), backgroundcolor = RGBf0(0.7, 0.8, 1))
    scene3D = LScene(scene, scenekw = (show_axis=false, limits = Rect(-7,-10.0,0, 11,10,11), resolution = (800, 800)), raw=false)
    create_coordinate_system(scene3D)
    cam = cameracontrols(scene3D.scene)
    FLYING[1] = false
    PLAYING[1] = false
    GUI_ACTIVE[1] = true

    reset_view(cam, scene3D)

    textsize[]  = TEXT_SIZE
    textsize2[] = AXIS_LABEL_SIZE
    text!(scene3D, "z", position = Point3f0(0, 0, 14.6), textsize = textsize2, align = (:center, :center), show_axis = false)
    text!(scene3D, "x", position = Point3f0(17, 0,0), textsize = textsize2, align = (:center, :center), show_axis = false)
    text!(scene3D, "y", position = Point3f0( 0, 14.5, 0), textsize = textsize2, align = (:center, :center), show_axis = false)

    text!(scene, status, position = Point2f0( 20, 0), textsize = TEXT_SIZE, align = (:left, :bottom), show_axis = false)
    status[]="Stopped"

    layout[1, 1] = scene3D
    layout[2, 1] = buttongrid = GridLayout(tellwidth = false)
    layout[1, 2] = ax1 = Axis(scene, xlabel = "time [s]", ylabel = y_label1)
    layout[2, 2] = ax2 = Axis(scene, xlabel = "time [s]", ylabel = y_label2)
    linkxaxes!(ax1, ax2)

    l_sublayout = GridLayout()
    layout[1:3, 1] = l_sublayout
    l_sublayout[:v] = [scene3D, buttongrid]

    log = demo_log("Launch test!")
    plot2d(se, ax1, y_label1, log, p1, :height)
    lines!(ax1, p1)
    vlines!(ax1, pos_x, color = :red)

    plot2d(se, ax2, y_label2, log, p2, :elevation, true)
    lines!(ax2, p2)
    vlines!(ax2, pos_x, color = :red)

    btn_RESET       = Button(scene, label = "RESET")
    btn_ZOOM_in     = Button(scene, label = "Zoom +")
    btn_ZOOM_out    = Button(scene, label = "Zoom -")
    btn_LAUNCH      = Button(scene, label = "LAUNCH")
    btn_PLAY_PAUSE  = Button(scene, label = @lift($running ? "PAUSE" : " PLAY  "))
    btn_STOP        = Button(scene, label = "STOP")
    sw = Toggle(scene, active = false)
    label = Label(scene, "repeat")

    buttongrid[1, 1:8] = [btn_PLAY_PAUSE, btn_LAUNCH, btn_ZOOM_in, btn_ZOOM_out, btn_RESET, btn_STOP, sw, label]

    gl_screen = display(scene)

    init_system(scene3D)
    update_system(scene3D, demo_state(INITIAL_HEIGHT, 0))

    camera = cameracontrols(scene3D.scene)
    reset_view(camera, scene3D)

    reset() = reset_and_zoom(camera, scene3D, zoom[1])
    layout[3, 2] = bg = GridLayout(tellwidth = false, default_colgap=10)
    buttons(scene, bg, se, ax1, ax2, y_label1, y_label2, reset)

    on(btn_LAUNCH.clicks) do c
        if ! PLAYING[1]
            FLYING[1] = true
            PLAYING[1] = false
            status[] = "Launching..."
            reset_and_zoom(camera, scene3D, zoom[1])
        end
    end

    @async begin
        while true
            if starting[1] == 1
                starting[1] = 0
                plot2d(se, ax1, y_label1, log, p1, :height)
                if PLAYING[]
                    plot2d(se, ax2, y_label2, log, p2, :power, true)
                else
                    plot2d(se, ax2, y_label2, log, p2, :elevation, true)
                end
                x2 = log.extlog.time
                xlims!(ax2, x2[1], x2[end])
                reset_and_zoom(camera, scene3D, zoom[1])
            else
                sleep(0.1)
            end
        end
    end

    @async begin
        logfile=se().log_file * ".arrow"
        if ! isfile(logfile)
            status[] = "The logfile $logfile is missing! Importing..."; sleep(0.1)
            include("src/Importer.jl"); sleep(0.1)
            if isfile(logfile)
                status[] = "Success!"
            end
        else
            sleep(0.1)
        end
    end

    on(btn_PLAY_PAUSE.clicks) do c
        if status[] != "Launching..."
            if ! running[]
                logfile=se().log_file * ".arrow"
                if isfile(logfile)
                    running[] = true
                    status[] = "Running"
                    FLYING[1] = true
                    PLAYING[1] = true
                else
                    status[] = "Failed to import $logfile !"
                end
            else
                running[] = false
                status[] = "Paused"
            end
            reset_and_zoom(camera, scene3D, zoom[1])
        end
    end

    on(btn_RESET.clicks) do c
        reset_view(camera, scene3D)
        zoom[1] = 1.0
    end

    on(btn_STOP.clicks) do c
        # if status[] != "Stopped"
            FLYING[1] = false
            PLAYING[1] = false
            running[] = false
            status[] = "Stopped"
            reset_and_zoom(camera, scene3D, zoom[1])
        # end
    end

    on(scene.px_area) do x
        textsize[] = round(x.widths[2]/900.0 * TEXT_SIZE)
        textsize2[] = round(x.widths[2]/900.0 * AXIS_LABEL_SIZE)
        reset_view(camera, scene3D)
        zoom[1] = 1.0
    end

    on(btn_ZOOM_in.clicks) do c
        zoom[1] *= 1.2
        reset_and_zoom(camera, scene3D, zoom[1])
    end

    on(btn_ZOOM_out.clicks) do c
        zoom[1] /= 1.2
        reset_and_zoom(camera, scene3D, zoom[1])
    end

    # launch the kite on button click
    delta_t = 1.0 / se().sample_freq
    active = false
    log = demo_log("Launch test!")
    reset_view(camera, scene3D)

    simulation = @async begin
        while GUI_ACTIVE[1]
            # wait for launch command
            while ! FLYING[1] && GUI_ACTIVE[1]
                active = false
                sleep(0.05)
            end
            # load log file
            if ! active && GUI_ACTIVE[1]
                if PLAYING[1]
                    logfile=basename(se().log_file)
                    if log.name != logfile
                        old =  status[]
                        status[] = "Loading log file..."
                        reset_and_zoom(camera, scene3D, zoom[1])
                        log = load_log(logfile)
                        status[] = old
                    end
                else
                    log = demo_log("Launch test!")
                end
                steps = length(log.syslog)
                starting[1] = 1
                active = true
            end
            i=0
            energy[1] = 0.0
            # fly...
            while FLYING[1]
                state = log.syslog[i+1]
                if running[] || ! PLAYING[1]
                    @sync update_system(scene3D, state, i)
                    pos_x[] = i*delta_t
                    i += 1
                end
                sleep(delta_t / se().time_lapse)
                if i >= steps
                    if ! sw.active[]
                        FLYING[1] = false
                        PLAYING[1] = false
                        running[] = false
                        status[] = "Stopped"
                        reset_and_zoom(camera, scene3D, zoom[1])
                    else
                        i = 0
                    end
                end
            end
        end
    end

    if gl_wait
        wait(gl_screen)
    end

    # terminate the simulation
    FLYING[1] = false
    GUI_ACTIVE[1] = false
    running[] = false
    energy[1] = 0.0
    return nothing
end


text:
type ╱ <: SingleArg #\diagup
	x
end
Div=╱
function print(io::IO,d::Div)
	print(io,"/($(d.x))")
end
/(x::X,ex::Ex)=expression(Factor[x,Div(ex)])
function /(ex::Expression,x::Ex)
	ap=dcterms(ex)
	for t in ap
		push!(t,Div(x))
	end
	return expression(ap)
end
function divbine!(term::Term)
	dis=indsin(term,Div)
	if isempty(dis)
		return term
	end
	df=Factor[]
	for i in dis
		push!(df,term[i].x)
	end
	deleteat!(term,dis)
	push!(term,Div(simplify(expression(df))))
	return term
end
divbine(term::Term)=divbine!(deepcopy(term))
function divbinedify!(term::Term)
	if !isa(term[end],Div)
		return term
	end
	x=term[end].x
	if isa(x,Expression)
		if length(x)==1
			x=x[1]
		else
			return term
		end
	elseif !isa(x,Array)
		x=[x]
	end
	remove=Integer[]
	for fi in 1:length(x)
		i=indin(term[1:end-1],x[fi])
		if i!=0
			push!(remove,fi)
			deleteat!(term,i)
		elseif isa(x[fi],Number)
			push!(remove,fi)
			unshift!(term,1/x[fi])
		end
	end
	deleteat!(x,remove)
	if isempty(x)
		deleteat!(term,length(term))
	elseif length(x)==1
		term[end].x=x[1]
	else
		term[end].x=simplify(expression(x))
	end
	return term
end
divbinedify(term::Term)=divbinedify!(deepcopy(term))
function divify!(term::Array)
	dis=indsin(term,Div)
	remove=Int64[]
	termrem=Integer[]
	for i in dis
		if term[i].x==1
			term[i]=1
		elseif isa(term[i].x,Div)
			term[i]=term[i].x.x
		elseif isa(term[i].x,Expression)
			ap=terms(term[i].x)
			if length(ap)==1
				aprem=Integer[]
				for fac in 1:length(ap[1])
					if isa(ap[1][fac],Div)
						push!(term,ap[1][fac].x)
						push!(aprem,fac)
					elseif has(term,ap[1][fac])
						ti=indsin(term,ap[1][fac])
						fi=0
						for tfi in ti
							if !in(tfi,termrem)
								fi=tfi
								break
							end
						end
						if fi!=0
							push!(termrem,fi)
							push!(aprem,fac)
						end
					end
				end
				deleteat!(ap[1],aprem)
				term[i]=Div(expression(ap))
				pushallunique!(remove,termrem)
			end
		else
			invinds=indsin(term,term[i].x)
			removed=findin(invinds,remove)
			deleteat!(invinds,removed)
			if !isempty(invinds)
				push!(remove,i)
				push!(remove,invinds[end])
			end
		end
	end
	if !isempty(remove)
		deleteat!(term,sort!(remove))
		if isempty(term)
			push!(term,1)
		end
		return term
	else
		return term
	end
end
divify(term::Array)=divify!(deepcopy(term))
divify(x::X)=x
function simplify(ex::Expression,t::Type{Div})
	ap=dcterms(ex)
	for term in 1:length(ap)
		ap[term]=divify!(ap[term]) #this should be phased out
		ap[term]=divbine!(ap[term])
		ap[term]=divbinedify!(ap[term]) #this needn't be called if divify succeeded however it will fail if divbine was called before and in that case divbine doesn't need be called here...
	end
	return simplify(expression(ap))
end
function simplify(d::Div)
	x=simplify(d.x)
	if isa(x,Number)
		return 1/x
	end
	return Div(x)
end
function matches(eq::Equation,t::Type{Div})
	lhs=terms(eq.lhs)
	rhs=terms(eq.rhs)
	m=Equation[]
	for term in lhs
		for fac in term
			nl=deepcopy(lhs)
			nr=deepcopy(rhs)
			for t in nl
				push!(t,Div(deepcopy(fac)))
			end
			for t in nr
				push!(t,Div(deepcopy(fac)))
			end
			nl=simplify(expression(nl))
			nr=simplify(expression(nr))
			teq=Equation(nl,nr,Any[fac])
			if !(teq∈m)&&!(isa(nl,Number)&&isa(nr,Number)&&nl!=nr) #prevents 1=0 from x=0
				push!(m,teq)
			end
		end
	end
	return m
end

ext: jl


text:
struct InformationFilter <: Updater
    up::Updater
    im::InformationMeasure
    initial_im::InformationMeasure
end
InformationFilter(up::Updater, im::InformationMeasure) = InformationFilter(up, im, im)

function POMDPs.update(up::InformationFilter, ib::InformationBelief{B}, a, o) where {B <: ParticleCollection}
    @debug("Start entropy filter update")
    b = ib.b
    bp = update(up.up, b, a, o)
    wbp = weighted_posterior_belief(up.up)
    bpred = predicted_belief(up.up)
    i = information(up.im, bp=bp, wbp=wbp, b=b, bpred=bpred, a=a, o=o)
    InformationBelief(bp, i)
end

function POMDPs.initialize_belief(up::InformationFilter, d::Any)
    b = initialize_belief(up.up, d)
    i = information(up.initial_im, bp=b)
    InformationBelief(b, i)
end


text:
struct Parameters
    ν::Float64
    σ::Float64
    tspan::Tuple{Float64,Float64}
end


text:
###########################  Blotter

function show(io::IO, ta::Blotter)
  # variables
  nrow          = size(ta.values, 1)
  ncol          = size(ta.values, 2)
  intcatcher    = falses(ncol)
  for c in 1:ncol
      rowcheck =  trunc(ta.values[:,c]) - ta.values[:,c] .== 0
      if sum(rowcheck) == length(rowcheck)
          intcatcher[c] = true
      end
  end
  spacetime     = strwidth(string(ta.timestamp[1])) + 3
  firstcolwidth = strwidth(ta.colnames[1])
  colwidth      = Int[]
      for m in 1:ncol
          push!(colwidth, max(strwidth(ta.colnames[m]), strwidth(@sprintf("%.2f", maximum(ta.values[:,m])))))
      end

  # summary line
  println(io,"")
  #print(io,@sprintf("%dx%d %s %s to %s", nrow, ncol, typeof(ta), string(ta.timestamp[1]), string(ta.timestamp[nrow])))
  print(io,@sprintf("%dx%d Financial Blotter for Ticker *** %s ***, %s to %s",
                    nrow, ncol, string(ta.meta), string(ta.timestamp[1]), string(ta.timestamp[nrow])))
  println(io,"")

  # row label line
   print_with_color(:magenta, io, ^(" ", spacetime), ta.colnames[1], ^(" ", colwidth[1] + 2 -firstcolwidth))

   for p in 2:length(colwidth)
     print_with_color(:magenta, io, ta.colnames[p], ^(" ", colwidth[p] - strwidth(ta.colnames[p]) + 2))
   end
   println(io,"")

  # timestamp and values line
    if nrow > 7
        for i in 1:4
            print(io, ta.timestamp[i], " | ")
        for j in 1:ncol
            intcatcher[j] && ta.values[i,j] > 0?
            print_with_color(:green, io, rpad(iround(ta.values[i,j]), colwidth[j] + 2, " ")) :
            intcatcher[j] && ta.values[i,j] < 0 ?
            print_with_color(:red, io, rpad(iround(ta.values[i,j]), colwidth[j] + 2, " ")) :
            print_with_color(:blue, io, rpad(round(ta.values[i,j], 2), colwidth[j] + 2, " "))
        end
        println(io,"")
        end
        println(io,'\u22EE')
        for i in nrow-3:nrow
            print(io, ta.timestamp[i], " | ")
        for j in 1:ncol
            intcatcher[j] && ta.values[i,j] > 0?
            print_with_color(:green, io, rpad(iround(ta.values[i,j]), colwidth[j] + 2, " ")) :
            intcatcher[j] && ta.values[i,j] < 0 ?
            print_with_color(:red, io, rpad(iround(ta.values[i,j]), colwidth[j] + 2, " ")) :
            print_with_color(:blue, io, rpad(round(ta.values[i,j], 2), colwidth[j] + 2, " "))
        end
        println(io,"")
        end
    else
        for i in 1:nrow
            print(io, ta.timestamp[i], " | ")
        for j in 1:ncol
            intcatcher[j] ?
            print(io, rpad(iround(ta.values[i,j]), colwidth[j] + 2, " ")) :
            print(io, rpad(round(ta.values[i,j], 2), colwidth[j] + 2, " "))
        end
        println(io,"")
        end
    end
end

########################### Trade

# function show(io::IO, t::Trade)
#     print_with_color(:blue, io, @sprintf("%s:", string(t.start)))
#     print_with_color(:blue, io, @sprintf("%s ", string(t.finish)))
#     print_with_color(:yellow, io, @sprintf("%s ~ ", string(t.open)))
#     print_with_color(:yellow, io, @sprintf("%s ", string(t.close)))
#     if t.side == "long" && t.close - t.open> 0
#         print_with_color(:green, io, @sprintf("qty = "))
#     else
#         print_with_color(:red, io, @sprintf("qty = "))
#     end
#     print_with_color(:blue, io, @sprintf("%s", string(t.quantity)))
# end

# This code needs to be reduced by several fold, storing it here for now

###########################  Account

# function show(io::IO, ta::Account)
#   # variables
#   nrow          = size(ta.values, 1)
#   ncol          = size(ta.values, 2)
#   intcatcher    = falses(ncol)
#   for c in 1:ncol
#       rowcheck =  trunc(ta.values[:,c]) - ta.values[:,c] .== 0
#       if sum(rowcheck) == length(rowcheck)
#           intcatcher[c] = true
#       end
#   end
#   spacetime     = strwidth(string(ta.timestamp[1])) + 3
#   firstcolwidth = strwidth(ta.colnames[1])
#   colwidth      = Int[]
#       for m in 1:ncol
#           push!(colwidth, max(strwidth(ta.colnames[m]), strwidth(@sprintf("%.2f", maximum(ta.values[:,m])))))
#       end

#   # summary line
#   print(io,@sprintf("%dx%d %s %s to %s", nrow, ncol, typeof(ta), string(ta.timestamp[1]), string(ta.timestamp[nrow])))
#   println(io,"")
#   println(io,"")

#   # row label line
#    print_with_color(:magenta, io, ^(" ", spacetime), ta.colnames[1], ^(" ", colwidth[1] + 2 -firstcolwidth))

#    for p in 2:length(colwidth)
#      print_with_color(:magenta, io, ta.colnames[p], ^(" ", colwidth[p] - strwidth(ta.colnames[p]) + 2))
#    end
#    println(io,"")

#   # timestamp and values line
#     if nrow > 7
#         for i in 1:4
#             print(io, ta.timestamp[i], " | ")
#         for j in 1:ncol
#             intcatcher[j] && ta.values[i,j] > 0?
#             print_with_color(:green, io, rpad(iround(ta.values[i,j]), colwidth[j] + 2, " ")) :
#             intcatcher[j] && ta.values[i,j] < 0 ?
#             print_with_color(:red, io, rpad(iround(ta.values[i,j]), colwidth[j] + 2, " ")) :
#             print_with_color(:blue, io, rpad(round(ta.values[i,j], 2), colwidth[j] + 2, " "))
#         end
#         println(io,"")
#         end
#         println(io,'\u22EE')
#         for i in nrow-3:nrow
#             print(io, ta.timestamp[i], " | ")
#         for j in 1:ncol
#             intcatcher[j] && ta.values[i,j] > 0?
#             print_with_color(:green, io, rpad(iround(ta.values[i,j]), colwidth[j] + 2, " ")) :
#             intcatcher[j] && ta.values[i,j] < 0 ?
#             print_with_color(:red, io, rpad(iround(ta.values[i,j]), colwidth[j] + 2, " ")) :
#             print_with_color(:blue, io, rpad(round(ta.values[i,j], 2), colwidth[j] + 2, " "))
#         end
#         println(io,"")
#         end
#     else
#         for i in 1:nrow
#             print(io, ta.timestamp[i], " | ")
#         for j in 1:ncol
#             intcatcher[j] ?
#             print(io, rpad(iround(ta.values[i,j]), colwidth[j] + 2, " ")) :
#             print(io, rpad(round(ta.values[i,j], 2), colwidth[j] + 2, " "))
#         end
#         println(io,"")
#         end
#     end
# end

###########################  Portfolio

# function show(io::IO, ta::Portfolio)
#   # variables
#   nrow          = size(ta.values, 1)
#   ncol          = size(ta.values, 2)
#   intcatcher    = falses(ncol)
#   for c in 1:ncol
#       rowcheck =  trunc(ta.values[:,c]) - ta.values[:,c] .== 0
#       if sum(rowcheck) == length(rowcheck)
#           intcatcher[c] = true
#       end
#   end
#   spacetime     = strwidth(string(ta.timestamp[1])) + 3
#   firstcolwidth = strwidth(ta.colnames[1])
#   colwidth      = Int[]
#       for m in 1:ncol
#           push!(colwidth, max(strwidth(ta.colnames[m]), strwidth(@sprintf("%.2f", maximum(ta.values[:,m])))))
#       end

#   # summary line
#   print(io,@sprintf("%dx%d %s %s to %s", nrow, ncol, typeof(ta), string(ta.timestamp[1]), string(ta.timestamp[nrow])))
#   println(io,"")
#   println(io,"")

#   # row label line
#    print_with_color(:magenta, io, ^(" ", spacetime), ta.colnames[1], ^(" ", colwidth[1] + 2 -firstcolwidth))

#    for p in 2:length(colwidth)
#      print_with_color(:magenta, io, ta.colnames[p], ^(" ", colwidth[p] - strwidth(ta.colnames[p]) + 2))
#    end
#    println(io,"")

#   # timestamp and values line
#     if nrow > 7
#         for i in 1:4
#             print(io, ta.timestamp[i], " | ")
#         for j in 1:ncol
#             intcatcher[j] && ta.values[i,j] > 0?
#             print_with_color(:green, io, rpad(iround(ta.values[i,j]), colwidth[j] + 2, " ")) :
#             intcatcher[j] && ta.values[i,j] < 0 ?
#             print_with_color(:red, io, rpad(iround(ta.values[i,j]), colwidth[j] + 2, " ")) :
#             print_with_color(:blue, io, rpad(round(ta.values[i,j], 2), colwidth[j] + 2, " "))
#         end
#         println(io,"")
#         end
#         println(io,'\u22EE')
#         for i in nrow-3:nrow
#             print(io, ta.timestamp[i], " | ")
#         for j in 1:ncol
#             intcatcher[j] && ta.values[i,j] > 0?
#             print_with_color(:green, io, rpad(iround(ta.values[i,j]), colwidth[j] + 2, " ")) :
#             intcatcher[j] && ta.values[i,j] < 0 ?
#             print_with_color(:red, io, rpad(iround(ta.values[i,j]), colwidth[j] + 2, " ")) :
#             print_with_color(:blue, io, rpad(round(ta.values[i,j], 2), colwidth[j] + 2, " "))
#         end
#         println(io,"")
#         end
#     else
#         for i in 1:nrow
#             print(io, ta.timestamp[i], " | ")
#         for j in 1:ncol
#             intcatcher[j] ?
#             print(io, rpad(iround(ta.values[i,j]), colwidth[j] + 2, " ")) :
#             print(io, rpad(round(ta.values[i,j], 2), colwidth[j] + 2, " "))
#         end
#         println(io,"")
#         end
#     end
# end


text:
# to get rid of eventually
const Columns = StructVector

# IndexedTable-like API

"""
    colnames(itr)

Returns the names of the "columns" in `itr`.

# Examples:

    colnames(1:3)
    colnames(Columns([1,2,3], [3,4,5]))
    colnames(table([1,2,3], [3,4,5]))
    colnames(Columns(x=[1,2,3], y=[3,4,5]))
    colnames(table([1,2,3], [3,4,5], names=[:x,:y]))
    colnames(ndsparse(Columns(x=[1,2,3]), Columns(y=[3,4,5])))
    colnames(ndsparse(Columns(x=[1,2,3]), [3,4,5]))
    colnames(ndsparse(Columns(x=[1,2,3]), [3,4,5]))
    colnames(ndsparse(Columns([1,2,3], [4,5,6]), Columns(x=[6,7,8])))
    colnames(ndsparse(Columns(x=[1,2,3]), Columns([3,4,5],[6,7,8])))

"""
function colnames end

Base.@pure colnames(t::AbstractVector) = (1,)
columns(v::AbstractVector) = v

Base.@pure colnames(t::Columns) = fieldnames(eltype(t))
Base.@pure colnames(t::Columns{<:Pair}) = colnames(t.first) => colnames(t.second)

"""
    columns(itr, select::Selection = All())

Select one or more columns from an iterable of rows as a tuple of vectors.

`select` specifies which columns to select. Refer to the [`select`](@ref) function for the
available selection options and syntax.

`itr` can be `NDSparse`, `Columns`, `AbstractVector`, or their distributed counterparts.

# Examples

    t = table(1:2, 3:4; names = [:x, :y])

    columns(t)
    columns(t, :x)
    columns(t, (:x,))
    columns(t, (:y, :x => -))
"""
function columns end

columns(c::Columns) = fieldarrays(c)
columns(c::Columns{<:Tuple}) = Tuple(fieldarrays(c))
columns(c::Columns{<:Pair}) = c.first => c.second

"""
    ncols(itr)

Returns the number of columns in `itr`.

# Examples

    ncols([1,2,3]) == 1
    ncols(rows(([1,2,3],[4,5,6]))) == 2
"""
function ncols end
ncols(c::Columns{T, C}) where {T, C} = fieldcount(C)
ncols(c::Columns{<:Pair}) = ncols(c.first) => ncols(c.second)
ncols(c::AbstractArray) = 1

summary(c::Columns{D}) where {D<:Tuple} = "$(length(c))-element Columns{$D}"

_sizehint!(c::Columns, n::Integer) = (foreachfield(x->_sizehint!(x,n), c); c)

function _strip_pair(c::Columns{<:Pair})
    f, s = map(columns, fieldarrays(c))
    (f isa AbstractVector) && (f = (f,))
    (s isa AbstractVector) && (s = (s,))
    Columns((f..., s...))
end

# fused indexing operations
# these can be implemented for custom vector types like PooledVector where
# you can get big speedups by doing indexing and an operation in one step.

@inline copyelt!(a, i, j) = (@inbounds a[i] = a[j])
@inline copyelt!(a, i, b, j) = (@inbounds a[i] = b[j])
@inline copyelt!(a::PooledArray, i, j) = (a.refs[i] = a.refs[j])

# row operations

@inline roweq(x::AbstractVector, i, j) = (@inbounds eq=isequal(x[i], x[j]); eq)
@inline roweq(a::PooledArray, i, j) = (@inbounds x=a.refs[i] == a.refs[j]; x)
@inline function roweq(a::StringArray{String}, i, j)
    weaksa = convert(StringArray{WeakRefString{UInt8}}, a)
    @inbounds isequal(weaksa[i], weaksa[j])
end

copyrow!(I::Columns, i, src) = foreachfield(c->copyelt!(c, i, src), I)
copyrow!(I::Columns, i, src::Columns, j) = foreachfield((c1,c2)->copyelt!(c1, i, c2, j), I, src)
copyrow!(I::AbstractArray, i, src::AbstractArray, j) = (@inbounds I[i] = src[j])
pushrow!(to::Columns, from::Columns, i) = foreachfield((a,b)->push!(a, b[i]), to, from)
pushrow!(to::AbstractArray, from::AbstractArray, i) = push!(to, from[i])

@generated function roweq(c::Columns{D,C}, i, j) where {D,C}
    N = fieldcount(C)
    ex = :(roweq(getfield(fieldarrays(c),1), i, j))
    for n in 2:N
        ex = :(($ex) && (roweq(getfield(fieldarrays(c),$n), i, j)))
    end
    ex
end

# uses number of columns from `d`, assuming `c` has more or equal
# dimensions, for broadcast joins.
@generated function rowcmp(c::Columns, i, d::Columns{D}, j) where D
    N = fieldcount(D)
    ex = :(cmp(getfield(fieldarrays(c),$N)[i], getfield(fieldarrays(d),$N)[j]))
    for n in N-1:-1:1
        ex = quote
            let k = rowcmp(getfield(fieldarrays(c),$n), i, getfield(fieldarrays(d),$n), j)
                (k == 0) ? ($ex) : k
            end
        end
    end
    ex
end

@inline function rowcmp(c::AbstractVector, i, d::AbstractVector, j)
    cmp(c[i], d[j])
end

@inline function rowcmp(c::StringArray{String}, i, d::StringArray{String}, j)
    wc = convert(StringArray{WeakRefString{UInt8}}, c)
    wd = convert(StringArray{WeakRefString{UInt8}}, d)
    cmp(wc[i], wd[j])
end

# test that the row on the right is "as of" the row on the left, i.e.
# all columns are equal except left >= right in last column.
# Could be generalized to some number of trailing columns, but I don't
# know whether that has applications.
@generated function row_asof(c::Columns{D,C}, i, d::Columns{D,C}, j) where {D,C}
    N = length(C.parameters)
    if N == 1
        ex = :(!isless(getfield(fieldarrays(c),1)[i], getfield(fieldarrays(d),1)[j]))
    else
        ex = :(isequal(getfield(fieldarrays(c),1)[i], getfield(fieldarrays(d),1)[j]))
    end
    for n in 2:N
        if N == n
            ex = :(($ex) && !isless(getfield(fieldarrays(c),$n)[i], getfield(fieldarrays(d),$n)[j]))
        else
            ex = :(($ex) && isequal(getfield(fieldarrays(c),$n)[i], getfield(fieldarrays(d),$n)[j]))
        end
    end
    ex
end

# map

"""
    map_rows(f, c...)

Transform collection `c` by applying `f` to each element. For multiple collection arguments, apply `f`
elementwise. Collect output as `Columns` if `f` returns
`Tuples` or `NamedTuples` with constant fields, as `Array` otherwise.

# Examples

    map_rows(i -> (exp = exp(i), log = log(i)), 1:5)
"""
function map_rows(f, iters...)
    collect_columns(f(i...) for i in zip(iters...))
end

# 1-arg case
map_rows(f, iter) = collect_columns(f(i) for i in iter)

## Special selectors to simplify column selector

"""
    All(cols::Union{Symbol, Int}...)

Select the union of the selections in `cols`. If `cols == ()`, select all columns.

# Examples

    t = table([1,1,2,2], [1,2,1,2], [1,2,3,4], [0, 0, 0, 0], names=[:a,:b,:c,:d])
    select(t, All(:a, (:b, :c)))
    select(t, All())
"""
struct All{T}
    cols::T
end

All(args...) = All(args)

"""
    Not(cols::Union{Symbol, Int}...)

Select the complementary of the selection in `cols`. `Not` can accept several arguments,
in which case it returns the complementary of the union of the selections.

# Examples

    t = table([1,1,2,2], [1,2,1,2], [1,2,3,4], names=[:a,:b,:c], pkey = (:a, :b))
    select(t, Not(:a))
    select(t, Not(:a, (:a, :b)))
"""
struct Not{T}
    cols::T
end

Not(args...) = Not(All(args))

"""
    Keys()

Select the primary keys.

# Examples

    t = table([1,1,2,2], [1,2,1,2], [1,2,3,4], names=[:a,:b,:c], pkey = (:a, :b))
    select(t, Keys())
"""
struct Keys; end

"""
    Between(first, last)

Select the columns between `first` and `last`.

# Examples

    t = table([1,1,2,2], [1,2,1,2], 1:4, 'a':'d', names=[:a,:b,:c,:d])
    select(t, Between(:b, :d))
"""
struct Between{T1 <: Union{Int, Symbol}, T2 <: Union{Int, Symbol}}
    first::T1
    last::T2
end

const SpecialSelector = Union{Not, All, Keys, Between, Function, Regex, Type}

hascolumns(t, s) = true
hascolumns(t, s::Symbol) = s in colnames(t)
hascolumns(t, s::Int) = s in 1:length(columns(t))
hascolumns(t, s::Tuple) = all(hascolumns(t, x) for x in s)
hascolumns(t, s::Not) = hascolumns(t, s.cols)
hascolumns(t, s::Between) = hascolumns(t, s.first) && hascolumns(t, s.last)
hascolumns(t, s::All) = all(hascolumns(t, x) for x in s.cols)
hascolumns(t, s::Type) = any(x -> eltype(x) <: s, columns(t))

lowerselection(t, s)                     = s
lowerselection(t, s::Union{Int, Symbol}) = colindex(t, s)
lowerselection(t, s::Tuple)              = map(x -> lowerselection(t, x), s)
lowerselection(t, s::Not)                = excludecols(t, lowerselection(t, s.cols))
lowerselection(t, s::Keys)               = lowerselection(t, IndexedTables.pkeynames(t))
lowerselection(t, s::Between)            = Tuple(colindex(t, s.first):colindex(t, s.last))
lowerselection(t, s::Function)           = colindex(t, Tuple(filter(s, collect(colnames(t)))))
lowerselection(t, s::Regex)              = lowerselection(t, x -> occursin(s, string(x)))
lowerselection(t, s::Type)               = Tuple(findall(x -> eltype(x) <: s, columns(t)))

function lowerselection(t, s::All)
    s.cols == () && return lowerselection(t, valuenames(t))
    ls = (isa(i, Tuple) ? i : (i,) for i in lowerselection(t, s.cols))
    ls |> Iterators.flatten |> union |> Tuple
end

### Iteration API

# For `columns(t, names)` and `rows(t, ...)` to work, `t`
# needs to support `colnames` and `columns(t)`

Base.@pure function colindex(t, col::Tuple)
    fns = colnames(t)
    map(x -> _colindex(fns, x), col)
end

Base.@pure function colindex(t, col)
    _colindex(colnames(t), col)
end

function colindex(t, col::SpecialSelector)
    colindex(t, lowerselection(t, col))
end

function _colindex(fnames::Union{Tuple, AbstractArray}, col, default=nothing)
    if isa(col, Int) && 1 <= col <= length(fnames)
        return col
    elseif isa(col, Symbol)
        idx = something(findfirst(isequal(col), fnames), 0)
        idx > 0 && return idx
    elseif isa(col, Pair{<:Any, <:AbstractArray})
        return 0
    elseif isa(col, Tuple)
        return 0
    elseif isa(col, Pair{Symbol, <:Pair}) # recursive pairs
        return _colindex(fnames, col[2])
    elseif isa(col, Pair{<:Any, <:Any})
        return _colindex(fnames, col[1])
    elseif isa(col, AbstractArray)
        return 0
    end
    default !== nothing ? default : error("column $col not found.")
end

# const ColPicker = Union{Int, Symbol, Pair{Symbol=>Function}, Pair{Symbol=>AbstractVector}, AbstractVector}
column(c, x) = columns(c)[colindex(c, x)]

# optimized method
@inline function column(c::Columns, x::Union{Int, Symbol})
    getfield(fieldarrays(c), x)
end

column(t, a::AbstractArray) = a
column(t, a::Pair{Symbol, <:AbstractArray}) = column(t, a[2])
column(t, a::Pair{Symbol, <:Pair}) = rows(t, a[2]) # renaming a selection
column(t, a::Pair{<:Any, <:Any}) = map(a[2], rows(t, a[1]))
column(t, s::SpecialSelector) = rows(t, lowerselection(t, s))

function columns(c, sel::Union{Tuple, SpecialSelector})
    which = lowerselection(c, sel)
    cnames = colnames(c, which)
    if all(x->isa(x, Symbol), cnames)
        tuplewrap = namedtuple(cnames...)∘tuple
    else
        tuplewrap = tuple
    end
    tuplewrap((rows(c, w) for w in which)...)
end

"""
`columns(itr, which)`

Returns a vector or a tuple of vectors from the iterator.

"""
columns(t, which) = column(t, which)

function colnames(c, cols::Union{Tuple, AbstractArray})
    map(x->colname(c, x), cols)
end

colnames(c, cols::SpecialSelector) = colnames(c, lowerselection(c, cols))

function colname(c, col)
    if isa(col, Union{Int, Symbol})
        col == 0 && return 0
        i = colindex(c, col)
        return colnames(c)[i]
    elseif isa(col, Pair{<:Any, <:Any})
        return col[1]
    elseif isa(col, Tuple)
        #ns = map(x->colname(c, x), col)
        return 0
    elseif isa(col, SpecialSelector)
        return 0
    elseif isa(col, AbstractVector)
        return 0
    end
    error("column named $col not found")
end

"""
    rows(itr, select = All())

Select one or more fields from an iterable of rows as a vector of their values.  Refer to
the [`select`](@ref) function for selection options and syntax.

`itr` can be [`NDSparse`](@ref), `StructArrays.StructVector`, `AbstractVector`, or their distributed counterparts.

# Examples

    t = table([1,2],[3,4], names=[:x,:y])
    rows(t)
    rows(t, :x)
    rows(t, (:x,))
    rows(t, (:y, :x => -))
"""
function rows end

rows(x::AbstractVector) = x
rows(cols::Tup) = Columns(cols)

rows(t, which...) = rows(columns(t, which...))

_cols_tuple(xs::Columns) = columns(xs)
_cols_tuple(xs::AbstractArray) = (xs,)
concat_cols(xs, ys) = rows(concat_tup(_cols_tuple(xs), _cols_tuple(ys)))

## Mutable Columns Dictionary

mutable struct ColDict{T}
    pkey::Vector{Int}
    src::T
    names::Vector
    columns::Vector
    copy::Union{Nothing, Bool}
end

"""
    d = ColDict(t)

Create a mutable dictionary of columns in `t`.

To get the immutable iterator of the same type as `t`
call `d[]`
"""
function ColDict(t; copy=nothing)
    cnames = colnames(t)
    if cnames isa AbstractArray
        cnames = Base.copy(cnames)
    end
    ColDict(Int[], t, convert(Array{Any}, collect(cnames)), Any[columns(t)...], copy)
end

Base.keys(d::ColDict) = d.names
Base.values(d::ColDict) = d.columns

function Base.getindex(d::ColDict{<:Columns})
    Columns(Tuple(d.columns); names=d.names)
end

Base.getindex(d::ColDict, key) = rows(d[], key)
Base.getindex(d::ColDict, key::AbstractArray) = key

function Base.setindex!(d::ColDict, x, key::Union{Symbol, Int})
    k = _colindex(d.names, key, 0)
    col = d[x]
    if k == 0
        push!(d.names, key)
        push!(d.columns, col)
    elseif k in d.pkey
        # primary key column has been modified.
        # copy the table as this results in a shuffle
        if d.copy === nothing
            d.copy = true
        end
        d.columns[k] = col
    else
        d.columns[k] = col
    end
end

set!(d::ColDict, key::Union{Symbol, Int}, x) = setindex!(d, x, key)

function Base.haskey(d::ColDict, key)
    _colindex(d.names, key, 0) != 0
end

function Base.insert!(d::ColDict, index, key, col)
    if haskey(d, key)
        error("Key $key already exists. Use dict[key] = col instead of inserting.")
    else
        insert!(d.names, index, key)
        insert!(d.columns, index, rows(d.src, col))
        for (i, pk) in enumerate(d.pkey)
            if pk >= index
                d.pkey[i] += 1 # moved right
            end
        end
    end
end

function insertafter!(d::ColDict, i, key, col)
    k = _colindex(d.names, i, 0)
    if k == 0
        error("$i not found. Cannot insert column after $i")
    end
    insert!(d, k+1, key, col)
end

function insertbefore!(d::ColDict, i, key, col)
    k = _colindex(d.names, i, 0)
    if k == 0
        error("$i not found. Cannot insert column after $i")
    end
    insert!(d, k, key, col)
end

function Base.pop!(d::ColDict, key::Union{Symbol, Int}=length(d.names))
    k = _colindex(d.names, key, 0)
    local col
    if k == 0
        error("Column $key not found")
    else
        col = d.columns[k]
        deleteat!(d.names, k)
        deleteat!(d.columns, k)
        idx = [pk[1] for pk in enumerate(d.pkey) if pk[2] == k]
        deleteat!(d.pkey, idx)
        for i in 1:length(d.pkey)
            if d.pkey[i] > k
                d.pkey[i] -= 1
            end
        end
        if !isempty(idx) && d.copy === nothing
            # set copy to true
            d.copy = true
        end
    end
    col
end

function rename!(d::ColDict, col::Union{Symbol, Int}, newname)
    k = _colindex(d.names, col, 0)
    if k == 0
        error("$col not found. Cannot rename it.")
    end
    d.names[k] = newname
end

Base.push!(d::ColDict, key::AbstractString, x) = push!(d, Symbol(key), x)
function Base.push!(d::ColDict, key::Union{Symbol, Int}, x)
    push!(d.names, key)
    push!(d.columns, rows(d.src, x))
end

for s in [:(Base.pop!), :(Base.push!), :(rename!), :(set!)]
    if s == :(Base.pop!)
        typ = :(Union{Symbol, Int})
    else
        typ = :Pair
        @eval $s(t::ColDict, x::Pair) = $s(t, x.first, x.second)
    end
    @eval begin
        function $s(t::ColDict, args)
            for i in args
                $s(t, i)
            end
        end
        $s(t::ColDict, args::Vararg{$typ}) = $s(t, args)
    end
end

function _cols(expr)
    if expr.head == :call
        dict = :(dict = ColDict($(expr.args[2])))
        expr.args[2] = :dict
        quote
            let $dict
                $expr
                dict[]
            end
        end |> esc
    else
        error("This form of @cols is not implemented. Use `@cols f(t,args...)` where `t` is the collection.")
    end
end

macro cols(expr)
    _cols(expr)
end

# Modifying a columns

"""
    setcol(t::Table, col::Union{Symbol, Int}, x::Selection)

Sets a `x` as the column identified by `col`. Returns a new table.

    setcol(t::Table, map::Pair{}...)

Set many columns at a time.

# Examples:

    t = table([1,2], [3,4], names=[:x, :y])

    # change second column to [5,6]
    setcol(t, 2 => [5,6])
    setcol(t, :y , :y => x -> x + 2)

    # add [5,6] as column :z
    setcol(t, :z => 5:6)
    setcol(t, :z, :y => x -> x + 2)

    # replacing the primary key results in a re-sorted copy
    t = table([0.01, 0.05], [1,2], [3,4], names=[:t, :x, :y], pkey=:t)
    t2 = setcol(t, :t, [0.1,0.05])
"""
setcol(t, args...) = @cols set!(t, args...)

"""
    pushcol(t, name, x)

Push a column `x` to the end of the table. `name` is the name for the new column. Returns a new table.

    pushcol(t, map::Pair...)

Push many columns at a time.

# Example

    t = table([0.01, 0.05], [2,1], [3,4], names=[:t, :x, :y], pkey=:t)
    pushcol(t, :z, [1//2, 3//4])
    pushcol(t, :z => [1//2, 3//4])
"""
pushcol(t, args...) = @cols push!(t, args...)

"""
    popcol(t, cols...)

Remove the column(s) `cols` from the table. Returns a new table.

# Example

    t = table([0.01, 0.05], [2,1], [3,4], names=[:t, :x, :y], pkey=:t)
    popcol(t, :x)
"""
popcol(t, args...) = @cols pop!(t, args...)

"""
    insertcol(t, position::Integer, name, x)

Insert a column `x` named `name` at `position`. Returns a new table.

# Example

    t = table([0.01, 0.05], [2,1], [3,4], names=[:t, :x, :y], pkey=:t)
    insertcol(t, 2, :w, [0,1])
"""
insertcol(t, i::Integer, name, x) = @cols insert!(t, i, name, x)

"""
    insertcolafter(t, after, name, col)

Insert a column `col` named `name` after `after`. Returns a new table.

# Example

    t = table([0.01, 0.05], [2,1], [3,4], names=[:t, :x, :y], pkey=:t)
    insertcolafter(t, :t, :w, [0,1])
"""
insertcolafter(t, after, name, x) = @cols insertafter!(t, after, name, x)

"""
    insertcolbefore(t, before, name, col)

Insert a column `col` named `name` before `before`. Returns a new table.

# Example

    t = table([0.01, 0.05], [2,1], [3,4], names=[:t, :x, :y], pkey=:t)
    insertcolbefore(t, :x, :w, [0,1])
"""
insertcolbefore(t, before, name, x) = @cols insertbefore!(t, before, name, x)

"""
    renamecol(t, col, newname)

Set `newname` as the new name for column `col` in `t`. Returns a new table.

    renamecol(t, map::Pair...)

Rename multiple columns at a time.

# Example

    t = table([0.01, 0.05], [2,1], names=[:t, :x])
    renamecol(t, :t, :time)
"""
renamecol(t, args...) = @cols rename!(t, args...)

## Utilities for mapping and reduction with many functions / OnlineStats

@inline _apply(f::OnlineStat, g, x) = (fit!(g, x); g)
@inline _apply(f::Tup, y::Tup, x::Tup) = _apply(astuple(f), astuple(y), astuple(x))
@inline _apply(f::Tuple, y::Tuple, x::Tuple) = map(_apply, f, y, x)
@inline _apply(f::NamedTuple, y::NamedTuple, x::NamedTuple) = map(_apply, f, y, x)
@inline _apply(f, y, x) = f(y, x)
@inline _apply(f::Tup, x::Tup) = _apply(astuple(f), astuple(x))
@inline _apply(f::NamedTuple, x::NamedTuple) = map(_apply, f, x)
@inline _apply(f::Tuple, x::Tuple) = map(_apply, f, x)
@inline _apply(f, x) = f(x)

@inline init_first(f, x) = x
@inline init_first(f::OnlineStat, x) = (g=copy(f); fit!(g, x); g)
@inline init_first(f::Tup, x::Tup) = map(init_first, f, x)

# Initialize type of output, functions to apply, input and output vectors

function reduced_type(f, x, isvec, key = nothing)
    if key !== nothing
        _promote_op(f, eltype(key), typeof(x))
    elseif isvec
        _promote_op(f, typeof(x))
    else
        _promote_op((a,b)->_apply(f, init_first(f, a), b),
                    eltype(x), eltype(x))
    end
end

function init_inputs(f, x, gettype, isvec) # normal functions
    f, x, gettype(f, x, isvec)
end

nicename(f::Function) = typeof(f).name.mt.name
nicename(f) = Symbol(last(split(string(f), ".")))
nicename(o::OnlineStat) = Symbol(typeof(o).name.name)

function mapped_type(f, x, isvec)
    _promote_op(f, eltype(x))
end

init_funcs(f, isvec) = init_funcs((f,), isvec)

function init_funcs(f::Tup, isvec)
    if isa(f, NamedTuple)
        return init_funcs((map(Pair, fieldnames(typeof(f)), f)...,), isvec)
    end

    funcmap = map(f) do g
        if isa(g, Pair)
            name = g[1]
            if isa(g[2], Pair)
                sel, fn = g[2]
            else
                sel = nothing
                fn = g[2]
            end
            (name, sel, fn)
        else
            (nicename(g), nothing, g)
        end
    end

    ns = map(x->x[1], funcmap)
    ss = map(x->x[2], funcmap)
    fs = map(map(x->x[3], funcmap)) do f
        f
    end

    NamedTuple{(ns...,)}((fs...,)), ss
end

function init_inputs(f::Tup, input, gettype, isvec)
    if isa(f, NamedTuple)
        return init_inputs((map(Pair, fieldnames(typeof(f)), f)...,), input, gettype, isvec)
    end
    fs, selectors = init_funcs(f, isvec)

    xs = map(s->s === nothing ? input : rows(input, s), selectors)

    output_eltypes = map((f,x) -> gettype(f, x, isvec), fs, xs)

    ns = fieldnames(typeof(fs))
    NT = namedtuple(ns...)

    # functions, input, output_eltype
    NT((fs...,)), rows(NT((xs...,))), NT{Tuple{output_eltypes...}}
end

### utils

compact_mem(x::Columns) = Columns(map(compact_mem, columns(x)))


text:
export transform_solutions

_parse_expression(exp) = exp isa String ? Num(eval(Meta.parse(exp))) : exp

"""
$(TYPEDSIGNATURES)

Takes a `Result` object and a string `f` representing a Symbolics.jl expression.
Returns an array with the values of `f` evaluated for the respective solutions.
Additional substitution rules can be specified in `rules` in the format `("a" => val)` or `(a => val)`
"""
function transform_solutions(res::Result, f::String; rules=Dict())
    # a string is used as input - a macro would not "see" the user's namespace while the user's namespace does not "see" the variables
    transformed = [Vector{ComplexF64}(undef, length(res.solutions[1])) for k in res.solutions] # preallocate

    # define variables in rules in this namespace
    new_keys = declare_variable.(string.(keys(Dict(rules))))
    expr = f isa String ? _parse_expression(f) : f

    fixed_subs = merge(res.fixed_parameters, Dict(zip(new_keys, values(Dict(rules)))))
    expr = substitute_all(expr, Dict(fixed_subs))

    vars = res.problem.variables
    all_symbols = cat(vars, collect(keys(res.swept_parameters)), dims=1)
    comp_func = build_function(expr, all_symbols)
    f = eval(comp_func)

    # preallocate an array for the numerical values, rewrite parts of it
    # when looping through the solutions
    vals = Vector{ComplexF64}(undef, length(all_symbols))
    n_vars = length(vars)
    n_pars = length(all_symbols) - n_vars

    for idx in CartesianIndices(res.solutions)
        params_values = res.swept_parameters[Tuple(idx)...]
        vals[end-n_pars+1:end] .= params_values # param values are common to all branches
        for (branch,soln) in enumerate(res.solutions[idx])
            vals[1:n_vars] .= soln
            transformed[idx][branch] = Base.invokelatest(f, vals)
        end
    end
    return transformed
end

# a simplified version meant to work with arrays of solutions
# cannot parse parameter values  mainly meant for time-dependent results
function transform_solutions(soln::Vector, f::String, harm_eq::HarmonicEquation)

    vars = _remove_brackets(get_variables(harm_eq))
    transformed = Vector{ComplexF64}(undef, length(soln))

    # parse the input with Symbolics
    expr = HarmonicBalance._parse_expression(f)

    rule(u) = Dict(zip(vars, u))

    transformed = map( x -> substitute_all(expr, rule(x)), soln)
    return convert(typeof(soln[1]), transformed)
end


text:
using Plots
include(joinpath(@__DIR__, "cartpole_friction.jl"))

# Nominal
X̄_nominal, Ū_nominal = unpack(Z̄_nominal, prob_nominal)
X̄_friction, Ū_friction = unpack(Z̄_friction, prob_friction)

# Plots results
S_friction_nominal = [Ū_friction[t][7] for t = 1:T-1]
@assert norm(S_friction_nominal, Inf) < 1.0e-4
b_friction_nominal = [(Ū_friction[t][2] - Ū_friction[t][3]) for t = 1:T-1]

t_nominal = range(0, stop = h * (T - 1), length = T)

# Control
plt = plot(t_nominal[1:T-1], hcat(Ū_nominal...)[1:1, :]',
    color = :cyan, width=2.0,
    title = "Cartpole", xlabel = "time (s)", ylabel = "control", label = "nominal",
    legend = :topright, linetype = :steppost)
plt = plot!(t_nominal[1:T-1], hcat(Ū_friction...)[1:1, :]', color = :magenta,
    width = 2.0, label = "nominal (friction)", linetype = :steppost)

# States
plt = plot(t_nominal, hcat(X̄_nominal...)[1:4, :]',
    color = :cyan, width = 2.0, xlabel = "time (s)",
    ylabel = "state", label = "", title = "Cartpole", legend = :topright)

plt = plot!(t_nominal, hcat(X̄_friction...)[1:4,:]',
    color = :magenta, width = 2.0, label = "")

# DPO
X̄_dpo, Ū_dpo = unpack(Z[prob_dpo.prob.idx.nom], prob_dpo.prob.prob.nom)
S_dpo = [Ū_dpo[t][7] for t = 1:T-1]
@assert norm(S_dpo, Inf) < 1.0e-4

X_sample = []
U_sample = []
for i = 1:N
	x, u = unpack(Z[prob_dpo.prob.idx.sample[i]],
		prob_dpo.prob.prob.sample[i])
	push!(X_sample, x)
	push!(U_sample, u)
end

# Control
plt = plot(t_nominal[1:T-1], hcat(Ū_dpo...)[1:1, :]',
    color = :orange, width=2.0,
    title = "Cartpole", xlabel = "time (s)", ylabel = "control", label = "nominal",
    legend = :topright, linetype = :steppost)

# States
plt = plot(t_nominal, hcat(X̄_nominal...)[1:4, :]',
    color = :orange, width = 2.0, xlabel = "time (s)",
    ylabel = "state", label = "", title = "Cartpole", legend = :topright)

# Simulation

plt_tvlqr_nom = plot(t_nominal,hcat(X_nominal...)[1:2,:]',legend=:topleft,color=:red,
    label=["nominal (no friction)" ""],
    width=2.0,xlabel="time",title="Cartpole",ylabel="state",ylims=(-1,5))
plt_tvlqr_nom = plot!(t_sim_nominal,hcat(z_tvlqr...)[1:2,:]',color=:purple,
    label=["tvlqr" ""],width=2.0)
savefig(plt_tvlqr_nom,joinpath(@__DIR__,"results/cartpole_friction_tvlqr_nom_sim.png"))

z_tvlqr_friction, u_tvlqr_friction, J_tvlqr_friction, Jx_tvlqr_friction, Ju_tvlqr_friction = simulate_cartpole_friction(K_friction_nominal,
    X_friction_nominal,U_friction_nominal,
    model_sim,Q_lqr,R_lqr,T_sim,Δt,X_friction_nominal[1],w,ul=ul_friction,uu=uu_friction,friction=true,
    μ=μ_sim)
plt_tvlqr_friction = plot(t_nominal,hcat(X_friction_nominal...)[1:2,:]',
    color=:red,label=["nominal (sample)" ""],
    width=2.0,xlabel="time",title="Cartpole",ylabel="state",legend=:topleft)
plt_tvlqr_friction = plot!(t_sim_nominal,hcat(z_tvlqr_friction...)[1:2,:]',
    color=:magenta,label=["tvlqr" ""],width=2.0)
savefig(plt_tvlqr_friction,joinpath(@__DIR__,"results/cartpole_friction_tvlqr_friction_sim.png"))

z_sample, u_sample, J_sample, Jx_sample, Ju_sample = simulate_cartpole_friction(K_sample,
    X_nom_sample,U_nom_sample,
    model_sim,Q_lqr,R_lqr,T_sim,Δt,X_nom_sample[1],w,ul=ul_friction,uu=uu_friction,friction=true,
    μ=μ_sim)

plt_sample = plot(t_sample,hcat(X_nom_sample...)[1:2,:]',legend=:bottom,color=:red,
    label=["nominal (sample)" ""],
    width=2.0,xlabel="time",title="Cartpole",ylabel="state")
plt_sample = plot!(t_sim_nominal,hcat(z_sample...)[1:2,:]',color=:orange,
    label=["sample" ""],width=2.0,legend=:topleft)
savefig(plt_sample,joinpath(@__DIR__,"results/cartpole_friction_sample_sim.png"))

# z_sample_c, u_sample_c, J_sample_c, Jx_sample_c, Ju_sample_c = simulate_cartpole_friction(K_sample_coefficients,
#     X_nom_sample_coefficients,U_nom_sample_coefficients,
#     model_sim,Q_lqr,R_lqr,T_sim,Δt,X_nom_sample_coefficients[1],w,ul=ul_friction,
#     uu=uu_friction,friction=true,
#     μ=μ_sim)

# plt_sample = plot(t_sample_c,hcat(X_nom_sample_coefficients...)[1:2,:]',
#     legend=:bottom,color=:red,label=["nominal (sample)" ""],
#     width=2.0,xlabel="time (s)",title="Cartpole",ylabel="state")
# plt_sample = plot!(t_sim_nominal,hcat(z_sample_c...)[1:2,:]',color=:magenta,
#     label=["sample" ""],width=2.0,legend=:topleft)
# savefig(plt_sample,joinpath(@__DIR__,"results/cartpole_friction_sample_sim.png"))

using PGFPlots
const PGF = PGFPlots

# nominal trajectory
psx_nom = PGF.Plots.Linear(t_nominal,hcat(X_nominal...)[1,:],mark="",
	style="color=cyan, line width=3pt")
psθ_nom = PGF.Plots.Linear(t_nominal,hcat(X_nominal...)[2,:],mark="",
	style="color=cyan, line width=3pt, densely dashed")

psx_nom_sim_tvlqr = PGF.Plots.Linear(t_sim_nominal,hcat(z_tvlqr...)[1,:],mark="",
	style="color=black, line width=2pt")
psθ_nom_sim_tvlqr = PGF.Plots.Linear(t_sim_nominal,hcat(z_tvlqr...)[2,:],mark="",
	style="color=black, line width=2pt, densely dashed")

a = Axis([psx_nom;psθ_nom;psx_nom_sim_tvlqr;psθ_nom_sim_tvlqr],
    xmin=0., ymin=-1, xmax=sum(H_nominal), ymax=3.5,
    axisEqualImage=false,
    hideAxis=false,
	ylabel="state",
	xlabel="time",
	legendStyle="{at={(0.01,0.99)},anchor=north west}",
	)

# Save to tikz format
dir = joinpath(@__DIR__,"results")
PGF.save(joinpath(dir,"cartpole_sim_tvlqr_no_friction.tikz"), a, include_preamble=false)

# nominal trajectory
psx_nom_friction = PGF.Plots.Linear(t_nominal,hcat(X_friction_nominal...)[1,:],mark="",
	style="color=purple, line width=3pt")
psθ_nom_friction = PGF.Plots.Linear(t_sample,hcat(X_friction_nominal...)[2,:],mark="",
	style="color=purple, line width=3pt, densely dashed")

psx_nom_sim_tvlqr_friction = PGF.Plots.Linear(t_sim_nominal,hcat(z_tvlqr_friction...)[1,:],mark="",
	style="color=black, line width=2pt",legendentry="pos.")
psθ_nom_sim_tvlqr_friction = PGF.Plots.Linear(t_sim_nominal,hcat(z_tvlqr_friction...)[2,:],mark="",
	style="color=black, line width=2pt, densely dashed",legendentry="ang.")

a_tvlqr_friction = Axis([psx_nom_friction;psθ_nom_friction;
    psx_nom_sim_tvlqr_friction;psθ_nom_sim_tvlqr_friction],
    xmin=0., ymin=-1, xmax=sum(H_nominal), ymax=3.5,
    axisEqualImage=false,
    hideAxis=false,
	ylabel="state",
	xlabel="time",
	legendStyle="{at={(0.01,0.99)},anchor=north west}",
	)

# Save to tikz format
dir = joinpath(@__DIR__,"results")
PGF.save(joinpath(dir,"cartpole_sim_tvlqr_friction.tikz"), a_tvlqr_friction, include_preamble=false)

# nominal trajectory
psx_dpo = PGF.Plots.Linear(t_sample,hcat(X_nom_sample...)[1,:],mark="",
	style="color=orange, line width=3pt")
psθ_dpo = PGF.Plots.Linear(t_sample,hcat(X_nom_sample...)[2,:],mark="",
	style="color=orange, line width=3pt, densely dashed")

psx_sim_dpo = PGF.Plots.Linear(t_sim_sample,hcat(z_sample...)[1,:],mark="",
	style="color=black, line width=2pt",legendentry="pos.")
psθ_sim_dpo = PGF.Plots.Linear(t_sim_sample,hcat(z_sample...)[2,:],mark="",
	style="color=black, line width=2pt, densely dashed",legendentry="ang.")

a_dpo = Axis([psx_dpo;psθ_dpo;psx_sim_dpo;
    psθ_sim_dpo],
    xmin=0., ymin=-1, xmax=sum(H_nom_sample), ymax=3.5,
    axisEqualImage=false,
    hideAxis=false,
	ylabel="state",
	xlabel="time",
	legendStyle="{at={(0.01,0.99)},anchor=north west}",
	)

# Save to tikz format
dir = joinpath(@__DIR__,"results")
PGF.save(joinpath(dir,"cartpole_sim_dpo.tikz"), a_dpo, include_preamble=false)

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
function bob(stimulus)

end


text:
using LinearAlgebra
@testset "error analyis" begin
    # Eval runerr
    # Same as Example 17 in manuscript
    graph = Compgraph(Float64);
    add_lincomb!(graph, :y, 1, :I, 1, :x)
    add_lincomb!(graph, :z, 1, :I, 1/2, :x)
    add_mult!(graph, :y2, :y, :y)
    add_mult!(graph, :z2, :z, :z)
    add_lincomb!(graph, :out, 1, :y2, -1, :z2)
    add_output!(graph, :out)
    x = 2. .^ -27
    running_error_bound = eval_runerr(graph, x, input = :x)
    @test running_error_bound < 3e-7;

    x = [1., 2.] .^ -27
    running_error_bound = eval_runerr(graph, x, input = :x)
    @test all(running_error_bound .< [2e-7, 3e-7])

    # compute theta
    (graph,_)=graph_ps([1.0,1.0, 1 / 2, 1 /6])
    q=compute_bnd_rel_bwd_err(:exp,graph,numterms=10,numdigits=30)
    @test q(0.1) < 1e-04

end


text:
## Various utility functions ##
export get_log_norm_weights, get_norm_weights
export effective_sample_size, get_ess
export mean, var

using Gen: effective_sample_size
using Statistics

@inline function update_refs!(state::ParticleFilterState)
    # Swap references
    tmp = state.traces
    state.traces = state.new_traces
    state.new_traces = tmp
end

@inline function update_refs!(state::ParticleFilterSubState)
    state.traces[:] = state.new_traces
end

lognorm(v::AbstractVector) = v .- logsumexp(v)

"""
    get_log_norm_weights(state::ParticleFilterState)

Return the vector of normalized log weights for the current state,
one for each particle.
"""
get_log_norm_weights(state::ParticleFilterView) = lognorm(state.log_weights)

"""
    get_norm_weights(state::ParticleFilterState)

Return the vector of normalized weights for the current state,
one for each particle.
"""
get_norm_weights(state::ParticleFilterView) = exp.(get_log_norm_weights(state))

"""
    effective_sample_size(state::ParticleFilterState)

Computes the effective sample size of the particles in the filter.
"""
Gen.effective_sample_size(state::ParticleFilterView) =
    Gen.effective_sample_size(get_log_norm_weights(state))

"""
    get_ess(state::ParticleFilterState)

Alias for `effective_sample_size`(@ref). Computes the effective sample size.
"""
get_ess(state::ParticleFilterView) = Gen.effective_sample_size(state)

"""
    mean(state::ParticleFilterState[, addr])

Returns the weighted empirical mean for a particular trace address `addr`.
If `addr` is not provided, returns the empirical mean of the return value.
"""
Statistics.mean(state::ParticleFilterView, addr) =
    sum(get_norm_weights(state) .* getindex.(state.traces, addr))

Statistics.mean(state::ParticleFilterView) =
    sum(get_norm_weights(state) .* get_retval.(state.traces))

"""
    var(state::ParticleFilterState[, addr])

Returns the empirical variance for a particular trace address `addr`.
If `addr` is not provided, returns the empirical variance of the return value.
"""
Statistics.var(state::ParticleFilterView, addr) =
    sum(get_norm_weights(state) .*
        (getindex.(state.traces, addr) .- mean(state, addr)).^2)

Statistics.var(state::ParticleFilterView) =
    sum(get_norm_weights(state) .*
        (get_retval.(state.traces) .- mean(state)).^2)


text:
using Test
using JSON
using MianalyzerBackend
const TimeSeriesController = MianalyzerBackend.TimeSeriesController

include("../TestSupport.jl")

@testset "TimeSeriesController" begin
  @testset ".call" begin
    request = Dict([
                    ("records", [[1.0, 100.0], [2.0, 200.0], [3.0, 300.0], [4.0, 200.0], [5.0, 100.0]]),
                    ("time_category", "number"),
                    ("grouping_type", "min")
                   ])

    actual = JSON.parse(IOBuffer(TimeSeriesController.call(request)))

    @test actual["correlogram"] == [[0.0, 1.0], [1.0, 0.05714285714285714], [2.0, -0.6714285714285714], [3.0, -0.11428571428571428], [4.0, 0.22857142857142856]]
    @test actual["fft"] == [["1/10", 1.0], ["1/4", 0.2908926654166549], ["1/2", 0.042440667916678346]]
    @test actual["hist"] == [[100.0, 2.0], [150.0, 0.0], [200.0, 2.0], [250.0, 0.0], [300.0, 1.0]]
  end
end


text:
"""
Likelihood-based calculations for ordinary least squares for the same right hand side (covariates) repeatedly.
"""
module BatchOLS

using ArgCheck

"""
    RHS(X, [Xfact])

Wrapper for covariates of an ordinary least squares (OLS) regression that precalculates the necessary factorizations. The latter can be supplied, but is not checked for consistency.
"""
struct RHS{TX, TXfact}
    "Regressors (covariates)."
    X::TX
    "Factorization of X that allows solving the least squares problem."
    Xfact::TXfact
end

function RHS(X::AbstractMatrix)
    N, K = size(X)
    @argcheck N ≥ K "Need more rows than columns."
    RHS(X, qrfact(X, Val{true}))
end

Base.size(rhs::RHS, dim...) = size(rhs.X, dim...)

Base.:*(rhs::RHS, β) = rhs.X * β

"""
    β, v = ML_βv(y, rhs)

Maximum likelihood estimate of

``y = X⋅β + ϵ``

where ``ϵ ∼ Normal(0, v)``, IID, and `v` is the variance. Use `RHS` to pre-calculate `rhs`.
"""
function ML_βv(y, rhs::RHS)
    β = rhs.Xfact \ y
    e = y-rhs.X*β
    β, mean(abs2, e)
end

"""
    loglikelihood(y, rhs, β, v)

Log likelihood of observations under given parameters `β`, `v` (not necessarily ML estimates). See `ML_βv`.
"""
loglikelihood(y, rhs::RHS, β, v) =
    -(size(rhs, 1)/2*(log(2*π)+log(v)) + sum(abs2, y-rhs.X*β)/(2*v))

end # module


text:
using Distributions, ForwardDiff, Plots, LinearAlgebra
using LogDensityProblems, Parameters, PositiveFactorizations

include(string(pwd(),"/hmc/rkhs.jl"))

include(string(pwd(),"/src/simulation.jl"))
include(string(pwd(),"/gali_bayesian.jl"))
include(string(pwd(),"/mcmc/priors_v2.jl"))
include(string(pwd(),"/mcmc/foos.jl"))

yy,shocks = simulate_dsge(GAMMA_0,GAMMA_1,PSI,PI,500)

prob = DSGE_Model(yy[:,2],1/3)

prob((bet = 0.99,epsilon = 6,theta=2/3,sig=1,s2=1,phi=1,phi_pi=1.5,phi_y=0.5/4,rho_v=0.5))

t = problem_transform(prob)

P = TransformedLogDensity(t,prob)

true_vals = TransformVariables.inverse(t,(bet = 0.99,epsilon = 6,theta=2/3,sig=1,s2=1,phi=1,phi_pi=1.5,phi_y=0.5/4,rho_v=0.5))

LogDensityProblems.logdensity(P,true_vals)

unif = Uniform(0,1)

M = I(9)
M_inv = inv(M)

replics = 3000
warm_up_size = 2000

moment_dist = MvNormal(M)

par = zeros(replics,9)

par[1,:] = [rand(prior_bet);rand(prior_epsilon);rand(prior_theta);rand(prior_sig);rand(prior_phi);rand(prior_phi_pi);rand(prior_phi_y);rand(prior_rho_v);rand(prior_s2)]

N = 300

# dados = zeros(N,10)
#
# for i in 1:N
#     xx = [rand(prior_bet);rand(prior_epsilon);rand(prior_theta);rand(prior_sig);rand(prior_phi);rand(prior_phi_pi);rand(prior_phi_y);rand(prior_rho_v);rand(prior_s2)]
#     dados[i,:] = [xx;LogDensityProblems.logdensity(P,xx)]
# end

par_warmup = zeros(warm_up_size,9)

par_warmup[1,:] = par[1,:]

coef_escala = 0.2#2.4/sqrt(9)

j = 2
rejec = 0

while j <= warm_up_size
    kernel_velho = MvNormal(par_warmup[j-1,:],coef_escala*M)
    novo_par = rand(kernel_velho)
    #pars_aceitos[j,2:10] = novo_par
    kernel_novo = MvNormal(novo_par,coef_escala*M)

    teste = LogDensityProblems.logdensity(P,novo_par)
    if isnan(teste)
        #pars_aceitos[j,2:10] = pars_aceitos[j-1,2:10]
        #global rejec += 1
        #global j += 1
        continue
    else
        num = teste + logpdf(kernel_novo,par_warmup[j-1,:])
        #println(string(teste," ",logpdf(kernel_novo, pars_aceitos[j-1,2:10])))
        dem = LogDensityProblems.logdensity(P,par_warmup[j-1,:]) + logpdf(kernel_velho,novo_par)
        alpha = exp(min(0,num - dem))
        p = rand(unif)
        if alpha > p
            par_warmup[j,:] = novo_par
        else
            par_warmup[j,:] = par_warmup[j-1,:]
            global rejec += 1
        end
        acc = 1 - rejec/j
        if j % 50 == 0
            println("Iteração ", j, " taxa de aceitação ", acc)
            sleep(0.4)
            cov_aux = cov(par_warmup[1:j,:])
            cov_aux = cholesky(Positive,cov_aux)
            global M = cov_aux.L*cov_aux.L'
            par_warmup[j,:] = rand(MvNormal(M))
        end
        global j += 1
    end
end

M_aux = cov(par_warmup[1000:2000,:])
M_aux = inv(M_aux)

M = diagm(diag(M_aux))*1e-9
M_inv = inv(M)

sigg = I(9)

# wp_idx = sample(1:size(dados,1),N)
#
# ll_val = [LogDensityProblems.logdensity(P,dados[wp_idx[j],1:9]) for j in 1:N]
#
# rkhs_approx = rkhs_gaussian(dados[wp_idx,1:9],ll_val,sigg,0)

wp_idx = sample(1000:size(par_warmup,1),N)

ll_val = [LogDensityProblems.logdensity(P,par_warmup[wp_idx[j],1:9]) for j in 1:N]

rkhs_approx = rkhs_gaussian(par_warmup[wp_idx,1:9],ll_val,sigg,0)

f_hat(x) = rkhs_approx(x)

dif_aprox(x) = ForwardDiff.gradient(f_hat,x)

moment_dist = MvNormal(M)

#moments = zeros(replics,5000,9)

i = 2
accep = 0

while i <= replics
    moment = rand(moment_dist)
    old_moment = moment
    par_aux = par[i-1,:]

    L = 1e4#ceil(rand(Uniform(1e4,6e4)))
    leapfrog_step = 1/L#rand(Uniform(1e-9,6e-9))

    for j in 1:L
        new_moment = moment + 1/2*leapfrog_step*ForwardDiff.gradient(pp->f_hat(pp),par_aux)
         par_aux = par_aux .+ leapfrog_step*M_inv*new_moment
         moment = new_moment + 1/2*leapfrog_step*ForwardDiff.gradient(pp->f_hat(pp),par_aux)
         #moments[i,Int64(j),:] = moment
    end

    #println(par_aux)
    num = LogDensityProblems.logdensity(P, par_aux) + logpdf(moment_dist,moment)#+ logpdf(kernel_novo,par[i-1])

    dem = LogDensityProblems.logdensity(P, par[i-1,:]) + logpdf(moment_dist,old_moment)#+ logpdf(kernel,novo_par)

    alpha = exp(min(0,num - dem))

    p = rand(unif)
    if alpha < p
        par[i,:] = par[i-1,:]
    else
        par[i,:] = par_aux
        global accep += 1
    end
    println("Step ",i, " Acceptance ", accep/i)
    # if i % 300 == 0
    #     idxs = sample(1:i,N)
    #     novos_dados = par[idxs,:]
    #     novas_ll = [LogDensityProblems.logdensity(P,novos_dados[j,:]) for j in 1:N]
    #     rkhs_approx = rkhs_gaussian(novos_dados,novas_ll,sig,0)
    #     global f_hat(x) = rkhs_approx(x)
    # end
    global i += 1
end

pars_convertidos = zeros(replics,9)

for i in 1:replics
    @unpack bet,epsilon,theta,sig,s2,phi,phi_pi,phi_y,rho_v = TransformVariables.transform(t,par[i,:])
    pars_convertidos[i,:] = [bet,epsilon,theta,sig,s2,phi,phi_pi,phi_y,rho_v]
end

using StatsPlots

StatsPlots.density(pars_convertidos[1000:3000,1])
StatsPlots.density(pars_convertidos[1000:5000,2])
StatsPlots.density(pars_convertidos[1000:5000,3])
StatsPlots.density(pars_convertidos[1000:5000,4])
StatsPlots.density(pars_convertidos[1000:5000,5])

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
# ---
# title: 1614. Maximum Nesting Depth of the Parentheses
# id: problem1614
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/maximum-nesting-depth-of-the-parentheses/description/>
# hidden: true
# ---
#
# A string is a **valid parentheses string** (denoted **VPS** ) if it meets one
# of the following:
#
#   * It is an empty string `""`, or a single character not equal to `"("` or `")"`,
#   * It can be written as `AB` (`A` concatenated with `B`), where `A` and `B` are **VPS** 's, or
#   * It can be written as `(A)`, where `A` is a **VPS**.
#
# We can similarly define the **nesting depth** `depth(S)` of any VPS `S` as
# follows:
#
#   * `depth("") = 0`
#   * `depth(C) = 0`, where `C` is a string with a single character not equal to `"("` or `")"`.
#   * `depth(A + B) = max(depth(A), depth(B))`, where `A` and `B` are **VPS** 's.
#   * `depth("(" + A + ")") = 1 + depth(A)`, where `A` is a **VPS**.
#
# For example, `""`, `"()()"`, and `"()(()())"` are **VPS** 's (with nesting
# depths 0, 1, and 2), and `")("` and `"(()"` are not **VPS** 's.
#
# Given a **VPS** represented as string `s`, return _the **nesting depth** of
# _`s`.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "(1+(2*3)+(( _8_ )/4))+1"
#     Output: 3
#     Explanation: Digit 8 is inside of 3 nested parentheses in the string.
#
#
# **Example 2:**
#
#
#
#     Input: s = "(1)+((2))+((( _3_ )))"
#     Output: 3
#
#
# **Example 3:**
#
#
#
#     Input: s = "1+( _2_ *3)/(2-1)"
#     Output: 1
#
#
# **Example 4:**
#
#
#
#     Input: s = " _1_ "
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 100`
#   * `s` consists of digits `0-9` and characters `'+'`, `'-'`, `'*'`, `'/'`, `'('`, and `')'`.
#   * It is guaranteed that parentheses expression `s` is a **VPS**.
#
#
## @lc code=start
using LeetCode

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


text:
abstract type AbstractAttitude end

@with_kw struct RandomWalkAttitude <: AbstractAttitude
    yaw::SVector{2, Float64} = SVector(0.0,0.0)
    pitch::SVector{2, Float64} = SVector(0.0,0.0)
    roll::SVector{2, Float64} = SVector(0.0,0.0)
    yaw_vel_std::Float64
    pitch_vel_std::Float64
    roll_vel_std::Float64
    pitch_limit::Float64
    roll_limit::Float64
end

@inline function propagate(attitude::RandomWalkAttitude, Δt, rng)
    P = get_process(Order(2), Δt / s)
    C = cholesky(get_process_covariance(Order(2), Δt / s)).L
    pitch_limit = attitude.pitch_limit
    roll_limit = attitude.roll_limit
    next_yaw = P * attitude.yaw + C * randn(rng, SVector{2,Float64}) * attitude.yaw_vel_std
    next_pitch = soft_bound(P * attitude.pitch, C * randn(rng, SVector{2,Float64}) * attitude.pitch_vel_std, pitch_limit, -pitch_limit)
    next_roll = soft_bound(P * attitude.roll, C * randn(rng, SVector{2,Float64}) * attitude.roll_vel_std, roll_limit, -roll_limit)
    RandomWalkAttitude(next_yaw, next_pitch, next_roll, attitude.yaw_vel_std, attitude.pitch_vel_std, attitude.roll_vel_std, pitch_limit, roll_limit)
end

function get_attitude(attitude::RandomWalkAttitude)
    RotXYZ(attitude.roll[1], attitude.pitch[1], attitude.yaw[1])
end

struct NoisyStaticAttitude <: AbstractAttitude
    attitude::RotXYZ{Float64}
    base_attitude::RotXYZ{Float64}
    roll_std::Float64
    pitch_std::Float64
    yaw_std::Float64
end

function NoisyStaticAttitude(attitude::RotXYZ, roll_std, pitch_std, yaw_std)
    NoisyStaticAttitude(attitude, attitude, roll_std, pitch_std, yaw_std)
end

struct DynamicAttitude <: AbstractAttitude
    attitudes::Vector{RotXYZ{Float64}}
    sample_freq::typeof(1.0Hz)
    time::typeof(1.0s)
end

function DynamicAttitude(attitudes, sample_freq; time = 0.0s)
    DynamicAttitude(attitudes, sample_freq, time)
end

@with_kw struct LinearDynamicAttitude <: AbstractAttitude
    attitude::RotXYZ{Float64}
    Δroll::typeof(1.0rad / 1.0s) = 0.0rad / 1.0s
    Δpitch::typeof(1.0rad / 1.0s) = 0.0rad / 1.0s
    Δyaw::typeof(1.0rad / 1.0s) = 0.0rad / 1.0s
end

struct NoisyDynamicAttitude <: AbstractAttitude
    attitude::RotXYZ{Float64}
    attitudes::Vector{RotXYZ{Float64}}
    sample_freq::typeof(1.0Hz)
    roll_std::Float64
    pitch_std::Float64
    yaw_std::Float64
    time::typeof(1.0s)
end

function NoisyDynamicAttitude(
    attitudes::Vector{<:RotXYZ},
    sample_freq,
    roll_std,
    pitch_std,
    yaw_std;
    time = 0.0s
)
    NoisyDynamicAttitude(
        first(attitudes),
        attitudes,
        sample_freq,
        roll_std,
        pitch_std,
        yaw_std,
        time
    )
end

@with_kw struct NoisyLinearDynamicAttitude <: AbstractAttitude
    attitude::RotXYZ{Float64}
    base_attitude::RotXYZ{Float64}
    Δroll::typeof(1.0rad / 1.0s) = 0.0rad / 1.0s
    Δpitch::typeof(1.0rad / 1.0s) = 0.0rad / 1.0s
    Δyaw::typeof(1.0rad / 1.0s) = 0.0rad / 1.0s
    roll_std::Float64
    pitch_std::Float64
    yaw_std::Float64
end

function NoisyLinearDynamicAttitude(
    attitude::RotXYZ,
    Δroll,
    Δpitch,
    Δyaw,
    roll_std,
    pitch_std,
    yaw_std
)
    NoisyLinearDynamicAttitude(
        attitude,
        attitude,
        Δroll,
        Δpitch,
        Δyaw,
        roll_std,
        pitch_std,
        yaw_std
    )
end

function get_attitude(attitude::A) where A <: AbstractAttitude
    attitude.attitude
end

"""
$(SIGNATURES)

Simulates a static noise-free attitude in radian.
"""
function propagate(attitude::RotXYZ, Δt, rng)
    attitude
end

function get_attitude(attitude::RotXYZ)
    attitude
end

"""
$(SIGNATURES)

Simulates a static noisy attitude in radian.
"""
function propagate(attitude::NoisyStaticAttitude, Δt, rng)
    T = attitude.base_attitude
    next_T = RotXYZ(
        T.theta1 + randn(rng) * attitude.roll_std,
        T.theta2 + randn(rng) * attitude.pitch_std,
        T.theta3 + randn(rng) * attitude.yaw_std
    )
    NoisyStaticAttitude(next_T, T, attitude.roll_std, attitude.pitch_std, attitude.yaw_std)
end

"""
$(SIGNATURES)

Simulates a time-varying noise-free attitude in radian.
If time index exceeds data length, last available value is returned.
"""
function propagate(attitude::DynamicAttitude, Δt, rng)
    next_time = attitude.time + Δt
    DynamicAttitude(attitude.attitudes, attitude.sample_freq, next_time)
end

function get_attitude(attitude::DynamicAttitude)
    get_sampled_value(attitude.time, attitude.sample_freq, attitude.attitudes)
end

"""
$(SIGNATURES)

Simulates a time-varying noisy attitude in radian.
If time index exceeds data length, last available value is returned.
"""
function propagate(attitude::NoisyDynamicAttitude, Δt, rng)
    next_time = attitude.time + Δt
    T = get_sampled_value(attitude.time, attitude.sample_freq, attitude.attitudes)
    next_T = RotXYZ(
        T.theta1 + randn(rng) * attitude.roll_std,
        T.theta2 + randn(rng) * attitude.pitch_std,
        T.theta3 + randn(rng) * attitude.yaw_std
    )
    NoisyDynamicAttitude(
        next_T,
        attitude.attitudes,
        attitude.sample_freq,
        attitude.roll_std,
        attitude.pitch_std,
        attitude.yaw_std,
        next_time
    )
end

"""
$(SIGNATURES)

Simulates a linearly time-varying noise-free attitude in radian.
"""
function propagate(attitude::LinearDynamicAttitude, Δt, rng)
    next_attitude = RotXYZ(
        attitude.attitude.theta1 + attitude.Δroll * Δt,
        attitude.attitude.theta2 + attitude.Δpitch * Δt,
        attitude.attitude.theta3 + attitude.Δyaw * Δt
    )
    LinearDynamicAttitude(next_attitude, attitude.Δroll, attitude.Δpitch, attitude.Δyaw)
end

"""
$(SIGNATURES)

Simulates a linearly time-varying noisy attitude in radian.
"""
function propagate(attitude::NoisyLinearDynamicAttitude, Δt, rng)
    T = attitude.base_attitude
    next_T = RotXYZ(
        T.theta1 + attitude.Δroll * Δt,
        T.theta2 + attitude.Δpitch * Δt,
        T.theta3 + attitude.Δyaw * Δt
    )
    noisy_next_T = RotXYZ(
        next_T.theta1 + randn(rng) * attitude.roll_std,
        next_T.theta2 + randn(rng) * attitude.pitch_std,
        next_T.theta3 + randn(rng) * attitude.yaw_std
    )
    NoisyLinearDynamicAttitude(
        noisy_next_T,
        next_T,
        attitude.Δroll,
        attitude.Δpitch,
        attitude.Δyaw,
        attitude.roll_std,
        attitude.pitch_std,
        attitude.yaw_std
    )
end


text:
@testset "MixSkewT" begin
  Random.seed!(0)
  println("Test MixSkewT")

  y = rand(SkewT(1, .6, 9, -5), 1000)

  model = MixSkewT(y, 5)
  init = MCMC.make_init_state(model)
  spl = make_sampler(model, init=init)

  nburn = 50
  nsamps = 100
  thin = 2

  chain, metrics = mcmc(spl, nsamps, init=init, nburn=nburn, thin=thin)
  @test length(unique(chain)) == nsamps

  # TODO:
  # - [ ] Test that the parameters are correct.

  syms = keys(init)
  for s in syms
    # println(s)
    samps = getindex.(chain, s)
    if s != :nu
      @test length(unique(samps)) > nsamps * 0.8
    else
      @test length(unique(samps)) > 2
    end
  end
end


text:
function sortieren(A,m,n)
    #n steht für 0 oder 1
    # 0 bedeutet aufsteigend und 1 absteigend
    # m steht für die Spalte nach der sortiert werden soll

    C = zeros(size(A,1),size(A,2));

    if n == 0
        vec = A[:,m];
        vec = sort!(vec);
        for i = 1:size(A,1)
            Ort = finde(A[:,m],vec[i]);
            if size(Ort,1) ==1
                x = Int(Ort[1,1]);
                for j = 1:size(A,2)
                    C[i,j] = A[x,j];
                end
            else
                for k = 1:size(Ort,1)
                    x = Int(Ort[k,1]);
                    for j = 1:size(A,2)
                        C[i,j] = A[x,j];
                    end
                    end
            end
        end
    elseif n == 1
        vec = A[:,m];
        vec = sort!(vec,rev = true);

        for i = 1:size(A,1)
            Ort = finde(A[:,m],vec[i]);
            if size(Ort,1) ==1
                x = Int(Ort[1,1]);
                for j = 1:size(A,2)
                    C[i,j] = A[x,j];
                end
            else
                for k = 1:size(Ort,1)
                    x = Int(Ort[k,1]);
                    for j = 1:size(A,2)
                        C[i,j] = A[x,j];
                    end
                end
            end
        end
    end
return C;

end


text:
@testset "util.iterator" begin

@testset "combinations" begin
    @test length(collect(combinations(1:10, 1))) == 10
    @test length(collect(combinations(1:10, 2))) == 45
    @test length(collect(combinations(1:10, 3))) == 120
    @test length(collect(combinations(1:10, 8))) == 45
    @test length(collect(combinations(1:10, 9))) == 10
    @test length(collect(combinations(1:10, 10))) == 1
    @test length(collect(combinations(1:10, 0))) == 1
    @test_throws DomainError combinations(1:10, -1)
    @test_throws ErrorException combinations(1:10, 11)
    @test_throws ErrorException combinations([], 1)
    @test_throws MethodError combinations(1:10, 1.5)
end

@testset "sortstablemax" begin
    @test sortstablemax([1, 4, 2, 7, 5]) == 7
    @test sortstablemax([8:10, 1:3, [], 4:23, 1:9], by=length) == 4:23
    @test sortstablemax([8:10, 1:3, 3:22, 4:23, 1:9], by=length) == 3:22
    @test_throws ErrorException sortstablemax([])
    @test sortstablemax([], init=0) == 0
end

@testset "sortstablemin" begin
    @test sortstablemin([1, 4, 2, 7, 5]) == 1
    @test sortstablemin([8:10, 1:3, [], 4:23, 1:9], by=length) == []
    @test sortstablemin([8:10, 1:3, 3:22, 4:23, 1:9], by=length) == 8:10
    @test_throws ErrorException sortstablemin([])
    @test sortstablemin([], init=0) == 0
end
end # util.iterator


text:
module FESpaceFactoriesTests

using Test

using Gridap.Geometry
using Gridap.Arrays
using Gridap.TensorValues
using Gridap.ReferenceFEs
using Gridap.FESpaces
using Gridap.CellData

domain = (0,1,0,1)
partition = (10,10)
model = CartesianDiscreteModel(domain,partition)
order = 3

# From a single reffe

V = FESpace(model,QUAD4)
@test isa(V,UnconstrainedFESpace)

V = FESpace(model,QUAD4,conformity=:H1)
@test isa(V,UnconstrainedFESpace)
@test num_free_dofs(V) == 121

V = FESpace(model,QUAD4,conformity=:L2)
@test isa(V,UnconstrainedFESpace)
@test num_free_dofs(V) == 400

V = FESpace(model,QUAD4,conformity=:L2,constraint=:zeromean)
@test isa(V,ZeroMeanFESpace)

# From parameter list describing the reffe

reffe = ReferenceFE(lagrangian,Float64,order)

V = FESpace(model,reffe)
@test isa(V,UnconstrainedFESpace)

V = FESpace(model,reffe,conformity=:H1)
@test isa(V,UnconstrainedFESpace)

V = FESpace(model,reffe,conformity=:L2)
@test isa(V,UnconstrainedFESpace)

V = FESpace(model,reffe,conformity=:L2,constraint=:zeromean)
@test isa(V,ZeroMeanFESpace)

# From a cell-wise vector of reffes

cell_reffe = ReferenceFE(model,lagrangian,Float64,order)

V = FESpace(model,cell_reffe)
@test isa(V,UnconstrainedFESpace)

V = FESpace(model,cell_reffe,conformity=:H1)
@test isa(V,UnconstrainedFESpace)

V = FESpace(model,cell_reffe,conformity=:L2)
@test isa(V,UnconstrainedFESpace)

V = FESpace(model,cell_reffe,conformity=:L2,constraint=:zeromean)
@test isa(V,ZeroMeanFESpace)

# From a CellFE

cell_fe = FiniteElements(PhysicalDomain(),model,lagrangian,Float64,order)

V = FESpace(model,cell_fe)
@test isa(V,UnconstrainedFESpace)

cell_fe = FiniteElements(PhysicalDomain(),model,lagrangian,Float64,order,conformity=:L2)
V = FESpace(model,cell_fe)
@test isa(V,UnconstrainedFESpace)

V = FESpace(model,cell_fe,constraint=:zeromean)
@test isa(V,ZeroMeanFESpace)

# Now from a restricted model

const R = 0.7

function is_in(coords)
  n = length(coords)
  x = (1/n)*sum(coords)
  d = x[1]^2 + x[2]^2 - R^2
  d < 0
end

Ω = Triangulation(model)
cell_coords = get_cell_coordinates(Ω)
cell_mask = lazy_map(is_in,cell_coords)
Ω_in = Triangulation(model,cell_mask)

# From a single reffe

V = FESpace(Ω_in,QUAD4)
@test isa(V,UnconstrainedFESpace)

V = FESpace(Ω_in,QUAD4,conformity=:L2)
@test isa(V,UnconstrainedFESpace)

V = FESpace(Ω_in,QUAD4,conformity=:L2,constraint=:zeromean)
@test isa(V,ZeroMeanFESpace)

# From parameter list describing the reffe

reffe = ReferenceFE(lagrangian,Float64,order)

V = FESpace(Ω_in,reffe)
@test isa(V,UnconstrainedFESpace)

V = FESpace(Ω_in,reffe,conformity=:L2)
@test isa(V,UnconstrainedFESpace)

V = FESpace(Ω_in,reffe,conformity=:L2,constraint=:zeromean)
@test isa(V,ZeroMeanFESpace)

# From a cell-wise vector of reffes

model_in = get_active_model(Ω_in)
cell_reffe = ReferenceFE(model_in,lagrangian,Float64,order)

V = FESpace(model_in,cell_reffe,trian=Ω_in)
@test isa(V,UnconstrainedFESpace)
@test Ω_in === get_triangulation(V)

V = FESpace(model_in,cell_reffe,trian=Ω_in,conformity=:L2)
@test isa(V,UnconstrainedFESpace)
@test Ω_in === get_triangulation(V)

V = FESpace(model_in,cell_reffe,trian=Ω_in,conformity=:L2,constraint=:zeromean)
@test isa(V,ZeroMeanFESpace)
@test Ω_in === get_triangulation(V)

# From a CellFE

cell_fe = FiniteElements(PhysicalDomain(),model_in,lagrangian,Float64,order)

V = FESpace(model_in,cell_fe,trian=Ω_in)
@test isa(V,UnconstrainedFESpace)
@test Ω_in === get_triangulation(V)

cell_fe = FiniteElements(PhysicalDomain(),model_in,lagrangian,Float64,order,conformity=:L2)
V = FESpace(model_in,cell_fe,trian=Ω_in)
@test isa(V,UnconstrainedFESpace)
@test Ω_in === get_triangulation(V)

V = FESpace(model_in,cell_fe,trian=Ω_in,constraint=:zeromean)
@test isa(V,ZeroMeanFESpace)
@test Ω_in === get_triangulation(V)

Γ = Boundary(model)
reffe = ReferenceFE(lagrangian,Float64,1)
V = FESpace(Γ,reffe)
@test get_triangulation(V) === Γ
V = FESpace(Γ,reffe,conformity=:H1)
@test get_triangulation(V) === Γ
V = FESpace(Γ,reffe,conformity=:L2)
@test get_triangulation(V) === Γ

end # module


text:
export tikz,
       tikz_output,
       display_tikz_output

tikz_e(out::IO, arg) =
  begin
    print(out, arg)
    println(out, ";")
  end

tikz_draw(out::IO, filled=false) = print(out, filled ? "\\fill " : "\\draw ")

tikz_number(out::IO, x::Real) =
  isinteger(x) ? print(out, x) : (abs(x) < 0.0001 ? print(out, 0) : print(out, round(x*10000.0)/10000.0))

tikz_cm(out::IO, x::Real) = begin
  tikz_number(out, x)
  print(out, "cm")
end

tikz_coord(out::IO, c::Loc) =
  begin
    print(out, "(")
    tikz_number(out, c.x)
    print(out, ",")
    tikz_number(out, c.y)
    if !iszero(c.z)
      print(out, ",")
      tikz_number(out, c.z)
    end
    print(out, ")")
  end

tikz_pgfpoint(out::IO, c::Loc) =
  begin
    if ! iszero(c.z)
      error("Can't handle 3D coords")
    end
    print(out, "\\pgfpoint{")
    tikz_cm(out, c.x)
    print(out, "}{")
    tikz_cm(out, c.y)
    print(out, "}")
  end

tikz_circle(out::IO, c::Loc, r::Real, filled::Bool=false) =
  begin
    tikz_draw(out, filled)
    tikz_coord(out, c)
    print(out, "circle(")
    tikz_cm(out, r)
    tikz_e(out, ")")
  end

tikz_point(out::IO, c::Loc) = tizk_circle(out, c, 0.01, true)

tikz_ellipse(out::IO, c::Loc, r0::Real, r1::Real, fi::Real, filled=false) =
  begin
    tikz_draw(out, filled)
    print(out, "[shift={")
    tikz_coord(out, c)
    print(out, "}]")
    print(out, "[rotate=")
    tikz_number(out, rad2deg(fi))
    print(out, "]")
    print(out, "(0,0)")
    print(out, "ellipse(")
    tikz_cm(out, r0)
    print(out, " and ")
    tikz_cm(out, r1)
    tikz_e(out, ")")
  end

tikz_arc(out::IO, c::Loc, r::Real, ai::Real, af::Real, filled=false) =
  begin
    tikz_draw(out, filled)
    if filled
      tikz_coord(out, c)
      print(out, "--")
    end
    tikz_coord(out, c+vpol(r, ai))
    print(out, "arc(")
    tikz_number(out, rad2deg(ai))
    print(out, ":")
    tikz_number(out, rad2deg(ai > af ? af+2*pi : af))
    print(out, ":")
    tikz_cm(out, r)
    print(out, ")")
    if filled
      tikz_e(out, "--cycle")    end
    println(out, ";")
  end

tikz_maybe_arc(out::IO, c::Loc, r::Real, ai::Real, da::Real, filled=false) =
  if iszero(r)
    tikz_point(out, c)
  elseif iszero(da)
    tikz_point(out, c + vpol(r, ai))
  elseif abs(da) >= 2*pi
    tikz_circle(out, c, r, filled)
  else
    let af = ai + da
      if af > ai
        tikz_arc(out, c, r, ai, af, filled)
      else
        tikz_arc(out, c, r, af, ai, filled)
      end
    end
  end

tikz_line(out::IO, pts::Locs) =
  begin
    tikz_draw(out, false)
    tikz_coord(out, first(pts))
    for pt in Iterators.drop(pts, 1)
      print(out, "--")
      tikz_coord(out, pt)
    end
    println(out, ";")
  end

tikz_closed_line(out::IO, pts::Locs, filled::Bool=false) =
  begin
    tikz_draw(out, filled)
    for pt in pts
      tikz_coord(out, pt)
      print(out, "--")
    end
    tikz_e(out, "cycle")
  end

tikz_spline(out::IO, pts::Locs, filled::Bool=false) =
  begin
    tikz_draw(out, filled)
    print(out, "plot [smooth,tension=1] coordinates {")
    for pt in pts
      tikz_coord(out, pt)
    end
    tikz_e(out, "}")
  end

tikz_closed_spline(out::IO, pts::Locs, filled::Bool=false) =
  begin
    tikz_draw(out, filled)
    print(out, "plot [smooth cycle,tension=1] coordinates {")
    for pt in pts
      tikz_coord(out, pt)
    end
    tikz_e(out, "}")
  end

# HACK we need to handle the starting and ending vectors
tikz_hobby_spline(out::IO, pts::Locs, filled::Bool=false) =
  begin
    tikz_draw(out, filled)
    print(out, "[hobby]")
    print(out, "plot coordinates {")
    for pt in pts
      tikz_coord(out, pt)
    end
    tikz_e(out, "}")
  end

# HACK we need to handle the starting and ending vectors
tikz_hobby_closed_spline(out::IO, pts::Locs, filled::Bool=false) =
  begin
    tikz_draw(out, filled)
    print(out, "[closed hobby]")
    print(out, "plot coordinates {")
    for pt in pts
      tikz_coord(out, pt)
    end
    tikz_e(out, "}")
  end

tikz_rectangle(out::IO, p::Loc, w::Real, h::Real, filled::Bool=false) =
  begin
    tikz_draw(out, filled)
    tikz_coord(out, p)
    print(out, "rectangle")
    tikz_coord(out, p+vxy(w, h))
    println(out, ";")
  end

# Assuming default Arial font for AutoCAD
tikz_text(out::IO, txt, p::Loc, h::Real) =
  let (scale_x, scale_y) = (3.7*h, 3.7*h)
    tikz_draw(out)
    print(out, "[anchor=base west]")
    tikz_coord(out, p)
    print(out, "node[font=\\fontfamily{phv}\\selectfont,outer sep=0pt,inner sep=0pt")
    print(out, ",xscale=")
    tikz_number(out, scale_x)
    print(out, ",yscale=")
    tikz_number(out, scale_y)
    print(out, "]{")
    print(out, txt)
    tikz_e(out, "}")
  end

# FINISH THIS
tikz_transform(out::IO, f::Function, c::Loc) =
  begin
    print(out, "\\begin{scope}")
    print(out, "[shift={")
    tikz_coord(out, c)
    print(out, "}]")
    print(out, "]")
    f(out)
    tikz_e(out, "\\end{scope}")
  end

tikz_set_view(out::IO, camera::Loc, target::Loc, lens::Real) =
  let (v, contents) = (camera-target, get_accumulated_tikz(out))
    print("\\tdplotsetmaincoords{")
    tikz_number(radians_>degrees(sph_psi(v)))
    print("}{")
    tikz_number(radians_>degrees(sph_phi(v))+90)
    println("}")
    println("\\begin{tikzpicture}[tdplot_main_coords]")
    write_bytes(contents, tikz_port)
    println("\\end{tikzpicture}")
  end

#

abstract type TikZKey end
const TikZId = Nothing
const TikZIds = Vector{TikZId}
const TikZRef = GenericRef{TikZKey, TikZId}
const TikZRefs = Vector{TikZRef}
const TikZEmptyRef = EmptyRef{TikZKey, TikZId}
const TikZUniversalRef = UniversalRef{TikZKey, TikZId}
const TikZNativeRef = NativeRef{TikZKey, TikZId}
const TikZUnionRef = UnionRef{TikZKey, TikZId}
const TikZSubtractionRef = SubtractionRef{TikZKey, TikZId}
const TikZ = IOBufferBackend{TikZKey, TikZId}

void_ref(b::TikZ) = TikZNativeRef(nothing)

create_TikZ_connection() = IOBuffer()

const tikz = TikZ(create_TikZ_connection())

get_accumulated_tikz(out::IO) = String(take!(out))

tikz_output(b::TikZ=current_backend()) = get_accumulated_tikz(connection(b))

display_tikz_output(b::TikZ=current_backend()) = print(tikz_output(b))

withTikZXForm(f, out, c) =
  if is_world_cs(c.cs)
    f(out, c)
  else
    tikz_transform(out,
      out -> f(out, u0(world_cs)),
      c)
   end

delete_all_shapes(b::TikZ) =
  truncate(connection(b), 0)

realize(b::TikZ, s::Point) =
  tikz_pgfpoint(connection(b), in_world(s.position))

realize(b::TikZ, s::Line) =
  tikz_line(connection(b), map(in_world, s.vertices))

realize(b::TikZ, s::Spline) =
  if (s.v0 == false) && (s.v1 == false)
    #TikZSpline(connection(b), s.points)
    tikz_hobby_spline(connection(b), map(in_world, s.points), false)
  elseif (s.v0 != false) && (s.v1 != false)
    TikZInterpSpline(connection(b), s.points, s.v0, s.v1)
  else
    TikZInterpSpline(connection(b),
                     s.points,
                     s.v0 == false ? s.points[2]-s.points[1] : s.v0,
                     s.v1 == false ? s.points[end-1]-s.points[end] : s.v1)
  end

realize(b::TikZ, s::ClosedSpline) =
  tikz_hobby_closed_spline(connection(b), map(in_world, s.points))

realize(b::TikZ, s::Circle) =
  withTikZXForm(connection(b), s.center) do out, c
    tikz_circle(out, c, s.radius)
  end
realize(b::TikZ, s::SurfaceCircle) =
  withTikZXForm(connection(b), s.center) do out, c
    tikz_circle(out, c, s.radius, true)
  end

realize(b::TikZ, s::Arc) =
  withTikZXForm(connection(b), s.center) do out, c
    tikz_maybe_arc(out, c, s.radius, s.start_angle, s.amplitude, false)
  end

realize(b::TikZ, s::SurfaceArc) =
  withTikZXForm(connection(b), s.center) do out, c
    tikz_maybe_arc(out, c, s.radius, s.start_angle, s.amplitude, true)
  end

realize(b::TikZ, s::Ellipse) =
  withTikZXForm(connection(b), s.center) do out, c
    tikz_ellipse(out, c, s.radius_x, s.radius_y, 0, false)
  end

realize(b::TikZ, s::SurfaceEllipse) =
  withTikZXForm(connection(b), s.center) do out, c
    tikz_ellipse(out, c, s.radius_x, s.radius_y, 0, true)
  end

realize(b::TikZ, s::EllipticArc) =
  error("Finish this")

#realize(b::TikZ, s::SurfaceElliptic_Arc) = TikZCircle(connection(b),

realize(b::TikZ, s::Polygon) =
  tikz_closed_line(connection(b), map(in_world, s.vertices))

realize(b::TikZ, s::SurfacePolygon) =
  tikz_closed_line(connection(b), map(in_world, s.vertices), true)

realize(b::TikZ, s::RegularPolygon) =
  tikz_closed_line(
    connection(b),
    map(in_world,
        regular_polygon_vertices(s.edges, s.center, s.radius, s.angle, s.inscribed)))

realize(b::TikZ, s::SurfaceRegularPolygon) =
  tikz_closed_line(
    connection(b),
    map(in_world,
        regular_polygon_vertices(s.edges, s.center, s.radius, s.angle, s.inscribed)),
    false)

realize(b::TikZ, s::Rectangle) =
  withTikZXForm(connection(b), s.corner) do out, c
    tikz_rectangle(out, c, s.dx, s.dy)
  end

realize(b::TikZ, s::SurfaceRectangle) =
  withTikZXForm(connection(b), s.corner) do out, c
    tikz_rectangle(out, c, s.dx, s.dy, true)
  end

realize(b::TikZ, s::Text) =
  withTikZXForm(connection(b), s.corner) do out, c
    tikz_text(out, s.str, c, s.height)
  end


text:
module FoldsCUDATests

using GPUArrays
using Test

include("utils.jl")

function include_tests(m = @__MODULE__, dir = @__DIR__)
    for file in readdir(dir)
        if match(r"^test_.*\.jl$", file) !== nothing
            Base.include(m, joinpath(dir, file))
        end
    end
end

include_tests()

module Generic
using ..FoldsCUDATests: include_tests
include_tests(@__MODULE__, joinpath(@__DIR__, "generic"))
end

should_test_gpu() =
    lowercase(get(ENV, "JULIA_PKGEVAL", "false")) != "true" &&
    lowercase(get(ENV, "CUDAFOLDS_JL_TEST_GPU", "true")) == "true"

requires_gpu(m::Module) = parentmodule(m) === @__MODULE__

function collect_modules(root::Module)
    modules = Module[]
    for n in names(root, all = true)
        m = getproperty(root, n)
        m isa Module || continue
        m === root && continue
        startswith(string(nameof(m)), "Test") || continue
        push!(modules, m)
    end
    return modules
end

function collect_modules()
    modules = collect_modules(@__MODULE__)
    append!(modules, collect_modules(Generic))
    return modules
end

this_project() = joinpath(dirname(@__DIR__), "Project.toml")

function is_in_path()
    project = this_project()
    paths = Base.load_path()
    project in paths && return true
    realproject = realpath(project)
    realproject in paths && return true
    matches(path) = path == project || path == realproject
    return any(paths) do path
        matches(path) || matches(realpath(path))
    end
end

function with_project(f)
    is_in_path() && return f()
    load_path = copy(LOAD_PATH)
    push!(LOAD_PATH, this_project())
    try
        f()
    finally
        append!(empty!(LOAD_PATH), load_path)
    end
end

function runtests(modules = collect_modules())
    with_project() do
        runtests_impl(modules)
    end
end

function runtests_impl(modules)
    GPUArrays.allowscalar(false)
    test_gpu = should_test_gpu()
    @testset "$(nameof(m))" for m in modules
        if !test_gpu && requires_gpu(m)
            continue  # branch inside `for` loop for printing skipped tests
        end
        if m in (
            Generic.TestReducePartitionBy,
            TestCoalescedCUDAEx,
            TestTypeChangingAccumulators,
        )
            VERSION < v"1.6-" && continue
        end
        tests = map(names(m, all = true)) do n
            n == :test || startswith(string(n), "test_") || return nothing
            f = getproperty(m, n)
            f !== m || return nothing
            parentmodule(f) === m || return nothing
            applicable(f) || return nothing  # removed by Revise?
            return f
        end
        filter!(!isnothing, tests)
        @testset "$f" for f in tests
            f()
        end
    end
end

""" A list of tests to be run from UnionArrays.jl """
const UNIONARRAYS_TESTS = [TestTypeChangingAccumulators]

function runtests_unionarrays()
    with_project() do
        runtests_impl(UNIONARRAYS_TESTS)
    end
end

end # module