text:
"""
    $SIGNATURES

Draw and show the factor graph `<:AbstractDFG` via system graphviz and pdf app.

Notes
- Should not be calling outside programs.
- Need long term solution
- DFG's `toDotFile` a better solution -- view with `xdot` application.
- also try `engine={"sfdp","fdp","dot","twopi","circo","neato"}`

Future:
- Might be kept with different call strategy since this function so VERY useful!
- Major issue that this function calls an external program such as "evince", which should be
   under user control only.
- Maybe solution is
- ```toDot(fg,file=...); @async run(`xdot file.dot`)```, or
  - ```toDot(fg,file=...); exportPdf(...); @async run(`evince ...pdf`)```.
"""
function writeGraphPdf(fgl::G;
                       viewerapp::String="evince",
                       filepath::AS="/tmp/fg.pdf",
                       engine::AS="neato", #sfdp
                       show::Bool=true ) where {G <: AbstractDFG, AS <: AbstractString}
  #
  @warn "writeGraphPdf is function changing to drawGraph, see DFG.toDotFile(dfg) as part of the long term solution."

  fgd = fgl
  @info "Writing factor graph file"
  fext = split(filepath, '.')[end]
  fpwoext = filepath[1:(end-length(fext)-1)] # split(filepath, '.')[end-1]
  dotfile = fpwoext*".dot"

  # create the dot file
  DFG.toDotFile(fgl, dotfile)

  try
    run(`$(engine) $(dotfile) -T$(fext) -o $(filepath)`)
    show ? (@async run(`$(viewerapp) $(filepath)`)) : nothing
  catch e
    @warn "not able to show $(filepath) with viewerapp=$(viewerapp). Exception e=$(e)"
  end
  nothing
end

# explore all shortest paths combinations in verts, add neighbors and reference subgraph
# Using unique index into graph data structure
function subgraphFromVerts(fgl::FactorGraph,
    verts::Array{Int,1};
    neighbors::Int=0  )

  vertdict = Dict{Int,Graphs.ExVertex}()
  for vert in verts
    vertdict[vert] = fgl.g.vertices[vert]
  end

  return subgraphFromVerts(fgl,vertdict,neighbors=neighbors)
end

"""
    $(SIGNATURES)

Explore all shortest paths combinations in verts, add neighbors and reference
subgraph using unique index into graph data structure.
"""
function subGraphFromVerts(fgl::FactorGraph,
                           verts::T;
                           neighbors::Int=0  ) where {T <: Union{Vector{String},Vector{Symbol}}}
  #
  vertdict = Dict{Int,Graphs.ExVertex}()
  for vert in verts
    id = -1
    vsym = Symbol(vert)
    if haskey(fgl.IDs, vsym)
      id = fgl.IDs[vsym]
    else
      error("FactorGraph01 only supports variable node subgraph search at this point")
    end
    vertdict[id] = getVert(fgl, vsym) # fgl.g.vertices[id]
  end

  return subgraphFromVerts(fgl,vertdict,neighbors=neighbors)
end

function subgraphFromVerts(fgl::FactorGraph,
                           verts::T;
                           neighbors::Int=0  ) where {T <: Union{Vector{String},Vector{Symbol}, Dict{Int,Graphs.ExVertex}}}
  #
  @warn "`subgraphFromVerts` deprecated, use `subGraphFromVerts` instead."
  subGraphFromVerts(fgl, verts, neighbors=neighbors)
end

"""
    $SIGNATURES

Return array of all variable vertices in a clique.
"""
function getCliqVars(subfg::FactorGraph, cliq::Graphs.ExVertex)
  verts = Graphs.ExVertex[]
  for vid in getCliqVars(subfg, cliq)
    push!(verts, getVert(subfg, vid, api=localapi))
  end
  return verts
end

### TODO: TO BE REFACTORED FOR DFG

# some plotting functions on the factor graph
function stackVertXY(fg::FactorGraph, lbl::String)
    v = dlapi.getvertex(fg,lbl)
    vals = getVal(v)
    X=vec(vals[1,:])
    Y=vec(vals[2,:])
    return X,Y
end

### TODO: TO BE REFACTORED FOR DFG

"""
    $(SIGNATURES)

Save mostly complete Factor Graph type by converting complicated FunctionNodeData
types to 'Packed' types using user supplied converters. Ground truth can also be
saved and recovered by the associated loadjld(file="tempfg.jld2") method.

Notes:
- Must use `.jld2` since Julia 1.0 (previous version was deprecated).
"""
function savejld(fgl::G;
                 file::AbstractString="tempfg.jld2"  ) where G <: AbstractDFG
  #
  @error "savejld has been deprecated, use saveDFG instead"
  # fgs = encodefg(fgl)
  @save file fgl
  return file
end

"""
    $(SIGNATURES)

Opposite of savejld(fg, gt=gt, file="tempfg.jl") to load data from file. This function
uses the unpacking converters for converting all PackedInferenceType to FunctorInferenceType.
"""
function loadjld(;file::AbstractString="tempfg.jld2")
  @error "loadjld has been deprecated, use loadDFG instead"
  fgd = @load file fgl
  return fgd
end

"""
    $(SIGNATURES)

Return the last up message stored in `cliq` of Bayes (Junction) tree.
"""
function upMsg(cliq::Graphs.ExVertex)
  @warn "deprecated upMsg, use getUpMsg instead"
  getData(cliq).upMsg
end
function upMsg(btl::BayesTree, sym::Symbol)
  @warn "deprecated upMsg, use getUpMsg instead"
  upMsg(whichCliq(btl, sym))
end

"""
    $(SIGNATURES)

Return the last down message stored in `cliq` of Bayes (Junction) tree.
"""
function dwnMsg(cliq::Graphs.ExVertex)
  @warn "deprecated dwnMsg, use getDwnMsgs instead"
  getData(cliq).dwnMsg
end
function dwnMsg(btl::BayesTree, sym::Symbol)
  @warn "deprecated dwnMsg, use getDwnMsgs instead"
  dwnMsg(whichCliq(btl, sym))
end

function getCliquePotentials!(fg::FactorGraph, bt::BayesTree, chkcliq::Int)
    @error "getCliquePotentials! deprecated, use setCliqPotentials! with DFG objects instead of FactorGraph"
    setCliqPotentials!(fg, bt.cliques[chkcliq])
end

"""
    $(SIGNATURES)

Pass NBPMessages back down the tree -- pre order tree traversal.
"""
function downMsgPassingRecursive(inp::ExploreTreeType{T}; N::Int=100, dbg::Bool=false, drawpdf::Bool=false) where {T}
  @info "====================== Clique $(inp.cliq.attributes["label"]) ============================="

  mcmciter = inp.prnt != nothing ? 3 : 0; # skip mcmc in root on dwn pass
  rDDT = downGibbsCliqueDensity(inp.fg, inp.cliq, inp.sendmsgs, N, mcmciter, dbg) #dwnMsg
  updateFGBT!(inp.fg, inp.bt, inp.cliq.index, rDDT, dbg=dbg, fillcolor="lightblue")
  setCliqStatus!(inp.cliq, :downsolved)
  drawpdf ? drawTree(inp.bt) : nothing

  # rr = Array{Future,1}()
  # pcs = procs()

  ddt=nothing
  for child in out_neighbors(inp.cliq, inp.bt.bt)
    ett = ExploreTreeType(inp.fg, inp.bt, child, inp.cliq, [rDDT.dwnMsg])#inp.fg
    ddt = downMsgPassingRecursive( ett , N=N, dbg=dbg )
    drawpdf ? drawTree(inp.bt) : nothing
  end

  # return modifications to factorgraph to calling process
  return ddt
end

# post order tree traversal and build potential functions
function upMsgPassingRecursive(inp::ExploreTreeType{T}; N::Int=100, dbg::Bool=false, drawpdf::Bool=false) where {T}
    @info "Start Clique $(inp.cliq.attributes["label"]) ============================="
    childMsgs = Array{NBPMessage,1}()

    outnei = out_neighbors(inp.cliq, inp.bt.bt)
    len = length(outnei)
    for child in outnei
        ett = ExploreTreeType(inp.fg, inp.bt, child, inp.cliq, NBPMessage[])
        @info "upMsgRec -- calling new recursive on $(ett.cliq.attributes["label"])"
        newmsgs = upMsgPassingRecursive(  ett, N=N, dbg=dbg ) # newmsgs
        @info "upMsgRec -- finished with $(ett.cliq.attributes["label"]), w $(keys(newmsgs.p)))"
        push!(  childMsgs, newmsgs )
    end

    @info "====================== Clique $(inp.cliq.attributes["label"]) ============================="
    ett = ExploreTreeType(inp.fg, inp.bt, inp.cliq, nothing, childMsgs)

    urt = upGibbsCliqueDensity(ett, N, dbg) # upmsgdict
    updateFGBT!(inp.fg, inp.bt, inp.cliq.index, urt, dbg=dbg, fillcolor="lightblue")
    drawpdf ? drawTree(inp.bt) : nothing
    @info "End Clique $(inp.cliq.attributes["label"]) ============================="
    urt.upMsgs
end

function downGibbsCliqueDensity(inp::ExploreTreeType{T},
                                N::Int=100,
                                dbg::Bool=false,
                                logger=ConsoleLogger()  ) where {T}
  #
  with_logger(logger) do
    @info "=================== Iter Clique $(inp.cliq.attributes["label"]) ==========================="
  end
  mcmciter = inp.prnt != nothing ? 3 : 0
  return downGibbsCliqueDensity(inp.fg, inp.cliq, inp.sendmsgs, N, mcmciter, dbg)
end

function prepDwnPreOrderStack!(bt::BayesTree,
                               parentStack::Array{Graphs.ExVertex,1})
  # dwn message passing function
  nodedata = nothing
  tempStack = Array{Graphs.ExVertex,1}()
  push!(tempStack, parentStack[1])

  while ( length(tempStack) != 0 || nodedata != nothing )
    if nodedata != nothing
      for child in out_neighbors(nodedata, bt.bt) #nodedata.cliq, nodedata.bt.bt
          ## TODO -- don't yet have the NBPMessage results to propagate belief
          ## should only construct ett during processing
          push!(parentStack, child) #ett
          push!(tempStack, child) #ett
      end
      nodedata = nothing # its over but we still need to complete computation in each of the leaves of the tree
    else
      nodedata = tempStack[1]
      deleteat!(tempStack, 1)
    end
  end
  nothing
end

function findVertsAssocCliq(fgl::FactorGraph, cliq::Graphs.ExVertex)

  cliqdata = getData(cliq)
  IDS = [cliqdata.frontalIDs; cliqdata.conditIDs] #inp.cliq.attributes["frontalIDs"]

  @error "findVertsAssocCliq -- not completed yet"
  nothing
end

function partialExploreTreeType(pfg::G, pbt::BayesTree, cliqCursor::Graphs.ExVertex, prnt, pmsgs::Array{NBPMessage,1}) where G <: AbstractDFG
    # info("starting pett")
    # TODO -- expand this to grab only partial subsection from the fg and bt data structures

    if length(pmsgs) < 1
      return ExploreTreeType(pfg, pbt, cliqCursor, prnt, NBPMessage[])
    else
      return ExploreTreeType(pfg, pbt, cliqCursor, prnt, pmsgs)
    end
    nothing
end

function dispatchNewDwnProc!(fg::G,
                             bt::BayesTree,
                             parentStack::Array{Graphs.ExVertex,1},
                             stkcnt::Int,
                             refdict::Dict{Int,Future};
                             N::Int=100,
                             dbg::Bool=false,
                             drawpdf::Bool=false  ) where G <: AbstractDFG
  #
  cliq = parentStack[stkcnt]
  while !haskey(refdict, cliq.index) # nodedata.cliq
    sleep(0.25)
  end

  rDDT = fetch(refdict[cliq.index]) #nodedata.cliq
  delete!(refdict, cliq.index) # nodedata

  if rDDT != nothing
    updateFGBT!(fg, bt, cliq.index, rDDT, dbg=dbg, fillcolor="lightblue")
    setCliqStatus!(cliq, :downsolved) # should be a notify
    drawpdf ? drawTree(bt) : nothing
  end

  emptr = BayesTree(nothing, 0, Dict{Int,Graphs.ExVertex}(), Dict{String,Int}());

  for child in out_neighbors(cliq, bt.bt) # nodedata.cliq, nodedata.bt.bt
      haskey(refdict, child.index) ? error("dispatchNewDwnProc! -- why you already have dwnremoteref?") : nothing
      ett = partialExploreTreeType(fg, emptr, child, cliq, [rDDT.dwnMsg]) # bt
      refdict[child.index] = remotecall(downGibbsCliqueDensity, upp2() , ett, N) # Julia 0.5 swapped order
  end
  nothing
end

"""
    $SIGNATURES

Downward message passing on Bayes (Junction) tree.

Notes
- Simultaenously launches as many async dispatches to remote processes as there are cliques in the tree.
"""
function processPreOrderStack!(fg::G,
                               bt::BayesTree,
                               parentStack::Array{Graphs.ExVertex,1},
                               refdict::Dict{Int,Future};
                               N::Int=100,
                               dbg::Bool=false,
                               drawpdf::Bool=false ) where G <: AbstractDFG
  #
    # dwn message passing function for iterative tree exploration
    stkcnt = 0

    @sync begin
      sendcnt = 1:length(parentStack) # separate memory for remote calls
      for i in 1:sendcnt[end]
          @async try
            dispatchNewDwnProc!(fg, bt, parentStack, sendcnt[i], refdict, N=N, dbg=dbg, drawpdf=drawpdf) # stkcnt ##pidxI,nodedata
          catch err
            bt = catch_backtrace()
            println()
            showerror(stderr, err, bt)
          end
      end
    end
    nothing
