text:
const COMPLEX_CHAR = 'K'

primitive type EdgeMotif 8 end

Base.convert( ::Type{EdgeMotif}, motif::UInt8 ) = reinterpret(EdgeMotif, motif )
Base.convert( ::Type{UInt8}, motif::EdgeMotif ) = reinterpret(UInt8, motif )

# OBLIGATE MOTIFS (provided the cognate nodes both have reads):

const TXST_MOTIF = convert(EdgeMotif, 0b000)
const TXEN_MOTIF = convert(EdgeMotif, 0b001)

const ALTF_MOTIF = convert(EdgeMotif, 0b010)
const ALTL_MOTIF = convert(EdgeMotif, 0b011)

# REQUIRES SPANNING EDGE:

const RETI_MOTIF = convert(EdgeMotif, 0b100)
const SKIP_MOTIF = convert(EdgeMotif, 0b101)

const ALTD_MOTIF = convert(EdgeMotif, 0b110)
const ALTA_MOTIF = convert(EdgeMotif, 0b111)

# NO MOTIF
const NONE_MOTIF = convert(EdgeMotif, 0b1000)

const MOTIF_STRING = fill( "", 9 )
      MOTIF_STRING[ UInt8(TXST_MOTIF) + 1 ] = "TS"
      MOTIF_STRING[ UInt8(TXEN_MOTIF) + 1 ] = "TE"
      MOTIF_STRING[ UInt8(ALTF_MOTIF) + 1 ] = "AF"
      MOTIF_STRING[ UInt8(ALTL_MOTIF) + 1 ] = "AL"
      MOTIF_STRING[ UInt8(RETI_MOTIF) + 1 ] = "RI"
      MOTIF_STRING[ UInt8(SKIP_MOTIF) + 1 ] = "CE"
      MOTIF_STRING[ UInt8(ALTD_MOTIF) + 1 ] = "AD"
      MOTIF_STRING[ UInt8(ALTA_MOTIF) + 1 ] = "AA"
      MOTIF_STRING[ UInt8(NONE_MOTIF) + 1 ] = "NA"

const MOTIF_TABLE = fill(NONE_MOTIF, 2^6 )
      # Alt TxStart  SL SL
      MOTIF_TABLE[ 0b000000 + 1 ] = TXST_MOTIF
      # Alt TxStart  LS SL
      MOTIF_TABLE[ 0b010000 + 1 ] = TXST_MOTIF

      # Alt PolyA    RS RS
      MOTIF_TABLE[ 0b011011 + 1 ] = TXEN_MOTIF
      # Alt PolyA    RS SR
      MOTIF_TABLE[ 0b011001 + 1 ] = TXEN_MOTIF

      # Alt FirstEx  LS LL
      MOTIF_TABLE[ 0b010101 + 1 ] = ALTF_MOTIF
      # Alt FirstEx  LS LR
      MOTIF_TABLE[ 0b010100 + 1 ] = ALTF_MOTIF
      # Alt FirstEx  LS LS
      MOTIF_TABLE[ 0b010010 + 1 ] = ALTF_MOTIF

      # Alt LastEx   RR SR
      MOTIF_TABLE[ 0b110001 + 1 ] = ALTL_MOTIF
      # Alt LastEx   LR SR
      MOTIF_TABLE[ 0b100001 + 1 ] = ALTL_MOTIF
      # Alt LastEx   SR SR
      MOTIF_TABLE[ 0b001001 + 1 ] = ALTL_MOTIF

      # RetainedInt  LL RR
      MOTIF_TABLE[ 0b101110 + 1 ] = RETI_MOTIF

      # Skipped      RR LL
      MOTIF_TABLE[ 0b110101 + 1 ] = SKIP_MOTIF
      #              RR LR
      MOTIF_TABLE[ 0b110100 + 1 ] = SKIP_MOTIF
      #              RR LS
      MOTIF_TABLE[ 0b110010 + 1 ] = SKIP_MOTIF

      #              LR LL
      MOTIF_TABLE[ 0b100101 + 1 ] = SKIP_MOTIF
      #              LR LR
      MOTIF_TABLE[ 0b100100 + 1 ] = SKIP_MOTIF
      #              LR LS
      MOTIF_TABLE[ 0b100010 + 1 ] = SKIP_MOTIF

      #              SR LL
      MOTIF_TABLE[ 0b001101 + 1 ] = SKIP_MOTIF
      #              SR LR
      MOTIF_TABLE[ 0b001100 + 1 ] = SKIP_MOTIF
      #              SR LS
      MOTIF_TABLE[ 0b001010 + 1 ] = SKIP_MOTIF

      # Alt Donor    LL LL
      MOTIF_TABLE[ 0b101101 + 1 ] = ALTD_MOTIF
      #              LL LR
      MOTIF_TABLE[ 0b101100 + 1 ] = ALTD_MOTIF
      #              LL LS
      MOTIF_TABLE[ 0b101010 + 1 ] = ALTD_MOTIF

      # Alt Acceptor RR RR
      MOTIF_TABLE[ 0b110110 + 1 ] = ALTA_MOTIF
      #              LR RR
      MOTIF_TABLE[ 0b100110 + 1 ] = ALTA_MOTIF
      #              SR RR
      MOTIF_TABLE[ 0b001110 + 1 ] = ALTA_MOTIF

Base.convert(::Type{EdgeMotif}, tup::Tuple{EdgeType,EdgeType}) = Base.convert(EdgeMotif, tup... )
Base.convert(::Type{EdgeMotif}, current::EdgeType, next::EdgeType) = MOTIF_TABLE[ (UInt8(current) << 3) | UInt8(next) + 1 ]

Base.convert(::Type{S}, edg::EdgeMotif ) where {S <: AbstractString} = MOTIF_STRING[ UInt8(edg) + 1 ]

#isobligate(  motif::EdgeMotif ) = motif != NONE_MOTIF && !( UInt8(motif) & 0b100 == 0b100 )
isobligate( motif::EdgeMotif ) = (0 <= UInt8(motif) <= 1)
isaltsplice( motif::EdgeMotif ) = (UInt8(motif) & 0b110) == 0b110

isspanning( edge::I, node::T ) where {I <: AbstractInterval, T <: Integer} = edge.first < node < edge.last ? true : false
isconnecting( edge::I, node::T ) where {I <: AbstractInterval, T <: Integer} = edge.first == node || edge.last == node ? true : false

# this is meant for short arrays when it is faster
# than using the overhead of a set
# DEPRECATED, IntSet is very efficient
function unique_push!( arr::Vector{T}, el::T ) where T
   if !( el in arr )
      push!( arr, el )
   end
end

# This holds one or many sets of connected
# nodes + the count of the reads and the eff_len
# We also store the min and max node of all sets
# Note: It did occur to me that a cleaner solution
# is PsiGraph holding paths::Vector{PsiPath} (deprecated)
# but the current implementation was choosen for
# faster memory access w.r.t. stride.
mutable struct PsiGraph
   psi::Vector{Float64}
   count::Vector{Float64}
   length::Vector{Float64}
   nodes::Vector{IntSet}
   min::NodeInt
   max::NodeInt
end

Base.sum( pgraph::PsiGraph ) = sum( pgraph.count )
Base.sum( ngraph::Nullable{PsiGraph} ) = isnull( ngraph ) ?
                                         0.0 : convert(Float64, sum( ngraph.value ))

function hasintersect( a::IntSet, b::IntSet )
   seta,setb = length(a) > length(b) ? (b,a) : (a,b)
   for elem in seta
      if elem in setb
         return true
      end
   end
   false
end

function hasintersect_terminal( a::IntSet, b::IntSet )
   if first(a) == last(b) || first(b) == last(a)
      return true
   else
      return false
   end
end

function hasintersect_terminal( a, b )
   if first(a) == last(b) || first(b) == last(a)
      return true
   else
      return false
   end
end

function Base.in( i::T, pgraph::PsiGraph ) where T <: Integer
   for nodeset in pgraph.nodes
      (i in nodeset) && return true
   end
   return false
end

complexity( one::PsiGraph, two::PsiGraph ) = complexity(length(one.nodes) + length(two.nodes))
complexity( one::PsiGraph ) = complexity(length(one.nodes))
# This function calculates complexity of a splicing event
# as the ceil log2 number of paths through the graph
complexity( num_paths::Int ) = @fastmath Int(ceil(log2( num_paths )))

# Calculate Shannon's Entropy, -Sum( Pi * log2( Pi ) )
function shannon_index( probs::Vector{Float64} )
   index = 0.0
   for i in 1:length(probs)
      prob = (probs[i] == 0.0 ? eps(Float64) : probs[i])
      index += prob * log2(prob)
   end
   index * -1
end

shannon_index( one::PsiGraph ) = shannon_index( one.psi )

# The Psi vectors in these two PsiGraphs need to sum to one!
# this does not assert this!
function shannon_index( one::PsiGraph, two::PsiGraph )
   index = 0.0
   for i in 1:length(one.psi)
      prob = (one.psi[i] == 0.0 ? eps(Float64) : one.psi[i])
      index += prob * log2(prob)
   end
   for j in 1:length(two.psi)
      prob = (two.psi[j] == 0.0 ? eps(Float64) : two.psi[j])
      index += prob * log2(prob)
   end
   index * -1
end

# Build minimal set of paths to explain graph
# Explanation: If we have 3 edges for example 1-2, 2-3, and 1-3
# We have to create 2 paths, 1-2-3, and 1-3, since 1-3 can't
# possibly contain the second node.
# This function will do this with any number of edges or complex
# paths through a graph and is the essence of how Whippet quantification
# works.. these graph paths are then assigned unambiguous counts, while
# ambiguous counts are then assigned to each graph path using the EM algorithm

function reduce_graph( edges::PsiGraph, graph::PsiGraph=deepcopy(edges), maxit::Int=100, maxiso::Int=1024 )
   newgraph = PsiGraph( Vector{IntSet}(),  Vector{Float64}(),
                        Vector{Float64}(), Vector{Float64}(),
                        graph.min, graph.max )
   it = 1
   while (newgraph.nodes != graph.nodes || it == 1) && it <= maxit
      if it > 1
         graph,newgraph = newgraph,graph
         empty!( newgraph.nodes  )
         empty!( newgraph.length )
         empty!( newgraph.count  )
         if length(graph.nodes) > maxiso
            ord = shuffle(1:length(graph.nodes))
            graph.nodes  = graph.nodes[ord[1:maxiso]]
            graph.length = graph.length[ord[1:maxiso]]
            graph.count  = graph.count[ord[1:maxiso]]
         end
      end
      @inbounds for i in 1:length(graph.nodes)
         push_i = false
         for j in 1:length(edges.nodes)
            if hasintersect_terminal( graph.nodes[i], edges.nodes[j] )
               push_i = true
               const jointset = union( graph.nodes[i], edges.nodes[j] )
               (jointset in newgraph.nodes) && continue
               push!( newgraph.nodes,  jointset )
               push!( newgraph.length, graph.length[i] + edges.length[j] )
               push!( newgraph.count,  graph.count[i] + edges.count[j] )
               newgraph.min = min( min(first(graph.nodes[i]), first(edges.nodes[j])), newgraph.min )
               newgraph.max = max( max(last(graph.nodes[i]),  last(edges.nodes[j])),  newgraph.max )
            end
         end
         if !push_i && !(graph.nodes[i] in newgraph.nodes)
            push!( newgraph.nodes, graph.nodes[i] )
            push!( newgraph.length, graph.length[i] )
            push!( newgraph.count, graph.count[i] )
            newgraph.min = min( first(graph.nodes[i]), newgraph.min )
            newgraph.max = max( last(graph.nodes[i]),  newgraph.max )
         end
      end
      it += 1
   end
   newgraph
end

@inline function reduce_graph_simple( edges::PsiGraph )
   if length(edges.nodes) == 2
      if hasintersect_terminal( edges.nodes[1], edges.nodes[2] )
         union!( edges.nodes[1], edges.nodes[2] )
         edges.length[1] += edges.length[2]
         edges.count[1]  += edges.count[2]
         pop!( edges.nodes )
         pop!( edges.length )
         pop!( edges.count )
      end
      return edges
   else
      return reduce_graph( edges )
   end
end

function Base.push!( pgraph::PsiGraph, edg::I;
                     value_bool=true, length=1.0 ) where I <: AbstractInterval
   push!( pgraph.count, value_bool ? get(edg.value) : 0.0 )
   push!( pgraph.length, value_bool ? length : 0.0 )
   push!( pgraph.nodes, IntSet([edg.first, edg.last]) )
   #reduce_graph_terminal!( pgraph )
   if edg.first < pgraph.min
      pgraph.min = edg.first
   end
   if edg.last > pgraph.max
      pgraph.max = edg.last
   end
end

mutable struct AmbigCounts
   paths::Vector{NodeInt}
   prop::Vector{Float64}
   prop_sum::Float64
   multiplier::Float64
end

Base.:(==)( a::AmbigCounts, b::AmbigCounts ) = a.paths == b.paths
Base.:(==)( a::AmbigCounts, b::IntSet ) = ( a.paths == b )
Base.:(==)( a::IntSet, b::AmbigCounts ) = ( a == b.paths )
Base.:(==)( a::Vector{I}, b::IntSet ) where {I <: Integer} = ( b == a )
function Base.:(==)( iset::IntSet, ivec::Vector{I} ) where I <: Integer
   length(iset) != length(ivec) && return false
   for intv in ivec
      !(intv in iset) && return false
   end
   true
end

Base.sum( nvec::Nullable{Vector{AmbigCounts}} ) = isnull( nvec ) ? 0.0 :
                                                  sum( nvec.value )
function Base.sum( vec::Vector{AmbigCounts} )
   sum = 0.0
   for i in 1:length(vec)
      sum += vec[i].multiplier
   end
   sum
end

