text:
facts("multi-path graphs 1") do
    a = Signal(0)
    b = Signal(0)

    c = map(+, a, b)
    d = merge(a, b)
    e = map(+, a, map(x->2x, a)) # Both depend on a
    f = map(+, a, b, c, e)

    @fact queue_size() --> 0

    for (av,bv) in [(1,2),(1,3),(7,7)]
        @show av bv
        push!(a, av)
        push!(b, bv-1)
        Reactive.run_till_now()
        push!(b, bv)
        Reactive.run_till_now()
        @fact value(c) --> av + bv
        @fact value(e) --> 3av
        @fact value(d) --> bv
        @fact value(f) --> 5av + 2bv
    end

    xv = 2
    x = Signal(xv)
    y = map(identity, x)
    x2 = map(x->2x, x)
    z = map(+, y, x2)
    Reactive.run_till_now()
    @fact value(y) --> xv
    @fact value(x2) --> 2xv
    @fact value(z) --> 3xv

    xv2 = xv + 1
    push!(x, xv2)
    Reactive.run_till_now()
    @fact value(z) --> 3xv2
end

facts("multi-path graphs 2") do
    a = Signal(0)
    b = Signal(0)

    c = map(+, a, b)
    d = merge(a, b)
    e = map(+, a, map(x->2x, a)) # Both depend on a
    f = map(+, a, b, c, e)

    for (av,bv) in [(1,2),(1,3),(7,7)]
        push!(a, av)
        push!(b, bv-1)
        Reactive.run_till_now()
        push!(b, bv)
        Reactive.run_till_now()
        @fact value(c) --> av + bv
        @fact value(e) --> 3av
        @fact value(d) --> bv
        @fact value(f) --> 5av + 2bv
    end

    xv = 2
    x = Signal(xv)
    y = map(identity, x)
    x2 = map(x->2x, x)
    z = map(+, y, x2)
    Reactive.run_till_now()
    @fact value(y) --> xv
    @fact value(x2) --> 2xv
    @fact value(z) --> 3xv

    xv2 = xv + 1
    push!(x, xv2)
    Reactive.run_till_now()
    @fact value(z) --> 3xv2
end

facts("multi-path graphs: dfs good, bfs bad") do
    # DFS good, BFS bad
    # s4x is initially 8 (4*2), after push!(sx,3), s4x should be 12 (4*3), but is instead 9.
    # initially: sx is 2, s2x is 4, s3x is 6, s4x is 8
    # after push!(sx, 3), what happens is:
    # BFS
        # order: sx, s2x, s4x, s3x (bad)
        # sx is updated to 3, s2x updates correctly (to 6), but then s4x is next in the queue,
        # and sees s3x's old value of 6, adds it to sx (3), and gets 9. Finally s3x updates
        # correctly to 9, but it's all too late
    # DFS
        # order: sx, s2x, s3x, s4x (good)
    xv = 2
    sx = Signal(xv; name="sx")
    s2x = map(x->2x, sx; name="s2x")
    s3x = map(x->x + x÷2, s2x; name="s3x")
    s4x = map(+, sx, s3x; name="s4x")
    @fact value.([sx, s2x, s3x, s4x]) --> [2, 4, 6, 8]
    push!(sx, 3)
    Reactive.run_till_now()
    @fact value.([sx, s2x, s3x, s4x]) --> [3, 6, 9, 12]
end

facts("multi-path graphs: bfs good, dfs bad") do
    #BFS good, DFS bad
    # s3x is initially 6 (3*2), after push!(sx, 3), s3x should be 9 (3*3), but is instead 7.
    # initially: sx and s1x1, s1x2 are 2, s2x is 4, s3x is 6
    # after push!(sx, 3), what happens is:
    # DFS:
        # order sx, s1x1, s3x, s1x2 (bad)
        # sx is updated to 3, s1x1 updates correctly (to 3), but then s3x is next in the queue,
        # and sees s1x2's old value of 2, adds it to sx (3), s1x1(3) and gets 8. Finally s1x2 updates
        # correctly to 3, but it's all too late
    # BFS:
        # order sx, s1x1, s1x2, sx3 (good)
    xv = 2
    sx = Signal(xv)
    s1x1 = map(identity, sx)
    s1x2 = map(identity, sx)
    s3x = map(+, s1x1, s1x2, sx)
    @fact value.([sx, s1x1, s1x2, s3x]) --> [2, 2, 2, 6]
    push!(sx, 3)
    Reactive.run_till_now()
    @fact value.([sx, s1x1, s1x2, s3x]) --> [3, 3, 3, 9]
end

facts("multi-path graphs: dfs bad, bfs bad") do
    #BFS bad, dfs bad
    # s3x is initially 6 (3*2), after push!(sx, 3), s3x should be 9 (3*3), but is instead 7.
    # what happens in BFS is:
    # initially: sx and s1x1, s1x2 are 2, s2x is 4, s3x is 6
    # after push!(sx, 3)
    # BFS:
        # order is x, s1x1, s1x2, s3x, s2x (bad)
        # correct: x, s1x1, s1x2, s2x, s3x)
        # sx is updated to 3, s1x1 and s1x2 update correctly (to 3), but s3x is next in the queue,
        # and sees s2x's old value of 4, adds it to sx (3), and gets 7. Finally s2x updates
        # correctly to 6 (2*3), but it's all too late
    # DFS:
        # order is x, s1x1, s2x, s3x, s1x2 (bad)
        # correct: x, s1x1, s1x2, s2x, s3x)
        # sx is updated to 3, s1x1 updates correctly to 3, but then:
        # s2x sees s1x2 as 2 and updates to 5 (3 + 2) - but should be 6, then
        # s3x sees s2x as 5, adds sx (3) and gets 8
        # s1x2 then updates to 3, but that's it, s2x is 5 (not 6), s3x is 8 (not 9)

    xv = 2
    sx = Signal(xv)
    s1x1 = map(identity, sx)
    s1x2 = map(identity, sx)
    s2x = map(+, s1x1, s1x2)
    s3x = map(+, sx, s2x)
    @fact value.([sx, s1x1, s1x2, s2x, s3x]) --> [2, 2, 2, 4, 6]
    push!(sx, 3)
    Reactive.run_till_now()
    @fact value.([sx, s1x1, s1x2, s2x, s3x]) --> [3, 3, 3, 6, 9]
end


text:
function real_time(Data)

    #return a power simulation model loaded with a power model

end

ext: jl


text:
global const jlConstantRegex 	= r"const (GL_)?([[:print:]]+)=([[:print:]]+)"
global const jlTypeAliasRegex = r"typealias ([[:graph:]]+) ([[:graph:]]+)"
global const jlFunctionRegex 	= r"@getCFun \"[[:graph:]]+\" [[:graph:]]+ ([[:graph:]]+)(\([[:print:]]*\)::[[:graph:]]+)"

function parseJLGLFile(path::ASCIIString)
	types 		= Dict{ASCIIString, ASCIIString}()
	functions 	= Dict{ASCIIString, ASCIIString}()
	constants 	= Dict{ASCIIString, ASCIIString}()

	fstream = open(path)
	lines = readlines(fstream)
	close(fstream)

	for line in lines
		@match line begin
			jlTypeAliasRegex(alias, _type) 		=> types[strip(alias)] 			= strip(_type)
			jlFunctionRegex(name, func) 		=> functions[strip(name)] 		= strip(func)
			jlConstantRegex(GL_, name, value)	=> constants["GL_"*strip(name)] = strip(value)
		end
	end
	lines = 0
	types, functions, constants
end


text:
# functiosn related to data about faces (edges in 2D)

"""
  This function gets the information about boundary and interior faces
  and stores it to the mesh object.  In particular, the following fields of
  the mesh are allocated and populatedby this function:
    numBoundaryFaces
    numInterfaces
    numPeriodicInterfaces
    numBC
    bndry_geo_nums
    bndryfaces
    interfaces

  The following fields are allocated but not populated:
    bndry_funcs
    bndry_funcs_revm

  Because the curvilinear metric calculation requires knowing about *all*
  the mesh faces (to calculate Eone), this function creates a new boundary
  condition and assigns all faces not used by one of the existing BCs to it.
  The BC_name for this new BC is defaultBC, which allows the solver to do the
  face integral for faces that don't have a regular BC.  The new boundary
  condition is only created if there are faces not previously assigned to a BC.

  Inputs/Outputs:
    mesh: a mesh object
    opts: the options dictionary, which contains BC info and is updated with the
           new BC info if needed

  Outputs:
    bndry_nums: array of [element number, global face number] for each
                face on a geometric face.  This is returned to make testing
                easier.
"""
function getAllFaceData(mesh::PumiMesh, opts)

  nbndryfaces, ninterfaces, npinterfaces, geo_edge_nums = countBoundaryEdges(mesh)
  mesh.numBoundaryFaces = nbndryfaces
  mesh.numInterfaces = ninterfaces
  mesh.numPeriodicInterfaces = npinterfaces

  if typeof(mesh) <: PumiMeshCG
    @assert mesh.numPeriodicInterfaces == 0
  end

  unused_geo_edge_nums = popBCEdges(geo_edge_nums, opts)
  # create an additional BC if needed to make sure the face integral gets
  # done for the unused edges
  add_bc = length(unused_geo_edge_nums) != 0
  if add_bc
    numBC = opts["numBC"] + 1

    opts["numBC"] = numBC
    opts[string("BC", numBC)] = unused_geo_edge_nums
    opts[string("BC", numBC, "_name")] = "defaultBC"

    if opts["error_undefined_bc"]
      error("dimension $(mesh.dim-1) geometric entities $unused_geo_edge_nums do not have assigned boundary condition")
    end
  end

  # populate mesh.bndry_faces from options dictionary
  mesh.numBC = opts["numBC"]
  allocateBoundaryAndInterfaceArrays(mesh, opts)
  boundary_nums = Array{Int}(mesh.numBoundaryFaces, 2)

  # get the element number, global face number of the mesh edges that have BCs
  offset = 1
  for i=1:mesh.numBC
    key_i = string("BC", i)
    model_edges = opts[key_i]
    # record geometric edges
    ngeo = length(model_edges)
    mesh.bndry_geo_nums[i] = Array{Int}(ngeo)
    copy!(mesh.bndry_geo_nums[i], model_edges)
    mesh.bndry_offsets[i] = offset
    offset, print_warning = getMeshEdgesFromModel(mesh, model_edges, offset, boundary_nums)  # get the mesh edges on the model edge
    if !(add_bc && i == mesh.numBC)  # don't error for the added BC
      if print_warning
        throw(ErrorException("Cannot apply boundary conditions to periodic boundary, model entity $model_edges"))
      end
    end
    # offset is incremented by getMeshEdgesFromModel
  end

  mesh.bndry_offsets[mesh.numBC + 1] = offset # = num boundary edges

  # get boundary information for entire mesh
  getBoundaryArray(mesh, boundary_nums)

  # get interface info for the entire mesh
  getInterfaceArray(mesh)
  sort!(mesh.interfaces)  #TODO: see if this is actually a good idea

  return boundary_nums
end

"""
  Smart allocator for the following fields of the mesh:

    bndry_offsets
    bndry_funcs
    bndry_funcs_revm
    bndry_geo_nums
    bndryfaces
    interfaces

  Because there is no reaonsable default value for bndry_funcs and
  bndry_funcs_revm, these fields are not modified if this function is called
  repeatedly.
"""
function allocateBoundaryAndInterfaceArrays(mesh, opts)

  if !isFieldDefined(mesh, :bndry_offsets, :bndry_funcs, :bndry_funcs_revm,
                     :bndry_geo_nums, :bndryfaces, :interfaces)

    mesh.bndry_offsets = Array{Int}(mesh.numBC + 1)
    mesh.bndry_funcs = Array{BCType}(mesh.numBC)
    mesh.bndry_funcs_revm = Array{BCType_revm}(mesh.numBC)
    mesh.bndry_funcs_revq = Array{BCType_revq}(mesh.numBC)
    mesh.bndry_geo_nums = Array{Array{Int, 1}}(mesh.numBC)

    mesh.bndryfaces = Array{Boundary}(mesh.numBoundaryFaces)
    mesh.interfaces = Array{Interface}(mesh.numInterfaces)
  else
    fill!(mesh.bndry_offsets, 0)
    fill!(mesh.bndry_geo_nums, 0)
    fill!(mesh.bndryfaces, Boundary(0, 0))
    fill!(mesh.interfaces, Interface(0, 0, 0, 0, 0))
  end

  return nothing
end

"""
  This function returns an array containing only the geometric edge numbers
  that do not have boundary conditions applied to them

  Inputs:
    geo_edge_nums: vector containing the numbers of all the geometric edges
                   that have non-periodic mesh edges on them
    opts: the options dictionary

  Outputs:
    unused_geo_edge_nums: numbers of the geometric edges (with mesh edges on
                          them) that do not have
                          boundary conditions applied to them

"""
function popBCEdges(geo_edge_nums::AbstractArray{I, 1}, opts) where I <: Integer

  numBC = opts["numBC"]

  # create array of all model edges that have a boundary condition
  bndry_edges_BC = Array{Int}(0)
  for i=1:numBC
    key_i = string("BC", i)
    bndry_edges_i = opts[key_i]
    for j=1:length(bndry_edges_i)
      push!(bndry_edges_BC, bndry_edges_i[j])
    end
  end

  # find the (first) duplicate BC, throw error
  for i=1:length(bndry_edges_BC)
    for j=(i+1):length(bndry_edges_BC)
      if bndry_edges_BC[i] == bndry_edges_BC[j]
        dup_edge = bndry_edges_BC[i]
        #TODO; find which BCs it is assigned to for an even better error
        #       message
        throw(ErrorException("Only 1 boundary condition may be applied to each geometric edge/face, edge/face $dup_edge appears twice"))
      end
    end
  end

  unused_geo_edge_nums = Array{Int}(0)

  for i=1:length(geo_edge_nums)
    edge_i = geo_edge_nums[i]
    if !(edge_i in bndry_edges_BC)
      push!(unused_geo_edge_nums, edge_i)
      continue
    end
  end

  @assert unused_geo_edge_nums == unique(unused_geo_edge_nums)

  return unused_geo_edge_nums
