text:
function init_hopping_matrices(mc::AbstractDQMC)
  p = mc.p
  l = mc.l

  if p.Bfield
    init_peirls_phases(mc)
    init_hopping_matrix_exp_Bfield(mc)
    p.chkr && init_checkerboard_matrices_Bfield(mc)
  else
    init_hopping_matrix_exp(mc)
    p.chkr && init_checkerboard_matrices(mc)
  end
  nothing
end

function init_hopping_matrix_exp(mc::AbstractDQMC)
  p = mc.p
  l = mc.l

  println("Initializing hopping exponentials")
  !p.Bfield || @warn("You should be using `init_hopping_matrix_exp_Bfield()` or set p.Bfield = false!")

  Tx = diagm(0 => fill(-p.mu1,l.sites))
  Ty = diagm(0 => fill(-p.mu2,l.sites))
  hor_nb = [2,4]
  ver_nb = [1,3]

  # Nearest neighbor hoppings
  if p.hoppings != "none"
    @inbounds @views begin
      for src in 1:l.sites
        for nb in hor_nb # horizontal neighbors
          trg = l.neighbors[nb,src]
          Tx[trg,src] += -l.t[1,1] # t_x_hor = tx1
          Ty[trg,src] += -l.t[1,2] # t_y_hor = tx2
        end

        for nb in ver_nb # vertical neighbors
          trg = l.neighbors[nb,src]
          Tx[trg,src] += -l.t[2,1] # t_x_ver = ty1
          Ty[trg,src] += -l.t[2,2] # t_y_ver = ty2
        end
      end
    end
  end

  # Next nearest neighbor hoppings (diagonal)
  if p.Nhoppings != "none"
    @inbounds @views begin
      for src in 1:l.sites
        for nb in 1:4
          trg = l.Nneighbors[nb,src]
          Tx[trg,src] += -l.tN[1,1] # t_x_diag = txy1
          Ty[trg,src] += -l.tN[1,2] # t_y_diag = txy2
        end
      end
    end
  end

  # Next-next nearest neighbor hoppings
  if p.NNhoppings != "none"
    @inbounds @views begin
      for src in 1:l.sites
        for nb in hor_nb # horizontal neighbors
          trg = l.NNneighbors[nb,src]
          Tx[trg,src] += -l.tNN[1,1] # t_x_hor = txx1
          Ty[trg,src] += -l.tNN[1,2] # t_y_hor = txx2
        end

        for nb in ver_nb # vertical neighbors
          trg = l.NNneighbors[nb,src]
          Tx[trg,src] += -l.tNN[2,1] # t_x_ver = tyy1
          Ty[trg,src] += -l.tNN[2,2] # t_y_ver = tyy2
        end
      end
    end
  end

  eTx_minus = exp(-0.5 * p.delta_tau * Tx)
  eTy_minus = exp(-0.5 * p.delta_tau * Ty)
  eTx_plus = exp(0.5 * p.delta_tau * Tx)
  eTy_plus = exp(0.5 * p.delta_tau * Ty)

  if p.opdim == 3
    l.hopping_matrix = cat(Tx, Ty, Tx, Ty, dims=(1,2))
    l.hopping_matrix_exp = cat(eTx_minus,eTy_minus,eTx_minus,eTy_minus, dims=(1,2))
    l.hopping_matrix_exp_inv = cat(eTx_plus,eTy_plus,eTx_plus,eTy_plus, dims=(1,2))
  else # O(1) and O(2)
    l.hopping_matrix = cat(Tx, Ty, dims=(1,2))
    l.hopping_matrix_exp = cat(eTx_minus,eTy_minus, dims=(1,2))
    l.hopping_matrix_exp_inv = cat(eTx_plus,eTy_plus, dims=(1,2))
  end

  return nothing
end

function init_peirls_phases(mc::AbstractDQMC)
  println("Initializing Peirls phases (Bfield)")
  l = mc.l
  L = mc.l.L
  p = mc.p

  B = zeros(2,2) # colidx = flavor, rowidx = spin up,down
  if p.Bfield
    B[1,1] = B[2,2] = 2 * pi / l.sites
    B[1,2] = B[2,1] = - 2 * pi / l.sites
  end

  l.peirls = Matrix{Matrix{Float64}}(undef, 2,2) # colidx = flavor, rowidx = spin up,down
  for f in 1:2 # flv
    for s in 1:2 # spin
      phis = fill(NaN, L, L, L, L)

      for x in 1:L
        for y in 1:L
          xp = mod1(x + 1, L)
          yp = mod1(y + 1, L)
          xm = mod1(x - 1, L)
          ym = mod1(y - 1, L)
          xp2 = mod1(x + 2, L)
          yp2 = mod1(y + 2, L)

          #nn
          if p.hoppings != "none"
            phis[x,y,x,yp] = 0
            phis[x,yp,x,y] = 0

            phis[x,y,xp,y] = - B[s,f] * (y - 1)
            phis[xp,y,x,y] = - phis[x,y,xp,y]
            if y == L
                phis[x,y,x,yp] = B[s,f] * L * (x -1)
                phis[x,yp,x,y] = - phis[x,y,x,yp]
            end
          end

          #nnn
          if p.Nhoppings != "none"
            phis[x,y,xp,yp] = -B[s,f] * (y -1 + 0.5)
            phis[xp,yp,x,y] = -phis[x,y,xp,yp]

            phis[x,y,xm,yp] = B[s,f] * (y -1 + 0.5)
            phis[xm,yp,x,y] = -phis[x,y,xm,yp]

            if y == L
                phis[x,y,xp,yp] = B[s,f] * (L * (x -1) + 0.5)
                phis[xp,yp,x,y] = - phis[x,y,xp,yp]

                phis[x,y,xm,yp] = B[s,f] * (L * (x-1) - 0.5)
                phis[xm,yp,x,y] = - phis[x,y,xm,yp]
            end
          end

          #nnnn
          if p.NNhoppings != "none"
            phis[x,y,x,yp2] = 0
            phis[x,yp2,x,y] = 0

            phis[x,y,xp2,y] = -2 * B[s,f] * (y-1)
            phis[xp2,y,x,y] = -phis[x,y,xp2,y]
            if (y == L) || (y == L - 1)
                phis[x,y,x,yp2] = B[s,f] * L * (x-1)
                phis[x,yp2,x,y] = - phis[x,y,x,yp2]
            end
          end

        end
      end
      l.peirls[s,f] = reshape(permutedims(phis, [2,1,4,3]), (l.sites,l.sites))

    end
  end
end

function init_hopping_matrix_exp_Bfield(mc::AbstractDQMC)
  p = mc.p
  l = mc.l
  H = heltype(mc)

  println("Initializing hopping exponentials (Bfield)")
  p.Bfield || @warn("You should be using `init_hopping_matrix_exp()` or set p.Bfield = true!")

  T = Matrix{Matrix{H}}(undef, 2,2) # colidx = flavor, rowidx = spin up,down

  hor_nb = [2,4]
  ver_nb = [1,3]

  for f in 1:2
    for s in 1:2
      T[s,f] = convert(Matrix{H}, diagm(0 => fill(f == 1 ? -p.mu1 : -p.mu2,l.sites)))

      # Nearest neighbor hoppings
      if p.hoppings != "none"
        @inbounds @views begin
          for src in 1:l.sites
            for nb in hor_nb # horizontal neighbors
              trg = l.neighbors[nb,src]
              T[s,f][trg,src] += - exp(im * l.peirls[s,f][trg,src]) * l.t[1,f] # horizontal
            end

            for nb in ver_nb # vertical neighbors
              trg = l.neighbors[nb,src]
              T[s,f][trg,src] += - exp(im * l.peirls[s,f][trg,src]) * l.t[2,f] # vertical
            end
          end
        end
      end

      # Next nearest neighbor hoppings (diagonal)
      if p.Nhoppings != "none"
        @inbounds @views begin
          for src in 1:l.sites
            for nb in 1:4
              trg = l.Nneighbors[nb,src]
              T[s,f][trg,src] += - exp(im * l.peirls[s,f][trg,src]) * l.tN[1,f]
            end
          end
        end
      end

      # Next-next nearest neighbor hoppings
      if p.NNhoppings != "none"
        @inbounds @views begin
          for src in 1:l.sites
            for nb in hor_nb # horizontal neighbors
              trg = l.NNneighbors[nb,src]
              T[s,f][trg,src] += - exp(im * l.peirls[s,f][trg,src]) * l.tNN[1,f] # horizontal
            end

            for nb in ver_nb # vertical neighbors
              trg = l.NNneighbors[nb,src]
              T[s,f][trg,src] += - exp(im * l.peirls[s,f][trg,src]) * l.tNN[2,f] # vertical
            end
          end
        end
      end

    end
  end

  eT_minus = map(Ti -> exp(-0.5 * p.delta_tau * Ti), T)
  eT_plus = map(Ti -> exp(0.5 * p.delta_tau * Ti), T)

  if p.opdim == 3
    l.hopping_matrix = cat(T[1,1], T[2,2], T[2,1], T[1,2], dims=(1,2))
    l.hopping_matrix_exp = cat(eT_minus[1,1], eT_minus[2,2], eT_minus[2,1], eT_minus[1,2], dims=(1,2))
    l.hopping_matrix_exp_inv = cat(eT_plus[1,1], eT_plus[2,2], eT_plus[2,1], eT_plus[1,2], dims=(1,2))
  else # O(1) and O(2)
    l.hopping_matrix = cat(T[1,1], T[2,2], dims=(1,2))
    l.hopping_matrix_exp = cat(eT_minus[1,1], eT_minus[2,2], dims=(1,2))
    l.hopping_matrix_exp_inv = cat(eT_plus[1,1], eT_plus[2,2], dims=(1,2))
  end

  return nothing
end

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
@testset "nograds" begin
  @test gradient(floor, 1) === nothing
  @test gradient(ceil, 1) === nothing
  @test gradient(round, 1) === nothing
  @test gradient(hash, 1) === nothing
  @test gradient(div, 1, 2) === nothing
end #testset

ext: jl


text:
@testset "STDOUT" begin
    # make sure it compiles and executes
    # TODO: replace this with VT100 tests
    imshow(colorview(RGB, rand(3,2,2))); println()
    imshow(colorview(RGB, rand(3,2,2)), (2,2)); println()
    imshow256(colorview(RGB, rand(3,2,2))); println()
    imshow256(colorview(RGB, rand(3,2,2)), (2,2)); println()
    imshow24bit(colorview(RGB, rand(3,2,2))); println()
    imshow24bit(colorview(RGB, rand(3,2,2)), (2,2)); println()
end

@testset "imshow" begin
    @testset "non colorant" begin
        @test_throws ArgumentError imshow(rand(5,5))
        @test_throws ArgumentError imshow(sprand(5,5,.5))
    end
    @testset "lena" begin
        img = imresize(lena, 10, 10)
        io = IOBuffer()
        ensurecolor(imshow256, io, img)
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))
        @test_reference "reference/lena_big_imshow.txt" res
    end
    @testset "ndarray" begin
        img = rgb_line_4d
        io = IOBuffer()
        ensurecolor(imshow256, io, img)
        res = readlines(seek(io,0))
        @test_reference "reference/ndarray_imshow.txt" res
    end
end

@testset "imshow256" begin
    @testset "non colorant" begin
        @test_throws ArgumentError imshow256(rand(5,5))
        @test_throws ArgumentError imshow256(sprand(5,5,.5))
    end
    @testset "rgb line" begin
        io = IOBuffer()
        ensurecolor(imshow256, io, rgb_line)
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))
        @test_reference "reference/rgbline_big_imshow256.txt" res
        io = IOBuffer()
        ensurecolor(imshow256, io, rgb_line, (1, 45))
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))
        @test_reference "reference/rgbline_small1_imshow256.txt" res
        io = IOBuffer()
        ensurecolor(imshow256, io, rgb_line, (1, 19))
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))
        @test_reference "reference/rgbline_small2_imshow256.txt" res
    end
    @testset "lena" begin
        img = imresize(lena, 10, 10)
        io = IOBuffer()
        ensurecolor(imshow256, io, img)
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))
        @test_reference "reference/lena_big_imshow256.txt" res
        io = IOBuffer()
        ensurecolor(imshow256, io, img, (10, 20))
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))
        @test_reference "reference/lena_small_imshow256.txt" res
    end
    @testset "ndarray" begin
        img = rgb_line_4d
        io = IOBuffer()
        ensurecolor(imshow256, io, img)
        res = readlines(seek(io,0))
        @test_reference "reference/ndarray_imshow256.txt" res
    end
end