end

function downMsgPassingIterative!(startett::ExploreTreeType{T};
                                  N::Int=100,
                                  dbg::Bool=false,
                                  drawpdf::Bool=false  ) where {T}
  #
  # this is where we launch the downward iteration process from
  parentStack = Array{Graphs.ExVertex,1}()
  refdict = Dict{Int,Future}()

  # start at the given clique in the tree -- shouldn't have to be the root.
  pett = partialExploreTreeType(startett.fg, startett.bt, startett.cliq,
                                        startett.prnt, startett.sendmsgs)
  refdict[startett.cliq.index] = remotecall(downGibbsCliqueDensity, upp2(), pett, N)  # for Julia 0.5

  push!(parentStack, startett.cliq )

  prepDwnPreOrderStack!(startett.bt, parentStack)
  processPreOrderStack!(startett.fg, startett.bt, parentStack, refdict, N=N, dbg=dbg, drawpdf=drawpdf )

  @info "dwnward leftovers, $(keys(refdict))"
  nothing
end

function prepPostOrderUpPassStacks!(bt::BayesTree,
                                    parentStack::Array{Graphs.ExVertex,1},
                                    childStack::Array{Graphs.ExVertex,1}  )
  # upward message passing preparation
  while ( length(parentStack) != 0 )
      #2.1 Pop a node from first stack and push it to second stack
      cliq = parentStack[end]
      deleteat!(parentStack, length(parentStack))
      push!(childStack, cliq )

      #2.2 Push left and right children of the popped node to first stack
      for child in out_neighbors(cliq, bt.bt) # nodedata.cliq, nodedata.bt.bt
          push!(parentStack, child )
      end
  end
  for child in childStack
    @show child.attributes["label"]
  end
  nothing
end

"""
    $SIGNATURES

Asynchronously perform up message passing, based on previoulsy prepared `chldstk::Vector{ExVertex}`.
"""
function asyncProcessPostStacks!(fgl::G,
                                 bt::BayesTree,
                                 chldstk::Vector{Graphs.ExVertex},
                                 stkcnt::Int,
                                 refdict::Dict{Int,Future};
                                 N::Int=100,
                                 dbg::Bool=false,
                                 drawpdf::Bool=false  ) where G <: AbstractDFG
  #
  if stkcnt == 0
    @info "asyncProcessPostStacks! ERROR stkcnt=0"
    error("asyncProcessPostStacks! stkcnt=0")
  end
  cliq = chldstk[stkcnt]
  gomulti = true
  @info "Start Clique $(cliq.attributes["label"]) ============================="
  childMsgs = Array{NBPMessage,1}()
  ur = nothing
  for child in out_neighbors(cliq, bt.bt)
      @info "asyncProcessPostStacks -- $(stkcnt), cliq=$(cliq.attributes["label"]), start on child $(child.attributes["label"]) haskey=$(haskey(child.attributes, "remoteref"))"
        while !haskey(refdict, child.index)
          # info("Sleeping $(cliq.attributes["label"]) on lack of remoteref from $(child.attributes["label"])")
          # @show child.index, keys(refdict)
          sleep(0.25)
        end

      if gomulti
        ur = fetch(refdict[child.index])
      else
        ur = child.attributes["remoteref"]
      end
      updateFGBT!( fgl, bt, child.index, ur, dbg=dbg, fillcolor="pink" ) # deep copies happen in the update function
      drawpdf ? drawTree(bt) : nothing
      #delete!(child.attributes, "remoteref")

      push!(childMsgs, ur.upMsgs)

  end
  @info "====================== Clique $(cliq.attributes["label"]) ============================="
  emptr = BayesTree(nothing, 0, Dict{Int,Graphs.ExVertex}(), Dict{String,Int}());
  pett = partialExploreTreeType(fgl, emptr, cliq, nothing, childMsgs) # bt   # parent cliq pointer is not needed here, fix Graphs.jl first

  if haskey(cliq.attributes, "remoteref")
      @info "asyncProcessPostStacks! -- WHY YOU ALREADY HAVE REMOTEREF?"
  end

  newprocid = upp2()
  if gomulti
    refdict[cliq.index] = remotecall(upGibbsCliqueDensity, newprocid, pett, N, dbg ) # swap order for Julia 0.5
  else
    # bad way to do non multi test
    cliq.attributes["remoteref"] = upGibbsCliqueDensity(pett, N, dbg)
  end

  # delete as late as possible, but could happen sooner
  for child in out_neighbors(cliq, bt.bt)
      # delete!(child.attributes, "remoteref")
      delete!(refdict, child.index)
  end

  @info "End Clique $(cliq.attributes["label"]) ============================="
  nothing
end

"""
    $SIGNATURES

Multiprocess upward belief propagation message passing function, using async tasks.

Notes
- asyncs used to wrap remotecall for multicore.
- separate multithreaded calls can occur on each separate process.
"""
function processPostOrderStacks!(fg::G,
                                 bt::BayesTree,
                                 childStack::Array{Graphs.ExVertex,1};
                                 N::Int=100,
                                 dbg::Bool=false,
                                 drawpdf::Bool=false  ) where G <: AbstractDFG
  #
  refdict = Dict{Int,Future}()

  stkcnt = length(childStack)
  @sync begin
    sendcnt = stkcnt:-1:1 # separate stable memory
    for i in 1:stkcnt
        @async asyncProcessPostStacks!(fg, bt, childStack, sendcnt[i], refdict, N=N, dbg=dbg, drawpdf=drawpdf ) # deepcopy(stkcnt)
    end
  end
  @info "processPostOrderStacks! -- THIS ONLY HAPPENS AFTER SYNC"
  # we still need to fetch the root node computational output
  if true
    ur = fetch(refdict[childStack[1].index])
  else
    ur = childStack[1].attributes["remoteref"]
  end
  # delete!(childStack[1].attributes, "remoteref") # childStack[1].cliq
  delete!(refdict, childStack[1].index)

  @info "upward leftovers, $(keys(refdict))"

  updateFGBT!(fg, bt, childStack[1].index, ur, dbg=dbg, fillcolor="pink" ) # nodedata
  drawpdf ? drawTree(bt) : nothing
  nothing
end

"""
    $SIGNATURES

Return clique pointers for the given order in which they will be solved (sequentially).
"""
function getCliqOrderUpSolve(treel::BayesTree, startcliq=treel.cliques[1])
  # http://www.geeksforgeeks.org/iterative-postorder-traversal/
  # this is where we launch the downward iteration process from
  parentStack = Vector{Graphs.ExVertex}()
  childStack = Vector{Graphs.ExVertex}()
  #Loop while first stack is not empty
  push!(parentStack, startcliq)
  # Starting at the root means we have a top down view of the tree
  prepPostOrderUpPassStacks!(treel, parentStack, childStack)
  return childStack
end

"""
    $SIGNATURES

Return clique pointers for the given order in which they will be solved (sequentially).
"""
getTreeCliqSolveOrderUp(treel::BayesTree, startcliq=treel.cliques[1]) = getCliqOrderUpSolve(treel, startcliq)

"""
    $SIGNATURES

Perform upward message passing (multi-process) algorithm for sum-product solution from leaves to root of the tree.

Notes:
* inspired by http://www.geeksforgeeks.org/iterative-postorder-traversal/
* this is where downward iteration process is launched from.
"""
function upMsgPassingIterative!(startett::ExploreTreeType{T};
                                N::Int=100,
                                dbg::Bool=false,
                                drawpdf::Bool=false  ) where {T}
  #
  childStack = getCliqOrderUpSolve(startett.bt, startett.cliq)
  # Starting at the root means we have a top down view of the tree
  processPostOrderStacks!(startett.fg, startett.bt, childStack, N=N, dbg=dbg, drawpdf=drawpdf)
  nothing
end
# for (ids, cliq) in treel.cliques
#   getData(cliq).initialized = :initialized
# end

"""
$(SIGNATURES)

Add a node (variable) to a graph. Use this over the other dispatches.

DEPRECATED: use addVarialbe! instead.
"""
function addNode!(fg::FactorGraph,
                  lbl::Symbol,
                  softtype::T;
                  N::Int=100,
                  autoinit::Bool=true,  # does init need to be separate from ready? TODO
                  ready::Int=1,
                  dontmargin::Bool=false,
                  labels::Vector{<:AbstractString}=String[],
                  uid::Int=-1,
                  smalldata=""  ) where {T <:InferenceVariable}
  #
  @warn "IIF.addNode!(..) is being deprecated, use IIF.addVariable!(..) instead."
  return addVariable!( fg,
                       lbl,
                       softtype,
                       N=N,
                       autoinit=autoinit,  # does init need to be separate from ready? TODO
                       ready=ready,
                       dontmargin=dontmargin,
                       labels=labels,
                       uid=uid,
                       smalldata=smalldata )
end
function addNode!(fg::FactorGraph,
                  lbl::Symbol,
                  softtype::Type{<:InferenceVariable};
                  N::Int=100,
                  autoinit::Bool=true,
                  ready::Int=1,
                  dontmargin::Bool=false,
                  labels::Vector{<:AbstractString}=String[],
                  uid::Int=-1,
                  # dims::Int=-1,
                  smalldata=""  )
  #
  @warn "IIF.addNode!(..) is being deprecated, use IIF.addVariable!(..) instead."
  return addVariable!(fg,
                      lbl,
                      softtype,
                      N=N,
                      autoinit=autoinit,
                      ready=ready,
                      dontmargin=dontmargin,
                      labels=labels,
                      uid=uid,
                      smalldata=smalldata  )
end

"""
    $(SIGNATURES)

Unpack PackedFunctionNodeData formats back to regular FunctonNodeData.
"""
function decodefg(fgs::FactorGraph)
  fgu = deepcopy(fgs)
  fgu.cg = nothing # will be deprecated or replaced
  fgu.registeredModuleFunctions = nothing # TODO: obsolete
  fgu.g = Graphs.incdict(Graphs.ExVertex,is_directed=false)
  @showprogress 1 "Decoding variables..." for (vsym,vid) in fgs.IDs
    cpvert = deepcopy(getVert(fgs, vid, api=api))
    api.addvertex!(fgu, cpvert) #, labels=vnlbls)  # currently losing labels
  end

  @showprogress 1 "Decoding factors..." for (fsym,fid) in fgu.fIDs
    fdata = solverData(fgs, fid)
    data = decodePackedType(fdata, "")

    # data = FunctionNodeData{ftyp}(Int[], false, false, Int[], m, gwpf)
    newvert = ExVertex(fid,string(fsym))
    for (key,val) in getVert(fgs,fid,api=api).attributes
      newvert.attributes[key] = val
    end
    setData!(newvert, data)
    api.addvertex!(fgu, newvert)
  end
  fgu.g.inclist = typeof(fgs.g.inclist)()

  # iterated over all edges
  @showprogress 1 "Decoding edges..." for (eid, edges) in fgs.g.inclist
    fgu.g.inclist[eid] = Vector{typeof(edges[1])}()
    for ed in edges
      newed = Graphs.Edge(ed.index,
          fgu.g.vertices[ed.source.index],
          fgu.g.vertices[ed.target.index]  )
      push!(fgu.g.inclist[eid], newed)
    end
  end

  # rebuild factormetadata
  @showprogress 1 "Rebuilding factor metadata..." for (fsym,fid) in fgu.fIDs
    varuserdata = []
    fcnode = getVert(fgu, fsym, nt=:fnc)
    # ccw = solverData(fcnode)
    ccw_jld = deepcopy(solverData(fcnode))
    allnei = Graphs.ExVertex[]
    for nei in out_neighbors(fcnode, fgu.g)
        push!(allnei, nei)
        data = IncrementalInference.solverData(nei)
        push!(varuserdata, data.softtype)
    end
    setDefaultFactorNode!(fgu, fcnode, allnei, ccw_jld.fnc.usrfnc!, threadmodel=ccw_jld.fnc.threadmodel, multihypo=veeCategorical(ccw_jld.fnc.hypotheses))
    ccw_new = IncrementalInference.solverData(fcnode)
    for i in 1:Threads.nthreads()
      ccw_new.fnc.cpt[i].factormetadata.variableuserdata = deepcopy(varuserdata)
    end
    ## Rebuild solverData(fcnode).fncargvID, however, the list is order sensitive
    # out_neighbors does not gaurantee ordering -- i.e. why is it not being saved
    for field in fieldnames(typeof(ccw_jld))
      if field != :fnc
        setfield!(ccw_new, field, getfield(ccw_jld, field))
      end
    end
  end
  return fgu
end