end

"""
  This function populates the mesh.bndryfaces field.

  Inputs:
    mesh: a mesh object
    boundary_nums: the 2 x mesh.numBoundaryFaces array containing the
                   [element number, global face number] for each
                   boundary faces

  Notes:
    the local face number is found using the *Pumi* topology.  It really should
    use the SBP topology, but they are the same for this purpose, so it
    hasn't caused a problem (yet)
"""
function getBoundaryArray(mesh::PumiMesh, boundary_nums::AbstractArray{Int, 2})
# get an array of type Boundary for SBP
# bnd_array is a vector of type Boundary with length equal to the number of edges on the boundary of the mesh

#  mesh.bndryfaces = Array{Boundary}(mesh.numBoundaryFaces)

  for i=1:mesh.numBoundaryFaces
    facenum = boundary_nums[i, 1]
    edgenum_global = boundary_nums[i, 2]
    edgenum_local = getBoundaryFaceLocalNum(mesh, edgenum_global)
    mesh.bndryfaces[i] = Boundary(facenum, edgenum_local)
  end

  return nothing
end

# deprecated
function getBoundaryFaceNormals(mesh::PumiMesh2D, sbp::AbstractSBP, bndry_faces::AbstractArray{Boundary, 1}, face_normals::Array{Tmsh, 3}) where Tmsh

  nfaces = length(bndry_faces)

  alpha = zeros(Tmsh, 2,2)
  dxidx = mesh.dxidx
  jac = mesh.jac
  for i=1:nfaces
    element_i = bndry_faces[i].element
    face_i = bndry_faces[i].face
    for j=1:sbp.numfacenodes
      node_index = sbp.facenodes[j, face_i]

      # calculate the mysterious parameter alpha
      for di1=1:2
	for di2=1:2
	  alpha[di1, di2] = (dxidx[di1,1, node_index, element_i].*dxidx[di2,1, node_index, element_i] + dxidx[di1,2, node_index, element_i].*dxidx[di2,2, node_index, element_i])*jac[node_index,element_i]
	end
      end

    # call SBP function
    smallmatvec!(alpha, sview(sbp.facenormal, :, face_i), sview(face_normals, :, j, i))

    end  # end loop over face nodes
  end  # end loop over faces

  return nothing
end

# deprecated
function getInternalFaceNormals(mesh::PumiMesh2D, sbp::AbstractSBP, internal_faces::AbstractArray{Interface, 1}, face_normals::Array{Tmsh, 4}) where Tmsh

  nfaces = length(internal_faces)

  alpha = zeros(Tmsh, 2,2)
  dxidx = mesh.dxidx
  jac = mesh.jac
  for i=1:nfaces
    element_iL = internal_faces[i].elementL
    face_iL = internal_faces[i].faceL
    element_iR = internal_faces[i].elementR
    face_iR = internal_faces[i].faceR
    for j=1:sbp.numfacenodes

      # calculate left face normal
      node_index = sbp.facenodes[j, face_iL]

      # calculate the mysterious parameter alpha
      for di1=1:2
	for di2=1:2
	  alpha[di1, di2] = (dxidx[di1,1, node_index, element_iL].*dxidx[di2,1, node_index, element_iL] + dxidx[di1,2, node_index, element_iL].*dxidx[di2,2, node_index, element_iL])*jac[node_index,element_iL]
	end
      end

    # call SBP function
    smallmatvec!(alpha, sview(sbp.facenormal, :, face_iL), sview(face_normals, :, 1, j, i))
      # calculate right fae normal
      node_index = sbp.facenodes[j, face_iR]

      # calculate the mysterious parameter alpha
      for di1=1:2
	for di2=1:2
	  alpha[di1, di2] = (dxidx[di1,1, node_index, element_iR].*dxidx[di2,1, node_index, element_iR] + dxidx[di1,2, node_index, element_iR].*dxidx[di2,2, node_index, element_iR])*jac[node_index,element_iR]
	end
      end

    # call SBP function
    smallmatvec!(alpha, sview(sbp.facenormal, :, face_iR), sview(face_normals, :, 2, j, i))

    end  # end loop over face nodes
  end  # end loop over faces

  return nothing
end

function getInterfaceArray(mesh::PumiMesh2D)
# get array of [elementL, elementR, edgeL, edgeR] for each internal edge,
# where elementL and R are the elements that use the edge, edgeL R are the
# local number of the edge within the element
# interfaces is the array to be populated with the data, it must be
# number of internal edges by 1
# the nodemap is used to determine the node ordering of elementR relative to
# elementL (ie. the relative orientation is determined between the nodes of
# the 2 element in SBP order)
# the edge is "owned" by elementL, ie the first node on the edge is the first
# node from the perspective of elementL

  # only need internal boundaries (not external)
  # get number of nodes affecting an edge
  num_edge_nodes = apf.countAllNodes(mesh.mshape_ptr, 1)

  # unused variable?
  nodemap = Array(num_edge_nodes:(-1):1)
  part_nums = Array{Cint}(1)
  matched_entities = Array{Ptr{Void}}(1)
  seen_entities = Set{Ptr{Void}}()
  sizehint!(seen_entities, mesh.numPeriodicInterfaces)

  pos = 1 # current position in interfaces
  for i=1:mesh.numEdge
    edge_i = mesh.edges[i]
    if in(edge_i, seen_entities)
      continue
    end
    # get number of elements using the edge
    adjacent_nums, num_adjacent = getAdjacentEntityNums(mesh, i, 1, 2)
    n = apf.countMatches(mesh.m_ptr, mesh.edges[i])
    apf.getMatches(part_nums, matched_entities)
    has_local_match = part_nums[1] == mesh.myrank && n > 0

    if num_adjacent > 1 || has_local_match # internal edge
      if !has_local_match
        element1 = adjacent_nums[1]
        element2 = adjacent_nums[2]
        edge1 = i
        edge2 = i
      else  # does have local match
        edge2_ptr = matched_entities[1]
        push!(seen_entities, edge2_ptr)
        element1 = adjacent_nums[1]

        # get the element of the other edge
        num_adjacent = apf.countAdjacent(mesh.m_ptr, edge2_ptr, mesh.dim)
        adjacent_entities = apf.getAdjacent(num_adjacent)
        element2 = apf.getNumberJ(mesh.el_Nptr, adjacent_entities[1], 0, 0) + 1
        edge1 = i
        edge2 = apf.getNumberJ(mesh.edge_Nptr, edge2_ptr, 0, 0) + 1
      end

      coords_1 = zeros(3,3)
      coords_2 = zeros(3,3)
      apf.getFaceCoords(mesh.m_ptr, mesh.elements[element1], coords_1, 3, 3)
      apf.getFaceCoords(mesh.m_ptr, mesh.elements[element2], coords_2, 3, 3)

      # calculate centroid
      centroid1 = sum(coords_1, 2)
      centroid2 = sum(coords_2, 2)

      if abs(centroid1[1] - centroid2[1]) > 1e-10  # if big enough difference
          if centroid1[1] < centroid2[1]
          elementL = element1
          elementR = element2
          edgeL = edge1
          edgeR = edge2
          centroidL = centroid1
          centroidR = centroid2
        else
          elementL = element2
          elementR = element1
          edgeL = edge2
          edgeR = edge1
          centroidL = centroid2
          centroidR = centroid1
        end
      else  # use y coordinate to decide
        if centroid1[2] < centroid2[2]
          elementL = element1
          elementR = element2
          edgeL = edge1
          edgeR = edge2
          centroidL = centroid1
          centroidR = centroid2
        else
          elementL = element2
          elementR = element1
          edgeL = edge2
          edgeR = edge1
          centroidL = centroid2
          centroidR = centroid1
        end
      end

      edgeL = getFaceLocalNum(mesh, edgeL, elementL)
      edgeR = getFaceLocalNum(mesh, edgeR, elementR)

      # here we set the orientation flag to 1 for all edges, because in 2D
      # the edges always have opposite orientation
      mesh.interfaces[pos] = Interface(elementL, elementR, edgeL, edgeR, UInt8(1))
      pos += 1

  #    print("\n")

    end  # end if internal edge

#     print("\n")
  end  # end loop over edges

  return nothing

end  # end function

"""
  This function populates the mesh.interfaces array

  Inputs/Outputs:
    mesh: a 3D DG mesh object.  The mesh.interfaces field should already have
          been allocated to be of length mesh.numInterfaces.  This function
          populates it with one Interface object for each face shared by
          two elements, including periodic faces.  The face numbers are the
          local face numbers, and the orientation is determined using the
          SBP topology.

    Notes:
      The local face numbers use the Pumi topology.  This face numbering
      (but not orientation) is consistent between Pumi and the SummationByParts
      package, but this could be a problem if that ever changes.

"""
function getInterfaceArray(mesh::PumiMesh3D)
# Interface.orient is calculated using the *SBP* topology of the faces

  adj_elements = Array{Ptr{Void}}(2)
  coords1 = zeros(3, 4)
  coords2 = zeros(3, 4)
  centroid1 = zeros(3)
  centroid2 = zeros(3)
  vertsL = Array{Ptr{Void}}(4)
  vertsR = Array{Ptr{Void}}(4)
  facevertsL = Array{Ptr{Void}}(3)
  facevertsR = Array{Ptr{Void}}(3)
  # there can be up to 400 elements using a vertex + 8
  # corners of the domain, so in the worst case a parallel partition could
  # generate 400 + 8 matches
  part_nums = Array{Cint}(400 + 8)
  matched_entities = Array{Ptr{Void}}(400 + 8)  # there can be up to 400

  pos = 1 # position in mesh.interfaces
  seen_entities = Set{Ptr{Void}}()
  sizehint!(seen_entities, mesh.numPeriodicInterfaces)

  for i=1:mesh.numFace
    face_i = mesh.faces[i]
    if face_i in seen_entities  # avoid counting matched entities twice
      continue
    end

    num_adjacent = apf.countAdjacent(mesh.m_ptr, face_i, mesh.dim)

    n = apf.countMatches(mesh.m_ptr, face_i)
    apf.getMatches(part_nums, matched_entities)
    has_local_match = part_nums[1] == mesh.myrank && n > 0

    if num_adjacent == 2 || has_local_match  # this is a shared interface

      apf.getAdjacent(adj_elements)
      if has_local_match
        # get the parent element
        el1 = adj_elements[1]
        edgenum1 = i
        elnum1 = apf.getNumberJ(mesh.el_Nptr, el1 , 0, 0) + 1

        # get the matched face
        other_face = matched_entities[1]
        # get the parent element of the matched face
        apf.countAdjacent(mesh.m_ptr, other_face, mesh.dim)
        apf.getAdjacent(adj_elements)
        el2 = adj_elements[1]
        edgenum2 = apf.getNumberJ(mesh.face_Nptr, other_face, 0, 0) + 1
        elnum2 = apf.getNumberJ(mesh.el_Nptr, el2, 0, 0) + 1
        push!(seen_entities, other_face)
      else
        # get both parent elements
        el1 = adj_elements[1]
        el2 = adj_elements[2]
        edgenum1 = i
        edgenum2 = i
        elnum1 = apf.getNumberJ(mesh.el_Nptr, el1, 0, 0) + 1
        elnum2 = apf.getNumberJ(mesh.el_Nptr, el2, 0, 0) + 1
      end

      # decide which one is elementL
      getElementCoords(mesh, el1, coords1)
      getElementCoords(mesh, el2, coords2)

      for j=1:4
        for k=1:3
          centroid1[k] += coords1[k, j]
          centroid2[k] += coords2[k, j]
        end
      end

      flag = getLR(centroid1, centroid2)
      if flag
        coordsL = coords2
        coordsR = coords1
        elnumL = elnum2
        elnumR = elnum1
        edgenumL = edgenum2
        edgenumR = edgenum1
        elL = el2
        elR = el1
      else
        coordsL = coords1
        coordsR = coords2
        elnumL = elnum1
        elnumR = elnum2
        edgenumL = edgenum1
        edgenumR = edgenum2
        elL = el1
        elR = el2
      end

      # Note: this must be consistent with the topo definition of face numbering
      #       (I think)
      localfacenumL = getFaceLocalNum(mesh, edgenumL, elnumL)
      localfacenumR = getFaceLocalNum(mesh, edgenumR, elnumR)
      fdata = FaceData(elnumL, elL, elnumR, elR, localfacenumL, localfacenumR,
                       vertsL, vertsR, facevertsL, facevertsR, part_nums,
                       matched_entities, has_local_match)

      rel_rotate = getRelativeOrientation(fdata, mesh)
      mesh.interfaces[pos] = Interface(elnumL, elnumR, localfacenumL, localfacenumR, UInt8(rel_rotate))
      pos += 1

      fill!(centroid1, 0.0)
      fill!(centroid2, 0.0)

    end  # end if
  end

  return nothing

end

function getLR(coordsL, coordsR)
# return true if coords should be reversed
  tol = 1e-10

  if abs(coordsL[1] - coordsR[1]) > tol
    if coordsL[1] > coordsR[1]
      return false
    else
      return true
    end
  elseif abs(coordsL[2] - coordsR[2]) > tol
    if coordsL[2] > coordsR[2]
      return false
    else
      return true
    end
  else  # if the first two coordinates are not decisive, use the last one
    if coordsL[3] > coordsR[3]
      return false
    else
      return true
    end
  end

  return false  # this should never be reached
end


text:
using AssetRegistry
using Sockets
using WebIO

struct IJuliaConnection <: AbstractConnection
    comm::IJulia.CommManager.Comm
end

function Sockets.send(c::IJuliaConnection, data)
    IJulia.send_comm(c.comm, data)
end

Base.isopen(c::IJuliaConnection) = haskey(IJulia.CommManager.comms, c.comm.id)

WebIO.register_renderable(T::Type, ::Val{:ijulia}) = nothing