@testset "imshow24bit" begin
    @testset "non colorant" begin
        @test_throws ArgumentError imshow24bit(rand(5,5))
        @test_throws ArgumentError imshow24bit(sprand(5,5,.5))
    end
    @testset "rgb line" begin
        io = IOBuffer()
        ensurecolor(imshow24bit, io, rgb_line)
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))
        @test_reference "reference/rgbline_big_imshow24bit.txt" res
        io = IOBuffer()
        ensurecolor(imshow24bit, io, rgb_line, (1, 45))
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))
        @test_reference "reference/rgbline_small1_imshow24bit.txt" res
        io = IOBuffer()
        ensurecolor(imshow24bit, io, rgb_line, (1, 19))
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))
        @test_reference "reference/rgbline_small2_imshow24bit.txt" res
    end
    @testset "lena" begin
        img = imresize(lena, 10, 10)
        io = IOBuffer()
        ensurecolor(imshow24bit, io, img)
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))
        @test_reference "reference/lena_big_imshow24bit.txt" res
        io = IOBuffer()
        ensurecolor(imshow24bit, io, img, (10, 20))
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))
        @test_reference "reference/lena_small_imshow24bit.txt" res
    end
    @testset "ndarray" begin
        img = rgb_line_4d
        io = IOBuffer()
        ensurecolor(imshow24bit, io, img)
        res = readlines(seek(io,0))
        @test_reference "reference/ndarray_imshow24bit.txt" res
    end
end

@testset "imshow non1" begin
    @testset "lena" begin
        img = OffsetArray(imresize(lena, 10, 10), (-10,5))
        io = IOBuffer()
        ensurecolor(imshow256, io, img)
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))
        @test_reference "reference/lena_big_imshow.txt" res
    end
    @testset "rotation" begin
        tfm = recenter(RotMatrix(-pi/4), center(lighthouse))
        lhr = ImageTransformations.warp(lighthouse, tfm)
        io = IOBuffer()
        ensurecolor(imshow256, io, lhr)
        res = replace.(readlines(seek(io,0)), Ref("\n" => ""))
        @test_reference "reference/lighthouse_rotated.txt" res
    end
end


text:
using View5D, Test, JavaCall
# JavaCall.addOpts("-Djava.awt.headless=true") # This is needed for Linux and Mac, but the Applet.init() causes an exception
# JavaCall.init()

data1 = rand(5,5,3,2,2);
data2 = rand(5,5,3,2,2);
data3 = rand(5,5,3,4,2); # more elements

@testset "start viewers" begin
    v = @vv data1; # start the viewer
    @ve data2; # append along element
    @vt data3; # append along time
    @vt data3; # append along time
    @vr data3; # replace first
    w = @vv data2
    @vr v data3; # replace first
    @test 6 == get_num_times();
    @test 4 == get_num_elements();
    sleep(0.1); hide_viewer() ;
end

@testset "interaction with markers" begin
    @vv data1; # start the viewer
    markers = empty_marker_list(2,2)
    markers[1][3]=1.0
    markers[2][4]=2.0
    markers[2][5]=1.0
    markers[3][3]=3.0
    markers[3][5]=2.0
    markers[4][3]=1.5
    markers[4][4]=1.5
    markers[4][5]=2.0
    sleep(0.1); import_marker_lists(markers)
    sleep(0.1); exported = export_marker_lists()
    mydiff = exported .- markers
    for d in 1:4
        sleep(0.1);
        @test mydiff[d][2:9] == zeros(8)
    end
    sleep(0.1);
    @test nothing != export_markers_string()
    sleep(0.1);
    delete_all_marker_lists()
end

data1 = rand(5,5,3,1,1) .+ 1im.*rand(5,5,3,1,1);
data2 = rand(5,5,3,1,1) .+ 1im.*rand(5,5,3,1,1);
data3 = rand(5,5,3,2,1) .+ 1im.*rand(5,5,3,2,1); # more elements

@testset "complex-valued display" begin
    @test nothing != @vp data1; # start a new viewer in phase mode
    @test nothing == set_gamma(1.0);
    @test nothing != @vep data2; # start a new viewer in phase mode
    @vp data3; # new phase viewer
    @vrp data3; # replace phase data
end

@testset "view5d datatypes" begin
    @test nothing != view5d(rand(Float32,2,2,2,2,2)) # -> float
    @test nothing != view5d(rand(Float64,2,2,2,2,2)) # -> double
    @test nothing != view5d(rand(ComplexF32,2,2,2,2,2)) # > Complex
    @test nothing != view5d(rand(ComplexF64,2,2,2,2,2)) # -> Complex
    @test nothing != view5d(rand(Int8,2,2,2,2,2)) # -> Byte
    @test nothing != view5d(rand(UInt8,2,2,2,2,2)) # -> Byte
    @test nothing != view5d(rand(Int16,2,2,2,2,2)) # -> Short (signed)
    @test nothing != view5d(rand(UInt16,2,2,2,2,2)) # -> Unsigned Short
    @test nothing != view5d(rand(Int32,2,2,2,2,2)) # -> int
    @test nothing != view5d(rand(UInt32,2,2,2,2,2)) # -> double
    @test nothing != view5d(rand(Int64,2,2,2,2,2)) # -> int
    @test nothing != view5d(rand(UInt64,2,2,2,2,2)) # -> int
    @test nothing == close_all()
    v = view5d(rand(Float32,2,2,2,2,2))
    @test nothing == set_gamma(1.0,v, element=0)
    @test nothing == set_display_size(400,500,v)
    @test nothing == set_fontsize(25,v)
    @test nothing == hide_viewer(v)
    @test nothing == to_front(v)
    @test nothing == set_value_unit("mm", v,element=1)
    @test nothing == set_min_max_thresh(0.0, 3.0, v,element=1)
    @test nothing == process_key_main_window('C', v)
    @test nothing == process_key_element_window('C', v)
    @test nothing == update_panels()
    @test nothing == set_time(0,v)
    @test nothing == set_element(0,v)
    @test nothing == set_times_linked(false,v)
    @test nothing == set_elements_linked(true,v)
    @test nothing == update_panels()
    @test nothing == close_all()
end


text:
"""
This is the replacement for Julia's `@threads` macro proposed in
https://github.com/JuliaLang/julia/pull/35003
"""
macro threads(args...)
    na = length(args)
    if na != 1
        throw(ArgumentError("wrong number of arguments in @threads"))
    end
    ex = args[1]
    if !isa(ex, Expr)
        throw(ArgumentError("need an expression argument to @threads"))
    end
    if ex.head === :for
        if ex.args[1] isa Expr && ex.args[1].head === :(=)
            return _threadsfor(ex.args[1], ex.args[2])
        else
            throw(ArgumentError("nested outer loops are not currently supported by @threads"))
        end
    else
        throw(ArgumentError("unrecognized argument to @threads"))
    end
end

function _threadsfor(iter_stmt, lbody)
    loopvar   = iter_stmt.args[1]
    iter      = iter_stmt.args[2]
    rng = gensym(:rng)
    out = quote
        Base.@sync for $rng in $(Iterators.partition)($iter, $(length)($iter) ÷ $(nthreads)())
            Base.Threads.@spawn begin
                Base.@sync for $loopvar in $rng
                    $lbody
                end
            end
        end
    end
    esc(out)
end


text:
mutable struct UnitsStruct

    style::String;
    boltz::Float64;
    hplanck::Float64;
    mvv2e::Float64;
    ftm2v::Float64;
    mv2d::Float64;
    nktv2p::Float64;
    qqr2e::Float64;
    qe2f::Float64;
    vxmu2f::Float64;
    xxt2kmu::Float64;
    e_mass::Float64;    # not yet set
    hhmrr2e::Float64;
    mvh2r::Float64;
    angstrom::Float64;
    femtosecond::Float64;
    qelectron::Float64;
    dt::Float64;
    skin::Float64;

    UnitsStruct() = new();
end

function setunits(style)

style = lowercase(style)
units = UnitsStruct();
units.style = style

if (style == "lj")  # lj
    units.boltz = 1.0;
    units.hplanck = 1.0;
    units.mvv2e = 1.0;
    units.ftm2v = 1.0;
    units.mv2d = 1.0;
    units.nktv2p = 1.0;
    units.qqr2e = 1.0;
    units.qe2f = 1.0;
    units.vxmu2f = 1.0;
    units.xxt2kmu = 1.0;
    units.e_mass = 0.0;    # not yet set
    units.hhmrr2e = 0.0;
    units.mvh2r = 0.0;
    units.angstrom = 1.0;
    units.femtosecond = 1.0;
    units.qelectron = 1.0;
    units.dt = 0.005;
    units.skin = 0.3;

elseif (style == "real")  # real
    units.boltz = 0.0019872067;
    units.hplanck = 95.306976368;
    units.mvv2e = 48.88821291 * 48.88821291;
    units.ftm2v = 1.0 / 48.88821291 / 48.88821291;
    units.mv2d = 1.0 / 0.602214129;
    units.nktv2p = 68568.415;
    units.qqr2e = 332.06371;     # see also units.qqr2d_lammps_real
    units.qe2f = 23.060549;
    units.vxmu2f = 1.4393264316e4;
    units.xxt2kmu = 0.1;
    units.e_mass = 1.0/1836.1527556560675;
    units.hhmrr2e = 0.0957018663603261;
    units.mvh2r = 1.5339009481951;
    units.angstrom = 1.0;
    units.femtosecond = 1.0;
    units.qelectron = 1.0;

    units.dt = 1.0;
    units.skin = 2.0;

elseif (style == "metal")  # metal
    units.boltz = 8.617343e-5;
    units.hplanck = 4.135667403e-3;
    units.mvv2e = 1.0364269e-4;
    units.ftm2v = 1.0 / 1.0364269e-4;
    units.mv2d = 1.0 / 0.602214129;
    units.nktv2p = 1.6021765e6;
    units.qqr2e = 14.399645;
    units.qe2f = 1.0;
    units.vxmu2f = 0.6241509647;
    units.xxt2kmu = 1.0e-4;
    units.e_mass = 0.0;    # not yet set
    units.hhmrr2e = 0.0;
    units.mvh2r = 0.0;
    units.angstrom = 1.0;
    units.femtosecond = 1.0e-3;
    units.qelectron = 1.0;

    units.dt = 0.001;
    units.skin = 2.0;

elseif (style == "si")  # si
    units.boltz = 1.3806504e-23;
    units.hplanck = 6.62606896e-34;
    units.mvv2e = 1.0;
    units.ftm2v = 1.0;
    units.mv2d = 1.0;
    units.nktv2p = 1.0;
    units.qqr2e = 8.9876e9;
    units.qe2f = 1.0;
    units.vxmu2f = 1.0;
    units.xxt2kmu = 1.0;
    units.e_mass = 0.0;    # not yet set
    units.hhmrr2e = 0.0;
    units.mvh2r = 0.0;
    units.angstrom = 1.0e-10;
    units.femtosecond = 1.0e-15;
    units.qelectron = 1.6021765e-19;
    units.dt = 1.0e-8;
    units.skin = 0.001;
elseif (style == "cgs")  # cgs
    units.boltz = 1.3806504e-16;
    units.hplanck = 6.62606896e-27;
    units.mvv2e = 1.0;
    units.ftm2v = 1.0;
    units.mv2d = 1.0;
    units.nktv2p = 1.0;
    units.qqr2e = 1.0;
    units.qe2f = 1.0;
    units.vxmu2f = 1.0;
    units.xxt2kmu = 1.0;
    units.e_mass = 0.0;    # not yet set
    units.hhmrr2e = 0.0;
    units.mvh2r = 0.0;
    units.angstrom = 1.0e-8;
    units.femtosecond = 1.0e-15;
    units.qelectron = 4.8032044e-10;
    units.dt = 1.0e-8;
    units.skin = 0.1;

elseif (style == "electron")  # electron
    units.boltz = 3.16681534e-6;
    units.hplanck = 0.1519829846;
    units.mvv2e = 1.06657236;
    units.ftm2v = 0.937582899;
    units.mv2d = 1.0;
    units.nktv2p = 2.94210108e13;
    units.qqr2e = 1.0;
    units.qe2f = 1.94469051e-10;
    units.vxmu2f = 3.39893149e1;
    units.xxt2kmu = 3.13796367e-2;
    units.e_mass = 0.0;    # not yet set
    units.hhmrr2e = 0.0;
    units.mvh2r = 0.0;
    units.angstrom = 1.88972612;
    units.femtosecond = 1.0;
    units.qelectron = 1.0;

    units.dt = 0.001;
    units.skin = 2.0;