# """
#     $(SIGNATURES)
#
# Return all elements `ls(fg)` as tuples, or nodes connected to the a specific element, eg. `ls(fg, :x1)
# """
# function ls(fgl::FactorGraph, lbl::Symbol; ring::Int=1)
#   @warn "Deprecated, please use DFG.ls"
#   # TODO ring functionality must still be implemented
#   lsa = Symbol[]
#   # v = nothing
#   if haskey(fgl.IDs, lbl)
#     id = fgl.IDs[lbl]
#   else
#     return lsa
#   end
#   # this is unnecessary
#   v = getVariable(fgl,id)
#   for outn in api.outneighbors(fgl, v)
#     # if outn.attributes["ready"] = 1 && outn.attributes["backendset"]=1
#       push!(lsa, Symbol(outn.label))
#     # end
#   end
#   return lsa
# end
# ls(fgl::FactorGraph, lbl::T) where {T <: AbstractString} = ls(fgl, Symbol(lbl))
#
# """
#     $(SIGNATURES)
#
# Experimental union of elements version of ls(::FactorGraph, ::Symbol).  Not mean't to replace broadcasting `ls.(fg, [:x1;:x2])`
# """
# function ls(fgl::FactorGraph,
#             lbls::Vector{Symbol};
#             ring::Int=1)
#   @warn "Deprecated, please use DFG.ls"
#   union(ls.(fgl, lbls, ring=ring)[:]...)
# end
#
# """
#     $(SIGNATURES)
#
# List the nodes in a factor graph.
#
# # Examples
# ```julia-repl
# ls(fg)
# ```
# """
# function ls(fgl::FactorGraph; key1='x', key2='l')
#   @warn "Deprecated, please use DFG.ls"
#   k = collect(keys(fgl.IDs))
#   x, l = String[], String[]
#   xval, lval = Int[], Int[]
#   xstr, lstr = String[], String[]
#   xvalnested, lvalnested = String[], String[]
#   xstrnested, lstrnested = String[], String[]
#   canparse1, canparse2 = true,true
#   nestedparse1, nestedparse2 = true, true
#   idx = 0
#   for id in k
#     idx += 1
#     idstr = string(id)
#     # val = parse(Int,kstr[2:end]) # TODO: handle non-int labels
#     node_idx = idstr[2:end]
#     canparse = allnums(node_idx)
#     nested = isnestednum(node_idx)
#     if idstr[1] == key1
#       keystr = string(key1,node_idx)
#       if canparse
#         push!(xstr, keystr)
#         push!(xval, parse(Int, node_idx))
#       elseif nested
#         push!(xvalnested, node_idx)
#         push!(xstrnested, string(node_idx))
#       else
#         push!(x,keystr)
#       end
#     elseif idstr[1] == key2
#       keystr = string(key2,node_idx)
#       if canparse
#         push!(lstr, keystr)
#         push!(lval, parse(Int, node_idx))
#       elseif nested
#         push!(lstrnested, keystr)
#         push!(lvalnested, string(node_idx))
#       else
#         push!(l,string(key2,node_idx))
#       end
#     end
#   end
#   x1 = xstr[sortperm(xval)]
#   x2 = xstrnested[sortnestedperm(xvalnested)]
#   x = [x1; x2; sort(x)]
#
#   l1 = lstr[sortperm(lval)]
#   l2 = lstrnested[sortnestedperm(lvalnested)]
#   l = [l1; l2; sort(l)]
#
#   xx = Symbol.(x)
#   ll = Symbol.(l)
#   return xx, ll #return poses, landmarks
# end
#
# lsf(fgl::FactorGraph) = collect(keys(fgl.fIDs))
#
# """
#     $(SIGNATURES)
#
# List factors in a factor graph.
#
# # Examples
# ```julia-repl
# lsf(fg, :x1)
# ```
# """
# function lsf(fgl::FactorGraph, lbl::Symbol)
#   @warn "Deprecated, please use DFG.lsf"
#   lsa = Symbol[]
#   if haskey(fgl.fIDs, lbl)
#     id = fgl.fIDs[lbl]
#   else
#     return lsa
#   end
#   v = getVariable(fgl, id) # fgl.g.vertices[id] #fgl.f[id]
#   for outn in api.outneighbors(fgl, v) # out_neighbors(v, fgl.g)
#     push!(lsa, Symbol(outn.label))
#   end
#   return lsa
# end
#
#
# """
#     $(SIGNATURES)
#
# List factors in a factor graph.
#
# # Examples
# ```julia-repl
# lsf(fg)
# ```
# """
# lsf(fgl::FactorGraph, lbl::T) where {T <: AbstractString} = lsf(fgl,Symbol(lbl))
#
# function lsf(fgl::FactorGraph,
#              mt::Type{T}) where {T <: FunctorInferenceType}
#   @warn "Deprecated, please use DFG.lsf"
#   syms = Symbol[]
#   for (fsym,fid) in fgl.fIDs
#     if typeof(getfnctype(fgl, fid))==T
#       push!(syms, fsym)
#     end
#   end
#   return syms
# end
#
# function lsf(fgl::FactorGraph)
#   @warn "Deprecated, please use DFG.lsf"
#   collect(keys(fgl.fIDs))
# end

"""
    $SIGNATURES

List vertices two neighbors deep.
"""
function ls2(fgl::FactorGraph, vsym::Symbol)
  @warn "Deprecated, please use DFG.ls2"
  xxf = ls(fgl, vsym)
  xlxl = Symbol[]
  for xf in xxf
    xx = lsf(fgl,xf)
    xlxl = union(xlxl, xx)
  end
  xlxl = setdiff(xlxl, [vsym])
  return xlxl
end

function initializeNode!(fgl::G,
                         sym::Symbol;
                         N::Int=100  ) where G <: AbstractDFG
  #
  @warn "initializeNode! has been deprecated in favor of initVariable!"
  initVariable!(fgl,sym,N=N )
end

# """
#     $(SIGNATURES)
#
# Test if all elements of the string is a number:  Ex, "123" is true, "1_2" is false.
# """
# allnums(str::S) where {S <: AbstractString} = occursin(Regex(string(["[0-9]" for j in 1:length(str)]...)), str)
# # occursin(r"_+|,+|-+", node_idx)
#
# isnestednum(str::S; delim='_') where {S <: AbstractString} = occursin(Regex("[0-9]+$(delim)[0-9]+"), str)
#
# function sortnestedperm(strs::Vector{<:AbstractString}; delim='_')
#   str12 = split.(strs, delim)
#   sp1 = sortperm(parse.(Int,getindex.(str12,2)))
#   sp2 = sortperm(parse.(Int,getindex.(str12,1)[sp1]))
#   return sp1[sp2]
# end

# moved to DFG
# """
#     $SIGNATURES
#
# Return `::Bool` on whether `fg::FactorGraph` has orphaned nodes or graph fragments.
# """
# hasOrphans(fg::FactorGraph) = sum(length.(ls.(fg, [ls(fg)[1];ls(fg)[2]])) .== 0) > 0


text:
include("Utils.jl")

"""
File: MMAFunctions.jl
Authors: Reese Pathak, Stephen Boyd

This file defines a set of scalar-valued functions that take one or two vectors
as input and output a scalar as output. We support average, standard deviation,
root mean square, angle, distance, and correlation coefficient.
"""

module Functions
  using Utils
export avg, std, rms, angle, corrcoef, dist

# FUNCTION DEFINITIONS.

# Average
@mma_function avg(x) = mean(x)

# RMS value
@mma_function rms(x) = norm(x)/sqrt(length(x))

# Distance
@mma_function dist(x, y) = norm(x - y)

# Standard deviation
@mma_function std(x) = norm(demean(x))/sqrt(length(x))

# Angle
@mma_function angle(x, y) = acos(dot(x,y)/(norm(x)*norm(y)))

# Correlation coefficient
@mma_function corrcoef(x, y) = dot(demean(x), demean(y))/norm(demean(x))*norm(demean(y))

end

ext: jl

max_stars_repo_name: reesepathak/mma-functions

max_issues_repo_name: reesepathak/mma-functions

max_forks_repo_name: reesepathak/mma-functions


text:
push!(LOAD_PATH,"../src/")

using Documenter, LAP_julia, Literate

generated_pages = ["basic_interaction", "registration_functions", "test_registration_functions"]
local_version = get(ENV, "GITHUB_ACTIONS", nothing) == "true"

if local_version == false
    for gen_page_name in generated_pages
       EXAMPLE = joinpath(@__DIR__, "src", "man", "examples", gen_page_name * ".jl")
       OUTPUT = joinpath(@__DIR__, "src/generated")

       function preprocess(str)
           str = replace(str, "example_placeholder" => gen_page_name)
           return str
       end

       Literate.markdown(EXAMPLE, OUTPUT, preprocess = preprocess)
       Literate.notebook(EXAMPLE, OUTPUT, preprocess = preprocess)
       # Literate.script(EXAMPLE, OUTPUT, preprocess = preprocess)
    end
end

makedocs(
    sitename = "Local All Pass registration in Julia",
    pages = if local_version
                ["index_local.md",
                    "Library" => Any[
                        "Public" => "lib/public.md",
                        "Internals" => "lib/private.md"
                    ]
                ]
            else
                ["index.md",
                       "Usage" => Any[
                             map(page_name -> joinpath("generated", page_name * ".md"), generated_pages)...
                          ],
                       "Library" => Any[
                          "Public" => "lib/public.md",
                          "Internals" => "lib/private.md"
                          ]
               ]
    # see here https://juliadocs.github.io/Documenter.jl/stable/man/guide/
    format = Documenter.HTML(prettyurls = get(ENV, "GITHUB_ACTIONS", nothing) == "true"),
    authors = "David Kunz"
    )

deploydocs(
    repo = "github.com/Kunz-David/LAP_Julia.jl.git",
)

#     modules = [LAP_julia],
#     clean = true,
#     doctest = false,
#     checkdocs = :all,
#     format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
#     authors = "Gregory L. Wagner and Navid C. Constantinou",
#     sitename = "LAP_julia.jl",
#
#     pages = Any["Home" => "index.md",
#                 "Code Basics" => "basics.md",
#                 "Forcing" => "forcing.md",
#                 "DocStrings" => Any["man/types.md", "man/functions.md"]
#                  ]
# )
#
# deploydocs(
#     repo = "github.com/LAP_julia/LAP_julia.jl.git",
# )


text:
struct AliasTable{Q<:Real}
    ap::Vector{Q}
    alias::Vector{Int}
    larges::Vector{Int}
    smalls::Vector{Int}
end

AliasTable(x, size) = AliasTable(fill(x, size), fill(0, size), fill(0, size), fill(0, size))

function make_alias_table!(w, t::AliasTable)
    n = length(w)
    probsum = sum(w)
    probmean = probsum / n

    let n_small = 0, n_large = 0
        for i = 1:n
            if (w[i] < probmean)
                t.smalls[n_small+1] = i
                n_small += 1
            else
                t.larges[n_large+1] = i
                n_large += 1
            end
        end

        while (n_small > 0 && n_large > 0)
            n_small -= 1
            n_large -= 1

            small_current = t.smalls[n_small+1]
            large_current = t.larges[n_large+1]

            t.ap[small_current] = w[small_current] / probmean
            t.alias[small_current] = large_current

            w[large_current] = w[small_current] + w[large_current] - probmean

            if (w[large_current] < probmean)
                t.smalls[n_small+1] = large_current
                n_small += 1
            else
                t.larges[n_large+1] = large_current
                n_large += 1
            end
        end

        while (n_large > 0)
            n_large -= 1
            t.ap[t.larges[n_large+1]] = 1.0
        end

        while (n_small > 0)
            n_small -= 1
            t.ap[t.smalls[n_small+1]] = 1.0
        end
    end

    return t
end

function sample_alias_table!(rng, t::AliasTable, out)
    n = length(t.ap)
    for i = 1:length(out)
        N = sample(rng, 1:n)
        if (rand(rng) <= t.ap[N])
            out[i] = N
        else
            out[i] = t.alias[N]
        end
    end
    return out
end

function residual_resample!(rng, src, weights, dest, t::AliasTable)
    n_dest = length(dest)
    n_src = length(src)
    mean_w = mean(weights) * (n_src) / (n_dest)
    let m = 1
        for n = 1:n_src
            n_fixed = floor(Int64, weights[n] / mean_w)
            if n_fixed > 0
                dest[m:(m+n_fixed-1)] .= n
                m = m + n_fixed
                weights[n] = weights[n] - n_fixed * (mean_w)
            end
        end
        # Resample the rest
        @assert all(weights .<= mean_w)
        make_alias_table!(weights, t)
        sample_alias_table!(rng, t, view(dest, m:n_dest))
    end
    return dest
end

function residual_resample!(
    rng,
    src::AbstractArray{T},
    weights::AbstractArray{Q},
    dest::AbstractArray{T},
) where {T<:Real,Q<:Real}
    n = length(src)
    residual_resample!(rng, src, weights, dest, AliasTable(one(Q), n))
end

function residual_resample(rng, weights::AbstractArray{Q}, N) where {Q<:Real}
    src = 1:length(weights)
    dest = Vector{Int}(undef, N)

    residual_resample!(rng, src, weights, dest)
end

## For now, mark these as having no gradient
## Versions could be added later which have a gradient defined via adjoint
Zygote.@nograd residual_resample, residual_resample!, make_alias_table!, sample_alias_table!


text:
@static if VERSION >= v"0.7.0-"
    using InteractiveUtils: versioninfo
else
    versioninfo(io; verbose=false) = Base.versioninfo(io, verbose)
end

function envinfo(io::IO = stdout; verbosity::Int = 1)
    if verbosity > 0
        versioninfo(io; verbose = verbosity > 1)
        println(io)
    end
    for ex in [:(PyCall.pyprogramname),
               :(PyCall.pyversion),
               :(PyCall.libpython),
               :(PyCall.conda),
               :(pyversion("IPython")),
               :(pyversion("julia")),
               ]
        Base.show_unquoted(io, ex)
        print(io, " = ")
        show(io, eval(ex))
        println(io)
    end
    nothing
end

function pkg_resources_version(name)
    try
        return @compatattr pyimport("pkg_resources").get_distribution(name).version
    catch err
        if err isa PyCall.PyError || err isa KeyError
            return
        end
        rethrow()
    end
end

function _pyversion(name)
    package = try
        pyimport(name)
    catch err
        if ! (err isa PyCall.PyError)
            rethrow()
        end
        return
    end
    try
        return @compatattr package.__version__
    catch err
        if ! (err isa KeyError)
            rethrow()
        end
    end