function IJulia.CommManager.register_comm(comm::IJulia.CommManager.Comm{:webio_comm}, x)
    conn = IJuliaConnection(comm)
    comm.on_msg = function (msg)
        data = msg.content["data"]
        WebIO.dispatch(conn, data)
    end
end

"""
A "dummy" type that is used to detect whether or not the Jupyter frontend has
the WebIO integration installed correctly.

This works by definining two Base.show methods:
* One for `WEBIO_NODE_MIME` that displays an empty div (which is effectively
  invisible in the browser).
* One for `text/html` which displays an error message and links to
  troubleshooting documentation.

This works since a properly configured Jupyter frontend should have a renderer
installed for `WEBIO_NODE_MIME` that is preferred over `text/html`. If the HTML
content is actually displayed, it means that the WebIO integration is not
correctly installed.
"""
struct _IJuliaInit end

function Base.show(io::IO, m::WEBIO_NODE_MIME, ::_IJuliaInit)
    Base.show(io, m, node(:div))
end

function Base.show(io::IO, m::MIME"text/html", ::_IJuliaInit)
    Base.print(
        io,
        """

The WebIO Jupyter extension was not detected. See the

            WebIO Jupyter integration documentation

        for more information.

        """,
    )
end

function main()
    if !IJulia.inited
        # If IJulia has not been initialized and connected to Jupyter itself,
        # then we have no way to display anything in the notebook and no way
        # to set up comms, so this function cannot run. That's OK, because
        # any IJulia kernels will start up with a fresh process and a fresh
        # copy of WebIO and IJulia.
        return
    end

    # https://github.com/JuliaLang/IJulia.jl/pull/755
    if isdefined(IJulia, :register_jsonmime)
        IJulia.register_jsonmime(WEBIO_NODE_MIME())
    else
        @warn "IJulia doesn't have register_mime; WebIO may not work as expected. Please upgrade to IJulia v1.13.0 or greater."
    end

    # See comment on _IJuliaInit for what this does
    display(_IJuliaInit())
end

WebIO.setup_provider(::Val{:ijulia}) = main() # calling setup_provider(Val(:ijulia)) will display the setup javascript
WebIO.setup(:ijulia)


text:
using Calendar

immutable Entry
	rate::ExchangeRate
	start::CalendarTime
	end_::CalendarTime
end

type ExchangeRateManager
	data::Dict{BigInt, Array{Entry}}
	ExchangeRateManager()=(x=new(Dict{BigInt, Array{Entry}}()); addKnownRates(x); x)
end

# Add an exchange rate
# The given rate is valid between the given dates
# Note:  If two rates are given between the same currencies
# and with overlapping date ranges, the latest one
# added takes precedence during lookup.
function add(ERM::ExchangeRateManager, ER::ExchangeRate, startdate::CalendarTime, enddate::CalendarTime)
	key=hash(ER.source, ER.target)
	if !haskey(ERM.data, key)
		ERM.data[key]=Array(Entry, 0)
	end
	push!(ERM.data[key], Entry(ER, startdate, enddate))
end

function addKnownRates(ERM::ExchangeRateManager)
	# currencies obsoleted by Euro
	add(ERM, ExchangeRate(EURCurrency(), ATSCurrency(), 13.7603), ymd_hms(1999,1,1,0,0,0), today())
	add(ERM, ExchangeRate(EURCurrency(), BEFCurrency(), 40.3399), ymd_hms(1999,1,1,0,0,0), today())
	add(ERM, ExchangeRate(EURCurrency(), DEMCurrency(), 1.95583), ymd_hms(1999,1,1,0,0,0), today())
	add(ERM, ExchangeRate(EURCurrency(), ESPCurrency(), 166.386), ymd_hms(1999,1,1,0,0,0), today())
	add(ERM, ExchangeRate(EURCurrency(), FIMCurrency(), 5.94573), ymd_hms(1999,1,1,0,0,0), today())
	add(ERM, ExchangeRate(EURCurrency(), FRFCurrency(), 6.55957), ymd_hms(1999,1,1,0,0,0), today())
	add(ERM, ExchangeRate(EURCurrency(), GRDCurrency(), 340.750), ymd_hms(2001,1,1,0,0,0), today())
	add(ERM, ExchangeRate(EURCurrency(), IEPCurrency(), 0.787564), ymd_hms(1999,1,1,0,0,0), today())
	add(ERM, ExchangeRate(EURCurrency(), ITLCurrency(), 1936.27), ymd_hms(1999,1,1,0,0,0), today())
	add(ERM, ExchangeRate(EURCurrency(), LUFCurrency(), 40.3399), ymd_hms(1999,1,1,0,0,0), today())
	add(ERM, ExchangeRate(EURCurrency(), NLGCurrency(), 2.20371), ymd_hms(1999,1,1,0,0,0), today())
	add(ERM, ExchangeRate(EURCurrency(), PTECurrency(), 200.482), ymd_hms(1999,1,1,0,0,0), today())
	add(ERM, ExchangeRate(EURCurrency(), CYPCurrency(), 0.585274), ymd_hms(2008,1,1,0,0,0), today())
	add(ERM, ExchangeRate(EURCurrency(), EEKCurrency(), 15.6466), ymd_hms(2011,1,1,0,0,0), today())
	add(ERM, ExchangeRate(EURCurrency(), LVLCurrency(), 1.42287), ymd_hms(2014,1,1,0,0,0), today())
	add(ERM, ExchangeRate(EURCurrency(), MTLCurrency(), 2.32937), ymd_hms(2008,1,1,0,0,0), today())
	add(ERM, ExchangeRate(EURCurrency(), SITCurrency(), 239.64), ymd_hms(2007,1,1,0,0,0), today())
	add(ERM, ExchangeRate(EURCurrency(), SKKCurrency(), 30.1260), ymd_hms(2009,1,1,0,0,0), today())

    # other obsoleted currencies
	add(ERM, ExchangeRate(TRYCurrency(), TRLCurrency(), 1000000.0), ymd_hms(2005,1,1,0,0,0), today())
	add(ERM, ExchangeRate(RONCurrency(), ROLCurrency(), 10000.0), ymd_hms(2005,7,1,0,0,0), today())
	add(ERM, ExchangeRate(PENCurrency(), PEICurrency(), 1000000.0), ymd_hms(1991,7,1,0,0,0), today())
	add(ERM, ExchangeRate(PEICurrency(), PEHCurrency(), 1000.0), ymd_hms(1985,2,1,0,0,0), today())
end

function clear(ERM::ExchangeRateManager)
	ERM=ExchangeRateManager()
end

hash(c1::Currency, c2::Currency)=BigInt(min(c1.numeric, c2.numeric))*1000+BigInt(max(c1.numeric, c2.numeric))
hashes(key::BigInt, c::Currency)=c.numeric==key%1000 || c.numeric==trunc(key/1000)

function fetch(ERM::ExchangeRateManager, source::Currency, target::Currency, date::CalendarTime)
	if !haskey(ERM.data, hash(source, target))
		return nothing
	end

	entries=ERM.data[hash(source, target)]
	for entry in entries
		if date>=entry.start && date<=entry.end_
			return entry.rate
		end
	end
	return nothing
end

function directLookup(ERM::ExchangeRateManager, source::Currency, target::Currency, date::CalendarTime)
	fetchedRate=fetch(ERM, source, target, date)

	if fetchedRate!=nothing
		return fetchedRate
	else
		error("no direct conversion available from $(source.code) to $(target.code) for $date")
	end
end

function smartLookup(ERM::ExchangeRateManager, source::Currency, target::Currency, date::CalendarTime, forbidden::Array{Int, 1})
	fetchedRate=fetch(ERM, source, target, date)

	# direct exchange rates are preferred.
	if fetchedRate!=nothing
		return fetchedRate
	end

	# if none is found, turn to smart lookup. The source currency
    # is forbidden to subsequent lookups in order to avoid cycles.
	push!(forbidden, source.numeric)

	for key in keys(ERM.data)
		# we look for exchange-rate data which involve our source currency
		if hashes(key, source) && length(ERM.data[key])!=0
			# ...whose other currency is not forbidden...
			e=ERM.data[BigInt(key)][1]
			other= source==e.rate.source ? e.rate.target : e.rate.source
			if !in(other.numeric, forbidden)
				# ...and which carries information for the requested date.
				head=fetch(ERM, source, other, date)

				if head!=nothing
					# if we can get to the target from here...
					try
						tail=smartLookup(ERM, other, target, date, forbidden)
						# ..we're done.
						return chain(head, tail)
					catch
						# otherwise, we just discard this rate.
					end
				end
			end
		end
	end

	# if the loop completed, we have no way to return the requested rate.
	error("no conversion available from $(source.code) to $(target.code) for $date")
end

function lookup(ERM::ExchangeRateManager, source::Currency, target::Currency, date::CalendarTime, ERType::Symbol=:Derived)

	if source==target
		return ExchangeRate(source, target, 1.0)
	end

	if date==ymd_hms(0, 0, 0, 0, 0, 0)
		date=today()
	end

	if ERType==:Direct
		return directLookup(source, target, date)
	elseif haskey(list_deprecated, source.code)
		link=list_deprecated[source.code]

		if link==target
			return directLookup(ERM, source, link, date)
		else
			return chain(directLookup(source, link, date), lookup(link, target, date))
		end

	elseif haskey(list_deprecated, target.code)
		link=list_deprecated[target.code]

		if link==source
			return directLookup(ERM, link, target, date)
		else
			return chain(lookup(source, link, date), directLookup(link, target, date))
		end

	else
		return smartLookup(ERM, source, target, date, Array(Int,0))
	end

end

# Singleton design pattern
function getInstance()
	global ERM_Instance
	if !isdefined(Currencies, :ERM_Instance)
		ERM_Instance=ExchangeRateManager()
	end
	return ERM_Instance
end

minDate=ymd_hms(1901, 1, 1, 0, 0, 0)
maxDate=ymd_hms(2199, 12, 31, 0, 0, 0)
add(ER::ExchangeRate, startdate::CalendarTime=minDate, enddate::CalendarTime=maxDate)=add(getInstance(), ER, startdate, enddate)
clear()=clear(getInstance())
lookup(source::Currency, target::Currency, date::CalendarTime=today(), ERType::Symbol=:Derived)=lookup(getInstance(), source, target, date, ERType)


text:
# A googol (10^100) is a massive number: one followed by one-hundred zeros;
# 100^100 is almost unimaginably large: one followed by two-hundred zeros.
# Despite their size, the sum of the digits in each number is only 1.
#
# Considering natural numbers of the form, a^b, where a, b < 100, what is the
# maximum digital sum?

using ProjectEulerSolutions

# Mostly an exercise in constraining the search space since the digits and
# sum built-in functions do most of the work when using BigInts.
function p056solution(n_max::Integer=10)::Integer

    # Assume that max value will be when a^b is many digits long, so select
    # a and b starting values to be 80% of their maxes.
    a_start = floor(BigInt, n_max * 0.8)
    b_start = floor(BigInt, n_max * 0.8)
    max_sum = 0
    for a in a_start:n_max
        for b in b_start:n_max
            max_sum = max(max_sum, sum(digits(a ^ b)))
        end
    end
    return max_sum
end

p056 = Problems.Problem(p056solution)

Problems.benchmark(p056, 100)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
# This file is a part of RRRMC.jl. License is MIT: http://github.com/carlobaldassi/RRRMC.jl/LICENCE.md

# A Real type which is actually an integer in disguise.
# Similar to FixedPoint, but more efficient for
# the purposes of this code.

module DFloats

export DFloat64, MAXDIGITS

const MAXDIGITS = 5

primitive type DFloat64 <: Real 64 end
const dfact = 10^MAXDIGITS

import Base: convert, ==, <, <=, *, /, +, -, round, typemin, show, promote_rule, decompose,
             zero, signbit, abs

i2d(x::Int64) = reinterpret(DFloat64, x)
d2i(x::DFloat64) = reinterpret(Int64, x)
convert(::Type{DFloat64}, x::DFloat64) = x
convert(::Type{DFloat64}, x::Real) = i2d(round(Int64, x * dfact))
convert(::Type{Bool}, x::DFloat64) = d2i(x) > 0
convert(::Type{Integer}, x::DFloat64) = Int64(d2i(x) / dfact)
convert(::Type{T}, x::DFloat64) where {T<:Real} = T(d2i(x) / dfact)
convert(::Type{DFloat64}, x::Integer) = i2d(Int64(x * dfact))
DFloat64(x::Number) = convert(DFloat64, x)

-(x::DFloat64) = i2d(-d2i(x))
-(x::DFloat64, y::DFloat64) = i2d(d2i(x) - d2i(y))
+(x::DFloat64, y::DFloat64) = i2d(d2i(x) + d2i(y))
*(a::Bool, b::DFloat64) = i2d(a * d2i(b))
*(a::DFloat64, b::Bool) = b * a
*(a::Integer, b::DFloat64) = i2d(Int64(a) * d2i(b))
*(a::DFloat64, b::Integer) = b * a
/(a::DFloat64, b::Integer) = i2d(d2i(a) ÷ Int64(b))
==(a::DFloat64, b::DFloat64) = d2i(a) == d2i(b)
<(a::DFloat64, b::DFloat64) = d2i(a) < d2i(b)
<=(a::DFloat64, b::DFloat64) = d2i(a) <= d2i(b)

zero(::Type{DFloat64}) = i2d(Int64(0))

promote_rule(::Type{Float64}, ::Type{DFloat64}) = Float64
promote_rule(::Type{DFloat64}, ::Type{Float64}) = Float64
promote_rule(::Type{Int}, ::Type{DFloat64}) = DFloat64
promote_rule(::Type{DFloat64}, ::Type{Int}) = DFloat64

round(x::DFloat64, digits::Integer) = (@assert digits == MAXDIGITS; x)
typemin(::Type{DFloat64}) = i2d(typemin(Int64))

signbit(x::DFloat64) = signbit(d2i(x))
abs(x::DFloat64) = i2d(abs(d2i(x)))

show(io::IO, x::DFloat64) = show(io, Float64(x))

