text:
export skip_complete

import Base: show

"""
    skip_complete()

Creates a `skip_complete` operator, which filters out `complete` event by the source Observable by emitting only
`next` and `error` messages.

# Producing

Stream of type `<: Subscribable{L}` where `L` refers to type of source stream

# Examples
```jldoctest
using Rocket

source = from([ 1, 2, 3 ])
subscribe!(source |> skip_complete(), logger())
;

# output
[LogActor] Data: 1
[LogActor] Data: 2
[LogActor] Data: 3

```

See also: [`AbstractOperator`](@ref), [`InferableOperator`](@ref), [`ProxyObservable`](@ref), [`skip_error`](@ref), [`skip_next`](@ref), [`logger`](@ref)
"""
skip_complete() = SkipCompleteOperator()

struct SkipCompleteOperator <: InferableOperator end

function on_call!(::Type{L}, ::Type{L}, operator::SkipCompleteOperator, source) where L
    return proxy(L, source, SkipCompleteProxy())
end

operator_right(operator::SkipCompleteOperator, ::Type{L}) where L = L

struct SkipCompleteProxy <: ActorProxy end

actor_proxy!(::Type{L}, proxy::SkipCompleteProxy, actor::A) where { L, A } = SkipCompleteActor{L, A}(actor)

struct SkipCompleteActor{L, A} <: Actor{L}
    actor :: A
end

on_next!(actor::SkipCompleteActor{L}, data::L) where L = next!(actor.actor, data)
on_error!(actor::SkipCompleteActor, err)               = error!(actor.actor, err)
on_complete!(actor::SkipCompleteActor)                 = begin end

Base.show(io::IO, ::SkipCompleteOperator)         = print(io, "SkipCompleteOperator()")
Base.show(io::IO, ::SkipCompleteProxy)            = print(io, "SkipCompleteProxy()")
Base.show(io::IO, ::SkipCompleteActor{L}) where L = print(io, "SkipCompleteActor($L)")


text:
push!(LOAD_PATH, joinpath(Pkg.dir("PumiInterface"), "src"))
using Base.Test
using PdePumiInterface
using ODLCommonTools
using SummationByParts
using PumiInterface
include("defs.jl")

@testset "----- Testing PdePumiInterface3DG -----" begin
  degree = 1
  Tsbp = Float64
  sbp = getTetSBPOmega(degree=degree, Tsbp=Tsbp)