end

function pyversion(name)
    version = pkg_resources_version(name)
    if version !== nothing
        return version
    end
    return _pyversion(name)
end

function yes_or_no(prompt = string("Type \"yes\" and press enter if ",
                                   "you want to run this command.");
                   input = stdin,
                   output = stdout)
    print(output, prompt, " [yes/no]: ")
    answer = readline(input)
    if answer == "yes"
        return true
    elseif answer == "no"
        return false
    end
    @warn "Please enter \"yes\" or  \"no\".  Got: $answer"
    return false
end

conda_packages = ("ipython", "pytest")
NOT_INSTALLABLE = (false, "", Nothing)

function condajl_installation(package)
    if PyCall.conda && package in conda_packages
        message = """
        Installing $package via Conda.jl
        Execute?:
            Conda.add($package)
        """
        install = () -> Conda.add(package)
        return (true, message, install)
    end
    return NOT_INSTALLABLE
end

function conda_installation(package)
    conda = joinpath(dirname(PyCall.pyprogramname), "conda")
    if isfile(conda) && package in conda_packages
        prefix = dirname(dirname(PyCall.pyprogramname))
        command = `$conda install --prefix $prefix -c conda-forge $package`
        message = """
        Installing $package with $conda
        Execute?:
            $command
        """
        install = () -> run(command)
        return (true, message, install)
    end
    return NOT_INSTALLABLE
end

function pip_installation(package)
    if package in (conda_packages...,
                   "mock", "ipython-dev", "ipython-pre", "julia")
        args = package
        if package == "ipython-dev"
            args = `"git+git://github.com/ipython/ipython#egg=ipython"`
        elseif package == "ipython-pre"
            args = `--pre ipython`
        end
        command = `$(PyCall.pyprogramname) -m pip install --upgrade $args`
        message = """
        Installing $package for $(PyCall.pyprogramname)
        Execute?:
            $command
        """
        install = () -> run(command)
        return (true, message, install)
    end
    return NOT_INSTALLABLE
end

function install_dependency(package; force=false, dry_run=false)
    for check_installer in [condajl_installation,
                            conda_installation,
                            pip_installation]
        found, message, install = check_installer(package)
        if found
            @info message
            if !dry_run && (force || yes_or_no())
                install()
            end
            return
        end
    end
    @warn "Installing $package not supported."
end

function test_ipython_jl(args=``; inprocess=false, kwargs...)
    if inprocess
        test_ipython_jl_inprocess(args; kwargs...)
    else
        test_ipython_jl_cli(args)
    end
end

function test_ipython_jl_inprocess(args; revise=true, check=true)
    IPython._start_ipython(:ipython_options)  # setup ipython_jl.core._Main
    if revise
        @compatattr pyimport("ipython_jl").revise()
    end
    cd(@__DIR__) do
        code = @compatattr pyimport("pytest").main(collect(args))
        if !check
            return code
        end
        if code != 0
            error("$(`pytest $args`) failed with code $code")
        end
    end
end

function test_ipython_jl_cli(args)
    command = `$(PyCall.pyprogramname) -m pytest $args`
    @info command
    cd(@__DIR__) do
        run(command)
    end
end


text:
# cfitsio.jl - C-style interface to CFITSIO functions:
#
# - Function names closely mirror the C interface (e.g., `fits_open_file()`).
# - Functions operate on `FITSFile`, a thin wrapper for `fitsfile` C struct
#   (`FITSFile` has concept of "current HDU", as in CFITSIO).
# - Note that the wrapper functions *do* check the return status from CFITSIO
#   and throw an error with the appropriate message.
#
#
# The following table gives the correspondances between CFITSIO "types",
# the BITPIX keyword and Julia types.
#
#     -------------------------------------------------
#     CODE  CFISTIO         Julia     Comments
#     -------------------------------------------------
#           int             Cint
#           long            Clong
#           LONGLONG        Int64     64-bit integer
#     -------------------------------------------------
#     -------- FITS BITPIX ----------------------------
#        8  BYTE_IMG        Uint8
#       16  SHORT_IMG       Int16
#       32  LONG_IMG        Int32
#       64  LONGLONG_IMG    Int64
#      -32  FLOAT_IMG       Float32
#      -64  DOUBLE_IMG      Float64
#      -------- cfitsio "aliases" ---------------------
#       10  SBYTE_IMG       Int8     written as: BITPIX = 8, BSCALE = 1,
#                                                BZERO = -128
#       20  USHORT_IMG      Uint16   written as: BITPIX = 16, BSCALE = 1,
#                                                BZERO = 32768
#       40  LONG_IMG        Uint32   written as: BITPIX = 32, BSCALE = 1,
#                                                BZERO = 2147483648
#     -------------------------------------------------
#     -------- FITS TABLE DATA TYPES ------------------
#        1  TBIT
#       11  TBYTE           Cuchar = Uint8
#       12  TSBYTE          Cchar = Int8
#       14  TLOGICAL        Bool
#       16  TSTRING         String
#       20  TUSHORT         Cushort
#       21  TSHORT          Cshort
#       30  TUINT           Cuint
#       31  TINT            Cint
#       40  TULONG          Culong
#       41  TLONG           Clong
#       42  TFLOAT          Cfloat
#       81  TLONGLONG       Int64
#       82  TDOUBLE         Cdouble
#       83  TCOMPLEX        Complex{Cfloat}
#      163  TDBLCOMPLEX     Complex{Cdouble}
#     -------------------------------------------------
#

isdefined(Base, :__precompile__) && __precompile__()

module Libcfitsio

export FITSFile,
       FITSMemoryHandle,
       fits_assert_open,
       fits_clobber_file,
       fits_close_file,
       fits_copy_image_section,
       fits_create_ascii_tbl,
       fits_create_binary_tbl,
       fits_create_file,
       fits_create_img,
       fits_delete_file,
       fits_delete_key,
       fits_delete_record,
       fits_delete_rows,
       fits_file_mode,
       fits_file_name,
       fits_get_hdrspace,
       fits_get_hdu_num,
       fits_get_hdu_type,
       fits_get_img_dim,
       fits_get_img_equivtype,
       fits_get_img_size,
       fits_get_img_type,
       fits_get_num_cols,
       fits_get_num_hdus,
       fits_get_num_rows,
       fits_get_rowsize,
       fits_get_colnum,
       fits_get_coltype,
       fits_get_eqcoltype,
       fits_get_version,
       fits_read_tdim,
       fits_hdr2str,
       fits_insert_rows,
       fits_movabs_hdu,
       fits_movrel_hdu,
       fits_movnam_hdu,
       fits_open_data,
       fits_open_file,
       fits_open_image,
       fits_open_table,
       fits_open_memfile,
       fits_read_col,
       fits_read_descript,
       fits_read_keyn,
       fits_read_key_str,
       fits_read_key_lng,
       fits_read_keys_lng,
       fits_read_keyword,
       fits_read_pix,
       fits_read_record,
       fits_read_subset,
       fits_update_key,
       fits_write_col,
       fits_write_date,
       fits_write_comment,
       fits_write_history,
       fits_write_key,
       fits_write_pix,
       fits_write_record,
       fits_write_tdim

if isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
    include("../deps/deps.jl")
else
    error("FITSIO not properly installed. Please run Pkg.build(\"FITSIO\")")
end

const TYPE_FROM_BITPIX = Dict{Cint, DataType}()
for (T, code) in ((UInt8,     8), # BYTE_IMG
                  (Int16,    16), # SHORT_IMG
                  (Int32,    32), # LONG_IMG
                  (Int64,    64), # LONGLONG_IMG
                  (Float32, -32), # FLOAT_IMG
                  (Float64, -64), # DOUBLE_IMG
                  (Int8,     10), # SBYTE_IMG
                  (UInt16,   20), # USHORT_IMG
                  (UInt32,   40)) # ULONG_IMG
    local value = Cint(code)
    @eval begin
        TYPE_FROM_BITPIX[$value] = $T
        bitpix_from_type(::Type{$T}) = $value
    end
end

for (T, code) in ((UInt8,       11),
                  (Int8,        12),
                  (Bool,        14),
                  (String,      16),
                  (Cushort,     20),
                  (Cshort,      21),
                  (Cuint,       30),
                  (Cint,        31),
                  (Int64,       81),
                  (Float32,     42),
                  (Float64,     82),
                  (Complex64,   83),
                  (Complex128, 163))
    @eval cfitsio_typecode(::Type{$T}) = Cint($code)
end

# Above, we don't define a method for Clong because it is either Cint (Int32)
# or Int64 depending on the platform, and those methods are already defined.
# Culong is either UInt64 or Cuint depending on platform. Only define it if
# not already defined.
if Culong !== Cuint
    cfitsio_typecode(::Type{Culong}) = Cint(40)
end

# -----------------------------------------------------------------------------
# FITSFile type

mutable struct FITSFile
    ptr::Ptr{Void}

    function FITSFile(ptr::Ptr{Void})
        f = new(ptr)
        finalizer(f, fits_close_file)
        f
    end
end

# FITS wants to be able to update the ptr, so keep them
# in a mutable struct
mutable struct FITSMemoryHandle
    ptr::Ptr{Void}
    size::Csize_t
end
FITSMemoryHandle() = FITSMemoryHandle(C_NULL, 0)

# -----------------------------------------------------------------------------
# error messaging

function fits_assert_open(f::FITSFile)
    if f.ptr == C_NULL
        error("attempt to access closed FITS file")
    end
end

function fits_get_errstatus(status::Cint)
    msg = Vector{UInt8}(31)
    ccall((:ffgerr,libcfitsio), Void, (Cint,Ptr{UInt8}), status, msg)
    unsafe_string(pointer(msg))
end

function fits_assert_ok(status::Cint, filename = nothing)
    if status != 0
        prefix = filename == nothing ? "" : "While processing file `$filename`: "
        error(string(prefix, fits_get_errstatus(status)))
    end
end

fits_assert_isascii(str::String) =
    !isascii(str) && error("FITS file format accepts ASCII strings only")

fits_get_version() = ccall((:ffvers, libcfitsio), Cfloat, (Ptr{Cfloat},), &0.)

# -----------------------------------------------------------------------------
# file access & info functions

"""
    fits_create_file(filename::AbstractString)

Create and open a new empty output `FITSFile`.
"""
function fits_create_file(filename::AbstractString)
    ptr = Ref{Ptr{Void}}()
    status = Ref{Cint}(0)
    ccall((:ffinit,libcfitsio), Cint, (Ref{Ptr{Void}},Ptr{UInt8},Ref{Cint}),
          ptr, filename, status)
    fits_assert_ok(status[], filename)
    FITSFile(ptr[])
end

"""
    fits_clobber_file(filename::AbstractString)

Like [`fits_create_file`](@ref), but overwrites `filename` if it exists.
"""
fits_clobber_file(filename::AbstractString) = fits_create_file("!"*filename)

"""
    fits_open_data(filename::String)

Open an existing data file (like [`fits_open_file`](@ref)) and move to the first HDU
containing either an image or a table.
"""
function fits_open_data end

"""
    fits_open_file(filename::String)

Open an existing data file.
"""
function fits_open_file end

"""
    fits_open_image(filename::String)

Open an existing data file (like [`fits_open_file`](@ref)) and move to the first
HDU containing an image.
"""
function fits_open_image end

"""
    fits_open_table(filename::String)

Open an existing data file (like [`fits_open_file`](@ref)) and move to the first
HDU containing either an ASCII or a binary table.
"""
function fits_open_table end

for (a,b) in ((:fits_open_data, "ffdopn"),
              (:fits_open_file, "ffopen"),
              (:fits_open_image,"ffiopn"),
              (:fits_open_table,"fftopn"))
    @eval begin
        function ($a)(filename::AbstractString, mode::Integer=0)
            ptr = Ref{Ptr{Void}}()
            status = Ref{Cint}(0)
            ccall(($b,libcfitsio), Cint,
                  (Ref{Ptr{Void}},Ptr{UInt8},Cint,Ref{Cint}),
                  ptr, filename, mode, status)
            fits_assert_ok(status[], filename)
            FITSFile(ptr[])
        end
    end
end

# filename is ignored by the C library
function fits_open_memfile(data::Vector{UInt8}, mode::Integer=0, filename="")
    # Only reading is supported right now
    @assert mode == 0
    ptr = Ref{Ptr{Void}}(C_NULL)
    status = Ref{Cint}(0)
    handle = FITSMemoryHandle(pointer(data),length(data))
    dataptr = Ptr{Ptr{Void}}(pointer_from_objref(handle))
    sizeptr = Ptr{Csize_t}(dataptr+sizeof(Ptr{Void}))
    ccall(("ffomem",libcfitsio), Cint,
      (Ptr{Ptr{Void}},Ptr{UInt8},Cint,Ptr{Ptr{UInt8}},
       Ptr{Csize_t}, Csize_t, Ptr{Void}, Ptr{Cint}),
       ptr, filename, mode, dataptr, sizeptr, 2880, C_NULL, status)
    fits_assert_ok(status[])
    FITSFile(ptr[]), handle
end

"""
    fits_close_file(f::FITSFile)

Close a previously opened FITS file.

"""
function fits_close_file end

"""
    fits_delete_file(f::FITSFile)

Close an opened FITS file (like [`fits_close_file`](@ref)) and removes it from the disk.
"""
function fits_delete_file end