# used for hashing
function decompose(x::DFloat64)
    i = d2i(x)
    m = gcd(i, dfact)
    return i ÷ m, 1, dfact ÷ m
end

end # module


text:
@doc raw"""
    ps_cost(psi_container::PSIContainer,
                var_name::Symbol,
                index::String,
                cost_component::Float64,
                dt::Float64,
                sign::Float64)

Returns linear cost terms for sum of variables with common factor to be used for cost expression for psi_container model.

# Equation

``` gen_cost = sum(variable)*cost_component ```

# LaTeX

`` cost = dt\times sign\sum_{i\in I} c x_i ``

Returns:

``` sign*gen_cost*dt ```

# Arguments

* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* var_name::Symbol: The variable name
* index::String: The index of the variable container
* cost_component::Float64 : cost to be associated with variable
* dt::Float64 : fraction of hour
* sign::Float64 : positive or negative sign to be associated cost term
"""
function ps_cost(
    psi_container::PSIContainer,
    var_name::Symbol,
    index::String,
    cost_component::Float64,
    dt::Float64,
    sign::Float64,
)
    variable = get_variable(psi_container, var_name)[index, :]
    gen_cost = sum(variable) * cost_component
    return sign * gen_cost * dt
end

@doc raw"""
    ps_cost(psi_container::PSIContainer,
                var_name::Symbol,
                index::String,
                cost_component::PSY.VariableCost{Float64},
                dt::Float64,
                sign::Float64)

Returns linear cost terms for sum of variables with common factor to be used for cost expression for psi_container model.
Does this by calling ```ps_cost``` that has Float64 cost component input.

Returns:

``` ps_cost(psi_container, variable, PSY.get_cost(cost_component), dt, sign) ```

# Arguments

* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* var_name::Symbol: The variable name
* index::String: The index of the variable container
* cost_component::PSY.VariableCost{Float64} : container for cost to be associated with variable
* dt::Float64 : fraction of hour
* sign::Float64 : positive or negative sign to be associated cost term
"""
function ps_cost(
    psi_container::PSIContainer,
    var_name::Symbol,
    index::String,
    cost_component::PSY.VariableCost{Float64},
    dt::Float64,
    sign::Float64,
)
    return ps_cost(psi_container, var_name, index, PSY.get_cost(cost_component), dt, sign)
end

@doc raw"""
    ps_cost(psi_container::PSIContainer,
                var_name::Symbol,
                index::String,
                cost_component::PSY.VariableCost{NTuple{2, Float64}}
                dt::Float64,
                sign::Float64)

Returns quadratic cost terms for sum of variables with common factor to be used for cost expression for psi_container model.

# Equation

``` gen_cost = dt*sign*(sum(variable.^2)*cost_component[1] + sum(variable)*cost_component[2]) ```

# LaTeX

`` cost = dt\times sign (sum_{i\in I} c_1 v_i^2 + sum_{i\in I} c_2 v_i ) ``

for quadratic factor large enough. Otherwise

``` return ps_cost(psi_container, variable, cost_component[2], dt, 1.0) ```

Returns ```gen_cost```

# Arguments

* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* var_name::Symbol: The variable name
* index::String: The index of the variable container
* cost_component::PSY.VariableCost{NTuple{2, Float64}} : container for quadratic and linear factors
* sign::Float64 : positive or negative sign to be associated cost term
"""
function ps_cost(
    psi_container::PSIContainer,
    var_name::Symbol,
    index::String,
    cost_component::PSY.VariableCost{NTuple{2, Float64}},
    dt::Float64,
    sign::Float64,
)
    if cost_component[1] >= eps()
        variable = get_variable(psi_container, var_name)[index, :]
        gen_cost =
            sum(variable .^ 2) * cost_component[1] + sum(variable) * cost_component[2]
        return sign * gen_cost * dt
    else
        return ps_cost(psi_container, var_name, index, cost_component[2], dt, 1.0)
    end
end

@doc raw"""
    _pwlparamcheck(cost_)

Returns True/False depending on compatibility of the cost data with the linear implementation method

Returns ```flag```

# Arguments

* cost_::PSY.VariableCost{NTuple{2, Float64}} : container for quadratic and linear factors
"""
function _pwlparamcheck(cost_)
    flag = true
    slopes = PSY.get_slopes(cost_)
    # First element of the array is the average cost at P_min
    for ix in 2:(length(slopes) - 1)
        if slopes[ix] > slopes[ix + 1]
            @debug slopes
            return flag = false
        end
    end
    return flag
end

@doc raw"""
    _pwlgencost_sos(psi_container::PSIContainer,
                variable::JuMP.Containers.DenseAxisArray{JV},
                cost_component::PSY.VariableCost{NTuple{2, Float64}}) where {JV <: JuMP.AbstractVariableRef}

Returns piecewise cost expression using SOS Type-2 implementation for psi_container model.

# Equations

``` variable = sum(sos_var[i]*cost_component[2][i])```

``` gen_cost = sum(sos_var[i]*cost_component[1][i]) ```

# LaTeX

`` variable = (sum_{i\in I} c_{2, i} sos_i) ``

`` gen_cost = (sum_{i\in I} c_{1, i} sos_i) ``

Returns ```gen_cost```

# Arguments

* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* variable::JuMP.Containers.DenseAxisArray{JV} : variable array
* cost_component::PSY.VariableCost{NTuple{2, Float64}} : container for quadratic and linear factors
"""
function _pwlgencost_sos(
    psi_container::PSIContainer,
    variable::JV,
    cost_component::Vector{NTuple{2, Float64}},
) where {JV <: JuMP.AbstractVariableRef}
    gen_cost = JuMP.GenericAffExpr{Float64, _variable_type(psi_container)}()
    pwlvars = JuMP.@variable(
        psi_container.JuMPmodel,
        [i in 1:length(cost_component)],
        base_name = "{$(variable)}_{sos}",
        start = 0.0,
        lower_bound = 0.0,
        upper_bound = 1.0
    )

    JuMP.@constraint(psi_container.JuMPmodel, sum(pwlvars) == 1.0)
    JuMP.@constraint(
        psi_container.JuMPmodel,
        pwlvars in MOI.SOS2(collect(1:length(pwlvars)))
    )
    for (ix, var) in enumerate(pwlvars)
        JuMP.add_to_expression!(gen_cost, cost_component[ix][1] * var)
    end

    JuMP.@constraint(
        psi_container.JuMPmodel,
        variable == sum([var * cost_component[ix][2] for (ix, var) in enumerate(pwlvars)])
    )

    return gen_cost, pwlvars
end

@doc raw"""
    _pwlgencost_linear(psi_container::PSIContainer,
                variable::JuMP.Containers.DenseAxisArray{JV},
                cost_component::PSY.VariableCost{NTuple{2, Float64}}) where {JV <: JuMP.AbstractVariableRef}

Returns piecewise cost expression using linear implementation for psi_container model.

# Equations

``` 0 <= pwl_var[i] <= (cost_component[2][i] - cost_component[2][i-1])```

``` variable = sum(pwl_var[i])```

``` gen_cost = sum(pwl_var[i]*cost_component[1][i]/cost_component[2][i]) ```

# LaTeX
`` 0 <= pwl_i <= (c_{2, i} - c_{2, i-1})``

`` variable = (sum_{i\in I} pwl_i) ``

`` gen_cost = (sum_{i\in I}  pwl_i) c_{1, i}/c_{2, i} ``

Returns ```gen_cost```

# Arguments

* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* variable::JuMP.Containers.DenseAxisArray{JV} : variable array
* cost_component::PSY.VariableCost{NTuple{2, Float64}} : container for quadratic and linear factors
"""
function _pwlgencost_linear(
    psi_container::PSIContainer,
    variable::JV,
    cost_component::Vector{NTuple{2, Float64}},
) where {JV <: JuMP.AbstractVariableRef}
    gen_cost = JuMP.GenericAffExpr{Float64, _variable_type(psi_container)}()
    pwlvars = JuMP.@variable(
        psi_container.JuMPmodel,
        [i = 1:length(cost_component)],
        base_name = "{$(variable)}_{pwl}",
        start = 0.0,
        lower_bound = 0.0,
        upper_bound = PSY.get_breakpoint_upperbounds(cost_component)[i]
    )

    for (ix, pwlvar) in enumerate(pwlvars)
        JuMP.add_to_expression!(gen_cost, PSY.get_slopes(cost_component)[ix] * pwlvar)
    end
    c = JuMP.@constraint(
        psi_container.JuMPmodel,
        variable == sum([pwlvar for (ix, pwlvar) in enumerate(pwlvars)])
    )
    return gen_cost, pwlvars
end

@doc raw"""
    _pwl_cost(cost)

Returns JuMP expression for a piecewise linear cost function depending on the data compatibility.

Returns ```gen_cost```

# Arguments

* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* variable::JuMP.Containers.DenseAxisArray{JV} : variable array
* cost_component::PSY.VariableCost{NTuple{2, Float64}} : container for quadratic and linear factors
"""
function _pwl_cost(
    psi_container::PSIContainer,
    variable::JV,
    cost_component::Vector{NTuple{2, Float64}},
) where {JV <: JuMP.AbstractVariableRef}
    if !_pwlparamcheck(cost_component)
        @warn("The cost function provided for $(variable) device is not compatible with a linear PWL cost function.
        An SOS-2 formulation will be added to the model.
        This will result in additional binary variables added to the model.")
        gen_cost, vars = _pwlgencost_sos(psi_container, variable, cost_component)
    else
        gen_cost, vars = _pwlgencost_linear(psi_container, variable, cost_component)
    end
    return gen_cost, vars
end

@doc raw"""
    ps_cost(psi_container::PSIContainer,
                var_name::Symbol,
                index::String,
                cost_component::PSY.VariableCost{Vector{NTuple{2, Float64}}},
                dt::Float64,
                sign::Float64)

Creates piecewise linear cost function using a sum of variables and expression with sign and time step included.

# Expression

```JuMP.add_to_expression!(gen_cost, c)```

Returns sign*gen_cost*dt

# LaTeX

``cost = sign\times dt \sum_{v\in V} c_v``

where ``c_v`` is given by

`` c_v = \sum_{i\in Ix} \frac{y_i - y_{i-1}}{x_i - x_{i-1}} v^{p.w.}_i ``

# Arguments

* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* var_name::Symbol: The variable name
* index::String: The index of the variable container
* cost_component::PSY.VariableCost{Vector{NTuple{2, Float64}}}
* dt::Float64 : fraction of hour
* sign::Float64 : positive or negative sign to be associated cost term
"""
function ps_cost(
    psi_container::PSIContainer,
    var_name::Symbol,
    index::String,
    cost_component::PSY.VariableCost{Vector{NTuple{2, Float64}}},
    dt::Float64,
    sign::Float64,
)
    cost_array = cost_component.cost
    # If array is full of tuples with zeros return 0.0
    all(iszero.(last.(cost_array))) && return JuMP.AffExpr(0.0)
    variable = get_variable(psi_container, var_name)[index, :]
    if !haskey(psi_container.variables, :PWL_cost_vars)
        time_steps = model_time_steps(psi_container)
        container = add_var_container!(
            psi_container,
            :PWL_cost_vars,
            [index],
            time_steps,
            1:length(cost_component);
            sparse = true,
        )
    else
        container = get_variable(psi_container, :PWL_cost_vars)
    end
    gen_cost = JuMP.GenericAffExpr{Float64, _variable_type(psi_container)}()
    for (t, var) in enumerate(variable)
        c, pwl_vars = _pwl_cost(psi_container, var, cost_array)
        for (ix, v) in enumerate(pwl_vars)
            container[(index, t, ix)] = v
        end
        JuMP.add_to_expression!(gen_cost, c)
    end

    return sign * gen_cost * dt
end

@doc raw"""
    add_to_cost(psi_container::PSIContainer,
                     devices::D,
                     var_name::Symbol,
                     cost_symbol::Symbol,
                     sign::Float64 = 1.0) where {D<:IS.FlattenIteratorWrapper{<:PSY.Device}}

Adds cost expression for each device using appropriate call to ```ps_cost```.

# Expression

for d in devices

```    cost_expression = ps_cost(psi_container,
                              variable[PSY.get_name(d), :],
                              getfield(PSY.get_op_cost(d), cost_symbol),
                              dt,
                              sign) ```
``` psi_container.cost_function += cost_expression ```

# LaTeX

`` COST = \sum_{d\in D} cost_d ``

# Arguments

* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* devices::D : set of devices
* var_name::Symbol : name of variable
* cost_symbol::Symbol : symbol associated with costx
"""
function add_to_cost(
    psi_container::PSIContainer,
    devices::D,
    var_name::Symbol,
    cost_symbol::Symbol,
    sign::Float64 = 1.0,
) where {D <: IS.FlattenIteratorWrapper{<:PSY.Device}}
    resolution = model_resolution(psi_container)
    dt = Dates.value(Dates.Minute(resolution)) / 60
    for d in devices
        cost_component = getfield(PSY.get_op_cost(d), cost_symbol)
        cost_expression =
            ps_cost(psi_container, var_name, PSY.get_name(d), cost_component, dt, sign)
        T_ce = typeof(cost_expression)
        T_cf = typeof(psi_container.cost_function)
        if T_cf <: JuMP.GenericAffExpr && T_ce <: JuMP.GenericQuadExpr
            psi_container.cost_function += cost_expression
        else
            JuMP.add_to_expression!(psi_container.cost_function, cost_expression)
        end
    end

    return
end


text:
# # Opinion spread

# ![](opinion.gif)

# This is a simple model of how an opinion spreads through a community.
# Each individual has a number of opinions as a list of integers.
# They can change their opinion by changing the numbers in the list.

# Agents can change their opinion at each step.
# They choose one of their neighbors randomly, and adopt one of the neighbor's opinion.
# They are more likely to adopt their neighbors opinion if the share more opinions with each other.

using Agents
using AgentsPlots
using Random

# ## Building the model
# ### 1. Model creation

mutable struct Citizen <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int}
    stabilized::Bool
    opinion::Array{Int,1}
    prev_opinion::Array{Int,1}