function extend_edges!( edges::IntervalMap{K,V}, pgraph::PsiGraph, igraph::PsiGraph,
                        ambig_edge::Nullable{Vector{IntervalValue}}, node::NodeInt ) where {K,V}

   minv = min( pgraph.min, igraph.min )
   maxv = max( pgraph.max, igraph.max )
   idx = node - 1
   shouldpush = false
   while idx > minv # <----- node iter left
      #iterate through local edges to the left
      for edg in intersect( edges, (idx,idx) )
         # if this is a new edge and connects to our idx from the left
         if isconnecting( edg, idx ) && edg.last == idx && edg.first >= minv
            shouldpush = false
            if edg.last in pgraph
               push!( pgraph, edg, value_bool=false )
               shouldpush = true
            end
            if edg.last in igraph
               push!( igraph, edg, value_bool=false )
               shouldpush = true
            end
            if shouldpush
               if edg.first < minv
                  minv = edg.first
               end
               if isnull( ambig_edge )
                  ambig_edge = Nullable( Vector{IntervalValue}() )
               end
               push!( ambig_edge.value, edg )
            end
         end
      end
      idx -= 1
   end

   idx = node + 1
   while idx < maxv # iter right node ------->
      # iterate through local edges to the right
      for edg in intersect( edges, (idx,idx) )
         # if this is a new edge and connects to our idx from the left
         if isconnecting( edg, idx ) && edg.first == idx && edg.last <= maxv
            shouldpush = false
            if edg.first in pgraph
               push!( pgraph, edg, value_bool=false )
               shouldpush = true
            end
            if edg.first in igraph
               push!( igraph, edg, value_bool=false )
               shouldpush = true
            end
            if shouldpush
               if edg.last > maxv
                  maxv = edg.last
               end
               if isnull( ambig_edge )
                  ambig_edge = Nullable( Vector{IntervalValue}() )
               end
               push!( ambig_edge.value, edg )
            end
         end
      end
      idx += 1
   end

   ambig_edge
end

# add_node_counts! for single pgraph (used in tandem_utrs)
function add_node_counts!( ambig::Vector{AmbigCounts}, pgraph::PsiGraph,
                           sgquant::SpliceGraphQuant, bias::Float64 )

   iset = IntSet()
   for n in pgraph.min:pgraph.max # lets go through all possible nodes
      for i in 1:length(pgraph.nodes)
         if n in pgraph.nodes[i]
            push!( iset, i )
            pgraph.length[i] += sgquant.leng[n]
         end
      end

      #=                                            =#
      length( iset ) == 0 && continue # unspliced node
      if length( iset ) == 1 # non-ambiguous node
         if first( iset ) <= length(pgraph.nodes) # belongs to inclusion
            idx = first( iset )
            @fastmath pgraph.count[idx] += get(sgquant.node[n]) #* bias / sgquant.leng[n]
         end
      else
         #check if there is already an entry for this set of paths
         # if so, just increment multiplier, if not, make new one
         exists = false
         for am in ambig
            if iset == am
               @fastmath am.multiplier += get(sgquant.node[n]) #* bias / sgquant.leng[n]
               exists = true
               break
            end
         end
         if !exists && get(sgquant.node[n]) > 0
            push!( ambig, AmbigCounts( collect(iset),
                                       ones( length(iset) ) / length(iset),
                                       1.0, get(sgquant.node[n]) ) ) #* bias / sgquant.leng[n] ) )
         end
      end
      empty!( iset ) # clean up
   end
end

# add_node_counts! for two graphs, inc_graph and exc_graph (used in spliced events)
function add_node_counts!( ambig::Vector{AmbigCounts}, igraph::PsiGraph,
                           egraph::PsiGraph, sgquant::SpliceGraphQuant,
                           bias::Float64 )

   minv = min( egraph.min, igraph.min )
   maxv = max( egraph.max, igraph.max )
   iset = IntSet()
   for n in minv:maxv # lets go through all possible nodes
      for i in 1:length(igraph.nodes)
         if n in igraph.nodes[i]
            push!( iset, i )
            igraph.length[i] += 1 # sgquant.leng[n]
         end
      end
      for i in 1:length(egraph.nodes)
         if n in egraph.nodes[i]
            push!( iset, i+length(igraph.nodes) )
            egraph.length[i] += 1 # sgquant.leng[n]
         end
      end
      #=                                            =#
      length( iset ) == 0 && continue # unspliced node
      if length( iset ) == 1 # non-ambiguous node
         if first( iset ) <= length(igraph.nodes) # belongs to inclusion
            idx = first( iset )
            @fastmath igraph.count[idx] += get(sgquant.node[n]) * bias / sgquant.leng[n]
         else # belongs to exclusion
            idx = first( iset ) - length(igraph.nodes)
            @fastmath egraph.count[idx] += get(sgquant.node[n]) * bias / sgquant.leng[n]
         end
      else
         #check if there is already an entry for this set of paths
         # if so, just increment multiplier, if not, make new one
         exists = false
         for am in ambig
            if iset == am
               @fastmath am.multiplier += get(sgquant.node[n]) * bias / sgquant.leng[n]
               exists = true
               break
            end
         end
         if !exists && get(sgquant.node[n]) > 0
            push!( ambig, AmbigCounts( collect(iset),
                                       ones( length(iset) ) / length(iset),
                                       1.0, get(sgquant.node[n]) * bias / sgquant.leng[n] ) )
         end
      end
      empty!( iset ) # clean up
   end
end

# add_edge_counts! for one psigraph, (used in tandem utr)
function add_edge_counts!( ambig::Vector{AmbigCounts}, pgraph::PsiGraph,
                           sgquant::SpliceGraphQuant )
   iset = IntSet()
   for edg in sgquant.edge
      for i in 1:length(pgraph.nodes)
         if edg in pgraph.nodes[i]
            push!(iset, i)
            #pgraph.length[i] += 1
         end
      end
      #=                                             =#
      length( iset ) == 0 && continue
      if length( iset ) == 1
         idx = first( iset )
         pgraph.count[idx] += get(edg.value)
      else
         exists = false
         for am in ambig
            if iset == am
               am.multiplier += get(edg.value)
               exists = true
               break
            end
         end
         if !exists && get(edg.value) > 0
            push!( ambig, AmbigCounts( collect(iset),
                                       ones( length(iset) ) / length(iset),
                                       1.0, get(edg.value) ) )
         end
      end
      empty!( iset )
   end
end

# add_edge_counts! for two graphs, inc_graph and exc_graph! (used in spliced events)
function add_edge_counts!( ambig::Vector{AmbigCounts}, igraph::PsiGraph,
                           egraph::PsiGraph, edges::Vector{IntervalValue} )
   iset = IntSet()
   for edg in edges
      for i in 1:length(igraph.nodes)
         if edg in igraph.nodes[i]
            push!(iset, i)
            igraph.length[i] += 1
         end
      end
      for i in 1:length(egraph.nodes)
         if edg in egraph.nodes[i]
            push!( iset, i+length(igraph.nodes) )
            egraph.length[i] += 1
         end
      end
      #=                                             =#
      length( iset ) == 0 && continue
      if length( iset ) == 1
         if first( iset ) <= length(igraph.nodes)
            idx = first( iset )
            igraph.count[idx] += get(edg.value)
         else
            idx = first( iset ) - length(igraph.nodes)
            egraph.count[idx] += get(edg.value)
         end
      else
         exists = false
         for am in ambig
            if iset == am
               am.multiplier += get(edg.value)
               exists = true
               break
            end
         end
         if !exists && get(edg.value) > 0
            push!( ambig, AmbigCounts( collect(iset),
                                       ones( length(iset) ) / length(iset),
                                       1.0, get(edg.value) ) )
         end
      end
      empty!( iset )
   end
end

function build_utr_graph( nodes::IntSet, motif::EdgeMotif, sgquant::SpliceGraphQuant )

   utr_graph = Nullable(PsiGraph( Vector{Float64}(), Vector{Float64}(),
                                  Vector{Float64}(), Vector{IntSet}(),
                                  first(nodes), last(nodes) ))
   curset = IntSet()

   while length(nodes) > 0
      if motif == TXST_MOTIF
         curval = pop!(nodes)
      else
         curval = shift!(nodes)
      end
      push!( curset, curval )
      push!( get(utr_graph).nodes, copy(curset) )
      push!( get(utr_graph).psi, 0.0 )
      push!( get(utr_graph).length, 0.0 )
      push!( get(utr_graph).count, 0.0 )
   end

   utr_graph
end

function _process_tandem_utr( sg::SpliceGraph, sgquant::SpliceGraphQuant{C,R},
                              node::NodeInt, motif::EdgeMotif, readlen::Int ) where {C <: SGAlignContainer, R <: ReadCounter}

   utr_graph  = Nullable{PsiGraph}()
   ambig_cnt  = Nullable{Vector{AmbigCounts}}()

   used_node  = IntSet()
   total_cnt  = 0.0

   if motif == TXEN_MOTIF
      # add node directly upstream (CE next to tandemUTR)
      push!( used_node, node-1 )
      total_cnt += get(sgquant.node[node-1])
   end

   i = node
   curmotif = motif
   while curmotif == motif
      push!( used_node, i )
      total_cnt += get(sgquant.node[i])
      if i > 1
         interv = Interval{ExonInt}( i-1, i )
         total_cnt += get(get( sgquant.edge, interv, IntervalValue(0,0,default(R)) ).value)
      end
      i += 1
      if i < length(sg.edgetype)
         curmotif = convert(EdgeMotif, sg.edgetype[i], sg.edgetype[i+1] )
      else
         break
      end
   end

   if motif == TXST_MOTIF
      # add node directly downstream (CE next to tandemUTR)
      push!( used_node, i )
      total_cnt += get(sgquant.node[i])
      interv = Interval{ExonInt}( i-1, i )
      total_cnt += get(get( sgquant.edge, interv, IntervalValue(0,0,default(R)) ).value)
   end

   psi = Nullable{Float64}()
   len = length(used_node)
   if total_cnt > 2
      ambig_cnt = Nullable( Vector{AmbigCounts}() )
      utr_graph = build_utr_graph( used_node, motif, sgquant )
      add_node_counts!( ambig_cnt.value, utr_graph.value, sgquant, 1.0 )
      add_edge_counts!( ambig_cnt.value, utr_graph.value, sgquant )
      it = rec_tandem_em!( utr_graph.value, ambig_cnt.value, sig=4, readlen=readlen )
      psi = Nullable( get(utr_graph).psi )
   end

   psi,utr_graph,ambig_cnt,len
end

function _process_spliced( sg::SpliceGraph, sgquant::SpliceGraphQuant,
                           node::NodeInt, motif::EdgeMotif, bias::Float64, isnodeok::Bool,
                           minedgeweight::Float64=0.02 )

   inc_graph  = Nullable{PsiGraph}()
   exc_graph  = Nullable{PsiGraph}()
   ambig_edge = Nullable{Vector{IntervalValue}}()
   ambig_cnt  = Nullable{Vector{AmbigCounts}}()

   connecting_val = 0.0
   spanning_val   = 0.0

   maxvalue = length(sgquant.edge) > 0 ? get(maximum( sgquant.edge ).value) : 0.0

   for edg in intersect( sgquant.edge, (node, node) )

      (get(edg.value) / maxvalue < minedgeweight) && continue

      if   isconnecting( edg, node )
         if isnull( inc_graph )
            inc_graph = Nullable(PsiGraph( Vector{Float64}(), Vector{Float64}(),
                                           Vector{Float64}(), Vector{IntSet}(),
                                           edg.first, edg.last ))
         end
         if isnull( ambig_edge )
            ambig_edge = Nullable( Vector{IntervalValue}() )
         end
         push!( inc_graph.value, edg, value_bool=false )
         push!( ambig_edge.value, edg )
         connecting_val += get(edg.value)
      elseif isspanning( edg, node )
         if isnull( exc_graph ) #don't allocate unless there is alt splicing
            exc_graph = Nullable(PsiGraph( Vector{Float64}(), Vector{Float64}(),
                                           Vector{Float64}(), Vector{IntSet}(),
                                           edg.first, edg.last ))
         end
         if isnull( ambig_edge )
            ambig_edge = Nullable( Vector{IntervalValue}() )
         end
         push!( exc_graph.value, edg, value_bool=false )
         push!( ambig_edge.value, edg )
         spanning_val += get(edg.value)
      else
         error("Edge has to be connecting or spanning!!!!" )
      end
   end

   if !isnull( exc_graph ) && !isnull( inc_graph )
      ambig_cnt = Nullable( Vector{AmbigCounts}() )

      # try to bridge nodes by extending with potentially ambiguous edges
      ambig_edge = extend_edges!( sgquant.edge, exc_graph.value, inc_graph.value, ambig_edge, node )

      inc_graph = Nullable( reduce_graph( inc_graph.value ) )
      exc_graph = Nullable( reduce_graph( exc_graph.value ) )

      add_edge_counts!( ambig_cnt.value, inc_graph.value,
                        exc_graph.value, get(ambig_edge) )

      if isnodeok
         add_node_counts!( ambig_cnt.value, inc_graph.value, exc_graph.value, sgquant, bias )
      end
   else
      !isnull( inc_graph ) && (inc_graph = Nullable( reduce_graph_simple( inc_graph.value ) ))
      !isnull( exc_graph ) && (exc_graph = Nullable( reduce_graph_simple( exc_graph.value ) ))
   end # end expanding module

  # lets finish up now.
   psi = Nullable{Float64}()
   if isnull( exc_graph ) # no spanning edge
      # check if we have both inclusion edges represented, or one if alt 5'/3'
      if !isnull( inc_graph ) && ((connecting_val >= 2) ||
                                  (connecting_val >= 1  &&
                                   isaltsplice(motif)))
          psi = Nullable( 1.0 ) #&& likelihood_ci( psi, inc_cnt, z=1.64 )
      else
         # NA is ignored.
      end
   else # there is skipping
      if isnull( inc_graph )
         if spanning_val >= 1
             psi = Nullable( 0.0 ) #&& likelihood_ci( psi, exc_graph.count, z=1.64 )
         else
            # NA
         end
      else
         get(exc_graph).psi = zeros( length( get(exc_graph).count ) )
         get(inc_graph).psi = zeros( length( get(inc_graph).count ) )
         calculate_psi!( inc_graph.value, exc_graph.value, [get(inc_graph).count; get(exc_graph).count] )
         it = spliced_em!( inc_graph.value, exc_graph.value, ambig_cnt.value, sig=4 )
         psi = Nullable( sum( get(inc_graph).psi ) )
      end
   end
   total_reads = connecting_val + spanning_val
   psi,inc_graph,exc_graph,ambig_cnt,total_reads
end

function process_events( outfile, lib::GraphLib, graphq::GraphLibQuant;
                         isnodeok=false, iscircok=false, readlen::Int=50 )
   io = open( outfile, "w" )
   stream = ZlibDeflateOutputStream( io )
   output_psi_header( stream )
   for g in 1:length(lib.graphs)
      name = lib.names[g]
      chr  = lib.info[g].name
      strand = lib.info[g].strand ? '+' : '-'
      _process_events( stream, lib.graphs[g], graphq.quant[g], (name,chr,strand),
                       isnodeok=isnodeok, iscircok=iscircok, readlen=readlen )
   end
   close(stream)
   close(io)
end