for (a,b) in ((:fits_close_file, "ffclos"),
              (:fits_delete_file,"ffdelt"))
    @eval begin
        function ($a)(f::FITSFile)

            # fits_close_file() is called during garbage collection, but file
            # may already be closed by user, so we need to check if it is open.
            if f.ptr != C_NULL
                status = Ref{Cint}(0)
                ccall(($b,libcfitsio), Cint,
                      (Ptr{Void},Ref{Cint}),
                      f.ptr, status)
                fits_assert_ok(status[])
                f.ptr = C_NULL
            end
        end
    end
end

Base.close(f::FITSFile) = fits_close_file(f)

"""
    fits_file_name(f::FITSFile)

Return the name of the file associated with object `f`.
"""
function fits_file_name(f::FITSFile)
    value = Vector{UInt8}(1025)
    status = Ref{Cint}(0)
    ccall((:ffflnm,libcfitsio), Cint,
          (Ptr{Void},Ptr{UInt8},Ref{Cint}),
          f.ptr, value, status)
    fits_assert_ok(status[])
    unsafe_string(pointer(value))
end

function fits_file_mode(f::FITSFile)
    result = Ref{Cint}(0)
    status = Ref{Cint}(0)
    ccall(("ffflmd", libcfitsio), Cint, (Ptr{Void}, Ref{Cint}, Ref{Cint}),
          f.ptr, result, status)
    fits_assert_ok(status[])
    result[]
end

# -----------------------------------------------------------------------------
# header access functions

"""
    fits_get_hdrspace(f::FITSFile) -> (keysexist, morekeys)

Return the number of existing keywords (not counting the END keyword)
and the amount of space currently available for more keywords.
"""
function fits_get_hdrspace(f::FITSFile)
    keysexist = Ref{Cint}(0)
    morekeys = Ref{Cint}(0)
    status = Ref{Cint}(0)
    ccall((:ffghsp,libcfitsio), Cint,
        (Ptr{Void},Ref{Cint},Ref{Cint},Ref{Cint}),
        f.ptr, keysexist, morekeys, status)
    fits_assert_ok(status[])
    (keysexist[], morekeys[])
end

function fits_read_key_str(f::FITSFile, keyname::String)
    value = Vector{UInt8}(71)
    comment = Vector{UInt8}(71)
    status = Ref{Cint}(0)
    ccall((:ffgkys, libcfitsio), Cint,
          (Ptr{Void}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ref{Cint}),
          f.ptr, keyname, value, comment, status)
    fits_assert_ok(status[])
    unsafe_string(pointer(value)), unsafe_string(pointer(comment))
end

function fits_read_key_lng(f::FITSFile, keyname::String)
    value = Ref{Clong}(0)
    comment = Vector{UInt8}(71)
    status = Ref{Cint}(0)
    ccall((:ffgkyj, libcfitsio), Cint,
          (Ptr{Void}, Ptr{UInt8}, Ref{Clong}, Ptr{UInt8}, Ref{Cint}),
          f.ptr, keyname, value, comment, status)
    fits_assert_ok(status[])
    value[], unsafe_string(pointer(comment))
end

function fits_read_keys_lng(f::FITSFile, keyname::String,
                            nstart::Integer, nmax::Integer)
    value = Vector{Clong}(nmax - nstart + 1)
    nfound = Ref{Cint}(0)
    status = Ref{Cint}(0)
    ccall((:ffgknj, libcfitsio), Cint,
          (Ptr{Void}, Ptr{UInt8}, Cint, Cint, Ptr{Clong}, Ref{Cint}, Ref{Cint}),
          f.ptr, keyname, nstart, nmax, value, nfound, status)
    fits_assert_ok(status[])
    value, nfound[]
end

"""
    fits_read_keyword(f::FITSFile, keyname::String) -> (value, comment)

Return the specified keyword.
"""
function fits_read_keyword(f::FITSFile, keyname::String)
    value = Vector{UInt8}(71)
    comment = Vector{UInt8}(71)
    status = Ref{Cint}(0)
    ccall((:ffgkey,libcfitsio), Cint,
        (Ptr{Void},Ptr{UInt8},Ptr{UInt8},Ptr{UInt8},Ref{Cint}),
        f.ptr, keyname, value, comment, status)
    fits_assert_ok(status[])
    unsafe_string(pointer(value)), unsafe_string(pointer(comment))
end

"""
    fits_read_record(f::FITSFile, keynum::Int) -> String

Return the nth header record in the CHU. The first keyword in the
header is at `keynum = 1`.
"""
function fits_read_record(f::FITSFile, keynum::Integer)
    card = Vector{UInt8}(81)
    status = Ref{Cint}(0)
    ccall((:ffgrec,libcfitsio), Cint,
        (Ptr{Void},Cint,Ptr{UInt8},Ref{Cint}),
        f.ptr, keynum, card, status)
    fits_assert_ok(status[])
    unsafe_string(pointer(card))
end

"""
    fits_read_keyn(f::FITSFile, keynum::Int) -> (name, value, comment)

Return the nth header record in the CHU. The first keyword in the header is at `keynum = 1`.
"""
function fits_read_keyn(f::FITSFile, keynum::Integer)
    keyname = Vector{UInt8}(9)
    value = Vector{UInt8}(71)
    comment = Vector{UInt8}(71)
    status = Ref{Cint}(0)
    ccall((:ffgkyn,libcfitsio), Cint,
        (Ptr{Void},Cint,Ptr{UInt8},Ptr{UInt8},Ptr{UInt8},Ref{Cint}),
        f.ptr, keynum, keyname, value, comment, status)
    fits_assert_ok(status[])
    (unsafe_string(pointer(keyname)), unsafe_string(pointer(value)),
     unsafe_string(pointer(comment)))
end

"""
    fits_write_key(f::FITSFile, keyname::String, value, comment::String)

Write a keyword of the appropriate data type into the CHU.
"""
function fits_write_key(f::FITSFile, keyname::String,
                        value::Union{Real,String}, comment::String)
    fits_assert_isascii(keyname)
    fits_assert_isascii(comment)
    cvalue = isa(value,String) ?  value :
             isa(value,Bool) ? Cint[value] : reinterpret(UInt8, [value])
    status = Ref{Cint}(0)
    ccall((:ffpky,libcfitsio), Cint,
        (Ptr{Void},Cint,Ptr{UInt8},Ptr{UInt8},Ptr{UInt8},Ref{Cint}),
        f.ptr, cfitsio_typecode(typeof(value)), keyname,
        cvalue, comment, status)
    fits_assert_ok(status[])
end

function fits_write_date(f::FITSFile)
    status = Ref{Cint}(0)
    ccall((:ffpdat, libcfitsio), Cint, (Ptr{Void}, Ref{Cint}), f.ptr, status)
    fits_assert_ok(status[])
end

function fits_write_comment(f::FITSFile, comment::String)
    fits_assert_isascii(comment)
    status = Ref{Cint}(0)
    ccall((:ffpcom, libcfitsio), Cint, (Ptr{Void}, Ptr{UInt8}, Ref{Cint}),
          f.ptr, comment, status)
    fits_assert_ok(status[])
end

function fits_write_history(f::FITSFile, history::String)
    fits_assert_isascii(history)
    status = Ref{Cint}(0)
    ccall((:ffphis, libcfitsio), Cint, (Ptr{Void}, Ptr{UInt8}, Ref{Cint}),
          f.ptr, history, status)
    fits_assert_ok(status[])
end

# update key: if already present, update it, otherwise add it.
for (a,T,S) in (("ffukys", :String, :(Ptr{UInt8})),
                ("ffukyl", :Bool,        :Cint),
                ("ffukyj", :Integer,     :Int64))
    @eval begin
        function fits_update_key(f::FITSFile, key::String, value::$T,
                                 comment::Union{String, Ptr{Void}}=C_NULL)
            isa(value, String) && fits_assert_isascii(value)
            isa(comment, String) && fits_assert_isascii(comment)
            status = Ref{Cint}(0)
            ccall(($a, libcfitsio), Cint,
                  (Ptr{Void}, Ptr{UInt8}, $S, Ptr{UInt8}, Ref{Cint}),
                  f.ptr, key, value, comment, status)
            fits_assert_ok(status[])
        end
    end
end

function fits_update_key(f::FITSFile, key::String, value::AbstractFloat,
                         comment::Union{String, Ptr{Void}}=C_NULL)
    isa(comment, String) && fits_assert_isascii(comment)
    status = Ref{Cint}(0)
    ccall(("ffukyd", libcfitsio), Cint,
          (Ptr{Void}, Ptr{UInt8}, Cdouble, Cint, Ptr{UInt8}, Ref{Cint}),
          f.ptr, key, value, -15, comment, status)
    fits_assert_ok(status[])
end

function fits_update_key(f::FITSFile, key::String, value::Void,
                         comment::Union{String, Ptr{Void}}=C_NULL)
    isa(comment, String) && fits_assert_isascii(comment)
    status = Ref{Cint}(0)
    ccall(("ffukyu", libcfitsio), Cint,
          (Ptr{Void}, Ptr{UInt8}, Ptr{UInt8}, Ref{Cint}),
          f.ptr, key, comment, status)
    fits_assert_ok(status[])
end

"""
    fits_write_record(f::FITSFile, card::String)

Write a user specified keyword record into the CHU.
"""
function fits_write_record(f::FITSFile, card::String)
    fits_assert_isascii(card)
    status = Ref{Cint}(0)
    ccall((:ffprec,libcfitsio), Cint,
        (Ptr{Void},Ptr{UInt8},Ref{Cint}),
        f.ptr, card, status)
    fits_assert_ok(status[])
end

"""
    fits_delete_record(f::FITSFile, keynum::Int)

Delete the keyword record at the specified index.
"""
function fits_delete_record(f::FITSFile, keynum::Integer)
    status = Ref{Cint}(0)
    ccall((:ffdrec,libcfitsio), Cint,
        (Ptr{Void},Cint,Ref{Cint}),
        f.ptr, keynum, status)
    fits_assert_ok(status[])
end

"""
    fits_delete_key(f::FITSFile, keyname::String)

Delete the keyword named `keyname`.
"""
function fits_delete_key(f::FITSFile, keyname::String)
    status = Ref{Cint}(0)
    ccall((:ffdkey,libcfitsio), Cint,
        (Ptr{Void},Ptr{UInt8},Ref{Cint}),
        f.ptr, keyname, status)
    fits_assert_ok(status[])
end

"""
    fits_hdr2str(f::FITSFile, nocomments::Bool=false)

Return the header of the CHDU as a string. If `nocomments` is `true`, comment
cards are stripped from the output.
"""
function fits_hdr2str(f::FITSFile, nocomments::Bool=false)
    status = Ref{Cint}(0)
    header = Ref{Ptr{UInt8}}()
    nkeys = Ref{Cint}(0)
    ccall((:ffhdr2str, libcfitsio), Cint,
          (Ptr{Void}, Cint, Ptr{Ptr{UInt8}}, Cint,
           Ptr{Ptr{UInt8}}, Ref{Cint}, Ref{Cint}),
          f.ptr, nocomments, &C_NULL, 0, header, nkeys, status)
    result = unsafe_string(header[])

    # free header pointer allocated by cfitsio (result is a copy)
    ccall((:fffree, libcfitsio), Ref{Cint},
          (Ptr{UInt8}, Ref{Cint}),
          header[], status)
    fits_assert_ok(status[])
    result
end

# -----------------------------------------------------------------------------
# HDU info functions and moving the current HDU

function hdu_int_to_type(hdu_type_int)
    if hdu_type_int == 0
        return :image_hdu
    elseif hdu_type_int == 1
        return :ascii_table
    elseif hdu_type_int == 2
        return :binary_table
    end

    :unknown
end

"""
    fits_movabs_hdu(f::FITSFile, hduNum::Integer)

Change the current HDU to the value specified by `hduNum`, and return a symbol
describing the type of the HDU.

Possible symbols are: `image_hdu`, `ascii_table`, or `binary_table`.
The value of `hduNum` must range between 1 and the value returned by
[`fits_get_num_hdus`](@ref).
"""
function fits_movabs_hdu end

"""
    fits_movrel_hdu(f::FITSFile, hduNum::Integer)

Change the current HDU by moving forward or backward by `hduNum` HDUs
(positive means forward), and return the same as [`fits_movabs_hdu`](@ref).
"""
function fits_movrel_hdu end
for (a,b) in ((:fits_movabs_hdu,"ffmahd"),
              (:fits_movrel_hdu,"ffmrhd"))
    @eval begin
        function ($a)(f::FITSFile, hduNum::Integer)
            hdu_type = Ref{Cint}(0)
            status = Ref{Cint}(0)
            ccall(($b,libcfitsio), Cint,
                  (Ptr{Void}, Cint, Ref{Cint}, Ref{Cint}),
                  f.ptr, hduNum, hdu_type, status)
            fits_assert_ok(status[])
            hdu_int_to_type(hdu_type[])
        end
    end
end

"""
    fits_movnam_hdu(f::FITSFile, extname::String, extver::Integer=0,
                    hdu_type_int::Integer=-1)

Change the current HDU by moving to the (first) HDU which has the specified
extension type and EXTNAME and EXTVER keyword values (or HDUNAME and HDUVER keywords).

If `extver` is 0 (the default) then the EXTVER keyword is ignored and the first HDU
with a matching EXTNAME (or HDUNAME) keyword will be found. If `hdu_type_int`
is -1 (the default) only the extname and extver values will be used to locate the
correct extension. If no matching HDU is found in the file, the current HDU will
remain unchanged.
"""
function fits_movnam_hdu(f::FITSFile, extname::String, extver::Integer=0,
                         hdu_type::Integer=-1)
    status = Ref{Cint}(0)
    ccall((:ffmnhd,libcfitsio), Cint,
          (Ptr{Void}, Cint, Ptr{UInt8}, Cint, Ref{Cint}),
          f.ptr, hdu_type, extname, extver, status)
    fits_assert_ok(status[])