end

function create_model(; dims = (10, 10), nopinions = 3, levels_per_opinion = 4)
    space = GridSpace(dims)
    properties = Dict(:nopinions => nopinions)
    model = AgentBasedModel(
        Citizen,
        space,
        scheduler = random_activation,
        properties = properties,
    )
    for pos in positions(model)
        add_agent!(
            pos,
            model,
            false,
            rand(1:levels_per_opinion, nopinions),
            rand(1:levels_per_opinion, nopinions),
        )
    end
    return model
end

# ### 2. Stepping functions

function adopt!(agent, model)
    neighbor = rand(collect(nearby_ids(agent, model)))
    matches = model[neighbor].opinion .== agent.opinion
    nmatches = count(matches)
    if nmatches < model.nopinions && rand() < nmatches / model.nopinions
        switchId = rand(findall(x -> x == false, matches))
        agent.opinion[switchId] = model[neighbor].opinion[switchId]
    end
end

function update_prev_opinion!(agent, model)
    for i in 1:(model.nopinions)
        agent.prev_opinion[i] = agent.opinion[i]
    end
end

function is_stabilized!(agent, model)
    if agent.prev_opinion == agent.opinion
        agent.stabilized = true
    else
        agent.stabilized = false
    end
end

function agent_step!(agent, model)
    update_prev_opinion!(agent, model)
    adopt!(agent, model)
    is_stabilized!(agent, model)
end

# ## Running the model

levels_per_opinion = 4
model = create_model(nopinions = 3, levels_per_opinion = levels_per_opinion)

"run the model until all agents stabilize"
rununtil(model, s) = count(a->a.stabilized, allagents(model)) == length(positions(model))

agentdata, _ = run!(model, agent_step!, dummystep, rununtil, adata = [(:stabilized, count)])

# ## Plotting

# The plot shows the number of stable agents, that is, number of agents whose opinions
# don't change from one step to the next. Note that the number of stable agents can
# fluctuate before the final convergence.

plot(
    1:size(agentdata, 1),
    agentdata.count_stabilized,
    legend = false,
    xlabel = "generation",
    ylabel = "# of stabilized agents",
)

# ### Animation

# Here is an animation that shows change of agent opinions over time.
# The first three opinions of an agent determines its color in RGB.
levels_per_opinion = 3
ac(agent) = RGB((agent.opinion[1:3] ./ levels_per_opinion)...)
model = create_model(nopinions = 3, levels_per_opinion = levels_per_opinion)
anim = @animate for sp in 1:500
    step!(model, agent_step!)
    p = plotabm(model, ac = ac, as = 12, am = :square)
    title!(p, "Step $(sp)")
    if rununtil(model, 1)
        break
    end
end
gif(anim, "opinion.gif")


text:
using DataStructures
using LightGraphs

# For plotting
using GraphRecipes
import Plots
Plots.pyplot()

function find_element_num(n::T) where {T}
    s = Stack{T}()
    push!(s, n)

    num = 1
    while(length(s) > 0)
        num += 1
        temp = pop!(s)
        for neighbor in temp.neighbors
            push!(s, neighbor)
        end
    end
    return num
end

# We'll talk about this later if we have time!
function GraphRecipes.graphplot(n::T) where {T}

    # Constructing a LightGraphs graph
    element_num = find_element_num(n)

    g = SimpleGraph(element_num)

    s = Stack{T}()
    push!(s, n)

    while(length(s) > 0)
        temp = pop!(s)

        for neighbor in temp.neighbors
            add_edge!(g,neighbor.ID, temp.ID)
            push!(s, neighbor)
        end
    end
    graphplot(g;curves=false, markersize=0.2, names=1:element_num)

end

function create_tree(num_row::Int64, num_neighbor::Int64; mode="Recursive")

    if mode == "recursive"
        return create_tree(num_row, num_neighbor, [1])
    elseif mode == "stack"
        return create_tree_stack(num_row, num_neighbor, [1])
    elseif mode == "queue"
        return create_tree_queue(num_row, num_neighbor)
    else
        println("error mode ", mode, " not found for tree construction!")
    end
end

function create_tree_queue(num_row::Int64, num_neighbor::Int64)

    # Creating tree
    tree = create_tree(num_row, num_neighbor, [1])

    # relabelling elements
    q = Queue{Node}()
    enqueue!(q, tree)

    ID = 1
    while(length(q) > 0)
        temp = dequeue!(q)
        temp.ID = ID
        ID += 1
        for neighbor in temp.neighbors
            enqueue!(q, neighbor)
        end
    end

    return tree

end

function create_tree_stack(num_row::Int64, num_neighbor::Int64, ID::Vector)
    final_tree = Node(ID[1])
    if (num_row == 0)
        return final_tree
    end

    #preallocating for stack ordering
    for i = 1:num_neighbor
        push!(final_tree.neighbors, Node(ID[1]))
    end

    for i = num_neighbor:-1:1
        ID[1] += 1
        neighbor = create_tree(num_row - 1, num_neighbor, ID)
        final_tree.neighbors[i] = neighbor
    end

    return final_tree
end

function create_tree(num_row::Int64, num_neighbor::Int64, ID::Vector)
    final_tree = Node(ID[1])
    if (num_row == 0)
        return final_tree
    end

    for i = 1:num_neighbor
        ID[1] += 1
        neighbor = create_tree(num_row - 1, num_neighbor, ID)
        push!(final_tree.neighbors, neighbor)
    end

    return final_tree
end


text:
struct BypassIndex{T}
    index::T
end

bypass(val::T) where {T} = BypassIndex{T}(val)

"""
    dim_to_indices(arr::IndexedDimsArray{T, N}, inds::Tuple{Vararg{Any, N}}) where {T, N}

Convert indexes into the IndexedDimsArray into indexes into the underlying array.

```jldoctest; setup = :(using IndexedDims)
julia> arr = IndexedDimsArray([4.0, 5.0], ["a", "b"]);

julia> IndexedDims.dim_to_indices(arr, (["b"],))
([2],)
```
"""
function dim_to_indices(arr::IndexedDimsArray{T, N}, inds::Tuple{Vararg{Any, N}}) where {T, N}
    return map(unwrap_index, arr.indexes, inds)
end

function unwrap_index(acc::AcceleratedVector{T, <:AbstractVector{T}, <:AcceleratedArrays.AbstractUniqueIndex}, ind) where T
    return _unwrap_unique_index(acc, ind)
end
_unwrap_unique_index(acc::AcceleratedVector, ind) = findfirst(isequal(ind), acc)
_unwrap_unique_index(acc::AcceleratedVector, ind::Union{AbstractArray, Colon, BypassIndex}) = _unwrap_index(acc, ind)
unwrap_index(acc::AcceleratedVector, ind) = _unwrap_index(acc, ind)
# NOTE: we don't guarantee that we will find the first match, just any single match
_unwrap_index(acc::AcceleratedVector, ind) = findfirst(isequal(ind), acc)
_unwrap_index(acc::AcceleratedVector, ind::AbstractArray) = findall(in(ind), acc)
_unwrap_index(acc::AcceleratedVector, ::Colon) = Colon()
_unwrap_index(acc::AcceleratedVector, ind::BypassIndex) = ind.index

# --- #
# TODO: move to AcceleratedArrays.jl

function _unwrap_index(acc::AcceleratedVector{T, <:AbstractUnitRange{T}}, ind::AbstractRange{<:Integer}) where T<:Integer
    return findall(in(ind), parent(acc))
end

# --- #

function dim_getindex(arr::IndexedDimsArray{T, N}, inds::Vararg{Any, N}) where {T, N}
    unwrapped_indices = dim_to_indices(arr, inds)
    new_indices = (reindex(arr_idx, idx) for (arr_idx, idx) in zip(arr.indexes, unwrapped_indices) if !(idx isa Integer))
    new_underlying_indices = (Base.OneTo(length(ind)) for ind in new_indices)
    new_array = IndexedDimsArray(similar(parent(arr), new_underlying_indices...), new_indices...)
    _unsafe_getindex!(parent(new_array), parent(arr), unwrapped_indices...)
    return new_array
end

function dim_setindex!(arr::IndexedDimsArray{T, N}, val, inds::Vararg{Any, N}) where {T, N}
    return Base.setindex!(parent(arr), val, dim_to_indices(arr, inds)...)
end

function Base.getindex(arr::IndexedDimsArray{T, N}, inds::Vararg{Any, N}) where {T, N}
    return dim_getindex(arr, inds...)
end

function Base.setindex!(arr::IndexedDimsArray{T, N}, val, inds::Vararg{Any, N}) where {T, N}
    return dim_setindex!(arr, inds...)
end


text:
#=
This file is auto-generated. Do not edit.
=#

mutable struct RenewableFix <: RenewableGen
    name::String
    available::Bool
    bus::Bus
    activepower::Float64
    reactivepower::Float64
    tech::TechRenewable
    _forecasts::InfrastructureSystems.Forecasts
    internal::InfrastructureSystemsInternal
end

function RenewableFix(name, available, bus, activepower, reactivepower, tech, _forecasts=InfrastructureSystems.Forecasts(), )
    RenewableFix(name, available, bus, activepower, reactivepower, tech, _forecasts, InfrastructureSystemsInternal())
end

function RenewableFix(; name, available, bus, activepower, reactivepower, tech, _forecasts=InfrastructureSystems.Forecasts(), )
    RenewableFix(name, available, bus, activepower, reactivepower, tech, _forecasts, )
end

# Constructor for demo purposes; non-functional.

function RenewableFix(::Nothing)
    RenewableFix(;
        name="init",
        available=false,
        bus=Bus(nothing),
        activepower=0.0,
        reactivepower=0.0,
        tech=TechRenewable(nothing),
        _forecasts=InfrastructureSystems.Forecasts(),
    )
end

"""Get RenewableFix name."""
get_name(value::RenewableFix) = value.name
"""Get RenewableFix available."""
get_available(value::RenewableFix) = value.available
"""Get RenewableFix bus."""
get_bus(value::RenewableFix) = value.bus
"""Get RenewableFix activepower."""
get_activepower(value::RenewableFix) = value.activepower
"""Get RenewableFix reactivepower."""
get_reactivepower(value::RenewableFix) = value.reactivepower
"""Get RenewableFix tech."""
get_tech(value::RenewableFix) = value.tech
"""Get RenewableFix _forecasts."""
get__forecasts(value::RenewableFix) = value._forecasts
"""Get RenewableFix internal."""
get_internal(value::RenewableFix) = value.internal


text:
"
Interface to PyStan to run Stan scripts
"
module PyStan

using MD5
using PyCall

export StanModel_cache, extract_samples

const pickle = PyNULL()
const pystan = PyNULL()
function __init__()
  copy!(pickle,pyimport_conda("pickle","pickle"))
  copy!(pystan,pyimport_conda("pystan","pystan"))
end

basedir = abspath(joinpath(@__DIR__,".."))

standir = joinpath(basedir,"stan")

"""StanModel_cache(model_name,model_code=nothing,cache_dir=joinpath(standir,"cache"),kwargs...) -> StanModel

Use just as you would `StanModel`. `kwargs` are optional keyword arguments for `StanModel`.

Returns a Python StanModel object.
"""
function StanModel_cache(model_name,model_code=nothing,cache_dir=joinpath(standir,"cache"),kwargs...)
  if isnothing(model_code)
    stanfile = joinpath(standir,"$model_name.stan")
    model_code = open(stanfile) do file
      read(file, String)
    end
  end
  code_hash = bytes2hex(md5(model_code))
  cache_fn = "cached-$model_name-$code_hash.pkl"
  if !isnothing(cache_dir)
    if !isdir(cache_dir) mkdir(cache_dir) end
    cache_fn = joinpath(cache_dir,cache_fn)
  end
  if isfile(cache_fn)
    sm = pickle.load(pybuiltin("open")(cache_fn,"rb"))
  else
    sm = pystan.StanModel(model_code=model_code,model_name=model_name,kwargs...)
    @pywith pybuiltin("open")(cache_fn,"wb") as f begin
      pickle.dump(sm,f)
  end
  end
  sm
end

# Although the documentation claims it is possible, I have had
# no luck in extracting a Dict in PyStan when permuted=false, hence this function
"""extract_dict(f,inc_warmup=false) -> Dict

Extract samples from a PyStan fit object. Also returns a dict if permuted is false
"""
function extract_samples(f;kwargs...)
  s = f.extract(;kwargs...)
  if !isa(s,Dict)
    n = [f.flatnames;"lp__"]
    s = Dict(k=>s[:,:,i] for (i,k) in enumerate(n))
  end
  s
end
end # module

using .PyStan

"""
$SIGNATURES

Input data for the stan script "sde.stan"
"""
function standata(problem::ProblemLike{S,M}) where {S,M}
  @unpack t,y,p = problem
  eststate = S isa Val{:est}
  S isa Val{:fixed} && @warn "SDE estimation with :fixed state not yet implemented"
  method = if M isa Val{:elimQ}
    error("Not implemented yet")
  elseif M isa Val{:matfrac}
    1
  end
  D = Dict(
    "N"=>length(t),
    "p"=>p,
    "t"=>t,
    "y"=>y,
    "eststate"=>Int64(eststate),
    "ref_prior_sigma_eps"=>Int64(false),
    "Tp"=>0,
    "method"=>method
  )
  if p>1
    D["p"] = p
  end
  return D
end

function standata(problem::ProblemProb)
  @unpack ref_prior_σϵ,Tp = problem
  D = standata(likeproblem(problem))
  D["ref_prior_sigma_eps"] = Int64(ref_prior_σϵ)
  D["Tp"] = isnothing(Tp) ? 0 : Tp
  return D
end

"""
$SIGNATURES
Stan parameters for process `proc`. `eststate` indicates whether
the state will be estimated.
"""
function stanpar(proc,emode,method)
  eq = proc.eq;
  spar =  Dict(
    "a" => eq.a,
    "sigma_eps" => eq.σϵ,
    "estinit" => Int64(emode isa Val{:est})
  )
  if emode isa Val{:est}
    proc.state isa State || error("State must be fixed Only fixed state (zero P) or no state is supported")
    spar["z2"] = proc.state.z[2:end]
  end
  return spar
