text:
import KernelDensityEstimatePlotting: getColorsByLength

# """
#     $(SIGNATURES)
#
# Standardize the length colors used by RoMEPlotting.
#
# Notes
# - Duplicated in KernelDensityEstimatePlotting
# """
# function getColorsByLength(len::Int=10)::Vector{String}
#   COLORS = String["red";"green";"blue";"black";"deepskyblue";"yellow";"magenta"]
#   morecyan = ["cyan" for i in (length(COLORS)+1):len]
#   retc = [COLORS; morecyan]
#   return retc[1:len]
# end

"""
    $(SIGNATURES)

A peneric KDE plotting function that allows marginals of higher dimensional beliefs and various keyword options.

Example:
```julia
p = kde!(randn(3,100))

plotKDE(p)
plotKDE(p, dims=[1;2], levels=3)
plotKDE(p, dims=[1])

q = kde!(5*randn(3,100))
plotKDE([p;q])
plotKDE([p;q], dims=[1;2], levels=3)
plotKDE([p;q], dims=[1])
```
"""
function plotKDE(fgl::G,
                 sym::Symbol;
                 dims=nothing,
                 title="",
                 levels::Int=5,
                 fill::Bool=false,
                 layers::Bool=false,
                 c=nothing  ) where G <: AbstractDFG
  #
  p = getKDE(getVariable(fgl,sym))
  # mmarg = length(marg) > 0 ? marg : collect(1:Ndim(p))
  # mp = marginal(p,mmarg)
  plotKDE(p, levels=levels, dims=dims, title=string(sym, "  ", title), fill=fill, layers=layers, c=c )
end
function plotKDE(fgl::G,
                 syms::Vector{Symbol};
                 addt::Vector{BallTreeDensity}=BallTreeDensity[],
                 dims=nothing,
                 title=nothing,
                 levels=3,
                 layers::Bool=false,
                 c=getColorsByLength(length(addt))  ) where G <: AbstractDFG
  #
  # TODO -- consider automated rotisary of color
  # colors = ["black";"red";"green";"blue";"cyan";"deepskyblue"; "yellow"]
  # COLORS = repmat(colors, 10)
  # COLORS = getColorsByLength(length(syms))
  MP = BallTreeDensity[]
  LEG = String[]
  # mmarg = Int[]
  for sym in syms
    p = getKDE(getVariable(fgl,sym))
    # mmarg = length(marg) > 0 ? marg : collect(1:Ndim(p))
    # mp = marginal(p,mmarg)
    push!(MP, p)
    push!(LEG, string(sym))
  end
  for p in addt
    # mp = marginal(p,mmarg)
    push!(MP, p)
    push!(LEG, "add")
  end
  plotKDE(MP, c=c, levels=levels, dims=dims, legend=LEG, title=title, layers=layers)
end

"""
    $(SIGNATURES)

Plot absolute values of variables and measurement surrounding fsym factor.
"""
function plotKDEofnc(fgl::G,
                     fsym::Symbol;
                     marg=nothing,
                     N::Int=100  ) where G <: AbstractDFG
  #
  fnc = nothing
  if haskey(fgl.fIDs, fsym)
    fnc = getfnctype( fgl, fgl.fIDs[fsym] )
  else
    error("fIDs doesn't have $(fsym)")
  end
  p = kde!( getSample(fnc, N)[1]  )
  # mmarg = length(marg) > 0 ? marg : collect(1:Ndim(p))
  plotKDE(fgl, lsf(fgl, fsym), addt=[p], marg=marg)
end

"""
    $(SIGNATURES)

Try plot relative values of variables and measurement surrounding fsym factor.
"""
function plotKDEresiduals(fgl::G,
                          fsym::Symbol;
                          N::Int=100,
                          levels::Int=3,
                          dims=nothing,
                          fill=false  ) where G <: AbstractDFG
  #
  # COLORS = ["black";"red";"green";"blue";"cyan";"deepskyblue"]
  fnc = getfnctype( fgl, fgl.fIDs[fsym] )
  # @show sxi = lsf(fgl, fsym)[1]
  # @show sxj = lsf(fgl, fsym)[2]
  fct = getFactor(fgl, fsym)
  @show sxi = getData(fct).fncargvID[1]
  @show sxj = getData(fct).fncargvID[2]
  xi = getVal(fgl, sxi)
  xj = getVal(fgl, sxj)
  measM = getSample(fnc, N)
  meas = length(measM) == 1 ? (0*measM[1], ) : (0*measM[1], measM[2])
  @show size(meas[1])
  d = size(measM[1],1)
  RES = zeros(d,N)
  for i in 1:N
    res = zeros(d)
    fnc(res, i, meas, xi, xj)
    RES[:,i] = res
    if length(measM) > 1
      if measM[2][i] == 0
        RES[:,i] = 0.5*randn(d)
      end
    end
  end
  pR = kde!(RES)
  pM = kde!(measM[1])
  fncvar = getfnctype(fct)

  COLORS = getColorsByLength()
  plotKDE([pR;pM], c=COLORS[1:2], dims=dims,levels=3, legend=["residual";"model"], fill=fill, title=string(fsym, ", ", fncvar))
end

"""
    $(SIGNATURES)

Plot the product of priors and incoming upward messages for variable in clique.

plotPriorsAtCliq(tree, :x2, :x1[, marg=[1;2]])
"""
function plotPriorsAtCliq(tree::BayesTree,
                          lb::Symbol,
                          cllb::Symbol;
                          dims::Vector{Int}=Int[],
                          levels::Int=1,
                          fill::Bool=false  )
  #
  # COLORS = ["black";"red";"green";"blue";"cyan";"deepskyblue"]

  cliq = whichCliq(tree, cllb)
  cliqprs = cliq.attributes["debug"].priorprods[1]

  vidx = 1
  for lbel in cliqprs.lbls
    if lbel == lb
      break;
    else
      vidx += 1
    end
  end
  dims = length(dims)>0 ? dims : collect(1:size(cliqprs.prods[vidx].prev,1))

  arr = BallTreeDensity[]
  push!(arr, marginal(kde!(cliqprs.prods[vidx].prev), dims)  )
  push!(arr, marginal(kde!(cliqprs.prods[vidx].product), dims)  )
  len = length(cliqprs.prods[vidx].potentials)
  lg = String["p";"n"]
  i=0
  for pot in cliqprs.prods[vidx].potentials
    push!(arr, marginal(pot, dims)  )
    i+=1
    push!(lg, cliqprs.prods[vidx].potentialfac[i])
  end

  COLORS = getColorsByLength(len+2)
  # cc = COLORS[1:(len+2)]
  plotKDE(arr, c=COLORS, legend=lg, levels=levels, fill=fill );
end

"""
    $(SIGNATURES)

Draw the up pass belief of lb from clique cllb.

plotUpMsgsAtCliq(tree, :x2, :x1)
"""
function plotUpMsgsAtCliq(treel::BayesTree, lb::Symbol, cllb::Symbol;
      show::Bool=true,
      w=20cm, h=15cm,
      levels::Int=1,
      marg::Vector{Int}=Int[] )
  #
  cliq = whichCliq(treel, string(cllb))
  cliqoutmsg = cliq.attributes["debug"].outmsg
  lbls = cliq.attributes["debug"].outmsglbls

  bel = kde!(cliqoutmsg.p[lbls[lb]])
  bel = length(marg)==0 ? bel : marginal(bel, marg)

  plotKDE(bel)
end

function plotMCMC(treel::BayesTree,
                  lbll::Symbol;
                  delay::Int=200,
                  show::Bool=true,
                  w=20cm, h=15cm,
                  levels::Int=1,
                  dims=nothing  )
  #
  cliq = whichCliq(treel, string(lbll))
  cliqdbg = cliq.attributes["debug"]

  vidx = 1
  for lb in cliqdbg.mcmc[1].lbls
    if lb == lbll
      break;
    else
      vidx += 1
    end
  end

  tmpfilepath = joinpath(dirname(@__FILE__),"tmpimgs")
  ARR = BallTreeDensity[]
  COLORS = getColorsByLength()
  # COLORS = ["black";"red";"green";"blue";"cyan";"deepskyblue";"magenta"]
  for i in 1:length(cliqdbg.mcmc)
    ppr = kde!(cliqdbg.mcmc[i].prods[vidx].prev)
    ppp = kde!(cliqdbg.mcmc[i].prods[vidx].product)
    ARR = [ARR;ppr;ppr;cliqdbg.mcmc[i].prods[vidx].potentials]
  end
  rangeV = getKDERange(ARR)
  ppp = nothing
  for i in 1:length(cliqdbg.mcmc)
    ppr = kde!(cliqdbg.mcmc[i].prods[vidx].prev)
    ppp = kde!(cliqdbg.mcmc[i].prods[vidx].product)
    arr = [ppr;ppp;cliqdbg.mcmc[i].prods[vidx].potentials]
    len = length(cliqdbg.mcmc[i].prods[vidx].potentials)
    lg = String["p";"n";cliqdbg.mcmc[i].prods[vidx].potentialfac] #map(string, 1:len)]
    COLORS_l = getColorsByLength(len+2)
    cc = plotKDE(arr, c=COLORS_l, legend=lg, levels=levels, fill=true, axis=rangeV, dims=dims );
    Gadfly.draw(PNG(joinpath(tmpfilepath,"$(string(lbll))mcmc$(i).png"),w,h),cc)
  end
  # draw initial and final result
  pp0 = kde!(cliqdbg.mcmc[1].prods[vidx].prev)
  i = 0
  cc = plotKDE([pp0], c=[COLORS[1]], legend=["0"], levels=levels, fill=true, axis=rangeV, dims=dims );
  Gadfly.draw(PNG(joinpath(tmpfilepath,"$(string(lbll))mcmc$(i).png"),w,h),cc)

  i = length(cliqdbg.mcmc)+1
  cc = plotKDE([pp0;ppp], c=COLORS[1:2], legend=["0";"n"], levels=levels, fill=true, axis=rangeV, dims=dims );
  Gadfly.draw(PNG(joinpath(tmpfilepath,"$(string(lbll))mcmc$(i).png"),w,h),cc)
  # generate output
  run(`convert -delay $(delay) $(tmpfilepath)/$(string(lbll))mcmc*.png $(tmpfilepath)/$(string(lbll))mcmc.gif`)
  !show ? nothing : (@async run(`eog $(tmpfilepath)/$(string(lbll))mcmc.gif`) )
  return "$(tmpfilepath)/$(string(lbll))mcmc.gif"
end

function drawOneMC!(cliqMC::CliqGibbsMC, minmax, mcmc=0; offs=2.0)

    if mcmc>minmax[4]
      minmax[4]=mcmc
    end

    i = 0.0
    for prod in cliqMC.prods
        prodVal = kde!(prod.product,"lcv")
        plotKDEProd!([prodVal;prod.potentials],minmax, h=-i*offs, mcmc=mcmc)
        i += 1.0
    end

end