elseif (style == "micro")  # micro
    units.boltz = 1.3806504e-8;
    units.hplanck = 6.62606896e-13;
    units.mvv2e = 1.0;
    units.ftm2v = 1.0;
    units.mv2d = 1.0;
    units.nktv2p = 1.0;
    units.qqr2e = 8.987556e6;
    units.qe2f = 1.0;
    units.vxmu2f = 1.0;
    units.xxt2kmu = 1.0;
    units.e_mass = 0.0;    # not yet set
    units.hhmrr2e = 0.0;
    units.mvh2r = 0.0;
    units.angstrom = 1.0e-4;
    units.femtosecond = 1.0e-9;
    units.qelectron = 1.6021765e-7;
    units.dt = 2.0;
    units.skin = 0.1;

elseif (style == "nano")  # nano
    units.boltz = 0.013806504;
    units.hplanck = 6.62606896e-4;
    units.mvv2e = 1.0;
    units.ftm2v = 1.0;
    units.mv2d = 1.0;
    units.nktv2p = 1.0;
    units.qqr2e = 230.7078669;
    units.qe2f = 1.0;
    units.vxmu2f = 1.0;
    units.xxt2kmu = 1.0;
    units.e_mass = 0.0;    # not yet set
    units.hhmrr2e = 0.0;
    units.mvh2r = 0.0;
    units.angstrom = 1.0e-1;
    units.femtosecond = 1.0e-6;
    units.qelectron = 1.0;
    units.dt = 0.00045;
    units.skin = 0.1;
else
    error("Illegal units command");

end

return units

end

function getunits(units)
    boltz = units.boltz
    hplanck = units.hplanck
    mvv2e = units.mvv2e
    ftm2v = units.ftm2v
    mv2d = units.mv2d
    nktv2p = units.nktv2p
    qqr2e = units.qqr2e
    qe2f = units.qe2f
    vxmu2f = units.vxmu2f
    xxt2kmu = units.xxt2kmu
    e_mass = units.e_mass    # not yet set
    hhmrr2e = units.hhmrr2e
    mvh2r = units.mvh2r
    angstrom = units.angstrom
    femtosecond = units.femtosecond
    qelectron = units.qelectron
    dt = units.dt
    skin = units.skin
end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
#
# Implementation of LDpred Gibbs sampler
#

#
# true effects μ distributed as
#
#        N(0, Nh^2/Mp)   with probability p
# μ_i ~
#        0 with probability 1-p
#
# define σ2 := Nh^2/M
# Nh²/M can be obtained from LD-score regression
#
# σ2_μ = σ2/p

#
# observed (marginal) z-scores distributed as
#
# z_i|μ_i = N(μ_i, 1)
#
# σ2_z = 1 + σ2_μ
# σ_z = sqrt(1 + σ2_μ)

#
# residualized marginal effect
# i.e. correct SNP for effect of all other SNPs in LD
#
# z̃ = ẑ - (Dμ - μ)
#

#
# probability that variant j is causal
# is probability that z̃_i ~ N(0, 1 + Nh^2/Mp) vs z̃_i ~  N(0, 1)
#
#
# p_j = 1 / (1 + (1-p)/p * sqrt(1+σ2_μ) * exp(-0.5 * z̃^2 * σ2_μ / (1 + σ2_μ)))
#

#
# Posterior for μ_j is then
#
#             N(σ2_μ/(1 + σ2_μ) z̃_j, σ2_μ / (1 + σ2_μ))  with probability p_j
# μ_j | z̃_j ~
#             0                                      with probability 1 - p_j
#

#
# posterior mean of μ_j is then
#
# ω_j = p_j * z̃_j * σ^2 / (1 + σ^2)
#

#
# Gibbs sampler
#
# for iter = 1:(n_burnin + n_iter)
#   for j = 1:n_variants
#      compute z̃_j
#      compute p_j
#      sample μ_j
#      compute ω_j
#      if iter > n_burnin
#        Ω_j += ω_j
#      end
#   end
# end
# μ = Ω / n_iter
#

using StatsBase

export ldpred_gibbs,
    ldpred_gibbs_EB,
    ldpred_gibbs_acc

"""
    ldpred_gibbs(z, D, p, σ2, μ0; n_burnin = 100, n_iter = 500, verbose = false)

Run LDpred Gibbs sampler.

# Arguments

- `z::Vector`: z-scores
- `D::LDMatrix`: LD matrix
- `p::Real` : proportion of variants deemed to be causal
- `σ2::Real`: ``Nh²/M`` as estimated from LDscore regression (`estimate_h2`, or `mean(z^2)/mean(lds)`)
              Note: prior variance of non-null component of ``μ``, ``σ2_μ = σ2/p``
- `μ0::Vector`: starting estimate (e.g. from infinite model)
"""
function ldpred_gibbs(z, D, p, σ2, μ0; n_burnin = 100, n_iter = 500, verbose = false)

    n_variants = length(z)
    μ = μ0[:]
    Ω = zeros(Float64, n_variants)

    σ2_μ = σ2 / p
    σ_z = sqrt(1.0 + σ2_μ)
    η = σ2_μ / (1.0 + σ2_μ)  # shrinkage factor

    for iter = 1:(n_burnin + n_iter)
        verbose && println("Iteration $iter/$(n_burnin + n_iter)")
        for j = 1:n_variants
            z̃_j = z[j] - μ'*D[:,j] + μ[j]
            p_j = 1.0 / (1.0 + (1.0-p)/p * σ_z * exp(-0.5 * z̃_j * z̃_j * η))
            μ[j] = (rand() <= p_j) ?  rand(Normal(η * z̃_j, sqrt(η))) : 0.0
            if iter > n_burnin
                Ω[j] += p_j * z̃_j * η
            end
        end
    end

    res = Ω ./ n_iter
    res[(!isfinite).(res)] .= 0.0 # DON'T DO THIS!!! THIS HIDES ERRORS
    res
end

# version that stores the simulated μ for each iteration
function ldpred_gibbs_acc(z, D, p, σ2, μ0; n_burnin = 100, n_iter = 500, verbose = false)

    n_variants = length(z)
    μ = μ0[:]
    Ω = zeros(Float64, n_variants)

    σ2_μ = σ2 / p
    σ_z = sqrt(1.0 + σ2_μ)
    η = σ2_μ / (1.0 + σ2_μ)  # shrinkage factor

    res = []

    for iter = 1:(n_burnin + n_iter)
        verbose && println("Iteration $iter/$(n_burnin + n_iter)")
        for j = 1:n_variants
            z̃_j = z[j] - μ'*D[:,j] + μ[j]
            p_j = 1.0 / (1.0 + (1.0-p)/p * σ_z * exp(-0.5 * z̃_j * z̃_j * η))
            μ[j] = (rand() <= p_j) ?  rand(Normal(η * z̃_j, sqrt(η))) : 0.0
        end

        if iter > n_burnin
            push!(res, copy(μ))
        end
    end

    res
end

#
# Empirical Bayes extension
#

# instead of using same prior for all μ_i, use prior as estimated from EB
#

function ldpred_gibbs_EB(z, D, π, σ_z, μ0, scale = 1.0; n_burnin = 100, n_iter = 500, verbose = false)

    K = length(π)
    n_variants = length(z)
    μ = μ0[:]
    Ω = zeros(Float64, n_variants)

    σ2_μ = max.(σ_z .^ 2 .- 1, 0.0)/scale
    σ_μ = sqrt.(σ2_μ)
    σ_z_scaled = sqrt.(1.0 .+ σ2_μ)
    η = σ2_μ ./ (1.0 .+ σ2_μ) # η[1,k] : shrinkage factor for component k

    for iter = 1:(n_burnin + n_iter)
        verbose && println("Iteration $iter/$(n_burnin + n_iter)")
        for j = 1:n_variants
            z̃_j = z[j] - μ'*D[:,j] + μ[j]
            p_j = estimate_component_prob([z̃_j], π, σ_z_scaled)

            k = sample(1:K, Weights(p_j[:]))
            μ[j] = rand(Normal(η[k] * z̃_j, sqrt(η[k])))
            #println("$j\t$(z[j])\t$(z̃_j)\t$(p_j)\t$(μ[j])")
            #if !isfinite(p_j[1])
            #    break
            #end
            if iter > n_burnin
                Ω[j] += sum(p_j .* η) * z̃_j
            end
        end
    end

    res = Ω ./ n_iter
    res[(!isfinite).(res)] .= 0.0 # DON'T DO THIS!!! THIS HIDES ERRORS
    res
end


text:
using SMC:
    dt_smc2_estimation,
    smc_pmcmc_proposal_logdens,
    smc_pmcmc_proposal,
    density_tempered_pmcmc
using SMC.Models: LinearGaussianInfer, rprior, dprior, loglik_fn
using Base.Iterators

y = [
    5.22896384735402,
    2.059848221254017,
    11.77696919584194,
    5.721925811242492,
    2.647894413351219,
    8.833040488128724,
    10.112474414603064,
    6.208953089526956,
    7.483218024927158,
    1.1450746725213778,
    2.3827281539488214,
    2.5154588929528154,
    7.6938119335857795,
    11.267170019825972,
    11.203900809640977,
    4.641402249798912,
    6.393560809006704,
    5.434277785253692,
    7.754107613915341,
    2.8018519282190573,
    8.844915959029999,
    6.179788030963477,
    8.986422350276268,
    8.129873884597485,
    10.464689898594969,
    7.43846896751153,
    6.195006607054435,
    7.684383628629019,
    6.432819127902091,
    9.784040586681458,
    3.9386281953008275,
    1.0524045623884053,
    5.2287370014479535,
    3.3793567409786296,
    9.109449328047496,
    2.4107833351064665,
    3.1442786924940815,
    8.694980615700167,
    9.34957815303634,
    1.9899759308134595,
    1.7676330607869204,
    6.164419273710474,
    7.378003270378701,
    3.730353213883039,
    8.15997706282842,
    3.912665621137947,
    2.8222087161027987,
    7.89108559137954,
    2.016432756238303,
    -0.7283326731242079,
    5.34191426139306,
    3.7422889485477566,
    6.879539434261266,
    1.697715484841611,
    5.7492808736295835,
    3.4683312150848638,
    7.897465524526299,
    9.358577329317843,
    9.01915967340802,
    5.098410508423983,
    10.161486865442141,
    12.905026331072838,
    11.193444598591357,
    5.284455585283422,
    5.403785641006497,
    2.984985087125228,
    -0.1108114821742765,
    3.8443055205297396,
    1.1335781598255243,
    7.313228382179326,
    1.1291030218984508,
    0.8188795644543472,
    13.044240581728545,
    4.551403331575773,
    5.182374381380807,
    4.064021944452551,
    13.385904690258663,
    2.7383489638189253,
    4.272321073129483,
    1.936449829003382,
    7.579915913135302,
    6.891186573209264,
    5.696714256363617,
    2.9401339639559136,
    4.045792031292042,
    -2.69331630340043,
    4.253305831682493,
    -6.440609695250053,
    1.2224659946193035,
    3.6232261847393623,
    5.212788618742227,
    -0.5932923334934714,
    6.036279124314067,
    6.862343503023183,
    -2.9726591794689665,
    4.290653945351619,
    -2.420184182764869,
    0.6612926367123761,
    0.6312700637171165,
    -0.041609972691161345,
    -0.24839219509891972,
];

Y = hcat(y, y, y)

m = LinearGaussianInfer(10, Y, 101)
ps = @inferred rprior(m, 10)
d1 = @inferred dprior(m, ps[1])
ds = @inferred dprior(m, ps)

l = @inferred loglik_fn(m, ps[2])
ls = @inferred broadcast(m, ps)

prior = p -> dprior(m, p)

lower_bound = (x, n) -> 0.0
upper_bound = (x, n) -> Inf
C = (2.38^2) / 2.0

rprop_pmh =
    (theta, mean, cov) -> begin
        thetavec = smc_pmcmc_proposal(
            vcat(theta...),
            mean,
            cov;
            lower_bound = lower_bound,
            upper_bound = upper_bound,
            C = C,
        )
        NamedTuple{keys(theta)}(Tuple(thetavec))
    end

dprop_pmh =
    (theta, theta_new, mean, cov) -> smc_pmcmc_proposal_logdens(
        vcat(theta...),
        vcat(theta_new...),
        mean,
        cov;
        lower_bound = lower_bound,
        upper_bound = upper_bound,
        C = C,
    )

# Currently type-unstable because of pmap
res = dt_smc2_estimation(ps, m, prior; rprop_pmh = rprop_pmh, dprop_pmh = dprop_pmh)
@test typeof(res[1]) == typeof(ps)
@test size(res[1]) == size(ps)
@test res[2] isa Vector{Float64}
@test size(res[2]) == (size(ps, 1),)
@test res[3] isa Array{Bool,3}
@test size(res[3], 1) == size(res[3], 2) == size(ps, 1)