end

function _modelsamples_stan_permuted(samples)
  S = length(samples["sigma_eps"])
  p = size(samples["a"],2)
  estinit = haskey(samples,"z")
  procsam = map(1:S) do i
    s = Dict(k=>v[i,:] for (k,v) in samples)
    eq = Sde(s["a"],s["sigma_eps"][1])
    state = estinit ? State(s["z"]) : nothing
    SdeProcess(eq,state)
  end
  return procsam
end

function _modelsamples_stan(samples)
  p = count(x->startswith(x,"a"),keys(samples))
  S = length(samples["sigma_eps"])
  procsam = map(1:S) do i
    s = Dict(k=>v[i] for (k,v) in samples)
    a = map(pc->s["a[$pc]"],1:p)
    eq = Sde(a,s["sigma_eps"])
    eststate = any(startswith.(keys(samples),"z")) # Is this correct?
    state = if eststate
      z = map(pc->s["a[$pc]"],1:p)
      State(z)
    else
      statstate(eq)
    end
    SdeProcess(eq,state)
  end
  return procsam
end

"""
$SIGNATURES

Converts `samples` returned by PyStan to an array of SdeProcess.
`permuted` is the parameter used with StanModel.extract (default `true`)
to extract the samples.
"""
modelsamples_stan(samples,permuted=true) = if permuted
    _modelsamples_stan_permuted(samples)
  else
    _modelsamples_stan(samples)
end

function process_from_stan(problem::Problem{S,M},spar) where {S,M}
  a = spar["a"][:] # [:] to convert 0-dim array that is returned for SDE(1)
  σϵ = spar["sigma_eps"][1]
  eq = Sde(a,σϵ)
  state = if S isa Val{:stat}
    statstate(eq)
  else
    state = S isa Val{:fixed} ? problem.state : State(spar["z"][:])
  end
  proc = SdeProcess(eq,state)
  return proc
end

function ml_estimate_init_stan(tor,problem::Problem{S,M},procstart) where {S,M}
  @info "Stan ML estimate is incorrect because of determinant adjustment that should not be used for ML" maxlog=1
  sm = StanModel_cache("sde")
  data = standata(problem)
  sparstart = stanpar(procstart,S,M)
  optim = sm.optimizing(data,init=sparstart,iter=tor.iter,as_vector=false)

  spar = optim["par"]
  est = process_from_stan(problem,spar)
  return est,optim
end

function posterior_samples_stan(tor,problem::ProblemProb{S,M},procstarts) where {S,M}
  sm = StanModel_cache("sde")
  Tav = mean(diff(problem.t))
  N = length(problem.y)
  Tp = isnothing(problem.Tp) ? nothing : Tav
  data = standata(problem)
  cstate = compatiblestate(problem.p,problem.y)
  eststate = S isa Val{:est}
  parinit = stanpar.(procstarts,eststate)
  seed = rand(tor.rng,UInt16)
  tor.chains>1 && error("tor.chains > 1 with Stan results in errors; currently not supported")
  fit = sm.sampling(data,seed=seed,init=parinit,iter=tor.iter,chains=tor.chains)
  info = (fit=fit,)
  samples = extract_samples(fit)
  procsam = modelsamples_stan(samples)
  return procsam,info
end


text:
import Gadfly

plot2!(::Val{:Gadfly}, p::Union{Plot,Nothing}, X, Ys; kind, title, xlab, ylab, legend, legendpos, names, colors, xlim, ylim, ycat, xunit, yunit, aspect) = begin
    n = length(Ys)
    Xs = [X for _ in 1:n]
    kinds = [kind for _ in 1:n]

    if kind == :line
        geoms = [Gadfly.Geom.line]
    elseif kind == :scatter
        geoms = [Gadfly.Geom.point]
    elseif kind == :scatterline
        geoms = [Gadfly.Geom.point, Gadfly.Geom.line]
    elseif kind == :step
        geoms = [Gadfly.Geom.step]
    else
        error("unrecognized plot kind = $kind")
    end

    #HACK: manual_color_key() expects [] while colorkey() expects nothing
    keypos = isnothing(legendpos) ? [] : legendpos .* [Gadfly.w, Gadfly.h]

    theme = Gadfly.Theme(
        background_color="white",
        plot_padding=[5*Gadfly.mm, 5*Gadfly.mm, 5*Gadfly.mm, 0*Gadfly.mm],
        major_label_font_size=10*Gadfly.pt,
        key_title_font_size=9*Gadfly.pt,
        key_position=isempty(keypos) ? :right : :inside,
        point_size=0.7*Gadfly.mm,
        discrete_highlight_color=_->Gadfly.RGBA(1, 1, 1, 0),
    )

    create_colors(colors; n0=0) = begin
        n = length(colors)
        C = Gadfly.Scale.default_discrete_colors(n0+n)[n0+begin:end]
        f(c::Int, _) = p.opt[:colors][c]
        f(c, _) = parse(Gadfly.Colorant, c)
        f(::Nothing, i) = C[i]
        [f(c, i) for (i, c) in enumerate(colors)]
    end
    colorkey(colors) = begin
        NC = filter(x -> let (n, c) = x; !isempty(n) end, collect(zip(names, colors)))
        if !isempty(NC)
            N, C = first.(NC), last.(NC)
            Gadfly.Guide.manual_color_key(legend, N, C; pos=keypos)
        end
    end
    colorkey!(key, colors) = begin
        k = colorkey(colors)
        if !isnothing(k)
            append!(key.labels, k.labels)
            append!(key.colors, k.colors)
        end
        k
    end
    update_color!(guides, colors) = begin
        #TODO: very hacky approach to append new plots... definitely need a better way
        keys = filter(x -> x isa Gadfly.Guide.ManualDiscreteKey, guides)
        if isempty(keys)
            key = colorkey(colors)
            !isnothing(key) && push!(guides, key)
        else
            key = only(keys)
            colorkey!(key, colors)
        end
    end
    create_layers(colors) = [Gadfly.layer(x=Xs[i], y=Ys[i], geoms..., Gadfly.Theme(theme; default_color=colors[i])) for i in 1:n]

    if isnothing(p)
        xmin, xmax = xlim
        ymin, ymax = ylim

        scales = if kind == :step
            [
                Gadfly.Scale.y_discrete(levels=ycat),
                Gadfly.Coord.cartesian(; xmin, xmax, ymin=1, ymax=length(ycat), aspect_ratio=aspect),
            ]
        #HACK: aesthetic adjustment for boolean (flag) plots
        #TODO: remove special adjustment in favor of new step plot
        elseif eltype(ylim) == Bool
            [
                #HACK: ensure correct level order (false low, true high)
                Gadfly.Scale.y_discrete(levels=[false, true]),
                #HACK: shift ylim to avoid clipping true values (discrete false=1, true=2)
                Gadfly.Coord.cartesian(; xmin, xmax, ymin=1, ymax=2, aspect_ratio=aspect),
            ]
        else
            [
                Gadfly.Coord.cartesian(; xmin, xmax, ymin, ymax, aspect_ratio=aspect),
            ]
        end
        guides = [
            Gadfly.Guide.title(title),
            Gadfly.Guide.xlabel(xlab),
            Gadfly.Guide.ylabel(ylab),
        ]
        colors = create_colors(colors)
        update_color!(guides, colors)
        layers = create_layers(colors)
        obj = Gadfly.plot(
            scales...,
            guides...,
            layers...,
            theme,
        )
        p = Plot(obj; Xs, Ys, kinds, colors, title, xlab, ylab, legend, names, xlim, ylim, xunit, yunit, aspect)
    else
        obj = p.obj
        n0 = length(obj.layers)
        colors = create_colors(colors; n0)
        update_color!(obj.guides, colors)
        foreach(l -> Gadfly.push!(obj, l), create_layers(colors))
        for l in create_layers(colors)
            Gadfly.push!(obj, l)
        end
        update!(p; Xs, Ys, kinds, colors, names)
    end
    p
end

plot3!(::Val{:Gadfly}, X, Y, Z; kind, title, legend, legendpos, xlab, ylab, zlab, xlim, ylim, zlim, zgap, zlabgap, aspect) = begin
    if kind == :heatmap
        geom = Gadfly.Geom.rectbin
        data = (x=X, y=Y, color=Z)
    elseif kind == :contour
        levels = isnothing(zgap) ? 100 : collect(zlim[1]:zgap:zlim[2])
        geom = Gadfly.Geom.contour(; levels)
        data = (x=X, y=Y, z=Z)
    else
        error("unrecognized plot kind = $kind")
    end

    #HACK: colorkey() expects nothing while manual_color_key() expects []
    keypos = isnothing(legendpos) ? nothing : legendpos .* [Gadfly.w, Gadfly.h]

    theme = Gadfly.Theme(
        background_color="white",
        plot_padding=[5*Gadfly.mm, 5*Gadfly.mm, 5*Gadfly.mm, 0*Gadfly.mm],
        major_label_font_size=10*Gadfly.pt,
        key_title_font_size=9*Gadfly.pt,
        key_position=legend ? isnothing(keypos) ? :right : :inside : :none,
    )

    label(z) = begin
        #TODO: remove redundant creation of Scale
        zmin, zmax = zlim
        zspan = zmax - zmin
        scale = Gadfly.Scale.color_continuous(minvalue=zmin, maxvalue=zmax)
        color = scale.f((z - zmin) / zspan)

        i = findmin(abs.(Z .- z))[2]
        #HACK: ignore lables presumably out of bound
        if i == firstindex(Z) || i == lastindex(Z)
            Gadfly.Guide.annotation(Gadfly.compose(Gadfly.context()))
        else
            x, y = X[i], Y[i]
            Gadfly.Guide.annotation(
                Gadfly.compose(
                    Gadfly.context(),
                    Gadfly.Compose.text(x, y, string(z), Gadfly.hcenter, Gadfly.vcenter),
                    Gadfly.font(theme.minor_label_font),
                    Gadfly.fontsize(theme.minor_label_font_size),
                    Gadfly.fill(color),
                    Gadfly.stroke("white"),
                    Gadfly.linewidth(0.1*Gadfly.pt),
                )
            )
        end
    end
    labels = isnothing(zlabgap) ? () : label.(zlim[1]:zlabgap:zlim[2])

    scales = [
        #TODO: fix performance regression with custom system image
        #Gadfly.Scale.x_continuous,
        #Gadfly.Scale.y_continuous,
        Gadfly.Scale.color_continuous(minvalue=zlim[1], maxvalue=zlim[2]),
        Gadfly.Coord.cartesian(xmin=xlim[1], ymin=ylim[1], xmax=xlim[2], ymax=ylim[2], aspect_ratio=aspect),
    ]
    guides = [
        Gadfly.Guide.title(title),
        Gadfly.Guide.xlabel(xlab),
        Gadfly.Guide.ylabel(ylab),
        Gadfly.Guide.colorkey(title=zlab; pos=keypos),
    ]

    obj = Gadfly.plot(
        scales...,
        guides...,
        geom,
        labels...,
        theme;
        data...,
    )
    Plot(obj; X, Y, Z, kind, title, xlab, ylab, zlab, xlim, ylim, zlim, aspect)
end

#HACK: use non-interactive SVG instead of SVGJS
_show(io::IO, m::MIME"text/html", p::Gadfly.Plot) = begin
    w = Gadfly.Compose.default_graphic_width
    h = Gadfly.Compose.default_graphic_height
    Gadfly.SVG(io, w, h, false)(p)
end

import Cairo
import ImageMagick
import FileIO
# https://github.com/tshort/SixelTerm.jl
sixel(p::Gadfly.Plot) = begin
    png = IOBuffer()
    #HACK: needs to set emit_on_finish false
    w = Gadfly.Compose.default_graphic_width
    h = Gadfly.Compose.default_graphic_height
    p |> Gadfly.PNG(png, w, h, false; dpi=144)
    im = ImageMagick.load(png)
    six = IOBuffer()
    st = FileIO.Stream(FileIO.format"six", six)
    ImageMagick.save(st, im)
    write(stdout, take!(six))
    nothing
end


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

"""
    add_fano_noise(E_dep::RealQuantity, E_ionisation::RealQuantity, f_fano::Real)::RealQuantity

Adds Fano noise to an energy deposition `E_dep`, assuming a detector material
ionisation energy `E_ionisation` and a Fano factor `f_fano`.

## Arguments
* `E_dep::RealQuantity`: Energy deposited in a semiconductor material.
* `E_ionisation`: Energy needed to create one electron-hole-pair in the semiconductor material.
* `f_fano`: Fano factor of the material.

## Example
```julia
add_fano_noise(100u"keV", 2.95u"eV", 0.129)
```

Some material properties are stored in `SolidStateDetectors.material_properties` and can be used here:
```julia
material = SolidStateDetectors.material_properties[:HPGe]
add_fano_noise(100u"keV", material.E_ionisation, material.f_fano)
```

!!! note
    Using values with units for `E_dep` or `E_ionisation` requires the package [Unitful.jl](https://github.com/PainterQubits/Unitful.jl).
"""
function add_fano_noise end
export add_fano_noise

function add_fano_noise(rng::AbstractRNG, E_dep::RealQuantity, E_ionisation::RealQuantity, f_fano::Real)
    target_unit = unit(E_dep)
    n_expected = uconvert(Unitful.NoUnits, E_dep/E_ionisation)

    # dist = Distributions.Poisson(n_expected)
    # n_observed = n_expected + f_fano * (rand(dist) - n_expected)

    # Use normal distribution instead of Poisson, makes no differenc for
    # typical numbers of charges:
    σ = sqrt(n_expected)
    n_observed = n_expected + f_fano * σ * randn(rng)
    n_observed_pos = max(0, n_observed)

    uconvert(target_unit, E_ionisation * n_observed_pos)
end