function _process_events( io::BufOut, sg::SpliceGraph, sgquant::SpliceGraphQuant, info::GeneMeta;
                          isnodeok=false, iscircok=false, readlen::Int=50.0 )
   # Heres the plan:
   # step through sets of edges, look for edge motifs, some are obligate calculations
   # others we only calculate psi if there is an alternative edge
   if length( sg.edgetype ) <= 2
      return 0
   end
   i = 1
   while i < length(sg.edgetype)
      motif = convert(EdgeMotif, sg.edgetype[i], sg.edgetype[i+1] )
      next_istandem = i+1 < length(sg.edgetype) ? isobligate( convert(EdgeMotif, sg.edgetype[i+1], sg.edgetype[i+2]) ) : false
      if motif == NONE_MOTIF
         if !next_istandem
            output_empty( io, motif, sg, i, info )
         end
      elseif isobligate( motif ) # is utr event
         psi,utr,ambig,len = _process_tandem_utr( sg, sgquant, convert(NodeInt, i), motif, readlen )
         if !isnull( psi ) && !any( map( isnan, psi.value ) )
            total_cnt = sum(utr) + sum(ambig)
            i = output_utr( io, round.(get(psi),4), utr, total_cnt, motif, sg, i , info )
         else
            # psi/utr/total_cnt ignored here.
            i = output_utr( io, zeros(len), utr, 0.0, motif, sg, i, info, empty=true )
         end
      else  # is a spliced node
         bias = calculate_bias!( sgquant )
         psi,inc,exc,ambig,total_cnt = _process_spliced( sg, sgquant, convert(NodeInt, i), motif, bias, isnodeok )
         #total_cnt = sum(inc) + sum(exc) + sum(ambig)
         if !isnull( psi ) && 0 <= psi.value <= 1 && total_cnt > 0
            conf_int  = beta_posterior_ci( psi.value, total_cnt, sig=3 )
            output_psi( io, signif(psi.value,3), inc, exc, total_cnt, conf_int, motif, sg, sgquant.edge, i, info, bias  ) # TODO bias
         else
            output_empty( io, motif, sg, i, info )
         end
      end
      i += 1
   end
   # process back-splicing
   iscircok && output_circular( io, sg, sgquant, info )
end

function divsignif!( arr::Vector{N}, divisor::D, sig::I ) where { N <: Number, D <: Number, I <: Integer }
   if sig > 0
      for i in 1:length( arr )
         @fastmath arr[i] = signif( arr[i] / divisor, sig )
      end
   else
      for i in 1:length( arr )
         @fastmath arr[i] /= divisor
      end
   end
end

# Beta(
@inline function beta_posterior_ci( p, n; ci=0.9, sig=0 )
   lo_q = (1 - ci)/2
   hi_q = 1 - lo_q
   beta = Beta( p*n + 1, (1-p)*n + 1 )
   if sig > 0
      lo = signif( quantile(beta, lo_q), sig )
      hi = signif( quantile(beta, hi_q), sig )
   else
      lo = quantile(beta, lo_q)
      hi = quantile(beta, hi_q)
   end
   lo,hi
end

# spliced psi
function calculate_psi!( igraph::PsiGraph, egraph::PsiGraph, counts::Vector{Float64}; sig=0 )
   for i in 1:length(igraph.psi)
      igraph.psi[i] = counts[i] / igraph.length[i]
   end
   for i in 1:length(egraph.psi)
      idx = i + length(igraph.psi)
      egraph.psi[i] = counts[idx] / egraph.length[i]
   end
   cnt_sum = sum( egraph.psi ) + sum( igraph.psi )
   divsignif!( igraph.psi, cnt_sum, sig )
   divsignif!( egraph.psi, cnt_sum, sig )
end

function calculate_psi!( pgraph::PsiGraph, counts::Vector{Float64}; sig=0, readlen::Int=50)
   for i in 1:length(pgraph.psi)
      pgraph.psi[i] = counts[i] / max( pgraph.length[i] - readlen, 1 )
   end
   cnt_sum = sum( pgraph.psi )
   divsignif!( pgraph.psi, cnt_sum, sig )
end

function spliced_em!( igraph::PsiGraph, egraph::PsiGraph, ambig::Vector{AmbigCounts};
                      it=1, maxit=1500, sig=0 )

   inc_temp   = zeros(length(igraph.count))
   exc_temp   = zeros(length(egraph.count))
   count_temp = ones(length(igraph.count)+length(egraph.count))

   while (inc_temp != igraph.psi || egraph.psi != exc_temp) && it < maxit

      unsafe_copy!( count_temp, igraph.count, indx_shift=0 )
      unsafe_copy!( count_temp, egraph.count, indx_shift=length(igraph.count) )
      unsafe_copy!( inc_temp, igraph.psi )
      unsafe_copy!( exc_temp, egraph.psi )

      for ac in ambig
         idx = 1
         if it > 1
            ac.prop_sum = 0.0
            for p in ac.paths
               prev_psi = p <= length(igraph.psi) ? igraph.psi[p] : egraph.psi[p-length(igraph.psi)]
               ac.prop[idx] = prev_psi
               ac.prop_sum += prev_psi
               idx += 1
            end
         end

         idx = 1
         for p in ac.paths
            @fastmath prop = ac.prop[idx] / ac.prop_sum
            ac.prop[idx] = prop
            @fastmath count_temp[p] += prop * ac.multiplier
            idx += 1
         end
      end

      calculate_psi!( igraph, egraph, count_temp, sig=sig )

      it += 1
   end
   it
end

function rec_tandem_em!( pgraph::PsiGraph, ambig::Vector{AmbigCounts};
                          utr_temp::Vector{Float64}=zeros(length(pgraph.count)),
                          count_temp::Vector{Float64}=ones(length(pgraph.count)),
                          it=1, maxit=1500, sig=0, readlen::Int=50 )

   unsafe_copy!( count_temp, pgraph.count, indx_shift=0 )
   unsafe_copy!( utr_temp, pgraph.psi )

   for ac in ambig
      idx = 1
      if it > 1
         ac.prop_sum = 0.0
         for p in ac.paths
            prev_psi = pgraph.psi[p]
            ac.prop[idx] = prev_psi
            ac.prop_sum += prev_psi
            idx += 1
         end
      end

      idx = 1
      for p in ac.paths
         @fastmath prop = ac.prop[idx] / ac.prop_sum
         ac.prop[idx] = prop
         @fastmath count_temp[p] += prop * ac.multiplier
         idx += 1
      end
   end

   calculate_psi!( pgraph, count_temp, sig=sig, readlen=readlen )

   if utr_temp != pgraph.psi && it < maxit
      it = rec_tandem_em!( pgraph, ambig,
                           utr_temp=utr_temp, count_temp=count_temp,
                           it=it+1, maxit=maxit, sig=sig, readlen=readlen )
   end
   it
end


text:
# This file includes an example file for 1D interpolation.

using FractalTools
using Makie

# Generate data
xi, dx, xf = 0., 0.1, 10.
x = collect(xi : dx : xf)
y = sin.(x)
pts = collect.(zip(x, y))
interp = interpolate(pts, Interp1D(0.1))

# Calculate interpolant
xt = collect(xi : 0.1dx : xf)
yt = interp.(xt)

# Plot interpolation
fig, ax, plt = lines(xt, yt, color=yt, linewidth=5)
scatter!(ax, x, y, markersize=5, color=:black)

ext: jl


text:
using ExecutableSpecifications.Gherkin

"""
The location of a step definition (the Julia code of a test) is a filename and a line number.
"""
struct StepDefinitionLocation
    filename::String
    lineno::Int
end

"Thrown if there is no step definition matching a scenario step."
struct NoMatchingStepDefinition <: Exception end

"Thrown when more than one step definition matches a scenario step."
struct NonUniqueStepDefinition <: Exception
    locations::Vector{StepDefinitionLocation}
end

"A step definition matcher takes a scenario step and finds a matching step definition to execute."
abstract type StepDefinitionMatcher end

# A StepDefinitionMatcher should define a method
# findstepdefinition(::StepDefinitionMatcher, ::Gherkin.ScenarioStep)

"A step definition has a description, which is used to find it, a function to execute, and a location."
struct StepDefinition
    description::String
    definition::Function
    location::StepDefinitionLocation
end

"""
The context in which a step definition executes. This context is used to share data between
different step definitions. It is created newly for each scenario. Thus, two scenarios cannot share
data.
"""
struct StepDefinitionContext
    variables::Dict{Symbol, Any}

    StepDefinitionContext() = new(Dict{Symbol, Any}())
end

"Find a variable value given a symbol name."
Base.getindex(context::StepDefinitionContext, sym::Symbol) = context.variables[sym]

"Set a variable value given a symbol name and a value."
Base.setindex!(context::StepDefinitionContext, value::Any, sym::Symbol) = context.variables[sym] = value

"Check for a mapping for a given key"
Base.haskey(context::StepDefinitionContext, sym::Symbol) = haskey(context.variables, sym)

#
# Global state
#

currentdefinitions = Vector{StepDefinition}()
currentfilename = ""

#
# Step definition macros
#

function step_definition_(description::String, definition::Expr)
    # The step definition function takes a context and executes a bit of code supplied by the
    # test writer. The bit of code is in $definition.
    definitionfunction = :(context -> $definition)
    quote
        # Push a given step definition to the global state so it can be found by the
        # `StepDefinitionMatcher`.
        push!(currentdefinitions, StepDefinition($description, (context) -> begin
                try
                    # Escape the step definition code so it gets the proper scope.
                    $(esc(definitionfunction))(context)
                    # Any step definition that does not throw an exception is successful.
                    SuccessfulStepExecution()
                catch ex
                    # StepAssertFailures are turned into a failed result here, but all other exceptions
                    # are propagated.
                    if ex isa StepAssertFailure
                        StepFailed(ex.assertion)
                    else
                        rethrow()
                    end
                end
            end,
            StepDefinitionLocation(currentfilename, 0)))
    end
end

"Provide a more user friendly @given macro for a step definition."
macro given(description, definition)
    step_definition_(description, definition)
end

"Provide a more user friendly @when macro for a step definition."
macro when(description, definition)
    step_definition_(description, definition)
end

"Provide a more user friendly @then macro for a step definition."
macro then(description, definition)
    step_definition_(description, definition)
end

#
# Step Definition Matcher implementation
#

"""
Finds step definitions defined in a Julia file with the @given, @when, @then macros defined above.
Takes a source text as input and reads the code defined in it.
"""
struct FromMacroStepDefinitionMatcher <: StepDefinitionMatcher
    stepdefinitions::Vector{StepDefinition}
    filename::String

    function FromMacroStepDefinitionMatcher(source::AbstractString; filename::String="<no filename>")
        global currentdefinitions
        global currentfilename
        currentfilename = filename
        # Read all step definitions as Julia code.
        include_string(Main, source)
        # Save the step definitions found in the global variable `currentdefinitions` into a local
        # variable, so that we can clear the global state and read another file.
        mydefinitions = currentdefinitions
        this = new(mydefinitions, filename)
        currentdefinitions = Vector{StepDefinition}()
        this
    end
end

"""
    findstepdefinition(matcher::FromMacroStepDefinitionMatcher, step::Gherkin.ScenarioStep)

Find a step definition that has a description that matches the provided scenario step.
If no such step definition is found, throw a `NoMatchingStepDefinition`.
If more than one such step definition is found, throw a `NonUniqueStepDefinition`.
"""
function findstepdefinition(matcher::FromMacroStepDefinitionMatcher, step::Gherkin.ScenarioStep)
    matchingindexes = findall(x -> x.description == step.text, matcher.stepdefinitions)
    matchingstepdefinitions = [matcher.stepdefinitions[i] for i in matchingindexes]
    if isempty(matchingstepdefinitions)
        throw(NoMatchingStepDefinition())
    end
    if length(matchingstepdefinitions) > 1
        locations = map(stepdefinition -> StepDefinitionLocation(matcher.filename, 0),
                        matchingstepdefinitions)
        throw(NonUniqueStepDefinition(locations))
    end
    matchingstepdefinitions[1]
end

#
# Composite matcher
#
"""
Find step definitions from multiple other step definition matchers.
"""
mutable struct CompositeStepDefinitionMatcher <: StepDefinitionMatcher
    matchers::Vector{StepDefinitionMatcher}

    CompositeStepDefinitionMatcher(matchers...) = new([matchers...])
end

function findstepdefinition(composite::CompositeStepDefinitionMatcher, step::Gherkin.ScenarioStep)
    matching = StepDefinition[]
    nonuniquesfound = StepDefinitionLocation[]
    # Recursively call `findstepdefinition(...)` on all sub-matchers.
    # When they throw a `NonUniqueStepDefinition`, record the location so it can be shown to the
    # user where the step definitions are.
    # Ignore `NonUniqueStepDefinition` exceptions, as normally all but one of the matchers will
    # throw it.
    for m in composite.matchers
        try
            stepdefinition = findstepdefinition(m, step)
            push!(matching, stepdefinition)
        catch ex
            if ex isa NonUniqueStepDefinition
                append!(nonuniquesfound, ex.locations)
            end
        end
    end
    if length(matching) > 1 || !isempty(nonuniquesfound)
        locations = vcat(nonuniquesfound, [d.location for d in matching])
        throw(NonUniqueStepDefinition(locations))
    end
    if isempty(matching)
        throw(NoMatchingStepDefinition())
    end
    matching[1]
end

function addmatcher!(composite::CompositeStepDefinitionMatcher, matcher::StepDefinitionMatcher)
    push!(composite.matchers, matcher)
end


text:
module Wavelets

include("util.jl")

include("wavelettypes.jl")
include("transforms.jl")

include("threshold.jl")
include("plot.jl")

using Reexport
@reexport using .Util, .WaveletTypes, .Transforms, .Threshold, .Plot

end


text:
"""

	X_1 ∩ X_2 ∩ ... ∩ X_n = ∅

"""

struct EmptyNIntersection <: Constraint

	X::Vector{Variable}

end

function EmptyNIntersection(m::Model, X::Vector{Variable})

	#push!(m.varsInter, Variable(1:m.p, "F[$(length(m.varsInter)+1)]"))
	Variable(model, 1:m.p, "F[$(length(m.varsInter)+1)]")
	Intersection(m, m.varsInter[end], X[1], X[2])

	for ind in 3:(length(X)-1)
		# push!(m.varsInter, Variable(1:m.p, "F[$(length(m.varsInter)+1)]"))
		Variable(model, 1:m.p, "F[$(length(m.varsInter)+1)]")
		Intersection(m, m.varsInter[end], m.varsInter[end-1], X[ind])
	end

	EmptyIntersection(m, m.varsInter[end], X[end])

	m.varsInter[end-(length(X)-2):end]

end


text:
# TODO:
#   1) better convergence testing