# Should be type-stable
ξ = 0.5
pmcmc_out = @inferred density_tempered_pmcmc(
    ps,
    ls,
    m,
    prior,
    ξ,
    ξ;
    rprop_pmh = rprop_pmh,
    dprop_pmh = dprop_pmh,
)


text:
#############################################################
# Download and parse every file from the NIST Matrix Market #
#############################################################

#Convenience function to emulate the behavior of gunzip
using GZip
function gunzip(fname)
    destname, ext = splitext(fname)
    if ext != ".gz"
        error("gunzip: $fname: unknown suffix -- ignored")
    end
    open(destname, "w") do f
        GZip.open(fname) do g
            write(f, readall(g))
        end
    end
    destname
end

#Download and parse master list of matrices
if !isfile("matrices.html")
    download("math.nist.gov/MatrixMarket/matrices.html", "matrices.html")
end

matrixmarketdata = Any[]
open("matrices.html") do f
   for line in readlines(f)
       if contains(line, """<A HREF="/MatrixMarket/data/""")
           collectionname, setname, matrixname = split(split(line, '"')[2], '/')[4:6]
           matrixname = split(matrixname, '.')[1]
           push!(matrixmarketdata, (collectionname, setname, matrixname) )
       end
   end
end

#Download one matrix at random
n = rand(1:length(matrixmarketdata))
for (collectionname, setname, matrixname) in matrixmarketdata[n:n]
    fn = string(collectionname, '_', setname, '_', matrixname)
    mtxfname = string(fn, ".mtx")
    if !isfile(mtxfname)
        url = "ftp://math.nist.gov/pub/MatrixMarket2/$collectionname/$setname/$matrixname.mtx.gz"
        gzfname = string(fn, ".mtx.gz")
        try
            download(url, gzfname)
        catch
            continue
        end
        gunzip(gzfname)
    end
end


text:
struct TransposeMap{T, A<:LinearMap{T}} <: LinearMap{T}
    lmap::A
end
struct AdjointMap{T, A<:LinearMap{T}} <: LinearMap{T}
    lmap::A
end

# transposition behavior of LinearMap objects
LinearAlgebra.transpose(A::TransposeMap) = A.lmap
LinearAlgebra.adjoint(A::AdjointMap) = A.lmap

LinearAlgebra.transpose(A::LinearMap) = TransposeMap(A)
LinearAlgebra.adjoint(A::LinearMap{<:Real}) = transpose(A)
LinearAlgebra.adjoint(A::LinearMap) = AdjointMap(A)

# properties
Base.size(A::Union{TransposeMap, AdjointMap}) = reverse(size(A.lmap))
LinearAlgebra.issymmetric(A::Union{TransposeMap, AdjointMap}) = issymmetric(A.lmap)
LinearAlgebra.ishermitian(A::Union{TransposeMap, AdjointMap}) = ishermitian(A.lmap)
LinearAlgebra.isposdef(A::Union{TransposeMap, AdjointMap}) = isposdef(A.lmap)

# comparison of TransposeMap objects
Base.:(==)(A::TransposeMap, B::TransposeMap) = A.lmap == B.lmap
Base.:(==)(A::AdjointMap, B::AdjointMap)     = A.lmap == B.lmap
Base.:(==)(A::TransposeMap, B::AdjointMap)   = false # isreal(B) && A.lmap == B.lmap # isreal(::AdjointMap) == false
Base.:(==)(A::AdjointMap, B::TransposeMap)   = false # isreal(A) && A.lmap == B.lmap # isreal(::AdjointMap) == false
Base.:(==)(A::TransposeMap, B::LinearMap)    = issymmetric(B) && A.lmap == B
Base.:(==)(A::AdjointMap, B::LinearMap)      = ishermitian(B) && A.lmap == B
Base.:(==)(A::LinearMap, B::TransposeMap)    = issymmetric(A) && B.lmap == A
Base.:(==)(A::LinearMap, B::AdjointMap)      = ishermitian(A) && B.lmap == A

# multiplication with vector
A_mul_B!(y::AbstractVector, A::TransposeMap, x::AbstractVector) =
    issymmetric(A.lmap) ? A_mul_B!(y, A.lmap, x) : At_mul_B!(y, A.lmap, x)

At_mul_B!(y::AbstractVector, A::TransposeMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x)

Ac_mul_B!(y::AbstractVector, A::TransposeMap, x::AbstractVector) =
    isreal(A.lmap) ? A_mul_B!(y, A.lmap, x) : (A_mul_B!(y, A.lmap, conj(x)); conj!(y))

A_mul_B!(y::AbstractVector, A::AdjointMap, x::AbstractVector) =
    ishermitian(A.lmap) ? A_mul_B!(y, A.lmap, x) : Ac_mul_B!(y, A.lmap, x)

At_mul_B!(y::AbstractVector, A::AdjointMap, x::AbstractVector) =
    isreal(A.lmap) ? A_mul_B!(y, A.lmap, x) : (A_mul_B!(y, A.lmap, conj(x)); conj!(y))

Ac_mul_B!(y::AbstractVector, A::AdjointMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x)


text:
# Side
# From https://bitbucket.org/ihmcrobotics/ihmc-open-robotics-software/src/575d33e1ab064f4e5957096414376fc011a98bce/IHMCRoboticsToolkit/src/us/ihmc/robotics/robotSide/RobotSide.java?at=develop&fileviewer=file-view-default
@enum Side left right
flipsign_if_right(x::Number, side::Side) = ifelse(side == right, -x, x)
Base.rand(rng::AbstractRNG, ::Type{Side}) = ifelse(rand(rng, Bool), left, right)
Base.:(-)(side::Side) = ifelse(side == right, left, right)


text:
correct = sqrt(π) * exp(-1/4) / 2
f(x) = sin(sqrt(x)) / sqrt(x)
α = 1/2

let
    x, w = laguerre(10, α)
    Q = sum(w .* f.(x))
    @test Q ≈ correct
end

let
    x, w = laguerre(10, α, left)
    X = x[2:end]
    W = w[2:end]
    Q = w[1] + sum(W .* f.(X))
    @test Q ≈ correct
end


text:
format_querypair_stem(stem::Any) = replace(string(stem), " " => "%20")

function format_querypair_stem(stem::Tuple{T,T}) where {T <: Number}
	m, M = stem
	if (M <= m )
		throw(ArgumentError("Range queries must be formatted as min,max"))
	end
	return replace(string(m)*","*string(M), " " => "%20")
end

function pairs_to_querystring(query::Pair...)
	if length(query) == 0
		return ""
	else
		# We start from an empty query string
		querystring = ""
		for (i, pair) in enumerate(query)
			# We build it pairwise for every
			delim = i == 1 ? "" : "&" # We use & to delimitate the queries
			# Let's be extra cautious and encode the spaces correctly
			root = replace(string(pair.first), " " => "%20")
			stem = format_querypair_stem(pair.second)
			# Then we can graft the pair string onto the query string
			pairstring = "$(delim)$(root)=$(stem)"
			querystring *= pairstring
		end
		return querystring
	end
end

"""
**Return an interpreted occurrence given its key**

	occurrence(key::Union{String, Integer})

The key can be given as a string or as an integer.
"""
function occurrence(key::Union{String, Integer})
  occ_url = gbifurl * "occurrence/" * string(key)
  occ_key_req = HTTP.get(occ_url)
  result = JSON.parse(String(occ_key_req.body))
  return GBIFRecord(result)
end

function _internal_occurrences_getter(query::Pair...)
	validate_occurrence_query.(query)
	occ_s_url = gbifurl * "occurrence/search"
	occ_s_req = HTTP.get(occ_s_url; query=pairs_to_querystring(query...))
	if occ_s_req.status == 200
		body = JSON.parse(String(occ_s_req.body))
		of_max = body["count"] > 100000 ? 100000 : body["count"]
		this_rec = GBIFRecord.(body["results"])
		return (this_rec, body["offset"], of_max)
	end
	return (nothing, nothing, nothing)
end

"""
**Retrieve latest occurrences based on a query**

    occurrences(query::Pair...)

This function will return the latest occurrences matching the queries -- usually
20, but this is entirely determined by the server default page size. The query
parametes must be given as pairs, and are optional. Omitting the query will
return the latest recorded occurrences.

The arguments accepted as queries are the one documented on the [GBIF
API](https://www.gbif.org/developer/occurrence) website.
"""
function occurrences(query::Pair...)
	retrieved, offset, of_max = _internal_occurrences_getter(query...)
	if !isnothing(retrieved)
		store = Vector{GBIFRecord}(undef, of_max)
		store[1:length(retrieved)] = retrieved
		return GBIFRecords(
			vcat(query...),
			store
		)
	else
		@error "Nothing retrieved"
	end
end

"""
**Retrieve latest occurrences for a taxon based on a query**

    occurrences(t::GBIFTaxon, query::Pair...)

Returns occurrences that correspond to a filter, given in `q` as a dictionary.
When first called, this function will return the latest 20 hits (or whichever
default page size GBIF uses). Future occurrences can be queried with `next!` or
`complete!`.
"""
function occurrences(t::GBIFTaxon, query::Pair...)
	levels = [:kingdom, :phylum, :class, :order, :family, :genus, :species]
	level = levels[findlast(l -> getfield(t, l) !== missing, levels)]
	taxon_query = String(level)*"Key" => getfield(t, level).second
	return occurrences(taxon_query, query...)
end


text:
using MPIMeasurements
#using Unitful
#using Compat
#using HDF5

error("Do NOT call this as a script, use commands manual in console each at a time")
# Init robot
configFile="HeadScanner.toml"
scanner = MPIScanner(configFile)
robot = getRobot(scanner)

# Reference robot und move to old teaching position
initRefZYX(robot)
#moveAbs(robot,0.0Unitful.mm,0.0Unitful.mm,0.0Unitful.mm)

# check if you at the right position
currentPos = getPos(robot,Unitful.mm)
@info "" currentPos robot.defCenterPos

# use moveAbs command to navigate to desired new teaching position
moveAbs(robot,robot.defCenterPos[1]+ 1.0Unitful.mm,robot.defCenterPos[2],robot.defCenterPos[3])
# or use moveRel
# moveRel(robot,1.0u"mm,0.0Unitful.mm,0.0Unitful.mm)
#
#
#...a few times moving manually with moveAbs or moveRel...
#
# if you have moved to your final new teaching position
TeachPosition(scanner, robot, configFile)

# change and test bacground position !!!


text:
export TodNoiseProperties, build_noise_properties, tod2map_mpi, baseline2map_mpi, destripe, baselines_covmat

using LinearAlgebra

try
    import MPI
catch
end

@doc raw"""
This structure holds a number of parameters relative to the noise level and
the 1/f baselines measured or simulated for a certain polarimeter.

Field              | Type           | Meaning
:----------------- |:-------------- |:----------------------------------------------
`pol_number`       | Int            | ID number of the polarimeter
`sigma`            | Float          | sigma of white noise
`num_of_samples`   | Int            | total number of samples
`baseline_lengths` | Array{Int}     | Array containing the length of each baseline

The structure provides a constructor that accepts the following
three keyword parameters:

- `pol`: number of the polarimeter
- `rms`: the sigma level
- `baseline_lengths`: array containing the elements of each baseline

# Example

The followng code creates a `TodNoiseProperties` object referring to polarimeter #1
(int this case, STRIP31) and specifying 0.002463 as the noise level. The data sample
contains 500 samples, and it is split into 5 baselines of equal length.

```julia
prop = TodNoiseProperties(
    pol = 1,
    rms = 0.002463,
    baselines = Int[100, 100, 100, 100, 100, 100],
)
```

"""
struct TodNoiseProperties
    polarimeter::Int
    sigma::Float64
    number_of_samples::Int
    baseline_lengths::Array{Int}

    TodNoiseProperties(;
        pol::Int,
        rms::Float64,
        baselines::Array{Int}) = new(pol, rms, sum(baselines), baselines)
end

@doc raw"""
    function build_noise_properties(detector_list, rms_list, num_of_baselines, num_of_samples) -> data_properties

This function builds a list of `TodNoiseProperties` object. These are
needed for desriping and calculation of `baselines_covmat`.

It requires in input 4 arrays containing:
- the ID number of the polarimeters that the current rank will simulate
- the white noise sigma for each polarimeter that the current rank will simulate
- the number of 1/f baselines for each polarimeter that the current rank will simulate
- the total number of samples for each polarimeter that the current rank will simulate

N.B. the ID numbers, the number of 1/f baselines and the total number of samples can be obtained using the function `get_chunk_properties`.

It returns an array of `TodNoiseProperties`, of length equal to the number of polarimeters simulated by current rank.

"""
function build_noise_properties(detector_list, rms_list, num_of_baselines, num_of_samples)
    @assert length(detector_list) == length(rms_list)
    @assert length(detector_list) == length(num_of_baselines)
    @assert length(detector_list) == length(num_of_samples)

    data_properties  = Array{TodNoiseProperties}(undef, length(detector_list))
    for i in 1:length(detector_list)
        # We use baselines of equal lengths
        baseline_lengths  = repeat([Int(num_of_samples[i] / num_of_baselines[i])], num_of_baselines[i])
        data_properties[i] = TodNoiseProperties(
            pol = detector_list[i],
            rms = rms_list[detector_list[i]],
            baselines = baseline_lengths,
        )
    end
    data_properties