end

function fits_get_hdu_num(f::FITSFile)
    hdunum = Ref{Cint}(0)
    ccall((:ffghdn,libcfitsio), Cint,
          (Ptr{Void}, Ref{Cint}),
          f.ptr, hdunum)
    hdunum[]
end

function fits_get_hdu_type(f::FITSFile)
    hdutype = Ref{Cint}(0)
    status = Ref{Cint}(0)
    ccall((:ffghdt, libcfitsio), Cint,
          (Ptr{Void}, Ref{Cint}, Ref{Cint}),
          f.ptr, hdutype, status)
    fits_assert_ok(status[])
    hdu_int_to_type(hdutype[])
end

# -----------------------------------------------------------------------------
# image HDU functions

"""
    fits_get_img_size(f::FITSFile)

Get the dimensions of the image.
"""
function fits_get_img_size end
for (a, b) in ((:fits_get_img_type,      "ffgidt"),
               (:fits_get_img_equivtype, "ffgiet"),
               (:fits_get_img_dim,       "ffgidm"))
    @eval function ($a)(f::FITSFile)
        result = Ref{Cint}(0)
        status = Ref{Cint}(0)
        ccall(($b, libcfitsio), Cint, (Ptr{Void}, Ref{Cint}, Ref{Cint}),
              f.ptr, result, status)
        fits_assert_ok(status[])
        result[]
    end
end

"""
    fits_create_img(f::FITSFile, t::Type, naxes::Vector{Int})

Create a new primary array or IMAGE extension with a specified data type and size.
"""
function fits_create_img{T, S<:Integer}(f::FITSFile, ::Type{T},
                                        naxes::Vector{S})
    status = Ref{Cint}(0)
    ccall((:ffcrimll, libcfitsio), Cint,
          (Ptr{Void}, Cint, Cint, Ptr{Int64}, Ref{Cint}),
          f.ptr, bitpix_from_type(T), length(naxes),
          Vector{Int64}(naxes), status)
    fits_assert_ok(status[])
end

"""
    fits_write_pix(f::FITSFile, fpixel::Vector{Int}, nelements::Int, data::Array)

Write pixels from `data` into the FITS file.
"""
function fits_write_pix{S<:Integer,T}(f::FITSFile, fpixel::Vector{S},
                                      nelements::Integer, data::Array{T})
    status = Ref{Cint}(0)
    ccall((:ffppxll, libcfitsio), Cint,
          (Ptr{Void}, Cint, Ptr{Int64}, Int64, Ptr{Void}, Ref{Cint}),
          f.ptr, cfitsio_typecode(T), Vector{Int64}(fpixel),
          nelements, data, status)
    fits_assert_ok(status[])
end

function fits_write_pix(f::FITSFile, data::Array)
    fits_write_pix(f, ones(Int64, length(size(data))), length(data), data)
end

function fits_read_pix{S<:Integer,T}(f::FITSFile, fpixel::Vector{S},
                                     nelements::Int, nullval::T,
                                     data::Array{T})
    anynull = Ref{Cint}(0)
    status = Ref{Cint}(0)
    ccall((:ffgpxvll, libcfitsio), Cint,
          (Ptr{Void}, Cint, Ptr{Int64}, Int64, Ptr{Void}, Ptr{Void},
           Ref{Cint}, Ref{Cint}),
          f.ptr, cfitsio_typecode(T), Vector{Int64}(fpixel),
          nelements, &nullval, data, anynull, status)
    fits_assert_ok(status[])
    anynull[]
end

"""
    fits_read_pix(f::FITSFile, fpixel::Vector{Int}, nelements::Int, data::Array)

Read pixels from the FITS file into `data`.
"""
function fits_read_pix{S<:Integer,T}(f::FITSFile, fpixel::Vector{S},
                                     nelements::Int, data::Array{T})
    anynull = Ref{Cint}(0)
    status = Ref{Cint}(0)
    ccall((:ffgpxvll, libcfitsio), Cint,
          (Ptr{Void}, Cint, Ptr{Int64}, Int64, Ptr{Void}, Ptr{Void},
           Ref{Cint}, Ref{Cint}),
          f.ptr, cfitsio_typecode(T), Vector{Int64}(fpixel),
          nelements, C_NULL, data, anynull, status)
    fits_assert_ok(status[])
    anynull[]
end

function fits_read_pix(f::FITSFile, data::Array)
    fits_read_pix(f, ones(Int64,length(size(data))), length(data), data)
end

function fits_read_subset{S1<:Integer,S2<:Integer,S3<:Integer,T}(
             f::FITSFile, fpixel::Vector{S1}, lpixel::Vector{S2},
             inc::Vector{S3}, data::Array{T})
    anynull = Ref{Cint}(0)
    status = Ref{Cint}(0)
    ccall((:ffgsv, libcfitsio), Cint,
          (Ptr{Void}, Cint, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Ptr{Void},
           Ptr{Void}, Ref{Cint}, Ref{Cint}),
          f.ptr, cfitsio_typecode(T),
          Vector{Clong}(fpixel),
          Vector{Clong}(lpixel),
          Vector{Clong}(inc),
          C_NULL, data, anynull, status)
    fits_assert_ok(status[])
    anynull[]
end

function fits_copy_image_section(fin::FITSFile, fout::FITSFile,
                                 section::String)
    status = Ref{Cint}(0)
    ccall((:fits_copy_image_section,libcfitsio), Cint,
          (Ptr{Void}, Ptr{Void}, Ptr{UInt8}, Ref{Cint}),
          fin.ptr, fout.ptr, section, status)
    fits_assert_ok(status[])
end

# -----------------------------------------------------------------------------
# ASCII/binary table HDU functions

# The three fields are: ttype, tform, tunit (CFITSIO's terminology)
const ColumnDef = Tuple{String, String, String}

"""
    fits_create_binary_tbl(f::FITSFile, numrows::Integer, coldefs::Array{ColumnDef},
                           extname::String)

Append a new HDU containing a binary table. The meaning of the parameters is the same
as in a call to [`fits_create_ascii_tbl`](@ref).

In general, one should pick this function for creating tables in a new HDU,
as binary tables require less space on the disk and are more efficient to read and write.
(Moreover, a few datatypes are not supported in ASCII tables).
"""
function fits_create_binary_tbl end

"""
    fits_create_ascii_tbl(f::FITSFile, numrows::Integer, coldefs::Array{ColumnDef},
                          extname::String)

Append a new HDU containing an ASCII table.

The table will have `numrows` rows (this parameter can be set to zero), each
initialized with the default value. In order to create a table, the programmer
must specify the characteristics of each column. The columns are specified by the
`coldefs` variable, which is an array of tuples.
Each tuple must have three string fields:

1. The name of the column.
2. The data type and the repetition count. It must be a string made by a number
   (the repetition count) followed by a letter specifying the type (in the example
   above, `D` stands for `Float64`, `E` stands for `Float32`, `A` stands for `Char`).
   Refer to the CFITSIO documentation for more information about the syntax of this
   parameter.
3. The measure unit of this field. This is used only as a comment.

The value of `extname` sets the "extended name" of the table, i.e., a string
that in some situations can be used to refer to the HDU itself.

Note that, unlike for binary tables, CFITSIO puts some limitations to the
types that can be used in an ASCII table column. Refer to the CFITSIO manual
for further information.

See also [`fits_create_binary_tbl`](@ref) for a similar function which
creates binary tables.
"""
function fits_create_ascii_tbl end
for (a,b) in ((:fits_create_binary_tbl, 2),
              (:fits_create_ascii_tbl,  1))
    @eval begin
        function ($a)(f::FITSFile, numrows::Integer,
                      coldefs::Array{ColumnDef}, extname::String)

            # Ensure that extension name, column names and units are
            # ASCII, as these get written to the file. We don't check
            # need to check that tform is ASCII because presumably
            # cfitsio will thrown an appropriate error if it doesn't
            # recognize the tform string.
            fits_assert_isascii(extname)
            for coldef in coldefs
                fits_assert_isascii(coldef[1])
                fits_assert_isascii(coldef[3])
            end

            # get length and convert coldefs to three arrays of Ptr{Uint8}
            ntype = length(coldefs)
            ttype = [pointer(x[1]) for x in coldefs]
            tform = [pointer(x[2]) for x in coldefs]
            tunit = [pointer(x[3]) for x in coldefs]
            status = Ref{Cint}(0)

            ccall(("ffcrtb", libcfitsio), Cint,
                  (Ptr{Void}, Cint, Int64, Cint,
                   Ptr{Ptr{UInt8}}, Ptr{Ptr{UInt8}},
                   Ptr{Ptr{UInt8}}, Ptr{UInt8}, Ref{Cint}),
                  f.ptr, $b, numrows, ntype,
                  ttype, tform, tunit, extname,
                  status)
            fits_assert_ok(status[])
        end
    end
end

"""
    fits_get_num_hdus(f::FITSFile)

Return the number of HDUs in the file.
"""
function fits_get_num_hdus end
for (a,b,T) in ((:fits_get_num_cols,  "ffgncl",  :Cint),
                (:fits_get_num_hdus,  "ffthdu",  :Cint),
                (:fits_get_rowsize,   "ffgrsz",  :Clong))
    @eval begin
        function ($a)(f::FITSFile)
            result = Ref{$T}(0)
            status = Ref{Cint}(0)
            ccall(($b,libcfitsio), Cint,
                  (Ptr{Void}, Ref{$T}, Ref{Cint}),
                  f.ptr, result, status)
            fits_assert_ok(status[])
            result[]
        end
    end
end

function fits_get_colnum(f::FITSFile, tmplt::String; case_sensitive::Bool=true)
    result = Ref{Cint}(0)
    status = Ref{Cint}(0)

    # Second argument is case-sensitivity of search: 0 = case-insensitive
    #                                                1 = case-sensitive
    ccall(("ffgcno", libcfitsio), Cint,
          (Ptr{Void}, Cint, Ptr{UInt8}, Ref{Cint}, Ref{Cint}),
          f.ptr, case_sensitive, tmplt, result, status)
    fits_assert_ok(status[])
    return result[]
end

# The following block are all functions that have separate variants for Clong
# and 64-bit integers in cfitsio. Rather than providing both of these, we
# provide only one according to the native integer type on the platform.
if  promote_type(Int, Clong) == Clong
    T = Clong
    ffgtdm = "ffgtdm"
    ffgnrw = "ffgnrw"
    ffptdm = "ffptdm"
    ffgtcl = "ffgtcl"
    ffeqty = "ffeqty"
    ffgdes = "ffgdes"
    ffgisz = "ffgisz"
else
    T = Int64
    ffgtdm = "ffgtdmll"
    ffgnrw = "ffgnrwll"
    ffptdm = "ffptdmll"
    ffgtcl = "ffgtclll"
    ffeqty = "ffeqtyll"
    ffgdes = "ffgdesll"
    ffgisz = "ffgiszll"
end

"""
    fits_get_coltype(f::FITSFile, colnum::Integer)

Provided that the current HDU contains either an ASCII or binary table, return
information about the column at position `colnum` (counting from 1).

Return is a tuple containing

- `typecode`: CFITSIO integer type code of the column.
- `repcount`: Repetition count for the column.
- `width`: Width of an individual element.
"""
function fits_get_coltype end

@eval begin
    function fits_get_coltype(ff::FITSFile, colnum::Integer)
        typecode = Ref{Cint}(0)
        repcnt = Ref{$T}(0)
        width = Ref{$T}(0)
        status = Ref{Cint}(0)
        ccall(($ffgtcl,libcfitsio), Cint,
              (Ptr{Void}, Cint, Ref{Cint}, Ref{$T}, Ref{$T}, Ref{Cint}),
              ff.ptr, colnum, typecode, repcnt, width, status)
        fits_assert_ok(status[])
        return Int(typecode[]), Int(repcnt[]), Int(width[])
    end

    function fits_get_eqcoltype(ff::FITSFile, colnum::Integer)
        typecode = Ref{Cint}(0)
        repcnt = Ref{$T}(0)
        width = Ref{$T}(0)
        status = Ref{Cint}(0)
        ccall(($ffeqty,libcfitsio), Cint,
              (Ptr{Void}, Cint, Ref{Cint}, Ref{$T}, Ref{$T}, Ref{Cint}),
              ff.ptr, colnum, typecode, repcnt, width, status)
        fits_assert_ok(status[])
        return Int(typecode[]), Int(repcnt[]), Int(width[])
    end

    function fits_get_img_size(f::FITSFile)
        ndim = fits_get_img_dim(f)
        naxes = Vector{$T}(ndim)
        status = Ref{Cint}(0)
        ccall(($ffgisz, libcfitsio), Cint,
              (Ptr{Void}, Cint, Ptr{$T}, Ref{Cint}),
              f.ptr, ndim, naxes, status)
        fits_assert_ok(status[])
        naxes
    end

    function fits_get_num_rows(f::FITSFile)
        result = Ref{$T}(0)
        status = Ref{Cint}(0)
        ccall(($ffgnrw, libcfitsio), Cint,
              (Ptr{Void}, Ref{$T}, Ref{Cint}),
              f.ptr, result, status)
        fits_assert_ok(status[])
        return Int(result[])
    end

    # `fits_read_tdim` returns the dimensions of a table column in a
    # binary table. Normally this information is given by the TDIMn
    # keyword, but if this keyword is not present then this routine
    # returns `[r]` with `r` equals to the repeat count in the TFORM
    # keyword.
    function fits_read_tdim(ff::FITSFile, colnum::Integer)
        naxes = Vector{$T}(99)  # 99 is the maximum allowed number of axes
        naxis = Ref{Cint}(0)
        status = Ref{Cint}(0)
        ccall(($ffgtdm,libcfitsio), Cint,
              (Ptr{Void}, Cint, Cint, Ref{Cint}, Ptr{$T}, Ref{Cint}),
              ff.ptr, colnum, length(naxes), naxis, naxes, status)
        fits_assert_ok(status[])
        return naxes[1:naxis[]]
    end

    function fits_write_tdim(ff::FITSFile, colnum::Integer,
                                 naxes::Array{$T})
        status = Ref{Cint}(0)
        ccall(($ffptdm, libcfitsio), Cint,
              (Ptr{Void}, Cint, Cint, Ptr{$T}, Ref{Cint}),
              ff.ptr, colnum, length(naxes), naxes, status)
        fits_assert_ok(status[])
    end

    function fits_read_descript(f::FITSFile, colnum::Integer, rownum::Integer)
        repeat = Ref{$T}(0)
        offset = Ref{$T}(0)
        status = Ref{Cint}(0)
        ccall(($ffgdes, libcfitsio), Cint,
              (Ptr{Void}, Cint, Int64, Ref{$T}, Ref{$T}, Ref{Cint}),
              f.ptr, colnum, rownum, repeat, offset, status)
        fits_assert_ok(status[])
        return Int(repeat[]), Int(offset[])
    end