#  sbp = TetSBP{Tsbp}(degree=degree, internal=true)
  ref_verts = sbp.vtx
  interp_op = SummationByParts.buildinterpolation(sbp, ref_verts.')
  face_verts = SummationByParts.SymCubatures.getfacevertexindices(sbp.cub)
  topo = ElementTopology{3}(face_verts)
  sbpface = TetFace{Tsbp}(degree, sbp.cub, ref_verts)
  shape_type = 2

  dmg_name = ".null"
#  smb_name = "tet1.smb"
  smb_name = "unitcube.smb"

  opts = PdePumiInterface.get_defaults()
  opts["smb_name"] = smb_name
  opts["dmg_name"] = dmg_name
  opts["order"] = degree
  opts["coloring_distance"] = 2
  opts["use_linear_metrics"] = true
  opts["numBC"] = 1
  opts["BC1"] = [0,1,2,3,4,5]

#  interp_op = eye(4)
  mesh = PumiMeshDG3(Float64, sbp, opts, sbpface, topo)
#  mesh = PumiMeshDG3{Float64, typeof(sbpface)}(dmg_name, smb_name, degree, sbp, opts, sbpface, topo)

  @test ( mesh.m_ptr )!=C_NULL
  @test ( mesh.mshape_ptr )!=C_NULL
  @test ( mesh.coordshape_ptr )!=mesh.mshape_ptr
  @test ( mesh.vert_Nptr )!=C_NULL
  @test ( mesh.edge_Nptr )!=C_NULL
  @test ( mesh.face_Nptr )!=C_NULL
  @test ( mesh.el_Nptr )!=C_NULL
  @test ( mesh.numBC )== 1
  @test ( mesh.bndry_geo_nums[1] )== opts["BC1"]

  function checkNumbering(nshape_ptr, dim, cnt)
    for i=1:4
      nodes = apf.countNodesOn(nshape_ptr, apf.simplexTypes[i])
      if i == dim
        @test ( nodes )== cnt
      else
        @test ( nodes )== 0
      end
    end
  end

  for i=1:4
    numbering_shape = apf.getNumberingShape(mesh.entity_Nptrs[i])
    checkNumbering(numbering_shape, i, 1)
  end

  @test ( mesh.numVert )== 8
  @test ( mesh.numEdge )== 19
  @test ( mesh.numFace )== 18
  @test ( mesh.numEl )== 6
  @test ( mesh.numBoundaryFaces )== 12
  @test ( mesh.numInterfaces )== 6
  @test ( mesh.numFacesPerElement )== 4
  @test isapprox( mesh.volume, 1.0) atol=1e-12
  @test ( mesh.numEntitiesPerType )== [mesh.numVert, mesh.numEdge, mesh.numFace, mesh.numEl]
  @test ( mesh.numTypePerElement )== [4, 6, 4, 1]
  @test ( mesh.el_type )== apf.TET
  @test ( mesh.face_type )== apf.TRIANGLE
  @test ( mesh.isDG )== true
  @test ( mesh.dim )== 3

  for i=1:min(length(mesh.faces), length(mesh.edges))
    @test ( mesh.edges[i] )!=mesh.faces[i]
  end

  nodeshape = apf.getNumberingShape(mesh.nodenums_Nptr)
  checkNumbering(nodeshape, 4, 4)
  dofshape = apf.getNumberingShape(mesh.dofnums_Nptr)
  checkNumbering(dofshape, 4, 4)

  @test ( mesh.bndry_offsets[1] )== 1
  @test ( mesh.bndry_offsets[2] )== 13

  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      for k=1:mesh.dim
        @test  mesh.coords[k, j, i]  > 0.99
        @test  mesh.coords[k, j, i]  < 3.01
      end
    end
  end
  dofs_sorted = reshape(mesh.dofs, mesh.numDof)
  sort!(dofs_sorted)
  @test ( unique(dofs_sorted) )== dofs_sorted

  @test ( maximum(mesh.dofs) )== mesh.numEl*mesh.numNodesPerElement*mesh.numDofPerNode
  @test ( minimum(mesh.dofs) )== 1

  for i=1:mesh.numEl
    @test  mesh.sparsity_counts[1, i]  > 0
    @test  mesh.sparsity_counts[2, i]  < 5
  end

  @test ( size(mesh.neighbor_colors, 1) )== 5
  @test ( size(mesh.neighbor_nums, 1) )== 5

  # check interface array
  for i=1:mesh.numInterfaces
    iface_i = mesh.interfaces[i]
    @test  iface_i.elementL  > 0
    @test  iface_i.elementL  < mesh.numEl+1
    @test  iface_i.elementR  > 0
    @test  iface_i.elementR  < mesh.numEl+1
    @test  iface_i.faceL  > 0
    @test  iface_i.faceL  < 5
    @test  iface_i.faceR  > 0
    @test  iface_i.faceR  < 5
    @test  iface_i.orient  > 0
    @test  iface_i.orient  < 4
  end

  # check boundary array
  for i=1:mesh.numBoundaryFaces
    bndry_i = mesh.bndryfaces[i]
    @test  bndry_i.element  > 0
    @test  bndry_i.element  < mesh.numEl + 1
    @test  bndry_i.face  > 0
    @test  bndry_i.face  < 5
  end

  testSurfaceNumbering(mesh, sbp, opts)

  # check mapping interpolation
  # should be constant within an element for straight-sided elements
  for i=1:mesh.numInterfaces
    iface_i = mesh.interfaces[i]
    el_i = iface_i.elementL
    dxidx_el = mesh.dxidx[:, :, 1, el_i]
    jac_el = mesh.jac[:, el_i]
    jac_face = mesh.jac_face[:, i]

    for j=1:mesh.numNodesPerFace
      dxidx_face = mesh.dxidx_face[:, :, j, i]
      for k=1:3
        for p=1:3
          @test isapprox( dxidx_face[p, k], dxidx_el[p, k]) atol=1e-13
        end
      end
      @test isapprox( jac_face[j], jac_el[j]) atol=1e-13
    end
  end  # end loop over interfaces

  for i=1:mesh.numBoundaryFaces
    bndry_i = mesh.bndryfaces[i]
    el_i = bndry_i.element
    dxidx_el = mesh.dxidx[:, :, 1, el_i]
    jac_el = mesh.jac[:, el_i]
    jac_face = mesh.jac_bndry[:, i]

    for j=1:mesh.numNodesPerFace
      dxidx_face = mesh.dxidx_bndry[:, :, j, i]
      for k=1:3
        for p=1:3
          @test isapprox( dxidx_face[p, k], dxidx_el[p, k]) atol=1e-13
        end
      end
      @test isapprox( jac_face[j], jac_el[j]) atol=1e-13
    end
  end  # end loop over interfaces

  # the 2D tests should work for this too
  test_interp_rev(mesh)

  # test update_coords
  coords_orig = copy(mesh.vert_coords)

  # double all coordinates
  coords_i = zeros(Float64, 3, 4)
  for i=1:mesh.numEl
    for j=1:4
      coords_i[1, j] = 2*coords_orig[1, j, i]
      coords_i[2, j] = 2*coords_orig[2, j, i]
      coords_i[3, j] = 2*coords_orig[3, j, i]
    end
    update_coords(mesh, i, coords_i)
  end

  commit_coords(mesh, sbp, opts)

#  PdePumiInterface.getCoordinatesAndMetrics(mesh, sbp)
  for i = 1:length(mesh.vert_coords)
    @test  abs(2*coords_orig[i] - mesh.vert_coords[i])  < 1e-10
  end

  # test vertmap

  nfaces_interior = 0
  for i=1:mesh.numEl
    vertnums_i = mesh.element_vertnums[:, i]
    for j=1:4
      @test  vertnums_i[j]  > 0
      @test  vertnums_i[j]  < mesh.numVert + 1
    end

    # test there are exactly 3 elements the current element shares 2 vertices
    # with
    nfaces = 0
    for j=1:mesh.numEl
      vertnums_j = mesh.element_vertnums[:, j]
      if length(intersect(vertnums_i, vertnums_j)) == 3
        nfaces += 1
        nfaces_interior += 1
      end
    end

    @test  nfaces  < 4
    @test  nfaces  > 0
  end  # end loop i

  @test ( nfaces_interior )== 2*(mesh.numFace - mesh.numBoundaryFaces)

  # test periodic
  @test ( mesh.numPeriodicInterfaces )== 0

  opts["smb_name"] = "tet3_pxz.smb"
  opts["BC1"] = [0,2,4,5,6]
  mesh = PumiMeshDG3(Float64, sbp, opts, sbpface, topo)
#  mesh = PumiMeshDG3{Float64, typeof(sbpface)}(dmg_name, smb_name, degree, sbp, opts, sbpface, topo)

  @test ( mesh.numPeriodicInterfaces )== 18
  for i=1:mesh.numInterfaces
    iface_i = mesh.interfaces[i]
    @test  iface_i.elementL  > 0
    @test  iface_i.elementL  < mesh.numEl + 1
    @test  iface_i.elementR  > 0
    @test  iface_i.elementR  < mesh.numEl + 1
    @test  iface_i.faceL  > 0
    @test  iface_i.faceL  < 5
    @test  iface_i.faceR  > 0
    @test  iface_i.faceR  < 5
  end

  # this mesh test that matched vertices are used for orientation
  # determiniation

  opts["smb_name"] = "tet5_periodic.smb"
  opts["numBC"] = 0
  delete!(opts, "BC1")

  mesh = PumiMeshDG3(Float64, sbp, opts, sbpface, topo)
#  mesh = PumiMeshDG3{Float64, typeof(sbpface)}(dmg_name, smb_name, degree, sbp, opts, sbpface, topo)

  @test ( mesh.numPeriodicInterfaces )== 3*(5*5*2)
  for i=1:mesh.numInterfaces
    @test ( mesh.interfaces[i].orient )!=0
    @test  mesh.interfaces[i].orient  < 4
  end

  for order=1:2
    fshape = apf.getFieldShape(0, order, 3)
    eshape = apf.getEntityShape(fshape, apf.TET)  # triangle
    nodexi = PdePumiInterface.getXiCoords(order, 3)
    numnodes = size(nodexi, 2)
    for i=1:numnodes
      vals = apf.getValues(mesh.m_ptr, eshape, nodexi[:, i], numnodes)
      for j=1:numnodes
        if i == j
          @test  abs(vals[j] - 1)  < 1e-12
        else
          @test  abs(vals[j])  < 1e-12
        end  # end if else
      end   # end loop j
    end  # end loop i

  end  # end loop order

  mesh = PumiMeshDG3(Float64, sbp, opts, sbpface, topo)
#  mesh = PumiMeshDG3{Float64, typeof(sbpface)}(dmg_name, smb_name, degree, sbp, opts, sbpface, topo)

  mesh_c = PumiMeshDG3(Complex128, sbp, opts, sbpface, topo)
#  mesh_c = PumiMeshDG3{Complex128, typeof(sbpface)}(dmg_name, smb_name, degree, sbp, opts, sbpface, topo)

  compare_meshes(mesh, mesh_c)

  # test curvilinear is same as linear for linear mesh
  println("testing curvilinear")
  degree = 2
  Tsbp = Float64
  sbp = getTetSBPOmega(degree=degree, Tsbp=Tsbp)
#  sbp = TetSBP{Tsbp}(degree=degree, internal=true)
#  ref_verts = sbp.vtx
  interp_op = SummationByParts.buildinterpolation(sbp, ref_verts.')
  face_verts = SummationByParts.SymCubatures.getfacevertexindices(sbp.cub)
  topo = ElementTopology{3}(face_verts)
  sbpface = TetFace{Tsbp}(degree, sbp.cub, ref_verts)

  opts = PdePumiInterface.get_defaults()
  opts["dmg_name"] = ".null"
#  smb_name = "tet1.smb"
  opts["smb_name"] = "unitcube.smb"
  opts["coloring_distance"] = 2

  opts["use_linear_metrics"] = true
  opts["order"] = degree
  opts["numBC"] = 1
  opts["BC1"] = [0,1,2,3,4,5]

#  interp_op = eye(4)
  mesh_c = PumiMeshDG3(Complex128, sbp, opts, sbpface, topo)
#  mesh_c = PumiMeshDG3{Complex128, typeof(sbpface)}(dmg_name, smb_name, degree, sbp, opts, sbpface, topo)
  mesh = PumiMeshDG3(Float64, sbp, opts, sbpface, topo)
#  mesh = PumiMeshDG3{Float64, typeof(sbpface)}(dmg_name, smb_name, degree, sbp, opts, sbpface, topo)

  dxidx_orig = copy(mesh.dxidx)
  jac_orig = copy(mesh.jac)
  nrm_bndry_orig = copy(mesh.nrm_bndry)
  nrm_face_orig = copy(mesh.nrm_face)
  coords_orig = copy(mesh.coords)
  coords_bndry_orig = copy(mesh.coords_bndry)
  coords_interface_orig = copy(mesh.coords_interface)

  # use curvilinear routines to do calculation
  PdePumiInterface.getMeshCoordinates(mesh, sbp)
  PdePumiInterface.getFaceCoordinatesAndNormals(mesh, sbp)
 # PdePumiInterface.getCurvilinearCoordinatesAndMetrics(mesh, sbp)

  # verify metrics are constant for a uniform linear mesh

  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      @test isapprox( norm(mesh.dxidx[:, :, j, i] - mesh.dxidx[:, :, 1, i]), 0.0) atol=1e-12
    end
  end

  for i=1:mesh.numBoundaryFaces
    for j=1:mesh.numNodesPerFace
      @test isapprox( norm(mesh.nrm_bndry[:, j, i] - mesh.nrm_bndry[:, 1, i]), 0.0) atol=1e-12
    end
  end

  for i=1:mesh.numInterfaces
    for j=1:mesh.numNodesPerFace
      @test isapprox( norm(mesh.nrm_face[:, j, i] - mesh.nrm_face[:, 1, i]), 0.0) atol=1e-12
    end
  end

  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      for k=1:mesh.dim
        for p=1:mesh.dim
          @test isapprox( mesh.dxidx[p, k, j, i], dxidx_orig[p, k, j, i]) atol=1e-12
        end
      end

      @test isapprox( mesh.jac[j, i], jac_orig[j, i]) atol=1e-12
    end
  end

  for i=1:mesh.numBoundaryFaces
    for j=1:mesh.numNodesPerFace
      for k=1:mesh.dim
        @test isapprox( mesh.nrm_bndry[k, j, i], nrm_bndry_orig[k, j, i]) atol=1e-12
        @test isapprox( mesh.coords_bndry[k, j, i], coords_bndry_orig[k, j, i]) atol=1e-12
      end
    end
  end

  for i=1:mesh.numInterfaces
    for j=1:mesh.numNodesPerFace
      for k=1:mesh.dim
        @test isapprox( mesh.nrm_face[k, j, i], nrm_face_orig[k, j, i]) atol=1e-12
        @test isapprox( mesh.coords_interface[k, j, i], coords_interface_orig[k, j, i]) atol=1e-12
      end
    end
  end

  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      for d=1:mesh.dim
        @test isapprox( mesh.coords[d, j, i], coords_orig[d, j, i]) atol=1e-13
      end
    end
  end

  testSurfaceNumbering(mesh, sbp, opts)

  # test reverse mode
  PdePumiInterface.copy_data!(mesh_c, mesh)
  test_metrics_rev(mesh, mesh_c, sbp, opts)

  # test with a diagonal E operator
  degree = 2
  Tsbp = Float64
  sbp = getTetSBPDiagE(degree=degree, Tsbp=Tsbp)
#  sbp = TetSBP{Tsbp}(degree=degree, internal=true)
#  ref_verts = sbp.vtx
  interp_op = SummationByParts.buildinterpolation(sbp, ref_verts.')
  face_verts = SummationByParts.SymCubatures.getfacevertexindices(sbp.cub)
  topo = ElementTopology{3}(face_verts)
  sbpface = getTetFaceForDiagE(degree, sbp.cub, ref_verts)
#  sbpface = TetFace{Tsbp}(degree, sbp.cub, ref_verts)

  opts["dmg_name"] = ".null"
#  smb_name = "tet1.smb"
  opts["smb_name"] = "unitcube.smb"

  opts["use_linear_metrics"] = false

#  interp_op = eye(4)

  mesh2_c = PumiMeshDG3(Complex128, sbp, opts, sbpface, topo; shape_type=4)
#  mesh2_c = PumiMeshDG3{Complex128, typeof(sbpface)}(dmg_name, smb_name, degree, sbp, opts, sbpface, topo, shape_type=4)

  mesh2 = PumiMeshDG3(Float64, sbp, opts, sbpface, topo, shape_type=4)
#  mesh2 = PumiMeshDG3{Float64, typeof(sbpface)}(dmg_name, smb_name, degree, sbp, opts, sbpface, topo, shape_type=4)

  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      @test isapprox( norm(mesh.dxidx[:, :, j, i] - mesh2.dxidx[:, :, 1, i]), 0.0) atol=1e-12
    end
  end

  for i=1:mesh.numBoundaryFaces
    for j=1:mesh.numNodesPerFace
      @test isapprox( norm(mesh.nrm_bndry[:, j, i] - mesh2.nrm_bndry[:, 1, i]), 0.0) atol=1e-12
    end
  end

  for i=1:mesh.numInterfaces
    for j=1:mesh.numNodesPerFace
      @test isapprox( norm(mesh.nrm_face[:, j, i] - mesh2.nrm_face[:, 1, i]), 0.0) atol=1e-12
    end
  end

  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      for k=1:mesh.dim
        for p=1:mesh.dim
          @test isapprox( mesh.dxidx[p, k, j, i], mesh2.dxidx[p, k, 1, i]) atol=1e-12
        end
      end

      @test isapprox( mesh.jac[j, i], mesh2.jac[j, i]) atol=1e-12
    end
  end

  for i=1:mesh.numBoundaryFaces
    for j=1:mesh.numNodesPerFace
      for k=1:mesh.dim
        @test isapprox( mesh.nrm_bndry[k, j, i], nrm_bndry_orig[k, 1, i]) atol=1e-12
      end
    end
  end

  for i=1:mesh.numInterfaces
    for j=1:mesh.numNodesPerFace
      for k=1:mesh.dim
        @test isapprox( mesh.nrm_face[k, j, i], nrm_face_orig[k, 1, i]) atol=1e-12
      end
    end
  end

  # test reverse mode
  PdePumiInterface.copy_data!(mesh2_c, mesh2)
  test_metrics_rev(mesh2, mesh2_c, sbp, opts)
  test_metrics_rev_1d(mesh2_c, sbp, opts)

  test_adapt_3d()

  test_geoNums(mesh)

  if HAVE_SIMMETRIX
    opts = PdePumiInterface.get_defaults()
    opts["smb_name"] = "./meshes/UnitCube/UnitCube.smb"
    opts["dmg_name"] = "./meshes/UnitCube/UnitCube.x_t"
    opts["order"] = degree
    opts["coloring_distance"] = 2
    opts["numBC"] = 1
    opts["BC1"] = [85, 83, 81, 79, 1, 43]

    edge_verts = [1 2 1 1 2 3;  # TODO: SBP should provide this
                  2 3 3 4 4 4]

    topo2 = ElementTopology2()
    topo = ElementTopology{3}(face_verts, edge_verts, topo2=topo2)

    mesh = PumiMeshDG3(Float64, sbp, opts, sbpface, topo, shape_type=4)
    test_setPoint(mesh, opts)

    mesh = PumiMeshDG3(Float64, sbp, opts, sbpface, topo, shape_type=4)
    test_geoNums(mesh)
    test_geoDerivative(mesh)

    opts["smb_name"] = "./meshes/UnitCubeCurve/UnitCubeCurve.smb"
    opts["dmg_name"] = "./meshes/UnitCubeCurve/UnitCubeCurve.x_t"

    mesh = PumiMeshDG3(Float64, sbp, opts, sbpface, topo, shape_type=4)
    test_setPoint(mesh, opts)

    mesh = PumiMeshDG3(Float64, sbp, opts, sbpface, topo, shape_type=4)
    test_geoDerivative(mesh)

    ctor = (T, sbp, opts, sbpface)-> PumiMeshDG3(T, sbp, opts, sbpface, topo, shape_type=4)
    test_coords_file(opts, sbp, sbpface, ctor)

  end

end

ext: jl


text:
portfolio_insts = [
    [(num_stocks, false, true)
    for num_stocks in vcat(10, 1000, 2000:2000:20000)] # includes compile run
    ]

insts = OrderedDict()
insts["nat"] = (nothing, portfolio_insts)
insts["ext"] = (:SOCExpPSD, portfolio_insts)
return (PortfolioJuMP, insts)


text:
@testset "`Samples` API" begin
    root = mktempdir()
    signals = Signal[]
    possible_recordings = (uuid4(), uuid4(), uuid4())
    expected_sample_types = (UInt8, UInt16, UInt32, UInt64, Int8, Int16, Int32, Int64, Float32, Float64)
    expected_parameters = [(recording=rand(possible_recordings),
                            file_path=joinpath(root, "x_$(i)_file"),
                            file_format=rand(("lpcm", "lpcm.zst")),
                            kind="x_$i",
                            channels=["a_$i", "b_$i", "c_$i"],
                            sample_unit="unit_$i",
                            sample_resolution_in_unit=rand((0.25, 1)),
                            sample_offset_in_unit=rand((-0.25, 0.25)),
                            sample_type=expected_sample_types[i],
                            sample_rate=rand((128, 50.5)),
                            start=Second(rand(0:30))) for i in 1:length(expected_sample_types)]
    for params in expected_parameters
        info = Onda.extract_samples_info(params)
        data = rand(params.sample_type, 3, sample_count(info, Second(rand(1:15))))
        samples = Samples(data, info, true)
        start = Second(rand(0:30))
        signal = store(params.file_path, params.file_format, samples, params.recording, params.start)
        push!(signals, signal)
    end
    for (expected, signal) in zip(expected_parameters, signals)
        expected_info = Onda.extract_samples_info(expected)
        for encoded in (false, true)
            s = @compat load(signal; encoded)
            s1 = @compat load(expected.file_path, expected.file_format, expected_info; encoded)
            s2 = @compat load(expected.file_path, Onda.format(expected.file_format, expected_info), expected_info; encoded)
            @test s == s1 == s2
            @test s.info == expected_info
            @test TimeSpans.istimespan(s)
            @test TimeSpans.duration(s) == TimeSpans.time_from_index(s.info.sample_rate, size(s.data, 2) + 1)
            @test channel_count(s) == channel_count(s.info) == length(s.info.channels)
            @test sample_count(s) == sample_count(s.info, TimeSpans.duration(s)) == size(s.data, 2)
            encoded || continue # everything after this assumes `encoded` is `true`
            if signal.file_format == "lpcm"
                s1_mmap = Onda.mmap(signal)
                s2_mmap = Onda.mmap(signal.file_path, s.info)
                s3_mmap = open(io -> Onda.mmap(io, s.info), signal.file_path)
                @test s == s1_mmap == s2_mmap == s3_mmap
            else
                @test_throws ArgumentError Onda.mmap(signal)
            end
            @test sizeof(Matrix(s.data)) == sizeof_samples(s.info, TimeSpans.duration(s))
            @test encode(s) === s
            tmp = similar(s.data)
            tmp_dither_storage = zeros(size(s.data))
            encode!(tmp, s.info.sample_resolution_in_unit,
                    s.info.sample_offset_in_unit, decode(s).data,
                    tmp_dither_storage)
            @test tmp == encode(sample_type(s.info), s.info.sample_resolution_in_unit,
                                s.info.sample_offset_in_unit, decode(s).data + tmp_dither_storage,
                                nothing)
            tmp = similar(s.data)
            encode!(tmp, s)
            @test tmp == s.data
            d = decode(s)
            @test decode(d) === d
            tmp = similar(d.data)
            decode!(tmp, d)
            @test tmp == d.data
            tmp = similar(d.data)
            decode!(tmp, s)
            @test tmp == d.data
            tmp = similar(d.data)
            decode!(tmp, s.info.sample_resolution_in_unit, s.info.sample_offset_in_unit, s.data)
            @test tmp == d.data
            @test d.data == (s.data .* s.info.sample_resolution_in_unit .+ s.info.sample_offset_in_unit)
            if sizeof(sample_type(s.info)) >= 8
                # decoding from 64-bit to floating point is fairly lossy
                tmp = similar(s.data)
                @test isapprox(encode(d, nothing).data, s.data, rtol = 10)
                encode!(tmp, d, nothing)
                @test isapprox(tmp, s.data, rtol = 10)
                @test isapprox(encode(d, missing).data, s.data, rtol = 10)
                encode!(tmp, d, missing)
                @test isapprox(tmp, s.data, rtol = 10)
            else
                tmp = similar(s.data)
                encode!(tmp, d, nothing)
                @test tmp == s.data
                @test encode(d, nothing).data == s.data
                encode!(tmp, d, missing)
                @test isapprox(tmp, s.data, rtol = 1)
                @test isapprox(encode(d, missing).data, s.data, rtol = 1)
            end
            chs = s.info.channels
            i = 27
            @test s[:, i].data == s.data[:, i:i]
            t = TimeSpans.time_from_index(s.info.sample_rate, i)
            t2 = TimeSpans.time_from_index(s.info.sample_rate, i + 15)
            j = TimeSpans.index_from_time(s.info.sample_rate, t2) - 1
            for (ch_inds, reg) in ((:, r""), (1:2, r"[ab]"), (2:3, r"[bc]"), (1:3, r"[abc]"), ([3,1], :skip), ([2,3,1], :skip), ([1,2,3], r"[abc]"))
                @test s[chs[ch_inds], t].data == s[ch_inds, i].data
                @test s[chs[ch_inds], TimeSpan(t, t2)].data == s.data[ch_inds, i:j]
                @test s[chs[ch_inds], i:j].data == s.data[ch_inds, i:j]
                @test s[ch_inds, t].data == s[ch_inds, i].data
                @test s[ch_inds, TimeSpan(t, t2)].data == s.data[ch_inds, i:j]
                @test s[ch_inds, i:j].data == s.data[ch_inds, i:j]
                reg === :skip && continue # can't represent the index with regex
                @test s[reg, t].data == s[ch_inds, i].data
                @test s[reg, TimeSpan(t, t2)].data == s.data[ch_inds, i:j]
                @test s[reg, i:j].data == s.data[ch_inds, i:j]
            end
            @test size(s[:, TimeSpan(0, Second(1))].data, 2) == floor(s.info.sample_rate)
            for i in 1:length(chs)
                @test channel(s, chs[i]) == i
                @test channel(s, i) == chs[i]
                @test channel(s.info, chs[i]) == i
                @test channel(s.info, i) == chs[i]
            end
            @test s[:, TimeSpan(0, TimeSpans.duration(s))].data == s.data
        end
    end
end

@testset "`Samples` pretty printing" begin
    info = SamplesInfo(kind="eeg",
                       channels=["a", "b", "c-d"],
                       sample_unit="unit",
                       sample_resolution_in_unit=0.25,
                       sample_offset_in_unit=-0.5,
                       sample_type=Int16,
                       sample_rate=50.2)
    samples = Samples(rand(Random.MersenneTwister(0), sample_type(info), 3, 5), info, true)
    M = VERSION >= v"1.6" ? "Matrix{Int16}" : "Array{Int16,2}"
    @test sprint(show, samples, context=(:compact => true)) == "Samples(3×5 $M)"
    @test sprint(show, samples) == """
                                   Samples (00:00:00.099601594):
                                     info.kind: "eeg"
                                     info.channels: ["a", "b", "c-d"]
                                     info.sample_unit: "unit"
                                     info.sample_resolution_in_unit: 0.25
                                     info.sample_offset_in_unit: -0.5
                                     sample_type(info): Int16
                                     info.sample_rate: 50.2 Hz
                                     encoded: true
                                     data:
                                   3×5 $M:
                                    20032  4760  27427  -20758   24287
                                    14240  5037   5598   -5888   21784
                                    16885   600  20880  -32493  -19305"""
end

@testset "Onda.VALIDATE_SAMPLES_DEFAULT" begin
    info = SamplesInfo(kind="kind",
                       channels=["a", "b", "c"],
                       sample_unit="microvolt",
                       sample_resolution_in_unit=1.0,
                       sample_offset_in_unit=0.0,
                       sample_type=Int16,
                       sample_rate=100.0)
    @test Onda.VALIDATE_SAMPLES_DEFAULT[]
    @test_throws ArgumentError Samples(rand(4, 10), info, false)
    @test_throws ArgumentError Samples(rand(Int32, 3, 10), info, true)
    Onda.VALIDATE_SAMPLES_DEFAULT[] = false
    @test Samples(rand(4, 10), info, false) isa Samples
    @test Samples(rand(Int32, 3, 10), info, true) isa Samples
    @test_throws ArgumentError Onda.validate(Samples(rand(4, 10), info, false))
    @test_throws ArgumentError Onda.validate(Samples(rand(Int32, 3, 10), info, true))
    Onda.VALIDATE_SAMPLES_DEFAULT[] = true
    @test_throws ArgumentError Samples(rand(4, 10), info, false)
    @test_throws ArgumentError Samples(rand(Int32, 3, 10), info, true)
end


text:
"""
Riemannian quadrature on the left.
"""
function RiemannQuadratureLeft(T=Float64)
    QuadratureRule(1, T[0], T[1])
end

"""
Riemannian quadrature on the right.
"""
function RiemannQuadratureRight(T=Float64)
    QuadratureRule(1, T[1], T[1])
end

"""
Midpoint quadrature.
"""
function MidpointQuadrature(T=Float64)
    QuadratureRule(2, T[1//2], T[1])
end

"""
Trapezoidal quadrature.
"""
function TrapezoidalQuadrature(T=Float64)
    QuadratureRule(2, T[0, 1], T[1//2, 1//2])
end


text:
using Drawing

with(Axes(centred=true), Pen(0.2)) do
    paint(Ink("yellow")) do
        circle(1)           # face background, painted yellow
    end
    draw() do
        circle(1)           # face outline, drawn black (default)
        circle(0.5; from=200deg, to=340deg)    # smile
    end
    paint() do
        move(0.3, 0.25)     # right eye position
        circle(0.2)         # paint eye
        move(-0.3, 0.25)    # left eye position
        circle(0.2)         # paint eye
    end
end


text:
"""
    compact_table!(table::CleanTable; empty_values::Vector=[])

Removes in-place from a `CleanTable` all rows and columns filled entirely by `missing` and
`empty_values`.
"""
function compact_table!(table::CleanTable; empty_values::Vector=[])
    return compact_rows!(
        compact_columns!(table; empty_values=empty_values); empty_values=empty_values
    )
end

"""
    compact_table(table; empty_values::Vector=[])

Creates a `CleanTable` with copied columns and removes from it all rows and columns filled
entirely by `missing` and `empty_values`.
"""
function compact_table(table; empty_values::Vector=[])
    return compact_table!(CleanTable(table); empty_values=empty_values)
end

"""
    compact_columns!(table::CleanTable; empty_values::Vector=[])

Removes in-place from a `CleanTable` all columns filled entirely by `missing` and
`empty_values`.
"""
function compact_columns!(table::CleanTable; empty_values::Vector=[])
    columns = cols(table)
    ndel = 0

    if Threads.nthreads() > 1 && length(columns) > 1 && length(columns[1]) >= 1_000_000
        col_state = replace!(Vector{Bool}(undef, length(columns)), true => false)

        Threads.@threads for i in 1:length(columns)
            if _is_empty_col(columns[i]; empty_values=empty_values)
                col_state[i] = true
            end
        end

        to_delete = findall(col_state)
        for j in to_delete
            deleteat!(columns, j - ndel)
            deleteat!(names(table), j - ndel)
            ndel += 1
        end
    else
        for i in 1:length(columns)
            if _is_empty_col(columns[i - ndel]; empty_values=empty_values)
                deleteat!(columns, i - ndel)
                deleteat!(names(table), i - ndel)
                ndel += 1
            end
        end
    end

    return table
end

"""
    compact_columns(table; empty_values::Vector=[])

Creates a `CleanTable` with copied columns and removes from it all columns filled entirely by
`missing` and `empty_values`.
"""
function compact_columns(table; empty_values::Vector=[])
    return compact_columns!(CleanTable(table); empty_values=empty_values)
end

function _is_empty_col(col; empty_values::Vector=[])
    for el in skipmissing(col)
        isempty(empty_values) && return false

        !in(el, empty_values) && return false
    end

    return true
end

"""
    compact_rows!(table::CleanTable; empty_values::Vector=[])

Removes in-place from a `CleanTable` all rows filled entirely by `missing` and
`empty_values`.
"""
function compact_rows!(table::CleanTable; empty_values::Vector=[])
    columns = cols(table)
    nrows = length(columns[1])
    row_state = replace!(Vector{Bool}(undef, nrows), true => false)

    if Threads.nthreads() > 1 && length(columns) > 1 && length(columns[1]) >= 1_000_000
        Threads.@threads for col in columns
            i = 1

            while i <= nrows
                if !row_state[i] && !ismissing(col[i]) && !in(col[i], empty_values)
                    row_state[i] = true
                end
                i += 1
            end
        end

        to_delete = findall(!, row_state)
        Threads.@threads for col in columns
            for num in 0:(length(to_delete) - 1)
                deleteat!(col, to_delete[num + 1] - num)
            end
        end
    else
        for col in columns
            i = 1

            while i <= nrows
                if !row_state[i] && !ismissing(col[i]) && !in(col[i], empty_values)
                    row_state[i] = true
                end
                i += 1
            end
        end

        to_delete = findall(!, row_state)
        for col in columns
            for num in 0:(length(to_delete) - 1)
                deleteat!(col, to_delete[num + 1] - num)
            end
        end
    end

    return table
end

"""
    compact_rows(table; empty_values::Vector=[])

Creates a `CleanTable` with copied columns and removes from it all rows filled entirely by
`missing` and `empty_values`.
"""
function compact_rows(table; empty_values::Vector=[])
    return compact_rows!(CleanTable(table); empty_values=empty_values)
end

max_stars_repo_name: TheRoniOne/Cleaner

max_issues_repo_name: TheRoniOne/Cleaner

max_forks_repo_name: TheRoniOne/Cleaner


text:
function (SL::SingleLayerKernel{T,S})(x,y)::T  where {T,S<:Maxwell}
    N = ambient_dimension(S)
    x==y && return zero(T)
    k = SL.op.k
    r = x - y
    d = norm(r)
    if N==2
        return error("Maxwell operator not implemented in 2d")
    elseif N==3
        g   = 1/(4π)/d * exp(im*k*d)
        gp  = im*k*g - g/d
        gpp = im*k*gp - gp/d + g/d^2
        ID    = SMatrix{3,3,Float64,9}(1,0,0,0,1,0,0,0,1)
        RRT   = r*transpose(r) # r ⊗ rᵗ
        return  g*ID + 1/k^2*(gp/d*ID + (gpp/d^2 - gp/d^3)*RRT)
    end
end

# Double Layer Kernel
# n × ∇ × G = γ₁ G
function (DL::DoubleLayerKernel{T,S})(x,y,ny)::T where {T,S<:Maxwell}
    N = ambient_dimension(S)
    x==y && return zero(T)
    k = DL.op.k
    r = x - y
    d = norm(r)
    g   = 1/(4π)/d * exp(im*k*d)
    gp  = im*k*g - g/d
    if N==2
        return error("Maxwell operator not yet defined in 2d")
    elseif N==3
        ID    = SMatrix{3,3,Float64,9}(1,0,0,0,1,0,0,0,1)
        ncross = transpose(SMatrix{3,3,Float64,9}(0,-ny[3],ny[2],
                                         ny[3],0,-ny[1],
                                         -ny[2],ny[1],0))
        rcross = transpose(SMatrix{3,3,Float64,9}(0,-r[3],r[2],
                                         r[3],0,-r[1],
                                         -r[2],r[1],0))
        # return -gp/d*ncross*rcross
        return gp/d*rcross*ncross
    end
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
# This file contains the implementation of various
# Monte Carlo samplers.

using SteinDistributions
using Distributions: MvNormal, TDist, ccdf

# This function generates a lambda function
# in the form of a power law. This is useful
# for computing the stepsize for algorithms
# that require some changing step size.
function epsilon_generator(;a=1.0, b=1.0, gamma=0.6)
    t::Int64 -> a * (b + t) ^ -gamma;
end

# A composition of the int and floor functions.
function intf(x::Float64)
    round(Int, floor(x))
end

# This function runs the approximate Metropolis-Hastings
# sampling procedure as described in Korattikara, Chen,
# and Welling.
#
# Inputs:
# d                - the target distribution to approximate a sample from
# theta0           - the starting point for the sampler
# numiter          - the number of iterations to run the sampler (i.e. number of points)
# numlikelihood    - the number of likelihood evaluations to make before halting
# batchsize        - the number of observations to evaluate before doing the significance test
# epsilon          - the threshold for assessing the significance test
# proposalvariance - the variance of the mean-zero Gaussian proposal distribution
#
# Returns:
# (samples, number of gradient evaluations)
function runapproxmh(d::SteinPosterior,
                     theta0::Array{Float64};
                     numiter=Inf,
                     numlikelihood=Inf,
                     batchsize=10,
                     epsilon=0.1,
                     proposalvariance=1e-2)

    if isinf(numiter) && isinf(numlikelihood)
        error("Need to specify numlikelihood or numiter")
    end

    p = length(theta0)
    t = 0
    noise = MvNormal(p, 1)
    thetas = Array(Float64, 0, p)
    theta_t = theta0
    numlikelihoods = Int64[]
    cumlikelihoods = 0

    while t < numiter && cumlikelihoods < numlikelihood
        delta = sqrt(proposalvariance) .* vec(rand(noise, 1))
        prop_t = theta_t + delta
        u = rand()

        (acceptprop, currbatchsize) =
            _speedyproposaltest(d, prop_t, theta_t, u, epsilon, batchsize)

        t += 1
        cumlikelihoods += currbatchsize
        push!(numlikelihoods, cumlikelihoods)

        if acceptprop
            theta_t = prop_t
        end
        thetas = [thetas; theta_t']
    end
    thetas, numlikelihoods
end

# This function creates samples via the method outlined in
# Teh and Welling called Stochastic Gradient Langevin Dynamics.
#
# Inputs:
# d                - the target distribution to approximate a sample from
# theta0           - the starting point for the sampler
# batchsize        - the size of the subsamples to evaluate the gradients
# numsweeps        - the number of passes to make over the entire dataset
# epsilonfunc      - the function generating epsilon as described in the algorithm
#
# Returns:
# (samples, number of gradient evaluations)
function runsgld(d::SteinPosterior,
                 theta0::Array{Float64};
                 batchsize=1,
                 numsweeps=100,
                 epsilonfunc=epsilon_generator())

    L = numdatapoints(d)
    p = length(theta0)
    noise = MvNormal(p, 1)
    t = 0
    batchratio = L / batchsize

    theta_t = theta0
    thetas = Array(Float64, numsweeps * intf(batchratio), p)

    for k = 1:numsweeps
        perm = randperm(L)
        for i = 1:intf(batchratio)
            j = batchsize * (i-1) + 1
            idx = perm[j:(j+batchsize-1)]

            epsilon_t = epsilonfunc(t)
            delta = (epsilon_t / 2) .* gradlogdensity(d, theta_t; idx=idx)
            delta += sqrt(epsilon_t) .* vec(rand(noise, 1))
            theta_t += delta
            t += 1
            thetas[t,:] = theta_t
        end
    end

    numgradients = numsweeps * intf(batchratio) * batchsize
    thetas, numgradients
end

# This function creates samples via a variant of SGLD that uses the
# Riemannian-Langevin operator as opposed to the Langevin one.
# The update is of the form
#
# Q(.|x) = N(x + eps * [a(x) grad log p(x) + <grad, a^T>(x)], eps a(x))
#
# Inputs:
# d                          - the target distribution to approximate a sample from
# theta0                     - the starting point for the sampler
# volatility_covariance      - the matrix-valued function a = (1/2) sigma sigma^T
# grad_volatility_covariance - the vector-valued function grad_a = <grad, a^t>
# batchsize                  - the size of the subsamples to evaluate the gradients
# numsweeps                  - the number of passes to make over the entire dataset
# epsilonfunc                - the function generating epsilon as described in the algorithm
# Returns:
# (samples, number of gradient evaluations)
function runsgrld(d::SteinPosterior,
                  theta0::Array{Float64,1},
                  volatility_covariance::Function,
                  grad_volatility_covariance::Function;
                  batchsize=1,
                  numsweeps=100,
                  epsilonfunc=epsilon_generator())
    L = numdatapoints(d)
    p = length(theta0)
    t = 0
    batchratio = L / batchsize

    theta_t = theta0
    thetas = Array(Float64, numsweeps * intf(batchratio), p)

    for k = 1:numsweeps
        perm = randperm(L)
        for i = 1:intf(batchratio)
            j = batchsize * (i-1) + 1
            idx = perm[j:(j+batchsize-1)]
            # get the mass matrix ready
            epsilon_t = epsilonfunc(t)
            a = volatility_covariance(theta_t)
            grad_a = grad_volatility_covariance(theta_t)
            # compute change!
            delta = a * gradlogdensity(d, theta_t; idx=idx) + grad_a
            delta *= epsilon_t
            # next compute the variance term
            noise = MvNormal(a)
            delta += sqrt(epsilon_t) * vec(rand(noise, 1))
            theta_t += delta
            t += 1
            thetas[t,:] = theta_t
        end
    end

    numgradients = numsweeps * intf(batchratio) * batchsize
    thetas, numgradients
end

# This function creates samples via Stochastic Gradient Fisher Scoring
# described in Korattikara, Welling and Ahn.
#
# Inputs:
# d                - the target distribution to approximate a sample from
# theta0           - the starting point for the sampler
# batchsize        - the size of the subsamples to evaluate the gradients
# numsweeps        - the number of passes to make over the entire dataset
# epsilon          - the parameter controlling the stepsizes and variance of noise
# diagonal         - if true, run SGFS-d (see the paper)
# kappa            - the function yielding kappa_t as described in the algorithm
#
# Returns:
# (samples, number of gradient evaluations)
function runsgfs(d::SteinPosterior,
                 theta0::Array{Float64};
                 batchsize=100,
                 numsweeps=100,
                 epsilon=1e-3,
                 diagonal=false,
                 kappa=(t -> 1/t))
    L = numdatapoints(d)
    p = length(theta0)
    gamma = (batchsize + L) / batchsize

    t = 1
    It = zeros(p, p)
    batchratio = L / batchsize

    theta_t = theta0;
    thetas = Array(Float64, numsweeps * intf(batchratio), p);

    for k = 1:numsweeps
        perm = randperm(L)
        for i = 1:intf(batchratio)
            j = batchsize * (i-1) + 1;
            idx = perm[j:(j+batchsize-1)]

            g = Array(Float64, 0, p)
            for yidx in idx
                g = vcat(g, gradloglikelihood(d, theta_t; idx=yidx))
            end
            V = cov(g)

            It = (1 - kappa(t)) .* It + kappa(t) .* V
            B = gamma * L .* It

            noise = MvNormal(B)
            eta = (2/sqrt(epsilon)) .* vec(rand(noise, 1))
            if diagonal
                Dvec = gamma * L * diag(It) + (4/epsilon) * diag(B)
                Dinv = diagm(1 ./ Dvec)
            else
                Dinv = inv(gamma * L * It + (4/epsilon) * B)
            end

            delta =
                2 .* Dinv *
                (gradlogdensity(d, theta_t; idx=idx) + eta)

            theta_t += delta
            thetas[t,:] = theta_t
            t += 1
        end
    end

    numgradients = numsweeps * intf(batchratio) * batchsize
    thetas, numgradients
end

# This function creates samples via Metropolis Adjusted Langevin Algorithm.
#
# Inputs:
# d           - the target distribution to approximate a sample from
# theta0      - the starting point for the sampler
# numiter     - the number of samples to create
# epsilonfunc - the function computing epsilon, which controls the
#               variance of the noise and step size
# Returns:
# (samples, number of gradient evaluations)
function runmala(d::SteinDistribution,
                 theta0::Array{Float64};
                 numiter=100,
                 epsilonfunc=epsilon_generator())

    L = numdatapoints(d)
    p = length(theta0)
    noise = MvNormal(p, 1)
    thetas = Array(Float64, numiter, p)
    theta_t = theta0

    for t = 1:numiter
        epsilon_t = epsilonfunc(t-1);
        delta = epsilon_t .* gradlogdensity(d, theta_t) ./ 2
        delta += sqrt(epsilon_t) .* vec(rand(noise, 1))

        prop_t = theta_t + delta
        logdensitydiff = logdensity(d, prop_t) - logdensity(d, theta_t)
        if (log(rand()) < logdensitydiff)
            theta_t = prop_t
        end
        thetas[t,:] = theta_t
    end
    thetas, numiter * L
end

# This function creates samples via Metropolis Adjusted Langevin Algorithm
# on a manifold. This can be thought of as a discretization of Riemannian
# Langevin diffusions. The update is of the form
#
# Q(.|x) = N(x + eps * [a(x) grad log p(x) + <grad, a^T>(x)], eps a(x))
#
# Inputs:
# d                          - the target distribution to approximate a sample from
# theta0                     - the starting point for the sampler
# volatility_covariance      - the matrix-valued function a = (1/2) sigma sigma^T
# grad_volatility_covariance - the vector-valued function grad_a = <grad, a^t>
# numiter                    - the number of samples to create
# epsilonfunc                - the function computing epsilon, which controls the
#                              variance of the noise and step size
# Returns:
# (samples, number of gradient evaluations, acceptance ratio)
function runmarla(d::SteinDistribution,
                  theta0::Array{Float64,1},
                  volatility_covariance::Function,
                  grad_volatility_covariance::Function;
                  numiter::Int=100,
                  epsilonfunc::Function=epsilon_generator(),
                  verbose::Bool=false)

    L = numdatapoints(d)
    p = length(theta0)
    thetas = Array(Float64, numiter, p)
    theta_t = theta0
    numaccepts = 0

    for t = 1:numiter
        epsilon_t = epsilonfunc(t-1)::Float64
        a = volatility_covariance(theta_t)::Array{Float64,2}
        grad_a = grad_volatility_covariance(theta_t)::Array{Float64,1}
        # first compute the drift part
        delta = a * gradlogdensity(d, theta_t)::Array{Float64,1} + grad_a
        delta *= epsilon_t
        # next compute the variance term
        noise = MvNormal(a)
        delta += sqrt(epsilon_t) * vec(rand(noise, 1))
        # then define the new proposal
        prop_t = theta_t + delta
        logdensitydiff = logdensity(d, prop_t) - logdensity(d, theta_t)
        if (log(rand()) < logdensitydiff)
            theta_t = prop_t
            numaccepts += 1
        end
        thetas[t,:] = theta_t

        if verbose && (t % 10000) == 0
            println("Finished iteration $t. Continuing.")
        end
    end
    acceptance_ratio = (numaccepts / numiter)
    thetas, 2 * numiter * L, acceptance_ratio
end

# This function creates samples via the approx slice sampler:
#
# http://jmlr.org/proceedings/papers/v33/dubois14.pdf
#
# Inputs:
# d                - the target distribution to approximate a sample from
# theta0           - the starting point for the sampler
# epsilon - the threshold for the sequential test to stop
# numiter - the number of points to sample before halting
# numlikelihood - the maximum number of likelihoods to use before halting
# batchsize - the size of each batch included during the sequential test
# w - the initial length of the line between L and R
# Vmax - the max number of times the initial [L, R] will be extended
# Returns:
# (samples, number of gradient evaluations)
function runapproxslice(d::SteinDistribution,
                        theta0::Array{Float64,1};
                        epsilon=0.1,
                        numiter=Inf,
                        numlikelihood=Inf,
                        batchsize=30,
                        w=1.0,
                        Vmax=100)
    if isinf(numiter) && isinf(numlikelihood)
        error("Need to specify numlikelihood or numiter")
    end

    p = length(theta0)
    thetas = Array(Float64, 0, p)
    theta_t = theta0
    numlikelihoods = Int64[]
    cumlikelihoods = 0
    t = 0

    while t < numiter && cumlikelihoods < numlikelihood
        u, wbar, vbar = [rand() for _ in 1:3]
        # pick some direction
        j = round(Int, ceil(p * rand()))
        z = zeros(p); z[j] = 1
        # now initialize the boundaries: L and R
        L = theta_t - w * wbar * z
        R = L + w * z
        # set VL and VR, needed for detailed balance
        VL = intf(Vmax * vbar)
        VR = Vmax - 1 - VL
        # search for the slice
        onslice = true
        while VL > 0 && onslice
            (onslice, likelihoods)  =
                _speedyproposaltest(d, L, theta_t, u, epsilon, batchsize)
            cumlikelihoods += likelihoods
            VL -= 1
            L -= (w * z)
        end
        onslice = true
        while VR > 0 && onslice
            (onslice, likelihoods)  =
                _speedyproposaltest(d, R, theta_t, u, epsilon, batchsize)
            cumlikelihoods += likelihoods
            VR -= 1
            R += (w * z)
        end
        # now pick the next candidate!
        while true
            eta = rand()
            thetaprop = L + eta * (R - L)
            (onslice, likelihoods) =
                _speedyproposaltest(d, thetaprop, theta_t, u, epsilon, batchsize)
            cumlikelihoods += likelihoods
            if onslice
                theta_t = thetaprop
                break
            end
            if dot(thetaprop, z) < dot(theta_t, z)
                L = thetaprop
            else
                R = thetaprop
            end
        end
        t += 1
        thetas = [thetas; theta_t']
        push!(numlikelihoods, cumlikelihoods)
    end

    thetas, numlikelihoods
end

# This does the sequential test to that runs until the new proposal
# (thetaprop) is significantly different than the acceptance threshold.  The
# sequential test to see if the probability P(thetaptop) / P(theta) > u,
# where u is drawn from Unif([0,1]).
#
# @params
# @d         - the distribution to be sampled from
# @thetaprop - the theta proposed to be on the slice
# @theta     - the current theta defining the slice
# @u         - the Unif([0,1]) that is the height in the sampler
# @epsilon   - the significance level of the test to cut short
#              the complete computation. If 0.0, this is exact.
# @batchsize - the number of samples to add at each test
# @returns: (is ratio larger, number of likelihood evals)
function _speedyproposaltest(d::SteinPosterior,
                             thetaprop::Array{Float64,1},
                             theta::Array{Float64,1},
                             u::Float64,
                             epsilon::Float64,
                             batchsize::Int)
    N = numdatapoints(d)
    logu = log(u)
    lprior = logprior(d, thetaprop) - logprior(d, theta)
    mu0 = (logu - lprior) / N

    significancelevel = 1.0
    prevbatchend = 1
    currbatchsize = 0
    ldiffmean = 0.0
    perm = randperm(N)
    ldiffs = Float64[]

    while (currbatchsize < N) && (significancelevel > epsilon)
        currbatchsize = min(currbatchsize + batchsize, N)
        for yidx in perm[prevbatchend:currbatchsize]
            push!(
                ldiffs, (
                    loglikelihood(d, thetaprop; idx=yidx:yidx) -
                    loglikelihood(d, theta; idx=yidx:yidx)
                )
            )
        end
        ldiffmean = mean(ldiffs)
        lstd = std(ldiffs) * sqrt(1.0/currbatchsize - 1.0/N)

        tstat = (ldiffmean - mu0) / lstd
        tdist = TDist(currbatchsize - 1)

        significancelevel = ccdf(tdist, abs(tstat))
        prevbatchend += batchsize
    end

    ldiffmean > mu0, currbatchsize
end


text:
using Documenter
using Literate
using ShallowWaters

EXAMPLE = joinpath(@__DIR__, "assets", "swm_equations.jl")
OUTPUT = joinpath(@__DIR__, "src")

# Generate markdown
binder_badge = "# [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/milankl/ShallowWaters.jl/gh-pages?labpath=dev%2Fswm_equations.ipynb)"
function preprocess_docs(content)
  return string(binder_badge, "\n\n", content)
end

Literate.markdown(EXAMPLE, OUTPUT; preprocess=preprocess_docs, codefence="```julia" => "```")
Literate.notebook(EXAMPLE, OUTPUT)

pages = [
  "Introduction" => "index.md",
  "Tutorial" => "swm_equations.md",
  "Reference" => "reference.md",
]

makedocs(
  sitename = "ShallowWaters.jl",
  format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
  modules = [ShallowWaters],
  pages = pages,
)

deploydocs(
  repo = "github.com/milankl/ShallowWaters.jl.git",
  push_preview = true,
  devbranch = "main",
)


text:
# Fractional Integration --
# Author: Ziyi Yin, ziyi.yin@gatech.edu
# Date: July 2020

export apply_frac_integral, frac_int

###############################################################################
# Fractional integration

function frac_int(nt::Int,order::Number;DDT=Float32)
	# fractional integration for a trace
	x    = zeros(Int64,nt,1);
	x[1] = 1;
	y    = convert(Array{DDT,2},fgl_deriv(-order,x,1));	# fractional integration
	C    = joConvolve(nt,1,y[:];DDT=DDT,RDT=DDT);
	return C
end

function apply_frac_integral(C::joLinearFunction{ADDT,ARDT},d_obs::judiVector{avDT}) where {ADDT,ARDT,avDT}
    # fractional integration based on C
    C.n == size(d_obs.data[1],1) || throw(judiVectorException("shape mismatch"))
    jo_check_type_match(ADDT,avDT,join(["DDT for *(joLinearFunction,judiWeights):",C.name,typeof(C),avDT]," / "))
	d_out = deepcopy(d_obs)
	for i = 1:d_obs.nsrc
		d_out.data[i] = C*d_obs.data[i]
	end
	return d_out
end


text:
# @testset "fit moments" begin
# end;

# @testset "fit to quantilepoint and mode - ignore mode" begin
# end;

# @testset "fit two quantiles same" begin
# end;

@testset "fit two quantiles" begin
    D = Weibull(1,1)
    #Plots.plot(D)
    median(D)
    qpl = @qp_m(median(D))
    qpu = @qp_u(quantile(D, 0.95))
    d = fit(Weibull, qpl, qpu);
    @test d == D
    #
    qpl = @qp_m(0.5)
    qpu = @qp_uu(5)
    d = fit(Weibull, qpl, qpu);
    @test median(d) == 0.5
    @test quantile(d, 0.975) ≈ 5
    #Plots.plot(d);
end;


text:
let
    using Random

    # Square matrix
    Random.seed!(0)
    p = 4; q = 5; m = 8; n = m;
    Ut = randn(p,m)
    Vt = randn(p,n)
    Pt = randn(q,m)
    Qt = randn(q,n)
    x = randn(n,1)
    A = tril(Ut'*Vt) + triu(Pt'*Qt,1)
    yref = A*x
    y = egrss.gemv(Ut,Vt,Pt,Qt,x)
    @test ( isapprox(y,yref) )

    # Fat matrix
    Random.seed!(0)
    p = 3; q = 4; m = 8; n = m+3;
    Ut = randn(p,m)
    Vt = randn(p,n)
    Pt = randn(q,m)
    Qt = randn(q,n)
    x = randn(n,1)
    A = tril(Ut'*Vt) + triu(Pt'*Qt,1)
    yref = A*x
    y = egrss.gemv(Ut,Vt,Pt,Qt,x)
    @test ( isapprox(y,yref) )
    y = egrss.gemv(Ut[:,1:m],Vt,Pt,Qt,x)
    @test ( isapprox(y,yref) )

    # Tall matrix
    Random.seed!(0)
    p = 3; q = 2; m = 8; n = m-2;
    Ut = randn(p,m)
    Vt = randn(p,n)
    Pt = randn(q,m)
    Qt = randn(q,n)
    x = randn(n,1)
    A = tril(Ut'*Vt) + triu(Pt'*Qt,1)
    yref = A*x
    y = egrss.gemv(Ut,Vt,Pt,Qt,x)
    @test ( isapprox(y,yref) )
    y = egrss.gemv(Ut,Vt,Pt[:,1:n],Qt,x);
    @test ( isapprox(y,yref) )

end

max_stars_repo_name: martinandersen/EGRSS

max_issues_repo_name: martinandersen/EGRSS

max_forks_repo_name: martinandersen/EGRSS


text:
function parametrized_kuramoto(du,u,p,t)
    # NumberOfOscillators, CouplingConstant, and NaturalFreqs
    N,K,w = p
    # Centroid
    u[1] = 0
    for i in 1:N
        u[1] += cos(u[i+1]) + sin(u[i+1])*im
    end
    u[1] = u[1]/N
    r = abs(u[1])
    ψ = angle(u[1])
    # Oscillators
    for i in 1:N
    	du[i+1] = w[i] + K*r*sin(ψ-u[i+1])
    end
end

function getparameters_kuramoto(N,K,γ,distribution)
    distribution == "Uniform" ? nothing : error("$(distribution) is not yet supported")
    j = collect(1:N)
    w = broadcast(j->((2j-N-1)/(N-1))*((pi*γ)/2),j)
    p = (N=N, K=K, w=w)
    return p
end

function getstartpoints_kuramoto(N)
    d = Uniform(-pi,pi)
    u0::Array{ComplexF64} = rand(d,N)
    pushfirst!(u0,mean(u0))
    return u0
end

function kuramoto(N,K,T,fs,γ,distribution)
    tspan = (0.0,T)
    u0 = getstartpoints_kuramoto(N)
    p = getparameters_kuramoto(N,K,γ,distribution)
    prob = ODEProblem(parametrized_kuramoto,u0,tspan,p)
    sol = solve(prob,RK4(),saveat=collect(0:1/fs:T),progress=true)
end


text:
struct Order_Report
    converged_K::Vector{Float64}
    converged_lh::Vector{Float64}
    failed_K::Vector{Float64}
    failed_lh::Vector{Float64}
end

struct Partition_Report
    partition_id::String
    naive_lh::Float64
    orddict::Dict{Integer,Order_Report}
    best_model::Tuple{Chain_ID,BHMM}
    best_repset::Vector{Chain_ID}
end

function Base.show(io::IO, report::Partition_Report)
    printstyled(io, "BACKGROUND HMM TRAINING REPORT\n", bold=true)
    printstyled(io, "Genome partition id: $(report.partition_id)\n", bold=true, color=:magenta)

    for (order, ordreport) in report.orddict
        if length(ordreport.failed_lh)>0
            ordplot=scatterplot(ordreport.converged_K, ordreport.converged_lh; name="Converged", title="Order $order HMMs", xlabel="K# States", ylabel="P(O|θ)", color=:green, xlim=(1,maximum(ordreport.converged_K)), ylim=(floor(min(minimum(ordreport.converged_lh),minimum(ordreport.failed_lh),report.naive_lh),sigdigits=2),ceil(max(maximum(ordreport.converged_lh),maximum(ordreport.failed_lh),report.naive_lh),sigdigits=2)))
            length(ordreport.failed_K)>=1 && scatterplot!(ordplot, ordreport.failed_K, ordreport.failed_lh; name="Unconverged",color=:red)
            lineplot!(ordplot, [report.naive_lh for i in 1:maximum(ordreport.converged_K)], name="Naive",color=:magenta)
            show(ordplot)
            println()
        else
            ordplot=scatterplot(ordreport.converged_K, ordreport.converged_lh; name="Converged", title="Order $order HMMs", xlabel="K# States", ylabel="P(O|θ)", color=:green, xlim=(1,maximum(ordreport.converged_K)), ylim=(floor(min(minimum(ordreport.converged_lh),report.naive_lh),sigdigits=2),ceil(max(maximum(ordreport.converged_lh),report.naive_lh),sigdigits=2)))
            lineplot!(ordplot, [report.naive_lh for i in 1:maximum(ordreport.converged_K)], name="Naive",color=:magenta)
            show(ordplot)
            println()
        end
    end
end

function report_partitions(chain_reports::Dict{Chain_ID,Chain_Report})
    partitions=Vector{String}()
    orders=Vector{Integer}()
    Ks=Vector{Integer}()
    reps=Vector{Integer}()
    for id in keys(chain_reports)
        !in(id.obs_id, partitions) && push!(partitions, id.obs_id)
        !in(id.order, orders) && push!(orders, id.order)
        !in(id.K, Ks) && push!(Ks, id.K)
        !in(id.replicate, reps) && push!(reps, id.replicate)
    end

    reports=Vector{Partition_Report}()
    for partition in partitions
        best_lh,best_rep=-Inf,Chain_ID("empty",1,0,1)
        naive_lh=1.
        orddict=Dict{Integer,Order_Report}()
        for order in orders
            conv_statevec=Vector{Float64}()
            fail_statevec=Vector{Float64}()
            conv_lh_vec=Vector{Float64}()
            fail_lh_vec=Vector{Float64}()
            for (id,report) in chain_reports
                if id.obs_id==partition && id.order==order
                    report.test_lh > best_lh && (best_lh=report.test_lh;best_rep=id)
                    naive_lh==1. && (naive_lh=chain_reports[id].naive_lh)
                    chain_reports[id].converged ? (push!(conv_statevec,id.K);
                    push!(conv_lh_vec,chain_reports[id].test_lh)) : (push!(fail_statevec,id.K);
                    push!(fail_lh_vec,chain_reports[id].test_lh))
                end
            end
            orddict[order]=Order_Report(conv_statevec,conv_lh_vec,fail_statevec,fail_lh_vec)
        end
        best_model=(best_rep,chain_reports[best_rep].final_hmm)
        best_repset=[Chain_ID(partition, best_rep.K, best_rep.order, rep) for rep in reps]
        push!(reports,Partition_Report(partition,naive_lh,orddict,best_model,best_repset))
    end
    return reports
end


text:
"""
Graphic outputs that display frames as RGB24 images.
"""
abstract type ImageOutput{T} <: GraphicOutput{T} end

(::Type{F})(o::T; kwargs...) where F <: ImageOutput where T <: ImageOutput =
    F(; frames=frames(o), starttime=starttime(o), endtime=endtime(o),
      fps=fps(o), showfps=showfps(o), timestamp=timestamp(o), stampframe=stampframe(o), store=store(o),
      processor=processor(o), minval=minval(o), maxval=maxval(o),
      kwargs...)

"""
Mixin for outputs that output images and can use an image processor.
"""
@premix @default_kw struct Image{P,Mi,Ma}
    processor::P | ColorProcessor()
    minval::Mi   | 0.0
    maxval::Ma   | 1.0
end

processor(o::ImageOutput) = o.processor
minval(o::ImageOutput) = o.minval
maxval(o::ImageOutput) = o.maxval

showframe(frame, o::ImageOutput, data::AbstractSimData, f) =
    showframe(frame, o, ruleset(data), f)
showframe(frame, o::ImageOutput, ruleset::AbstractRuleset, f) =
    showframe(frametoimage(o, ruleset, frame, f), o, f)

# Manual showframe without data/ruleset
showframe(o::ImageOutput, f=lastindex(o)) = showframe(o[f], o, Ruleset(), f)

"""
Default colorscheme. Better performance than using a Colorschemes.jl scheme.
"""
struct Greyscale{M1,M2}
    max::M1
    min::M2
end
Greyscale(; min=nothing, max=nothing) = Greyscale(max, min)

Base.get(scheme::Greyscale, x) = scale(x, scheme.min, scheme.max)

scale(x, ::Nothing, max) = x * max
scale(x, min, ::Nothing) = x * (one(min) - min) + min
scale(x, ::Nothing, ::Nothing) = x
scale(x, min, max) = x * (max - min) + min

const Grayscale = Greyscale

# Greyscale is the default processor
processor(x) = Greyscale()

"""
Frame processors convert arrays into RGB24 images for display.
"""
abstract type FrameProcessor end

abstract type MultiFrameProcessor <: FrameProcessor end

"""
Convert frame matrix to RGB24, using an FrameProcessor
"""
function frametoimage end

frametoimage(o::ImageOutput, i::Integer) = frametoimage(o, o[i], i)
frametoimage(o::ImageOutput, frame, i::Integer) =
    frametoimage(processor(o), o, Ruleset(), o[i], i)
frametoimage(o::ImageOutput, args...) = frametoimage(processor(o), o, args...)

""""
Converts output frames to a colorsheme.
# Arguments
`scheme`: a ColorSchemes.jl colorscheme.
"""
struct ColorProcessor{S,Z,M} <: FrameProcessor
    scheme::S
    zerocolor::Z
    maskcolor::M
end
ColorProcessor(; scheme=Greyscale(), zerocolor=nothing, maskcolor=nothing) =
    ColorProcessor(scheme, zerocolor, maskcolor)

scheme(processor::ColorProcessor) = processor.scheme
zerocolor(processor::ColorProcessor) = processor.zerocolor
maskcolor(processor::ColorProcessor) = processor.maskcolor

frametoimage(p::ColorProcessor, o::Output,
             ruleset::AbstractRuleset, frame::AbstractArray, t) = begin
    img = fill(RGB24(0), size(frame))
    for i in CartesianIndices(frame)
        img[i] = if !(maskcolor(p) isa Nothing) && ismasked(mask(ruleset), i)
            maskcolor(p)
        else
            x = frame[i]
            if hasminmax(o) isa HasMinMax
                x = normalise(x, minval(o), maxval(o))
            end
            if !(zerocolor(p) isa Nothing) && x == zero(x)
                zerocolor(p)
            else
                rgb(scheme(p), x)
            end
        end
    end
    img
end

abstract type BandColor end

struct Red <: BandColor end
struct Green <: BandColor end
struct Blue <: BandColor end

struct ThreeColor{C<:Tuple,Z,M} <: MultiFrameProcessor
    colors::C
    zerocolor::Z
    maskcolor::M
end
ThreeColor(; colors=(Red(), Green(), Blue()), zerocolor=nothing, maskcolor=nothing) =
    ThreeColor(colors, zerocolor, maskcolor)

colors(processor::ThreeColor) = processor.colors
zerocolor(processor::ThreeColor) = processor.zerocolor
maskcolor(processor::ThreeColor) = processor.maskcolor

frametoimage(p::ThreeColor, o::Output, ruleset, bands::NamedTuple, t) = begin
    img = fill(RGB24(0), size(first(bands)))
    ncols = length(colors(p))
    nbands = length(bands)
    nbands == ncols || throw(ArgumentError("$nbands layers in model but $ncols colors"))

    for i in CartesianIndices(first(bands))
        img[i] = if !(maskcolor(p) isa Nothing) && ismasked(mask(ruleset), i)
            maskcolor(p)
        else
            if hasminmax(o) isa HasMinMax
                xs = map((f, mi, ma) -> normalise(f[i], mi, ma), values(bands), minval(o), maxval(o))
            else
                xs = map(f -> f[i], values(bands))
            end
            if !(zerocolor(p) isa Nothing) && all(map(x -> x .== zero(x), xs))
                zerocolor(p)
            else
                acc = (0.0, 0.0, 0.0)
                combine(colors(p), acc, xs)
            end
        end
    end
    img
end

combine(c::Tuple{Red,Vararg}, acc, xs) = combine(tail(c), (acc[1] + xs[1], acc[2], acc[3]), tail(xs))
combine(c::Tuple{Green,Vararg}, acc, xs) = combine(tail(c), (acc[1], acc[2] + xs[1], acc[3]), tail(xs))
combine(c::Tuple{Blue,Vararg}, acc, xs) = combine(tail(c), (acc[1], acc[2], acc[3] + xs[1]), tail(xs))
combine(c::Tuple{Nothing,Vararg}, acc, xs) = combine(tail(c), acc, tail(xs))
combine(c::Tuple{}, acc, xs) = RGB24(acc...)

rgb(scheme, x) = RGB24(get(scheme, x))
rgb(c::RGB24) = c
rgb(c::Tuple) = RGB24(c...)
rgb(c::Number) = RGB24(c)

"""
    savegif(filename::String, o::Output, ruleset; [processor=processor(o)], [kwargs...])

Write the output array to a gif. You must pass a processor keyword argument for any
`Output` objects not in `ImageOutput` (which allready have a processor attached).

Saving very large gifs may trigger a bug in Imagemagick.
"""
savegif(filename::String, o::Output, ruleset=Ruleset();
        processor=processor(o), minval=0, maxval=1, kwargs...) = begin
    # fr = normaliseframe.(frames(o), minval, maxval)
    images = frametoimage.(Ref(processor), Ref(o), Ref(ruleset), frames(o), collect(firstindex(o):lastindex(o)))
    array = cat(images..., dims=3)
    FileIO.save(filename, array; kwargs...)
end

struct HasMinMax end
struct NoMinMax end

hasminmax(output::T) where T = (:minval in fieldnames(T)) ? HasMinMax() : NoMinMax()

normaliseframe(output::Output, a) =
    normaliseframe(hasminmax(output), output, a)
normaliseframe(::HasMinMax, output, a::NamedTuple) =
    map(normaliseframe, values(a), minval(output), maxval(output))
normaliseframe(::HasMinMax, output, a::AbstractArray) =
    normaliseframe(a, minval(output), maxval(output))
normaliseframe(::NoMinMax, output, a) = a
normaliseframe(a::AbstractArray, minval::Number, maxval::Number) = normalise.(a, minval, maxval)
normaliseframe(a::AbstractArray, minval, maxval) = a

normalise(x::Number, minval::Number, maxval::Number) =
    min((x - minval) / (maxval - minval), oneunit(eltype(x)))


text:
#v.2
function inputStr(mensagem::String)::String # armazenar resultado em txt
    println(mensagem)
    return readline()
end

function gerar_grafo()::Dict{String,Dict{String,Float64}}
    arquivo = open("C:\\Users\\Carlos Eduardo\\arquivo.txt", "r");
    arquivo = split(read(arquivo, String), keepempty = true)

    grafo::Dict{String,Dict{String,Float64}} = Dict()
    i::UInt64, j::UInt64, contador::UInt64 = 1, 1, 0
    aux_no = Array[]
    FloatOuString = Union{Float64,AbstractString}
    push!(aux_no,FloatOuString[])

    for linha in arquivo
        contador += 1

        if linha == ""
            continue
        elseif j % 3 != 0
            push!(aux_no[i], linha)
        else
            linha = parse(Float64, linha)
            push!(aux_no[i], linha)
            if contador < length(arquivo) # limitar a quantidade de array
                push!(aux_no,FloatOuString[])
            end
            i += 1
        end

        j += 1
    end

    qtd_no::UInt64 = i-1

    for i = 1:qtd_no
        grafo[aux_no[i][1]] = Dict()
        grafo[aux_no[i][2]] = Dict()
    end

    for i = 1:qtd_no
        grafo[aux_no[i][1]][aux_no[i][2]] = aux_no[i][3]
        grafo[aux_no[i][2]][aux_no[i][1]] = aux_no[i][3]
    end

    #close(arquivo);
    return grafo
end

function dijkstra(no_inicial::String, no_final::String, grafo::Dict{String,Dict{String,Float64}})::Nothing
    nao_visitado = grafo
    menor_distancia::Dict{String,Float64} = Dict()
    no_anterior::Dict{String,String} = Dict()
    caminho = String[]
    infinito = Inf

    [menor_distancia[no] = infinito for (no,) in nao_visitado]

    menor_distancia[no_inicial] = 0.0
    no = no_inicial

    while length(nao_visitado) != 0
        aux_grafo = nao_visitado[no]

        for (aux_no,peso) = aux_grafo
            if peso + menor_distancia[no] < menor_distancia[aux_no]
                menor_distancia[aux_no] = peso + menor_distancia[no]
                no_anterior[aux_no] = no
            end
        end

        pop!(nao_visitado,no)
        no = nothing

        for (aux_no,peso) in nao_visitado
            if no == nothing
                no = aux_no
            elseif menor_distancia[aux_no] < menor_distancia[no]
                no = aux_no
            end
        end

        if no == no_final
            break
        end
    end

    no_atual = no_final

    while no_atual != no_inicial
        pushfirst!(caminho,no_atual)
        no_atual = no_anterior[no_atual]
    end

    pushfirst!(caminho,no_inicial)

    println("Peso mínimo: $(menor_distancia[no_final])")
    println("Caminho mínimo: $caminho")
    return nothing
end

function main()
    no_inicial = inputStr("Digite o nó inicial: ")
    no_final = inputStr("Digite o nó final: ")

    grafo = gerar_grafo()

    println("")

    dijkstra(no_inicial, no_final, grafo)
end

@time main()


text:
@testset "Implicit" begin
   function h!(out::A,x::B,p::C) where {A,B,C}
       out[1] = x[1]^2 + x[1]*p[1] + 4.0
       return
   end
   function hj!(out::A,x::B,p::C) where {A,B,C}
       out[1,1] = 2.0*x[1]+p[1]
       return
   end
   nx = 1
   np = 1
   P = [Interval{Float64}(6.0,9.0)]
   X = [Interval{Float64}(-0.78,-0.4)]
   pmid = mid.(P)
   pref_mc = [MC{1,NS}(pmid[1],P[1],1)]
   p_mc = [MC{1,NS}(pmid[1],P[1],1)]

   mccall! = MCCallback(h!, hj!, nx, np)
   mccall!.P[:] = P
   mccall!.X[:] = X
   mccall!.p_mc[1] = MC{1,NS}(pmid[1],P[1],1)
   mccall!.pref_mc[1] = MC{1,NS}(pmid[1],P[1],1)

   gen_expansion_params!(mccall!)
   implicit_relax_h!(mccall!)

   @test isapprox(mccall!.x_mc[1].cv, -0.6307841, atol = 1E-5)
   @test isapprox(mccall!.x_mc[1].cc, -0.5209112, atol = 1E-5)
   @test isapprox(mccall!.x_mc[1].Intv.lo, -0.7720045045045048, atol = 1E-5)
   @test isapprox(mccall!.x_mc[1].cv_grad[1], 0.09414689079323225, atol = 1E-5)
   @test isapprox(mccall!.x_mc[1].cc_grad[1], 0.0455312, atol = 1E-5)

   dense = McCormick.DenseMidInv(nx, np)
   @test dense.Y == zeros(Float64,nx,nx)
   @test dense.YInterval == zeros(Interval{Float64},1)
   @test dense.nx == nx
   @test dense.np == np

   t = KrawczykCW()
   dmc = MCCallback(h!, hj!, nx, np, KrawczykCW())
   McCormick.contract!(t, dmc)
end


text:
import Base: parse

global LVTYPES = Dict(
    "Cluster" => Dict,
    "Array"   => Array,
    "Refnum"  => String,
    "DBL"     => Float64,
    "I64"     => Int64,
    "I32"     => Int32,
    "I16"     => Int16,
    "I8"      => Int8,
    "U64"     => UInt64,
    "U32"     => UInt32,
    "U16"     => UInt16,
    "U8"      => UInt8,
    "String"  => String,
    "Boolean" => Bool,
)

Base.parse(::Type{T}, x::AbstractString) where T <: AbstractString = x

function parse_array(element)
    # find size of the array and parse as Int. This only works for 1D arrays!
    dimsize = LightXML.find_element(element, "Dimsize") |> LightXML.content |> x -> parse(Int, x)
    lvtype = get_array_eltype(element)
    jltype = LVTYPES[lvtype]
    # preallocate Array
    array = Array{jltype}(undef, dimsize)
    if dimsize > 0
        for (i, _el) in enumerate(LightXML.child_elements(element))
            i < 3 && continue
            value_str = LightXML.find_element(_el, "Val") |> LightXML.content
            value     = parse(jltype, value_str)
            array[i-2] = value
        end
    end
    array
end

function get_array_eltype(element)
    for (i,_e) in enumerate(LightXML.child_elements(element))
        i < 3 && continue
        if i == 3
            return LightXML.name(_e)
        end
    end
    error("No type information found in array $(LightXML.name(element))")
end

"""
    readxml(path::String) -> Dict{String,Any}

Reads a LabVIEW generated XML file and returns a `Dict`.
"""
function readxml(path::String)
    function loopelements(element, dict)
        for e in LightXML.child_elements(element)
            if haskey(LVTYPES, LightXML.name(e))
                if LightXML.name(e) == "Cluster"
                    # Go deeper if we hit a cluster
                    clustername = LightXML.find_element(e, "Name") |> LightXML.content
                    dict[clustername] = Dict{String,Any}()
                    loopelements(e, dict[clustername])
                elseif LightXML.name(e) == "Array"
                    arrayname = LightXML.find_element(e, "Name") |> LightXML.content
                    dict[arrayname] = parse_array(e)
                else
                    key, value = parseelement(e)
                    # Put out a warning if the key already exists
                    haskey(dict, key) && @warn "Overwriting key '$key'."
                    dict[key] = value
                end
            elseif any(LightXML.name(e) .== ["NumElts", "Name", "Dimsize"])
                # Ignore entry
                continue
            else
                # Put out a warning if the type can't be found
                @warn "Could not find type for $(LightXML.name(e))"
            end
        end
    end

    function parseelement(element)
        lvtype = LightXML.name(element)
        jltype = LVTYPES[lvtype]
        element_name = LightXML.find_element(element, "Name") |> LightXML.content
        value_str    = LightXML.find_element(element, "Val") |> LightXML.content
        value = parse(jltype, value_str)
        element_name, value
    end

    dict = Dict{String,Any}()

    xdoc = LightXML.parse_file(path)
    xroot = LightXML.root(xdoc)

    loopelements(xroot, dict)

    LightXML.free(xdoc)

    dict
end


text:
function fsc_fastmulti{T}(X::Array{T, 2}, dictSize::Int, optParams::FSCOpts{T})
# Input:
#   X: H*WxN column of image patches
#
# Alg. parameters:
#   mu_max: 1e5
#   tau: [1.01 1.5] larger means faster convergence but sacrifices primal feasibility,
#        smaller means slower convergence

  const mu_max = 1e5;
  const tau = 1.01;

  const Dinit    = optParams.Dinit;
  const beta::T  = optParams.beta;
        mu_s::T  = optParams.mu_s;
        mu_t::T  = optParams.mu_t;
  const max_iter = optParams.max_iter;
  const tol      = optParams.tol;
  const FIXED_D  = optParams.FIXED_D;
  const DEBUG    = optParams.DEBUG;

  (filtLen, numPatches) = size(X);
  assert(!FIXED_D || !isempty(Dinit), "Dinit must be provided when FIXED_D is enabled.");
  assert(isempty(Dinit) || size(X, 1)==size(Dinit, 1), "Number of rows in Dinit and X do not match.");
  assert(isempty(Dinit) || size(Dinit, 2)==dictSize, "Dinit does not have the correct number of columns.");

  # Initialization
  if !FIXED_D;
    if isempty(Dinit);
      Dinit = reshape(mod([convert(T, 1):filtLen*dictSize], 13), (filtLen, dictSize));
    end
    d = Dinit';
    s = copy(d);
    normalize!(s);
    lambda_s = zeros(T, dictSize, filtLen);
  else;
    d = Dinit';
    s = d;
  end

  x = X;
  z = zeros(T, dictSize, numPatches);
  t = zeros(T, dictSize, numPatches);
  lambda_t = zeros(T, dictSize, numPatches);

  @printf("iter=00000");
  if DEBUG;
    debugData = {"x"=>X,
                 "d"=>d, "s"=>s,
                 "lambda_s"=>lambda_s,
                 "mu_s"=>mu_s,
                 "z"=>z, "t"=>t,
                 "lambda_t"=>lambda_t,
                 "mu_t"=>mu_t,
                 "beta"=>beta}
    (prev_obj,rerr,spar) = compute_lagrangian(debugData);
    @printf(": obj=%.10f (rerr=%.10f, spar=%.10f)\n", prev_obj, rerr, spar);
  end

  iter = 0;
  converged = false;
  if !FIXED_D;
    learnData = LearningData(zeros(T, dictSize, filtLen), zeros(T, dictSize, dictSize));
    learnParams = LearningOpts(T); learnParams.max_iter = 1; learnParams.tol = tol;
    if DEBUG;
      learnParams.DEBUG = DEBUG;
      learnParams.debugData = debugData;
    end
  end
  codeData = FastCodingData(zeros(T, dictSize, numPatches), zeros(T, dictSize, dictSize));
  codeParams = FastCodingOpts(T); codeParams.beta = beta; codeParams.max_iter = 1; codeParams.tol = tol;
  if DEBUG;
    codeParams.DEBUG = DEBUG;
    codeParams.debugData = debugData;
  end
  while ~converged;
    iter += 1;

    @printf("\b\b\b\b\b\b\b\b\b\biter=%05d", iter);
    if DEBUG;
      debugData["z"] = z; debugData["t"] = t;
      debugData["lambda_t"] = lambda_t;
      debugData["mu_t"] = mu_t;
      debugData["d"] = d; debugData["s"] = s;
      debugData["lambda_s"] = lambda_s;
      debugData["mu_s"] = mu_s;
      (prev_obj,rerr,spar) = compute_lagrangian(debugData);
      @printf(": obj=%.10f (rerr=%.10f, spar=%.10f)\n", prev_obj, rerr, spar);
    end

    #
    # CODING
    #
    Dt = d;
    codeData.DtD = Dt*Dt';
    codeData.Dtx = Dt*x;
    codeParams.z = z;
    codeParams.t = t;
    codeParams.lambda_t = lambda_t;
    codeParams.mu_t = mu_t;
    (z, t) = fsc_coding_fastmulti(codeData, codeParams);

    #
    # LEARNING
    #
    if !FIXED_D;
      fill!(learnData.ZtZ, 0);
      fill!(learnData.Ztx, 0);
      for j=1:numPatches;
        learnData.ZtZ += z[:, j]*z[:, j]';
        learnData.Ztx += z[:, j]*x[:, j]';
      end
      learnParams.Dinit = d;
      learnParams.s = s;
      learnParams.mu_s = mu_s;
      learnParams.lambda_s = lambda_s;
      (d, s) = fsc_learning_multi(learnData, learnParams);
    end

    # Update the Lagrangian multiplier
    for i=1:dictSize*numPatches; lambda_t[i] += mu_t*(z[i]-t[i]); end
    if !FIXED_D;
      for i=1:dictSize*filtLen; lambda_s[i] += mu_s*(d[i]-s[i]); end
    end
#    debugData["lambda_t"] = lambda_t;
#    debugData["lambda_s"] = lambda_s;
#    (obj,rerr,spar) = compute_lagrangian(debugData);
#    is_mono = (prev_obj>=obj) ? 'O' : 'X';
#    @printf("          L: %c obj=%.10f (rerr=%.10f, spar=%.10f)\n", is_mono, obj, rerr, spar);
#    prev_obj = obj;

    # Update penalty term
    mu_t = min(tau*mu_t, mu_max);
    mu_s = min(tau*mu_s, mu_max);
#    debugData["mu_t"] = mu_t;
#    debugData["mu_s"] = mu_s;
#    (obj,rerr,spar) = compute_lagrangian(debugData);
#    is_mono = (prev_obj>=obj) ? 'O' : 'X';
#    @printf("          M: %c obj=%.10f (rerr=%.10f, spar=%.10f)\n", is_mono, obj, rerr, spar);
#    prev_obj = obj;

    learnConverged = FIXED_D || mean((d-s).^2)<=tol;
    codeConverged = true;
    for j=1:numPatches;
      if mu_t>0 && mean((z[:, j]-t[:, j]).^2)>tol;
        codeConverged = false;
        break
      end
    end
    if iter>=max_iter || (learnConverged && codeConverged);
      converged = true;
    end
  end
  return (s', z, t);
end

max_stars_repo_name: bkong/FSC

max_issues_repo_name: bkong/FSC

max_forks_repo_name: bkong/FSC


text:
@testset "Combined Mean and Variance" begin

    @testset "Two samples" begin
        test_two = (x1,x2) -> begin
            xc = vcat(x1,x2)
            meanc, varc = combined_mean_and_var(x1,x2)
            @test abs(mean(xc) - meanc) < 1e-12
            @test abs(var(xc) - varc) < 1e-12
        end

        x1 = rand(30_000)
        x2 = rand(20_000)
        test_two(x1,x2)

        x1 = rand(ComplexF64, 30_000)
        x2 = rand(ComplexF64, 20_000)
        test_two(x1,x2)
    end

    @testset "Three samples" begin
        test_three = (x1,x2,x3) -> begin
            xc = vcat(x1,x2,x3)
            meanc, varc = combined_mean_and_var(x1,x2,x3)
            @test abs(mean(xc) - meanc) < 1e-12
            @test abs(var(xc) - varc) < 1e-12

            meanc_three, varc_three = combined_mean_and_var_three(x1,x2,x3)
            @test isapprox(meanc, meanc_three)
            @test isapprox(varc, varc_three)
        end

        x1 = rand(30_000)
        x2 = rand(20_000)
        x3 = rand(40_000)
        test_three(x1, x2, x3)

        x1 = rand(ComplexF64, 30_000)
        x2 = rand(ComplexF64, 20_000)
        x3 = rand(ComplexF64, 40_000)
        test_three(x1, x2, x3)
    end

    @testset "N samples" begin
        test_N = (xs...) -> begin
            xc = vcat(xs...)
            meanc, varc = combined_mean_and_var(xs...)
            @test abs(mean(xc) - meanc) < 1e-12
            @test abs(var(xc) - varc) < 1e-12
        end

        lengths = [30_000, 20_000, 40_000]
        N = 5

        xs = [rand(Float64, rand(lengths)) for _ in 1:N]
        test_N(xs...)

        xs = [rand(ComplexF64, rand(lengths)) for _ in 1:N]
        test_N(xs...)
    end

    @testset "N samples (ns, μs, vs)" begin
        test_N_moments = (ns, μs, vs, mean_exact, var_exact) -> begin
            meanc, varc = combined_mean_and_var(ns, μs, vs)
            @test abs(mean_exact - meanc) < 1e-12
            @test abs(var_exact - varc) < 1e-12
        end

        lengths = [30_000, 20_000, 40_000]
        N = 5

        xs = [rand(Float64, rand(lengths)) for _ in 1:N]
        xc = vcat(xs...)
        test_N_moments(length.(xs), mean.(xs), var.(xs), mean(xc), var(xc))

        xs = [rand(ComplexF64, rand(lengths)) for _ in 1:N]
        xc = vcat(xs...)
        test_N_moments(length.(xs), mean.(xs), var.(xs), mean(xc), var(xc))
    end
end

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
## Functions for mean feature map
function apply_mean(X::Array{Float64}, rowcounts::Int64, mom::StatisticalMoments,options)
    # X is the matrix formed by all data inputs. All "apply" functions for feature maps need to fill in X  [rowcount,:] with the mapped features.
   lenresponse = length(mom.moments[1]);
   X[rowcounts,1:lenresponse] = mom.moments[1];
   # NOTE: scientifically it's important to see how well mean(real.(responses),2) performs, and not to use abs. If someone could measure the abs of each response, then they could use any statistic they wanted really. However, many experiments (and theoretical developments) use only mean(real.(responses),2).   #Works better with abs.(responses)
end

function get_num_mean_features(mom::StatisticalMoments, options)
   lenresponse = length(mom.moments[1]);
    return lenresponse;
end

# Defines the feature map
fmmean = FeatureMap(apply_mean,get_num_mean_features,"meanmap");
##end functions for identity feature map.