end

@doc raw"""
    tod2map(pix_idx, tod, num_of_pixels; comm = nothing) -> binned_map
    tod2map(pix_idx, tod, num_of_pixels, data_properties; comm = nothing) -> binned_map

This function creates a binned map from the time-ordered data kept in
the array `tod`, assuming that each sample is observing the pixel
whose index is in `pix_idx`. The parameter `num_of_pixels` contains
the number of pixels in the Healpix map, and it is used as an upper
bound for the values in `pix_idx`. The parameter `comm` must be a MPI
communicator, or `nothing` if you are not using MPI.

If `comm` is not `nothing`, the function is parallelized using MPI, and
each process computes a map from its available data.  All partial maps
are then combined together with `MPI.Allreduce`. The function returns
an array containing the binned map.

If Array of structures `TodNoiseProperties`is passed to the function, the
output binned map will be a weighted binned map. Each sample is weighted
according to the inverse white noise variance sigma^2 of the corrispondingù
polarimeter. In this way, the less noisy polarimeters will count more in the
estimation of the map.

# Requirements
- The length of the arrays `pix_idx` and `tod` must be the same

"""
function tod2map_mpi(pix_idx, tod, num_of_pixels; comm = nothing, unseen = NaN)
    @assert length(pix_idx) == length(tod)

    T = eltype(tod)
    N = eltype(pix_idx)

    partial_map = zeros(T, num_of_pixels)
    partial_hits = zeros(N, num_of_pixels)
    binned_map = zeros(T, num_of_pixels)
    hits = zeros(N, num_of_pixels)

    @inbounds for i in eachindex(pix_idx)
        partial_map[pix_idx[i]] += tod[i]
        partial_hits[pix_idx[i]] += 1
    end

    if comm != nothing
        binned_map = MPI.Allreduce(partial_map, MPI.SUM, comm)
        hits = MPI.Allreduce(partial_hits, MPI.SUM, comm)
    else
        binned_map .= partial_map
        hits .= partial_hits
    end

    @inbounds for i in eachindex(binned_map)
        if (hits[i] > 0)
            binned_map[i] = binned_map[i] / hits[i]
        else
            binned_map[i] = unseen
        end
    end

    binned_map
end

function tod2map_mpi(pix_idx, tod, num_of_pixels, data_properties; comm = nothing, unseen = NaN)
    @assert length(pix_idx) == length(tod)

    T = eltype(tod)

    partial_map = zeros(T, num_of_pixels)
    partial_hits = zeros(T, num_of_pixels)
    binned_map = zeros(T, num_of_pixels)
    hits = zeros(T, num_of_pixels)

    start_idx = 1

    for j in eachindex(data_properties)                 #loop on detectors
        end_idx = start_idx + data_properties[j].number_of_samples - 1
        for i in start_idx:end_idx             #loop on samples
            partial_map[pix_idx[i]] += tod[i] * 1 / (data_properties[j].sigma)^2
            partial_hits[pix_idx[i]] += 1 / (data_properties[j].sigma)^2
        end
        start_idx += data_properties[j].number_of_samples
    end

    if comm != nothing
        binned_map = MPI.Allreduce(partial_map, MPI.SUM, comm)
        hits = MPI.Allreduce(partial_hits, MPI.SUM, comm)
    else
        binned_map .= partial_map
        hits .= partial_hits
    end

    @inbounds for i in eachindex(binned_map)
        if (hits[i] > 0)
            binned_map[i] = binned_map[i] / hits[i]
        else
            binned_map[i] = unseen
        end
    end

    binned_map
end

@doc raw"""
    baseline2map_mpi(pix_idx, baselines, baseline_lengths, num_of_pixels; comm = nothing)-> noise_map
    baseline2map_mpi(pix_idx, baselines, num_of_pixels, data_properties; comm = nothing) -> noise_map

This is a MPI based function: each MPI process computes a map from its
available data.  All partial maps are then combined together with
MPI.Allreduce.
The function returns an array containing the binned map.

If Array of structures `TodNoiseProperties`is passed to the function (instead of `baseline_lengths`) the output binned map will be a weighted binned map.
Each sample will be weighted according to the inverse white noise variance sigma^2 of the corrisponding polarimeter.
In this way, the less noisy polarimeters will count more in the estimation of the map.

# Requirements
- The length of `baselines` and `baseline_lengths` must be the same;
- The value `sum(baseline_lengths)` must be the same as the length of `pix_idx`.
"""
baseline2map_mpi

function baseline2map_mpi(pix_idx, baselines, baseline_lengths, num_of_pixels; comm = nothing,
                          unseen = NaN)

    T = eltype(baselines)
    N = eltype(pix_idx)

    partial_map = zeros(T, num_of_pixels)
    partial_hits = zeros(N, num_of_pixels)
    noise_map = zeros(T, num_of_pixels)
    hits = zeros(N, num_of_pixels)

    startidx = 1

    for i in eachindex(baseline_lengths)
        endidx = baseline_lengths[i] + startidx - 1

        for j in startidx:endidx
            partial_map[pix_idx[j]] += baselines[i]
            partial_hits[pix_idx[j]] += 1
        end
        startidx += baseline_lengths[i]
    end

    if comm != nothing
        noise_map = MPI.Allreduce(partial_map, MPI.SUM, comm)
        hits = MPI.Allreduce(partial_hits, MPI.SUM, comm)
    else

        noise_map .= partial_map
        hits .= partial_hits
    end

    @inbounds for i in eachindex(noise_map)
        if (hits[i] > 0)
            noise_map[i] /= hits[i]
        else
            noise_map[i] = unseen
        end
    end

    noise_map
end

function baseline2map_mpi(pix_idx, baselines, num_of_pixels, data_properties, num_of_baselines; comm = nothing, unseen = NaN)

    T = eltype(baselines)

    partial_map = zeros(T, num_of_pixels)
    partial_hits = zeros(T, num_of_pixels)
    noise_map = zeros(T, num_of_pixels)
    hits = zeros(T, num_of_pixels)

    startidx = 1
    baseline_idx = 1

    for det_idx in eachindex(data_properties)  #loop on detectors
        for sample_idx in eachindex(data_properties[det_idx].baseline_lengths)
            endidx = data_properties[det_idx].baseline_lengths[sample_idx] + startidx - 1

            @inbounds for j in startidx:endidx
                partial_map[pix_idx[j]] += baselines[baseline_idx] / (data_properties[det_idx].sigma)^2
                partial_hits[pix_idx[j]] += 1 / (data_properties[det_idx].sigma)^2
            end
            startidx += data_properties[det_idx].baseline_lengths[sample_idx]
            baseline_idx += 1
        end
    end

    if comm != nothing
        noise_map = MPI.Allreduce(partial_map, MPI.SUM, comm)
        hits = MPI.Allreduce(partial_hits, MPI.SUM, comm)
    else

        noise_map .= partial_map
        hits .= partial_hits
    end

    @inbounds for i in eachindex(noise_map)
        if (hits[i] > 0)
            noise_map[i] /= hits[i]
        else
            noise_map[i] = unseen
        end
    end

    noise_map
end

function applyz_and_sum(pix_idx, tod, num_of_pixels, data_properties, num_of_baselines; comm = nothing, unseen = NaN)

    @assert length(tod) == length(pix_idx)

    baselines_sum = zeros(eltype(tod), num_of_baselines)

    binned_map = tod2map_mpi(pix_idx,
        tod,
        num_of_pixels,
        data_properties,
        comm = comm,
        unseen = unseen)

    baseline_idx = 1
    startidx = 1
    for det_idx in eachindex(data_properties)  #loop on detectors
        for baseline_idx in eachindex(data_properties[det_idx].baseline_lengths)
            endidx = data_properties[det_idx].baseline_lengths[baseline_idx] + startidx - 1

            for j in startidx:endidx
                baselines_sum[baseline_idx] += (tod[j] - binned_map[pix_idx[j]]) * 1 / (data_properties[det_idx].sigma)^2
            end

            startidx += data_properties[det_idx].baseline_lengths[baseline_idx]
            baseline_idx += 1
        end
    end

    baselines_sum
end

function applya(baselines, pix_idx, num_of_baselines, num_of_pixels, data_properties; comm = nothing, unseen = NaN)
    @assert length(baselines) == num_of_baselines

    baselines_sum = zeros(eltype(baselines), num_of_baselines)
    total_sum = zero(eltype(baselines))

    binned_map = baseline2map_mpi(pix_idx,
        baselines,
        num_of_pixels,
        data_properties,
        num_of_baselines,
        comm = comm,
        unseen = unseen)

    startidx = 1
    baseline_idx = 1

    for det_idx in eachindex(data_properties)
        for baseline_idx in eachindex(data_properties[det_idx].baseline_lengths)
            endidx = data_properties[det_idx].baseline_lengths[baseline_idx] + startidx - 1

            for j in startidx:endidx
                baselines_sum[baseline_idx] += (baselines[baseline_idx] - binned_map[pix_idx[j]]) * 1 / (data_properties[det_idx].sigma)^2
            end

            startidx += data_properties[det_idx].baseline_lengths[baseline_idx]
            baseline_idx += 1
        end
    end

    #needed to assure that sum(baselines)==0

    if comm != nothing
        total_sum = MPI.Allreduce([sum(baselines)], MPI.SUM, comm)[1]
    else
        total_sum = sum(baselines)
    end

    baselines_sum .+= total_sum
end

function mpi_dot_prod(x, y; comm = nothing)

    @assert eltype(x) == eltype(y)
    result = zero(eltype(x))
    local_sum::eltype(x) = dot(x, y)

    if comm != nothing
        result = MPI.Allreduce([local_sum], MPI.SUM, comm)[1]
    else
        result = local_sum
    end

    return result
end

mutable struct DestripingResults
    threshold::Float64
    max_iter::Int
    convergence_param_list::Array{Float64}
    best_iteration::Int
    best_baselines::Array{Float64}
    best_sky_map::Any
    baseline_history::Any

    DestripingResults() = new(
        1e-9,
        1000,
        Float64[],
        -1,
        Float64[],
        nothing,
        nothing,
    )
end

function conj_grad(
    results::DestripingResults,
    baselines_sum,
    pix_idx,
    tod,
    num_of_pixels,
    data_properties,
    num_of_baselines,
    rank;
    comm = nothing,
    save_baseline_history = false,
    callback = nothing,
)

    T = eltype(tod)
    N = eltype(pix_idx)

    baselines = Array{T}(undef, num_of_baselines)
    r = Array{T}(undef, num_of_baselines)
    r_next = Array{T}(undef, num_of_baselines)
    p = Array{T}(undef, num_of_baselines)
    Ap = Array{T}(undef, num_of_baselines)

    baselines .= 0    #starting baselines
    convergence_parameter = zero(T)
    rdotr = zero(T)
    rdotr_next = zero(T)

    best_convergence_parameter = zero(T)
    best_baselines = zeros(T, num_of_baselines)
    results.best_iteration = 0

    r = baselines_sum - applya(baselines, pix_idx, num_of_baselines, num_of_pixels, data_properties, comm = comm)  #residual
    p .= r
    rdotr = mpi_dot_prod(r, r, comm = comm)
    best_convergence_parameter = sqrt(rdotr)

    results.baseline_history = save_baseline_history ? [] : nothing
    results.best_baselines = Array{Float64}(undef, num_of_baselines)
    results.best_baselines .= baselines

    rdotr == 0 && return results

    results.convergence_param_list = Float64[]
    iter_idx = 0
    while true
        Ap = applya(p, pix_idx,  num_of_baselines, num_of_pixels, data_properties, comm = comm)

        rdotr = mpi_dot_prod(r, r, comm = comm)
        pdotAp = mpi_dot_prod(p, Ap, comm = comm)

        alpha = rdotr / pdotAp
        @. baselines += alpha * p
        save_baseline_history && push!(results.baseline_history, copy(baselines))

        @. r_next = r - alpha * Ap
        rdotr_next = mpi_dot_prod(r_next, r_next, comm = comm)
        convergence_parameter = sqrt(rdotr_next)
        push!(results.convergence_param_list, convergence_parameter)

        if (convergence_parameter < best_convergence_parameter)
            best_convergence_parameter = convergence_parameter
            results.best_baselines .= baselines
            results.best_iteration = iter_idx
        end

        isnothing(callback) || callback(
            iter_idx = iter_idx,
            max_iter = results.max_iter,
            convergence_parameter = convergence_parameter,
            convergence_threshold = results.threshold,
        )

        ((convergence_parameter < results.threshold) || (iter_idx > results.max_iter)) && break

        beta = rdotr_next / rdotr

        @. p = r_next + beta * p
        r .= r_next
        iter_idx += 1
    end