function drawMCMCDebug(cliq; offs=2.0)
    println("$(cliq.attributes["label"])")
    cliqDbg = cliq.attributes["data"].debug
    MCd = 100.0/length(cliqDbg.mcmc)
    MCo=0.0
    minmax=[99999,-99999,-10,-99999.0]
    for mc in cliqDbg.mcmc
        drawOneMC!(mc, minmax, MCo, offs=offs)
        MCo+=MCd
    end

    n = 3
    y = linspace(minmax[1], minmax[2], n)
    x = linspace(minmax[3],minmax[4],n)
    xgrid = repmat(x',n,1)
    ygrid = repmat(y,1,n)
    z = zeros(n,n)

    for i in 1:length(cliqDbg.mcmc[1].prods)
      surf(xgrid,ygrid,z-(i-1)*offs,alpha=0.04, linewidth=0.0)
    end
end

function drawTreeUpwardMsgs(fgl::G,
                            bt::BayesTree;
                            N=300 ) where G <: AbstractDFG
  #
  len = length(bt.cliques)-1
  vv = Array{Gadfly.Compose.Context,1}(undef, len)
  #r = Array{RemoteRef,1}(len)
  i = 0
  for cliq in bt.cliques
      if cliq[1] == 1 println("No upward msgs from root."); continue; end
      @show cliq[2].attributes["label"]
      i+=1
      vv[i] = drawHorDens(fgl, cliq[2].attributes["debug"].outmsg.p, N)
  end
  #[r[j] = @spawn drawCliqueMsgs(bt.cliques[j+1]) for j in 1:len]
  #[vv[j] = fetch(r[j]) for j in 1:len]
  vv
end

function drawFrontalDens(fg::G,
                         bt::BayesTree;
                         N=300,
                         gt=Union{} ) where G <: AbstractDFG
    #
    len = length(bt.cliques)
    vv = Array{Gadfly.Compose.Context,1}(len)
    i = 0
    for cliq in bt.cliques
        #@show cliq[2].attributes["label"]
        lenfr = length(cliq[2].attributes["data"].frontalIDs)

        p = Array{BallTreeDensity,1}(lenfr)
        j=0
        #pvals = Array{Array{Float64,2},1}(lenfr)
        gtvals = Dict{Int,Array{Float64,2}}()
        lbls = String[]

        for frid in cliq[2].attributes["data"].frontalIDs
            j+=1
            p[j] = getKDE(getVariable(fg, frid)) # getKDE(fg.v[frid])
            # p[j] = kde!(fg.v[frid].attributes["val"])

            #pvals[j] = fg.v[frid].attributes["val"]

            if gt!=Union{}
              gtvals[j] = gt[getVariable(fg,frid).label] # fg.v[frid].
              #push!(gtvals, gt[fg.v[frid].attributes["label"]][1])
              #push!(gtvals, gt[fg.v[frid].attributes["label"]][2])
            end
            push!(lbls, getVariable(fg,frid).label) # fg.v[frid].

        end

        #r = Array{RemoteRef,1}(lenfr)
        #[r[j] = @spawn kde!(pvals[j]) for j in 1:lenfr]
        #[p[j] = fetch(r[j]) for j in 1:lenfr]

        i+=1
        if length(gtvals) > 0
          #gtvals = reshape(gtvals,2,round(Int,length(gtvals)/2))'
          vv[i] = drawHorDens(p, N, gt=gtvals, lbls=lbls)
        else
          vv[i] = drawHorDens(p, N,lbls=lbls)
        end
    end

    #r = Array{RemoteRef,1}(lenfr)
    #[r[j] = @spawn kde!(pvals[j]) for j in 1:lenfr]
    #[p[j] = fetch(r[j]) for j in 1:lenfr]

    i+=1
    if length(gtvals) > 0
      #gtvals = reshape(gtvals,2,round(Int,length(gtvals)/2))'
      vv[i] = drawHorDens(p, N, gt=gtvals, lbls=lbls)
    else
      vv[i] = drawHorDens(p, N,lbls=lbls)
    end
  #
  return vv
end

# for some reason we still need msgPlots of right size in the global for these functions to work.
# precall drawTreeUpwardMsgs or drawFrontalDens to make this work properly TODO
function vstackedDensities(msgPlots)
    #msgPlots = f(fg, bt) # drawTreeUpwardMsgs
    evalstr = ""
    for i in 1:length(msgPlots)
        evalstr = string(evalstr, ",msgPlots[$(i)]")
    end
    eval(parse(string("vstack(",evalstr[2:end],")")))
end

function investigateMultidimKDE(p::BallTreeDensity, p0::BallTreeDensity)
    co = ["black"; "blue"]
    h = Union{}
    x = plotKDE([marginal(p,[1]); marginal(p0,[1])], c=co )
    y = plotKDE([marginal(p,[2]); marginal(p0,[2])], c=co )
    if p.bt.dims >= 3
      th = plotKDE([marginal(p,[3]); marginal(p0,[3])], c=co )
      h = hstack(x,y,th)
    else
      h = hstack(x,y)
    end

    return h
end

function investigateMultidimKDE(p::Array{BallTreeDensity,1})
    co = ["black"; "blue"; "green"; "red"; "magenta"; "cyan"; "cyan1"; "cyan2";
    "magenta"; "cyan"; "cyan1"; "cyan2"; "magenta"; "cyan"; "cyan1"; "cyan2"; "magenta";
    "cyan"; "cyan1"; "cyan2"; "magenta"; "cyan"; "cyan1"; "cyan2"]
    # compute all the marginals
    Pm = Array{Array{BallTreeDensity,1},1}()
    push!(Pm,stackMarginals(p,1)) #[marginal(p[1],[1]); marginal(p[2],[1])]
    push!(Pm,stackMarginals(p,2)) #[marginal(p[1],[2]); marginal(p[2],[2])]

    h = Union{}
    x = plotKDE(Pm[1], c=co )
    y = plotKDE(Pm[2], c=co )
    if p[1].bt.dims >= 3
      #Pm3 = [marginal(p[1],[3]); marginal(p[2],[3])]
      push!(Pm,stackMarginals(p,3)) # [marginal(p[1],[3]); marginal(p[2],[3])]
      th = plotKDE(Pm[3], c=co )
      h = hstack(x,y,th)
    else
      h = hstack(x,y)
    end
    return h
end

function investigateMultidimKDE(p::BallTreeDensity)
    x = plotKDE(marginal(p,[1]) )
    y = plotKDE(marginal(p,[2]) )
    if p.bt.dims >= 3
      th = plotKDE(marginal(p,[3]) )
      return hstack(x,y,th)
    end
    return hstack(x,y)
end

function vArrPotentials(potens::Dict{Symbol,EasyMessage})
  vv = Array{Gadfly.Compose.Context,1}(length(potens))
  i = 0
  oned=false
  for p in potens
      i+=1
      pb = kde!(p[2].pts, p[2].bws)
      if size(p[2].pts,1) > 3
        # vv[i] = plotKDE(pb)
        error("can't handle higher dimensional plots here yet")
      elseif size(p[2].pts,1) > 1
        vv[i] = investigateMultidimKDE(pb)
      else
        vv[i] = plotKDE(pb)
      end
  end
  return vv
end

function draw(em::EasyMessage;xlbl="X")
  p = Union{}
  if size(em.pts,1) == 1
    p=plotKDE(kde!(em),xlbl=xlbl)
  else
    p=plotKDE(kde!(em))
  end
  return p
end

function whosWith(cliq::Graphs.ExVertex)
  println("$(cliq.attributes["label"])")
  for pot in cliq.attributes["data"].potentials
      println("$(pot)")
  end
  nothing
end

function whosWith(bt::BayesTree, frt::String)
    whosWith(whichCliq(bt,frt))
end

function drawUpMsgAtCliq(fg::G,
                         cliq::Graphs.ExVertex  ) where G <: AbstractDFG
  #
  for id in keys(cliq.attributes["data"].debug.outmsg.p)
      print("$(getVariable(fg,id).label), ") #fg.v[id].
  end
  println("")
  sleep(0.1)
  potens = getData(cliq).debug.outmsg.p
  vArrPotentials(potens)
end

function drawUpMsgAtCliq(fg::G,
                         bt::BayesTree,
                         lbl::Union{String,Symbol}  ) where G <: AbstractDFG
  #
  drawUpMsgAtCliq(fg, whichCliq(bt, Symbol(lbl)) )
end

# function drawDwnMsgAtCliq(fg::FactorGraph, cliq::Graphs.ExVertex)
function dwnMsgsAtCliq(fg::G,
                       cliq::Graphs.ExVertex  ) where G <: AbstractDFG
  #
  for id in keys(cliq.attributes["data"].debugDwn.outmsg.p)
      print("$(getVariable(fg,id).label), ") # fg.v[id].
  end
  println("")
  sleep(0.1)
  potens = cliq.attributes["data"].debugDwn.outmsg.p
  potens
end

function dwnMsgsAtCliq(fg::G,
                       bt::BayesTree,
                       lbl::Union{String,Symbol}  ) where G <: AbstractDFG
  #
  dwnMsgsAtCliq(fg, whichCliq(bt, Symbol(lbl)) )
end

function drawPose2DMC!(plots::Array{Gadfly.Compose.Context,1},
                       cliqMC::CliqGibbsMC )
  #
  for prod in cliqMC.prods
    prodVal = kde!(prod.product, "lcv") #cliqMC.prods[1]
    push!(plots, plotKDE([prodVal;prod.potentials]) )
  end
  vstackedPlots(plots)
end

function mcmcPose2D!(plots::Array{Gadfly.Compose.Context,1},
                     cliqDbg::DebugCliqMCMC,
                     iter::Int=1  )
    # for mc in cliqDbg.mcmc
    mc = cliqDbg.mcmc[iter]
    v = drawPose2DMC!(plots, mc)
    # end
    return v
end

function drawUpMCMCPose2D!(plots::Array{Gadfly.Compose.Context,1},
                           cliq::Graphs.ExVertex,
                           iter::Int=1 )
  #
  whosWith(cliq)
  cliqDbg = cliq.attributes["data"].debug
  sleep(0.1)
  mcmcPose2D!(plots, cliqDbg, iter)
end

function drawUpMCMCPose2D!(plots::Array{Gadfly.Compose.Context,1},
                           bt::BayesTree,
                           frt::String,
                           iter::Int=1 )
  #
  drawUpMCMCPose2D!(plots, whichCliq(bt,frt), iter)
end

function drawDwnMCMCPose2D!(plots::Array{Gadfly.Compose.Context,1},
                            cliq::Graphs.ExVertex,
                            iter::Int=1  )
  #
  whosWith(cliq)
  cliqDbg = cliq.attributes["data"].debugDwn
  sleep(0.1)
  mcmcPose2D!(plots, cliqDbg, iter)
end

function drawDwnMCMCPose2D!(plots::Array{Gadfly.Compose.Context,1},
                            bt::BayesTree,
                            frt::String,
                            iter::Int=1 )
  #
  drawDwnMCMCPose2D!(plots, whichCliq(bt,frt), iter)
end

function drawLbl(fgl::G, lbl::Symbol) where G <: AbstractDFG
    plotKDE(getKDE(getVariable(fgl,lbl)))
end
drawLbl(fgl::G, lbl::T) where {G <: AbstractDFG, T <: AbstractString} = drawLbl(fgl, Symbol(lbl))

function predCurrFactorBeliefs(fgl::G,
                               fc::Graphs.ExVertex  ) where G <: AbstractDFG
  #
  # TODO update to use ls and lsv functions
  prjcurvals = Dict{String, Array{BallTreeDensity,1}}()
  for v in getNeighbors(fgl, fc)
    pred = kde!(evalFactor2(fgl, fc, v.index))
    curr = kde!(getVal(v))
    prjcurvals[v.attributes["label"]] = [curr; pred]
  end
  return prjcurvals, collect(keys(prjcurvals))
end

function drawHorDens(fgl::G,
                     pDens::Dict{Int,EasyMessage},
                     N=200 ) where G <: AbstractDFG
  #
  p = BallTreeDensity[]
  lbls = String[]
  for pd in pDens
    push!(p, kde!(pd[2].pts,pd[2].bws))
    push!(lbls, getVariable(fgl,pd[1]).label)
  end
  drawHorDens(p,N=N,lbls=lbls)
end

function drawHorBeliefsList(fgl::G,
                            lbls::Array{Symbol,1};
                            nhor::Int=-1,
                            gt=nothing,
                            N::Int=200,
                            extend=0.1  ) where G <: AbstractDFG
  #
  len = length(lbls)
  pDens = BallTreeDensity[]
  for lb in lbls
    ptkde = getKDE(getVariable(fgl,lb))
    push!(pDens, ptkde )
  end

  if nhor<1
    nhor = round(Int,sqrt(len))
  end
  vlen = ceil(Int, len/nhor)
  vv = Array{Gadfly.Compose.Context,1}(vlen)
  conslb = deepcopy(lbls)
  vidx = 0
  for i in 1:nhor:len
    pH = BallTreeDensity[]
    gtvals = Dict{Int,Array{Float64,2}}()
    labels = String[]
    for j in 0:(nhor-1)
      if i+j <= len
        push!(pH, pDens[i+j])
        push!(labels, string(lbls[i+j]))
        if gt != nothing gtvals[j+1] = gt[lbls[i+j]]  end
      end
    end
    vidx+=1
    if gt !=nothing
      vv[vidx] = KernelDensityEstimatePlotting.drawHorDens(pH, N=N, gt=gtvals, lbls=labels, extend=extend)
    else
      vv[vidx] = KernelDensityEstimatePlotting.drawHorDens(pH, N=N, lbls=labels, extend=extend)
    end
  end
  vv
end

function drawFactorBeliefs(fgl::G,
                           flbl::Symbol ) where G <: AbstractDFG
  #
  if !haskey(fgl.fIDs, flbl)
    println("no key $(flbl)")
    return nothing
  end
  # for fc in fgl.f
    # if fc[2].attributes["label"] == flbl

    fc = fgl.g.vertices[fgl.fIDs[flbl]]  # fc = fgl.f[fgl.fIDs[flbl]]
      prjcurvals, lbls = predCurrFactorBeliefs(fgl, fc)
      if length(lbls) == 3
        return vstack(
        plotKDE(prjcurvals[lbls[1]]),
        plotKDE(prjcurvals[lbls[2]]),
        plotKDE(prjcurvals[lbls[3]]),
        )
      elseif length(lbls) == 2
        return vstack(
        plotKDE(prjcurvals[lbls[1]]),
        plotKDE(prjcurvals[lbls[2]]),
        )
      elseif length(lbls) == 1
        return plotKDE(prjcurvals[lbls[1]])
      end

    # end
  # end
  nothing
end
drawFactorBeliefs(fgl::G, flbl::T) where {G <: AbstractDFG, T <: AbstractString} = drawFactorBeliefs(fgl, Symbol(flbl))

"""
    $(SIGNATURES)

Plot the proposal belief from neighboring factors to `lbl` in the factor graph (ignoring Bayes tree representation),
and show with new product approximation for reference.
"""
function plotLocalProduct(fgl::G,
                          lbl::Symbol;
                          N::Int=100,
                          dims::Vector{Int}=Int[],
                          levels::Int=1,
                          show::Bool=false,
                          dirpath="/tmp/",
                          mimetype::AbstractString="svg",
                          sidelength=20cm,
                          title::String="Local product, "  ) where G <: AbstractDFG
  #
  @warn "not showing partial constraints, but included in the product"
  arr = Array{BallTreeDensity,1}()
  lbls = String[]
  push!(arr, getKDE(getVariable(fgl, lbl)))
  push!(lbls, "curr")
  pl = nothing
  pp, parr, partials, lb = IncrementalInference.localProduct(fgl, lbl, N=N)
  if length(parr) > 0 && length(partials) == 0
    if pp != parr[1]
      push!(arr,pp)
      push!(lbls, "prod")
      for a in parr
        push!(arr, a)
      end
      @show lb, lbls
      lbls = union(lbls, string.(lb))
    end
    dims = length(dims) > 0 ? dims : collect(1:Ndim(pp))
    colors = getColorsByLength(length(arr))
    pl = plotKDE(arr, dims=dims, levels=levels, c=colors, title=string(title,lbl), legend=string.(lbls)) #
  elseif length(parr) == 0 && length(partials) > 0
    # stack 1d plots to accomodate all the partials
    PL = []
    lbls = String["prod";"curr";string.(lb)]
    pdims = sort(collect(keys(partials)))
    for dimn in pdims
      vals = partials[dimn]
      proddim = marginal(pp, [dimn])
      colors = getColorsByLength(length(vals)+2)
      pl = plotKDE([proddim;getKDE(getVariable(fgl, lbl));vals], dims=[1;], levels=levels, c=colors, title=string("Local product, dim=$(dimn), ",lbl))
      push!(PL, pl)
    end
    pl = Gadfly.vstack(PL...)
  else
    return error("plotLocalProduct not built for lengths parr, partials = $(length(parr)), $(length(partials)) yet.")
  end

  # now let's export:
  backend = getfield(Gadfly, Symbol(uppercase(mimetype)))
  Gadfly.draw(backend(dirpath*"test_$(lbl).$(mimetype)",sidelength,0.75*sidelength), pl)
  driver = mimetype in ["pdf"] ? "evince" : "eog"
  show ? (@async run(`$driver $(dirpath)test_$(lbl).$(mimetype)`)) : nothing

  return pl
end

"""
    $(SIGNATURES)

Plot---for the cylindrical manifold only---the proposal belief from neighboring factors to `lbl` in the factor graph (ignoring Bayes tree representation),
and show with new product approximation for reference.  The linear and circular belief products are returned as a tuple.
"""
function plotLocalProductCylinder(fgl::G,
                                  lbl::Symbol;
                                  N::Int=100,
                                  levels::Int=1,
                                  show=false,
                                  dirpath="/tmp/",
                                  mimetype::AbstractString="svg",
                                  sidelength=30cm,
                                  scale::Float64=0.2 ) where G <: AbstractDFG
  #
  @warn "not showing partial constraints, but included in the product"
  arr = Array{BallTreeDensity,1}()
  carr = Array{BallTreeDensity,1}()
  lbls = String[]
  push!(arr, getKDE(getVariable(fgl, lbl)))
  push!(lbls, "curr")
  pl = nothing
  plc = nothing
  pp, parr, partials, lb = IncrementalInference.localProduct(fgl, lbl, N=N)
  if length(parr) > 0 && length(partials) == 0
    if pp != parr[1]
      push!(arr, marginal(pp,[1]))
      push!(lbls, "prod")
      push!(carr, marginal(pp, [2]))
      for a in parr
        push!(arr, marginal(a,[1]))
        push!(carr, marginal(a,[2]))
      end
      lbls = union(lbls, lb)
    end
    colors = getColorsByLength(length(arr))
    pll = plotKDE(arr, dims=[1;], levels=levels, c=colors, legend=lbls, title=string("Local product, ",lbl))
    plc = AMP.plotKDECircular(carr, c=colors, scale=scale)
    # (pll, plc)
  elseif length(parr) == 0 && length(partials) > 0
    # stack 1d plots to accomodate all the partials
    PL = []
    PLC = []
    lbls = ["prod";"curr";lb]
    pdims = sort(collect(keys(partials)))
    for dimn in pdims
      vals = partials[dimn]
      proddim = marginal(pp, [dimn])
      colors = getColorsByLength(length(vals)+2)
      if dimn == 1
        pll = plotKDE([proddim;getKDE(getVariable(fgl, lbl));vals], dims=[1;], levels=levels, c=colors, title=string("Local product, dim=$(dimn), ",lbl))
        push!(PL, pl)
      else
        plc = AMP.plotKDECircular([proddim; marginal(getKDE(getVariable(fgl, lbl)),[2]);vals], c=colors, scale=scale)
        push!(PLC, plc)
      end
    end
    pl = Gadfly.vstack(PL..., PLC...)
  else
    return error("plotLocalProduct not built for lengths parr, partials = $(length(parr)), $(length(partials)) yet.")
  end

  # now let's export:
  # backend = getfield(Gadfly, Symbol(uppercase(mimetype)))
  # Gadfly.draw(backend(dirpath*"test_$(lbl).$(mimetype)",sidelength,sidelength), pl)
  # driver = mimetype in ["pdf"] ? "evince" : "eog"
  # show ? (@async run(`$driver $(dirpath)test_$(lbl).$(mimetype)`)) : nothing

  return pll, plc
end

"""
    $SIGNATURES

Plot the proposal belief from neighboring factors to `lbl` in the factor graph (ignoring Bayes tree representation),
and show with new product approximation for reference. String version is obsolete and will be deprecated.
"""
plotLocalProduct(fgl::G, lbl::T; N::Int=100, dims::Vector{Int}=Int[]) where {G <: AbstractDFG, T <: AbstractString} = plotLocalProduct(fgl, Symbol(lbl), N=N, dims=dims)

"""
    $SIGNATURES

Project (convolve) to and take product of variable in Bayes/Junction tree.

Notes
- assume cliq and var sym the same, unless both specified.
- `cliqsym` defines a frontal variable of a clique.
"""
function plotTreeProductUp(fgl::G,
                           treel::BayesTree,
                           cliqsym::Symbol,
                           varsym::Symbol=cliqsym  ) where G <: AbstractDFG
  #
  # build a subgraph copy of clique
  cliq = whichCliq(treel, cliqsym)
  syms = getCliqAllVarIds(cliq)
  subfg = buildSubgraphFromLabels(fgl,syms)

  # add upward messages to subgraph
  msgs = getCliqChildMsgsUp(treel,cliq, BallTreeDensity)
  # @show typeof(msgs)
  addMsgFactors!(subfg, msgs)

  # predictBelief
  # stuff = treeProductUp(fgl, treel, cliqsym, varsym)
  # plotKDE(manikde!(stuff[1], getManifolds(fgl, varsym)))
  cllbl = cliq.attributes["label"]
  return plotLocalProduct(subfg, varsym, title="Tree Up $(cllbl) | ")
end

function plotTreeProductDown(fgl::G,
                             treel::BayesTree,
                             cliqsym::Symbol,
                             varsym::Symbol=cliqsym;
                             levels::Int=1  ) where G <: AbstractDFG
  #
  # build a subgraph copy of clique
  cliq = whichCliq(treel, cliqsym)
  syms = getCliqAllVarIds(cliq)
  subfg = buildSubgraphFromLabels(fgl,syms)

  # add upward messages to subgraph
  msgs = getCliqParentMsgDown(treel,cliq)
  addMsgFactors!(subfg, msgs)

  # predictBelief
  # stuff = treeProductUp(fgl, treel, cliqsym, varsym)
  # plotKDE(manikde!(stuff[1], getManifolds(fgl, varsym)))
  cllbl = cliq.attributes["label"]
  return plotLocalProduct(subfg, varsym, title="Tree Dwn $(cllbl) | ", levels=levels)
end

function saveplot(pl;name="pl",frt=:png,w=25cm,h=25cm,nw=false,fill=true)
  if frt==:png
    Gadfly.draw(PNG(string(name,".png"),w,h),pl)
    # if fill run(`composite $(name).png plB.png $(name).png`) end
    if !nw run(`eog $(name).png`) end
  end
  if frt==:pdf
    Gadfly.draw(PDF(string(name,".pdf"),w,h),pl)
    if !nw run(`evince $(name).pdf`) end
  end
  nothing
end

function animateVertexBelief(FGL::Array{<:AbstractDFG,1}, lbl; nw=false)
  len = length(FGL)
  [saveplot(plotLocalProduct(FG[i],lbl),h=15cm,w=30cm,name="gifs/pl$(i)",nw=true) for i=1:len];
  run(`convert -delay 100 gifs/pl'*'.png result.gif`)
  if !nw run(`eog result.gif`) end
  nothing
end

function fixRotWrapErr!(RT::Array{Float64,1})

  for i in 1:length(RT)
    if RT[i] > pi
      RT[i] = abs(RT[i]-2.0*pi)
    end
  end
  nothing
end

function asyncUniComp(fgl::G, isamdict::Dict{Int,Array{Float64,1}}) where G <: AbstractDFG
  r,rt,lb = computeGraphResiduals(fgl,isamdict);
  fixRotWrapErr!(rt)
  return [sqrt(mean(r.^2));maximum(abs(r));sqrt(mean(rt.^2));maximum(rt)]
end

function unimodalCompare(FGL::Array{<:AbstractDFG,1},isamdict::Dict{Int,Array{Float64,1}})
  len = length(FGL)
  RMS = Float64[]
  MAX = Float64[]
  RMSth = Float64[]
  MAXth = Float64[]

  rr = Future[] #RemoteRef[]

  for fgl in FGL
    push!(rr, remotecall(uppA(),asyncUniComp, fgl, isamdict))
  end

  for r in rr
    err = fetch(r)
    push!(RMS, err[1])
    push!(MAX, err[2])
    push!(RMSth, err[3])
    push!(MAXth, err[4])
  end

  x=0:(len-1)
  df1 = DataFrame(x=x, y=RMS, label="rms")
  df2 = DataFrame(x=x, y=MAX, label="max")
  df3 = DataFrame(x=x, y=RMSth*180.0/pi, label="rmsth")
  df4 = DataFrame(x=x, y=MAXth*180.0/pi, label="maxth")
  df = vcat(df1, df2)
  dfth = vcat(df3,df4)

  return df,dfth
end

function asyncAnalyzeSolution(fgl::G, sym::Symbol) where G <: AbstractDFG
  lbl = string(sym)
  pp, arr, partials = IncrementalInference.localProduct(fgl, lbl)
  lpm = getKDEMax(pp)
  em = getKDEMax(getKDE(getVariable(fgl,lbl)))
  err1 = norm(lpm[1:2]-em[1:2])
  err2 = 0.0
  if lbl[1]=='x'
    err2 = abs(lpm[3]-em[3])
  end
  return [err1;err2]
end

function analyzeSolution(FGL::Array{<: AbstractDFG,1},fggt=Union{})
  len = length(FGL)
  RMS = Float64[]
  MAX = Float64[]
  RMSth = Float64[]
  MAXth = Float64[]
  for fgl in FGL
    xLB, lLB = ls(fgl)
    ERR = Float64[]
    ERRth = Float64[]
    ALB = [xLB;lLB]
    rr = Future[] #RemoteRef[]
    for lbl in ALB
      push!(rr, remotecall(uppA(),asyncAnalyzeSolution, fgl, lbl))
      # err
      # push!(ERR, err[1])
      # if lbl[1]=='x'
      #   push!(ERRth, err[2])
      # end
    end

    idx = 1
    for r in rr
      err = fetch(r)
      push!(ERR, err[1])
      if ALB[idx][1]=='x'
        push!(ERRth, err[2])
      end
      idx += 1
    end
    push!(RMS, sqrt(mean(ERR.^2)))
    push!(MAX, maximum(abs(ERR)))
    push!(RMSth, sqrt(mean(ERRth.^2)))
    push!(MAXth, maximum(ERRth))
  end

  x=0:(len-1)
  df1 = DataFrame(x=x, y=RMS, label="rms")
  df2 = DataFrame(x=x, y=MAX, label="max")
  df3 = DataFrame(x=x, y=RMSth*180.0/pi, label="rmsth")
  df4 = DataFrame(x=x, y=MAXth*180.0/pi, label="maxth")
  df = vcat(df1, df2)
  dfth = vcat(df3,df4)
  return df,dfth
end
# discrete_color_manual(colors...; levels=nothing,order=nothing) is deprecated, use color_discrete_manual(colors...; levels=levels,order=order) instead.

function drawAnalysis(df,dfth)
  return vstack(
  Gadfly.plot(df, x="x", y="y", color="label", Geom.line,
         Scale.discrete_color_manual("red","black")),
  Gadfly.plot(dfth, x="x", y="y", color="label", Geom.line,
        Scale.discrete_color_manual("red","black"))
        )
end

function getAllFGsKDEs(fgD::Array{<: AbstractDFG,1}, vertsym::Symbol)
  ret = Array{BallTreeDensity,1}()
  for i in 1:length(fgD)
    push!(ret, getKDE(getVariable(fgD[i],vertsym)) )
  end
  return ret
end

function drawAllPose2DBeliefs(plots::Array{Gadfly.Compose.Context,1}, fgD::Array{<: AbstractDFG,1})
    ids = sort(collect(keys(fgD[1].v)))
    co = ["black"; "blue"; "green"; "red"; "magenta"; "cyan"; "cyan1"; "cyan2"]
    println(co[1:length(fgD)])
    for i in ids
        getVariable(fgD[1],i).label
        kdes = getAllFGsKDEs(fgD, i)
        push!(plots, plotKDE(  kdes  )) # [kde!(getVal(V)); kde!(getVal(V0))]
    end
    vstackedPlots(plots)
end

# legacy function -- use the array version instead
function drawAllPose2DBeliefs(plots::Array{Gadfly.Compose.Context,1}, fgD::G, fgD0=Union{}) where G <: AbstractDFG
  println("WARNING: drawAllPose2DBeliefs -- legacy function -- use the array version instead.")
  if fgD0 != Union{}
    drawAllPose2DBeliefs(plots, [fgD;fgD0])
  else
    drawAllPose2DBeliefs(plots, [fgD])
  end
end

function drawComicStripLM(fgD::Array{<: AbstractDFG,1})
    comicA = Array{Gadfly.Plot,1}()
    for fgd in fgD
        cv = drawPosesLandms(fgd)
        # cv = drawPoses(fgd)
        push!(comicA,cv)
    end
    hstack(comicA)
end

function drawComicStrip(fgD::Array{<: AbstractDFG,1})
    comicA = Array{Gadfly.Plot,1}()
    for fgd in fgD
        cv = drawPoses(fgd)
        push!(comicA,cv)
    end
    hstack(comicA)
end

function compositeComic(fnc::Function, fgGT, fgA::Array{<: AbstractDFG,1})
    v = Union{}
    @show length(fgA)
    if length(fgA) == 2
        Gadfly.set_default_plot_size(25cm, 10cm)
        v = fnc([fgA[1:2];fgGT])
    elseif length(fgA) == 3
        Gadfly.set_default_plot_size(25cm, 20cm)
        v = vstack(fnc(fgA[1:2])
        ,fnc([fgA[3];fgGT])    )
    elseif length(fgA) == 4
        Gadfly.set_default_plot_size(25cm, 20cm)
        v = vstack(fnc(fgA[1:3])
        ,fnc([fgA[4];fgGT])    )
    elseif length(fgA) == 7
        Gadfly.set_default_plot_size(25cm, 25cm)
        v = vstack(fnc(fgA[1:3])
        ,fnc(fgA[4:6])
        ,fnc([fgA[7];fgGT])    )
    elseif length(fgA) == 10
        Gadfly.set_default_plot_size(25cm, 25cm)
        v = vstack(fnc(fgA[1:3])
        ,fnc(fgA[4:6])
        ,fnc(fgA[7:9])
        ,fnc([fgA[10];fgGT])    )
    elseif length(fgA) == 13
        Gadfly.set_default_plot_size(25cm, 30cm)
        v = vstack(fnc(fgA[1:3])
        ,fnc(fgA[4:6])
        ,fnc(fgA[7:9])
        ,fnc(fgA[10:12])
        ,fnc([fgA[13];fgGT])    )
    end
    v
end

function compositeComic(fnc::Function, fgA::Array{<: AbstractDFG,1})
    v = Union{}
    @show length(fgA)
    if length(fgA) == 2
        Gadfly.set_default_plot_size(25cm, 10cm)
        v = fnc(fgA[1:2])
    elseif length(fgA) == 3
        Gadfly.set_default_plot_size(25cm, 20cm)
        v = fnc(fgA[1:3])
    elseif length(fgA) == 4
        Gadfly.set_default_plot_size(25cm, 25cm)
        v = vstack(fnc(fgA[1:2])
        ,fnc(fgA[3:4]))
    elseif length(fgA) == 6
        Gadfly.set_default_plot_size(25cm, 25cm)
        v = vstack(fnc(fgA[1:3])
        ,fnc(fgA[4:6])  )
    elseif length(fgA) == 9
        Gadfly.set_default_plot_size(25cm, 25cm)
        v = vstack(fnc(fgA[1:3])
        ,fnc(fgA[4:6])
        ,fnc(fgA[7:9])    )
    elseif length(fgA) == 12
        Gadfly.set_default_plot_size(25cm, 30cm)
        v = vstack(fnc(fgA[1:3])
        ,fnc(fgA[4:6])
        ,fnc(fgA[7:9])
        ,fnc(fgA[10:12])    )
    end
    v
end

#
#
#
# function spyCliqMat(cliq::Graphs.ExVertex; showmsg=true)
#   mat = deepcopy(getCliqMat(cliq, showmsg=showmsg))
#   # TODO -- add improved visualization here, iter vs skip
#   mat = map(Float64, mat)*2.0.-1.0
#   numlcl = size(IIF.getCliqAssocMat(cliq),1)
#   mat[(numlcl+1):end,:] *= 0.9
#   mat[(numlcl+1):end,:] .-= 0.1
#   numfrtl1 = floor(Int,length(getData(cliq).frontalIDs) + 1)
#   mat[:,numfrtl1:end] *= 0.9
#   mat[:,numfrtl1:end] .-= 0.1
#   @show getData(cliq).itervarIDs
#   @show getData(cliq).directvarIDs
#   @show getData(cliq).msgskipIDs
#   @show getData(cliq).directFrtlMsgIDs
#   @show getData(cliq).directPriorMsgIDs
#   sp = Gadfly.spy(mat)
#   push!(sp.guides, Gadfly.Guide.title("$(cliq.attributes["label"]) || $(cliq.attributes["data"].frontalIDs) :$(cliq.attributes["data"].conditIDs)"))
#   push!(sp.guides, Gadfly.Guide.xlabel("fmcmcs $(cliq.attributes["data"].itervarIDs)"))
#   push!(sp.guides, Gadfly.Guide.ylabel("lcl=$(numlcl) || msg=$(size(getCliqMsgMat(cliq),1))" ))
#   return sp
# end
# function spyCliqMat(bt::BayesTree, lbl::Symbol; showmsg=true)
#   spyCliqMat(whichCliq(bt,lbl), showmsg=showmsg)
# end

function plotPose2Vels(fgl::G, sym::Symbol; coord=nothing) where G <: AbstractDFG
  X = getKDE(getVariable(fgl, sym))
  px = plotKDE(X, dims=[4], title="Velx")
  coord != nothing ? (px.coord = coord) : nothing
  py = plotKDE(X, dims=[5], title="Vely")
  coord != nothing ? (py.coord = coord) : nothing
  hstack(px, py)
end

"""
    $(SIGNATURES)

Analysis function to compare KDE plots between the factor graph centric product of a variable with
current value stored in the factor graph object.
"""
function plotProductVsKDE(fgl::G,
                          sym::Symbol;
                          levels::Int=3,
                          c::Vector{String}=["red";"black"] ) where  G <: AbstractDFG
    #
    plotKDE([IIF.localProduct(fgl, sym)[1], getKDE(getVariable(fgl, sym))], levels=3, c=c)
end

"""
    $SIGNATURES

Overlay plot all upward messages from cliques.
"""
function plotCliqUpMsgs(fg::G,
                        tree::BayesTree,
                        sym::Symbol;
                        show::Bool=true,
                        dims::Vector{Int}=Int[],
                        levels::Int=1,
                        c=nothing,
                        title="up msgs on $(sym)"    ) where G <: AbstractDFG
  #
  # get all msgs
  allmsgs = getTreeCliqUpMsgsAll(tree)
  # stack messages by variable
  sckmsgs = stackCliqUpMsgsByVariable(tree, allmsgs)

  if !haskey(sckmsgs, sym)
    @warn "plotCliqUpMsgs -- tree does not have up messages for $sym."
    return nothing
  end

  # prepend current estimate too
  Xs = getKDE(fg, sym)
  # vectorize beliefs
  beliefs = BallTreeDensity[Xs;]
  lbls = String["curr,-1";]
  for (frt,dep,bel,infd) in sckmsgs[sym]
    push!(beliefs, bel)
    push!(lbls, "$frt,$dep")
  end

  # ignoring legend and color information
  cc = c == nothing ? getColorsByLength(length(beliefs)) : c

  # plot and return
  plotKDE(beliefs, levels=levels, c=cc, title=title, legend=lbls)
end

"""
    $SIGNATURES

Plot new proposal (convolution) for factor x of variable y given all other -- i.e. y|X

Notes
- plot new result on `towards` as first color.
- Plot other variables in factor on 2nd to m colors.
"""
function plotVariableGivenFactor(dfg::G,
                                 fct::Symbol,
                                 towards::Symbol;
                                 levels::Int=2,
                                 dims=nothing  ) where G <: AbstractDFG
  #
  pts = approxConv(dfg,fct,towards)
  mani = getManifolds(dfg, towards)
  res = manikde!(pts,mani)

  lie = ls(dfg, fct)
  setdiff!(lie, [towards])

  otr = map(x->getKDE(dfg,x),lie)
  lbls = string.([towards;lie])

  pl = plotKDE([res;otr],dims=dims,levels=levels,legend=lbls)

  return pl
end

"""
    $SIGNATURES

Plot the downward messages currently stored in a clique.
"""
function plotCliqDownMsgs(tree::BayesTree,
                          frnt::Symbol;
                          show::Bool=true,
                          levels::Int=2,
                          dims=nothing,
                          existing=nothing  )
  #
  cliq = getCliq(tree,frnt)
  msgs = getCliqMsgsDown(cliq)

  PL = []

  for (key, beldim) in msgs
    npl = plotKDE(beldim[1], levels=levels, title="dwn msg $key", dims=dims)
    existing == nothing ? nothing : union!(npl.layers, existing.layers)
    push!(PL, npl)
  end

  existing == nothing ? nothing : push!(PL, existing)
  pl = vstack(PL...)

  folderpath = "/tmp/caesar/random/"
  filepath = folderpath*"downmsgs_cliq$(cliq.index).pdf"
  Base.mkpath(folderpath)
  pl |> PDF(filepath, 20cm, length(PL)*12cm)
  @async run(`evince $filepath`)
  pl
end


text:
dir_list = [
"0_5-1_25_050-075_bin",
"0_5-1_25_075-100_bin",
"0_5-1_25_100-125_bin",
"0_5-1_25_125-150_bin",
"0_5-1_25_150-175_bin",
"0_5-1_25_175-200_bin",
"0_5-1_25_2-2_5_bin",
"0_5-1_25_2_5-3_bin",
"0_5-1_25_3-4_bin",
"0_5-1_25_4-6_bin",
"0_5-1_25_6-8_bin",
"0_5-1_25_8-12_bin",
"0_5-1_25_12-16_bin",
"1_25-2_5_050-075_bin",
"1_25-2_5_075-100_bin",
"1_25-2_5_100-125_bin",
"1_25-2_5_125-150_bin",
"1_25-2_5_150-175_bin",
"1_25-2_5_175-200_bin",
"1_25-2_5_2-2_5_bin",
"1_25-2_5_2_5-3_bin",
"1_25-2_5_3-4_bin",
"1_25-2_5_4-6_bin",
"1_25-2_5_6-8_bin",
"1_25-2_5_8-12_bin",
"1_25-2_5_12-16_bin",
"2_5-5_050-075_bin",
"2_5-5_075-100_bin",
"2_5-5_100-125_bin",
"2_5-5_125-150_bin",
"2_5-5_150-175_bin",
"2_5-5_175-200_bin",
"2_5-5_2-2_5_bin",
"2_5-5_2_5-3_bin",
"2_5-5_3-4_bin",
"2_5-5_4-6_bin",
"2_5-5_6-8_bin",
"2_5-5_8-12_bin",
"2_5-5_12-16_bin",
"5-10_050-075_bin",
"5-10_075-100_bin",
"5-10_100-125_bin",
"5-10_125-150_bin",
"5-10_150-175_bin",
"5-10_175-200_bin",
"5-10_2-2_5_bin",
"5-10_2_5-3_bin",
"5-10_3-4_bin",
"5-10_4-6_bin",
"5-10_6-8_bin",
"5-10_8-12_bin",
"5-10_12-16_bin",
"10-20_050-075_bin",
"10-20_075-100_bin",
"10-20_100-125_bin",
"10-20_125-150_bin",
"10-20_150-175_bin",
"10-20_175-200_bin",
"10-20_2-2_5_bin",
"10-20_2_5-3_bin",
"10-20_3-4_bin",
"10-20_4-6_bin",
"10-20_6-8_bin",
"10-20_8-12_bin",
"10-20_12-16_bin",
"20-40_050-075_bin",
"20-40_075-100_bin",
"20-40_100-125_bin",
"20-40_125-150_bin",
"20-40_150-175_bin",
"20-40_175-200_bin",
"20-40_2-2_5_bin",
"20-40_2_5-3_bin",
"20-40_3-4_bin",
"20-40_4-6_bin",
"20-40_6-8_bin",
"20-40_8-12_bin",
"20-40_12-16_bin",
"40-80_050-075_bin",
"40-80_075-100_bin",
"40-80_100-125_bin",
"40-80_125-150_bin",
"40-80_150-175_bin",
"40-80_175-200_bin",
"40-80_2-2_5_bin",
"40-80_2_5-3_bin",
"40-80_3-4_bin",
"40-80_4-6_bin",
"40-80_6-8_bin",
"40-80_8-12_bin",
"40-80_12-16_bin",
"80-160_050-075_bin",
"80-160_075-100_bin",
"80-160_100-125_bin",
"80-160_125-150_bin",
"80-160_150-175_bin",
"80-160_175-200_bin",
"80-160_2-2_5_bin",
"80-160_2_5-3_bin",
"80-160_3-4_bin",
"80-160_4-6_bin",
"80-160_6-8_bin",
"80-160_8-12_bin",
"80-160_12-16_bin",
"160-320_050-075_bin",
"160-320_075-100_bin",
"160-320_100-125_bin",
"160-320_125-150_bin",
"160-320_150-175_bin",
"160-320_175-200_bin",
"160-320_2-2_5_bin",
"160-320_2_5-3_bin",
"160-320_3-4_bin",
"160-320_4-6_bin",
"160-320_6-8_bin",
"160-320_8-12_bin",
"160-320_12-16_bin"
]

p_lim = [0.5, 1.25, 2.5, 5., 10., 20., 40., 80., 160., 320.]
r_lim = [0.5, 0.75, 1., 1.25, 1.5, 1.75, 2., 2.5, 3., 4., 6., 8., 12., 16.]

file_ind = 1
for p_ind = 1:(length(p_lim)-1)
    for r_ind = 1:(length(r_lim)-1)
        if !isdir(dir_list[file_ind])
            mkdir(dir_list[file_ind])
        end
        cd(dir_list[file_ind])
        f_param = open("param.in","w")
        write(f_param,"""stellar_catalog = "q1q17_dr25_gaia_fgk.jld2"\n""")
        write(f_param,"""koi_catalog = "q1_q17_dr25_koi.csv"\n""")
        write(f_param,"num_targ_sim = 5000\n")
        write(f_param,string("p_bin_lim = [",p_lim[p_ind],", ",p_lim[p_ind+1],"]\n"))
        write(f_param,string("r_bin_lim = [",r_lim[r_ind],", ",r_lim[r_ind+1],"]\n"))
        write(f_param,"#rate_init = 1.0")
        file_ind += 1
        cd("..")
    end
end

ext: jl


text:
"""
`Basics.tabulate()`  ... tabulates the data from various objects due to different criteria.

  + `("multiplet: energies", multiplet::Multiplet; stream::IO=stdout)`
    ... to tabulate the energies of all levels of the given multiplet into a neat format; nothing is returned.
  + `("multiplet: energy relative to immediately lower level", multiplet::Multiplet; stream::IO=stdout)`
    ... to tabulate the energy splitting between neighboured levels of all levels of the given multiplet into a neat
        format; nothing is returned.
  + `("multiplet: energy of each level relative to lowest level", multiplet::Multiplet; stream::IO=stdout)`
    ... to tabulate the energy splitting of all levels with regard to the lowest level of the given multiplet into
        a neat format; nothing is returned.
"""
function Basics.tabulate(sa::String, multiplet::Multiplet; stream::IO=stdout)
    if        sa == "multiplet: energies"
        println(stream, "\n  Eigenenergies:")
        sb = "  Level  J Parity          Hartrees       " * "             eV                   " *  JAC.TableStrings.inUnits("energy")
        println(stream, "\n", sb, "\n")
        for  i = 1:length(multiplet.levels)
            lev = multiplet.levels[i]
            en  = lev.energy;    en_eV = Defaults.convertUnits("energy: from atomic to eV", en);    en_requested = Defaults.convertUnits("energy: from atomic", en)
            sc  = " " * JAC.TableStrings.level(i) * "    " * string(LevelSymmetry(lev.J, lev.parity)) * "    "
            @printf(stream, "%s %.15e %s %.15e %s %.15e %s", sc, en, "  ", en_eV, "  ", en_requested, "\n")
        end

    elseif    sa == "multiplet: energy relative to immediately lower level"
        println(stream, "\n  Energy of each level relative to immediately lower level:")
        sb = "  Level  J Parity          Hartrees       " * "             eV                   " * JAC.TableStrings.inUnits("energy")
        println(stream, "\n", sb, "\n")
        for  i = 2:length(multiplet.levels)
            lev = multiplet.levels[i]
            en    = lev.energy - multiplet.levels[i-1].energy;
            en_eV = Defaults.convertUnits("energy: from atomic to eV", en);    en_requested = Defaults.convertUnits("energy: from atomic", en)
            sc  = " " * JAC.TableStrings.level(i) * "    " * string(LevelSymmetry(lev.J, lev.parity)) * "    "
            @printf(stream, "%s %.15e %s %.15e %s %.15e %s", sc, en, "  ", en_eV, "  ", en_requested, "\n")
        end

    elseif    sa == "multiplet: energy of each level relative to lowest level"
        println(stream, "\n  Energy of each level relative to lowest level:")
        sb = "  Level  J Parity          Hartrees       " * "             eV                   " * JAC.TableStrings.inUnits("energy")
        println(stream, "\n", sb, "\n")
        for  i = 2:length(multiplet.levels)
            lev = multiplet.levels[i]
            en    = lev.energy - multiplet.levels[1].energy;
            en_eV = Defaults.convertUnits("energy: from atomic to eV", en);    en_requested = Defaults.convertUnits("energy: from atomic", en)
            sc  = " " * JAC.TableStrings.level(i) * "    " * string(LevelSymmetry(lev.J, lev.parity)) * "    "
            @printf(stream, "%s %.15e %s %.15e %s %.15e %s", sc, en, "  ", en_eV, "  ", en_requested, "\n")
        end
    else
        error("Unsupported keystring.")
    end

    return( nothing )
end

"""
  + `("multiplet: energies", multiplet::JAC.Hfs.IJF_Multiplet; stream::IO=stdout)`
    ... to tabulate the energies of all hyperfine levels of the given multiplet into a neat format; nothing is returned.
  + `("multiplet: energy of each level relative to lowest level", multiplet::JAC.Hfs.IJF_Multiplet; stream::IO=stdout)`
    ... to tabulate the energy splitting of all levels with regard to the lowest level of the given multiplet into
        a neat format; nothing is returned.
"""
function Basics.tabulate(sa::String, multiplet::JAC.Hfs.IJF_Multiplet; stream::IO=stdout)
    if        sa == "multiplet: energies"
        println(stream, "\n  Eigenenergies for nuclear spin I = $(multiplet.levelFs[1].I):")
        sb = "  Level  F Parity          Hartrees       " * "             eV                   " *  JAC.TableStrings.inUnits("energy")
        println(stream, "\n", sb, "\n")
        for  i = 1:length(multiplet.levelFs)
            lev = multiplet.levelFs[i]
            en  = lev.energy;    en_eV = Defaults.convertUnits("energy: from atomic to eV", en);    en_requested = Defaults.convertUnits("energy: from atomic", en)
            sc  = " " * JAC.TableStrings.level(i) * "    " * string(LevelSymmetry(lev.F, lev.parity)) * "    "
            @printf(stream, "%s %.15e %s %.15e %s %.15e %s", sc, en, "  ", en_eV, "  ", en_requested, "\n")
        end

    elseif    sa == "multiplet: energy of each level relative to lowest level"
        println(stream, "\n  Energy of each level relative to lowest level for nuclear spin I = $(multiplet.levelFs[1].I):")
        sb = "  Level  F Parity          Hartrees       " * "             eV                   " * JAC.TableStrings.inUnits("energy")
        println(stream, "\n", sb, "\n")
        for  i = 2:length(multiplet.levelFs)
            lev = multiplet.levelFs[i]
            en    = lev.energy - multiplet.levelFs[1].energy;
            en_eV = Defaults.convertUnits("energy: from atomic to eV", en);    en_requested = Defaults.convertUnits("energy: from atomic", en)
            sc  = " " * JAC.TableStrings.level(i) * "    " * string(LevelSymmetry(lev.F, lev.parity))  * "    "
            @printf(stream, "%s %.15e %s %.15e %s %.15e %s", sc, en, "  ", en_eV, "  ", en_requested, "\n")
        end
    else
        error("Unsupported keystring.")
    end

    return( nothing )
end

"""
`Basics.tabulateKappaSymmetryEnergiesDirac(kappa::Int64, evalues::Array{Float64,1}, ns::Int64, nuclearModel::Nuclear.Model)`
    ... tabulates the eigenenergies for a given symmetry block kappa together with the corresponding Dirac energies for a
        point-like nucleus. The index ns tells the number of 'negative-continnum' energies in the given evalues. nothing is
        returned.
"""
function Basics.tabulateKappaSymmetryEnergiesDirac(kappa::Int64, evalues::Array{Float64,1}, ns::Int64, nuclearModel::Nuclear.Model)
    Z = nuclearModel.Z
    # Determine the allowed principal quantum numbers n
    l = Basics.subshell_l(Subshell(101,kappa))
    println("  ", JAC.TableStrings.hLine(77))
    sa = "  "
    sa = sa * JAC.TableStrings.center( 7, "Index";            na=2)
    sa = sa * JAC.TableStrings.center(10, "Subshell";         na=3)
    sa = sa * JAC.TableStrings.center(17, "Energies [a.u.]";  na=2)
    sa = sa * JAC.TableStrings.center(17, "Dirac-E  [a.u.]";  na=2)
    sa = sa * JAC.TableStrings.center(17, "Delta-E / |E|";    na=2)
    println(sa)
    println("  ", JAC.TableStrings.hLine(77))
    for  i = ns+1:ns+8
        sa = " " * JAC.TableStrings.center( 6, JAC.TableStrings.level(i-ns); na=2)
        sa = sa *  JAC.TableStrings.flushright(10, string(Subshell(i-ns+l, kappa)); na=6)
        en = Basics.computeDiracEnergy(Subshell(i-ns+l, kappa), Z)
        if  evalues[i] >= 0.      sa = sa * "+"   end;    sa = sa * @sprintf("%.8e", evalues[i])                       * "    "
        if  en         >= 0.      sa = sa * "+"   end;    sa = sa * @sprintf("%.8e", en)                               * "    "
        if  evalues[i]-en >= 0.   sa = sa * "+"   end;    sa = sa * @sprintf("%.8e", (evalues[i]-en)/abs(evalues[i])) * "    "
        println(sa)
    end
    println("      :       :    ")
    for  i = length(evalues)-5:length(evalues)
        sa = " " * JAC.TableStrings.center( 6, JAC.TableStrings.level(i-ns); na=2)
        sa = sa *  JAC.TableStrings.flushright(10, string(Subshell(i-ns+l, kappa)); na=6)
        en = Basics.computeDiracEnergy(Subshell(i-ns+l, kappa), Z)
        if  evalues[i] >= 0.      sa = sa * "+"   end;    sa = sa * @sprintf("%.8e", evalues[i])                       * "    "
        if  en         >= 0.      sa = sa * "+"   end;    sa = sa * @sprintf("%.8e", en)                               * "    "
        if  evalues[i]-en >= 0.   sa = sa * "+"   end;    sa = sa * @sprintf("%.8e", (evalues[i]-en)/abs(evalues[i])) * "    "
        println(sa)
    end
    println("  ", JAC.TableStrings.hLine(77))

    return( nothing )
end


text:
mutable struct Release <: GitHubType
    url::Union{Nothing, HTTP.URI}
    html_url::Union{Nothing, HTTP.URI}
    assets_url::Union{Nothing, HTTP.URI}
    upload_url::Union{Nothing, HTTP.URI}
    tarball_url::Union{Nothing, HTTP.URI}
    zipball_url::Union{Nothing, HTTP.URI}
    id::Union{Nothing, Int}
    node_id::Union{Nothing, String}
    tag_name::Union{Nothing, String}
    target_commitish::Union{Nothing, String}
    name::Union{Nothing, String}
    body::Union{Nothing, String}
    draft::Union{Nothing, Bool}
    prerelease::Union{Nothing, Bool}
    created_at::Union{Nothing, String}
    published_at::Union{Nothing, String}
    author::Union{Nothing, Dict{String, Any}}
    assets::Union{Nothing, Array{Any, 1}}
end

Release(data::Dict) = json2github(Release, data)
namefield(r::Release) = r.id

@api_default function create_release(api::GitHubAPI, repo; options...)
    result = gh_post_json(api, "/repos/$(name(repo))/releases"; options...)
    return Release(result)
end

@api_default function releases(api::GitHubAPI, repo; options...)
    result, paged_data = gh_get_paged_json(api, "/repos/$(name(repo))/releases"; options...)
    return map(Release, result), paged_data
end


text:
"""
    source_paths(mod, file)

Returns list of all source files in the module with matching (partial) filename,
the filename does not have to be relative to the base directory.
It can be any partial path.
if no matching path is found then an empty list is returned
"""
function source_paths(mod, file)
    mdata = pkgfiles(mod)
    if mdata === nothing
        is_stdlib = CodeTracking.PkgId(mod).uuid === nothing
        is_stdlib && Revise.track(mod) # need to track stdlibs explictly

        Revise.revise()   # Or maybe a Revise was pending
        mdata = pkgfiles(mod)
        if mdata === nothing
            error(
                "According to CodeTracking.jl $mod is not loaded. " *
                "Revise.jl might be malfuncting."
            )
        end
    end
    mfiles = joinpath.(mdata.basedir, mdata.files)

    file = expanduser(file)
    if isfile(file)
        # This was something point to a real file that exists at a location
        # not just filename **somewhere**
        # It may have been a relative path, if so that is hard to deal with
        # So we will make it absolute.
        file=abspath(file)
    end

    if !isabspath(file)
        # We do not need it to be absolute
        # But we do want it to start with a "/" if it isn't
        # So we can use `endswith` to check that it matchs one of the paths
        # that we know about, without worry of matching part of a filename
        file = "/" * file
    end
    matched_inds = map(mfile->endswith(mfile, file), mfiles)
    return mfiles[matched_inds]
end

"""
    containing_methods([module], filename, linenum)

Returns the methods within which the provided line number
in the given file, occurs.
Returns an empty list on no match.

If the module is not provided, then all modules loaded
will be searched for a file with that name that has a function over that line.

Filenames can be absolute, or partial.
E.g. `/home/user1/dev/MyMod/src/helpers/utils.jl`,
`helpers/utils.jl` or `utils.jl` are all acceptable.

If it is ambigious, then all matching methods in all file will be returned.
"""
function containing_methods(mod, file, linenum)
    paths = source_paths(mod, file)
    isempty(paths) && return Method[]
    sigs = mapreduce(vcat, paths) do path
        signatures_at(path, linenum)
    end
    # TODO: workout  a way to round-forward as the signatures_at
    # only matches to
    # statement within the functions body, not from the line it is "declared"
    # or any intervening whitespace. It also doesn't match to the `end` line.

    ret = map(sigt2meth, sigs)
    return ret
end

function containing_methods(file, linenum)
    #TODO: raise issue on CodeTracking.jl to expose public way
    # to get list of all loaded modules
    pkgids = keys(CodeTracking._pkgfiles)
    isempty(pkgids) && return Method[]
    mapreduce(vcat, pkgids) do pkgid
        containing_methods(pkgid, file, linenum)
    end
end

function sigt2meth(::Type{SIGT})::Method where SIGT
    params = parameter_typeof(SIGT)
    func_t = params[1]
    func = func_t.instance

    args_t = Tuple{params[2:end]...}

    meths = methods(func, args_t)
    if length(meths) > 1
        @info "Multiple possible methods. Falling back to first" func args_t meths
    end
    return first(meths)
end

##############
"""
    src_line2ir_statement_ind(method|ir, src_line)

Given a CodeIR, and line number from source code
determines the index of the last IR statement that occurs on that line.
"""
function src_line2ir_statement_ind(ir, src_line)
    linetable_ind = findlast(ir.linetable) do lineinfo
        lineinfo.line == src_line
    end
    statement_ind = findlast(isequal(linetable_ind), ir.codelocs)
    return statement_ind
end

function src_line2ir_statement_ind(meth::Method, src_line)
    ir = Base.uncompressed_ast(meth)
    return src_line2ir_statement_ind(ir, src_line)
end

"""
    ir_statement_ind2src_linenum(method|ir, statement_ind)

Given a CodeIR, and a statement index (code loc) within that method's IR
determines the line in the source code wher that occurs.
"""
function statement_ind2src_linenum(ir, statement_ind)
    # handle boundry cases
    statement_ind < 0 && throw(ArgumentError("Can not have negative statement index"))
    # We use 0 for "before first statement".
    statement_ind == 0 && return before_body_linenum(ir)

    code_loc = ir.codelocs[statement_ind]
    return ir.linetable[code_loc].line
end

function statement_ind2src_linenum(meth::Method, statement_ind)
    ir = Base.uncompressed_ast(meth)
    return statement_ind2src_linenum(ir, statement_ind)
end

"""
    before_body_linenum(ir)
Returns the line number of the line before the first line in the body of the ir.
"""
function before_body_linenum(ir)
    first_lineinfo = ir.linetable[first(ir.codelocs)]
    return first_lineinfo.line - 1
end

######################################################

"""
    loc_for_file(file)
Returns a vector of lines of code for the file.
This special cases the name `"REPL[\\d]"`, as being REPL history
and so returns that REPL history cell instead.
If for some reason the file can not be found,
this returns a vector with a single line and a message explaining as such.
"""
function loc_for_file(file::AbstractString)
    if isfile(file)
        return readlines(file)
    elseif startswith(file, "REPL[") && isdefined(Base, :active_repl)
        # extract the number from "REPL[123]"
        hist_idx = parse(Int,string(file)[6:end-1])
        hist = Base.active_repl.interface.modes[1].hist
        source_code = hist.history[hist.start_idx+hist_idx]
        return split(source_code, "\n")
    else
        return ["-- source for $file not found --"]
    end
end


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

"""
    ValueShapes.default_datatype(T::Type)

Return a default specific type U that is more specific than T, with U <: T.

e.g.

    ValueShapes.default_datatype(Real) == Float64
    ValueShapes.default_datatype(Complex) == Complex{Float64}
"""
function default_datatype end

@inline default_datatype(::Type{>:Int}) = Int
@inline default_datatype(::Type{>:Float64}) = Float64
@inline default_datatype(::Type{>:Real}) = Float64
@inline default_datatype(::Type{>:Complex{Float64}}) = Complex{Float64}
@inline default_datatype(T::Type) = T

"""
    abstract type AbstractValueShape

An `AbstractValueShape` combines type and size information.

Subtypes are defined for shapes of scalars (see [`ScalarShape`](@ref)),
arrays (see [`ArrayShape`](@ref)), constant values
(see [`ConstValueShape`](@ref)) and `NamedTuple`s (see
[`NamedTupleShape`](@ref)).

Subtypes of `AbstractValueShape` must support `eltype`, `size` and
[`totalndof`](@ref).

Value shapes can be used as constructors to generate values of the given
shape with undefined content. If the element type of the shape is an abstract
or union type, a suitable concrete type will be chosen automatically, if
possible (see [`ValueShapes.default_datatype`](@ref)):

```julia
shape = ArrayShape{Real}(2,3)
A = shape(undef)
typeof(A) == Array{Float64,2}
size(A) == (2, 3)
valshape(A) == ArrayShape{Float64}(2,3)
```

Use

    (shape::AbstractValueShape)(data::AbstractVector{<:Real})::eltype(shape)

to view a flat vector of anonymous real values
as a value of the given shape:

```julia
data = [1, 2, 3, 4, 5, 6]
shape(data) == [1 3 5; 2 4 6]
```

In return,

    Base.Vector{<:Real}(undef, shape::AbstractValueShape)

will create a suitable uninitialized vector of the right length to hold such
flat data for the given shape. If no type `T` is given, a suitable data
type will be chosen automatically.

When dealing with multiple vectors of flattened data, use

    shape.(data::ArraysOfArrays.AbstractVectorOfSimilarVectors)

ValueShapes supports this via specialized broadcasting.

In return,

    ArraysOfArrays.VectorOfSimilarVectors{<:Real}(shape::AbstractValueShape)

will create a suitable vector (of length zero) of vectors that can hold
flattened data for the given shape. The result will be a
`VectorOfSimilarVectors` wrapped around a 2-dimensional `ElasticArray`.
This way, all data is stored in a single contiguous chunk of memory.

`AbstractValueShape`s can be compared with `<=` and `>=`, with semantics that
are similar to compare type with `<:` and `>:`:

```julia
a::AbstractValueShape <= b::AbstractValueShape == true
```

implies that values of shape `a` are can be used in contexts that expect
values of shape `b`. E.g.:

```julia
(ArrayShape{Float64}(4,5) <= ArrayShape{Real}(4,5)) == true
(ArrayShape{Float64}(4,5) <= ArrayShape{Integer}(4,5)) == false
(ArrayShape{Float64}(2,2) <= ArrayShape{Float64}(3,3)) == false
(ScalarShape{Real}() >= ScalarShape{Int}()) == true
```
"""
abstract type AbstractValueShape end

export AbstractValueShape

import Base.>=
@inline >=(a::AbstractValueShape, b::AbstractValueShape) = b <= a

# Reserve broadcasting semantics for value shapes:
@inline Base.Broadcast.broadcastable(shape::AbstractValueShape) =
    throw(ArgumentError("broadcasting over `AbstractValueShape`s is reserved"))

function _valshapeoftype end

"""
    ValueShapes.default_unshaped_eltype(shape::AbstractValueShape)

Returns the default real array element type to use for unshaped
representations of data with shape `shape`.

Subtypes of `AbstractValueShape` must implemenent
`ValueShapes.default_unshaped_eltype`.
"""
function default_unshaped_eltype end

"""
    ValueShapes.shaped_type(shape::AbstractValueShape, ::Type{T}) where {T<:Real}
    ValueShapes.shaped_type(shape::AbstractValueShape)

Returns the type the will result from reshaping a real-valued vector (of
element type `T`, if specified) with `shape`.

Subtypes of `AbstractValueShape` must implement

    ValueShapes.shaped_type(shape::AbstractValueShape, ::Type{T}) where {T<:Real}
"""
function shaped_type end

shaped_type(shape::AbstractValueShape) = shaped_type(shape, default_unshaped_eltype(shape))

"""
    valshape(x)::AbstractValueShape
    valshape(acc::ValueAccessor)::AbstractValueShape

Get the value shape of an arbitrary value, resp. the shape a `ValueAccessor`
is based on, or the shape of the variates for a `Distribution`.
"""
function valshape end
export valshape

@inline valshape(x::T) where T = _valshapeoftype(T)

"""
    elshape(x)::AbstractValueShape

Get the shape of the elements of x
"""
function elshape end
export elshape

@inline elshape(x::T) where T = _valshapeoftype(eltype(T))

@inline elshape(A::AbstractArray{<:AbstractArray}) = ArrayShape{eltype(eltype(A))}(innersize(A)...)

"""
    totalndof(shape::AbstractValueShape)

Get the total number of degrees of freedom of values of the given shape.

Equivalent to the length of a vector that would result from flattening the
data into a sequence of real numbers, excluding any constant values.
"""
function totalndof end
export totalndof

"""
    unshaped(x)::AbstractVector{<:Real}
    unshaped(x, shape::AbstractValueShape)::AbstractVector{<:Real}

Retrieve the unshaped underlying data of x, assuming x is a structured view
(based on some [`AbstractValueShape`](@ref)) of a flat/unstructured
real-valued data vector.

If `shape` is given, ensures that the shape of `x` is compatible with it.
Specifying a shape may be necessary if the correct shape of `x` cannot be
inferred from `x`, e.g. because `x` is assumed to have fewer degrees of
freedom (because of constant components) than would be inferred from
the plain value of `x`.

Example:

```julia
shape = NamedTupleShape(
    a = ScalarShape{Real}(),
    b = ArrayShape{Real}(2, 3)
)
data = [1, 2, 3, 4, 5, 6, 7]
x = shape(data)
@assert unshaped(x) == data
@assert unshaped(x.a) == view(data, 1:1)
@assert unshaped(x.b) == view(data, 2:7)
```
"""
function unshaped end
export unshaped

unshaped(x::Real) = Fill(x, 1)
unshaped(x::AbstractArray{<:Real,0}) = view(x, firstindex(x):firstindex(x))
unshaped(x::SubArray{<:Real,0}) = view(parent(x), x.indices[1]:x.indices[1])
unshaped(x::AbstractArray{<:Real,1}) = x
unshaped(x::Base.ReshapedArray{T,N,<:AbstractArray{T,1}}) where {T<:Real,N} = parent(x)

"""
    stripscalar(x)

Dereference value `x`.

If x is a scalar-like object, like a 0-dimensional array or a `Ref`,
`stripscalar` returns it's inner value. Otherwise, `x` is returned unchanged.

Useful to strip shaped scalar-like views of their 0-dim array semantics
(if present), but leave array-like views unchanged.

Example:

```julia
data = [1, 2, 3]
shape1 = NamedTupleShape(a = ScalarShape{Real}(), b = ArrayShape{Real}(2))
x1 = shape1(data)
@assert x1 isa AbstractArray{<:NamedTuple,0}
@assert stripscalar(x) isa NamedTuple

shape2 = ArrayShape{Real}(3)
x2 = shape2(data)
@assert x2 isa AbstractArray{Int,1}
@assert ref(x2) isa AbstractArray{Int,1}
```
"""
function stripscalar end
export stripscalar

stripscalar(x::Any) = x
stripscalar(x::Ref) = x[]
stripscalar(x::AbstractArray{T,0}) where T = x[]

function _checkcompat(shape::AbstractValueShape, data::AbstractVector{<:Real})
    n_shape = totalndof(shape)
    n_data = length(eachindex(data))
    if n_shape != n_data
        throw(ArgumentError("Data vector of length $(n_data) incompatible with value shape with $(n_shape) degrees of freedom"))
    end
    nothing
end

function _checkcompat_inner(shape::AbstractValueShape, data::AbstractArray{<:AbstractVector{<:Real}})
    n_shape = totalndof(shape)
    n_data = prod(innersize(data))
    if n_shape != n_data
        throw(ArgumentError("Data vector of length $(n_data) incompatible with value shape with $(n_shape) degrees of freedom"))
    end
    nothing
end

@inline function _apply_shape_to_data(shape::AbstractValueShape, data::AbstractVector{<:Real})
    @boundscheck _checkcompat(shape, data)
    view(data, ValueAccessor(shape, 0))
end

@inline function (shape::AbstractValueShape)(data::AbstractVector{<:Real})
    _apply_shape_to_data(shape, data)
end

Base.Vector{T}(::UndefInitializer, shape::AbstractValueShape) where {T <: Real} =
    Vector{T}(undef, totalndof(shape))

Base.Vector{<:Real}(::UndefInitializer, shape::AbstractValueShape) =
    Vector{default_unshaped_eltype(shape)}(undef, shape)

ArraysOfArrays.VectorOfSimilarVectors{T}(shape::AbstractValueShape) where {T<:Real} =
    VectorOfSimilarVectors(ElasticArray{T}(undef, totalndof(shape), 0))

const VSBroadcasted1{N,F,T} = Base.Broadcast.Broadcasted{
    <:Base.Broadcast.AbstractArrayStyle{N},
    <:Any,
    F,
    <:Tuple{T}
}

const VSBroadcasted2{N,F,T1,T2} = Base.Broadcast.Broadcasted{
    <:Base.Broadcast.AbstractArrayStyle{N},
    <:Any,
    F,
    <:Tuple{T1,T2}
}

# Specialize (::AbstractValueShape).(::AbstractVector{<:AbstractVector{<:Real}}):
Base.copy(instance::VSBroadcasted1{N,<:AbstractValueShape,AbstractArray{<:AbstractVector{<:Real},N}}) where N =
    broadcast(view, instance.args[1], Ref(ValueAccessor(instance.f, 0)))

# Specialize unshaped for real vectors (semantically vectors of scalar-shaped values)
function Base.copy(instance::VSBroadcasted1{1,typeof(unshaped),AbstractVector{<:Real}})
    x = instance.args[1]
    nestedview(reshape(view(x, :), 1, length(eachindex(x))))
end

# Specialize unshaped for real vectors that are array slices:
const _MatrixSliceFirstDim{T} = SubArray{T,1,<:AbstractArray{T,2},<:Tuple{Int,AbstractArray{Int}}}
function Base.copy(instance::VSBroadcasted1{1,typeof(unshaped),<:_MatrixSliceFirstDim{<:Real}})
    instance
    x = instance.args[1]
    nestedview(view(parent(x), x.indices[1]:x.indices[1], x.indices[2]))
end

function _zerodim_array(x::T) where T
    A = Array{T,0}(undef)
    A[] = x
end

"""
    const_zero(x::Any)

Get the equivalent of a constant zero for values the same type as .
"""
function const_zero end

const_zero(x::Number) = zero(x)
const_zero(A::AbstractArray{T}) where T <: Number = Fill(zero(T), size(A)...)

"""
    replace_const_shapes(f::Function, shape::AbstractValueShape)

If `shape` is a, or contains, [`ConstValueShape`](@ref) shape(s), recursively
replace it/them with the result of `f(s::Shape)`.
"""
function replace_const_shapes end
export replace_const_shapes

"""
    gradient_shape(argshape::AbstractValueShape)

Return the value shape of the gradient of functions that take values of
shape `argshape` as an input.
"""
function gradient_shape end

gradient_shape(vs::AbstractValueShape) = replace_const_shapes(const_zero_shape, vs)
export gradient_shape


text:
# # Arrays and Loops
#
# We'll cover:
# 1. Array literals
# 2. Concatenation
# 3. For loops
# 4. Comprehensions
# 5. Element types
# 6. Dequeues
# 7. The `bang!` convention
# 8. Variable names vs. copies
#
# # Arrays
#
# Julia has highly efficient multidimensional arrays, both constructed and indexed with square brackets.
#
# Syntax: <br>
# ```julia
# [item1, item2, ...]
# ```

squares = [1, 4, 9, 15, 25, 36, 49, 64]

#-

squares[1]

#-

squares[1:3]

#-

squares[end]

#-

squares[4] = 16

#-

squares

#-

typeof(squares)

# ## Concatenation
#
# If, instead of commas, you just use spaces, then the values are concatenated horizontally.

cubes = [1, 8, 27, 64, 125, 216, 343, 512]

#-

powers = [1:8 squares cubes]

#-

powers[4, 2]

#-

powers[:, 3]

#-

powers[7, :]

#-

typeof(powers)

# Semicolon separators perform vertical concatenation:

[squares; cubes]

# Whereas commas would simply create an array of arrays:

nested_powers = [[1,2,3,4,5,6,7,8], squares, cubes]

#-

nested_powers[2]

# Horizontal and vertical concatenation can be used together to as a simple syntax for matrix literals:

[1 3 5; 2 4 6]

# # Loops

#-

# Of course, we could construct this programmatically with a for-loop. The syntax for a `for` loop is
#
# ```julia
# for *var* in *loop iterable*
#     *loop body*
# end
# ```
#

A = fill(0, (8, 3)) # Allocate an 8x3 matrix to store the values into
for pow in 1:3
    for value in 1:8
        A[value, pow] = value ^ pow
    end
end
A

#-

A == powers

# ## Array Comprehensions

squares = [value^2 for value in 1:8]

#-

cubes = [value^3 for value in 1:8]

#-

powers = [value^pow for value in 1:8, pow in 1:3]

# # The element type
#
# Note that every time an array prints out, it is displaying its element type and dimensionality, for example `Array{Int64, 2}`. This describes what it can store — and thus what it can return upon indexing.

typeof(powers)

#-

typeof(powers[1, 1])

# Further, the array will try to convert any new values assigned into it to its element type:

powers[1, 1] = 1.6

#-

powers[1, 1] = -5.0 # This can be losslessly converted to an integer

#-

powers

# Arrays that have an exact and concrete element type are generally significantly faster, so Julia will try to find an amenable element type for you in its literal construction syntax:

fortytwosarray = [42, 42.0, 4.20e1, 4.20f1, 84//2, 0x2a]

#-

for x in fortytwosarray
    show(x)
    println("\tisa $(typeof(x))")
end

# The `Any` array can be helpful for disabling these behaviors and allowing all kinds of different objects:

anyfortytwos = Any[42, 42.0, 4.20e1, 4.20f1, 84//2, 0x2a]

#-

anyfortytwos[1] = "FORTY TWO"
anyfortytwos

# # Vectors as dequeues

#-

# One-dimensional arrays can be appended to and have items removed from them:

fib = [1, 1, 2, 3, 5, 8, 13]

#-

push!(fib, 21)

#-

push!(fib, sum(fib[end-1:end]))

#-

pop!(fib)

#-

pushfirst!(fib, 0)

#-

popfirst!(fib)

# ## Aside: why so shouty!?
#
# Why are there exclamations in all the above function calls? `push!(fib, ...)`, `pop!(fib)`, etc.?

push!

# This is entirely a convention. It's a signal to the caller that one of the argument is going to be _mutated_. It's perhaps easiest to demonstrate with an example:

A = rand(0:10, 10)

#-

sort(A)

#-

A

#-

sort!(A)

#-

A # That changed A!

# ## Aside: names vs. copies
#
# Remember that variables are just names we give our objects. So watch what happens if we give the same object two different names:

fibonacci = [1, 1, 2, 3, 5, 8, 13]

#-

some_numbers = fibonacci
some_numbers[1] = 404
some_numbers

#-

fibonacci

#-

fibonacci[1] = 1
some_numbers = copy(fibonacci)
some_numbers[2] = 404
fibonacci

max_stars_repo_path: Courses/Introduction to Julia/02. Arrays and Loops.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
"""
    insert_location(ex::Expr, location)

Insert a symbolic representation of `location` into the arguments of an `expression`.

Used in the `@at` macro for specifying the location of an `AbstractOperation`.
"""
function insert_location!(ex::Expr, location)
    if ex.head === :call && ex.args[1] ∈ operators
        push!(ex.args, ex.args[end])
        ex.args[3:end-1] .= ex.args[2:end-2]
        ex.args[2] = location
    end

    for arg in ex.args
        insert_location!(arg, location)
    end

    return nothing
end

"Fallback for when `insert_location` is called on objects other than expressions."
insert_location!(anything, location) = nothing

"""
    @at location abstract_operation

Modify the `abstract_operation` so that it returns values at
`location`, where `location` is a 3-tuple of `Face`s and `Cell`s.
"""
macro at(location, abstract_operation)
    insert_location!(abstract_operation, location)
    return esc(abstract_operation)
end


text:
function Ellipse(arg0::Vector3D, arg1::Vector3D, arg2::Vector3D, arg3::jdouble, arg4::jdouble, arg5::Frame)
    return Ellipse((Vector3D, Vector3D, Vector3D, jdouble, jdouble, Frame), arg0, arg1, arg2, arg3, arg4, arg5)
end

function get_a(obj::Ellipse)
    return jcall(obj, "getA", jdouble, ())
end

function get_b(obj::Ellipse)
    return jcall(obj, "getB", jdouble, ())
end

function get_center(obj::Ellipse)
    return jcall(obj, "getCenter", Vector3D, ())
end

function get_center_of_curvature(obj::Ellipse, arg0::Vector2D)
    return jcall(obj, "getCenterOfCurvature", Vector2D, (Vector2D,), arg0)
end

function get_frame(obj::Ellipse)
    return jcall(obj, "getFrame", Frame, ())
end

function get_u(obj::Ellipse)
    return jcall(obj, "getU", Vector3D, ())
end

function get_v(obj::Ellipse)
    return jcall(obj, "getV", Vector3D, ())
end

function point_at(obj::Ellipse, arg0::jdouble)
    return jcall(obj, "pointAt", Vector3D, (jdouble,), arg0)
end

function project_to_ellipse(obj::Ellipse, arg0::TimeStampedPVCoordinates)
    return jcall(obj, "projectToEllipse", TimeStampedPVCoordinates, (TimeStampedPVCoordinates,), arg0)
end

function project_to_ellipse(obj::Ellipse, arg0::Vector2D)
    return jcall(obj, "projectToEllipse", Vector2D, (Vector2D,), arg0)
end

function to_plane(obj::Ellipse, arg0::Vector3D)
    return jcall(obj, "toPlane", Vector2D, (Vector3D,), arg0)
end

function to_space(obj::Ellipse, arg0::Vector2D)
    return jcall(obj, "toSpace", Vector3D, (Vector2D,), arg0)
end


text:
# a rather literal translation of the pseudocode at https://en.wikipedia.org/wiki/MD5

const s = UInt32[7, 12, 17, 22, 7, 12, 17, 22, 7, 12, 17, 22, 7, 12, 17, 22,
                 5, 9, 14, 20, 5, 9, 14, 20, 5, 9, 14, 20, 5, 9, 14, 20,
                 4, 11, 16, 23, 4, 11, 16, 23, 4, 11, 16, 23,  4, 11, 16, 23,
                 6, 10, 15, 21, 6, 10, 15, 21, 6, 10, 15, 21, 6, 10, 15, 21]

const K = UInt32[0xd76aa478, 0xe8c7b756, 0x242070db, 0xc1bdceee, 0xf57c0faf,
    0x4787c62a, 0xa8304613, 0xfd469501, 0x698098d8, 0x8b44f7af, 0xffff5bb1,
    0x895cd7be, 0x6b901122, 0xfd987193, 0xa679438e, 0x49b40821, 0xf61e2562,
    0xc040b340, 0x265e5a51, 0xe9b6c7aa, 0xd62f105d, 0x02441453, 0xd8a1e681,
    0xe7d3fbc8, 0x21e1cde6, 0xc33707d6, 0xf4d50d87, 0x455a14ed, 0xa9e3e905,
    0xfcefa3f8, 0x676f02d9, 0x8d2a4c8a, 0xfffa3942, 0x8771f681, 0x6d9d6122,
    0xfde5380c, 0xa4beea44, 0x4bdecfa9, 0xf6bb4b60, 0xbebfbc70, 0x289b7ec6,
    0xeaa127fa, 0xd4ef3085, 0x04881d05, 0xd9d4d039, 0xe6db99e5, 0x1fa27cf8,
    0xc4ac5665, 0xf4292244, 0x432aff97, 0xab9423a7, 0xfc93a039, 0x655b59c3,
    0x8f0ccc92, 0xffeff47d, 0x85845dd1, 0x6fa87e4f, 0xfe2ce6e0, 0xa3014314,
    0x4e0811a1, 0xf7537e82, 0xbd3af235, 0x2ad7d2bb, 0xeb86d391]

function md5(msgbytes)
    a0::UInt32 = 0x67452301  # A
    b0::UInt32 = 0xefcdab89  # B
    c0::UInt32 = 0x98badcfe  # C
    d0::UInt32 = 0x10325476  # D

    oldlen = length(msgbytes)
    umsg = push!([UInt8(b) for b in msgbytes], UInt8(0x80))
    while length(umsg) % 64 != 56
        push!(umsg, UInt8(0))
    end
    append!(umsg, reinterpret(UInt8, [htol(UInt64(oldlen) * 8)]))

    for j in 1:64:length(umsg)-1
        arr = view(umsg, j:j+63)
        M = [reinterpret(UInt32, arr[k:k+3])[1] for k in 1:4:62]
        A = a0
        B = b0
        C = c0
        D = d0

        for i in 0:63
            if 0 ≤ i ≤ 15
                F = D ⊻ (B & (C ⊻ D))
                g = i
            elseif 16 ≤ i ≤ 31
                F = C ⊻ (D & (B ⊻ C))
                g = (5 * i + 1) % 16
            elseif 32 ≤ i ≤ 47
                F = B ⊻ C ⊻ D
                g = (3 * i + 5) % 16
            elseif 48 ≤ i ≤ 63
                F = C ⊻ (B | (~D))
                g = (7 * i) % 16
            end
            F += A + K[i+1] + M[g+1]
            A = D
            D = C
            C = B
            B += ((F) << s[i+1]) | (F >> (32 - s[i+1]))
        end

        a0 += A
        b0 += B
        c0 += C
        d0 += D
    end
    digest = join(map(x -> lpad(string(x, base=16), 2, '0'), reinterpret(UInt8, [a0, b0, c0, d0])), "") # Output is in little-endian
end

for pair in [0xd41d8cd98f00b204e9800998ecf8427e => "", 0x0cc175b9c0f1b6a831c399e269772661 => "a",
   0x900150983cd24fb0d6963f7d28e17f72 => "abc", 0xf96b697d7cb7938d525a2f31aaf161d0 => "message digest",
   0xc3fcd3d76192e4007dfb496cca67e13b => "abcdefghijklmnopqrstuvwxyz",
   0xd174ab98d277d9f5a5611c2c9f419d9f => "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789",
   0x57edf4a22be3c955ac49da2e2107b67a => "12345678901234567890123456789012345678901234567890123456789012345678901234567890"]
   println("MD5 of $(pair[2]) is $(md5(pair[2])), which checks with $(string(pair[1], base=16)).")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
struct BagModel{T <: MillModel, U <: MillModel} <: MillModel
im::T
a::Aggregation
bm::U
end

use a `im` model on data in `BagNode`, the uses function `a` to aggregate individual bags,
and finally it uses `bm` model on the output
"""
struct BagModel{T <: MillModel, A, U <: MillModel} <: MillModel
    im::T
    a::A
    bm::U
end

Flux.@treelike BagModel
# Flux.@functor BagModel

BagModel(im::MillFunction, a, bm::MillFunction) = BagModel(ArrayModel(im), a, ArrayModel(bm))
BagModel(im::MillModel, a, bm::MillFunction) = BagModel(im, a, ArrayModel(bm))
BagModel(im::MillFunction, a, bm::MillModel) = BagModel(ArrayModel(im), a, bm)
BagModel(im::MillFunction, a) = BagModel(im, a, identity)
BagModel(im::MillModel, a) = BagModel(im, a, ArrayModel(identity))

(m::BagModel)(x::WeightedBagNode{<: AbstractNode}) = m.bm(m.a(m.im(x.data), x.bags, x.weights))
# if the data is missing, we do not use the mapping
(m::BagModel)(x::WeightedBagNode{<: Missing}) = m.bm(ArrayNode(m.a(x.data, x.bags, x.weights)))

(m::BagModel)(x::BagNode{<: AbstractNode}) = m.bm(m.a(m.im(x.data), x.bags))
# if the data is missing, we do not use the mapping
(m::BagModel)(x::BagNode{<: Missing}) = m.bm(ArrayNode(m.a(x.data, x.bags)))

function modelprint(io::IO, m::BagModel{T, A}; pad=[], s="", tr=false) where {T, A}
    c = COLORS[(length(pad)%length(COLORS))+1]
    paddedprint(io, "BagModel$(tr_repr(s, tr))\n", color=c)
    paddedprint(io, "  ├── ", color=c, pad=pad)
    modelprint(io, m.im, pad=[pad; (c, "  │   ")], s=s * encode(1, 1), tr=tr)
    println(io);
    paddedprint(io, "  ├── ", color=c, pad=pad)
    modelprint(io, m.a, pad=[pad; (c, "  │   ")])
    # A == Aggregation || paddedprint(io, '\n', color=c)
    println(io);
    paddedprint(io, "  └── ", color=c, pad=pad)
    modelprint(io, m.bm, pad=[pad; (c, "  │   ")])
end


text:
InputFolder = './Images/Mitochondria/';
OutputFolder = './Results/Images/Mitochondria/';

@iA = '*.tif';

@fxg_mRadStdThr [iA] > [M];
params.GRad = 1;
params.RSet = [2 4 6 8];
params.Fracstd = 0.55;
params.SegThr = 0.55;
/endf

/show iA > M;
/keep M > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
# This file was generated, do not modify it. # hide
dtr = @load DecisionTreeRegressor
forest = EnsembleModel(atom=dtr)

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
using Pkg
Pkg.activate(".")

verbose = true
using Revise
 if verbose   println("# Loading RvSpecML")    end
 using RvSpectML
 if verbose   println("# Loading other packages")    end
 using DataFrames, Query, Statistics

# USER:  The default paths that specify where datafiles can be entered here or overridden in examples/data_paths.jl
target_subdir = "101501"   # USER: Replace with directory of your choice
 fits_target_str = "101501"
 output_dir = "examples/output"
 default_paths_to_search = [pwd(),"examples",joinpath(pkgdir(RvSpectML),"examples"),"/gpfs/group/ebf11/default/ebf11/expres/inputs"]

 make_plots = false
 write_ccf_to_csv = false
 write_rvs_to_csv = false
 write_template_to_csv = false
 write_spectral_grid_to_jld2 = false
 write_dcpca_to_csv = false

global has_loaded_data = false
global has_computed_ccfs = false
global has_computed_rvs = false
global has_computed_tempalte = false
global has_computed_dcpca = false

using Profile
for i in 1:2
global has_loaded_data
global has_computed_ccfs
global has_computed_rvs
global has_computed_tempalte
global has_computed_dcpca

   Profile.clear_malloc_data()
if !has_loaded_data
  global df_files = make_manifest(target_subdir, EXPRES )
  eval(code_to_include_param_jl())

  if verbose println("# Reading in FITS files.")  end
  @time global expres_data = map(EXPRES.read_data,eachrow(df_files_use))
  global has_loaded_data = true
end

if !has_computed_ccfs
 if verbose println("# Extracting order list timeseries from spectra.")  end
   order_list_timeseries = RvSpectML.make_order_list_timeseries(expres_data)
   order_list_timeseries = RvSpectML.filter_bad_chunks(order_list_timeseries,verbose=true)
   RvSpectML.normalize_spectra!(order_list_timeseries,expres_data);
 end

 if verbose println("# Reading line list for CCF: ", linelist_for_ccf_filename, ".")  end
 lambda_range_with_good_data = get_λ_range(expres_data)
 espresso_filename = joinpath(pkgdir(RvSpectML),"data","masks",linelist_for_ccf_filename)
 espresso_df = RvSpectML.read_linelist_espresso(espresso_filename)
 line_list_df = EXPRES.filter_line_list(espresso_df,first(expres_data).inst)

    # Compute CCF's & measure RVs
 if verbose println("# Computing CCF.")  end
 mask_shape = RvSpectML.CCF.TopHatCCFMask(order_list_timeseries.inst, scale_factor=tophap_ccf_mask_scale_factor)
 line_list = RvSpectML.CCF.BasicLineList(line_list_df.lambda, line_list_df.weight)
 ccf_plan = RvSpectML.CCF.BasicCCFPlan(mask_shape = mask_shape, line_list=line_list, midpoint=ccf_mid_velocity)
 v_grid = RvSpectML.CCF.calc_ccf_v_grid(ccf_plan)
 @time ccfs = RvSpectML.CCF.calc_ccf_chunklist_timeseries(order_list_timeseries, ccf_plan)
 # Write CCFs to file
 if write_ccf_to_csv
    using CSV
    CSV.write(target_subdir * "_ccfs.csv",Tables.table(ccfs',header=Symbol.(v_grid)))
    global has_computed_ccfs = true
 end

if make_plots
   using Plots
   t_idx = 1
   using Plots
   plot(v_grid,ccfs[:,t_idx],label=:none)
   xlabel!("v (m/s)")
   ylabel!("CCF")
 end

if !has_computed_rvs
 if verbose println("# Measuring RVs from CCF.")  end
 #rvs_ccf_gauss = [ RvSpectML.RVFromCCF.measure_rv_from_ccf(v_grid,ccfs[:,i],fit_type = "gaussian") for i in 1:length(order_list_timeseries) ]
 rvs_ccf_gauss = RvSpectML.RVFromCCF.measure_rv_from_ccf(v_grid,ccfs,fit_type = "gaussian")
 # Store estimated RVs in metadata
 map(i->order_list_timeseries.metadata[i][:rv_est] = rvs_ccf_gauss[i]-mean(rvs_ccf_gauss), 1:length(order_list_timeseries) )
 if write_rvs_to_csv
    using CSV
    CSV.write(target_subdir * "_rvs_ccf.csv",DataFrame("Time [MJD]"=>order_list_timeseries.times,"CCF RV [m/s]"=>rvs_ccf_gauss))
 end
 global has_computed_rvs = true
end

if !has_computed_tempalte
   if verbose println("# Making template spectra.")  end
   @time ( spectral_orders_matrix, f_mean, var_mean, deriv, deriv2, order_grids )  = RvSpectML.make_template_spectra(order_list_timeseries)

   if write_template_to_csv
      using CSV
      CSV.write(target_subdir * "_template.csv",DataFrame("λ"=>spectral_orders_matrix.λ,"flux_template"=>f_mean,"var"=>var_mean, "dfluxdlnλ_template"=>deriv,"d²fluxdlnλ²_template"=>deriv2))
   end
   if write_spectral_grid_to_jld2
      using JLD2, FileIO
      save(target_subdir * "_matrix.jld2", Dict("λ"=>spectral_orders_matrix.λ,"spectra"=>spectral_orders_matrix.flux,"var_spectra"=>spectral_orders_matrix.var,"flux_template"=>f_mean,"var"=>var_mean, "dfluxdlnλ_template"=>deriv,"d²fluxdlnλ²_template"=>deriv2))
   end
end

if !has_computed_dcpca
   if verbose println("# Performing Doppler constrained PCA analysis.")  end
   using MultivariateStats
   dcpca_out, M = RvSpectML.DCPCA.doppler_constrained_pca(spectral_orders_matrix.flux, deriv, rvs_ccf_gauss)
   frac_var_explained = 1.0.-cumsum(principalvars(M))./tvar(M)
   println("# Fraction of variance explained = ", frac_var_explained[1:min(5,length(frac_var_explained))])
   if write_dcpca_to_csv
      using CSV
      CSV.write(target_subdir * "_dcpca_basis.csv",  Tables.table(M.proj) )
      CSV.write(target_subdir * "_dcpca_scores.csv", Tables.table(dcpca_out) )
   end
end

#=
using MultivariateStats
function doppler_constrained_pca(flux::AbstractArray{T1,2}, deriv::AbstractVector{T2}, rvs::AbstractVector{T3} )  where { T1<:Real, T2<:Real, T3<:Real  }#, times::AbstractVector, rvs::AbstractVector)
  fm_perp = RvSpectML.compute_spectra_perp_doppler_shift(flux,deriv, rvs )
  idx_good_obs = 1:length(rvs)
  M = fit(PCA, fm_perp[:,idx_good_obs]; maxoutdim=10)
  pca_out = MultivariateStats.transform(M,fm_perp[:,idx_good_obs])
  return pca_out, M
end

=#
if make_plots
  using Plots
  # Set parameters for plotting analysis
  plt_order = 42
  plt_order_pix = 4500:5000
  plt = scatter(frac_var_explained, xlabel="Number of PCs", ylabel="Frac Variance Unexplained")
  display(plt)
end

if make_plots
  plt_order = 13
  RvSpectML.plot_basis_vectors(order_grids, f_mean, deriv, M.proj, idx_plt = spectral_orders_matrix.chunk_map[plt_order], num_basis=3 )
  #xlims!(5761.5,5766)
end

if make_plots
  RvSpectML.plot_basis_scores(order_list_timeseries.times, rvs_ccf_gauss, pca_out, num_basis=3 )
end

if make_plots
  RvSpectML.plot_basis_scores_cor( rvs_ccf_gauss, pca_out, num_basis=3)
end

global has_loaded_data = false
global has_computed_ccfs = false
global has_computed_rvs = false
global has_computed_tempalte = false
global has_computed_dcpca = false

end # for


text:
using Statistics: mean

# Code adapted from Flux.jl
# Ref: https://github.com/FluxML/Flux.jl/blob/master/src/layers/normalise.jl#L93-L177
# License: https://github.com/FluxML/Flux.jl/blob/master/LICENSE.md

istraining() = false

mutable struct InvertibleBatchNorm{T1,T2,T3} <: Bijector{1}
    b       ::  T1  # bias
    logs    ::  T1  # log-scale
    m       ::  T2  # moving mean
    v       ::  T2  # moving variance
    eps     ::  T3
    mtm     ::  T3  # momentum
end
function Base.:(==)(b1::InvertibleBatchNorm, b2::InvertibleBatchNorm)
    return  b1.b == b2.b &&
            b1.logs == b2.logs &&
            b1.m == b2.m &&
            b1.v == b2.v &&
            b1.eps == b2.eps &&
            b1.mtm == b2.mtm
end

function InvertibleBatchNorm(
    chs::Int;
    eps::T=1f-5,
    mtm::T=1f-1,
) where {T<:AbstractFloat}
    return InvertibleBatchNorm(
        zeros(T, chs),
        zeros(T, chs),  # logs = 0 means s = 1
        zeros(T, chs),
         ones(T, chs),
        eps,
        mtm,
    )
end

# define numerical parameters
# TODO: replace with `Functors.@functor InvertibleBatchNorm (b, logs)` when
# https://github.com/FluxML/Functors.jl/pull/7 is merged
function Functors.functor(::Type{<:InvertibleBatchNorm}, x)
    function reconstruct_invertiblebatchnorm(xs)
        return InvertibleBatchNorm(xs.b, xs.logs, x.m, x.v, x.eps, x.mtm)
    end
    return (b = x.b, logs = x.logs), reconstruct_invertiblebatchnorm
end

function forward(bn::InvertibleBatchNorm, x)
    dims = ndims(x)
    size(x, dims - 1) == length(bn.b) ||
        error("InvertibleBatchNorm expected $(length(bn.b)) channels, got $(size(x, dims - 1))")
    channels = size(x, dims - 1)
    as = ntuple(i -> i == ndims(x) - 1 ? size(x, i) : 1, dims)
    logs = reshape(bn.logs, as...)
    s = exp.(logs)
    b = reshape(bn.b, as...)
    if istraining()
        n = div(prod(size(x)), channels)
        axes = [1:dims-2; dims] # axes to reduce along (all but channels axis)
        m = mean(x, dims = axes)
        v = sum((x .- m) .^ 2, dims = axes) ./ n
        # Update moving mean and variance
        mtm = bn.mtm
        T = eltype(bn.m)
        bn.m = (1 - mtm) .* bn.m .+ mtm .* T.(reshape(m, :))
        bn.v = (1 - mtm) .* bn.v .+ (mtm * n / (n - 1)) .* T.(reshape(v, :))
    else
        m = reshape(bn.m, as...)
        v = reshape(bn.v, as...)
    end

    rv = s .* (x .- m) ./ sqrt.(v .+ bn.eps) .+ b
    logabsdetjac = (
        fill(sum(logs - log.(v .+ bn.eps) / 2), size(x, dims))
    )
    return (rv=rv, logabsdetjac=logabsdetjac)
end

logabsdetjac(bn::InvertibleBatchNorm, x) = forward(bn, x).logabsdetjac

(bn::InvertibleBatchNorm)(x) = forward(bn, x).rv

function forward(invbn::Inverse{<:InvertibleBatchNorm}, y)
    @assert !istraining() "`forward(::Inverse{InvertibleBatchNorm})` is only available in test mode."
    dims = ndims(y)
    as = ntuple(i -> i == ndims(y) - 1 ? size(y, i) : 1, dims)
    bn = inv(invbn)
    s = reshape(exp.(bn.logs), as...)
    b = reshape(bn.b, as...)
    m = reshape(bn.m, as...)
    v = reshape(bn.v, as...)

    x = (y .- b) ./ s .* sqrt.(v .+ bn.eps) .+ m
    return (rv=x, logabsdetjac=-logabsdetjac(bn, x))
end

(bn::Inverse{<:InvertibleBatchNorm})(y) = forward(bn, y).rv

function Base.show(io::IO, l::InvertibleBatchNorm)
    print(io, "InvertibleBatchNorm($(join(size(l.b), ", ")))")
end


text:
import Base.-, Base.+, Base.*, Base./
import Base.getindex, Base.setindex!

struct Field{T}
    values :: Array{T,1}
    mesh
    boundaries
end

struct BoundaryCondition{T}
    c1    :: Array{T,1}
    c2    :: Array{Scalar,1}
    field :: Field{T}
    patch :: Label
end

getindex(f::Field{T}, i::Int) where {T} = f.values[i]
setindex!(f::Field{T}, v::T, i::Int) where {T} = begin f.values[i] = v end

const ScalarField = Field{Scalar}
ScalarField(m::Mesh) = ScalarField( zeros(Scalar,length(m.volume)), m, Dict{String,Any}());

const VectorField = Field{Vector}
VectorField(m::Mesh) = VectorField( zeros(Vector,length(m.volume)), m, Dict{String,Any}());

for op ∈ [:+, :-, :*, :/]
    @eval begin
        function ($op)(a::Field{T}, b::Field{T}) where {T}
            assert(a.mesh == b.mesh)
            Field{T}(broadcast($op,a.values,b.values), a.mesh, Dict{String,Any}());
        end
    end
end

for op ∈ [:+, :-, :*, :/]
    @eval begin
        function ($op)(a::Field{T}, b::T) where {T}
            Field{T}(broadcast($op,a.values,b), a.mesh, Dict{String,Any}());
        end
    end
end

for op ∈ [:+, :-, :*, :/]
    @eval begin
        function ($op)(a::Field{T}, b::Array{T,1}) where {T}
            @assert size(a.values) == size(b)
            Field(broadcast($op,a.values,b), a.mesh, Dict{String,Any}());
        end
    end
end

*(a::Number, f::Field{T}) where {T} = Field(a*f.values, f.mesh, Dict{String,Any}());
-(a::Field{T}) where {T} = Field{T}(-a.values, a.mesh);

function ←(a::Field{T}, b::Field{T}) where {T}
    @assert a.mesh == b.mesh
    for i in eachindex(a.values)
        @inbounds a.values[i] = b.values[i]
    end
end

function ←(a::Field{T}, b::Float64) where {T}
    a.values[:] .= b
end

function ←(a::VectorField, b::Vector) where {T}
    for i in eachindex(a.values)
        @inbounds a.values[i] = b
    end
end

function ←(a::Field{T}, b::Array{T,1}) where {T}
    @assert length(a.values) == length(b)
    for i in eachindex(a.values)
        @inbounds a.values[i] = b[i]
    end
end

function boundary_coeffs(bc::BoundaryCondition{T}, f::Label) where {T}
    field = bc.field
    mesh = field.mesh
    start = first(mesh.patch[bc.patch].range)
    i = f + 1 - start

    c1 = bc.c1[i]
    c2 = bc.c2[i]

    return c1, c2
end

function boundary_value(bc::BoundaryCondition{T}, f::Label) where {T}
    field = bc.field
    mesh = field.mesh
    u1 = field.values[mesh.owner[f]]
    c1, c2 = boundary_coeffs(bc, f)
    return c1 + c2*u1
end

function set_dirichlet_patch!(field::Field{T}, name::String, val::T) where {T}
    mesh = field.mesh
    p = patch_by_name(mesh, name)
    c1 = [val for f in patch_faces(mesh, p)]
    c2 = zeros(Scalar, length(c1))
    field.boundaries[name] = BoundaryCondition(c1, c2, field, p)
end

function set_neumann_patch!(field::Field{T}, name::String) where {T}
    mesh = field.mesh
    p = patch_by_name(mesh, name)
    c2 = ones(length(patch_faces(mesh, p)))
    c1 = zeros(T, length(c2))
    field.boundaries[name] = BoundaryCondition(c1, c2, field, p)
end

function interpolate_to_points(f::Field{T}) where{T}
    mesh = f.mesh
    fn = zeros(T, length(mesh.point))
    vn = zeros(Scalar, length(mesh.point))
    for c in cells(mesh)
        for pid in mesh.cell2point[c]
            fn[pid] += f[c] * mesh.volume[c]
            vn[pid] += mesh.volume[c]
        end
    end
    return fn ./ vn
end

# Vykresleni pomoci matplotlib
function plot_contour(f::ScalarField, kwargs...)
    mesh = f.mesh
    x = [p[1] for p in mesh.point]
    y = [p[2] for p in mesh.point]
    tri  = triangles(mesh) .- 1
    fn = interpolate_to_points(f)

    PyPlot.tricontour(x, y, tri, fn; kwargs...)
end

function plot_contour(f::VectorField; kwargs...)
    mag = ScalarField(f.mesh)
    mag ← [norm(v) for v in f.values]
    plot_contour(mag; kwargs...)
end

function plot_contourf(f::ScalarField; kwargs...)
    mesh = f.mesh
    x = [p[1] for p in mesh.point]
    y = [p[2] for p in mesh.point]
    tri  = triangles(mesh)
    for i=1:length(tri)
        tri[i] .-= 1
    end
    fn = interpolate_to_points(f)

    PyPlot.tricontourf(x, y, tri, fn; kwargs...)
end

function plot_contourf(f::VectorField; kwargs...)
    mag = ScalarField(f.mesh)
    mag ← [norm(v) for v in f.values]
    plot_contourf(mag; kwargs...)
end

function plot_arrows(f::VectorField; kwargs...)
    mesh = f.mesh
    x = [c[1] for c in mesh.centre]
    y = [c[2] for c in mesh.centre]
    u = [v[1] for v in f.values]
    v = [v[2] for v in f.values]

    PyPlot.quiver(x, y, u, v; kwargs...)
end

max_stars_repo_name: furstj/NSP

max_issues_repo_name: furstj/NSP

max_forks_repo_name: furstj/NSP


text:
# not generating a trajectory yet, let's just test ergodic manager

n = 28
K = 17

d = Domain([0,0,0], [1,1,2*pi], [n,n,36])
em = ErgodicManagerSE2(d, K)

phi = gaussian(d, [0.5,0.5,pi], 0.03*eye(3))
copy!(em.phi, phi)
normalize!(em)
decompose!(em)

phi2 = reconstruct(em, em.phik)

x0 = [0.49, 0.01, 0.01]
dt = 0.5
N = 80
tm = TrajectoryManager(x0, dt, N, ConstantInitializer([0.0,0.0,0.0]))
dynamics!(tm, SingleIntegrator(3,dt))

xd, ud = pto_trajectory(em, tm, max_iters=20)


text:
using Dates
using Random

"""
    ea1p1(instance, solution, Mutation, Fitness, Stop, max_iteration, db)

    This is exactly (1+1)EA algorithm.
"""
function ea1p1(instance:: Dict{String, }, candidate, mutation, fitness, stop, db:: dbEngine)

    best = deepcopy(candidate)

    for iter in 1:BigInt(typemax(Int64))

        candidate = mutation(instance, best, fitness)

        if candidate.fitness >= best.fitness
            best = candidate

            record(db, (id=db.id, best=best.fitness, step=iter))

            if stop(instance, best, iter)
                result(db, (id=db.id,step=iter, value=best.fitness))
                break
            end
        end
    end
    best
end

max_stars_repo_name: ahrajabi/DiscO

max_issues_repo_name: ahrajabi/DiscO

max_forks_repo_name: ahrajabi/DiscO