function learn!(
        model::Model,
        opt::Optimizer;
        iterations=1000,
        ftol::Float64=1e-6
    )
    x = paramvec(model) # returns view of all parameters as a vector
    ∇ = gradvec(model)  # returns view of gradient as a vector
    f = value(model)
    fhist = [f]

    converged = false
    iter = 0
    while !converged && iter < iterations
        iter += 1
        f = update!(x,∇,f,model,opt)    # uodates x
        grad!(model)                    # updates ∇
        converged = (fhist[end]-f)/fhist[end] < ftol
        push!(fhist,f)
    end

    return fhist,converged
end

max_stars_repo_name: JuliaML/StochasticOptimization

max_issues_repo_name: JuliaML/StochasticOptimization

max_forks_repo_name: JuliaML/StochasticOptimization


text:
function sampleFromDiscrete(probs,Pw,nonzeroLength)
# draw a sample from unnormalized discrete distribution, probs
# Pw is the total of discrete distributions
    temp = rand()*Pw;
    total = 0.0;
    for i = 1:nonzeroLength
        total = total + probs[i];
        if temp<total
            return i;
            break
        end
    end
end

function generateAlias(p, k)
# creates a 2xnumTopics alias table
    # p is the unnormalized probability distribution of size 1xk
    # k is the number of outcomes
    # returns a generated alias table of size kx2

    # initialization...
    AT = zeros(k,2); # generated alias table of size kx2
    # temporary array to store low and high probabilities
    L_probInd = Int64[]; L_prob = Float64[];H_probInd = Int64[]; H_prob = Float64[];
    Lind = 1; Hind = 1; # current index of AT and L or H
    kInv = 1/k; # average probability of the distribution
    sum_p = sum(p);

    # generation alias table...
    for i = 1:k
        p[i] = p[i] ./ sum_p; # normalizing the previously unnormalized discrete distribution
        if p[i]<=kInv
            push!(L_probInd,i);
            push!(L_prob,p[i]);
            Lind+= 1;
        else
            push!(H_probInd,i);
            push!(H_prob,p[i]);
            Hind+= 1;
        end
        AT[i,1]=i; AT[i,2]=i*kInv; # initial formation of alias table
    end
    Lind-= 1; Hind-= 1; # this points to the last added element
    while Lind>=1 && Hind>=1
        ATind = L_probInd[Lind]; # extract last index from L_probInd
        AT[ATind,1] = H_probInd[Hind]; # extract new index from H_probInd and map it to AT
        AT[ATind,2] = (ATind-1)*kInv+L_prob[Lind]; # place low probability value in AT
        Ph = H_prob[Hind]-(kInv-L_prob[Lind]); # cut alias probability from high value probabilities
        Lind-= 1; Hind-= 1; # points to next entry we can extract
        if Ph>kInv
            Hind+= 1;
            H_probInd[Hind] = AT[ATind,1]; H_prob[Hind] = Ph;
        else
            Lind+= 1;
            L_probInd[Lind] = AT[ATind,1];L_prob[Lind] = Ph;
        end
    end
    return AT;
end

function sampleAlias(AT,k)
# generates a sample from alias table
    bin = rand(1:k);
    if AT[bin,2]>rand()
        s = bin;
    else
        s = convert(Int64,AT[bin,1]); # to make sure we are getting integer samples from alias table
    end
    return s;
end

function aliasMethodOnWord(word,Alpha,wordTopicCounts,topicCounts,Beta,BetaSum,numTopics,numCachedSamples,cachedSamples,cached_qw)
    qw = Array{Float64}(numTopics); # unnormalized discrete distribution of dense bucket
    for topic = 1:numTopics
        qw[topic] = Alpha[topic] .* (wordTopicCounts[word,topic] + Beta[word]) ./ (topicCounts[topic] + BetaSum);
        cached_qw[word,topic] = qw[topic];
    end
    aliasTable = generateAlias(qw,numTopics);
    for sample = 1:numCachedSamples
        cachedSamples[word,sample] = sampleAlias(aliasTable,numTopics);
    end

    return cachedSamples,cached_qw;
end


text:
# Julia wrapper for header: ImfCRgbaFile.h
# Automatically generated using Clang.jl

function ImfFloatToHalf(f, h)
    ccall((:ImfFloatToHalf, libOpenEXR), Cvoid, (Cfloat, Ptr{ImfHalf}), f, h)
end

function ImfFloatToHalfArray(n, f, h)
    ccall((:ImfFloatToHalfArray, libOpenEXR), Cvoid, (Cint, Ptr{Cfloat}, Ptr{ImfHalf}), n, f, h)
end

function ImfHalfToFloat(h)
    ccall((:ImfHalfToFloat, libOpenEXR), Cfloat, (ImfHalf,), h)
end

function ImfHalfToFloatArray(n, h, f)
    ccall((:ImfHalfToFloatArray, libOpenEXR), Cvoid, (Cint, Ptr{ImfHalf}, Ptr{Cfloat}), n, h, f)
end

function ImfNewHeader()
    ccall((:ImfNewHeader, libOpenEXR), Ptr{ImfHeader}, ())
end

function ImfDeleteHeader(hdr)
    ccall((:ImfDeleteHeader, libOpenEXR), Cvoid, (Ptr{ImfHeader},), hdr)
end

function ImfCopyHeader(hdr)
    ccall((:ImfCopyHeader, libOpenEXR), Ptr{ImfHeader}, (Ptr{ImfHeader},), hdr)
end

function ImfHeaderSetDisplayWindow(hdr, xMin, yMin, xMax, yMax)
    ccall((:ImfHeaderSetDisplayWindow, libOpenEXR), Cvoid, (Ptr{ImfHeader}, Cint, Cint, Cint, Cint), hdr, xMin, yMin, xMax, yMax)
end

function ImfHeaderDisplayWindow(hdr, xMin, yMin, xMax, yMax)
    ccall((:ImfHeaderDisplayWindow, libOpenEXR), Cvoid, (Ptr{ImfHeader}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), hdr, xMin, yMin, xMax, yMax)
end

function ImfHeaderSetDataWindow(hdr, xMin, yMin, xMax, yMax)
    ccall((:ImfHeaderSetDataWindow, libOpenEXR), Cvoid, (Ptr{ImfHeader}, Cint, Cint, Cint, Cint), hdr, xMin, yMin, xMax, yMax)
end

function ImfHeaderDataWindow(hdr, xMin, yMin, xMax, yMax)
    ccall((:ImfHeaderDataWindow, libOpenEXR), Cvoid, (Ptr{ImfHeader}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), hdr, xMin, yMin, xMax, yMax)
end

function ImfHeaderSetPixelAspectRatio(hdr, pixelAspectRatio)
    ccall((:ImfHeaderSetPixelAspectRatio, libOpenEXR), Cvoid, (Ptr{ImfHeader}, Cfloat), hdr, pixelAspectRatio)
end

function ImfHeaderPixelAspectRatio(hdr)
    ccall((:ImfHeaderPixelAspectRatio, libOpenEXR), Cfloat, (Ptr{ImfHeader},), hdr)
end

function ImfHeaderSetScreenWindowCenter(hdr, x, y)
    ccall((:ImfHeaderSetScreenWindowCenter, libOpenEXR), Cvoid, (Ptr{ImfHeader}, Cfloat, Cfloat), hdr, x, y)
end

function ImfHeaderScreenWindowCenter(hdr, x, y)
    ccall((:ImfHeaderScreenWindowCenter, libOpenEXR), Cvoid, (Ptr{ImfHeader}, Ptr{Cfloat}, Ptr{Cfloat}), hdr, x, y)
end

function ImfHeaderSetScreenWindowWidth(hdr, width)
    ccall((:ImfHeaderSetScreenWindowWidth, libOpenEXR), Cvoid, (Ptr{ImfHeader}, Cfloat), hdr, width)
end

function ImfHeaderScreenWindowWidth(hdr)
    ccall((:ImfHeaderScreenWindowWidth, libOpenEXR), Cfloat, (Ptr{ImfHeader},), hdr)
end

function ImfHeaderSetLineOrder(hdr, lineOrder)
    ccall((:ImfHeaderSetLineOrder, libOpenEXR), Cvoid, (Ptr{ImfHeader}, Cint), hdr, lineOrder)
end

function ImfHeaderLineOrder(hdr)
    ccall((:ImfHeaderLineOrder, libOpenEXR), Cint, (Ptr{ImfHeader},), hdr)
end

function ImfHeaderSetCompression(hdr, compression)
    ccall((:ImfHeaderSetCompression, libOpenEXR), Cvoid, (Ptr{ImfHeader}, Cint), hdr, compression)
end

function ImfHeaderCompression(hdr)
    ccall((:ImfHeaderCompression, libOpenEXR), Cint, (Ptr{ImfHeader},), hdr)
end

function ImfHeaderSetIntAttribute(hdr, name, value)
    ccall((:ImfHeaderSetIntAttribute, libOpenEXR), Cint, (Ptr{ImfHeader}, Ptr{UInt8}, Cint), hdr, name, value)
end

function ImfHeaderIntAttribute(hdr, name, value)
    ccall((:ImfHeaderIntAttribute, libOpenEXR), Cint, (Ptr{ImfHeader}, Ptr{UInt8}, Ptr{Cint}), hdr, name, value)
end

function ImfHeaderSetFloatAttribute(hdr, name, value)
    ccall((:ImfHeaderSetFloatAttribute, libOpenEXR), Cint, (Ptr{ImfHeader}, Ptr{UInt8}, Cfloat), hdr, name, value)
end

function ImfHeaderSetDoubleAttribute(hdr, name, value)
    ccall((:ImfHeaderSetDoubleAttribute, libOpenEXR), Cint, (Ptr{ImfHeader}, Ptr{UInt8}, Cdouble), hdr, name, value)
end

function ImfHeaderFloatAttribute(hdr, name, value)
    ccall((:ImfHeaderFloatAttribute, libOpenEXR), Cint, (Ptr{ImfHeader}, Ptr{UInt8}, Ptr{Cfloat}), hdr, name, value)
end

function ImfHeaderDoubleAttribute(hdr, name, value)
    ccall((:ImfHeaderDoubleAttribute, libOpenEXR), Cint, (Ptr{ImfHeader}, Ptr{UInt8}, Ptr{Cdouble}), hdr, name, value)
end

function ImfHeaderSetStringAttribute(hdr, name, value)
    ccall((:ImfHeaderSetStringAttribute, libOpenEXR), Cint, (Ptr{ImfHeader}, Ptr{UInt8}, Ptr{UInt8}), hdr, name, value)
end

function ImfHeaderStringAttribute(hdr, name, value)
    ccall((:ImfHeaderStringAttribute, libOpenEXR), Cint, (Ptr{ImfHeader}, Ptr{UInt8}, Ptr{Cstring}), hdr, name, value)
end

function ImfHeaderSetBox2iAttribute(hdr, name, xMin, yMin, xMax, yMax)
    ccall((:ImfHeaderSetBox2iAttribute, libOpenEXR), Cint, (Ptr{ImfHeader}, Ptr{UInt8}, Cint, Cint, Cint, Cint), hdr, name, xMin, yMin, xMax, yMax)
end

function ImfHeaderBox2iAttribute(hdr, name, xMin, yMin, xMax, yMax)
    ccall((:ImfHeaderBox2iAttribute, libOpenEXR), Cint, (Ptr{ImfHeader}, Ptr{UInt8}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), hdr, name, xMin, yMin, xMax, yMax)
end

function ImfHeaderSetBox2fAttribute(hdr, name, xMin, yMin, xMax, yMax)
    ccall((:ImfHeaderSetBox2fAttribute, libOpenEXR), Cint, (Ptr{ImfHeader}, Ptr{UInt8}, Cfloat, Cfloat, Cfloat, Cfloat), hdr, name, xMin, yMin, xMax, yMax)
end

function ImfHeaderBox2fAttribute(hdr, name, xMin, yMin, xMax, yMax)
    ccall((:ImfHeaderBox2fAttribute, libOpenEXR), Cint, (Ptr{ImfHeader}, Ptr{UInt8}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}), hdr, name, xMin, yMin, xMax, yMax)
end

function ImfHeaderSetV2iAttribute(hdr, name, x, y)
    ccall((:ImfHeaderSetV2iAttribute, libOpenEXR), Cint, (Ptr{ImfHeader}, Ptr{UInt8}, Cint, Cint), hdr, name, x, y)
end

function ImfHeaderV2iAttribute(hdr, name, x, y)
    ccall((:ImfHeaderV2iAttribute, libOpenEXR), Cint, (Ptr{ImfHeader}, Ptr{UInt8}, Ptr{Cint}, Ptr{Cint}), hdr, name, x, y)
end

function ImfHeaderSetV2fAttribute(hdr, name, x, y)
    ccall((:ImfHeaderSetV2fAttribute, libOpenEXR), Cint, (Ptr{ImfHeader}, Ptr{UInt8}, Cfloat, Cfloat), hdr, name, x, y)
end

function ImfHeaderV2fAttribute(hdr, name, x, y)
    ccall((:ImfHeaderV2fAttribute, libOpenEXR), Cint, (Ptr{ImfHeader}, Ptr{UInt8}, Ptr{Cfloat}, Ptr{Cfloat}), hdr, name, x, y)
end

function ImfHeaderSetV3iAttribute(hdr, name, x, y, z)
    ccall((:ImfHeaderSetV3iAttribute, libOpenEXR), Cint, (Ptr{ImfHeader}, Ptr{UInt8}, Cint, Cint, Cint), hdr, name, x, y, z)
end

function ImfHeaderV3iAttribute(hdr, name, x, y, z)
    ccall((:ImfHeaderV3iAttribute, libOpenEXR), Cint, (Ptr{ImfHeader}, Ptr{UInt8}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), hdr, name, x, y, z)
end

function ImfHeaderSetV3fAttribute(hdr, name, x, y, z)
    ccall((:ImfHeaderSetV3fAttribute, libOpenEXR), Cint, (Ptr{ImfHeader}, Ptr{UInt8}, Cfloat, Cfloat, Cfloat), hdr, name, x, y, z)
end

function ImfHeaderV3fAttribute(hdr, name, x, y, z)
    ccall((:ImfHeaderV3fAttribute, libOpenEXR), Cint, (Ptr{ImfHeader}, Ptr{UInt8}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}), hdr, name, x, y, z)
end

function ImfHeaderSetM33fAttribute(hdr, name, m)
    ccall((:ImfHeaderSetM33fAttribute, libOpenEXR), Cint, (Ptr{ImfHeader}, Ptr{UInt8}, Ptr{NTuple{3, Cfloat}}), hdr, name, m)
end

function ImfHeaderM33fAttribute(hdr, name, m)
    ccall((:ImfHeaderM33fAttribute, libOpenEXR), Cint, (Ptr{ImfHeader}, Ptr{UInt8}, Ptr{NTuple{3, Cfloat}}), hdr, name, m)