end

function destriped_map(baselines, pix_idx, tod, data_properties, num_of_pixels, num_of_baselines; comm = nothing, unseen = NaN)
    @assert length(tod) == length(pix_idx)
    tod2map_mpi(pix_idx, tod, num_of_pixels, data_properties, comm = comm) - baseline2map_mpi(pix_idx, baselines, num_of_pixels, data_properties, num_of_baselines, comm = comm)
end

@doc raw"""
    destripe(pix_idx, tod, num_of_pixels, data_properties, rank; comm = nothing, threshold = 1e-9, max_iter = 10000, save_baseline_history = false) -> DestripingResults

This MPI based function creates a map from a TOD and removes both 1/f
and white noise, using the destriping technique.

The parameters passed to the function have the following meaning:

- `pix_idx`: array containing the indices of the pixels visited by the instrument

- `tod`: the values measured by the polarimeters for each pixel
  (either I, Q, or U)

- `num_of_pixels`: the number of pixels in the map to be
   produced. This is used as an upper limit for the values in `pix_idx`

- `data_properties`: an array of structures `TodNoiseProperties`,
   holding information on each simulated polarimeter noise level,
   number and length of 1/f baselines and total number of samples.
   It can be obtained by using function `build_noise_properties`.

- `rank`: the rank of the current MPI process

- `comm`: the MPI communicator object.

The following arguments are optional:

- `threshold` is used by the conjugated-gradient algorithm. When the
   residual error of the iteration goes below this value, the iteration
   stops. The smaller the value, the more accurate the solution.

- `max_iter` is the maximum number of iterations to be executed in the
   conjugated-gradient algorithm. If the algorithm does not converge
   within this number of iterations, the process will quit without
   having reached the convergence threshold (see the `threshold`
   keyword above).

- If `save_baseline_history` is `true`, the return value will contain the
  sequence of baselines tested by the CG algorithm. Each MPI process will
  hold its own baselines.

The function returns a `DestripingResults` object containings the destriped
map, the sequence of baselines, and other information describing the
convergence of the CG algorithm.

# Remarks
- The length of the arrays `pix_idx` and `tod` must be the same;
- If you do not specify `comm`, no MPI will be used
"""
function destripe(
    pix_idx,
    tod,
    num_of_pixels,
    data_properties,
    rank;
    comm = nothing,
    threshold = 1e-9,
    max_iter = 1_000,
    save_baseline_history = false,
    unseen = NaN,
    callback = nothing,
)

    @assert length(pix_idx) == length(tod)

    num_of_baselines = 0
    for i in 1:length(data_properties)
        num_of_baselines += length(data_properties[i].baseline_lengths)
    end

    baselines_sum = applyz_and_sum(pix_idx,
        tod,
        num_of_pixels,
        data_properties,
        num_of_baselines,
        comm = comm,
        unseen = unseen,
    )

    results = DestripingResults()
    results.threshold = threshold
    results.max_iter = max_iter

    conj_grad(
        results,
        baselines_sum,
        pix_idx,
        tod,
        num_of_pixels,
        data_properties,
        num_of_baselines,
        rank,
        comm = comm,
        save_baseline_history = save_baseline_history,
        callback = callback,
    )

    # once we have an estimate of the baselines, we can build the destriped map
    results.best_sky_map = destriped_map(results.best_baselines,
        pix_idx,
        tod,
        data_properties,
        num_of_pixels,
        num_of_baselines,
        comm = comm,
        unseen = unseen,
    )

    results
end

@doc raw"""
    baselines_covmat(polarimeters, sigma_k, baseline_length_s, fsamp_hz, total_time) -> covariance_matrix

This function produces the covariance matrix of the 1/f baselines computed by the destriper.
As an approximation, we ignore nondiagonal terms (i.e. the correlations between different baselines).
The function output is therefore an array corresponding to the covariance matrix diagonal.
Another assumption is that the white noise variance stays constant over a given baseline.

The baseline error is computed according to equation 29 in https://arxiv.org/abs/0904.3623

The parameters passed to the function have the following meaning:
    -`polarimeters`: the array of polarimeters ID numbers
    -`sigma_k`: the array of corresponding white noise sigma (in K)
    -`baseline_length_s` = the length (in s) of each 1/f baseline
    -`fsamp_hz` = the sampling frequency (in Hz)
    -`total_time` = the duration (in s) of the observation

The function firstly computes an Array of structures `TodNoiseProperties`,
which matches the white noise variance sigma with the number of 1/f baselines and their lengths, for each polarimeter.
Then it computes the covariance matrix according to these matches.
"""
function baselines_covmat(polarimeters, sigma_k, baseline_length_s, fsamp_hz, total_time)

    ALL_data_properties = Array{TodNoiseProperties}(undef, length(polarimeters))
    baselines_per_pol =  Int64(total_time / baseline_length_s)
    for i in 1:length(polarimeters)
        baseline_lengths  = repeat([baseline_length_s * fsamp_hz], baselines_per_pol)
        ALL_data_properties[i]  = TodNoiseProperties(pol = polarimeters[i],
            rms = sigma_k[i],
            baselines = baseline_lengths,
        )
    end

    covariance_matrix = []
    for i in 1:length(ALL_data_properties)
        partial_covmat = Array{Float64}(undef, length(ALL_data_properties[i].baseline_lengths))

        for j in 1:length(ALL_data_properties[i].baseline_lengths)
            sigma = ALL_data_properties[i].sigma
            this_baseline_length = ALL_data_properties[i].baseline_lengths[j]

            partial_covmat[j] = sigma^2 / this_baseline_length
        end
        covariance_matrix = append!(covariance_matrix, partial_covmat)
    end
    return covariance_matrix
end


text:
module PkgC

g = Int64(2)

greet(x::Int64) = x * 5

precompile(greet, (Int64,))

end # module


text:
import RowEchelon

#=
Homework 5
Section 2.1
Problem 1
=#

# Initialize all Matrices
A = [2 0 -1 ; 4 -5 2]
B = [7 -5 1 ; 1 -4 -3]
C = [1 2 ;-2 1]
D = [3 5 ; -1 4]
E = [-5 ; 3]

println("-2A:")
println(-2*A)

println("B-2A:")
println(B-(2*A))

# After this point, products may not exist
# Use try/catch to find nonexistent products.
println("AC:")
try
	println(A*C)
catch
	println("AC cannot be computed.")
end

println("CD:")
try
	println(C*D)
catch
	println("CD cannot be computed.")
end