end

"""
    fits_read_col(f, colnum, firstrow, firstelem, data)

Read data from one column of an ASCII/binary table and convert the data into the
specified type `T`.

### Arguments ###

* `f::FITSFile`: the file to be read.
* `colnum::Integer`: the column number, where the value of the first column is `1`.
* `firstrow::Integer`: the elements to be read start from this row.
* `firstelem::Integer`: specifies which is the first element to be read, when each
  cell contains more than one element (i.e., the "repetition count" of the field is
  greater than one).
* `data::Array`: at the end of the call, this will be filled with the elements read
  from the column. The length of the array gives the overall number of elements.
"""
function fits_read_col(f::FITSFile,
                       colnum::Integer,
                       firstrow::Integer,
                       firstelem::Integer,
                       data::Array{String})

    # get width: number of characters in each string
    typecode, repcount, width = fits_get_eqcoltype(f, colnum)

    # ensure that data are strings, otherwise cfitsio will try to write
    # formatted strings, which have widths given by fits_get_col_display_width
    # not by the repeat value from fits_get_coltype.
    abs(typecode) == 16 || error("not a string column")

    # create an array of character buffers of the correct width
    buffers = [Vector{UInt8}(width) for i in 1:length(data)]

    # Call the CFITSIO function
    anynull = Ref{Cint}(0)
    status = Ref{Cint}(0)
    ccall((:ffgcvs, libcfitsio), Cint,
          (Ptr{Void}, Cint, Int64, Int64, Int64,
           Ptr{UInt8}, Ptr{Ptr{UInt8}}, Ref{Cint}, Ref{Cint}),
          f.ptr, colnum, firstrow, firstelem, length(data),
          " ", buffers, anynull, status)
    fits_assert_ok(status[])

    # Create strings out of the buffers, terminating at null characters.
    # Note that `String(x)` does not copy the buffer x.
    for i in 1:length(data)
        zeropos = search(buffers[i], 0x00)
        data[i] = (zeropos >= 1) ? String(buffers[i][1:(zeropos-1)]) :
                                   String(buffers[i])
    end
end

function fits_read_col{T}(f::FITSFile,
                          colnum::Integer,
                          firstrow::Integer,
                          firstelem::Integer,
                          data::Array{T})
    anynull = Ref{Cint}(0)
    status = Ref{Cint}(0)
    ccall((:ffgcv,libcfitsio), Cint,
          (Ptr{Void}, Cint, Cint, Int64, Int64, Int64,
           Ptr{Void}, Ptr{Void}, Ref{Cint}, Ref{Cint}),
          f.ptr, cfitsio_typecode(T), colnum,
          firstrow, firstelem, length(data),
          C_NULL, data, anynull, status)
    fits_assert_ok(status[])
end

"""
    fits_write_col(f, colnum, firstrow, firstelem, data)

Write some data in one column of a ASCII/binary table.

If there is no room for the elements, new rows will be created. (It is therefore
useless to call [`fits_insert_rows`](@ref) if you only need to *append* elements
to the end of a table.)

* `f::FITSFile`: the file in which data will be written.
* `colnum::Integer`: the column number, where the value of the first column is `1`.
* `firstrow::Integer`: the data wil be written from this row onwards.
* `firstelem::Integer`: specifies the position in the row where the first element
  will be written.
* `data::Array`: contains the elements that are to be written to the column of the table.
"""
function fits_write_col(f::FITSFile,
                        colnum::Integer,
                        firstrow::Integer,
                        firstelem::Integer,
                        data::Array{String})
    for el in data; fits_assert_isascii(el); end
    status = Ref{Cint}(0)
    ccall((:ffpcls, libcfitsio), Cint,
          (Ptr{Void}, Cint, Int64, Int64, Int64,
           Ptr{Ptr{UInt8}}, Ref{Cint}),
          f.ptr, colnum, firstrow, firstelem, length(data),
          data, status)
    fits_assert_ok(status[])
end

function fits_write_col{T}(f::FITSFile,
                           colnum::Integer,
                           firstrow::Integer,
                           firstelem::Integer,
                           data::Array{T})
    status = Ref{Cint}(0)
    ccall((:ffpcl, libcfitsio), Cint,
          (Ptr{Void}, Cint, Cint, Int64, Int64, Int64,
           Ptr{Void}, Ref{Cint}),
          f.ptr, cfitsio_typecode(T), colnum,
          firstrow, firstelem, length(data),
          data, status)
    fits_assert_ok(status[])
end

"""
    fits_insert_rows(f::FITSFile, firstrow::Integer, nrows::Integer)

Insert a number of rows equal to `nrows` after the row number `firstrow`.

The elements in each row are initialized to their default value: you can
modify them later using [`fits_write_col`](@ref).

Since the first row is at position 1, in order to insert rows *before*
the first one `firstrow` must be equal to zero.
"""
function fits_insert_rows end

"""
    fits_delete_rows(f::FITSFile, firstrow::integer, nrows::Integer)

Delete `nrows` rows, starting from the one at position `firstrow`. The index of
the first row is 1.
"""
function fits_delete_rows end

for (a,b) in ((:fits_insert_rows, "ffirow"),
              (:fits_delete_rows, "ffdrow"))
    @eval begin
        function ($a)(f::FITSFile, firstrow::Integer, nrows::Integer)
            status = Ref{Cint}(0)
            ccall(($b,libcfitsio), Cint,
                  (Ptr{Void}, Int64, Int64, Ref{Cint}),
                  f.ptr, firstrow, nrows, status)
            fits_assert_ok(status[])
        end
    end
end

end # module


text:
using Test
using Dates, DataFrames, CSV, YAML
using Streamfall

DATA_PATH = joinpath(@__DIR__, "data/hymod")

# Load and generate stream network
network = YAML.load_file(joinpath(DATA_PATH, "hymod_network.yml"))
sn = create_network("HyMod Network", network)

@testset "Single node network" begin
    inlets, outlets = find_inlets_and_outlets(sn)

    @test length(inlets) == 0
    @test length(outlets) == 1
end

@testset "Saving a network spec" begin
    spec = Streamfall.extract_network_spec(sn)
    @test spec isa Dict
    @test haskey(spec, "leaf_river")

    spec["02473000"] = spec["leaf_river"]

    tmp_sn = create_network("Test", spec)

    tmp_fn, io = mktemp()
    @test Streamfall.save_network_spec(tmp_sn, tmp_fn) isa Any

    # Ensure reloaded network is usable
    tmp_network = YAML.load_file(tmp_fn)
    new_sn = create_network("Test", tmp_network)

    @test new_sn isa Streamfall.StreamfallNetwork

    # There is currently a bug in YAML.jl which parses string integers as Ints
    # which causes issues when the value starts with a 0
    # (incorrectly parses as an octal integer)
    # If this test does not throw an error, then it means the issue has been resolved
    # and can be removed after checking.
    # see: https://github.com/JuliaData/YAML.jl/pull/45
    @test_throws ErrorException new_sn["02473000"]
end

@testset "Reloading a network spec (HyMod)" begin
    spec::Dict{String, Union{Dict{String, Any}, Any}} = Streamfall.extract_network_spec(sn)
    @test spec isa Dict
    @test haskey(spec, "leaf_river")

    tmp_sn = create_network("Test", spec)

    tmp_fn, io = mktemp()
    @test Streamfall.save_network_spec(tmp_sn, tmp_fn) isa Any

    # Ensure reloaded network is usable
    tmp_network = YAML.load_file(tmp_fn)
    new_sn = create_network("Test", tmp_network)

    @test new_sn isa Streamfall.StreamfallNetwork

    # Load climate data
    date_format = "YYYY-mm-dd"
    obs_data = CSV.File(joinpath(DATA_PATH, "leaf_river_data.csv"),
                        comment="#",
                        dateformat=date_format) |> DataFrame

    hist_streamflow = obs_data[:, "leaf_river_outflow"]
    climate_data = obs_data[:, ["Date", "leaf_river_P", "leaf_river_ET"]]
    climate = Climate(climate_data, "_P", "_ET")

    # Ensure reloaded spec still calibrates
    metric = (obs, sim) -> 1.0 - Streamfall.NNSE(obs, sim)
    @test calibrate!(new_sn, 1, climate, hist_streamflow; metric=metric, MaxTime=10) isa Any
end

@testset "Recursing IHACRESNode upstream" begin

    begin
        include("../examples/run_nodes.jl")
        # Ensure example does not error out

        reset!(sn)
        run_basin!(sn, climate; extraction=hist_dam_releases)
        @test Streamfall.RMSE(n_data, sn[dam_id].level) == 0.0

        reset!(sn)
        Streamfall.run_node!(sn, dam_id, climate; extraction=hist_dam_releases)
        @test Streamfall.RMSE(n_data, sn[dam_id].level) == 0.0

        # Ensure example results haven't changed much...
        # @test nse_score >= 0.95 && nse_score < 1.0 && rmse < 2.0
    end
end


text:
function ldlt!(a::Matrix{Float64}, d = Vector{Float64}; small=1.0e-10)
  n = size(a, 1)
  for k in 1:n-1
    d[1] = a[1,1]
    if abs(a[k,k]) > small
      for i in k+1:n
        x = a[i,k] / a[k,k]
        for j = k+1:n
          a[i,j] -= a[k,j] * x
        end
        d[i] = a[i,i]
      end
    else
      println("Zero pivot found in row $k")
    end
  end
end

function ldlt(a::Matrix{Float64}, d = Vector{Float64}; small=1.0e-10)
  at = deepcopy(a)
  dt = deepcopy(d)
  n = size(at, 1)
  for k in 1:n-1
    dt[1] = at[1,1]
    if abs(at[k,k]) > small
      for i in k+1:n
        x = at[i,k] / at[k,k]
        for j = k+1:n
          at[i,j] -= at[k,j] * x
        end
        dt[i] = at[i,i]
      end
    else
      println("Zero pivot found in row $k")
    end
  end
  (at, dt)
end

export
  ldlt!,
  ldlt


text:
"""
Solve MPC controller

    #Parameters

    -A MPC regulator method with x_t, x_ref and u_ref values.

    #Example

    u, x = solve_mpc(controller)

"""

function solve_mpc(controller)
    if controller.method == "l"
        # The variable are retrieved from the model
        x = controller.formulation[:x]
        x_tilde = controller.formulation[:x_tilde]
        x_ref = controller.formulation[:x_ref]
        u = controller.formulation[:u]
        u_tilde = controller.formulation[:u_tilde]
        u_ref = controller.formulation[:u_ref]
        x_ref_update = controller.x_ref
        u_ref_update = controller.u_ref

        # The JuMP model is updated
        JuMP.fix.(x[:,1], controller.x_t[:,1]; force = true)
        for k in 1 : 1 : size(u,2)
            JuMP.fix.(u_ref[:,k], u_ref_update[:,k]; force = true)
        end
        for k in 1 : 1 : size(x,2)
            JuMP.fix.(x_ref[:,k], x_ref_update[:,k]; force = true)
        end

        # The model is optimized
        JuMP.optimize!(controller.formulation)

        # The optimised values are retrieved
        u_optimised = JuMP.value.(u[:,:])
        u_tilde_optimised = JuMP.value.(u_tilde[:,:])
        x_optimised = JuMP.value.(x[:,:])
        x_tilde_optimised = JuMP.value.(x_tilde[:,:])
        obj_optimised = JuMP.objective_value(controller.formulation)
    end
    return u_optimised, x_optimised, u_tilde_optimised, x_tilde_optimised, obj_optimised
end


text:
@testset "Hat basis" begin

using ValidatedNumerics
using InvariantMeasures
using InvariantMeasures: HatFunction, HatFunctionOnTorus, IntervalOnTorus, nonzero_on

f = HatFunction(1., 2, 3)
@test f(1.5) == 0.5
@test f(1..1.5) == 0..0.5

f = HatFunctionOnTorus(0.125, 0.25, 0.375)
x = IntervalOnTorus(0.375..1.1875)
@test f(x) == 0..0.5

@test f(IntervalOnTorus(0..1)) == 0..1