end

function ImfHeaderSetM44fAttribute(hdr, name, m)
    ccall((:ImfHeaderSetM44fAttribute, libOpenEXR), Cint, (Ptr{ImfHeader}, Ptr{UInt8}, Ptr{NTuple{4, Cfloat}}), hdr, name, m)
end

function ImfHeaderM44fAttribute(hdr, name, m)
    ccall((:ImfHeaderM44fAttribute, libOpenEXR), Cint, (Ptr{ImfHeader}, Ptr{UInt8}, Ptr{NTuple{4, Cfloat}}), hdr, name, m)
end

function ImfOpenOutputFile(name, hdr, channels)
    ccall((:ImfOpenOutputFile, libOpenEXR), Ptr{ImfOutputFile}, (Ptr{UInt8}, Ptr{ImfHeader}, Cint), name, hdr, channels)
end

function ImfCloseOutputFile(out)
    ccall((:ImfCloseOutputFile, libOpenEXR), Cint, (Ptr{ImfOutputFile},), out)
end

function ImfOutputSetFrameBuffer(out, base, xStride, yStride)
    ccall((:ImfOutputSetFrameBuffer, libOpenEXR), Cint, (Ptr{ImfOutputFile}, Ptr{ImfRgba}, Csize_t, Csize_t), out, base, xStride, yStride)
end

function ImfOutputWritePixels(out, numScanLines)
    ccall((:ImfOutputWritePixels, libOpenEXR), Cint, (Ptr{ImfOutputFile}, Cint), out, numScanLines)
end

function ImfOutputCurrentScanLine(out)
    ccall((:ImfOutputCurrentScanLine, libOpenEXR), Cint, (Ptr{ImfOutputFile},), out)
end

function ImfOutputHeader(out)
    ccall((:ImfOutputHeader, libOpenEXR), Ptr{ImfHeader}, (Ptr{ImfOutputFile},), out)
end

function ImfOutputChannels(out)
    ccall((:ImfOutputChannels, libOpenEXR), Cint, (Ptr{ImfOutputFile},), out)
end

function ImfOpenTiledOutputFile(name, hdr, channels, xSize, ySize, mode, rmode)
    ccall((:ImfOpenTiledOutputFile, libOpenEXR), Ptr{ImfTiledOutputFile}, (Ptr{UInt8}, Ptr{ImfHeader}, Cint, Cint, Cint, Cint, Cint), name, hdr, channels, xSize, ySize, mode, rmode)
end

function ImfCloseTiledOutputFile(out)
    ccall((:ImfCloseTiledOutputFile, libOpenEXR), Cint, (Ptr{ImfTiledOutputFile},), out)
end

function ImfTiledOutputSetFrameBuffer(out, base, xStride, yStride)
    ccall((:ImfTiledOutputSetFrameBuffer, libOpenEXR), Cint, (Ptr{ImfTiledOutputFile}, Ptr{ImfRgba}, Csize_t, Csize_t), out, base, xStride, yStride)
end

function ImfTiledOutputWriteTile(out, dx, dy, lx, ly)
    ccall((:ImfTiledOutputWriteTile, libOpenEXR), Cint, (Ptr{ImfTiledOutputFile}, Cint, Cint, Cint, Cint), out, dx, dy, lx, ly)
end

function ImfTiledOutputWriteTiles(out, dxMin, dxMax, dyMin, dyMax, lx, ly)
    ccall((:ImfTiledOutputWriteTiles, libOpenEXR), Cint, (Ptr{ImfTiledOutputFile}, Cint, Cint, Cint, Cint, Cint, Cint), out, dxMin, dxMax, dyMin, dyMax, lx, ly)
end

function ImfTiledOutputHeader(out)
    ccall((:ImfTiledOutputHeader, libOpenEXR), Ptr{ImfHeader}, (Ptr{ImfTiledOutputFile},), out)
end

function ImfTiledOutputChannels(out)
    ccall((:ImfTiledOutputChannels, libOpenEXR), Cint, (Ptr{ImfTiledOutputFile},), out)
end

function ImfTiledOutputTileXSize(out)
    ccall((:ImfTiledOutputTileXSize, libOpenEXR), Cint, (Ptr{ImfTiledOutputFile},), out)
end

function ImfTiledOutputTileYSize(out)
    ccall((:ImfTiledOutputTileYSize, libOpenEXR), Cint, (Ptr{ImfTiledOutputFile},), out)
end

function ImfTiledOutputLevelMode(out)
    ccall((:ImfTiledOutputLevelMode, libOpenEXR), Cint, (Ptr{ImfTiledOutputFile},), out)
end

function ImfTiledOutputLevelRoundingMode(out)
    ccall((:ImfTiledOutputLevelRoundingMode, libOpenEXR), Cint, (Ptr{ImfTiledOutputFile},), out)
end

function ImfOpenInputFile(name)
    ccall((:ImfOpenInputFile, libOpenEXR), Ptr{ImfInputFile}, (Ptr{UInt8},), name)
end

function ImfCloseInputFile(in)
    ccall((:ImfCloseInputFile, libOpenEXR), Cint, (Ptr{ImfInputFile},), in)
end

function ImfInputSetFrameBuffer(in, base, xStride, yStride)
    ccall((:ImfInputSetFrameBuffer, libOpenEXR), Cint, (Ptr{ImfInputFile}, Ptr{ImfRgba}, Csize_t, Csize_t), in, base, xStride, yStride)
end

function ImfInputReadPixels(in, scanLine1, scanLine2)
    ccall((:ImfInputReadPixels, libOpenEXR), Cint, (Ptr{ImfInputFile}, Cint, Cint), in, scanLine1, scanLine2)
end

function ImfInputHeader(in)
    ccall((:ImfInputHeader, libOpenEXR), Ptr{ImfHeader}, (Ptr{ImfInputFile},), in)
end

function ImfInputChannels(in)
    ccall((:ImfInputChannels, libOpenEXR), Cint, (Ptr{ImfInputFile},), in)
end

function ImfInputFileName(in)
    ccall((:ImfInputFileName, libOpenEXR), Cstring, (Ptr{ImfInputFile},), in)
end

function ImfOpenTiledInputFile(name)
    ccall((:ImfOpenTiledInputFile, libOpenEXR), Ptr{ImfTiledInputFile}, (Ptr{UInt8},), name)
end

function ImfCloseTiledInputFile(in)
    ccall((:ImfCloseTiledInputFile, libOpenEXR), Cint, (Ptr{ImfTiledInputFile},), in)
end

function ImfTiledInputSetFrameBuffer(in, base, xStride, yStride)
    ccall((:ImfTiledInputSetFrameBuffer, libOpenEXR), Cint, (Ptr{ImfTiledInputFile}, Ptr{ImfRgba}, Csize_t, Csize_t), in, base, xStride, yStride)
end

function ImfTiledInputReadTile(in, dx, dy, lx, ly)
    ccall((:ImfTiledInputReadTile, libOpenEXR), Cint, (Ptr{ImfTiledInputFile}, Cint, Cint, Cint, Cint), in, dx, dy, lx, ly)
end

function ImfTiledInputReadTiles(in, dxMin, dxMax, dyMin, dyMax, lx, ly)
    ccall((:ImfTiledInputReadTiles, libOpenEXR), Cint, (Ptr{ImfTiledInputFile}, Cint, Cint, Cint, Cint, Cint, Cint), in, dxMin, dxMax, dyMin, dyMax, lx, ly)
end

function ImfTiledInputHeader(in)
    ccall((:ImfTiledInputHeader, libOpenEXR), Ptr{ImfHeader}, (Ptr{ImfTiledInputFile},), in)
end

function ImfTiledInputChannels(in)
    ccall((:ImfTiledInputChannels, libOpenEXR), Cint, (Ptr{ImfTiledInputFile},), in)
end

function ImfTiledInputFileName(in)
    ccall((:ImfTiledInputFileName, libOpenEXR), Cstring, (Ptr{ImfTiledInputFile},), in)
end

function ImfTiledInputTileXSize(in)
    ccall((:ImfTiledInputTileXSize, libOpenEXR), Cint, (Ptr{ImfTiledInputFile},), in)
end

function ImfTiledInputTileYSize(in)
    ccall((:ImfTiledInputTileYSize, libOpenEXR), Cint, (Ptr{ImfTiledInputFile},), in)
end

function ImfTiledInputLevelMode(in)
    ccall((:ImfTiledInputLevelMode, libOpenEXR), Cint, (Ptr{ImfTiledInputFile},), in)
end

function ImfTiledInputLevelRoundingMode(in)
    ccall((:ImfTiledInputLevelRoundingMode, libOpenEXR), Cint, (Ptr{ImfTiledInputFile},), in)
end

function ImfNewRound12logLut(channels)
    ccall((:ImfNewRound12logLut, libOpenEXR), Ptr{ImfLut}, (Cint,), channels)
end

function ImfNewRoundNBitLut(n, channels)
    ccall((:ImfNewRoundNBitLut, libOpenEXR), Ptr{ImfLut}, (UInt32, Cint), n, channels)
end

function ImfDeleteLut(lut)
    ccall((:ImfDeleteLut, libOpenEXR), Cvoid, (Ptr{ImfLut},), lut)
end

function ImfApplyLut(lut, data, nData, stride)
    ccall((:ImfApplyLut, libOpenEXR), Cvoid, (Ptr{ImfLut}, Ptr{ImfRgba}, Cint, Cint), lut, data, nData, stride)
end

function ImfErrorMessage()
    ccall((:ImfErrorMessage, libOpenEXR), Cstring, ())
end


text:
module test_julia_splat

using Test

@test Base.splat(+)((1,2,3)) == .+(1,2,3) == 6

end # module test_julia_splat


text:
export savename, @dict, @ntuple, @strdict
export ntuple2dict, dict2ntuple

"""
    allaccess(c)
Return all the keys `c` can be accessed using [`access`](@ref).
For dictionaries/named tuples this is `keys(c)`,
for everything else it is `fieldnames(typeof(c))`.
"""
allaccess(c::AbstractDict) = collect(keys(c))
allaccess(c::NamedTuple) = keys(c)
allaccess(c::Any) = fieldnames(typeof(c))
allaccess(c::DataType) = fieldnames(c)
allaccess(c::String) = error("`c` must be a container, not a string!")

"""
    access(c, key)
Access `c` with given key. For `AbstractDict` this is `getindex`,
for anything else it is `getproperty`.

    access(c, keys...)
When given multiple keys, `access` is called recursively, i.e.
`access(c, key1, key2) = access(access(c, key1), key2)` and so on.
For example, if `c, c.k1` are `NamedTuple`s then
`access(c, k1, k2) == c.k1.k2`.

!!! note
    Please only extend the single key method when customizing `access`
    for your own Types.
"""
access(c, keys...) = access(access(c, keys[1]), Base.tail(keys)...)
access(c::AbstractDict, key) = getindex(c, key)
access(c, key) = getproperty(c, key)

"""
    default_allowed(c) = (Real, String, Symbol)
Return the (super-)Types that will be used as `allowedtypes`
in [`savename`](@ref) or other similar functions.
"""
default_allowed(c) = (Real, String, Symbol)

"""
    default_prefix(c) = ""
Return the `prefix` that will be used by default
in [`savename`](@ref) or other similar functions.
"""
default_prefix(c) = ""

"""
    savename([prefix,], c [, suffix]; kwargs...)
Create a shorthand name, commonly used for saving a file, based on the
parameters in the container `c` (`Dict`, `NamedTuple` or any other Julia
composite type, e.g. created with Parameters.jl). If provided use
the `prefix` and end the name with `.suffix` (i.e. you don't have to include
the `.` in your `suffix`).

The function chains keys and values into a string of the form:
```julia
key1=val1_key2=val2_key3=val3
```
while the keys are **always sorted alphabetically.** If you provide
the prefix/suffix the function will do:
```julia
prefix_key1=val1_key2=val2_key3=val3.suffix
```
assuming you chose the default `connector`, see below. Notice
that `prefix` can be any path and in addition if
it ends as a path (`/` or `\\`) then the `connector` is ommited.

`savename` can be very conveniently combined with
[`@dict`](@ref) or [`@ntuple`](@ref).

## Keywords
* `allowedtypes = default_allowed(c)` : Only values of type subtyping
  anything in `allowedtypes` are used in the name. By default
  this is `(Real, String, Symbol)`.
* `accesses = allaccess(c)` : You can also specify which specific keys you want
  to use with the keyword `accesses`. By default this is all possible
  keys `c` can be accessed with, see [`allaccess`](@ref).
* `digits = 3` : Floating point values are rounded to `digits`.
  In addition if the following holds:
  ```julia
  round(val; digits = digits) == round(Int, val)
  ```
  then the integer value is used in the name instead.
* `connector = "_"` : string used to connect the various entries.

## Examples
```julia
d = (a = 0.153456453, b = 5.0, mode = "double")
savename(d; digits = 4) == "a=0.1535_b=5_mode=double"
savename("n", d) == "n_a=0.153_b=5_mode=double"
savename("n/", d) == "n/a=0.153_b=5_mode=double"
savename(d, "n") == "a=0.153_b=5_mode=double.n"
savename("data/n", d, "n") == "data/n_a=0.153_b=5_mode=double.n"
savename("n", d, "n"; connector = "-") == "n-a=0.153-b=5-mode=double.n"
savename(d, allowedtypes = (String,)) == "mode=double"

rick = (never = "gonna", give = "you", up = "!");
savename(rick) == "give=you_never=gonna_up=!" # keys are sorted!
```
"""
savename(c; kwargs...) = savename(default_prefix(c), c, ""; kwargs...)
savename(c::Any, suffix::String; kwargs...) =
    savename(default_prefix(c), c, suffix; kwargs...)
savename(prefix::String, c::Any; kwargs...) = savename(prefix, c, ""; kwargs...)
function savename(prefix::String, c, suffix::String;
                  allowedtypes = default_allowed(c),
                  accesses = allaccess(c), digits = 3,
                  connector = "_")

    labels = vecstring(accesses) # make it vector of strings
    p = sortperm(labels)
    first = prefix == "" || prefix[end] == '\\' || prefix[end] == '/'
    s = prefix
    for j ∈ p
        val = access(c, accesses[j])
        label = labels[j]
        t = typeof(val)
        if any(x -> (t <: x), allowedtypes)
            !first && (s *= connector)
            if t <: AbstractFloat
                if round(val; digits = digits) == round(Int, val)
                    val = round(Int, val)
                else
                    val = round(val; digits = digits)
                end
            end
            s *= label*"="*string(val);
            first = false
        end
    end
    suffix != "" && (s *= "."*suffix)
    return s
end