text:
#=

  Zebra puzzle in Julia ConstraintSolver.jl

  Lewis Carrol's classical puzzle with five houses and a zebra:

  Five men with different nationalities live in the first five houses
  of a street.  They practise five distinct professions, and each of
  them has a favourite animal and a favourite drink, all of them
  different.  The five houses are painted in different colours.

  The Englishman lives in a red house.
  The Spaniard owns a dog.
  The Japanese is a painter.
  The Italian drinks tea.
  The Norwegian lives in the first house on the left.
  The owner of the green house drinks coffee.
  The green house is on the right of the white one.
  The sculptor breeds snails.
  The diplomat lives in the yellow house.
  Milk is drunk in the middle house.
  The Norwegian's house is next to the blue one.
  The violinist drinks fruit juice.
  The fox is in a house next to that of the doctor.
  The horse is in a house next to that of the diplomat.

  Who owns a Zebra, and who drinks water?

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function zebra(print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    n = 5
    @variable(model, 1 <= nationality[1:n] <= n, Int)
    English, Spaniard, Japanese, Italian, Norwegian = nationality
    nationality_str = ["English", "Spaniard", "Japanese", "Italian", "Norwegian"]
    @variable(model, 1 <= color[1:n] <= n, Int)
    Red, Green, White, Yellow, Blue = color
    @variable(model, 1 <= profession[1:n] <= n, Int)
    Painter, Sculptor, Diplomat, Violinist, Doctor = profession
    @variable(model, 1 <= pet[1:n] <= n, Int)
    Dog, Snails, Fox, Horse, Zebra = pet
    @variable(model, 1 <= drink[1:n] <= n, Int)
    Tea, Coffee, Milk, Juice, Water = drink

    @constraint(model, nationality in CS.AllDifferent())
    @constraint(model, color in CS.AllDifferent())
    @constraint(model, profession in CS.AllDifferent())
    @constraint(model, pet in CS.AllDifferent())
    @constraint(model, drink in CS.AllDifferent())

    @constraint(model, English == Red)
    @constraint(model, Spaniard == Dog)
    @constraint(model, Japanese == Painter)
    @constraint(model, Italian == Tea)
    @constraint(model, Norwegian == 1)
    @constraint(model, Green == Coffee)
    @constraint(model, Green == White + 1)
    @constraint(model, Sculptor == Snails)
    @constraint(model, Diplomat == Yellow)
    @constraint(model, Milk == 3)
    my_abs(model,Norwegian,Blue,1)
    @constraint(model, Violinist == Juice)
    my_abs(model,Fox,Doctor,1)
    my_abs(model, Horse,Diplomat, 1)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                nationality_val = convert.(Integer,JuMP.value.(nationality; result=sol))
                color_val = convert.(Integer,JuMP.value.(color; result=sol))
                profession_val = convert.(Integer,JuMP.value.(profession; result=sol))
                pet_val = convert.(Integer,JuMP.value.(pet; result=sol))
                drink_val = convert.(Integer,JuMP.value.(drink; result=sol))
                println("nationality:$nationality_val")
                println("color      :$color_val")
                println("profession :$profession_val")
                println("pet        :$pet_val")
                println("drink      :$drink_val")
                zebra_val = convert.(Integer,JuMP.value.(Zebra; result=sol))
                water_val = convert.(Integer,JuMP.value.(Water; result=sol))
                println("zebra:$zebra_val water:$water_val")
                z = ["$(nationality_str[i]) owns the zebra\n" for i in 1:n if nationality_val[i] == zebra_val][1]
                w = ["$(nationality_str[i]) drinks water\n" for i in 1:n if nationality_val[i] == water_val][1]
                println("$z$w")
                println()
            end
        end
    else
        println("status:$status")
    end

    return status
end

@time zebra()

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
@benchset "kernel" begin

    @benchset "ZeroKernel{Float64}()" begin
        create_benchmarks(ZeroKernel{Float64}(); grads=false)
    end
    # @benchset "OneKernel()" begin
    #     create_benchmarks(OneKernel(); grads=false)
    # end
    # @benchset "EQ" begin
    #     @benchset "Real CPU" create_benchmarks(EQ())
        # @benchset "Real GPU" create_benchmarks(
        #     EQ();
        #     x=randn(Float32), x′=randn(Float32),
        #     x̄s=[CuArray{Float32}(randn(N)) for N in Ns()],
        #     x̄′s=[CuArray{Float32}(randn(N)) for N in Ns()],
        # )

    #     for D in Ds()
    #         @benchset "ColVecs (D=$D) CPU" create_benchmarks(
    #             EQ();
    #             x=randn(D), x′=randn(D),
    #             x̄s=[ColVecs(randn(D, N)) for N in Ns()],
    #             x̄′s=[ColVecs(randn(D, N)) for N in Ns()],
    #         )
    #     end

    #     # See https://github.com/FluxML/Zygote.jl/issues/44
    #     @benchset "EQ Almost-Toeplitz" create_benchmarks(
    #         EQ();
    #         x=randn(), x′=randn(),
    #         x̄s=[range(-randn(), step=randn(), length=N) for N in Ns()],
    #         x̄′s=[range(-randn(), step=randn(), length=N) for N in Ns()],
    #     )

    #     δ = randn()
    #     @benchset "EQ Toeplitz" create_benchmarks(
    #         EQ();
    #         x=randn(), x′=randn(),
    #         x̄s=[range(-randn(), step=δ, length=N) for N in Ns()],
    #         x̄′s=[range(-randn(), step=δ, length=N) for N in Ns()],
    #     )
    # end

    # @benchset "PerEQ" begin
    #     @benchset "Real CPU" create_benchmarks(PerEQ())
    #     # @benchset "Real GPU" create_benchmarks(PerEQ();
    #     #     x = randn(Float32), x′=randn(Float32),
    #     #     x̄s=[CuArray{Float32}(randn(N)) for N in Ns()],
    #     #     x̄′s=[CuArray{Float32}(randn(N)) for N in Ns()],
    #     # )
    # end
end


text:
module PiTemperature

export
    sensorids,
    readtemp

include("tempsensor.jl")

end # module


text:
## ExoplanetsSysSim/examples/hsu_etal_2018/abc_run_quick.jl
## (c) 2018 Danley C. Hsu & Eric B. Ford
# Script for producing Q1-Q16 FGK planet candidate occurrence rate estimates
#    using smaller simulated catalogs

include(joinpath(pwd(), "abc_setup.jl"))

using SysSimABC
using ExoplanetsSysSim
using JLD
using StatsBase

out2txt = false # Write occurrence rates & densities to text files

abc_plan = setup_abc()
EvalSysSimModel.add_param_fixed(EvalSysSimModel.sim_param_closure,"num_targets_sim_pass_one",1000)  # Set universal simulated catalog size
@time output = run_abc(abc_plan)
println("")

if out2txt
    file_rate = open("rate_output.txt", "w")
    file_dens = open("dens_output.txt", "w")
end

limitP = get_any(EvalSysSimModel.sim_param_closure, "p_lim_arr", Array{Float64,1})
limitR = get_any(EvalSysSimModel.sim_param_closure, "r_lim_arr", Array{Float64,1})

weight_vec = pweights(output.weights)

for p_ind = 1:(length(limitP)-1)
    for r_ind = 1:(length(limitR)-1)
        dens_denom = 1.0/(log2(limitP[p_ind+1])-log2(limitP[p_ind]))/(log2(limitR[r_ind+1])-log2(limitR[r_ind]))

        bin_ind = (p_ind-1)*(length(limitR)-1)+r_ind
        quant_arr = quantile(output.theta[bin_ind,:], weight_vec, [0.1587, 0.5, 0.8413])

        println("-----------------------------")
        println("Orbital Period (day) = ", string(limitP[p_ind:p_ind+1]), " / Planet Radius (R_earth) = ", string(limitR[r_ind:r_ind+1]/ExoplanetsSysSim.earth_radius))
        println("")
        println("Rate = ", string(quant_arr[2], " + ", quant_arr[3]-quant_arr[2], " - ", quant_arr[2]-quant_arr[1]))
        println("Density = ", string(quant_arr[2]*dens_denom, " + ", (quant_arr[3]-quant_arr[2])*dens_denom, " - ", (quant_arr[2]-quant_arr[1])*dens_denom))

        if out2txt
            write(file_rate, string(quant_arr[2], " + ", quant_arr[3]-quant_arr[2], " - ", quant_arr[2]-quant_arr[1], "\n"))
            write(file_dens, string(quant_arr[2]*dens_denom, " + ", (quant_arr[3]-quant_arr[2])*dens_denom, " - ", (quant_arr[2]-quant_arr[1])*dens_denom, "\n"))
        end
    end
end

if out2txt
    close(file_rate)
    close(file_dens)
end

println("-----------------------------")
println("")
println(EvalSysSimModel.get_ss_obs())
save(string("test-pop-out.jld"), "output", output, "ss_true", EvalSysSimModel.get_ss_obs())


text:
using Base.Test
using MathProgBase
using MathProgBase.SolverInterface
using Xpress

function quadprogtest(solver=MathProgBase.defaultQPsolver)
    println("Testing quadprog with solver ", string(typeof(solver)))

    sol = quadprog([0., 0., 0.],[2. 1. 0.; 1. 2. 1.; 0. 1. 2.],[1. 2. 3.; 1. 1. 0.],'>',[4., 1.],-Inf,Inf,solver)

    @test sol.status == :Optimal
    @test_approx_eq_eps sol.objval 130/70 1e-5
    @test_approx_eq_eps norm(sol.sol[1:3] - [0.5714285714285715,0.4285714285714285,0.8571428571428572]) 0.0 1e-4
    println("prs 1")
    let
    m = LinearQuadraticModel(solver)
    loadproblem!(m, [1. 2. 3.; 1. 1. 0.],[-Inf,-Inf,-Inf],[Inf,Inf,Inf],[0.,0.,0.],[4., 1.],[Inf,Inf], :Min)

    setquadobj!(m,diagm([10.0,10.0,10.0]))
    rows = [1, 2, 2, 2, 3, 3, 3]
    cols = [1, 1, 1, 2, 2, 3, 3]
    vals = Float64[2, 0.5, 0.5, 2, 1, 1, 1]
    setquadobj!(m,rows,cols,vals)
    optimize!(m)
    stat = status(m)
    @test stat == :Optimal
    @test_approx_eq_eps getobjval(m) 130/70 1e-5
    @test_approx_eq_eps norm(getsolution(m) - [0.5714285714285715,0.4285714285714285,0.8571428571428572]) 0.0 1e-4
    end
    println("prb 2")

    let
    m = LinearQuadraticModel(solver)
    loadproblem!(m, [-1. 1.; 1. 1.], [0.,0.], [Inf,Inf], [1.,1.], [0.,0.], [Inf,Inf], :Max)
    addquadconstr!(m, [2], [1.], [1], [1], [1.], '<', 2)
    optimize!(m)
    stat = status(m)
    @test stat == :Optimal
    @test_approx_eq_eps getobjval(m) 2.25 1e-4
    @test_approx_eq_eps norm(getsolution(m) - [0.5,1.75]) 0.0 1e-3
    end
    println("Done")
end
quadprogtest(XpressSolver(OPTIMALITYTOL = 1e-8))


text:
const D3=3
const POINT=15
const UNSET = 0

function GmshDiscreteModel(mshfile; renumber=true)
  @check_if_loaded
  if !isfile(mshfile)
    error("Msh file not found: $mshfile")
  end

  gmsh.initialize()
  gmsh.option.setNumber("General.Terminal", 1)
  gmsh.option.setNumber("Mesh.SaveAll", 1)
  gmsh.option.setNumber("Mesh.MedImportGroupsOfNodes", 1)
  gmsh.open(mshfile)

  renumber && gmsh.model.mesh.renumberNodes()
  renumber && gmsh.model.mesh.renumberElements()

  grid, cell_to_entity = _setup_grid(gmsh)
  grid_topology = UnstructuredGridTopology(grid)
  labeling = _setup_labeling(gmsh,grid,grid_topology,cell_to_entity)

  gmsh.finalize()

  UnstructuredDiscreteModel(grid,grid_topology,labeling)
end

function _setup_grid(gmsh)

  D = _setup_dim(gmsh)

  node_to_coords = _setup_node_coords(gmsh,D)

  cell_to_nodes, nminD = _setup_connectivity(gmsh,D)

  cell_to_type, reffes, orientation = _setup_reffes(gmsh,D)

  cell_to_entity = _setup_cell_to_entity(
    gmsh,D,length(cell_to_nodes),nminD)

  grid = UnstructuredGrid(
    node_to_coords,
    cell_to_nodes,
    reffes,
    cell_to_type,
    orientation)

  (grid, cell_to_entity)

end

function _setup_dim(gmsh)
  entities = gmsh.model.getEntities()
  D = -1
  for e in entities
    D = max(D,e[1])
  end
  if D == -1
    gmsh.finalize()
    error("No entities in the msh file.")
  end
  D
end

function _setup_node_coords(gmsh,D)
  nodeTags, coord, parametricCoord = gmsh.model.mesh.getNodes()
  nmin = minimum(nodeTags)
  nmax = maximum(nodeTags)
  nnodes = length(nodeTags)
  if !(nmax == nnodes && nmin == 1)
    gmsh.finalize()
    error("Only consecutive node tags allowed.")
  end
  node_to_coords = zeros(Point{D,Float64},nnodes)
  _fill_node_coords!(node_to_coords,nodeTags,coord,D)
  node_to_coords
end

function _fill_node_coords!(node_to_coords,nodeTags,coord,D)
  m = zero(mutable(Point{D,Float64}))
  for node in nodeTags
    for j in 1:D
      k = (node-1)*D3 + j
      xj = coord[k]
      m[j] = xj
    end
    node_to_coords[node] = m
  end
end

function _setup_connectivity(gmsh,d)

  elemTypes, elemTags, nodeTags = gmsh.model.mesh.getElements(d)

  if length(elemTypes) == 0
    ncells = 0
    ndata = 0
    nmin = 1
    cell_to_nodes_prts = zeros(Int,ncells+1)
    cell_to_nodes_data = zeros(Int32,ndata)
    cell_to_nodes = Table(cell_to_nodes_data,cell_to_nodes_prts)
    return (cell_to_nodes, nmin)
  end

  ncells, nmin, nmax = _check_cell_tags(elemTags)

  etype_to_nlnodes = _setup_etype_to_nlnodes(elemTypes,gmsh)

  ndata = sum([length(t) for t in nodeTags])

  cell_to_nodes_data = zeros(Int,ndata)
  cell_to_nodes_prts = zeros(Int32,ncells+1)

  _fill_connectivity!(
    cell_to_nodes_data,
    cell_to_nodes_prts,
    nmin-1,
    etype_to_nlnodes,
    elemTypes,
    elemTags,
    nodeTags,
    d)

  cell_to_nodes = Table(cell_to_nodes_data,cell_to_nodes_prts)

  (cell_to_nodes, nmin)

end

function _check_cell_tags(elemTags)
  nmin::Int = minimum([minimum(t) for t in elemTags])
  nmax::Int = maximum([maximum(t) for t in elemTags])
  ncells = sum([length(t) for t in elemTags])
  if !( (nmax-nmin+1) == ncells)
    gmsh.finalize()
    error("Only consecutive elem tags allowed.")
  end
  (ncells,nmin,nmax)
end

function _setup_etype_to_nlnodes(elemTypes,gmsh)
  netypes = maximum(elemTypes)
  etype_to_nlnodes = zeros(Int,netypes)
  for etype in elemTypes
    name, dim, order, numv, parv =
      gmsh.model.mesh.getElementProperties(etype)
    etype_to_nlnodes[etype] = numv
  end
  etype_to_nlnodes
end

function  _fill_connectivity!(
    cell_to_nodes_data,
    cell_to_nodes_prts,
    o,
    etype_to_nlnodes,
    elemTypes,
    elemTags,
    nodeTags,
    d)

  for (j,etype) in enumerate(elemTypes)
    nlnodes = etype_to_nlnodes[etype]
    i_to_cell = elemTags[j]
    for cell in i_to_cell
      cell_to_nodes_prts[cell+1-o] = nlnodes
    end
  end

  length_to_ptrs!(cell_to_nodes_prts)

  for (j,etype) in enumerate(elemTypes)
    nlnodes = etype_to_nlnodes[etype]
    i_to_cell = elemTags[j]
    i_lnode_to_node = nodeTags[j]
    if (nlnodes == d+1)
      _orient_simplex_connectivities!(nlnodes,i_lnode_to_node)
    elseif (nlnodes == 4)
      _sort_quad_connectivites!(nlnodes,i_lnode_to_node)
    elseif (nlnodes == 8)
      _sort_hex_connectivites!(nlnodes,i_lnode_to_node)
    end
    for (i,cell) in enumerate(i_to_cell)
      a = cell_to_nodes_prts[cell-o]-1
      for lnode in 1:nlnodes
        node = i_lnode_to_node[(i-1)*nlnodes+lnode]
        cell_to_nodes_data[a+lnode] = node
      end
    end
  end

end

function _orient_simplex_connectivities!(nlnodes,i_lnode_to_node)
  aux = zeros(eltype(i_lnode_to_node),nlnodes)
  offset = nlnodes-1
  for i in 1:nlnodes:length(i_lnode_to_node)
    aux = i_lnode_to_node[i:i+offset]
    sort!(aux)
    i_lnode_to_node[i:i+offset] = aux
  end
end

function _sort_quad_connectivites!(nlnodes,i_lnode_to_node)
  aux = zero(eltype(i_lnode_to_node))
  offset = nlnodes-1
  for i in 1:nlnodes:length(i_lnode_to_node)
    aux = i_lnode_to_node[i+offset-1]
    i_lnode_to_node[i+offset-1] = i_lnode_to_node[i+offset]
    i_lnode_to_node[i+offset] = aux
  end
end

function _sort_hex_connectivites!(nlnodes,i_lnode_to_node)
  perm = [1, 2, 4, 3, 5, 6, 8, 7]
  aux = zeros(eltype(i_lnode_to_node),nlnodes)
  offset = nlnodes-1
  for i in 1:nlnodes:length(i_lnode_to_node)
    aux = i_lnode_to_node[i:i+offset]
    i_lnode_to_node[i:i+offset] = aux[perm]
  end
end

function _setup_reffes(gmsh,d)

  elemTypes, elemTags, nodeTags = gmsh.model.mesh.getElements(d)

  if !(length(elemTypes)==1)
    gmsh.finalize()
    s = """
    Only one element type per dimension allowed for the moment.
    Dimension $d has $(length(elemTypes)) different element types
    """
    error(s)
  end

  ncells, nmin, nmax = _check_cell_tags(elemTags)
  cell_to_type = fill(Int8(1),ncells)

  etype::Int = first(elemTypes)
  name, dim, order::Int, numv, parv = gmsh.model.mesh.getElementProperties(etype)

  if order == 0 && etype == POINT
    order = 1
  end

  if order != 1
    gmsh.finalize()
    error("For the moment only for first-order elements")
  end

  reffe = _reffe_from_etype(etype)
  reffes = [reffe,]

  orientation = Val(is_simplex(get_polytope(reffe)))

  (cell_to_type, reffes, orientation)
end

function _reffe_from_etype(eltype)
  if eltype == 1
    SEG2
  elseif eltype == 2
    TRI3
  elseif eltype == 3
    QUAD4
  elseif eltype == 4
    TET4
  elseif eltype == 5
    HEX8
  elseif eltype == POINT
    VERTEX1
  else
    gmsh.finalize()
    error("Unsupported element type. elemType: $etype")
  end
end

function _setup_cell_to_entity(gmsh,d,ncells,nmin)
  cell_to_entity = fill(UNSET,ncells)
  entities = gmsh.model.getEntities(d)
  for e in entities
    _, elemTags, _ = gmsh.model.mesh.getElements(e[1], e[2])
    _fill_cell_to_entity!(cell_to_entity,elemTags,e[2],nmin-1)
  end
  cell_to_entity
end

function _fill_cell_to_entity!(cell_to_entity,elemTags,entity,o)
  for i_to_cell in elemTags
    for cell in i_to_cell
      cell_to_entity[cell-o] = entity
    end
  end
end

function _setup_labeling(gmsh,grid,grid_topology,cell_to_entity)

  D = num_cell_dims(grid)
  dim_to_gface_to_nodes, dim_gface_to_entity = _setup_faces(gmsh,D)
  push!(dim_to_gface_to_nodes,get_cell_nodes(grid))
  push!(dim_gface_to_entity,cell_to_entity)

  dim_to_group_to_entities = _setup_dim_to_group_to_entities(gmsh)
  dim_to_group_to_name = _setup_dim_to_group_to_name(gmsh)
  dim_to_offset = _setup_dim_to_offset(gmsh)

  nnodes = num_nodes(grid)
  node_to_label = _setup_node_to_label(gmsh,dim_to_offset,nnodes)

  labeling = _setup_face_labels(
    grid_topology,
    dim_to_gface_to_nodes,
    dim_gface_to_entity,
    dim_to_offset,
    dim_to_group_to_entities,
    dim_to_group_to_name,
    node_to_label)

  labeling

end

function _setup_faces(gmsh,D)
  dim_to_gface_to_nodes = []
  dim_gface_to_entity = []
  for d in 0:(D-1)
    face_to_nodes, nmin  = _setup_connectivity(gmsh,d)
    face_to_entity = _setup_cell_to_entity(gmsh,d,length(face_to_nodes),nmin)
    push!(dim_to_gface_to_nodes,face_to_nodes)
    push!(dim_gface_to_entity,face_to_entity)
  end
  (dim_to_gface_to_nodes, dim_gface_to_entity)
end

function _setup_dim_to_group_to_entities(gmsh)
  dim_to_group_to_entities = Vector{Vector{Int}}[]
  for d in 0:D3
    pgs = gmsh.model.getPhysicalGroups(d)
    n = length(pgs)
    group_to_entities = Vector{Int}[]
    for pg in pgs
      es = gmsh.model.getEntitiesForPhysicalGroup(pg...)
      entities = Int[]
      for e in es
        push!(entities,e)
      end
      push!(group_to_entities,entities)
    end
    push!(dim_to_group_to_entities,group_to_entities)
  end
  dim_to_group_to_entities
end

function _setup_dim_to_group_to_name(gmsh)
  u = 1
  dim_to_group_to_name = Vector{String}[]
  for d in 0:D3
    pgs = gmsh.model.getPhysicalGroups(d)
    n = length(pgs)
    names = String[]
    for pg in pgs
      _name = gmsh.model.getPhysicalName(pg...)
      if _name == ""
        name = "untitled$u"
        u += 1
      else
        name = _name
      end
      push!(names,name)
    end
    push!(dim_to_group_to_name,names)
  end
  dim_to_group_to_name
end

function _setup_dim_to_offset(gmsh)
  entities = gmsh.model.getEntities()
  dim_to_nentities = zeros(Int,D3+1)
  for e in entities
    d = e[1]
    id = e[2]
    _id = dim_to_nentities[d+1]
    dim_to_nentities[d+1] = max(id,_id)
  end
  dim_to_offset = zeros(Int,D3+1)
  for d=1:D3
    dim_to_offset[d+1] += dim_to_nentities[d-1+1] + dim_to_offset[d-1+1]
  end
  dim_to_offset
end

function _setup_node_to_label(gmsh,dim_to_offset,nnodes)
  node_to_label = zeros(Int,nnodes)
  entities = gmsh.model.getEntities()
  for (dim, entity) in entities
    nodes, _, _ = gmsh.model.mesh.getNodes(dim, entity)
    node_to_label[nodes] .= (entity + dim_to_offset[dim+1])
  end
  node_to_label
end

function _setup_face_labels(
  grid_topology,
  dim_to_gface_to_nodes,
  dim_gface_to_entity,
  dim_to_offset,
  dim_to_group_to_entities,
  dim_to_group_to_name,
  node_to_label)

  D = num_cell_dims(grid_topology)

  dim_to_face_to_label = [node_to_label,]

  for d in 1:D
    z = fill(UNSET,num_faces(grid_topology,d))
    push!(dim_to_face_to_label,z)
  end

  _fill_dim_to_face_to_label!(
    dim_to_face_to_label,
    grid_topology,
    dim_to_gface_to_nodes,
    dim_gface_to_entity,
    dim_to_offset)

  tag_to_name, tag_to_groups, _ = _setup_tag_to_name(dim_to_group_to_name)

  tag_to_labels = _setup_tag_to_labels(
    tag_to_groups,dim_to_group_to_entities,dim_to_offset)

  FaceLabeling(dim_to_face_to_label,tag_to_labels,tag_to_name)

end

function _fill_dim_to_face_to_label!(
  dim_to_face_to_label,
  grid_topology,
  dim_to_gface_to_nodes,
  dim_gface_to_entity,
  dim_to_offset)

  D = length(dim_to_face_to_label)-1
  d = D
  cell_to_entity =  dim_gface_to_entity[d+1]
  cell_to_label = dim_to_face_to_label[d+1]
  offset = dim_to_offset[d+1]
  _apply_offset_for_cells!(cell_to_label,cell_to_entity,offset)

  for d in 0:(D-1)

    gface_to_nodes = dim_to_gface_to_nodes[d+1]
    face_to_nodes = get_faces(grid_topology,d,0)
    node_to_faces = get_faces(grid_topology,0,d)
    gface_to_face = _setup_gface_to_face(
      face_to_nodes,
      node_to_faces,
      gface_to_nodes)
    face_to_label = dim_to_face_to_label[d+1]
    gface_to_entity = dim_gface_to_entity[d+1]
    offset = dim_to_offset[d+1]
    _apply_offset_for_faces!(face_to_label,gface_to_entity,gface_to_face,offset)

  end

  for d = 0:(D-1)
    for j in (d+1):D
      dface_to_jfaces = get_faces(grid_topology,d,j)
      dface_to_label = dim_to_face_to_label[d+1]
      jface_to_label = dim_to_face_to_label[j+1]
      _fix_dface_to_label!(dface_to_label,jface_to_label,dface_to_jfaces)
    end
  end

end

function _apply_offset_for_cells!(cell_to_label,cell_to_entity,offset)
  ncells = length(cell_to_label)
  for cell in 1:ncells
    entity = cell_to_entity[cell]
    cell_to_label[cell] = entity+offset
  end
end

function _apply_offset_for_faces!(face_to_label,gface_to_entity,gface_to_face,offset)
  ngfaces = length(gface_to_face)
  for gface in 1:ngfaces
    entity = gface_to_entity[gface]
    face = gface_to_face[gface]
    face_to_label[face] = entity+offset
  end
end

function _setup_gface_to_face(
  face_to_nodes,
  node_to_faces,
  gface_to_nodes)

  gface_to_face = _find_gface_to_face(
    face_to_nodes.data,
    face_to_nodes.ptrs,
    node_to_faces.data,
    node_to_faces.ptrs,
    gface_to_nodes.data,
    gface_to_nodes.ptrs)

  gface_to_face

end

function _find_gface_to_face(
  face_to_nodes_data,
  face_to_nodes_ptrs,
  node_to_faces_data::AbstractVector{T},
  node_to_faces_ptrs,
  gface_to_nodes_data,
  gface_to_nodes_ptrs) where T

  ngfaces = length(gface_to_nodes_ptrs) - 1
  gface_to_face = zeros(T,ngfaces)
  n = max_cells_arround_vertex(node_to_faces_ptrs)
  faces_around = fill(UNSET,n)
  faces_around_scratch = fill(UNSET,n)

  _fill_gface_to_face!(
    gface_to_face,
    face_to_nodes_data,
    face_to_nodes_ptrs,
    node_to_faces_data,
    node_to_faces_ptrs,
    gface_to_nodes_data,
    gface_to_nodes_ptrs,
    faces_around,
    faces_around_scratch)

  gface_to_face

end

function  _fill_gface_to_face!(
  gface_to_face,
  face_to_nodes_data,
  face_to_nodes_ptrs,
  node_to_faces_data,
  node_to_faces_ptrs,
  gface_to_nodes_data,
  gface_to_nodes_ptrs,
  faces_around,
  faces_around_scratch)

  ngfaces = length(gface_to_nodes_ptrs) - 1

  nfaces_around = UNSET
  nfaces_around_scratch = UNSET

  for gface in 1:ngfaces

    a = gface_to_nodes_ptrs[gface]-1
    b = gface_to_nodes_ptrs[gface+1]
    nlnodes = b-(a+1)

    for lnode in 1:nlnodes
      node = gface_to_nodes_data[lnode+a]
      if lnode == 1
        nfaces_around = _fill_cells_around_scratch!(
          faces_around,
          node,
          node_to_faces_data,
          node_to_faces_ptrs)
      else
        nfaces_around_scratch = _fill_cells_around_scratch!(
          faces_around_scratch,
          node,
          node_to_faces_data,
          node_to_faces_ptrs)
        _set_intersection!(
          faces_around,faces_around_scratch,
          nfaces_around,nfaces_around_scratch)
      end
    end

    for face in faces_around
      if face != UNSET
        gface_to_face[gface] = face
        break
      end
    end

  end

end

function _fix_dface_to_label!(dface_to_label,jface_to_label,dface_to_jfaces)

  ndfaces = length(dface_to_label)
  @assert ndfaces == length(dface_to_jfaces)

  for dface in 1:ndfaces

    dlabel = dface_to_label[dface]
    if dlabel != UNSET
      continue
    end

    jfaces = dface_to_jfaces[dface]
    for jface in jfaces
      jlabel = jface_to_label[jface]
      if jlabel != UNSET
        dface_to_label[dface] = jlabel
        break
      end
    end

  end

end

function _setup_tag_to_name(dim_to_group_to_name)
  tag_to_name = String[]
  tag_to_groups = Vector{Tuple{Int,Int}}[]
  name_to_tag = Dict{String,Int}()
  tag = 1
  for (dim, group_to_name) in enumerate(dim_to_group_to_name)
    for (id,name) in enumerate(group_to_name)
      group = (dim,id)
      if !haskey(name_to_tag,name)
        push!(tag_to_name,name)
        groups = [group,]
        push!(tag_to_groups,groups)
        name_to_tag[name] = tag
        tag += 1
      else
        _tag = name_to_tag[name]
        groups = tag_to_groups[_tag]
        push!(groups,group)
      end
    end
  end
  (tag_to_name, tag_to_groups, name_to_tag)
end

function _setup_tag_to_labels(
  tag_to_groups,dim_to_group_to_entities,dim_to_offset)

  tag_to_labels = Vector{Int}[]
  for groups in tag_to_groups

    labels = Int[]
    for group in groups
      dim, id = group
      offset = dim_to_offset[dim]
      entities = dim_to_group_to_entities[dim][id]
      for entity in entities
        label = entity + offset
        push!(labels,label)
      end
    end
    push!(tag_to_labels,labels)

  end

  tag_to_labels

end