function add_fano_noise(rng::AbstractRNG, evts::DetectorHitEvents, E_ionisation::RealQuantity, f_fano::Real)
    noisy_edep = deepmap(x -> add_fano_noise(x, E_ionisation, f_fano), evts.edep)

    TypedTables.Table(merge(
        Tables.columns(evts),
        (edep = noisy_edep,)
    ))
end

add_fano_noise(x, E_ionisation::RealQuantity, f_fano::Real) =
    add_fano_noise(Random.GLOBAL_RNG, x, E_ionisation, f_fano)


text:
function calculateLogProb(x,p::PrefEl)
	out = 0.0

	# println("x = $x")

	for i in 1:size(p.strict,1)
		out += getLogStrictProb(x,p.data[p.strict[i,1],:],p.data[p.strict[i,2],:],p.Sigma)
	end

	for i in 1:size(p.indif,1)
		out += getLogIndifProb(x,p.data[p.indif[i,1],:],p.data[p.indif[i,2],:],p.Sigma)
	end

	out += logPrior(x,p.priors.dists)

	println(out)

	return out
end

function getLogStrictProb{R<:Real}(x::Array{R,1}, a::Array{R,2}, b::Array{R,2}, Sigma::Array{R,2})
	d = a .- b
	varAlongD = d * Sigma * d'
	meanAlongD = dot(vec(d),vec(x))
	temp = Normal(0,sqrt(varAlongD[1]))
	return logcdf(temp,meanAlongD)
end