f = HatFunctionOnTorus(0.125, 0.25, 0.375)
x = IntervalOnTorus(3.1875..3.25)
@test f(x) == 0.5..1

f = HatFunctionOnTorus(0, 0.125, 0.25)
x = IntervalOnTorus(0..0.0625)
@test f(x) == 0..0.5

f = HatFunctionOnTorus(0.875, 0, 0.125)
x = IntervalOnTorus(0..0.0625)
@test f(x) == 0.5..1

f = HatFunctionOnTorus(0.875, 0, 0.125)
x = IntervalOnTorus(0.9375..1)
@test f(x) == 0.5..1

B = Hat(4)
@test nonzero_on(B, (0.1..0.3, NaN)) == (1,3)
@test nonzero_on(B, (0..1, NaN)) == (1,4)
@test nonzero_on(B, (0.3..0.31, NaN)) == (2,3)

@test is_refinement(Hat(8), Hat(4))
@test is_refinement(Hat(8), Hat(8))
@test !is_refinement(Hat(4), Hat(8))

end


text:
#=
Created on August 2019
Last update: Tuesday 01 October 2019

@author: Michiel Stock
michielfmstock@gmail.com

Methods for fitting a utility matrix to an observed interaction network.
=#

using Zygote, Optim

# regularization methods
abstract type Regularization end

struct L1 <: Regularization end
struct L2 <: Regularization end
struct ElasticNet <: Regularization
    α::Real
    function ElasticNet(α::Real=0.5)
        @assert α ≥ 0 && α ≤ 1
        return new(α)
    end
end
struct Exp <: Regularization end
struct Entropic <: Regularization
    offset::Real
end

Entropic() = Entropic(0.0)

r(x, reg::L1) = sum(abs.(x))
r(x, reg::L2) = sum(x.^2)
r(x, reg::ElasticNet) = reg.α * sum(abs.(x)) + (1 - reg.α) * sum(x.^2)
r(x, reg::Exp) = sum(exp.(x))
r(x, reg::Entropic) = sum(x .* log.(x .+ reg.offset))

function _fitM(P::AbstractMatrix; γ=1e-2, reg::Regularization=L2())

    l = M -> relative_entropy(P, optimaltransport(M)) + γ * r(M, reg)
    dM = similar(P)
    ∇l!(dM, M) = dM .= l'(M)
    return Optim.minimizer(optimize(l, ∇l!, zero(P), BFGS()))
end

function _fitM(P::AbstractMatrix, a::AbstractVector; γ=1e-2, reg::Regularization=L2())

    l = M -> relative_entropy(P, optimaltransport(M, a)) + γ * r(M, reg)
    dM = similar(P)
    ∇l!(dM, M) = dM .= l'(M)
    return Optim.minimizer(optimize(l, ∇l!, zero(P), BFGS()))
end

function _fitM(P, a::AbstractVector, b::AbstractVector; γ=1e-2,
                                    reg::Regularization=L2(),
                                    maxitter::Int=100, ϵ::Real=1e-5)

    l = M -> relative_entropy(P, optimaltransport(M, a, b; maxitter=maxitter,
                                        ϵ=ϵ)) + γ * r(M, reg)
    dM = similar(P)
    ∇l!(dM, M) = dM .= l'(M)
    return Optim.minimizer(optimize(l, ∇l!, zero(P), BFGS()))
end

function _fitMnograd(P, a::AbstractVector, b::AbstractVector; γ=1e-2,
                                    reg::Regularization=L2(),
                                    maxitter::Int=100, ϵ::Real=1e-5)

    l = M -> relative_entropy(P, optimaltransport(M, a, b; maxitter=maxitter,
                                        ϵ=ϵ)) + γ * r(M, reg)
    dM = similar(P)
    return Optim.minimizer(optimize(l, zero(P), BFGS()))
end

"""
    fitM(P::AbstractMatrix; fix_a::Bool, fix_b::Bool,
                    γ=1e-2,
                    reg::Regularization=L2(),
                    maxitter::Int=100, ϵ::Real=1e-5)

Fit a utility matrix `M` to an observed interaction frequency matrix `P`.
The parameters `fix_a` and `fix_b` determine whether the marginals should be
fixed.
"""
function fitM(P::AbstractMatrix; fix_a::Bool=true, fix_b::Bool=true,
                    γ=1e-2,
                    reg::Regularization=L2(),
                    maxitter::Int=100, ϵ::Real=1e-5)
    a, b = marginals(P)
    if fix_a && fix_b
        return _fitM(P, a, b, reg=reg, γ=γ, maxitter=maxitter, ϵ=ϵ)
    elseif fix_a
        return _fitM(P, a, reg=reg, γ=γ)
    elseif fix_b
        return _fitM(P', b, reg=reg, γ=γ)'
    else
        return _fitM(P, reg=reg, γ=γ)
    end
end

# loss over several observed frequency matrices
function global_loss(M, Ps::Array{M, 1} where M <: AbstractMatrix,
                fix_a::Bool, fix_b::Bool;
                γ=1e-2,
                reg::Regularization=L2(),
                maxitter::Int=10,
                ϵ::Real=1e-5)
    l = 0.0
    for P in Ps
        a, b = marginals(P)
        Q = optimaltransport(M, fix_a ? a : nothing,
                                fix_b ? b : nothing,
                                ϵ=ϵ, maxitter=maxitter)
        l += relative_entropy(P, Q)
    end
    l += γ * r(M, reg)
    return l
end

"""
    parallel_cross_entropy(M, A, B, Ps; maxitter::Integer=20)

Compute the cross-entropy compared to optimal transport on a series of interaction
distributions.
"""
function parallel_cross_entropy(M::AbstractMatrix,
            Ps::Array{M, 1} where M <: AbstractMatrix;
            A::AbstractMatrix,
            B::AbstractMatrix,
            maxitter::Integer=20)
    K = exp.(M)
    # perform Sinkhorn iterations
    U, V = A, B
    for itter in 1:maxitter
        U = A ./ (K * V)  # fix row marginals
        V = B ./ (K' * U)  # fix column marginals
    end
    # compute cross-entropy
    c = 0
    for Pₒ in Ps
        c -= sum(Pₒ .* M)
    end
    c -= sum(A .* log.(U))
    c -= sum(B .* log.(V))
    return c
end

function parallel_cross_entropy_A(M::AbstractMatrix,
            Ps::Array{M, 1} where M <: AbstractMatrix;
            A::AbstractMatrix)
    K = exp.(M)
    U = A ./ sum(K, dims=2)  # fix row marginals
    # compute cross-entropy
    c = 0.0
    for Pₒ in Ps
        c -= sum(Pₒ .* M)
    end
    c -= sum(A .* log.(U))
    return c
end

function parallel_cross_entropy_B(M::AbstractMatrix,
            Ps::Array{M, 1} where M <: AbstractMatrix;
            B::AbstractMatrix)
    K = exp.(M)
    V = B ./ sum(K', dims=2)  # fix row marginals
    # compute cross-entropy
    c = 0.0
    for Pₒ in Ps
        c -= sum(Pₒ .* M)
    end
    c -= sum(B .* log.(V))
    return c
end

function parallel_cross_entropy_free(M::AbstractMatrix,
            Ps::Array{M, 1} where M <: AbstractMatrix)
    logδ = exp.(M) |> sum |> log
    c = 0.0
    for Pₒ in Ps
        c -= sum(Pₒ .* M) - logδ
    end
    return c
end

"""
    fitM(Ps::Array{M, 1} where M <: AbstractMatrix,
                    fix_a::Bool, fix_b::Bool,
                    γ=1e-2,
                    reg::Regularization=L2(),
                    maxitter::Int=100, ϵ::Real=1e-5)

Fit a utility matrix `M` to a collection of observed probability matrices `Ps=[Pi]`.
"""
function fitM(Ps::Array{M, 1} where M <: AbstractMatrix; fix_a=true, fix_b=true,
                γ=1e-2,
                reg::Regularization=L2(),
                maxitter::Int=20)
    # compute marginals
    A, B = marginals(Ps)
    if fix_a && fix_b
        l = M -> parallel_cross_entropy(M, Ps, A=A, B=B, maxitter=maxitter) + γ * r(M, reg)
    elseif fix_a
        l = M -> parallel_cross_entropy_A(M, Ps, A=A) + γ * r(M, reg)
    elseif fix_b
        l = M -> parallel_cross_entropy_B(M, Ps, B=B) + γ * r(M, reg)
    else
        l = M -> parallel_cross_entropy_free(M, Ps) + γ * r(M, reg)
    end
    dM = similar(first(Ps))
    ∇l!(dM, M) = dM .= l'(M)
    return Optim.minimizer(optimize(l, ∇l!, 0.01randn(size(first(Ps))), BFGS()))
end

using StatsBase: sample

"""
    fitMsgd(Ps::Array{M, 1} where M <: AbstractMatrix;
                    fix_a::Bool, fix_b::Bool,
                    γ=1e-2,
                    reg::Regularization=L2(),
                    maxitter::Int=20,
                    ϵ::Real=1e-4,
                    α::Real=0.1,
                    β::Real=0.8,
                    nsteps::Real=100)

Fit a utility matrix `M` to a collection of observed probability matrices
`Ps=[Pi]` using stochastic gradient descent with momentum.
The parameters `fix_a` and `fix_b` determine whether the marginals should be
fixed.
"""
function fitMsgd(Ps::Array{M, 1} where M <: AbstractMatrix;
                fix_a::Bool, fix_b::Bool,
                γ=1e-2,
                reg::Regularization=L2(),
                maxitter::Int=20,
                ϵ::Real=1e-4,
                α::Real=0.1,
                β::Real=0.8,
                nsteps::Real=100)  # momentum parameter)
    global_loss = M -> global_loss(M, Ps, fix_a, fix_b, γ=γ, reg=reg,
                                    maxitter=maxitter, ϵ=ϵ)
    losses = zeros(nsteps)
    ΔM = similar(first(Ps))
    fill!(ΔM, 0.0)
    M = zeros(size(ΔM)...)
    for i in 1:nsteps
        losses[i] = global_loss(M)
        for rep in 1:length(Ps)
            P = sample(Ps)
            a, b = marginals(P)
            l = M -> relative_entropy(P,
                    optimaltransport(M, fix_a ? a : nothing,
                            fix_b ? b : nothing,
                            maxitter=maxitter, ϵ=ϵ)) + (γ / length(Ps)) * r(M, reg)
            ΔM .*= 1 - β
            ΔM .+= α * l'(M)
            M .-= ΔM
        end
    end
    return M, losses
end

max_stars_repo_name: PoisotLab/OTSIN

max_issues_repo_name: PoisotLab/OTSIN

max_forks_repo_name: PoisotLab/OTSIN


text:
using VennABERS
using Test

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

max_stars_repo_name: PatrikDurdevic/VennABERS

max_issues_repo_name: PatrikDurdevic/VennABERS

max_forks_repo_name: PatrikDurdevic/VennABERS


text:
"""
    ncatt_get(f::NCfiles, key)
    ncatt_put(f::AbstractString, attrib = Dict())
    ncatt_del(f::AbstractString, keys::Vector{<:AbstractString})

add or delete global attributes
"""
function ncatt_put(f::NCfiles, atts = Dict())
    nc_open(f, "a") do nc
        names = keys(atts) |> collect
        vals = values(atts) |> collect
        for i = 1:length(atts)
            nc.attrib[names[i]] = vals[i]
        end
    end
    # close(nc)
end

function ncatt_get(f::NCfiles)
    nc_open(f) do nc
        nc.attrib |> collect
    end
end

function ncatt_get(f::NCfiles, keys::Vector)
    nc = nc_open(f)
    res = map(key -> begin
        if haskey(nc.attrib, key)
            nc.attrib[key]
        else
            @warn("`$key` not exits in file '$(basename(f))'")
            nothing
        end
    end, keys)
    nc_close(nc)
    # filter(not_empty, res)
    res
end

function ncatt_get(f::NCfiles, key::AbstractString)
    ncatt_get(f, [key])[1]
end

function ncatt_del(f::NCfiles, keys::Vector{<:AbstractString})
    nc_open(f, "a") do nc
        for i = 1:length(keys)
            delete!(nc.attrib, keys[i])
        end
    end
end

function ncatt_del(f::NCfiles)
    keys = map(x -> x.first, ncatt_get(f))
    ncatt_del(f, keys)
end


text:
using BinDeps

@BinDeps.setup

deps = [
  libkyotocabinet = library_dependency("libkyotocabinet", aliases = ["libkyotocabinet", "libkyotocabinet16"])
]

@osx_only begin
  using Homebrew
  provides(Homebrew.HB, "kyotocabinet", libkyotocabinet, os = :Darwin )
end

# System Package Managers
provides(AptGet, {
  "libkyotocabinet-dev" => libkyotocabinet
})

const lib_kc_ver = "1.2.76"

provides(Sources, {
  URI("http://fallabs.com/kyotocabinet/pkg/kyotocabinet-$(lib_kc_ver).tar.gz") => libkyotocabinet
})

builddir = joinpath(BinDeps.depsdir(libkyotocabinet), "usr")
srcdir = joinpath(BinDeps.depsdir(libkyotocabinet), "src", "kyotocabinet-$(lib_kc_ver)")
provides(SimpleBuild,
  (@build_steps begin
    GetSources(libkyotocabinet)
    CreateDirectory(builddir)
    @build_steps begin
      ChangeDirectory(srcdir)
      (@build_steps begin
        `./configure --prefix=$(builddir)`
        `make`
        `make install`
      end)
    end
  end), libkyotocabinet, os = :Unix)

@BinDeps.install [:libkyotocabinet => :libkyotocabinet]