text:
import KernelDensityEstimate.kde!

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") #cliqMC.prods[1]
        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["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::FactorGraph, bt::BayesTree; N=300)
    len = length(bt.cliques)-1
    vv = Array{Gadfly.Compose.Context,1}(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::FactorGraph, bt::BayesTree;
                        N=300, gt=Union{})
    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["frontalIDs"])

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

        for frid in cliq[2].attributes["frontalIDs"]
            j+=1
            p[j] = getVertKDE(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[dlapi.getvertex(fg,frid).attributes["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, dlapi.getvertex(fg,frid).attributes["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
    #
    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
#
# ## TODO -- you were here with port starboard lines
# function stbPrtLineLayers!(pl, Xpp, Ypp, Thpp)
#
#     l = 5.0
#     lnstpr = [0.0;l;0.0]
#     lnstpg = [0.0;-l;0.0]
#
#     Rd  =SE2(lnstpr)
#     Gr = SE2(lnstpg)
#
#     for i in 1:length(Xpp)
#       lnstt = [Xpp[i];Ypp[i];Thpp[i]]
#       Ps = SE2(lnstt)
#       lnr = se2vee(Ps*Rd)
#       lng = se2vee(Ps*Gr)
#       xsr = [Xpp[i];lnr[1]]
#       ysr = [Ypp[i];lnr[2]]
#       xsg = [Xpp[i];lng[1]]
#       ysg = [Ypp[i];lng[2]]
#
#       push!(pl.layers, layer(x=xsr, y=ysr, Geom.path(), Gadfly.Theme(default_color=colorant"red", line_width=1.5pt))[1] )
#       push!(pl.layers, layer(x=xsg, y=ysg, Geom.path(), Gadfly.Theme(default_color=colorant"green", line_width=1.5pt))[1] )
#     end
#     nothing
# end
#
# # function lblsFromTo(from,to)
# #   lbls=ASCIIString[]
# #   [push!(lbls, "$(i)") for i in from:to]
# #   return lbls
# # end
#
# function drawPoses(fg::FactorGraph; from::Int64=0,to::Int64=99999999,
#                     meanmax=:max, lbls=true, drawhist=true)
#     #Gadfly.set_default_plot_size(20cm, 30cm)
#     Xp,Yp = get2DPoseSamples(fg, from=from, to=to)
#     Xpp = Float64[]; Ypp=Float64[]; Thpp=Float64[]; LBLS=ASCIIString[];
#     if meanmax == :mean
#       Xpp,Ypp, Thpp, LBLS = get2DPoseMeans(fg, from=from, to=to)
#     elseif meanmax == :max
#       Xpp,Ypp, Thpp, LBLS = get2DPoseMax(fg, from=from, to=to)
#     end
#
#     # lbls = lblsFromTo(1,length(Xpp))
#     psplt = Union{}
#     if lbls
#       psplt = Gadfly.plot(
#       Gadfly.layer(x=Xpp,y=Ypp,label=LBLS,Geom.path(), Theme(line_width=2pt), Geom.label)
#       )
#     else
#       psplt = Gadfly.plot(
#       Gadfly.layer(x=Xpp,y=Ypp,Geom.path(), Theme(line_width=2pt))
#       )
#     end
#     stbPrtLineLayers!(psplt, Xpp, Ypp, Thpp)
#     if drawhist
#       push!(psplt.layers,  Gadfly.layer(x=Xp, y=Yp, Geom.histogram2d)[1] )#(xbincount=100, ybincount=100))
#     end
#     psplt
# end
#
# function drawLandms(fg::FactorGraph;
#               from::Int64=0, to::Int64=99999999, minnei::Int64=0,
#               meanmax=:max,lbls=true,showmm=false,drawhist=true,c="red",MM=Union{})
#     #Gadfly.set_default_plot_size(20cm, 30cm)
#     Xp,Yp = get2DLandmSamples(fg, from=from, to=to)
#     Xpp = Float64[]; Ypp=Float64[]; Thpp=Float64[]; lblstags=ASCIIString[];
#     if meanmax==:mean
#       Xpp,Ypp, t, lbltags = get2DLandmMeans(fg, from=from, to=to)
#     elseif meanmax==:max
#       Xpp,Ypp, t, lbltags = get2DLandmMax(fg, from=from, to=to,showmm=showmm,MM=MM)
#     end
#
#     if lbls
#       psplt = Gadfly.plot(
#       Gadfly.layer(x=Xpp,y=Ypp, label=lbltags, Geom.point, Theme(line_width=2pt, default_color=parse(Colorant,c)), Geom.label)
#       # ,Gadfly.layer(x=Xp, y=Yp, Geom.histogram2d)#(xbincount=100, ybincount=100)
#       )
#     else
#       psplt = Gadfly.plot(
#       Gadfly.layer(x=Xpp,y=Ypp, Geom.point, Theme(line_width=2pt, default_color=parse(Colorant,c)))
#       )
#     end
#
#     if drawhist
#       push!(psplt.layers, Gadfly.layer(x=Xp, y=Yp, Geom.histogram2d)[1])#(xbincount=100, ybincount=100)
#     end
#
#     psplt
# end
#
# function drawPosesLandms(fg::FactorGraph;
#                     from::Int64=0, to::Int64=99999999, minnei::Int64=0,
#                     meanmax=:max,lbls=true,drawhist=true, MM=Union{})
#   p = drawPoses(fg, from=from,to=to,meanmax=meanmax,lbls=lbls,drawhist=drawhist)
#   pl = drawLandms(fg, from=from, to=to, minnei=minnei,lbls=lbls,drawhist=drawhist, MM=MM)
#   for l in pl.layers
#     push!(p.layers, l)
#   end
#   return p
# end
#
# function drawSubmaps(fgl::FactorGraph, fromto::Array{Int,2};
#                     m1hist=false,m2hist=false,m3hist=false, showmm=false)
#   p = drawLandms(fgl, from=fromto[1,1], to=fromto[1,2], drawhist=m1hist, showmm=showmm)
#   if size(fromto,1) >1
#     p2 = drawLandms(fgl, from=fromto[2,1], to=fromto[2,2], drawhist=m2hist,c="blue", showmm=showmm)
#     for l in p2.layers
#       push!(p.layers, l)
#     end
#   end
#   if size(fromto,1) >2
#     p3 = drawLandms(fgl, from=fromto[3,1], to=fromto[3,2], drawhist=m3hist,c="magenta", showmm=showmm)
#     for l in p3.layers
#       push!(p.layers, l)
#     end
#   end
#   return p
# end
#
# function drawSubmaps(fgl::FactorGraph, fromto::Array{Int,1}; spread::Int=25,
#                     m1hist=false,m2hist=false,m3hist=false, showmm=false)
#   ft = zeros(Int,length(fromto),2)
#   for i in 1:length(fromto)
#     ft[i,1] = fromto[i]-spread; ft[i,2] = fromto[i]+spread;
#   end
#   drawSubmaps(fgl, ft, m1hist=m1hist, m2hist=m2hist, m3hist=m3hist, showmm=showmm)
# end
#
# # function getKDEMax(p::BallTreeDensity;N=200)
# #   m = zeros(p.bt.dims)
# #   for i in 1:p.bt.dims
# #     mm = marginal(p,[i])
# #     rangeV = getKDERange(mm)
# #     X = linspace(rangeV[1],rangeV[2],N)
# #     yV = evaluateDualTree(mm,X)
# #     m[i] = X[findfirst(yV,maximum(yV))]
# #   end
# #   return m
# # 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{Int,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] = investigatePoseKDE(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 kde!(em::EasyMessage)
  return kde!(em.pts,em.bws)
end

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

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

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

function drawUpMsgAtCliq(fg::FactorGraph, cliq::Graphs.ExVertex)
    for id in keys(cliq.attributes["debug"].outmsg.p)
        print("$(dlapi.getvertex(fg,id).attributes["label"]), ") #fg.v[id].
    end
    println("")
    sleep(0.1)
    potens = cliq.attributes["debug"].outmsg.p
    vArrPotentials(potens)
end

function drawUpMsgAtCliq(fg::FactorGraph, bt::BayesTree, lbl::ASCIIString)
    drawUpMsgAtCliq(fg, whichCliq(bt, lbl) )
end

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

function dwnMsgsAtCliq(fg::FactorGraph, bt::BayesTree, lbl::ASCIIString)
    dwnMsgsAtCliq(fg, whichCliq(bt, 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, investigatePoseKDE([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["debug"]
    sleep(0.1)
    mcmcPose2D!(plots, cliqDbg, iter)
end

function drawUpMCMCPose2D!(plots::Array{Gadfly.Compose.Context,1}, bt::BayesTree, frt::ASCIIString, 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["debugDwn"]
    sleep(0.1)
    mcmcPose2D!(plots, cliqDbg, iter)
end

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

function drawLbl(fg::FactorGraph, lbl::ASCIIString)
    # v = dlapi.getvertex(fg,lbl) # fg.v[fg.IDs[lbl]]
    # investigatePoseKDE(kde!(getVal(v)))
    investigatePoseKDE(getVertKDE(fg,lbl))
end

function predCurrFactorBeliefs(fgl::FactorGraph, fc::Graphs.ExVertex)
  # TODO update to use ls and lsv functions
  prjcurvals = Dict{ASCIIString, Array{BallTreeDensity,1}}()
  for v in dlapi.outneighbors(fgl, fc) #out_neighbors(fc, fgl.g)
    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::FactorGraph, pDens::Dict{Int,EasyMessage}, N=200)
  p = BallTreeDensity[]
  lbls = ASCIIString[]
  for pd in pDens
    push!(p, kde!(pd[2].pts,pd[2].bws))
    push!(lbls, dlapi.getvertex(fgl,pd[1]).attributes["label"]) # fgl.v[pd[1]].
  end
  @show lbls
  drawHorDens(p,N,lbls=lbls)
end

function drawHorBeliefsList(fgl::FactorGraph, lbls::Array{ASCIIString,1};
                        nhor::Int=-1,gt=Union{},N::Int=200, extend=0.1)
  len = length(lbls)
  pDens = BallTreeDensity[]
  for lb in lbls
    ptkde = getVertKDE(fgl,lb)
    push!(pDens, ptkde )
    # pt = getVal(fgl.v[fgl.IDs[lb]])
    # push!(pDens, kde!(pt,"lcv"))
  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 = ASCIIString[]
    for j in 0:(nhor-1)
      if i+j <= len
        push!(pH, pDens[i+j])
        push!(labels, lbls[i+j])
        if gt != Union{} gtvals[j+1] = gt[lbls[i+j]]  end
      end
    end
    vidx+=1
    if gt !=Union{}
      vv[vidx] = KernelDensityEstimate.drawHorDens(pH, N=N, gt=gtvals, lbls=labels, extend=extend)
    else
      vv[vidx] = KernelDensityEstimate.drawHorDens(pH, N=N, lbls=labels, extend=extend)
    end
  end
  vv
end

function drawFactorBeliefs(fgl::FactorGraph, flbl::ASCIIString)
  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.f[fgl.fIDs[flbl]]
      prjcurvals, lbls = predCurrFactorBeliefs(fgl, fc)
      if length(lbls) == 3
        return vstack(
        investigatePoseKDE(prjcurvals[lbls[1]]),
        investigatePoseKDE(prjcurvals[lbls[2]]),
        investigatePoseKDE(prjcurvals[lbls[3]]),
        )
      elseif length(lbls) == 2
        return vstack(
        investigatePoseKDE(prjcurvals[lbls[1]]),
        investigatePoseKDE(prjcurvals[lbls[2]]),
        )
      elseif length(lbls) == 1
        return investigatePoseKDE(prjcurvals[lbls[1]])
      end

    # end
  # end
  nothing
end

function localProduct(fgl::FactorGraph, lbl::ASCIIString;N=300)
  arr = Array{BallTreeDensity,1}()
  cf = ls(fgl, lbl)
  pp = Union{}
  for f in cf
    fgf = dlapi.getvertex(fgl,f,:fnc) #fgl.f[fgl.fIDs[f]]
    push!(arr, kde!(evalFactor2(fgl, fgf, fgl.IDs[lbl])))
  end
  if length(arr)>1
    Ndims = arr[1].bt.dims
    dummy = kde!(rand(Ndims,N),[1.0]);
    pGM, = prodAppxMSGibbsS(dummy, arr, Union{}, Union{}, 15)
    pp = kde!(pGM,"lcv")
  end
  return pp,arr
end

function drawLocalProduct(fgl::FactorGraph, lbl::ASCIIString;N=300)
  arr = Array{BallTreeDensity,1}()
  push!(arr, getVertKDE(fgl, lbl)) #kde!(getVal(fgl.v[fgl.IDs[lbl]]))
  pp, parr = localProduct(fgl, lbl, N=N)
  if pp != Union{}
    push!(arr,pp)
    for a in parr
      push!(arr, a)
    end
  end
  return investigatePoseKDE(arr)
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{FactorGraph,1}, lbl;nw=false)
  len = length(FGL)
  [saveplot(drawLocalProduct(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 ls(fgl::FactorGraph, lbl::ASCIIString)
  ls = ASCIIString[]
  v = Union{}
  if haskey(fgl.IDs, lbl)
    id = fgl.IDs[lbl]
  else
    return ls
  end
  v = dlapi.getvertex(fgl,id) #fgl.v[id]
  for outn in dlapi.outneighbors(fgl, v) # out_neighbors(v, fgl.g)
    push!(ls, outn.label)
  end
  return ls
end

function ls(fgl::FactorGraph)
  k = collect(keys(fgl.IDs))
  l = Int[]
  x = Int[]
  for kk in k
    val = parse(Int,kk[2:end])
    if kk[1] == 'l'
      push!(l,val)
    elseif kk[1] == 'x'
      push!(x,val)
    end
  end
  l = sort(l)
  x = sort(x)
  ll = Array{ASCIIString,1}(length(l))
  xx = Array{ASCIIString,1}(length(x))
  for i in 1:length(l)
    ll[i] = string("l",l[i])
  end
  for i in 1:length(x)
    xx[i] = string("x",x[i])
  end
  return xx,ll
end

function lsv(fgl::FactorGraph, lbl::ASCIIString)
  ls = ASCIIString[]
  v = Union{}
  if haskey(fgl.fIDs, lbl)
    id = fgl.fIDs[lbl]
  else
    return ls
  end
  v = fgl.f[id]
  for outn in dlapi.outneighbors(fgl, v) # out_neighbors(v, fgl.g)
    push!(ls, outn.label)
  end
  return ls
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::FactorGraph, isamdict::Dict{Int,Array{Float64,1}})
  r,rt,lb = computeGraphResiduals(fgl,isamdict);
  fixRotWrapErr!(rt)
  return [sqrt(Base.mean(r.^2));maximum(abs(r));sqrt(Base.mean(rt.^2));maximum(rt)]
end

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

  rr = 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::FactorGraph, lbl::ASCIIString)
  pp, arr = localProduct(fgl, lbl)
  lpm = getKDEMax(pp)
  em = getKDEMax(getVertKDE(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{FactorGraph,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 = 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(Base.mean(ERR.^2)))
    push!(MAX, maximum(abs(ERR)))
    push!(RMSth, sqrt(Base.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 Gadfly.plot(
  # layer(,y=RMS,Geom.path(),Geom.point, Gadfly.Theme(default_color=parse(Colorant,"red"))),
  # layer(x=1:len,y=MAX,Geom.path(),Geom.point, Gadfly.Theme(default_color=parse(Colorant,"black"))),
  # layer(x=1:len,y=RMSth*180.0/pi,Geom.path(),Geom.point, Gadfly.Theme(default_color=parse(Colorant,"blue"))),
  # layer(x=1:len,y=MAXth*180.0/pi,Geom.path(),Geom.point, Gadfly.Theme(default_color=parse(Colorant,"green")),
  # )
  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{FactorGraph,1}, vertid::Int64)
  ret = Array{BallTreeDensity,1}()
  for i in 1:length(fgD)
    push!(ret, getVertKDE(fgD[i],vertid) )
    # V = dlapi.getvertex(fgD[i],vertid)# fgD[i].v[vertid]
    # push!(ret, kde!(getVal(V)))
  end
  return ret
end

function drawAllPose2DBeliefs(plots::Array{Gadfly.Compose.Context,1}, fgD::Array{FactorGraph,1})
    ids = sort(collect(keys(fgD[1].v)))
    co = ["black"; "blue"; "green"; "red"; "magenta"; "cyan"; "cyan1"; "cyan2"]
    println(co[1:length(fgD)])
    # V = fg.v[7]
    for i in ids
        # V = fgD[1].v[i]
        @show dlapi.getvertex(fgD[1],i).attributes["label"] #fgD[1].v[i].
        # if length(fgD) >= 2
            # V0 = fgD[2].v[i]
            kdes = getAllFGsKDEs(fgD, i)
            push!(plots, investigatePoseKDE(  kdes  )) # [kde!(getVal(V)); kde!(getVal(V0))]
        # else
        #     push!(plots, investigatePoseKDE(  [kde!(getVal(V))]  ))
        # end
    end
    vstackedPlots(plots)
end

# legacy function -- use the array version instead
function drawAllPose2DBeliefs(plots::Array{Gadfly.Compose.Context,1}, fgD::FactorGraph, fgD0=Union{})
  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{FactorGraph,1})
    comicA = Array{Gadfly.Plot,1}()
    for fgd in fgD
        cv = drawPosesLandms(fgd)
        # cv = drawPoses(fgd)
        push!(comicA,cv)
    end
    hstack(comicA)
end

# function drawComicStripLM(fgD::Array{FactorGraph,1})
#   len = length(fgD)
#   comicA = Array{Gadfly.Plot,1}(len)
#   r = Array{RemoteRef,1}(len)
#   # for fgd in fgD
#   #     cv = drawPosesLandms(fgd)
#   #     # cv = drawPoses(fgd)
#   #     push!(comicA,cv)
#   # end
#   [r[i] = @spawn drawPosesLandms(fgD[i]) for i in 1:len]
#   [comicA[i] = fetch(r[i]) for i in 1:len]
#   hstack(comicA)
# end

function drawComicStrip(fgD::Array{FactorGraph,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{FactorGraph,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{FactorGraph,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


text:
import PTChain
import ProposalFunctions

f(x) = sum(x.^2)

@timev x = PTChain.sampleTempering(randn(3), f, [0.1, 1], ProposalFunctions.widescaleRandomNoise, 5, 10000)

ext: jl

max_stars_repo_name: jacob-ruth/pt-optimization

max_issues_repo_name: jacob-ruth/pt-optimization

max_forks_repo_name: jacob-ruth/pt-optimization


text:
include("building_tree.jl")
include("building_weighted_tree.jl")
include("utilities.jl")
using LinearAlgebra

function fusion_points(X, Y, proportion)
  n = length(Y)
  m = n*proportion
  m ÷= 1
  for i in 1:m
    # trouver points d'indice a et b les plus proches de meme classe
    distance = 1
    a, b = 0, 0
    for k in 1:size(X, 1)
      for l in 1:size(X, 1)
        if Y[k] == Y[l] && k != l
          distance_provisoire = norm(X[k,:] - X[l,:]) # distance euclidienne
          if distance_provisoire < distance && distance_provisoire > 0
            distance = distance_provisoire
            a = k
            b = l
          end
        end
      end
    end
    # les fusionner dans X et dans Y
    X[a,:] = 0.5*(X[a,:] + X[b,:])
    # supprimer X[b] et Y[b]
    X = X[1:end .!= b,:]
    splice!(Y, b)
  end
  return X,Y
end

function main()

    # Pour chaque jeu de données
    for dataSetName in ["iris", "seeds", "wine", "haberman", "banknote"]

        print("=== Dataset ", dataSetName)

        # Préparation des données
        include("../data/" * dataSetName * ".txt")
        train, test = train_test_indexes(length(Y))
        X_train = X[train, :] # modifier X_train pour fusionner gamma pourcent des points de meme classe
        Y_train = Y[train] # idem sur Y_train
        X_train, Y_train = fusion_points(X_train, Y_train, 0.2)
        X_test = X[test, :] # mais pas sur les donnees de test !
        Y_test = Y[test]

        println(" (train size ", size(X_train, 1), ", test size ", size(X_test, 1), ", ", size(X_train, 2), ", features count: ", size(X_train, 2), ")")

        # Temps limite de la méthode de résolution en secondes
        time_limit = 30

        # On ajoute ou pas des poids aléatoires aux points/aux clusters
        weigthed = true
        if weigthed
          println("avec poids : ")
        end

        # Pour chaque profondeur considérée
        for D in 2:4

            println("  D = ", D)

            ## 1 - Univarié (séparation sur une seule variable à la fois)
            # Création de l'arbre
            print("    Univarié...  \t")
            if weigthed
              T, obj, resolution_time, gap = build_weighted_tree(X_train, Y_train, D,  multivariate = false, time_limit = time_limit)
            else
              T, obj, resolution_time, gap = build_tree(X_train, Y_train, D,  multivariate = false, time_limit = time_limit)
            end

            # Test de la performance de l'arbre
            print(round(resolution_time, digits = 1), "s\t")
            print("gap ", round(gap, digits = 1), "%\t")
            if T != nothing
                print("Erreurs train/test ", prediction_errors(T,X_train,Y_train))
                print("/", prediction_errors(T,X_test,Y_test), "\t")
            end
            println()

            ## 2 - Multivarié
            print("    Multivarié...\t")
            if weigthed
              T, obj, resolution_time, gap = build_weighted_tree(X_train, Y_train, D, multivariate = true, time_limit = time_limit)
            else
              T, obj, resolution_time, gap = build_tree(X_train, Y_train, D, multivariate = true, time_limit = time_limit)
            end
            print(round(resolution_time, digits = 1), "s\t")
            print("gap ", round(gap, digits = 1), "%\t")
            if T != nothing
                print("Erreurs train/test ", prediction_errors(T,X_train,Y_train))
                print("/", prediction_errors(T,X_test,Y_test), "\t")
            end
            println("\n")
        end
    end
end

max_stars_repo_name: boumpteryx/RODM

max_issues_repo_name: boumpteryx/RODM

max_forks_repo_name: boumpteryx/RODM


text:
colscitypes = Dict(
    :family_id => Union{Missing, Multiclass},
    :person_id => Union{Missing, Multiclass},
    :centre => Union{Missing, Multiclass},
    :status => Union{Missing, Multiclass},
    :dob => Union{Missing, Multiclass},
    :age => Union{Missing, MLJ.Continuous},
    :sex => Union{Missing, Multiclass},
    :assessment_date => Union{Missing, Multiclass},
    :most_recent_update => Union{Missing, Multiclass},
    :affected_status => Union{Missing, Multiclass},
    :diagnosis => Union{Missing, Multiclass},
    :gaf_last_assessment => Union{Missing, MLJ.Continuous},
    :marital_status => Union{Missing, Multiclass},
    :ses => Union{Missing, Multiclass},
    :father_ethnicity => Union{Missing, Multiclass},
    :mother_ethnicity => Union{Missing, Multiclass},
    :primary_diagnosis => Union{Missing, Multiclass},
    :age_of_onset => Union{Missing, MLJ.Continuous},
    :onset_d => Union{Missing, MLJ.Continuous},
    :onset_M => Union{Missing, MLJ.Continuous},
    :onset_m => Union{Missing, MLJ.Continuous},
    :clinical_course => Union{Missing, Multiclass},
    :lifetime_manias => Union{Missing, MLJ.Continuous},
    :lifetime_depressions => Union{Missing, MLJ.Continuous},
    :lifetime_mixed => Union{Missing, MLJ.Continuous},
    :lifetime_multiphasic => Union{Missing, MLJ.Continuous},
    :total_lifetime_episodes => Union{Missing, MLJ.Continuous},
    :polarity_first_episode => Union{Missing, Multiclass},
    :rapid_cycling => Union{Missing, Multiclass},
    :lifetime_psychosis => Union{Missing, Multiclass},
    :rapid_mood_switching => Union{Missing, Multiclass},
    :first_degree_rel_mood_do => Union{Missing, Multiclass},
    :first_degree_rel_bd => Union{Missing, Multiclass},
    :n_first_degree_rel_bd1 => Union{Missing, MLJ.Continuous},
    :n_first_degree_rel_bd2 => Union{Missing, MLJ.Continuous},
    :n_first_degree_rel_unipolar_depression => Union{Missing, MLJ.Continuous},
    :n_first_degree_rel_sza => Union{Missing, MLJ.Continuous},
    :n_first_degree_rel_scz => Union{Missing, MLJ.Continuous},
    :n_first_degree_rel_anx => Union{Missing, MLJ.Continuous},
    :n_first_degree_rel_unaffected => Union{Missing, MLJ.Continuous},
    :n_first_degree_rel_completed_suicide => Union{Missing, MLJ.Continuous},
    :n_first_degree_rel_suicide_attempt => Union{Missing, MLJ.Continuous},
    :n_second_degree_rel_completed_suicide => Union{Missing, MLJ.Continuous},
    :n_second_degree_rel_suicide_attempt => Union{Missing, MLJ.Continuous},
    :total_n_first_degree_rel => Union{Missing, MLJ.Continuous},
    :n_episodes_pre_li => Union{Missing, MLJ.Continuous},
    :n_episodes_on_li => Union{Missing, MLJ.Continuous},
    :li_response_narrow => Union{Missing, Multiclass},
    :li_a_score => Union{Missing, MLJ.Continuous},
    :li_b1_score => Union{Missing, MLJ.Continuous},
    :li_b2_score => Union{Missing, MLJ.Continuous},
    :li_b3_score => Union{Missing, MLJ.Continuous},
    :li_b4_score => Union{Missing, MLJ.Continuous},
    :li_b5_score => Union{Missing, MLJ.Continuous},
    :li_total_score => Union{Missing, MLJ.Continuous},
    :n_suicide_attempts => Union{Missing, MLJ.Continuous},
    :n_potentially_lethal_suicide_attempts => Union{Missing, MLJ.Continuous},
    :age_first_suicide_attempt => Union{Missing, MLJ.Continuous},
    :mood_polarity => Union{Missing, Multiclass},
    :lifetime_hx_suicide_ideation => Union{Missing, Multiclass},
    :suicide_ideation_related_to_mood_episode => Union{Missing, Multiclass},
    :social_anxiety => Union{Missing, Multiclass},
    :age_onset_social_anxiety => Union{Missing, MLJ.Continuous},
    :panic_disorder => Union{Missing, Multiclass},
    :age_onset_panic_disorder => Union{Missing, MLJ.Continuous},
    :generalized_anxiety => Union{Missing, Multiclass},
    :age_onset_generalized_anxiety => Union{Missing, MLJ.Continuous},
    :ocd => Union{Missing, Multiclass},
    :age_onset_ocd => Union{Missing, MLJ.Continuous},
    :substance_abuse => Union{Missing, Multiclass},
    :age_onset_substance_abuse => Union{Missing, MLJ.Continuous},
    :adhd => Union{Missing, Multiclass},
    :age_onset_adhd => Union{Missing, MLJ.Continuous},
    :learning_disability => Union{Missing, Multiclass},
    :age_onset_learning_disability => Union{Missing, MLJ.Continuous},
    :primary_insomnia => Union{Missing, Multiclass},
    :age_onset_primary_insomnia => Union{Missing, MLJ.Continuous},
    :personality_disorder => Union{Missing, Multiclass},
    :age_onset_personality_disorder => Union{Missing, MLJ.Continuous},
    :diabetes_mellitus => Union{Missing, Multiclass},
    :age_onset_diabetes_mellitus => Union{Missing, MLJ.Continuous},
    :diabetes_mellitus_type => Union{Missing, Multiclass},
    :hypertension => Union{Missing, Multiclass},
    :age_onset_hypertension => Union{Missing, MLJ.Continuous},
    :menstrual_abnormality => Union{Missing, Multiclass},
    :age_onset_menstrual_abnormality => Union{Missing, MLJ.Continuous},
    :thyroid_disease => Union{Missing, Multiclass},
    :age_onset_thyroid_disease => Union{Missing, MLJ.Continuous},
    :head_injury => Union{Missing, Multiclass},
    :age_onset_head_injury => Union{Missing, MLJ.Continuous},
    :migraine => Union{Missing, Multiclass},
    :age_onset_migraine => Union{Missing, MLJ.Continuous},
    :opcrit_source_of_rating => Union{Missing, Multiclass},
    :opcrit_time_frame => Union{Missing, Multiclass},
    :opcrit_gender => Union{Missing, Multiclass},
    :opcrit_age_of_onset => Union{Missing, MLJ.Continuous},
    :opcrit_mode_of_onset => Union{Missing, Multiclass},
    :opcrit_single => Union{Missing, Multiclass},
    :opcrit_unemployed_at_onset => Union{Missing, Multiclass},
    :opcrit_illness_duration => Union{Missing, MLJ.Continuous},
    :opcrit_poor_work_adjustment => Union{Missing, Multiclass},
    :opcrit_poor_premorbid_social_adjustment => Union{Missing, Multiclass},
    :opcrit_premorbid_personality_disorder => Union{Missing, Multiclass},
    :opcrit_alcohol_drug_abuse => Union{Missing, Multiclass},
    :opcrit_famhx_scz => Union{Missing, Multiclass},
    :opcrit_famhx_other_psychiatric_disorder => Union{Missing, Multiclass},
    :opcrit_premorbid_coarse_brain_disease => Union{Missing, Multiclass},
    :opcrit_premorbid_definite_psychosocial_stressor => Union{Missing, Multiclass},
    :opcrit_bizarre_behaviour => Union{Missing, Multiclass},
    :opcrit_catatonia => Union{Missing, Multiclass},
    :opcrit_excessive_activity => Union{Missing, Multiclass},
    :opcrit_reckless_activity => Union{Missing, Multiclass},
    :opcrit_distractibility => Union{Missing, Multiclass},
    :opcrit_reduced_need_for_sleep => Union{Missing, Multiclass},
    :opcrit_agitated_activity => Union{Missing, Multiclass},
    :opcrit_slowed_activity => Union{Missing, Multiclass},
    :opcrit_loss_energy_tiredness => Union{Missing, Multiclass},
    :opcrit_speech_difficult_to_understand => Union{Missing, Multiclass},
    :opcrit_incoherent => Union{Missing, Multiclass},
    :opcrit_positive_formal_thought_disorder => Union{Missing, Multiclass},
    :opcrit_negative_formal_thought_disorder => Union{Missing, Multiclass},
    :opcrit_pressured_speech => Union{Missing, Multiclass},
    :opcrit_thoughts_racing => Union{Missing, Multiclass},
    :opcrit_restricted_affect => Union{Missing, Multiclass},
    :opcrit_blunted_affect => Union{Missing, Multiclass},
    :opcrit_inappropriate_affect => Union{Missing, Multiclass},
    :opcrit_elevated_mood => Union{Missing, Multiclass},
    :opcrit_irritable_mood => Union{Missing, Multiclass},
    :opcrit_dysphoria => Union{Missing, Multiclass},
    :opcrit_diurnal_variation => Union{Missing, Multiclass},
    :opcrit_loss_of_pleasure => Union{Missing, Multiclass},
    :opcrit_altered_libido => Union{Missing, Multiclass},
    :opcrit_poor_concentration => Union{Missing, Multiclass},
    :opcrit_excessive_self_reproach => Union{Missing, Multiclass},
    :opcrit_suicidal_ideation => Union{Missing, Multiclass},
    :opcrit_initial_insomnia => Union{Missing, Multiclass},
    :opcrit_middle_insomnia => Union{Missing, Multiclass},
    :opcrit_early_morning_waking => Union{Missing, Multiclass},
    :opcrit_excessive_sleep => Union{Missing, Multiclass},
    :opcrit_poor_appetite => Union{Missing, Multiclass},
    :opcrit_weight_loss => Union{Missing, Multiclass},
    :opcrit_increased_appetite => Union{Missing, Multiclass},
    :opcrit_weight_gain => Union{Missing, Multiclass},
    :opcrit_relationship_between_psychotic_and_affective_symptoms => Union{Missing, Multiclass},
    :opcrit_increased_sociability => Union{Missing, Multiclass},
    :opcrit_persecutory_delusions => Union{Missing, Multiclass},
    :opcrit_well_organized_delusions => Union{Missing, Multiclass},
    :opcrit_increased_self_esteem => Union{Missing, Multiclass},
    :opcrit_grandiose_delusions => Union{Missing, Multiclass},
    :opcrit_delusions_of_influence => Union{Missing, Multiclass},
    :opcrit_bizarre_delusions => Union{Missing, Multiclass},
    :opcrit_widespread_delusions => Union{Missing, Multiclass},
    :opcrit_delusions_of_passivity => Union{Missing, Multiclass},
    :opcrit_primary_delusional_perception => Union{Missing, Multiclass},
    :opcrit_other_primary_delusions => Union{Missing, Multiclass},
    :opcrit_delusions_hallucinations_one_week => Union{Missing, Multiclass},
    :opcrit_persecutory_jealous_delusion_hallucination => Union{Missing, Multiclass},
    :opcrit_thought_insertion => Union{Missing, Multiclass},
    :opcrit_thought_withdrawal => Union{Missing, Multiclass},
    :opcrit_thought_broadcast => Union{Missing, Multiclass},
    :opcrit_delusions_of_guilt => Union{Missing, Multiclass},
    :opcrit_delusions_of_poverty => Union{Missing, Multiclass},
    :opcrit_nihilistic_delusions => Union{Missing, Multiclass},
    :opcrit_thought_echo => Union{Missing, Multiclass},
    :opcrit_third_person_auditory_hallucinations => Union{Missing, Multiclass},
    :opcrit_running_commentary_voices => Union{Missing, Multiclass},
    :opcrit_abusive_accusatory_persecutory_voices => Union{Missing, Multiclass},
    :opcrit_other_nonaffective_auditory_hallucinations => Union{Missing, Multiclass},
    :opcrit_nonaffective_hallucination_any_modality => Union{Missing, Multiclass},
    :opcrit_lifetime_dx_alcohol_abuse_dependence => Union{Missing, Multiclass},
    :opcrit_lifetime_diagnosis_cannabis_dependence => Union{Missing, Multiclass},
    :opcrit_lifetime_diagnosis_other_abuse_dependence => Union{Missing, Multiclass},
    :opcrit_alcohol_abuse_dependence_with_psychopathology => Union{Missing, Multiclass},
    :opcrit_cannabis_abuse_dependence_with_psychopathology => Union{Missing, Multiclass},
    :opcrit_other_abuse_dependence_with_psychopathology => Union{Missing, Multiclass},
    :opcrit_information_not_credible => Union{Missing, Multiclass},
    :opcrit_lack_of_insight => Union{Missing, Multiclass},
    :opcrit_rapport_difficult => Union{Missing, Multiclass},
    :opcrit_impairment_incapacity_during_disorder => Union{Missing, Multiclass},
    :opcrit_deterioration_from_premorbid_level_functioning => Union{Missing, Multiclass},
    :opcrit_psychotic_symptoms_respond_to_neuroleptics => Union{Missing, Multiclass},
    :opcrit_course_of_disorder => Union{Missing, Multiclass},
    :vpa_response_narrow => Union{Missing, Multiclass},
    :vpa_a_score => Union{Missing, MLJ.Continuous},
    :vpa_b1_score => Union{Missing, MLJ.Continuous},
    :vpa_b2_score => Union{Missing, MLJ.Continuous},
    :vpa_b3_score => Union{Missing, MLJ.Continuous},
    :vpa_b4_score => Union{Missing, MLJ.Continuous},
    :vpa_b5_score => Union{Missing, MLJ.Continuous},
    :vpa_total_score => Union{Missing, MLJ.Continuous},
    :ltg_response_narrow => Union{Missing, Multiclass},
    :ltg_a_score => Union{Missing, MLJ.Continuous},
    :ltg_b1_score => Union{Missing, MLJ.Continuous},
    :ltg_b2_score => Union{Missing, MLJ.Continuous},
    :ltg_b3_score => Union{Missing, MLJ.Continuous},
    :ltg_b4_score => Union{Missing, MLJ.Continuous},
    :ltg_b5_score => Union{Missing, MLJ.Continuous},
    :ltg_total_score => Union{Missing, MLJ.Continuous}
);

"""
    coerce_scitypes(data)

Function to coerce to appropriate scientific types, such that it can be piped.
"""
coerce_scitypes(data::DataFrame) = coerce(data, colscitypes)

"""
    remove_columns(data, colregex[, fuzzy])

Selects all columns in `data` except those matching the `colregex` pattern. If `fuzzy=false`
then only columns matching `colregex` exactly will be removed.
"""
function remove_columns(data::DataFrame, colregex::String, fuzzy::Bool=true)
    if fuzzy == true
        out = select(data, Regex("^((?!" * colregex *").)*\$"))
    else
        out = select(data, Regex("^((?!" * colregex *"))"))
    end
    return out
end

remove_ltg_targets(data::DataFrame) = remove_columns(data, "ltg_")
remove_vpa_targets(data::DataFrame) = remove_columns(data, "vpa_")
remove_opcrit(data::DataFrame) = remove_columns(data, "opcrit_")

"""
    remove_n_episode_outliers(data, cols)

Removes values specified as >=99, which in the dataset for many variables means "too high to count"
"""
function remove_n_episode_outliers(data::DataFrame)
    data[:, :n_episodes_pre_li] = passmissing(x -> x >=99 ? missing : x).(data[:, :n_episodes_pre_li])
    return data
end

"""
    probands_only(data)

Excludes any relatives from the dataset.
"""
function probands_only(data::DataFrame)
    data = data[.!ismissing.(data.status), :]
    data = data[data.status .== "Proband", :]
    return select(data, Not(:status))
end

"""
    limit_to_bd_diagnosis(data)

Includes only subjects with BDI or II
"""
function limit_to_bd_diagnosis(data::DataFrame)
    data[ismissing.(data.diagnosis), :diagnosis] .= "FILLER"
    filter!(:diagnosis => x -> (x == "BD I") || (x == "BD II"), data)
    return data
end

"""
    merge_mania_hypomania_ageofonset(data)

Combines the ages of mania and hypomania onset into a single variable where the minimum is
used as the age of onset for the union of these phases
"""
function merge_mania_hypomania_ageofonset(data::DataFrame)
    data[ismissing.(data.onset_M), "onset_M"] .= 10000
    data[ismissing.(data.onset_m), "onset_m"] .= 10000
    insertcols!(data, :onset_Mm=> min.(data.onset_M, data.onset_m))
    allowmissing!(data)
    data[data.onset_Mm .== 10000,"onset_Mm"] .= missing
    data = select(data, Not(:onset_M))
    data = select(data, Not(:onset_m))
    return data
end

"""
    normalize_family_history(data)

Normalizes family history counts to [0,1]. If no first degree relatives, value is 0.
"""
function normalize_family_history(data::DataFrame)
    # Dictionary mapping count variable names to the proportion
    fdr_dict = Dict(
        "n_first_degree_rel_anx"=>"p_fdr_anx",
        "n_first_degree_rel_bd1"=>"p_fdr_bd1",
        "n_first_degree_rel_bd2"=>"p_fdr_bd2",
        "n_first_degree_rel_completed_suicide"=>"p_fdr_suic",
        "n_first_degree_rel_suicide_attempt"=>"p_fdr_satt",
        "n_first_degree_rel_scz"=>"p_fdr_scz",
        "n_first_degree_rel_sza"=>"p_fdr_sza",
        "n_first_degree_rel_unaffected"=>"p_fdr_unaff",
        "n_first_degree_rel_unipolar_depression"=>"p_fdr_mdd")
    nfdr = data[:,"total_n_first_degree_rel"]
    for k ∈ keys(fdr_dict)
        p_fdr = passmissing((x,y) -> y == 0 ? 0.0 : x/y).(data[:,k], nfdr)
        insertcols!(data, fdr_dict[k]=> p_fdr)
        data = select(data, Not(k))
    end
    return select(data, Not("total_n_first_degree_rel"))
end

"""
    narrow_comorbidity_definition(data)

For the comorbidity fields `"social_anxiety", "panic_disorder", "generalized_anxiety",
"ocd", "substance_abuse", "adhd","learning_disability", "primary_insomnia",
"personality_disorder"` we set a binary indicator as positive even if they are subsyndromal.
"""
function narrow_comorbidity_definition(data::DataFrame)
    allowmissing(data)
    comorb = ["social_anxiety", "panic_disorder", "generalized_anxiety",
              "ocd", "substance_abuse", "adhd","learning_disability",
              "primary_insomnia", "personality_disorder"]
    for c ∈ comorb
        data[:,c] = passmissing(x -> replace(x, "Subsyndromal"=>"Yes")).(data[:,c])
        data[:,c] = passmissing(x -> replace(x, "Full criteria"=>"Yes")).(data[:,c])
    end

    # Remove comorbidity ages of onset
    data = remove_columns(data, "age_onset_")
    # Remove diabetes mellitus type
    data = remove_columns(data, "diabetes_mellitus_type")
    return data
end

"""
    narrow_comorbidity_definition(data)

Removes variables that are excluded from the analysis based on being (A) after medication
initiation, or deemed irrelevant a priori. These variables are "dob", "assessment_date",
"most_recent_update", "affected_status", "gaf_last_assessment", "marital_status", "ses",
"primary_diagnosis", "n_episodes_on_li", "father_ethnicity", "mother_ethnicity",
"diagnosis", "total_lifetime_episodes", "n_second_degree_rel_completed_suicide",
"n_second_degree_rel_suicide_attempt", "n_suicide_attempts", "lifetime_manias",
"lifetime_depressions", "lifetime_mixed", "lifetime_multiphasic", "lifetime_psychosis", and
"lifetime_hx_suicide_ideation".
"""
function remove_excluded_variables(data::DataFrame)
    excluded_vars = ["dob",
                     "assessment_date",
                     "most_recent_update",
                     "affected_status",
                     "gaf_last_assessment",
                     "marital_status",
                     "ses",
                     "primary_diagnosis",
                     "n_episodes_on_li",
                     "father_ethnicity",
                     "mother_ethnicity",
                     "diagnosis",
                     "total_lifetime_episodes",
                     "n_second_degree_rel_completed_suicide",
                     "n_second_degree_rel_suicide_attempt",
                     "n_suicide_attempts",
                     "lifetime_manias",
                     "lifetime_depressions",
                     "lifetime_mixed",
                     "lifetime_multiphasic",
                     "lifetime_psychosis",
                     "lifetime_hx_suicide_ideation"]
    for v ∈ excluded_vars
        data = select(data, Not(v))
    end
    return data
end

"""
    remove_ids(data)

Removes family and person ID's.
"""
function remove_ids(data::DataFrame)
    for v ∈ ["family_id", "person_id"]
        data = select(data, Not(v))
    end
    return data
end

"""
    replace_centre_names(data)

Replace centre names with their acronyms.
"""
function replace_centre_names(data::DataFrame)
    data[:,:centre] = replace.(data[:,:centre], "Ottawa/Hamilton"=>"ON")
    data[:,:centre] = replace.(data[:,:centre], "Montreal"=>"MTL")
    data[:,:centre] = replace.(data[:,:centre], "Maritimes"=>"MAR")
    data[:,:centre] = replace.(data[:,:centre], "Cagliari (Centro Bini)"=>"CB")
    data[:,:centre] = replace.(data[:,:centre], "Cagliari (University)"=>"CU")
    data[:,:centre] = replace.(data[:,:centre], "Cagliari (PLOT-BD)"=>"CP")
    data[:,:centre] = replace.(data[:,:centre], "Poznan"=>"POZ")
    data[:,:centre] = replace.(data[:,:centre], "Austria"=>"AUT")
    data[:,:centre] = replace.(data[:,:centre], "Barcelona (Colom)"=>"BARC")
    data[:,:centre] = replace.(data[:,:centre], "Barcelona (Murru)"=>"BARM")
    return data
end

"""
    remove_incomplete_targets(data, regex_str)

Includes only those rows that have complete records for the specified column
"""
function remove_incomplete_targets(data::DataFrame, regex_str::Regex)
    return data[completecases(select(data, regex_str)),:]
end

remove_incomplete_li_target_binary(data::DataFrame) = remove_incomplete_targets(data, r"li_response_narrow")
remove_incomplete_li_target_numerical(data::DataFrame) = remove_incomplete_targets(data, r"li_[b|a]")

"""
    remove_lowvar(X, v)

    Removes variables whose variance is <v or NaN.
"""
Base.isnan(x::Nothing) = false
function remove_lowvar(X::DataFrame, v::Float64=1e-3)
    feat = describe(X, :std)
    return select(X, filter(:std => x -> !any(f -> f(x), (y-> !isnothing(y) ? y < v : false, ismissing, Base.isnan)), feat).variable)
end

"""
    compute_alda_b_score(data[, med])

Sum over Alda B scores for `med ∈ ["li", "vpa", "ltg"]`
"""
function compute_alda_b_score(data::DataFrame, med::String="li")
    regex_str = Regex(med*"_b")
    return sum(eachcol(select(data, regex_str)))
end

"""
    compute_alda_total_score(data[, med])

Sum over Alda total scores.
"""
function compute_alda_total_score(data::DataFrame, med::String="li")
    b_score = compute_alda_b_score(data, med)
    data[:, med*"_total_score"] .= data[:, med*"_a_score"] - b_score
    return data
end

# Splits the data into features, targets and confounds (centre)
function unpack_features_targets_covariates(
        data::Union{DataFrame, SubDataFrame},
        target::String,
        thresh::Union{Nothing, Float64, Int64}=nothing
    )
    if target == "li_response_binary"
        Y = data[:, "li_total_score"] |> Vector{Float64}
        Y = map(x -> x >= thresh ? 1. : 0., Y)
    elseif target == "li_nonresponse_binary"
        Y = data[:, "li_total_score"] |> Vector{Float64}
        Y = map(x -> thresh >= x ? 1. : 0., Y)
    else
        Y = data[:, target] |> Vector{Float64}
    end
    Z = select(data, "centre")
    X = remove_columns(remove_columns(data, "li_"), "centre")
    return X, Y, Z
end

"""
    prune_missingness(data, p)

Function that identifies which features have less than p (proportion) of values missing.
"""
function prune_missingness(data::DataFrame, p::Float64=0.6)
    return describe(data, :nmissing) |> x -> x[x[:,:nmissing]./size(data,1) .< p,:variable]
end

"""
    prune_sparse_observations(X, prune_p)

Removes subjects with more than prune_p amount of variables missing
"""
function prune_sparse_observations(X::DataFrame, prune_p::Float64=1.0)
    # Count missingness by subject
    M = ismissing.(data) |> Matrix{Bool}
    v = sum(M, dims=2)./size(M, 2) |> vec
    return X[v .<= prune_p, :]
end


text:
################################################################################
############################# QUANTUM  CIRCUIT #################################
################################################################################
abstract type QuantumCircuit end
abstract type ParametricQuantumCircuit <: QuantumCircuit end

#  QuantumCircuit
struct QulacsQuantumCircuit <: QuantumCircuit
    pyobj::PyObject
end

function QulacsQuantumCircuit(n_qubit::Int)
    QulacsQuantumCircuit(qulacs.QuantumCircuit(n_qubit))
end

function Base.copy(circuit::QulacsQuantumCircuit)
    QulacsQuantumCircuit(circuit.pyobj.copy())
end

function add_X_gate!(circuit::QulacsQuantumCircuit, idx_qubit::Int)
    circuit.pyobj.add_X_gate(idx_qubit-1)
end

function add_H_gate!(circuit::QulacsQuantumCircuit, idx_qubit::Int)
    circuit.pyobj.add_H_gate(idx_qubit-1)
end

function add_S_gate!(circuit::QulacsQuantumCircuit, idx_qubit::Int)
    circuit.pyobj.add_S_gate(idx_qubit-1)
end

function add_Sdag_gate!(circuit::QulacsQuantumCircuit, idx_qubit::Int)
    circuit.pyobj.add_Sdag_gate(idx_qubit-1)
end

function add_CNOT_gate!(circuit::QulacsQuantumCircuit, control::Int, target::Int)
    circuit.pyobj.add_CNOT_gate(control-1, target-1)
end

function add_RY_gate!(circuit::QulacsQuantumCircuit, idx_qubit::Int, angle::Float64)
    circuit.pyobj.add_RY_gate(idx_qubit-1, angle)
end

function add_RZ_gate!(circuit::QulacsQuantumCircuit, idx_qubit::Int, angle::Float64)
    circuit.pyobj.add_RZ_gate(idx_qubit-1, angle)
end

function add_SWAP_gate!(circuit::QulacsQuantumCircuit, idx_qubit_1::Int, idx_qubit_2::Int)
    circuit.pyobj.add_SWAP_gate(idx_qubit_1-1, idx_qubit_2-1)
end

function add_U1_gate!(circuit::QulacsQuantumCircuit,idx_qubit::Int,angle::Float64)
    circuit.pyobj.add_U1_gate(idx_qubit-1, angle)
end

################################################################################
############################# Quantum Parametric CIRCUIT #######################
################################################################################
struct QulacsParametricQuantumCircuit <: ParametricQuantumCircuit　#QulacsParametricQuantumCircuitはjulia
    pyobj::PyObject
end

function QulacsParametricQuantumCircuit(n_qubit::Int)
    QulacsParametricQuantumCircuit(qulacs.ParametricQuantumCircuit(n_qubit))
end

function Base.copy(circuit::QulacsParametricQuantumCircuit)
    QulacsParametricQuantumCircuit(circuit.pyobj.copy())
end

function add_parametric_multi_Pauli_rotation_gate!(circuit::QulacsParametricQuantumCircuit,
    pauli_indices::Vector{Int}, pauli_ids::Vector{PauliID}, theta_init::Float64=0.0)
    if !all(pauli_indices .>= 1)
        error("pauli indices are out of range!")
    end
    circuit.pyobj.add_parametric_multi_Pauli_rotation_gate(
        pauli_indices .- 1, Int.(pauli_ids), theta_init)
end

function set_parameter!(circuit::QulacsParametricQuantumCircuit, index, theta)
    if index <= 0 || index > get_parameter_count(circuit)
        error("index is out of range!")
    end
    circuit.pyobj.set_parameter(index-1, theta)
end

#=
ispython(circuit::ParametricQuantumCircuit) = isa(circuit, PyObject)

function set_parameter!(circuit::ParametricQuantumCircuit, index, theta)
    if index <= 0 || index > get_parameter_count(circuit)
        error("index is out of range!")
    end
    if ispython(circuit)
        circuit.pyobj.set_parameter(index, theta)
    else
        # ここは Julia の実装
    end
end
=#

function get_parameter_count(circuit::QulacsParametricQuantumCircuit)::Int64
    circuit.pyobj.get_parameter_count()
end

function get_parameter(circuit::QulacsParametricQuantumCircuit, idx::Int)
    @assert idx >= 1
    circuit.pyobj.get_parameter(idx-1)
end

function add_parametric_RX_gate!(circuit::QulacsParametricQuantumCircuit, i::Int, angle::Float64)
    circuit.pyobj.add_parametric_RX_gate(i-1, angle)
end

function add_parametric_RY_gate!(circuit::QulacsParametricQuantumCircuit, i::Int, angle::Float64)
    circuit.pyobj.add_parametric_RY_gate(i-1, angle)
end

function add_parametric_RZ_gate!(circuit::QulacsParametricQuantumCircuit, i::Int, angle::Float64)
    circuit.pyobj.add_parametric_RZ_gate(i-1, angle)
end

function add_S_gate!(circuit::QulacsParametricQuantumCircuit, idx_qubit::Int)
    circuit.pyobj.add_S_gate(idx_qubit-1)
end

function add_Sdag_gate!(circuit::QulacsParametricQuantumCircuit, idx_qubit::Int)
    circuit.pyobj.add_Sdag_gate(idx_qubit-1)
end

function add_Z_gate!(circuit::QulacsParametricQuantumCircuit, idx_qubit::Int)
    circuit.pyobj.add_Z_gate(idx_qubit-1)
end

function add_X_gate!(circuit::QulacsParametricQuantumCircuit, idx_qubit::Int)
    circuit.pyobj.add_X_gate(idx_qubit-1)
end

function add_CNOT_gate!(circuit::QulacsParametricQuantumCircuit, control::Int, target::Int)
    circuit.pyobj.add_CNOT_gate(control-1, target-1)
end

function add_RY_gate!(circuit::QulacsParametricQuantumCircuit, idx_qubit::Int, angle::Float64)
    circuit.pyobj.add_RY_gate(idx_qubit-1, angle)
end

function add_RZ_gate!(circuit::QulacsParametricQuantumCircuit, idx_qubit::Int, angle::Float64)
    circuit.pyobj.add_RZ_gate(idx_qubit-1, angle)
end

function add_U1_gate!(circuit::QulacsParametricQuantumCircuit,idx_qubit::Int,angle::Float64)
    circuit.pyobj.add_U1_gate(idx_qubit-1, angle)
end

#function add_SWAP_gate!(circuit::QulacsPatametricQuantumCircuit, idx_qubit_1::Int, idx_qubit_2::Int)
#    circuit.pyobj.add_SWAP_gate(idx_qubit_1-1, idx_qubit_2-1)
#end

#function create_simul_gate(circuit::QulacsParametricQuantumCircuit,idx_qubit::Int, thetas)
#  theta1 = thetas[1]
#  theta2 = thetas[2]
#  #CNOT02 = CNOT(first, second)
#  add_parametric_RZ_gate!(circuit,idx_qubit, -theta[2])
#  add_parametric_RY_gate!(circuit,idx_qubit, theta[1])
#  add_parametric_RZ_gate!(circuit,idx_qubit, theta[2])
#  add_Z_gate!(circuit, idx_qubit)
#  return merge(merge(merge(merge(CNOT02, CRY0), CNOT20), CRY0_minus), CNOT02)

################################################################################
##################  Variational QUANTUM  CIRCUIT ###############################
################################################################################
abstract type VariationalQuantumCircuit end

"""
Return the number of independent variational parameters
"""
function num_theta(circuit::VariationalQuantumCircuit)::Int64
    0
end

"""
Return a copy of variational parameters
"""
function get_thetas(circuit::VariationalQuantumCircuit)::Vector{Float64}
    Float64[]
end

"""
Update the values of the independent variational parameters
"""
function update_circuit_param!(circuit::VariationalQuantumCircuit, thetas::Vector{Float64})
    # update parameters
end

"""
Update a state using a circuit
"""
function update_quantum_state!(ucccirc::VariationalQuantumCircuit, state::QuantumState)
    # do something
end

"""
Wrap a QulacsParametricQuantumCircuit object, which will not be copied.
"""
struct QulacsVariationalQuantumCircuit <: VariationalQuantumCircuit
    qcircuit::QulacsParametricQuantumCircuit #qcircuit=>field, QulacsParametricQuantumCircuit::型(具体型で良い)
end

#=
abstract type Animal end
struct Goma <: Animal
      nakigoe::String
end

struct MyFamily
    members::Vector{Goma}
end

=#

"""
Copy a QulacsVariationalQuantumCircuit object.
This makes a copy of the underlying QulacsParametricQuantumCircuit object.
"""
function Base.copy(circuit::QulacsVariationalQuantumCircuit)
    QulacsVariationalQuantumCircuit(copy(circuit.qcircuit))
end

function num_theta(circuit::QulacsVariationalQuantumCircuit)
    get_parameter_count(circuit.qcircuit)
end

function get_thetas(circuit::QulacsVariationalQuantumCircuit)::Vector{Float64}
    return Float64[get_parameter(circuit.qcircuit, i)
        for i in 1:get_parameter_count(circuit.qcircuit)]
end

function update_circuit_param!(circuit::QulacsVariationalQuantumCircuit, thetas::Vector{Float64})
    for i in eachindex(thetas)
        set_parameter!(circuit.qcircuit, i, thetas[i])
    end
end

function update_quantum_state!(circuit::QulacsVariationalQuantumCircuit, state::QulacsQuantumState)
    update_quantum_state!(circuit.qcircuit, state)
end

function add_parametric_RY_gate!(circuit::QulacsVariationalQuantumCircuit, i::Int, angle::Float64)
    circuit.qcircuit.add_parametric_RY_gate(i-1, angle)
end

#function add_parametric_RZ_gate!(circuit::QulacsParametricQuantumCircuit, i::Int, angle::Float64)
#    circuit.pyobj.add_parametric_RZ_gate(i-1, angle)
#end

max_stars_repo_name: Sakurairihito/QCMaterialNew

max_issues_repo_name: sakurairihito/QCMaterialNew

max_forks_repo_name: sakurairihito/QCMaterialNew


text:
"Multiplication"
*(op::AbstractLinearTransform, x) = op.prod(x)

"Adjoint operator"
adjoint(op::AbstractLinearTransform) = LinearTransform(op.dimsRange, op.dimsDomain, op.tprod, op.prod)

"Unary addition"
+(op::AbstractLinearTransform) = op

"Unary negation"
function -(op::AbstractLinearTransform)
    prod = v -> -op.prod(v)
    tprod = w -> -op.tprod(w)
    LinearTransform(op.dimsRange, op.dimsDomain, prod, tprod)
end


text:
"""
Compare if two `TD`'s  are equal
"""
function Base.isequal(dat1::TD, dat2::TD)
	fnames = fieldnames(TD)
	vec=[(isequal(getfield(dat1, name),getfield(dat2, name))) for name in fnames]
	return all(vec)
end

"""
Return if two `TD`'s have same dimensions and bounds.
"""
function Base.isapprox(dat1::TD, dat2::TD)
	vec=([
		isequal(dat1.tgrid, dat2.tgrid),
		isequal(dat1.fields, dat2.fields),
		isequal(dat1.acqgeom, dat2.acqgeom, :receivers), # only receivers have to be the same
       		(size(dat1.d)==size(dat2.d)),
		])
	vec2=[size(dat1.d[iss,ifield])==size(dat2.d[iss,ifield]) for iss=1:dat1.acqgeom.nss, ifield=1:length(dat1.fields)]
	return (all(vec) & all(vec2))
end

function Base.length(data::TD)
	nr = data.acqgeom.nr;	nss = data.acqgeom.nss;	nt = length(data.tgrid);
	l=0
	for ifield = 1:length(data.fields), iss = 1:nss
		for ir = 1:nr[iss]
			for it in 1:nt
				l+=1
			end
		end
	end
	return l
end

"""
Return a vec of data object sorted in the order
time, receivers, supersource, fields
"""
function Base.vec(data::TD)
	nr = data.acqgeom.nr;	nss = data.acqgeom.nss;	nt = length(data.tgrid);
	v=Vector{Float64}()
	for ifield = 1:length(data.fields), iss = 1:nss
		dd=data.d[iss, ifield]
		append!(v,dd)
	end
	return v
end

"""
Method to Devectorize data!
No memory allocations
"""
function Base.copyto!(dataout::TD, v::AbstractVector{Float64})
	nr = dataout.acqgeom.nr;	nss = dataout.acqgeom.nss;	nt = length(dataout.tgrid);
	dout=getfield(dataout, :d)
	i0=0
	for iss=1:nss, ifield=1:length(dataout.fields)
		ddout=dout[iss,ifield]
		for ir = 1:nr[iss]
			for it in 1:nt
				ddout[it,ir]=v[i0+it]
			end
			i0+=nt
		end
	end
	return dataout
end

"""
Method to vectorize data!
No memory allocations
"""
function Base.copyto!(v::AbstractVector{Float64}, data::TD)
	nr = data.acqgeom.nr;	nss = data.acqgeom.nss;	nt = length(data.tgrid);
	d=getfield(data, :d)
	i0=0
	for iss=1:nss, ifield=1:length(data.fields)
		dd=d[iss,ifield]
		for ir = 1:nr[iss]
			for it in 1:nt
				v[i0+it]=dd[it,ir]
			end
			i0+=nt
		end
	end
	return v
end

"""
Fill with randn values
"""
function Random.randn!(data::TD)
	nr = data.acqgeom.nr;	nss = data.acqgeom.nss;	nt = length(data.tgrid);
	for ifield = 1:length(data.fields), iss = 1:nss
		dd=data.d[iss, ifield]
		Random.randn!(dd)
	end
	return data
end

"""
Copy `TD`'s, which are similar.
"""
function Base.copyto!(dataout::TD, data::TD)
	if(isapprox(dataout, data))
		dout=getfield(dataout, :d)
		din=getfield(data, :d)
		for iss=1:data.acqgeom.nss, ifield=1:length(data.fields)
			ddout=dout[iss,ifield]
			ddin=din[iss,ifield]
			copyto!(ddout,ddin)
		end
		return dataout
	else
		error("attempt to copy dissimilar data")
	end
end

"""
Returns bool depending on if input `data::TD` has all zeros or not.
"""
function Base.iszero(data::TD)
	return maximum(broadcast(maximum,abs,data.d)) == 0.0 ? true : false
end

"""
Returns dot product of data.

# Arguments

* `data1::TD` : data 1
* `data2::TD` : data 2

# Return

* dot product as `Float64`
"""
function LinearAlgebra.dot(data1::TD, data2::TD)
	if(isapprox(data1, data2))
		dotd = 0.0;
		for ifield = 1:length(data1.fields), iss = 1:data1.acqgeom.nss
			for ir = 1:data1.acqgeom.nr[iss], it = 1:length(data1.tgrid)
				dotd += data1.d[iss, ifield][it, ir] * data2.d[iss, ifield][it, ir]
			end
		end
		return dotd
	else
		error("cannot dot dissimilar datasets")
	end
end


text:
# Bogumił Kamiński, 2021

# Codes for chapter 1

# Code from section 1.2.1

function sum_n(n)
    s = 0
    for i in 1:n
        s += i
    end
    return s
end

@time sum_n(1_000_000_000)

# Code from section 1.4

@time using Plots
@time plot(1:10)
@time plot(1:10)

# Code allowing to reproduce the data frame presented in section 1.3

using DataFrames

DataFrame(a=1:3, name=["Alice", "Bob", "Clyde"],
          age=[19, 24, 21], friends=[[2], [1, 3], [2]],
          location=[(city="Atlanta", state="GA"),
                    (city="Boston", state="MA"),
                    (city="Austin", state="TX")])

# Code for comparison of languages in figure 1.1

# Source of data:
#   * https://benchmarksgame-team.pages.debian.net/benchmarksgame/fastest/julia-gcc.html
#   * https://benchmarksgame-team.pages.debian.net/benchmarksgame/fastest/julia-python3.html
#   * https://benchmarksgame-team.pages.debian.net/benchmarksgame/fastest/julia.html
# Download date:
#   29.01.2022
# Programming language versions used:
#   * Julia:    julia version 1.7.1
#   * Python 3: Python 3.10.1
#   * C:        gcc (Ubuntu 11.2.0-7ubuntu2) 11.2.0
#   * Java      openjdk 17.0.1 2021-10-19
#               OpenJDK Runtime Environment (build 17.0.1+12-39)
#               OpenJDK 64-Bit Server VM (build 17.0.1+12-39, mixed mode, sharing)

data = """
problem,language,time,size
n-body,C,2.13,1633
mandelbrot,C,1.3,1135
spectral norm,C,0.41,1197
fannkuch-redux,C,7.58,910
fasta,C,0.78,1463
k-nucleotide,C,3.96,1506
binary-trees,C,1.58,809
reverse-complement,C,0.41,1965
pidigits,C,0.56,1090
regex-redux,C,0.8,1397
n-body,Java,6.77,1489
mandelbrot,Java,4.1,796
spectral norm,Java,1.55,756
fannkuch-redux,Java,10.48,1282
fasta,Java,1.2,2543
k-nucleotide,Java,4.83,1812
binary-trees,Java,2.51,835
reverse-complement,Java,1.57,2183
pidigits,Java,0.79,764
regex-redux,Java,5.34,929
n-body,Python,541.34,1196
mandelbrot,Python,177.35,688
spectral norm,Python,112.97,407
fannkuch-redux,Python,341.45,950
fasta,Python,36.9,1947
k-nucleotide,Python,46.31,1967
binary-trees,Python,44.7,660
reverse-complement,Python,6.62,814
pidigits,Python,1.16,567
regex-redux,Python,1.34,1403
n-body,Julia,4.21,1111
mandelbrot,Julia,1.42,619
spectral norm,Julia,1.11,429
fannkuch-redux,Julia,7.83,1067
fasta,Julia,1.13,1082
k-nucleotide,Julia,4.94,951
binary-trees,Julia,7.28,634
reverse-complement,Julia,1.44,522
pidigits,Julia,0.97,506
regex-redux,Julia,1.74,759
"""

using CSV
using DataFrames
using Plots

df = CSV.read(IOBuffer(data), DataFrame)

plot(map([:time, :size],
         ["execution time (relative to C)",
          "code size (relative to C)"]) do col, title
    df_plot = unstack(df, :problem, :language, col)
    df_plot[!, Not(:problem)] ./= df_plot.C
    select!(df_plot, Not(:C))
    scatter(df_plot.problem, Matrix(select(df_plot, Not(:problem)));
            labels=permutedims(names(df_plot, Not(:problem))),
            ylabel=title,
            yaxis = col == :time ? :log : :none,
            xrotation=20,
            markershape=[:rect :diamond :circle],
            markersize=[4 5 5],
            markercolor=[:lightgray :lightgray :gold],
            xtickfontsize=7, ytickfontsize=7,
            legendfontsize=7, ylabelfontsize=7)
            hline!([1.0]; color="orange", labels="C")
end...)

max_stars_repo_name: Mo-Gul/JuliaForDataAnalysis

max_issues_repo_name: Mo-Gul/JuliaForDataAnalysis

max_forks_repo_name: Mo-Gul/JuliaForDataAnalysis


text:
module LogisticRegression2
using Revise
using Flux
using Flux: binarycrossentropy, throttle, @epochs
using Flux: Descent, TrackedArray, param, data
using Distributions
using LinearAlgebra
using Random
using Plots
using Base.Iterators: repeated

theme(:juno)

# Generate data
D = 2;   # Data dimension
n = 100; # Number of observations
x1 = rand(MvNormal([1., -5.], I), n);
x2 = rand(MvNormal([-2., -1.5], I), n);
y1 = repeat([1], n);
y2 = repeat([0], n);

classOneColor = :orange
classTwoColor = RGBA(0,.2,.7,1);
p = scatter(x1[1, :], x1[2, :], label="Class 1", color=classOneColor)
scatter!(p, x2[1, :], x2[2, :], label="Class 2", color=classTwoColor)
X = transpose([x1 x2]);
y = [y1; y2];

# Create model
struct LogisticRegression
    α :: TrackedArray # bias
    β :: TrackedArray # weights
    LogisticRegression(α :: Array, β :: Array) = new(param(α), param(β))
end

function (m :: LogisticRegression)(X)
    σ.(m.α .+ X*m.β)
end

function plot_decision_surface(m, xlim :: Tuple{Int, Int}, ylim :: Tuple{Int, Int})
    x = xlim[1]:0.5:xlim[2]
    y = ylim[1]:0.5:ylim[2]
    f = (x, y) -> (first ∘ data ∘ m)([x y])
    contour(x, y, f, fill = true, color = :pu_or)
end

function plot_fit(m)
    xlim = ylim = (-8, 8)
    p = plot_decision_surface(m, xlim, ylim)
    scatter!(p, x1[1, :], x1[2, :], label="Class 1", color=classOneColor)
    scatter!(p, x2[1, :], x2[2, :], label="Class 2", color=classTwoColor)
    p
end

function plot_loss(losses)
    plot(1:length(losses), losses)
end

function create_animation(αs :: Array, βs :: Array)
    anim = @animate for (α, β) in zip(αs, βs)
        m = LogisticRegression(α, β)
        display(plot_fit(m))
    end

    gif(anim, "plots/transition.gif", fps = 60)
end

m = LogisticRegression(rand(1), rand(D, 1)); display(plot_fit(m));
θ = Flux.params(m.α, m.β);
numEpochs = 100;
loss = (x, y) -> mean(binarycrossentropy.(m(x), y));

losses, αs, βs = [], [], []
cb = () -> begin
    l = Flux.data(loss(X, y));
    push!(losses, l);
    push!(αs, Flux.data(m.α));
    push!(βs, Flux.data(m.β));
end

dataset = repeated((X, y), numEpochs)
Flux.train!(loss, θ, dataset, Descent(0.1), cb = cb)

plot_loss(losses)
plot_fit(m)

create_animation(αs, βs)

end

max_stars_repo_name: SebastianCallh/julia-baby-steps

max_issues_repo_name: SebastianCallh/julia-baby-steps

max_forks_repo_name: SebastianCallh/julia-baby-steps


text:
using Behavior
using Niancat
using Test

#
# TestingGameService is a mock interface between
# a game and the Niancat service.
#
using Niancat.Gameface
using Niancat.Users
using Niancat.Scores
import Niancat.Gameface: score!, notify!, event!

struct TestingGameService <: GameService
    scores::Vector{Tuple{User, Score}}
    notifications::Vector{GameNotification}
    userevents::Vector{GameUserEvent}
    gameevents::Vector{GameEvent}
    instanceid::Int
    dictionaryrequests::Vector{String}

    TestingGameService(instanceid::Int = 1) = new([], [], [], [], instanceid, [])
end

Gameface.score!(tgs::TestingGameService, user::User, score::Score) = push!(tgs.scores, (user, score))
Gameface.notify!(tgs::TestingGameService, notification::GameNotification) = push!(tgs.notifications, notification)
Gameface.event!(tgs::TestingGameService, event::GameUserEvent) = push!(tgs.userevents, event)
Gameface.event!(tgs::TestingGameService, event::GameEvent) = push!(tgs.gameevents, event)
Gameface.gameinstanceid(tgs::TestingGameService) = tgs.instanceid
function Gameface.finddictionary(tgs::TestingGameService, dictionaryid::AbstractString) :: LanguageDictionary
    push!(tgs.dictionaryrequests, dictionaryid)
    SwedishDictionary(["DATORSPEL", "LEDARPOST", "ORDPUSSEL", "PUSSGURKA"])
end

lastnotification(tgs::TestingGameService) = tgs.notifications[end]
lastuserevent(tgs::TestingGameService) = tgs.userevents[end]
lastgameevent(tgs::TestingGameService) = tgs.gameevents[end]

#
# Include test files
#

include("games/niancat/runtests.jl")

@test runspec(pkgdir(Niancat), tags = "not @wip", presenter = TerseRealTimePresenter(ColorConsolePresenter()))


text:
include("kernels.jl")
include("gaussian_process_regressor.jl")

max_stars_repo_name: soblin/prml-julia

max_issues_repo_name: soblin/prml-julia

max_forks_repo_name: soblin/prml-julia


text:
## ==================================================================================================================='
# Should be made Common CSM functions
## ==================================================================================================================='

"""
    $SIGNATURES

Notify possible parent if clique is upsolved and exit the state machine.

Notes
- State machine function nr.0
- can recycle if two checks:
  - previous clique was identically downsolved
  - all children are also :uprecycled
"""
function canCliqMargRecycle_StateMachine(csmc::CliqStateMachineContainer)
  # @show getCliqFrontalVarIds(csmc.oldcliqdata), getCliqueStatus(csmc.oldcliqdata)
  infocsm(csmc, "0., $(csmc.incremental) ? :uprecycled => getCliqueStatus(csmc.oldcliqdata)=$(getCliqueStatus(csmc.oldcliqdata))")

  if areCliqVariablesAllMarginalized(csmc.dfg, csmc.cliq)

    # no work required other than assembling upward message
    if getSolverParams(csmc.cliqSubFg).dbg
      tmnow = now()
      tmpst = getCliqueStatus(csmc.cliq)
      @async begin
        mkpath(joinLogPath(csmc.cliqSubFg,"logs","cliq$(csmc.cliq.index)"))
        open(joinLogPath(csmc.cliqSubFg,"logs","cliq$(csmc.cliq.index)","marginalization.log"), "w") do f
          println(f, tmnow, ", marginalized from previous status ", tmpst)
        end
      end
    end
    prepPutCliqueStatusMsgUp!(csmc, :marginalized, dfg=csmc.dfg)

    # set marginalized color
    setCliqueDrawColor!(csmc.cliq, "blue")

    # set flag, looks to be previously unused???
    getCliqueData(csmc.cliq).allmarginalized = true

    # FIXME divert to rapid CSM exit
	# GUESSING THIS THE RIGHT WAY go to 4
	# return canCliqMargSkipUpSolve_StateMachine
  end

  # go to 0c.
  return canCliqIncrRecycle_StateMachine
end

"""
    $SIGNATURES

Final determination on whether can promote clique to `:uprecycled`.

Notes
- State machine function nr.0b
- Assume children clique status is available
- Will return to regular init-solve if new information in children -- ie not uprecycle or marginalized
"""
function checkChildrenAllUpRecycled_StateMachine(csmc::CliqStateMachineContainer)
  count = Int[]
  chldr = getChildren(csmc.tree, csmc.cliq)
  for ch in chldr
    chst = getCliqueStatus(ch)
    if chst in [:uprecycled; :marginalized]
      push!(count, 1)
    end
  end
  infocsm(csmc, "0b, checkChildrenAllUpRecycled_StateMachine -- length(chldr)=$(length(chldr)), sum(count)=$(sum(count))")

  # all children can be used for uprecycled -- i.e. no children have new information
  if sum(count) == length(chldr)
    # set up msg and exit go to 1
    sdims = Dict{Symbol,Float64}()
    for varid in getCliqAllVarIds(csmc.cliq)
      sdims[varid] = 0.0
    end

    # NOTE busy consolidating #459
    updateCliqSolvableDims!(csmc.cliq, sdims, csmc.logger)
        # setCliqueStatus!(csmc.cliq, :uprecycled)
        # replacing similar functionality from CSM 1.
    if getSolverParams(csmc.cliqSubFg).dbg
      tmnow = now()
      tmpst = getCliqueStatus(csmc.cliq)
      @async begin
        mkpath(joinLogPath(csmc.cliqSubFg,"logs","cliq$(csmc.cliq.index)"))
        open(joinLogPath(csmc.cliqSubFg,"logs","cliq$(csmc.cliq.index)","incremental.log"), "w") do f
          println(f, tmnow, ", marginalized from previous status ", tmpst)
        end
      end
    end
    prepPutCliqueStatusMsgUp!(csmc, :uprecycled, dfg=csmc.dfg)
    setCliqueDrawColor!(csmc.cliq, "orange")
    #go to 10
    return canCliqDownSolve_StateMachine
        # # go to 1
        # return isCliqUpSolved_StateMachine
  end

  # return to regular solve, go to 2
  return buildCliqSubgraph_StateMachine
end

"""
    $SIGNATURES

Determine if clique is upsolved by incremental update and exit the state machine.

Notes
- State machine function nr.0c
- can recycle if two checks:
  - previous clique was identically downsolved
  - all children are also :uprecycled
"""
function canCliqIncrRecycle_StateMachine(csmc::CliqStateMachineContainer)
	# check if should be trying and can recycle clique computations
    if csmc.incremental && getCliqueStatus(csmc.oldcliqdata) == :downsolved
      csmc.cliq.data.isCliqReused = true
      # check if a subgraph will be needed later
      if csmc.dodownsolve
        # yes need subgraph and need more checks, so go to 2
        return buildCliqSubgraph_StateMachine
      else
        # one or two checks say yes, so go to 4
        return canCliqMargSkipUpSolve_StateMachine
      end
    end

    # nope, regular clique init-solve, go to 1
    return isCliqUpSolved_StateMachine
end

"""
    $SIGNATURES

Either construct and notify of a new upward initialization message and progress to downsolve checks,
or circle back and start building the local clique subgraph.

Notes
- State machine function nr.1
- Root clique message should be empty since it has an empty separator.
"""
function isCliqUpSolved_StateMachine(csmc::CliqStateMachineContainer)

  infocsm(csmc, "1, isCliqUpSolved_StateMachine")
  cliqst = getCliqueStatus(csmc.cliq)

  # if upward complete for any reason, prepare and send new upward message
  if cliqst in [:upsolved; :downsolved; :marginalized; :uprecycled]
    # construct init's up msg from initialized separator variables
    # NOTE cliqSubFg has not been copied yet
    prepPutCliqueStatusMsgUp!(csmc, cliqst, dfg=csmc.dfg)
    #go to 10
    return canCliqDownSolve_StateMachine
  end
  # go to 2
  return buildCliqSubgraph_StateMachine
end

"""
    $SIGNATURES

Build a sub factor graph for clique variables from the larger factor graph.

Notes
- State machine function nr.2
"""
function buildCliqSubgraph_StateMachine(csmc::CliqStateMachineContainer)
  # build a local subgraph for inference operations
  infocsm(csmc, "2, build subgraph syms=$(getCliqAllVarIds(csmc.cliq))")
  buildCliqSubgraph!(csmc.cliqSubFg, csmc.dfg, csmc.cliq)

  # if dfg, store the cliqSubFg for later debugging
  _dbgCSMSaveSubFG(csmc, "fg_build")

  # go to 4
  return canCliqMargSkipUpSolve_StateMachine
end

"""
$SIGNATURES

Quick redirection of out-marginalized cliques to downsolve path, or wait on children cliques to get a csm status.

Notes
- State machine function nr.4
"""
function canCliqMargSkipUpSolve_StateMachine(csmc::CliqStateMachineContainer)

  cliqst = getCliqueStatus(csmc.oldcliqdata)
  infocsm(csmc, "4, canCliqMargSkipUpSolve_StateMachine, $cliqst, csmc.incremental=$(csmc.incremental)")

  # if clique is out-marginalized, then no reason to continue with upsolve
  # marginalized state is set in `canCliqMargRecycle_StateMachine`
  if cliqst == :marginalized
    # go to 10 -- Add case for IIF issue #474
    return canCliqDownSolve_StateMachine
  end

  # go to 4e
  return blockUntilChildrenHaveStatus_StateMachine
end

## ==================================================================================================================='
## Does this have a place
## ==================================================================================================================='

"""
    $SIGNATURES

Build a sub factor graph for clique variables from the larger factor graph.

Notes
- State machine function nr.2r
"""
function buildCliqSubgraphForDown_StateMachine(csmc::CliqStateMachineContainer)
  # build a local subgraph for inference operations
  syms = getCliqAllVarIds(csmc.cliq)
  infocsm(csmc, "2r, build subgraph syms=$(syms)")
  csmc.cliqSubFg = buildSubgraph(csmc.dfg, syms, 1; verbose=false)

  opts = getSolverParams(csmc.dfg)
  # store the cliqSubFg for later debugging
  _dbgCSMSaveSubFG(csmc, "fg_build_down")

  # go to 10
  return canCliqDownSolve_StateMachine
end

## ==================================================================================================================='
## Split and use
## ==================================================================================================================='

"""
    $SIGNATURES

One of the last steps in CSM to clean up after a down solve.

Notes
- CSM function 11b.
"""
function cleanupAfterDownSolve_StateMachine(csmc::CliqStateMachineContainer)
  # RECENT split from 11 (using #760 solution for deleteMsgFactors)
  opts = getSolverParams(csmc.cliqSubFg)

  # set PPE and solved for all frontals
  for sym in getCliqFrontalVarIds(csmc.cliq)
    # set PPE in cliqSubFg
    setVariablePosteriorEstimates!(csmc.cliqSubFg, sym)
    # set solved flag
    vari = getVariable(csmc.cliqSubFg, sym)
    setSolvedCount!(vari, getSolvedCount(vari, :default)+1, :default )
  end

  # store the cliqSubFg for later debugging
  _dbgCSMSaveSubFG(csmc, "fg_afterdownsolve")

  # transfer results to main factor graph
  frsyms = getCliqFrontalVarIds(csmc.cliq)
  infocsm(csmc, "11, finishingCliq -- going for transferUpdateSubGraph! on $frsyms")
  transferUpdateSubGraph!(csmc.dfg, csmc.cliqSubFg, frsyms, csmc.logger, updatePPE=true)

  infocsm(csmc, "11, doCliqDownSolve_StateMachine -- before prepPutCliqueStatusMsgDwn!")
  cliqst = prepPutCliqueStatusMsgDwn!(csmc, :downsolved)
  infocsm(csmc, "11, doCliqDownSolve_StateMachine -- just notified prepPutCliqueStatusMsgDwn!")

  # remove msg factors that were added to the subfg
  rmFcts = deleteMsgFactors!(csmc.cliqSubFg)
  infocsm(csmc, "11, doCliqDownSolve_StateMachine -- removing all up/dwn message factors, length=$(length(rmFcts))")

  infocsm(csmc, "11, doCliqDownSolve_StateMachine -- finished, exiting CSM on clique=$(csmc.cliq.index)")
  # and finished
  return IncrementalInference.exitStateMachine
end

"""
$SIGNATURES

Root clique upsolve and downsolve are equivalent, so skip a repeat downsolve, just set messages and just exit directly.

Notes
- State machine function nr. 10b
- Separate out during #459 dwnMsg consolidation.

DevNotes
- TODO should this consolidate some work with 11b?
"""
function specialCaseRootDownSolve_StateMachine(csmc::CliqStateMachineContainer)
  # this is the root clique, so assume already downsolved -- only special case
  dwnmsgs = getCliqDownMsgsAfterDownSolve(csmc.cliqSubFg, csmc.cliq)
  setCliqueDrawColor!(csmc.cliq, "lightblue")

  # this part looks like a pull model
  # JT 459 putMsgDwnThis!(csmc.cliq, dwnmsgs)
  putDwnMsgConsolidated!(csmc.cliq.data, dwnmsgs) # , from=:putMsgDwnThis!  putCliqueMsgDown!
  setCliqueStatus!(csmc.cliq, :downsolved)
  csmc.dodownsolve = false

  # Update estimates and transfer back to the graph
  frsyms = getCliqFrontalVarIds(csmc.cliq)
  # set PPE and solved for all frontals
  for sym in frsyms
    # set PPE in cliqSubFg
    setVariablePosteriorEstimates!(csmc.cliqSubFg, sym)
    # set solved flag
    vari = getVariable(csmc.cliqSubFg, sym)
    setSolvedCount!(vari, getSolvedCount(vari, :default)+1, :default )
  end
  # Transfer to parent graph
  transferUpdateSubGraph!(csmc.dfg, csmc.cliqSubFg, frsyms, updatePPE=true)

  prepPutCliqueStatusMsgDwn!(csmc, :downsolved)
  # notifyCliqDownInitStatus!(csmc.cliq, :downsolved, logger=csmc.logger)

  # bye
  return IncrementalInference.exitStateMachine
end

## ==================================================================================================================='
## Can be consolidated/used (mostly used already as copies in X)
## ==================================================================================================================='

"""
    $SIGNATURES

Do cliq downward inference

Notes:
- State machine function nr. 11
"""
function doCliqDownSolve_StateMachine(csmc::CliqStateMachineContainer)
  infocsm(csmc, "11, doCliqDownSolve_StateMachine")
  setCliqueDrawColor!(csmc.cliq, "red")
  # get down msg from parent (assuming root clique CSM wont make it here)
  # this looks like a pull model #674
  prnt = getParent(csmc.tree, csmc.cliq)
  dwnmsgs = fetchDwnMsgConsolidated(prnt[1])
  infocsm(csmc, "11, doCliqDownSolve_StateMachine -- dwnmsgs=$(collect(keys(dwnmsgs.belief)))")

  __doCliqDownSolve!(csmc, dwnmsgs)

  # compute new down messages
  infocsm(csmc, "11, doCliqDownSolve_StateMachine -- going to set new down msgs.")
  newDwnMsgs = getSetDownMessagesComplete!(csmc.cliqSubFg, csmc.cliq, dwnmsgs, csmc.logger)

  prepPutCliqueStatusMsgDwn!(csmc, :downsolved)

  # update clique subgraph with new status
  setCliqueDrawColor!(csmc.cliq, "lightblue")

  infocsm(csmc, "11, doCliqDownSolve_StateMachine -- finished with downGibbsCliqueDensity, now update csmc")

  # go to 11b.
  return cleanupAfterDownSolve_StateMachine
end

# XXX only use skip down part
"""
    $SIGNATURES

Direct state machine to continue with downward solve or exit.

Notes
- State machine function nr. 10
"""
function canCliqDownSolve_StateMachine(csmc::CliqStateMachineContainer)
  infocsm(csmc, "10, canCliqDownSolve_StateMachine, csmc.dodownsolve=$(csmc.dodownsolve).")
  # finished and exit downsolve
  if !csmc.dodownsolve
    infocsm(csmc, "10, canCliqDownSolve_StateMachine -- shortcut exit since downsolve not required.")
    return IncrementalInference.exitStateMachine
  end

  # assume separate down solve via solveCliq! call, but need a csmc.cliqSubFg
  # could be dedicated downsolve that was skipped during previous upsolve only call
  # e.g. federated solving case (or debug)
  if length(ls(csmc.cliqSubFg)) == 0
    # first need to fetch cliq sub graph
    infocsm(csmc, "10, canCliqDownSolve_StateMachine, oops no cliqSubFg detected, lets go fetch a copy first.")

    # go to 2b
    return buildCliqSubgraphForDown_StateMachine
  end

  # both parent or otherwise might start by immediately doing downsolve, so likely need cliqSubFg in both cases
  # e.g. federated solving case (or debug)
  prnt = getParent(csmc.tree, csmc.cliq)
  if 0 == length(prnt) # check if have parent
    # go to 10b
    return specialCaseRootDownSolve_StateMachine
  end

  # go to 8c
  return waitChangeOnParentCondition_StateMachine
  # # go to 10a
  # return wipRedirect459Dwn_StateMachine
end

# XXX does not look like it has a place
"""
    $SIGNATURES

Is upsolve complete or should the CSM solving process be repeated.

Notes
- State machine function nr.9

DevNotes
- FIXME FIXME FIXME ensure init worked
"""
function checkUpsolveFinished_StateMachine(csmc::CliqStateMachineContainer)
  cliqst = getCliqueStatus(csmc.cliq)
  infocsm(csmc, "9, checkUpsolveFinished_StateMachine")
  if cliqst == :upsolved
      frsyms = getCliqFrontalVarIds(csmc.cliq)
    infocsm(csmc, "9, checkUpsolveFinished_StateMachine -- going for transferUpdateSubGraph! on $frsyms")
    # TODO what about down solve??
    transferUpdateSubGraph!(csmc.dfg, csmc.cliqSubFg, frsyms, csmc.logger, updatePPE=false)

    # remove any solvable upward cached data -- TODO will have to be changed for long down partial chains
    # assuming maximally complte up solved cliq at this point
    # lockUpStatus!(csmc.cliq, csmc.cliq.index, true, csmc.logger, true, "9.finishCliqSolveCheck")
    sdims = Dict{Symbol,Float64}()
    for varid in getCliqAllVarIds(csmc.cliq)
      sdims[varid] = 0.0
    end
    updateCliqSolvableDims!(csmc.cliq, sdims, csmc.logger)
    # unlockUpStatus!(csmc.cliq)

    # go to 10
    return canCliqDownSolve_StateMachine # IncrementalInference.exitStateMachine
  elseif cliqst == :initialized
    # setCliqueDrawColor!(csmc.cliq, "sienna")

    # go to 7
    return determineCliqNeedDownMsg_StateMachine
  else
    infocsm(csmc, "9, checkUpsolveFinished_StateMachine -- init not complete and should wait on init down message.")
    # setCliqueDrawColor!(csmc.cliq, "coral")
    # TODO, potential problem with trying to downsolve
    # return canCliqMargSkipUpSolve_StateMachine
  end

  # go to 4b (redirected here during #459 dwnMsg effort)
  return trafficRedirectConsolidate459_StateMachine
    # # go to 4
    # return canCliqMargSkipUpSolve_StateMachine # whileCliqNotSolved_StateMachine
end

# XXX does not look like it has a place
"""
    $SIGNATURES

Do up initialization calculations,  preparation for solving Chapman-Kolmogorov
transit integral in upward direction.

Notes
- State machine function nr. 8f
- Includes initialization routines.
- Adds `:__LIKELIHOODMESSAGE__` factors but does not remove.
- gets msg likelihoods from cliqSubFg, see #760

DevNotes
- TODO: Make multi-core
"""
function prepInitUp_StateMachine(csmc::CliqStateMachineContainer)
  setCliqueDrawColor!(csmc.cliq, "green")

  # check if init is required and possible
  infocsm(csmc, "8f, prepInitUp_StateMachine -- going for doCliqAutoInitUpPart1!.")
  # get incoming clique up messages
  upmsgs = getMsgsUpInitChildren(csmc, skip=[csmc.cliq.index;])
  # Filter for usable messages
  ## FIXME joint decomposition as differential likelihoods conversion must still be done for init
  dellist = []
  for (chid, lm) in upmsgs
    if !(lm.status in [:initialized;:upsolved;:marginalized;:downsolved;:uprecycled])
      push!(dellist, chid)
    end
  end
  dellist .|> x->delete!(upmsgs, x)

  # remove all lingering upmessage likelihoods
  oldTags = deleteMsgFactors!(csmc.cliqSubFg)
  0 < length(oldTags) ? @warn("stale LIKELIHOODMESSAGE tags present in prepInitUp_StateMachine") : nothing

  # add incoming up messages as priors to subfg
  infocsm(csmc, "8f, prepInitUp_StateMachine -- adding up message factors")
  # interally adds :__LIKELIHOODMESSAGE__, :__UPWARD_DIFFERENTIAL__, :__UPWARD_COMMON__ to each of the factors
  msgfcts = addMsgFactors!(csmc.cliqSubFg, upmsgs, UpwardPass)

  # store the cliqSubFg for later debugging
  _dbgCSMSaveSubFG(csmc, "fg_beforeupsolve")

  # go to 8m
  return tryUpInitCliq_StateMachine
end

#XXX uses same functions, can be split at message
"""
    $SIGNATURES

Calculate the full upward Chapman-Kolmogorov transit integral solution approximation (i.e. upsolve).

Notes
- State machine function nr. 8g
- Assumes LIKELIHOODMESSAGE factors are in csmc.cliqSubFg but does not remove them.
- TODO: Make multi-core

DevNotes
- NEEDS DFG v0.8.1, see IIF #760
"""
function doCliqUpSolveInitialized_StateMachine(csmc::CliqStateMachineContainer)

  __doCliqUpSolveInitialized!(csmc)

  # Send upward message, NOTE consolidation WIP #459
  infocsm(csmc, "8g, doCliqUpSolveInitialized_StateMachine -- setting up messages with status = :upsolved")
  prepPutCliqueStatusMsgUp!(csmc, :upsolved)

  # go to 8h
  return rmUpLikeliSaveSubFg_StateMachine
end

# XXX can be changed to utility function
"""
    $SIGNATURES

Close out up solve attempt by removing any LIKELIHOODMESSAGE and save a debug cliqSubFg.

Notes
- State machine function nr. 8h
- Assumes LIKELIHOODMESSAGE factors are in csmc.cliqSubFg and also removes them.
- TODO: Make multi-core

DevNotes
- NEEDS DFG v0.8.1, see IIF #760
"""
function rmUpLikeliSaveSubFg_StateMachine(csmc::CliqStateMachineContainer)
  #
  status = getCliqueStatus(csmc.cliq)
  opts = getSolverParams(csmc.dfg)

  # remove msg factors that were added to the subfg
  tags__ = opts.useMsgLikelihoods ? [:__UPWARD_COMMON__;] : [:__LIKELIHOODMESSAGE__;]
  msgfcts = deleteMsgFactors!(csmc.cliqSubFg, tags__)
  infocsm(csmc, "8g, doCliqUpsSolveInit.! -- status = $(status), removing $(tags__) factors, length=$(length(msgfcts))")

  # store the cliqSubFg for later debugging
  _dbgCSMSaveSubFG(csmc, "fg_afterupsolve")

  # go to 9
  return checkUpsolveFinished_StateMachine
end

# XXX does not look like it has a place, replaced by take!
"""
    $SIGNATURES

When this clique needs information from parent to continue but parent is still busy.
Reasons are either downsolve or need down init message information (which are similar).

Notes
- State machine function nr. 8c
- bad idea to injectDelayBefore this function, because it will delay waiting on the parent past the event.
"""
function waitChangeOnParentCondition_StateMachine(csmc::CliqStateMachineContainer)
  #
  # setCliqueDrawColor!(csmc.cliq, "coral")

  prnt = getParent(csmc.tree, csmc.cliq)
  if 0 < length(prnt)
    infocsm(csmc, "8c, waitChangeOnParentCondition_StateMachine, wait on parent=$(prnt[1].index) for condition notify.")

    prntst = fetchDwnMsgConsolidated(prnt[1]).status
    if prntst != :downsolved
      wait(getSolveCondition(prnt[1]))
    end

    # wait(getSolveCondition(prnt[1]))
  else
    infocsm(csmc, "8c, waitChangeOnParentCondition_StateMachine, cannot wait on parent for condition notify.")
    @warn "no parent!"
  end

  # Listing most likely status values that might lead to 4b (TODO needs to be validated)
  if getCliqueStatus(csmc.cliq) in [:needdownmsg; :initialized; :null]
    # go to 4b
    return trafficRedirectConsolidate459_StateMachine
  end

  ## consolidation from CSM 10a
  # yes, continue with downsolve
  infocsm(csmc, "10a, wipRedirect459Dwn_StateMachine, parent status=$prntst.")
  if prntst != :downsolved
    infocsm(csmc, "10a, wipRedirect459Dwn_StateMachine, going around again.")
    return canCliqDownSolve_StateMachine
  end

  infocsm(csmc, "10a, wipRedirect459Dwn_StateMachine, going for down solve.")
  # go to 11
  return doCliqDownSolve_StateMachine
end

# XXX functions replaced in take! structure
"""
    $SIGNATURES

Do down solve calculations, loosely translates to solving Chapman-Kolmogorov
transit integral in downward direction.

Notes
- State machine function nr. 8e.ii.
  - Follows routines in 8c.
- Pretty major repeat of functionality, FIXME
- TODO: Make multi-core

DevNotes
- TODO Lots of cleanup required, especially from calling function.
- TODO move directly into a CSM state function
- TODO figure out the diference, redirection, ie relation to 8m??

CONSOLIDATED OLDER FUNCTION DOCS:
Initialization requires down message passing of more specialized down init msgs.
This function performs any possible initialization of variables and retriggers
children cliques that have not yet initialized.

Notes:
- Assumed this function is only called after status from child clique up inits completed.
- Assumes cliq has parent.
  - will fetch message from parent
- Will perform down initialization if status == `:needdownmsg`.
- might be necessary to pass furhter down messges to child cliques that also `:needdownmsg`.
- Will not complete cliq solve unless all children are `:upsolved` (upward is priority).
- `dwinmsgs` assumed to come from parent initialization process.
- assume `subfg` as a subgraph that can be modified by this function (add message factors)
  - should remove message prior factors from subgraph before returning.
- May modify `cliq` values.
  - `putMsgUpInit!(cliq, msg)`
  - `setCliqueStatus!(cliq, status)`
  - `setCliqueDrawColor!(cliq, "sienna")`
  - `notifyCliqDownInitStatus!(cliq, status)`

Algorithm:
- determine which downward messages influence initialization order
- initialize from singletons to most connected non-singletons
- revert back to needdownmsg if cycleInit does nothing
- can only ever return :initialized or :needdownmsg status
"""
function tryDwnInitCliq_StateMachine(csmc::CliqStateMachineContainer)
  setCliqueDrawColor!(csmc.cliq, "green")
  opt = getSolverParams(csmc.cliqSubFg)
  dwnkeys_ = lsf(csmc.cliqSubFg, tags=[:__DOWNWARD_COMMON__;]) .|> x->ls(csmc.cliqSubFg, x)[1]

  ## TODO deal with partial inits only, either delay or continue at end...
  # find intersect between downinitmsgs and local clique variables
  # if only partials available, then
  infocsm(csmc, "8e.ii, tryDwnInitCliq_StateMachine, do cliq init down dwinmsgs=$(dwnkeys_)")
  # get down variable initialization order
  initorder = getCliqInitVarOrderDown(csmc.cliqSubFg, csmc.cliq, dwnkeys_)
  infocsm(csmc, "8e.ii, tryDwnInitCliq_StateMachine,  initorder=$(initorder)")

  # store the cliqSubFg for later debugging
  if opt.dbg
    DFG.saveDFG(csmc.cliqSubFg, joinpath(opt.logpath,"logs/cliq$(csmc.cliq.index)/fg_beforedowninit"))
  end

  # cycle through vars and attempt init
  infocsm(csmc, "8e.ii, tryDwnInitCliq_StateMachine, cycle through vars and attempt init")
  # cliqst = :needdownmsg
  if cycleInitByVarOrder!(csmc.cliqSubFg, initorder, logger=csmc.logger)
    # cliqst = :initialized
    # TODO: transfer values changed in the cliques should be transfered to the tree in proc 1 here.
    # # TODO: is status of notify required here? either up or down msg??
    setCliqueDrawColor!(csmc.cliq, "sienna")
    setCliqueStatus!(csmc.cliq, :initialized)
  end

  # go to 8l
  return rmMsgLikelihoodsAfterDwn_StateMachine
end

# XXX functions replaced in take! structure
"""
    $SIGNATURES

Check if the clique is fully initialized.

Notes
- State machine function nr. 8m

DevNotes
- TODO figure out the relation or possible consolidation with 8e.ii
"""
function tryUpInitCliq_StateMachine(csmc::CliqStateMachineContainer)
  # attempt initialize if necessary
  setCliqueDrawColor!(csmc.cliq, "green")
  someInit = false
  if !areCliqVariablesAllInitialized(csmc.cliqSubFg, csmc.cliq)
    # structure for all up message densities computed during this initialization procedure.
    varorder = getCliqVarInitOrderUp(csmc.cliqSubFg)
    someInit = cycleInitByVarOrder!(csmc.cliqSubFg, varorder, logger=csmc.logger)
    # is clique fully upsolved or only partially?
    # print out the partial init status of all vars in clique
    printCliqInitPartialInfo(csmc.cliqSubFg, csmc.cliq, csmc.logger)
    infocsm(csmc, "8m, tryUpInitCliq_StateMachine -- someInit=$someInit, varorder=$varorder")
  end

  chldneed = doAnyChildrenNeedDwnMsg(getChildren(csmc.tree, csmc.cliq))
  allvarinit = areCliqVariablesAllInitialized(csmc.cliqSubFg, csmc.cliq)
  infocsm(csmc, "8m, tryUpInitCliq_StateMachine -- someInit=$someInit, chldneed=$chldneed, allvarinit=$allvarinit")

  upmessages = fetchMsgsUpChildrenDict(csmc)
  all_child_status = map(msg -> msg.status, values(upmessages))

  # redirect if any children needdownmsg
  if someInit || chldneed
    # Calculate and share the children sum solvableDim information for priority initialization
    totSolDims = Dict{Int, Float64}()
    for (clid, upmsg) in upmessages
      totSolDims[clid] = 0
      for (varsym, tbup) in upmsg.belief
        totSolDims[clid] += tbup.solvableDim
      end
    end
    infocsm(csmc, "8m, tryUpInitCliq_StateMachine -- totSolDims=$totSolDims")

    # prep and put down init message
    setCliqueDrawColor!(csmc.cliq, "sienna")
    prepPutCliqueStatusMsgDwn!(csmc, :initialized, childSolvDims=totSolDims)

    # go to 7e
    return slowWhileInit_StateMachine

    # (short cut) check again if all cliq vars have been initialized so that full inference can occur on clique
  # clique should be initialized and all children upsolved, uprecycled, or marginalized
  elseif allvarinit && all(in.(all_child_status, Ref([:upsolved; :uprecycled; :marginalized])))
    infocsm(csmc, "8m, tryUpInitCliq_StateMachine -- all initialized")
    setCliqueDrawColor!(csmc.cliq, "sienna")
    # don't send a message yet since the upsolve is about to occur too
    setCliqueStatus!(csmc.cliq, :initialized)

    # go to 8g.
    return doCliqUpSolveInitialized_StateMachine
  end

  infocsm(csmc, "8m, tryUpInitCliq_StateMachine -- not able to init all")
  # TODO Simplify this
  status = getCliqueStatus(csmc.cliq)

  if !(status == :initialized || length(getParent(csmc.tree, csmc.cliq)) == 0)
    # notify of results (big part of #459 consolidation effort)
    setCliqueDrawColor!(csmc.cliq, "orchid")
    prepPutCliqueStatusMsgUp!(csmc, :needdownmsg)
  end

  # go to 8h
  return rmUpLikeliSaveSubFg_StateMachine
end

# XXX maybe change to utility function
"""
    $SIGNATURES

Remove any `:__LIKELIHOODMESSAGE__` from `cliqSubFg`.

Notes
- State machine function nr.8l
"""
function rmMsgLikelihoodsAfterDwn_StateMachine(csmc::CliqStateMachineContainer)
  ## TODO only remove :__DOWNWARD_COMMON__ messages here
  #

  _dbgCSMSaveSubFG(csmc, "fg_afterdowninit")

  ## FIXME move this to separate state in CSM.
  # remove all message factors
  # remove msg factors previously added
  fctstorm = deleteMsgFactors!(csmc.cliqSubFg)
  infocsm(csmc, "8e.ii., tryDwnInitCliq_StateMachine, removing factors $fctstorm")

  # go to 8d
  return decideUpMsgOrInit_StateMachine
end