function getLogIndifProb{R<:Real}(x::Array{R,1}, a::Array{R,2}, b::Array{R,2}, Sigma::Array{R,2})
	d = a .- b
	varAlongD = (d * Sigma * d')[1]
	meanAlongD = dot(vec(d),x)
	temp = Normal(0,sqrt(varAlongD))
	return log(cdf(temp,meanAlongD + 0.5) - cdf(temp,meanAlongD - 0.5))
end

function logPrior{R <: Real}(x::Vector{R}, priors::Vector{Distribution})
	out = 0.0
	@simd for i in 1:length(x)
		if typeof(priors[i]) != Uniform || priors[i].a != -Inf
			@inbounds out += logpdf(priors[i],x[i])
		end
	end
	if out == -Inf
		out = realmin()
	end
	return out
end

using NLopt
function infer(p::PrefEl; method = "MAP")

	if method == "MAP"

		# Set up optimization
		n = size(p.data,2)
		opt = Opt(:LN_NEWUOA_BOUND,n) # slow
		# opt = Opt(:LN_BOBYQA, n) # really slow
		# opt = Opt(:LN_COBYLA, n) # slow
		# opt = Opt(:LD_LBFGS,n)
		xtol_rel!(opt,1e-4)
		# ftol_rel!(opt,1e-8)

		# Create dummy function
		f(x,grad) = calculateLogProb(x,p)

		# Want MAP estimate
		max_objective!(opt,f)

		# Restrict optimization to be within support of prior
		# necessary for priors like the Exponential dist
		lb = zeros(n)
		ub = zeros(n)
		for i in 1:n
			s = support(p.priors.dists[i])
			lb[i] = max(s.lb,p.priors.lb[i])
			ub[i] = min(s.ub,p.priors.ub[i])
		end
		lower_bounds!(opt,lb)
		upper_bounds!(opt,ub)

		# Set the starting point to be the prior mode
		means = zeros(F,n)
		for i in 1:n
			means[i] = mean(p.priors.dists[i]) # Would like to use mode, but exponential will kill autodiff because it'll check negative values
			if isnan(means[i])
				means[i] = 0.0
			end
		end

		minF, maxX, ret = optimize(opt,means)

		println(size(maxX))

		println("Optimization returned with code $ret")

		return maxX

	elseif method == "MCMC"
		samples = getMCMCSamples(p)
		return vec(mean(samples,1))
	else
		error("Unknown method $method")
	end
end

function suggest(p::PrefEl)
	n = size(p.data,1)

	# Set up all possible recommendations
	recs = cell(div(n * (n-1),2)) # to maintain int-ness
	count = 1
	for i in 1:n
		for j in i+1:n
			recs[count] = (i,j)
			count += 1
		end
	end

	# Create dummy function for pmap
	f(x) = estimateEntropy(p,x[1],x[2])

 	# Map f to our candidates
	results = pmap(f,recs)

	# Find the optimum
	_, id = findmin(results)

	return recs[id]
end

using KernelDensity
function estimateEntropy(p::PrefEl, a::Int,b::Int)
	pNew = deepcopy(p) # need pass-by-value, or else procs will talk over each other and die

	# Get the probability that a > b
	est = vec(infer(p))
	abProb = exp(getLogStrictProb(est,p.data[a,:],p.data[b,:],p.Sigma))

	# Get samples from posterior if a > b
	pNew.strict = [pNew.strict; a b]
	abSamples = getMCMCSamples(pNew)

	# Get samples from posterior if b > a
	pNew.strict = [pNew.strict[1:end-1,:]; b a]
	baSamples = getMCMCSamples(pNew)

	# Calculate entropies for each solution
	aEntropy = 0.0
	bEntropy = 0.0
	for dim in 1:size(abSamples,2)
		aKernel = KernelDensity.kde(abSamples[:,dim])
		bKernel = KernelDensity.kde(baSamples[:,dim])
		aEntropy += dot(aKernel.density[:], log(aKernel.density[:]))
		bEntropy += dot(bKernel.density[:], log(bKernel.density[:]))
	end

	# Return expected entropy
	return -(abProb * aEntropy + (1-abProb) * bEntropy)
end

function getMCMCSamples(p::PrefEl, samples = 1000)
	n = size(p.data,2)

	sigma = 1.0
	C = eye(n) * sigma
	propDist = MvNormal(C)

	samps = zeros(samples,n)

	samps[1,:] = (infer(p))'
	zOld = calculateLogProb(vec(samps[1,:]),p)

	for i in 2:samples
		prop = samps[i-1,:] .+ rand(propDist)'

		zNew = calculateLogProb(vec(samps[1,:]),p)

		if zNew > zOld || log(rand()) < zNew - zOld
			zOld = zNew
			samps[i,:] = prop
		else
			samps[i,:] = samps[i-1,:]
		end
	end

	return samps
end

max_stars_repo_name: jlepird/prefEl

max_issues_repo_name: jlepird/prefEl

max_forks_repo_name: jlepird/prefEl


text:
using Test
using Rectangle

import Rectangle: _inorder, isnil

const RUN_BENCHMARK = false

function tree_get_data(t::AbstractBST{K, V}) where {K, V}
    karr = Vector{K}()
    varr = Vector{V}()

    for p in collect(Iterator(t))
        push!(karr, p[1])
        push!(varr, p[2])
    end
    return karr, varr
end

function bstvalidity(t::AbstractBST)
    valid = true
    _inorder(t.root, t) do n
        !valid && return valid
        r1 = (!isnil(t, n.l) && !(n < n.l)) || isnil(t, n.l)
        r2 = (!isnil(t, n.r) && !(n.r < n)) || isnil(t, n.r)
        valid &= (r1 & r2)
        return valid
    end
    return valid
end

function countvalidity(t::RandomizedBinarySearchTree{K,V}) where {K,V}
    valid = true
    _inorder(t.root, t) do n
        !valid && return valid
        valid &= n.l.n + n.r.n + 1 == n.n
    end
end

function parentvalidity(t::AbstractBST)
    valid = true
    _inorder(t.root, t) do n
        @assert !isnil(t, n)
        !valid && return valid
        r1 = ( isnil(t, n.p) && (t.root === n))
        r2 = (!isnil(t, n.p) && (n.p.l === n || n.p.r === n))
        @assert isnil(t, n.p) || (n.p.l !== n.p.r)
        valid &= (r1 | r2)
        return valid
    end
    return valid
end

#=
CLRS 3rd Ed. Chapter 13

1. Every node is either red or black.
2. The root is black.
3. Every leaf (NIL ) is black.
4. If a node is red, then both its children are black.
5. For each node, all simple paths from the node to descendant leaves contain the
same number of black nodes.

=#
function rbvalidity(t::RBTree{K, V}) where {K, V}
    valid = (t.root.red == false)
    !valid && error("Invalid value at root: $(t.root.k)")
    valid = (valid && (t.nil.red == false))
    !valid && error("Invalid value at nil")
    _inorder(t.root, t) do n
        valid = (valid && ((n.red && !n.l.red && !n.r.red && !n.p.red) || !n.red))
        !valid && println("Invalid value at $(n.k)")
        cl = cr = 0
        tn = n
        while !isnil(t, tn)
            cl += (!tn.red ? 1 : 0)
            tn = tn.l
        end
        tn = n
        while !isnil(t, tn)
            cr += (!tn.red ? 1 : 0)
            tn = tn.r
        end
        valid = (valid && cl == cr)
        !valid && println("Invalid balance value at $(n.k): $(cl) and $(cr)")
        return valid
    end
    return valid
end

@testset "Binary Search Trees" begin
    a = [4, 5, 3, 1, 10, 2, 7, 6, 8, 9]
    t = BinarySearchTree{Int, Int}()

    @test get!(t, 1, 2) == get(t, 1, 3)
    @test endswith(string(t), "Tree with 1 nodes.\nRoot at: 1.\n")
    @test get!(t, 2, 3) == get(t, 2, 4)
    @test get!(t, 2, 4) == 3
    @test get(t, 2, 4) == 3
    @test delete!(t, 1) == (1 => 2)
    @test delete!(t, 2) == (2 => 3)
    @test isempty(t)

    for i in a
        insert!(t, i, i*10)
    end
    @test collect(Iterator(t)) == [1=>10, 2=>20, 3=>30, 4=>40, 5=>50,
                                   6=>60, 7=>70, 8=>80, 9=>90, 10=>100]

    b = [i for i in 1:100]
    for i = 1:50
        n1 = rand(b)
        n2 = rand(b)
        b[n1], b[n2] = b[n2], b[n1]
    end
    s = BinarySearchTree{Int, Int}()
    for i in b
       insert!(s, i, 10i)
    end
    bb = [i for i in 1:10]
    @test collect(Iterator(s, -3, 10)) == Pair.(bb, 10bb)
    bb = [i for i in 11:20]
    @test collect(Iterator(s, bb[1], bb[end])) == Pair.(bb, 10bb)
    bb = [i for i in 101:103]
    @test collect(Iterator(s, bb[1], bb[end])) == []
    bb = [i for i in 99:103]
    @test collect(Iterator(s, bb[1], bb[end])) == [99 => 990, 100 => 1000]
    bb = [i for i in -3:2]
    @test collect(Iterator(s, bb[1], bb[end])) == [1 => 10, 2 => 20]
    bb = [i for i in -3:0]
    @test collect(Iterator(s, bb[1], bb[end])) == []
    for i = 1:10
        insert!(s, i, 10i)
    end
    @test collect(Iterator(s, 1, 2)) == [1 => 10, 1=>10, 2=>20, 2 => 20]
    @test length(t) == 10
    @test !isempty(t)
    @test maximum(t) == (10 =>100)
    @test minimum(t) == (1 => 10)
    data = tree_get_data(t)
    @test data[1] == sort(a)
    @test delete!(t, 5) == (5 => 50)
    @test delete!(t, 25) === nothing
    data = tree_get_data(t)
    @test data[2] == [10, 20, 30, 40, 60, 70, 80, 90, 100]
    empty!(t)
    @test isempty(t)
    @testset "Non-unique" begin
        t.unique = false
        n = 1000
        a = []
        for i = 1:n
            j = rand(1:10000)
            push!(a, j, i)
        end
        for i = 1:n
            insert!(t, a[i], i)
        end
        @test parentvalidity(t)
        @test bstvalidity(t)
        for i = 1:div(n, 2)
            delete!(t, a[i])
        end
        @test parentvalidity(t)
        @test bstvalidity(t)
        for i = 1:div(n, 2)
            insert!(t, a[i], i)
        end
        @test parentvalidity(t)
        @test bstvalidity(t)
    end
    @testset "Unique" begin
        t = BinarySearchTree{Int, Int}()
        t.unique = true
        n = 1000
        b = [i for i in 1:n]
        for i = 1:div(n, 2)
            x = rand(1:n)
            y = rand(1:n)
            b[x], b[y] = b[y], b[x]
        end
        a = Tuple{Int, Int}[]
        for i = 1:n
            insert!(t, b[i], i)
            push!(a, (b[i], i))
        end
        ka, kv = tree_get_data(t)
        aa = sort(a, lt= (x, y) -> x[1] < y[1])
        @test aa == collect(zip(ka, kv))
        @test bstvalidity(t)
        @test parentvalidity(t)
        @test_throws ErrorException insert!(t, 10, 10)

        # Randomly delete 50% of the data. Data in a are not sorted

        n_2 = div(n, 2)
        for i = 1:n_2
            Base.delete!(t, b[i])
        end
        bb = b[n_2+1:end]
        a = a[n_2+1:end]
        ka, kv = tree_get_data(t)
        aa = sort(a, lt= (x, y) -> x[1] < y[1])
        @test aa == collect(zip(ka, kv))
        @test bstvalidity(t)
        @test parentvalidity(t)
    end
end

@testset "Red and Black Trees" begin

    a = [4, 5, 3, 1, 10, 2, 7, 6, 8, 9]
    t = RBTree{Int, Int}()

    @test isempty(t)

    for i in a
        insert!(t, i, i*10)
    end
    @test collect(Iterator(t)) == [1=>10, 2=>20, 3=>30, 4=>40, 5=>50,
                                   6=>60, 7=>70, 8=>80, 9=>90, 10=>100]
    @test collect(Iterator(t, 3, 8)) == [3=>30, 4=>40, 5=>50,
                                         6=>60, 7=>70, 8=>80]
    @test parentvalidity(t)
    @test bstvalidity(t)
    @test rbvalidity(t)
    @test length(t) == 10
    @test !isempty(t)
    @test maximum(t) == (10 =>100)
    @test minimum(t) == (1 => 10)
    data = tree_get_data(t)
    @test data[1] == sort(a)
    @test delete!(t, 5) == (5 => 50)
    @test parentvalidity(t)
    @test bstvalidity(t)
    @test rbvalidity(t)
    data = tree_get_data(t)
    @test data[2] == [10, 20, 30, 40, 60, 70, 80, 90, 100]
    empty!(t)
    @test isempty(t)

    @testset "Non-unique" begin
        t.unique = false
        n = 10000
        a = Int[]
        for i = 1:n
            j = rand(1:1000)
            push!(a, j)
        end
        for i = 1:n
            insert!(t, a[i], i)
        end
        @test parentvalidity(t)
        @test bstvalidity(t)
        @test rbvalidity(t)
        delete!(t, t.root.k)
        for i = 1:div(n, 2)
            delete!(t, a[i])
        end
        @test rbvalidity(t)
        @test parentvalidity(t)
        @test bstvalidity(t)

        for i = 1:div(n, 2)
            insert!(t, a[i], i)
        end
        @test rbvalidity(t)
        @test parentvalidity(t)
        @test bstvalidity(t)
    end
    empty!(t)
    @testset "Unique" begin
        t.unique = true
        n = 10000
        b = [i for i in 1:n]

        for i = div(n, 2):-1:1
            x = rand(1:n)
            y = rand(1:n)
            b[x], b[y] = b[y], b[x]
        end
        a = Tuple{Int, Int}[]
        pa = Pair{Int, Int}[]
        for i = 1:n
            insert!(t, b[i], i)
            push!(a, (b[i], i))
            push!(pa, (b[i] => i))
        end
        ka, kv = tree_get_data(t)
        aa = sort(a, lt= (x, y) -> x[1] < y[1])
        pa = sort(pa, lt= (x, y) -> x[1] < y[1])
        @test aa == collect(zip(ka, kv))
        @test pa == collect(Iterator(t))
        @test rbvalidity(t)
        @test bstvalidity(t)
        @test parentvalidity(t)
        #@test_throws ErrorException insert!(t, 10, 10)

        # Randomly delete 50% of the data. Data in a are not sorted
        n_2 = div(n, 2)
        for i = 1:n_2
            Base.delete!(t, b[i])
            ka, kv = tree_get_data(t)
        end
        bb = b[n_2+1:end]
        a  = a[n_2+1:end]
        ka, kv = tree_get_data(t)
        aa = sort(a, lt= (x, y) -> x[1] < y[1])
        @test aa == collect(zip(ka, kv))
        @test rbvalidity(t)
        @test bstvalidity(t)
        @test parentvalidity(t)
    end
end

@testset "Random BST" begin
    a = [4, 5, 3, 1, 10, 2, 7, 6, 8, 9]
    t = RandomizedBinarySearchTree{Int, Int}()

    @test get!(t, 1, 2) == get(t, 1, 3)
    @test endswith(string(t), "Tree with 1 nodes.\nRoot at: 1.\n")
    @test get!(t, 2, 3) == get(t, 2, 4)
    @test get!(t, 2, 4) == 3
    @test get(t, 2, 4) == 3
    @test delete!(t, 1) == (1 => 2)
    @test delete!(t, 2) == (2 => 3)
    @test isempty(t)
    for i in a
        insert!(t, i, i*10)
    end
    @test collect(Iterator(t)) == [1=>10, 2=>20, 3=>30, 4=>40, 5=>50,
                                   6=>60, 7=>70, 8=>80, 9=>90, 10=>100]
    b = [i for i in 1:100]
    for i = 1:50
        n1 = rand(b)
        n2 = rand(b)
        b[n1], b[n2] = b[n2], b[n1]
    end
    s = RandomizedBinarySearchTree{Int, Int}()
    for i in b
       insert!(s, i, 10i)
    end
    bb = [i for i in 1:10]
    @test collect(Iterator(s, -3, 10)) == Pair.(bb, 10bb)
    bb = [i for i in 11:20]
    @test collect(Iterator(s, bb[1], bb[end])) == Pair.(bb, 10bb)
    bb = [i for i in 101:103]
    @test collect(Iterator(s, bb[1], bb[end])) == []
    bb = [i for i in 99:103]
    @test collect(Iterator(s, bb[1], bb[end])) == [99 => 990, 100 => 1000]
    bb = [i for i in -3:2]
    @test collect(Iterator(s, bb[1], bb[end])) == [1 => 10, 2 => 20]
    bb = [i for i in -3:0]
    @test collect(Iterator(s, bb[1], bb[end])) == []
    for i = 1:10
        insert!(s, i, 10i)
    end
    @test collect(Iterator(s, 1, 2)) == [1 => 10, 1=>10, 2=>20, 2 => 20]
    @test length(t) == 10
    @test !isempty(t)
    @test maximum(t) == (10 =>100)
    @test minimum(t) == (1 => 10)
    data = tree_get_data(t)
    @test data[1] == sort(a)
    @test delete!(t, 5) == (5 => 50)
    @test delete!(t, 25) === nothing
    data = tree_get_data(t)
    @test data[2] == [10, 20, 30, 40, 60, 70, 80, 90, 100]
    empty!(t)
    @test isempty(t)
    @testset "Non-unique" begin
        t.unique = false
        n = 10000
        a = []
        for i = 1:n
            j = rand(1:10000)
            push!(a, j, i)
        end
        for i = 1:n
            insert!(t, a[i], i)
        end
        @test countvalidity(t)
        @test parentvalidity(t)
        @test bstvalidity(t)
        x = delete!(t, t.root.k)
        @assert t.root.n == 9999
        @test countvalidity(t)
        @test parentvalidity(t)
        @test bstvalidity(t)
        insert!(t, x.first, x.second)
        @test countvalidity(t)
        @test parentvalidity(t)
        @test bstvalidity(t)
        @assert t.root.n == 10000
        for i = 1:div(n, 2)
            delete!(t, a[i])
        end
        @test countvalidity(t)
        @test parentvalidity(t)
        @test bstvalidity(t)
        for i = 1:div(n, 2)
            insert!(t, a[i], i)
        end
        @test countvalidity(t)
        @test parentvalidity(t)
        @test bstvalidity(t)
    end
    @testset "Unique" begin
        t = RandomizedBinarySearchTree{Int, Int}()
        t.unique = true
        n = 10000
        b = [i for i in 1:n]
        for i = 1:div(n, 2)
            x = rand(1:n)
            y = rand(1:n)
            b[x], b[y] = b[y], b[x]
        end
        a = Tuple{Int, Int}[]
        for i = 1:n
            insert!(t, b[i], i)
            push!(a, (b[i], i))
        end
        ka, kv = tree_get_data(t)
        aa = sort(a, lt= (x, y) -> x[1] < y[1])
        @test aa == collect(zip(ka, kv))
        @test countvalidity(t)
        @test bstvalidity(t)
        @test parentvalidity(t)
        @test_throws ErrorException insert!(t, 10, 10)

        # Randomly delete 50% of the data. Data in a are not sorted

        n_2 = div(n, 2)
        for i = 1:n_2
            Base.delete!(t, b[i])
        end
        bb = b[n_2+1:end]
        a = a[n_2+1:end]
        ka, kv = tree_get_data(t)
        aa = sort(a, lt= (x, y) -> x[1] < y[1])
        @test aa == collect(zip(ka, kv))
        @test bstvalidity(t)
        @test parentvalidity(t)
    end
end

RUN_BENCHMARK && using BenchmarkTools

@static if RUN_BENCHMARK

    begin
        t = BinarySearchTree{Int, Int}()
        t.unique = false
        n = 1000
        a = []
        for i = 1:n
            j = rand(1:10000000000)
            push!(a, j, i)
        end
        print("Inserting random $n values into BinarySearchTree: \t")
        @btime for i = 1:n
            insert!(t, a[i], i)
        end
        print("Base.find call to find a number in array: \t\t")
        @btime begin
            find(x -> x == 15, a)
        end
        @assert parentvalidity(t)
        @assert bstvalidity(t)
        print("delete 50% from BinarySearchTree and add again: \t")
        @btime begin
            for i = 1:div(n, 2)
                delete!(t, a[i])
            end
            for i = 1:div(n, 2)
                insert!(t, a[i], i)
            end
        end
    end

    begin
        t = RBTree{Int, Int}()
        t.unique = false
        n = 10000
        a = []
        for i = 1:n
            j = rand(1:10000000000)
            push!(a, j, i)
        end
        print("Inserting $n values into RBTree: \t\t\t")
        @btime for i = 1:n
            insert!(t, a[i], i)
        end
        print("Base.find call to find a number in array: \t\t")
        @btime find(x -> x == 15, a)
        @assert parentvalidity(t)
        @assert bstvalidity(t)
        print("delete 50% from RBTree and add again: \t\t\t")
        @btime begin
            for i = 1:div(n, 2)
                delete!(t, a[i])
            end
            for i = 1:div(n, 2)
                insert!(t, a[i], i)
            end
        end
    end

    begin
        t = RandomizedBinarySearchTree{Int, Int}()
        t.unique = false
        n = 1000
        a = []
        for i = 1:n
            j = rand(1:10000000000)
            push!(a, j, i)
        end
        print("Inserting random $n values into RandomizedBinarySearchTree: \t")
        @btime for i = 1:n
            insert!(t, a[i], i)
        end
        print("Base.find call to find a number in array: \t\t")
        @btime begin
            find(x -> x == 15, a)
        end
        @assert countvalidity(t)
        @assert parentvalidity(t)
        @assert bstvalidity(t)
        print("delete 50% from BinarySearchTree and add again: \t")
        @btime begin
            for i = 1:div(n, 2)
                delete!(t, a[i])
            end
            for i = 1:div(n, 2)
                insert!(t, a[i], i)
            end
        end
        @assert countvalidity(t)
        @assert parentvalidity(t)
        @assert bstvalidity(t)
    end
end


text:
function sim_power!(gpT::GPE, gpC::GPE, gpNull::GPE, τ::Float64,
            treat::BitVector, Xb::AbstractMatrix,
                chi_null::Vector{Float64},
                mll_null::Vector{Float64},
                pval_invvar_null::Vector{Float64},
                Σcliff::PDMat,
                cK_T::AbstractMatrix,
                cK_C::AbstractMatrix,
                KCT::AbstractMatrix
               )
    # simulate data
    Ysim = prior_rand(gpNull)
    Ysim[treat] .+= τ

    # Modify data in GP objects
    gpT.y[:] = Ysim[treat]
    gpC.y[:] = Ysim[.!treat]
    gpNull.y[:] = Ysim

    # Update GPs
    update_alpha!(gpT)
    update_alpha!(gpC)
    update_alpha!(gpNull)

    gpT.mll = mll(gpT)
    gpC.mll = mll(gpC)
    gpNull.mll = mll(gpNull)

    # mll
    Δmll = gpT.mll + gpC.mll - gpNull.mll
    pval_mll = mean(mll_null .> Δmll)

    # χ2
    χ2 = chistat(gpT, gpC, Xb, Σcliff, cK_T, cK_C)
    pval_χ2 = mean(chi_null .> χ2)

    # inverse-var
    pval_invvar_obs = pval_invvar_uncalib(gpT, gpC, Xb, Σcliff, cK_T, cK_C)

    invvar_bootcalib = mean(pval_invvar_null .< pval_invvar_obs)
    invvar_calib = pval_invvar_calib(gpT, gpC, Xb, Σcliff, cK_T, cK_C, KCT)

    return (pval_mll,pval_χ2,pval_invvar_obs,invvar_bootcalib,invvar_calib)
end

function nsim_power(gpT::GPE, gpC::GPE, τ::Float64,
                Xb::AbstractMatrix,
                chi_null::Vector{Float64},
                mll_null::Vector{Float64},
                pval_invvar_null::Vector{Float64},
                nsim::Int
               )
    gpT_mod = modifiable(gpT)
    gpC_mod = modifiable(gpC)
    gpNull = make_null(gpT_mod, gpC_mod)

    _, Σcliff = cliff_face(gpT, gpC, Xb)
    Σraw, chol = make_posdef!(copy(Σcliff))
    PDcliff = PDMat(Σraw, chol)
    cK_T = cov(gpT.kernel, Xb, gpT.x)
    cK_C = cov(gpC.kernel, Xb, gpC.x)
    KCT = cov(gpC.kernel, gpC.x, gpT.x)

    treat = BitVector(undef, gpNull.nobs)
    treat[:] .= false
    treat[1:gpT.nobs] .= true
    power_sims = [sim_power!(gpT_mod, gpC_mod, gpNull, τ, treat, Xb,
                            chi_null, mll_null, pval_invvar_null,
                            PDcliff, cK_T, cK_C, KCT
                           )
                  for _ in 1:nsim];
    return power_sims
end


text:
"""
Chain a series of transform, here `t1` will be called first
```
    t1 = ScaleTransform()
    t2 = LowRankTransform(rand(3,4))
    ct = ChainTransform([t1,t2]) #t1 will be called first
    ct == t2∘t1
```
"""
struct ChainTransform <: Transform
    transforms::Vector{Transform}
end

Base.length(t::ChainTransform) = length(t.transforms) #TODO Add test

function ChainTransform(v::AbstractVector{<:Transform})
    ChainTransform(v)
end

function ChainTransform(v::AbstractVector{<:Type{<:Transform}},θ::AbstractVector)
    @assert length(v) == length(θ)
    ChainTransform(v.(θ))
end

function transform(t::ChainTransform,X::T,obsdim::Int=defaultobs) where {T}
    Xtr = copy(X)
    for tr in t.transforms
        Xtr = transform(tr,Xtr,obsdim)
    end
    return Xtr
end

set!(t::ChainTransform,θ) = set!.(t.transforms,θ)
params(t::ChainTransform) = (params.(t.transforms))
duplicate(t::ChainTransform,θ) = ChainTransform(duplicate.(t.transforms,θ))

Base.:∘(t₁::Transform,t₂::Transform) = ChainTransform([t₂,t₁])
Base.:∘(t::Transform,tc::ChainTransform) = ChainTransform(vcat(tc.transforms,t)) #TODO add test
Base.:∘(tc::ChainTransform,t::Transform) = ChainTransform(vcat(t,tc.transforms))