"""
    @dict vars...
Create a dictionary out of the given variables that has as keys the variable
names and as values their values.

## Examples
```jldoctest; setup = :(using DrWatson)
julia> ω = 5; χ = "test"; ζ = π/3;

julia> @dict ω χ ζ
Dict{Symbol,Any} with 3 entries:
  :ω => 5
  :χ => "test"
  :ζ => 1.0472
```
"""
macro dict(vars...)
    expr = Expr(:call, :Dict)
    for i in 1:length(vars)
        push!(expr.args, :($(QuoteNode(vars[i])) => $(esc(vars[i]))))
    end
    return expr
end

"""
    @strdict vars...
Same as [`@dict`](@ref) but the key type is `String`.
"""
macro strdict(vars...)
    expr = Expr(:call, :Dict)
    for i in 1:length(vars)
        push!(expr.args, :(string($(QuoteNode(vars[i]))) => $(esc(vars[i]))))
    end
    return expr
end

"""
    @ntuple vars...
Create a `NamedTuple` out of the given variables that has as keys the variable
names and as values their values.

## Examples
```jldoctest; setup = :(using DrWatson)
julia> ω = 5; χ = "test"; ζ = 3.14;

julia> @ntuple ω χ ζ
(ω = 5, χ = "test", ζ = 3.14)
```
"""
macro ntuple(vars...)
   args = Any[]
   for i in 1:length(vars)
       push!(args, Expr(:(=), esc(vars[i]), :($(esc(vars[i])))))
   end
   expr = Expr(:tuple, args...)
   return expr
end
# Credit of `ntuple` macro goes to Sebastian Pfitzner, @pfitzseb

"""
    ntuple2dict(nt) -> dict
Convert a `NamedTuple` to a dictionary.
"""
ntuple2dict(nt::NamedTuple) = Dict(k => nt[k] for k in keys(nt))

"""
    dict2ntuple(dict) -> ntuple
Convert a dictionary (with `Symbol` or `String` as key type) to
a `NamedTuple`.
"""
function dict2ntuple(dict::Dict{String, T}) where T
    NamedTuple{Tuple(Symbol.(keys(dict)))}(values(dict))
end
function dict2ntuple(dict::Dict{Symbol, T}) where T
    NamedTuple{Tuple(keys(dict))}(values(dict))
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule PulseAudio_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("PulseAudio")
JLLWrappers.@generate_main_file("PulseAudio", UUID("02771fc1-bdb7-5db5-8d11-300768e00fbd"))
end  # module PulseAudio_jll


text:
# finiteDMRG.jl
function finiteDMRG(mps::FiniteMPS, mpo::MPO, final_position::Int,
                    DMRG_parameters::Dict, truncation_parameters::Dict)
    L = length(mps.Ts)

    if mps.canonical_form[1] != :right_canonical
        make_right_canonical(mps)
    end

    max_sweep = get(DMRG_parameters,"max_sweep",50)
    min_sweep = get(DMRG_parameters,"min_sweep",5)
    max_E_error = get(DMRG_parameters,"max_E_error",10^-8)
    max_S_error = get(DMRG_parameters,"max_S_error",10^-3)

    env = initialize_FiniteEnvironment(mps,mpo,final_position)

    println("Sweep 1")
    E = one_sweep(mps,mpo,env,truncation_parameters)
    chi = get_chi(mps)
    println("chi = ", chi)
    println("E = ", E)
    E_list = [E]
    S = entanglement_entropy(mps)
    println("S = ", S)
    S_list = [S]

    for i in 2:max_sweep
        println("Sweep ", i)
        E = one_sweep(mps,mpo,env,truncation_parameters)
        chi = get_chi(mps)
        println("chi = ", chi)
        println("E = ", E)
        push!(E_list,E)
        delta_E = E_list[end]-E_list[end-1]
        println("delta_E = ", delta_E)
        delta_E > max_E_error*0.1 && error("the energy increases, stop DMRG!")
        -delta_E < max_E_error && println("The energy has converged!")
        S = entanglement_entropy(mps)
        push!(S_list,S)
        delta_S = mapreduce(x->abs(x), +, S_list[end]-S_list[end-1])/L
        println("delta_S = ", delta_S)
        delta_S < max_S_error && println("The entanglement entropy has converged!")
        if (-delta_E < max_E_error) && (delta_S < max_S_error) && (i >= min_sweep)
            println("Both energy and entropy has converged! DMRG finishes!")
            break
        end
        if i == max_sweep
            println("NOTE!!! The maximal number of sweep has been reached, but the
                simulation has not been converged to the given precision!!!")
        end
    end
end

function one_sweep(mps::FiniteMPS,mpo::MPO,env::FiniteEnvironment,
                    truncation_parameters::Dict)
    E_list_one_sweep = [ ]
    truncation_error_list_one_sweep = [ ]

    # sweep from left to right
    for i in 2:length(mps.Ts)
        U, Σ, Vʰ, ϵ, E, info = one_update(i, mps, mpo, env, truncation_parameters)
        push!(E_list_one_sweep, E)
        push!(truncation_error_list_one_sweep, ϵ)
        mps.Ts[i-1] = U
        mps.Ss[i] = Σ
        mps.Ts[i] = permute(Σ*Vʰ, (1,3),(2,))
        update_left_env(i,mps,mpo,env)
    end

    # sweep from right to left
    for i in length(mps.Ts):-1:2
        U, Σ, Vʰ, ϵ, E, info = one_update(i, mps, mpo, env, truncation_parameters)
        push!(E_list_one_sweep, E)
        push!(truncation_error_list_one_sweep, ϵ)
        mps.Ts[i-1] = U*Σ
        mps.Ss[i] = Σ
        mps.Ts[i] = permute(Vʰ,(1,3),(2,))
        update_right_env(i,mps,mpo,env)
    end
    println("E_list_one_sweep = ", E_list_one_sweep)
    println("truncation_error_list_one_sweep = ", truncation_error_list_one_sweep)
    return E_list_one_sweep[end]
end

function one_update(i::Int, mps::FiniteMPS, mpo::MPO, env::FiniteEnvironment,
                        truncation_parameters::Dict)
    chi_max = get(truncation_parameters,"chi_max",50)
    svd_min = get(truncation_parameters,"svd_min",10^-10)
    truncation_error = get(truncation_parameters,"truncation_error", 10^-8)

    i > 1 || error("bond should be larger than 1")
    i < length(mps.Ts) + 1 || error("bond should be smaller than the length of mps")

    theta_ini = TensorMap(undef, codomain(mps.Ts[i-1])*codomain(mps.Ts[i],2),
                            domain(mps.Ts[i]))
    @tensor theta_ini[a b d;e] = mps.Ts[i-1][a b;c]*mps.Ts[i][c d;e]

    E, theta_new, info = eigsolve(theta_ini, 1, :SR, Lanczos()) do x
        build_two_site_H_eff(x, mpo.Ws[i-1], mpo.Ws[i],
                                env.left_env[i-1], env.right_env[i+1])
    end
    U, Σ, Vʰ, ϵ = tsvd(theta_new[1], (1,2), (4,3),
               trunc = truncdim(chi_max) & truncbelow(svd_min) & truncerr(truncation_error),
               p = 2, alg= SDD())
    return U, Σ, Vʰ, ϵ, E[1], info
end

function build_two_site_H_eff(theta::AbstractTensorMap{S,3,1},
                              W1::AbstractTensorMap{S,2,2}, W2::AbstractTensorMap{S,2,2},
                              left_env::AbstractTensorMap{S,1,2},
                              right_env::AbstractTensorMap{S,2,1}) where S
    T_tmp1 = TensorMap(undef, codomain(left_env)*codomain(theta,2)*codomain(theta,3),
                        domain(left_env,2)*domain(theta))
    @tensor T_tmp1[a d e;c f] = left_env[a;b c]*theta[b d e;f]

    T_tmp2 = TensorMap(undef,codomain(left_env)*codomain(theta,3)*codomain(W1,2),
                        domain(theta)*domain(W1,2))
    @tensor T_tmp2[e f b;g d] = T_tmp1[e c f;a g]*W1[a b;c d]

    T_tmp3 = TensorMap(undef,codomain(left_env)*codomain(W1,2)*codomain(W2,2),
                        domain(theta)*domain(W2,2))
    @tensor T_tmp3[e f b;g d] = T_tmp2[e c f;g a]*W2[a b;c d]

    theta_H = TensorMap(undef, codomain(left_env)*codomain(W1,2)*codomain(W2,2),
                            domain(right_env))
    @tensor theta_H[d e f;c] = T_tmp3[d e f;a b]*right_env[a b;c]
    return theta_H
end


text:
""" ISSNeuralODE

    Dynamics models of the form
        dxdt = f(x,u)
           y = h(x) + g(u)

    * dynamics: the function f
    * links: the function h
    * feedthrough: the function g
    * modifiers: a parametrized function that takes parameters of f and of h,
      and returns modified parameters of f and h. For example, implements
      stability constraints.

    Catalog:

        dynamics
            affine
                AffineOde
                LinearOde
            ctrnn0
                CTRNN0Elman
                CTRNN0Canonical
                CTRNN0Zero
            ctrnn1
                CTRNN1Nominal
                CTRNN1Canonical
                CTRNN1Zero
        links
            affine_canonical
                AffineCanonical
                LinearCanonical
        feedthrough
            linear_feedthrough
                LinearFeedthrough
                ResistiveFeedthrough
        modifiers
            Nominal
            ISSCTRNN0
            ISSCTRNN1

            regularization
                no_regularization
                lds_condition_V
                lds_condition_ctrnn0
                lds_condition_ctrnn1
"""

using DiffEqFlux, Flux
import DiffEqFlux: paramlength, initial_params
using DifferentialEquations,
    DiffEqFlux, DiffEqSensitivity, SteadyStateDiffEq, Flux, Zygote, Random
using IterTools: ncycle
using ParameterSchedulers
using ParameterSchedulers: Stateful, next!
using Printf

include("dynamics/dynamics.jl")
include("links/links.jl")
include("modifiers/modifiers.jl")
include("modifiers/regularization.jl")
include("feedthrough/feedthrough.jl")
include("interpolator.jl")
include("trainer.jl")
include("forward_solve.jl")
include("initial_condition.jl")
include("configurations.jl")

max_stars_repo_name: syanga/ISSNeuralODE

max_issues_repo_name: syanga/ISSNeuralODE

max_forks_repo_name: syanga/ISSNeuralODE


text:
using Genie

import Genie.Router: route, @params
import Base.convert

convert(::Type{Int}, s::SubString{String}) = parse(Int, s)

route("/sum/:x::Int/:y::Int") do
    @params(:x) + @params(:y)
end

route("/hello") do
    "Welcome to Genie!"
end

route("/") do
    "root"
end

Genie.AppServer.startup(8080, "0.0.0.0", async=true)

max_stars_repo_name: kallooa-bhr/dockerfiles

max_issues_repo_name: kallooa-bhr/dockerfiles

max_forks_repo_name: kallooa-bhr/dockerfiles


text:
function sum_vector(x::AbstractArray)
s = 0.0
for i = 1:length(x)
s = s + x[i]
end
return s
end
function sum_cols_matrix_views(x::Array{Float64, 2})
num_cols = size(x, 2); num_rows = size(x, 1)
s = zeros(num_cols)
for i = 1:num_cols
s[i] = sum_vector(sub(x, 1:num_rows, i))
end
return s
end

max_stars_repo_path: Module 2/Chapter06/array_views_ad.jl

max_issues_repo_path: Module 2/Chapter06/array_views_ad.jl

max_forks_repo_path: Module 2/Chapter06/array_views_ad.jl


text:
####
#### expmv utils
####
diagshift!(A, mu) = A - mu * I # fallback out-of-place (for sparse) or lazy (for LinearMaps) operator
diagshift!(A::StridedArray, mu) = A .-= mu * I # modify in-place

####
#### norm1est utils
####
function A_pow_n_B!(res, A, n::Integer, v)
    tmp = similar(v)
    mul!(res, A, v)
    for i in 1:n-1
        mul!(tmp, A, res)
        copyto!(res, tmp)
    end
end

function At_pow_n_B!(res, A, n::Integer, v)
    tmp = similar(v)
    mul!(res, adjoint(A), v)
    for i in 1:n-1
        mul!(tmp, adjoint(A), res)
        copyto!(res, tmp)
    end
end

####
#### normest1 utils
####
function lazyfindmax(f,x)
    idx = 1 # default
    maxval = f(first(x)) # default
    for (i,xi) in enumerate(x)
        (i == 1) && continue
        fval = f(xi)
        if fval > maxval
            idx, maxval = i, fval
        end
    end
    return maxval, idx
end
lazyfindmax(x) = lazyfindmax(identity,x)

# extension of Base.sign where x==0
zerosign(x) = x == zero(x) ? one(x) : Base.sign(x)


text:
function uhlig05()
#
#  test polynomial used by F. Uhlig
#
    p = poly([-1*ones(6); 2; 2]);
    z = [-1.0 6; 2 2];

    p, PolyZeros(z)
end


text:
using GeneralizedGenerated
using TupleVectors: chainvec
import MeasureTheory: testvalue

export testvalue
EmptyNTtype = NamedTuple{(),Tuple{}} where T<:Tuple

@inline function testvalue(mc::ModelClosure; cfg = NamedTuple(), ctx=NamedTuple())
    gg_call(mc, testvalue, cfg, ctx, KeepReturn())
end

###############################################################################
# ctx::NamedTuple

@inline function tilde(::typeof(testvalue), lens::typeof(identity), xname, x, d, cfg, ctx::NamedTuple)
    xnew = testvalue(d)
    ctx′ = merge(ctx, NamedTuple{(dynamic(xname),)}((xnew,)))
    (xnew, ctx′, ctx′)
end

@inline function tilde(::typeof(testvalue), lens, xname, x, d, cfg, ctx::NamedTuple)
    xnew = set(x, Lens!!(lens), testvalue(d))
    ctx′ = merge(ctx, NamedTuple{(dynamic(xname),)}((xnew,)))
    (xnew, ctx′, ctx′)
end


text:
using tempfile: NamedTemporaryFile
using textwrap: wrap
abstract type AbstractFileOp end
mutable struct FileOp <: AbstractFileOp
    file_name
    method
    file
end
function __enter__(self::FileOp)
    self.file = open(self.file_name, self.method)
    return self.file
end

function __exit__(self::FileOp, type_ = nothing, value = nothing, traceback = nothing)
    close(self.file)
end

if abspath(PROGRAM_FILE) == @__FILE__
    NamedTempFile::new() do temp_file
        file_path = temp_file.name
        FileOp(file_path, "w") do file
            write(file, "test")
        end
        open(file_path, "r") do f
            @assert(read(f) == "test")
        end
    end
end


text:
using CSV
using MortalityTables

# load the table
table = MortalityTables.table(3277)

select_filepath = joinpath(@__DIR__, "3277_sel.csv")
select_csv = CSV.File(select_filepath)
# test the select values
for row in select_csv
    age, dur, val = row
    if !(table.select[age][age + dur - 1] == val)
        error("The values don't match")
    end
end

ultimate_filepath = joinpath(@__DIR__, "3277_ult.csv")
ultimate_csv = CSV.File(ultimate_filepath)
for row in ultimate_csv
    age, val = row
    if !(table.ultimate[age] == val)
        error("The values don't match")
    end
end

max_stars_repo_name: actuarialopensource/validations

max_issues_repo_name: actuarialopensource/validations

max_forks_repo_name: actuarialopensource/validations


text:
@testset "Vector-vector projection" begin

    x_unit = [1, 0]
    y_unit = [0, 1]

    project(x_unit, x_unit) == x_unit
    project(y_unit, y_unit) == y_unit

    # The projection is the zero vector if u and v are perpendicular.
    project(x_unit, y_unit) == [0, 0]
    project(y_unit, x_unit) == [0, 0]
    project([5, 3], [-3, 5]) == [0, 0]

    @test project([1, 0], x_unit) == [1, 0]
    @test project([1, 1], x_unit) == [1, 0]
    @test project([5, 5], x_unit) == [5, 0]
    @test project([5, -5], x_unit) == [5, 0]
    @test project([-5, 5], x_unit) == [-5, 0]

    # The magnitude of the second vector shouldn't matter.
    @test project([1, 0], 2 * x_unit) == [1, 0]

    @test project([1, 0], [1, 1]) == [0.5, 0.5]
    @test project([1, 2, 4], [5, 1, -1]) == [5, 1, -1] / 9
end

@testset "Point-line projection" begin

    line_x = Line([0, 0], [1, 0])
    @test project([0, 0], line_x) == [0, 0]
    @test project([10, 0], line_x) == [10, 0]
    @test project([0, 10], line_x) == [0, 0]
    @test project([10, 10], line_x) == [10, 0]

    line_diag = Line([0, 0], [1, 1])
    @test project([0, 0], line_diag) == [0, 0]
    @test project([1, 0], line_diag) == 0.5 * ones(2)

    line_diag_3d = Line([0, 0, 0], [1, 1, 1])
    @test project([1, 0, 0], line_diag_3d) ≈ 1/3 * ones(3)

    # The magnitude of the direction vector should not matter.
    @test project([1, 0, 0], Line([0, 0, 0], -2 * ones(Int, 3))) ≈ 1/3 * ones(3)

    @test project([50, 10], Line([1, -5], [0, 3])) == [1, 10]
end

@testset "Point-plane projection" begin

    plane_xy = Plane([0, 0, 0], [0, 0, 1])

    @test project([0, 0, 0], plane_xy) == [0, 0, 0]
    @test project([1, 0, 0], plane_xy) == [1, 0, 0]
    @test project([50, -34, 0], plane_xy) == [50, -34, 0]
    @test project([50, -34, 50], plane_xy) == [50, -34, 0]
    @test project([50, -34, -150], plane_xy) == [50, -34, 0]

    # The magnitude of the plane normal should not matter.
    @test project([50, -34, -150], Plane([0, 0, 0], [0, 0, 5])) == [50, -34, 0]
    @test project([50, -34, -150], Plane([0, 0, 0], [0, 0, -5])) == [50, -34, 0]

    @test project([50, -34, -150], Plane([0, 0, -5], [0, 0, 1])) == [50, -34, -5]
    @test project([50, -34, -150], Plane([18, 23, -5], [0, 0, 1])) == [50, -34, -5]

    @test project([50, -34, -150], Plane([18, 23, -5], [0, 0, 1])) == [50, -34, -5]

    @test project([7, -3, 5], Plane([1, 2, 3], [1, 1, 1])) == [6, -4, 4]
end


text:
################################################################################
#
# WRI dual formulation (frequency domain version)
#
################################################################################

module WRIdual

using LinearAlgebra, SparseArrays
using ImageFiltering
using PyPlot

# const R = Float32
# const C = ComplexF32
const R = Float64
const C = ComplexF64
const RuC = Union{R, C}
export R, C, RuC

include("utils.jl")
include("RestrictionLinAlg.jl")
include("HelmholtzLinAlg.jl")
include("Objectives_FWI.jl")
include("Objectives_WRIdual_eps.jl")
include("Objectives_WRIdual_lambda.jl")
include("PDEWeightFun.jl")

end


text:
using  ForwardDiff,GenericLinearAlgebra,LinearAlgebra

mutable struct svdmat
    U::Any
    S::Any
    V::Any
end
export svd_lowrank

function svd_lowrank(A::Matrix, l::Int, q::Int)
    #= note:: The implementation is based on the Algorithm 5.1 from
    Halko et al, 2009.=#
    Q = get_approximate_basis(A, l, q)
    Q = conj(Q')
    B = Q * A
    U, S, V = svd(B)
    V = conj(V')
    U = Q * U
    return svdmat(U, S, V)
end

function get_approximate_basis(A::Matrix, l::Int, q::Int)
    #= note:: The implementation is based on the Algorithm 4.4 from
    Halko et al, 2009. =#
    m, n = size(A)
    omega = randn(n, l)
    Y = A * omega
    Q, R = qr(Y)
    Y_H = conj(Y')

    for i = 1:q
        Q, R = qr(Y_H * Q)
        Q, R = qr(Y * Q)
    end
    return Q
end

max_stars_repo_name: SamanthaWangdl/ReSparse

max_issues_repo_name: SamanthaWangdl/ReSparse

max_forks_repo_name: SamanthaWangdl/ReSparse


text:
module SBMLCounterfactual

using Omega
using StatsBase
using Random
using Plots
using Distributions
using PyCall
using DataStructures

te = pyimport("tellurium")

function get_parameters(sbml_file_path)
    """
    Function to extract model parameters from sbml. Returns dictionary with
    a transition matrix, species list, and rates list.
    """
    ## Load SBML
    sbml = te.loadSBMLModel(sbml_file_path)

    ## Transition Matrix
    transition = sbml.getFullStoichiometryMatrix()
    transition = transpose(transition)

    ## Extract Species
    species_ids = sbml.getFloatingSpeciesIds()
    species_starting_vals = sbml.getFloatingSpeciesAmounts()
    species = Dict(species_ids . => species_starting_vals)

    ## Extract Rates
    rates = Dict(sbml.getGlobalParameterIds() . => sbml.getGlobalParameterValues())

    return Dict("transition" => transition, "species" => species, "rates" => rates)

end

end # module


text:
function mutaterandom(individual, mutationRate)
    if rand() < mutationRate
        s = size(individual, 1)
        mutatepoints = rand(1:s,2)
        individual[mutatepoints[1]], individual[mutatepoints[2]] = individual[mutatepoints[2]], individual[mutatepoints[1]]
    end
    return individual
end

max_stars_repo_name: ariririri/heuristicopt

max_issues_repo_name: ariririri/heuristicopt


text:
# This will take a text file of space-separated RGB values and create the
# proper `register_gradient_colors(:magma, sample_evenly([...]))` calls

function build_grad_string(fn::AbstractString, libname, separator = ' ')
    name = splitext(basename(fn))[1]
    if length(name) > 4 && name[end-3:end] == "-rgb"
        name = name[1:end-4]
    end
    @show name
    colors = [begin
        "    RGB($(join(split(strip(l), separator), ", "))),\n"
    end for l in eachline(open(fn))]
    "\nregister_gradient_colors(:$(name), sample_evenly([\n$(join(colors))], $(min(64, length(colors)))), :$(libname))\n"
end

function build_grad_file(dir::AbstractString, outname::AbstractString, separator = ' ')
    outfn = joinpath(dirname(@__FILE__()), outname*".jl")
    @show outfn
    outfile = open(outfn, "w")
    write(outfile, "register_color_library(:$(outname), ColorLibrary(Dict(:default => :sequential)))\n")
    for (root, dirs, files) in walkdir(dir)
        @show root, dirs, files
        for fn in files
            #if length(fn) > 8 && fn[end-7:end] == "-rgb.txt" #specific to cmocean
            if fn[end-3:end] == ".csv"
                write(outfile, build_grad_string(joinpath(root, fn), outname, separator))
            end
        end
    end
    close(outfile)
end


text:
# function to sort hoppings
"""
    sort_hoppings!(sys::TBHamiltonian)

Sorts the hoppings of the TBHamiltonian (sys.hoppings), such that they are ordered first by the start orbital, then by the destination orbital and finally by the direction indices. This is necessary to work with sparce Hamiltonians.
"""
function sort_hoppings!(sys::TBHamiltonian)

    sort!(sys.hoppings, by = hop -> (hop.from, hop.to, hop.dir...))
    sys.sorted[1] = true

end

# utility functions to extract information from TBHamiltonian
"""
    dimension(sys::TBHamiltonian{D, T})  where {D, T}

Returns the spacial dimensional of a tight-binding Hamiltonian.
"""
function dimension(sys::TBHamiltonian{D, T}) where {D, T}
    return D
end

"""
    hopping_type(sys::TBHamiltonian{D, T}) where {D, T}

Returns the type of hopping of a tight-binding Hamiltonian.
"""
function hopping_type(sys::TBHamiltonian{D, T})  where {D, T}
    return T
end

"""
    norbs(sys::TBHamiltonian)

Returns the number of orbitals within a unit cell of a tight-binding TBHamiltonian.
"""
function norbs(sys::TBHamiltonian)
    return length(sys.orbitals)
end

"""
    hopping_ranges(sys::TBHamiltonian)

Returns a vector with the range of hoppings in the TBHamiltonian, i.e., for a list of hopping in 3D along the directions n*a1 + m*a2+ k*a3, returns maximum(n), maximum(m), maximum(k)
"""
function hopping_ranges(sys::TBHamiltonian)

    dim = dimension(sys)
    ranges = [[0, 0] for n=1:dim]

    for hop in sys.hoppings
        for i=1:dim
            if ranges[i][1] > hop.dir[i] #min of direction i
                ranges[i][1] = hop.dir[i]
            end
            if ranges[i][2] < hop.dir[i] # max of direction i
                ranges[i][2] = hop.dir[i]
            end
        end
    end

    return [ranges[i][1]:ranges[i][2] for i=1:dim]
end

"""
    hopping_array(sys::TBHamiltonian)

Returns the hoppings t^{to, from}_{n1, n2, ...} store in a array t[n1, n2, .., to, from]. n1, ... run from -range1, ..., range1
"""
function hopping_array(sys::TBHamiltonian)
    dim = dimension(sys)

    ranges = hopping_ranges(sys)
    Norbs = norbs(sys)

    t = OffsetArray{Float64, dim + 2}(ranges..., 1:Norbs, 1:Norbs)
    for i in eachindex(t)
        t[i] = zero(hopping_type(sys))
    end

    for hop in sys.hoppings
        t[hop.dir..., hop.to, hop.from] = hop.val
    end

    return t # array t^{t,j}
end


text:
###########
# This gives an example of Gauss–Seidel
#
###########

using BlockBandedMatrices, BandedMatrices, ArrayLayouts, FillArrays

function finitedifference_2d(n)
    h = 1/n
    D² = BandedMatrix(0 => Fill(-2,n), 1 => Fill(1,n-1), -1 => Fill(1,n-1))/h^2
    D_xx = BandedBlockBandedMatrix(Kron(D², Eye(n)))
    D_yy = BandedBlockBandedMatrix(Kron(Eye(n), D²))
    D_xx + D_yy
end

function _gaussseidel(L, U, b, x=copy(b), y=copy(b), M=5)
    for _=1:M
        @view(y[1:end-1]) .= Mul(U , @view(x[2:end]))
        y[end] = 0
        x .= b .- y
        x .= Ldiv(L, x)
    end
    x
end

function _gaussseidel2(L, U, b, x=copy(b), y=copy(b), M=5)
    for _=1:M
        mul!(@view(y[1:end-1]) , U , @view(x[2:end]))
        y[end] = 0
        y .= b .- y
        x .= L\ y
    end
    x
end

function gaussseidel(A, b, M)
    n = Int(sqrt(length(b)))
    L = LowerTriangular(A)
    U = BandedBlockBandedMatrix(UpperTriangular(@view(A[1:end-1,2:end])), ([fill(n,n-1); n-1], [n-1; fill(n,n-1)]),
                                            (0,1), (0,1))
    x = copy(b)
    _gaussseidel(L,U, b, x, copy(x), M) # 1.6s
    x
end

n = 1000;
    Δt = (1/n^2)/4; @time Δ = finitedifference_2d(n);
    @time A = I - Δt*Δ  # 1m x 1m discretization
    b = randn(n^2)
    L = LowerTriangular(A)
    @time U = BandedBlockBandedMatrix(UpperTriangular(@view(A[1:end-1,2:end])), ([fill(n,n-1); n-1], [n-1; fill(n,n-1)]),
                                            (0,1), (0,1))
    x = copy(b)
    y = copy(x)
    @time _gaussseidel(L,U, b, x, y, 20) # 0.4s

    @time L̃ = sparse(BandedBlockBandedMatrix(L))
    @time Ũ = sparse(U)
    x = copy(b)
    y = copy(x)
    @time _gaussseidel2(L̃,Ũ, b, x, y, 20)

Ã = sparse(A)
@time qr(Ã);
@time Ã \b ;

A = randn(n,n)

n = 1000;
    As = [randn(2k+1,k) for k=1:n];

    @time for A in As
        qr(A)
    end
n = 1000;
    # As = [randn(2k+1,k) for k=1:n];

    @time for A in As
        lu(A)
    end

0.02*n

A*x - b

x̃ = copy(x)
ỹ = copy(y)

@view(y[1:end-1]) .= Mul(U , @view(x[2:end]))

mul!(@view(ỹ[1:end-1]) , Ũ , @view(x̃[2:end]))
y[end] = 0
ỹ[end] = 0
x = copy(b)
    y = copy(x)
    _gaussseidel(L,U, b, x, y, 1) - x_bb

x .= b .- y
ỹ .= b .- ỹ
x .= Ldiv(L, x)

x̃ .= L̃ \ ỹ

x - x_bb

x̃-x

(L̃ \ x) - (L \ x)

ỹ - y |> norm
x - x_bb

A*x - b

L̃ == L
Ũ == U

norm(A*x - b)
@time

@time copy(x)
L̃ = sparse(L)
Ũ = sparse(U)

y = x[1:end-1]
z = similar(y)
@time z .= Mul(U,y)
y = similar(x)
@time (@view(y[1:end-1]) .= Mul(U , @view(x[2:end])))
y[end] = 0
@time x .= b .- y
@time x .= Ldiv(L, x)
A = randn(9,40_000);
    @time A*x;

A
S = sparse(A);

400^2

@time L\x

h = 1/n
    @time D² = BandedMatrix(0 => Fill(-2,n), 1 => Fill(1,n-1), -1 => Fill(1,n-1))/h^2
    @time D_xx = BandedBlockBandedMatrix(Kron(D², Eye(n)))
    @time D_yy = BandedBlockBandedMatrix(Kron(Eye(n), D²))
    @time D_xx .+ D_yy
    4

using Profile
dest = Δ;
    @time BlockBandedMatrices.blockbanded_copyto!(dest, D_yy)
    view(dest, Block(1,1)) .= view(dest, Block(1,1)) .+ view(D_yy, Block(1,1))
D_yy
@time copyto!(dest, Broadcast.broadcasted(+, D_xx , D_yy))

using FFTW
A = randn(400, 400);
    @time fft(A);
dest[Block(1,1)]


text:
# low-level wrappers of the CUSPARSE library

#helper functions
function cusparseCreate()
  handle = Ref{cusparseHandle_t}()
  @check ccall( (:cusparseCreate, libcusparse), cusparseStatus_t, (Ptr{cusparseHandle_t},), handle)
  handle[]
end

function cusparseDestroy(handle)
  @check ccall( (:cusparseDestroy, libcusparse), cusparseStatus_t, (cusparseHandle_t,), handle)
end

function cusparseGetVersion(handle, version)
  @check ccall( (:cusparseGetVersion, libcusparse), cusparseStatus_t, (cusparseHandle_t, Ptr{Cint}), handle, version)
end

function cusparseSetStream(handle, streamId)
  @check ccall( (:cusparseSetStream, libcusparse), cusparseStatus_t, (cusparseHandle_t, CuStream_t), handle, streamId)
end

function cusparseGetStream(handle, streamId)
  @check ccall( (:cusparseGetStream, libcusparse), cusparseStatus_t, (cusparseHandle_t, Ptr{CuStream_t}), handle, streamId)
end

function cusparseGetPointerMode(handle, mode)
  @check ccall( (:cusparseGetPointerMode, libcusparse), cusparseStatus_t, (cusparseHandle_t, Ptr{cusparsePointerMode_t}), handle, mode)
end

function cusparseSetPointerMode(handle, mode)
  @check ccall( (:cusparseSetPointerMode, libcusparse), cusparseStatus_t, (cusparseHandle_t, cusparsePointerMode_t), handle, mode)
end

function cusparseCreateHybMat(hybA)
  @check ccall( (:cusparseCreateHybMat, libcusparse), cusparseStatus_t, (Ptr{cusparseHybMat_t},), hybA)
end

function cusparseDestroyHybMat(hybA)
  @check ccall( (:cusparseDestroyHybMat, libcusparse), cusparseStatus_t, (cusparseHybMat_t,), hybA)
end

function cusparseCreateSolveAnalysisInfo(info)
  @check ccall( (:cusparseCreateSolveAnalysisInfo, libcusparse), cusparseStatus_t, (Ptr{cusparseSolveAnalysisInfo_t},), info)
end

function cusparseDestroySolveAnalysisInfo(info)
  @check ccall( (:cusparseDestroySolveAnalysisInfo, libcusparse), cusparseStatus_t, (cusparseSolveAnalysisInfo_t,), info)
end

function cusparseCreateBsrsm2Info(info)
  @check ccall( (:cusparseCreateBsrsm2Info, libcusparse), cusparseStatus_t, (Ptr{bsrsm2Info_t},), info)
end

function cusparseDestroyBsrsm2Info(info)
  @check ccall( (:cusparseDestroyBsrsm2Info, libcusparse), cusparseStatus_t, (bsrsm2Info_t,), info)
end

function cusparseCreateBsrsv2Info(info)
  @check ccall( (:cusparseCreateBsrsv2Info, libcusparse), cusparseStatus_t, (Ptr{bsrsv2Info_t},), info)
end

function cusparseDestroyBsrsv2Info(info)
  @check ccall( (:cusparseDestroyBsrsv2Info, libcusparse), cusparseStatus_t, (bsrsv2Info_t,), info)
end

function cusparseCreateCsrsv2Info(info)
  @check ccall( (:cusparseCreateCsrsv2Info, libcusparse), cusparseStatus_t, (Ptr{csrsv2Info_t},), info)
end

function cusparseDestroyCsrsv2Info(info)
  @check ccall( (:cusparseDestroyCsrsv2Info, libcusparse), cusparseStatus_t, (csrsv2Info_t,), info)
end

function cusparseCreateCsric02Info(info)
  @check ccall( (:cusparseCreateCsric02Info, libcusparse), cusparseStatus_t, (Ptr{csric02Info_t},), info)
end

function cusparseDestroyCsric02Info(info)
  @check ccall( (:cusparseDestroyCsric02Info, libcusparse), cusparseStatus_t, (csric02Info_t,), info)
end

function cusparseCreateCsrilu02Info(info)
  @check ccall( (:cusparseCreateCsrilu02Info, libcusparse), cusparseStatus_t, (Ptr{csrilu02Info_t},), info)
end

function cusparseDestroyCsrilu02Info(info)
  @check ccall( (:cusparseDestroyCsrilu02Info, libcusparse), cusparseStatus_t, (csrilu02Info_t,), info)
end

function cusparseCreateBsric02Info(info)
  @check ccall( (:cusparseCreateBsric02Info, libcusparse), cusparseStatus_t, (Ptr{bsric02Info_t},), info)
end

function cusparseDestroyBsric02Info(info)
  @check ccall( (:cusparseDestroyBsric02Info, libcusparse), cusparseStatus_t, (bsric02Info_t,), info)
end

function cusparseCreateBsrilu02Info(info)
  @check ccall( (:cusparseCreateBsrilu02Info, libcusparse), cusparseStatus_t, (Ptr{bsrilu02Info_t},), info)
end

function cusparseDestroyBsrilu02Info(info)
  @check ccall( (:cusparseDestroyBsrilu02Info, libcusparse), cusparseStatus_t, (bsrilu02Info_t,), info)
end

function cusparseGetProperty(property::CUDAapi.libraryPropertyType)
  value_ref = Ref{Cint}()
  @check ccall((:cusparseGetProperty, libcusparse),
               cusparseStatus_t,
               (Cint, Ptr{Cint}),
               property, value_ref)
  value_ref[]
end


text:
using Documenter, BOMoD

makedocs(sitename="BoMoD.jl",
        authors = "Sam Dierickx")

deploydocs(
        repo = "https://github.com/SamDier/BOMoD.jl",
        )


text:
@testset "conversion" begin
    img = rand(Gray{Float32}, 4, 4)
    img_sp = [SuperPixel(img, (1:2, 1:2)),
              SuperPixel(img, (1:2, 3:4)),
              SuperPixel(img, (3:4, 1:2)),
              SuperPixel(img, (3:4, 3:4))]
    @test img == synthesize(img_sp)
    @test img == synthesize(img_sp, Raw())
    @test img == synthesize(Gray{Float32}, img_sp)
    @test img == synthesize(Gray{Float32}, img_sp, Raw())
    @test RGB{Float32}.(img) == synthesize(RGB{Float32}, img_sp)

    # order independent
    img_sp = [SuperPixel(img, (1:2, 1:2)),
              SuperPixel(img, (3:4, 1:2)),
              SuperPixel(img, (1:2, 3:4)),
              SuperPixel(img, (3:4, 3:4))]
    @test img == synthesize(img_sp, Raw())
    @test_nowarn synthesize(img_sp, Average())

    # overlap
    img_sp = [SuperPixel(img, (1:2, 1:2)),
              SuperPixel(img, (1:2, 3:4)),
              SuperPixel(img, (2:4, 1:3)),
              SuperPixel(img, (3:4, 3:4))]
    @test img == synthesize(img_sp, Raw())
    @test_nowarn synthesize(img_sp, Average())
end


text:
module MultiFieldSparseMatrixAssemblersTests

using Gridap.Arrays
using Gridap.Algebra
using Gridap.Geometry
using Gridap.FESpaces
using Gridap.Fields
using Gridap.Integration
using SparseArrays
using Gridap.CellData
using Gridap.MultiField

using Test

order = 2

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

degree = order
trian = get_triangulation(model)
quad = CellQuadrature(trian,degree)

trian_Γ = SkeletonTriangulation(model)
quad_Γ = CellQuadrature(trian_Γ,degree)

V = TestFESpace(model=model,order=order,reffe=:Lagrangian,conformity=:H1,valuetype=Float64)
Q = TestFESpace(model=model,order=order-1,reffe=:Lagrangian,conformity=:L2,valuetype=Float64)

U = TrialFESpace(V)
P = TrialFESpace(Q)

Y = MultiFieldFESpace([V,Q])
X = MultiFieldFESpace([U,P])

free_values = rand(num_free_dofs(X))
xh = FEFunction(X,free_values)
uh, ph = xh

dy = get_cell_basis(Y)
dx = get_cell_basis(X)
dv, dq = dy
du, dp = dx

cellmat = integrate(dv*du,trian,quad)
cellvec = integrate(dv*2,trian,quad)
cellids = get_cell_id(trian)
cellmatvec = pair_arrays(cellmat,cellvec)

dv_Γ, dq_Γ = restrict(get_cell_basis(Y), trian_Γ)
du_Γ, dp_Γ = restrict(get_cell_basis(X), trian_Γ)
cellmat_Γ = integrate(  jump(dv_Γ)*dp_Γ.⁺ + mean(dq_Γ)*jump(dp_Γ), trian_Γ,quad_Γ)
cellvec_Γ = integrate(  jump(dv_Γ) + mean(dq_Γ), trian_Γ,quad_Γ)
cellmatvec_Γ = pair_arrays(cellmat_Γ,cellvec_Γ)
cellids_Γ = get_cell_id(trian_Γ)

assem = SparseMatrixAssembler(SparseMatrixCSR{0},X,Y)

matvecdata = ([cellmatvec,cellmatvec_Γ],[cellids,cellids_Γ],[cellids,cellids_Γ])
matdata = ([cellmat,cellmat_Γ],[cellids,cellids_Γ],[cellids,cellids_Γ])
vecdata = ([cellvec,cellvec_Γ],[cellids,cellids_Γ])
data = (matvecdata,matdata,vecdata)

test_assembler(assem,matdata,vecdata,data)

A = assemble_matrix(assem,matdata)

A = allocate_matrix(assem,matdata)

assem = SparseMatrixAssembler(X,Y)
test_assembler(assem,matdata,vecdata,data)

end # module


text:
module Pmm
using ...Types, ...Utils
using DifferentialEquations
#include("aerodynamic_coeff_fragment.jl")
export timeOfFlight!, trajectoryPMM!

    function pmm(du,u,p,t)
        A,proj,coeffAero = p
        v = sqrt(u[3]^2 + u[4]^2)
        ma = machNumber(v,u[2])
        Cd = coeffAero(proj,ma)
        ρ = density(u[2])
        θ = atan(u[4]/u[3])
        #A = projectedArea()

        du[1] = u[3]#vx
        du[2] = u[4]#vy
        #du[3] =  -1/2*pi*calibre^2/4*Cd*ρ*v^2*cos(θ)/m
        du[3] =  -A*Cd*ρ*v^2*cos(θ)/m
        #du[4] = gravity(u[2])-1/2*pi*calibre^2/4*Cd*ρ*v^2*sin(θ)/m
        du[4] = gravity(u[2])-A*Cd*ρ*v^2*sin(θ)/m
    end

    function timeOfFlight!(p::Projectile2D, target::Target2D) #PmmZEroDrag
        #t = 0.0
        range = target.position[1]
        flight_time = 0.0

    x0 = p.position[1]
    y0 = p.position[2]
    v0x = p.velocity_m[1]
    v0y = p.velocity_m[2]
    u0 = [x0,y0,v0x,v0y]

    tspan = (0.0,tmax)
    prob = ODEProblem(pmm,u0,tspan)
    #sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
    integrator = init(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
    oldx = integrator.u[1]
    oldy = integrator.u[2]
    oldu = integrator.u[3]
    oldv = integrator.u[4]
    while range > p.position[1]
        oldx = integrator.u[1]
        oldy = integrator.u[2]
        oldu = integrator.u[3]
        oldv = integrator.u[4]
    step!(integrator,dt,true)
    p.position[1]= integrator.u[1]
    #println(integrator.t, " ",p.position[1], " ", range)
    end
    p.position[1] = oldx
    p.position[2] = oldy
    p.velocity_m[1] = oldu
    p.velocity_m[2] = oldv
    flight_time = integrator.t-dt
    p.yaw = atan(p.velocity_m[2]/p.velocity_m[1])

    return flight_time
end

function condition(u,t,integrator) # Event when event_f(u,t) == 0
  #u[2]
  #target.x-u[1]
  dist-u[1]
end

affect!(integrator) = terminate!(integrator)

function trajectoryPMM!(u0, tspan, p, proj::AbstractPenetrator, target::Target2D)
    global dist = target.position[1]
    prob = ODEProblem(pmm,u0,tspan, p)
#cb = DiscreteCallback(condition,affect!)
#println(target.x)
    cb = ContinuousCallback(condition,affect!)
    sol = solve(prob,Tsit5(),callback=cb)
#sol = solve(prob,Tsit5(),callback=cb)
#sol = solve(prob,BS5(),callback=cb)
#sol = solve(prob,AutoTsit5(Rosenbrock23()),callback=cb)
 proj.position[1] = sol.u[end][1]
 proj.position[2] = sol.u[end][2]
 #proj.z = sol.u[end][3]
 proj.velocity[1] = sol.u[end][4]
 proj.velocity[2] = sol.u[end][5]
 #proj.vz = sol.u[end][6]
 proj.tof = sol.t[end]
end

end


text:
using PkgBenchmark

results = benchmarkpkg("CovarianceEstimation")
export_markdown("benchmark/results.md", results)


text:
# gco types
const GCObool = Cuchar
const GCOnode_id = Cint
const GCOLabelID = Cint
const GCOEnergyTermType = Cint
const GCOEnergyType = Union{Cint, Clonglong}
const GCOVarID = GCOnode_id
const GCOSiteID = GCOVarID

# gco structs
struct GCOSparseDataCost
    site::GCOSiteID
    cost::GCOEnergyTermType
end


text:
module QuandlAccess

export Quandl, TimeSeries, Table
export eq, gt, lt, gte, lte

using CSV
using DataFrames
using HTTP

using Dates

include("types.jl")
include("api.jl")

end