text:
# ===============================
# Written by AAE
# TU Delft, Winter 2014
# simulkade.com
# ===============================

# =============================== SOLVERS ===================================
function solveLinearPDE(m::MeshStructure, M::SparseMatrixCSC{Float64, Int64}, RHS::Array{Float64,1})
N=m.dims
x=M\RHS # until the problem is solved with Julia "\" solver
phi = CellValue(m, reshape(full(x), tuple(N+2...)))
phi
end

function solvePDE(m::MeshStructure, M::SparseMatrixCSC{Float64, Int64}, RHS::Array{Float64,1})
N=m.dims
x=M\RHS # until the problem is solved with Julia "\" solver
phi = CellValue(m, reshape(full(x), tuple(N+2...)))
phi
end

function solveMUMPSLinearPDE(m::MeshStructure, M::SparseMatrixCSC{Float64, Int64}, RHS::Array{Float64,1})
N=m.dims
x=solveMUMPS(M,RHS) # until the problem is solved with Julia "\" solver
phi = CellValue(m, reshape(full(x), tuple(N+2...)))
phi
end

function solveExplicitPDE(phi_old::CellValue, dt::Real, RHS::Array{Float64,1},
  BC::BoundaryCondition)
  d = phi_old.domain.dimension
  N = phi_old.domain.dims
  phi_val=reshape(phi_old.value[:]+dt*RHS, tuple(N+2...))
  if (d==1) || (d==1.5)
  	phi_val= phi_val[2:N[1]+1]
  elseif (d==2) || (d==2.5) || (d==2.8)
  	phi_val= phi_val[2:N[1]+1, 2:N[2]+1]
  elseif (d==3) || (d==3.2)
    phi_val= phi_val[2:N[1]+1, 2:N[2]+1, 2:N[3]+1]
  end
  return createCellVariable(phi_old.domain, phi_val, BC)
end

function solveExplicitPDE(phi_old::CellValue, dt::Real, RHS::Array{Float64,1},
  BC::BoundaryCondition, alfa::CellValue)
  d = phi_old.domain.dimension
  N = phi_old.domain.dims
  phi_val=reshape(phi_old.value[:]+dt*RHS./alfa.value[:], tuple(N+2...))
  if (d==1) || (d==1.5)
  	phi_val= phi_val[2:N[1]+1]
  elseif (d==2) || (d==2.5) || (d==2.8)
  	phi_val= phi_val[2:N[1]+1, 2:N[2]+1]
  elseif (d==3) || (d==3.2)
    phi_val= phi_val[2:N[1]+1, 2:N[2]+1, 2:N[3]+1]
  end
  return createCellVariable(phi_old.domain, phi_val, BC)
end

# =========================== Visualization =================================
function visualizeCells(phi::CellValue)
d=phi.domain.dimension
if d==1 || d==1.5
  x = [phi.domain.facecenters.x[1]; phi.domain.cellcenters.x; phi.domain.facecenters.x[end]]
  phi = [0.5*(phi.value[1]+phi.value[2]); phi.value[2:end-1]; 0.5*(phi.value[end-1]+phi.value[end])]
  plot(x, phi)
elseif d==2 || d==2.5
  x = [phi.domain.facecenters.x[1]; phi.domain.cellcenters.x; phi.domain.facecenters.x[end]]
  y = [phi.domain.facecenters.y[1]; phi.domain.cellcenters.y; phi.domain.facecenters.y[end]]
  phi0 = Base.copy(phi.value)
  phi0[:,1] = 0.5*(phi0[:,1]+phi0[:,2])
  phi0[1,:] = 0.5*(phi0[1,:]+phi0[2,:])
  phi0[:,end] = 0.5*(phi0[:,end]+phi0[:,end-1])
  phi0[end,:] = 0.5*(phi0[end,:]+phi0[end-1,:])
  phi0[1,1] = phi0[1,2]
  phi0[1,end] = phi0[1,end-1]
  phi0[end,1] = phi0[end,2]
  phi0[end,end] = phi0[end,end-1]
  pcolor(x, y, phi0')
elseif d==2.8
  x = [phi.domain.facecenters.x[1]; phi.domain.cellcenters.x; phi.domain.facecenters.x[end]]
  y = [phi.domain.facecenters.y[1]; phi.domain.cellcenters.y; phi.domain.facecenters.y[end]]
  phi0 = Base.copy(phi.value)
  phi0[:,1] = 0.5*(phi0[:,1]+phi0[:,2])
  phi0[1,:] = 0.5*(phi0[1,:]+phi0[2,:])
  phi0[:,end] = 0.5*(phi0[:,end]+phi0[:,end-1])
  phi0[end,:] = 0.5*(phi0[end,:]+phi0[end-1,:])
  phi0[1,1] = phi0[1,2]
  phi0[1,end] = phi0[1,end-1]
  phi0[end,1] = phi0[end,2]
  phi0[end,end] = phi0[end,end-1]
  subplot(111, polar="true")
  pcolor(y, x, phi0)
elseif d==3
  Nx = phi.domain.dims[1]
  Ny = phi.domain.dims[2]
  Nz = phi.domain.dims[3]
  x=[phi.domain.facecenters.x[1]; phi.domain.cellcenters.x; phi.domain.facecenters.x[end]]
  y=Array(Float64,1,Ny+2)
  y[:]=[phi.domain.facecenters.y[1]; phi.domain.cellcenters.y; phi.domain.facecenters.y[end]]
  z=Array(Float64,1,1,Nz+2)
  z[:]=[phi.domain.facecenters.z[1]; phi.domain.cellcenters.z; phi.domain.facecenters.z[end]]

  phi0 = Base.copy(phi.value)
  phi0[:,1,:]=0.5*(phi0[:,1,:]+phi0[:,2,:])
  phi0[:,end,:]=0.5*(phi0[:,end-1,:]+phi0[:,end,:])
  phi0[:,:,1]=0.5*(phi0[:,:,1]+phi0[:,:,1])
  phi0[:,:,end]=0.5*(phi0[:,:,end-1]+phi0[:,:,end])
  phi0[1,:,:]=0.5*(phi0[1,:,:]+phi0[2,:,:])
  phi0[end,:,:]=0.5*(phi0[end-1,:,:]+phi0[end,:,:])

  vmin = minimum(phi0)
  vmax = maximum(phi0)

  a=ones(Nx+2,Ny+2,Nz+2)
  X = x.*a
  Y = y.*a
  Z = z.*a

  s=mayavis.pipeline[:scalar_field](X,Y,Z,phi0)

  mayavis.pipeline[:image_plane_widget](s, plane_orientation="x_axes", slice_index=0, vmin=vmin, vmax=vmax)
  mayavis.pipeline[:image_plane_widget](s, plane_orientation="y_axes", slice_index=0, vmin=vmin, vmax=vmax)
  mayavis.pipeline[:image_plane_widget](s, plane_orientation="z_axes", slice_index=0, vmin=vmin, vmax=vmax)
  mayavis.pipeline[:image_plane_widget](s, plane_orientation="z_axes", slice_index=floor(Integer,Nz/2.0), vmin=vmin, vmax=vmax)
  mayavis.outline()

#   # 6 surfaces
#   # surfaces 1,2 (x=x[1], x=x[end])
#   Y=repmat(y,1,Nz)
#   Z=repmat(z,1,Ny)
#   mayavis.mesh(x[1]*ones(Ny,Nz),Y,Z',scalars=squeeze(phi.value[2,2:end-1,2:end-1],1), vmin=vmin, vmax=vmax, opacity=0.8)
#   mayavis.mesh(x[end]*ones(Ny,Nz),Y,Z',scalars=squeeze(phi.value[end-1,2:end-1,2:end-1],1), vmin=vmin, vmax=vmax, opacity=0.8)
#
#   # surfaces 3,4 (y=y[1], y=y[end]
#   X = repmat(x,1,Nz)
#   Z = repmat(z,1,Nx)
#   mayavis.mesh(X,y[1]*ones(Nx,Nz),Z',scalars=squeeze(phi.value[2:end-1,2,2:end-1],2), vmin=vmin, vmax=vmax, opacity=0.8)
#   mayavis.mesh(X,y[end]*ones(Nx,Nz),Z',scalars=squeeze(phi.value[2:end-1,end-1,2:end-1],2), vmin=vmin, vmax=vmax, opacity=0.8)
#   mayavis.axes()
#
#   # surfaces 5,6 (z=z[1], z=z[end]
#   X = repmat(x,1,Ny)
#   Y = repmat(y,1,Nx)
#   mayavis.mesh(X,Y',z[1]*ones(Nx,Ny),scalars=phi.value[2:end-1,2:end-1,2], vmin=vmin, vmax=vmax, opacity=0.8)
#   mayavis.mesh(X,Y',z[end]*ones(Nx,Ny),scalars=phi.value[2:end-1,2:end-1,end-1], vmin=vmin, vmax=vmax, opacity=0.8)

  mayavis.colorbar()
  mayavis.axes()
  mshot= mayavis.screenshot()
  mayavis.show()
  return mshot

elseif d==3.2
  Nx = phi.domain.dims[1]
  Ny = phi.domain.dims[2]
  Nz = phi.domain.dims[3]
  r=[phi.domain.facecenters.x[1]; phi.domain.cellcenters.x; phi.domain.facecenters.x[end]]
  theta = Array(Float64,1,Ny+2)
  theta[:]=[phi.domain.facecenters.y[1]; phi.domain.cellcenters.y; phi.domain.facecenters.y[end]]
  z=Array(Float64,1,1,Nz+2)
  z[:]=[phi.domain.facecenters.z[1]; phi.domain.cellcenters.z; phi.domain.facecenters.z[end]]
  a=ones(Nx+2,Ny+2,Nz+2)
  R=r.*a
  TH = theta.*a
  Z = z.*a

  X=R.*cos(TH)
  Y=R.*sin(TH)

  phi0 = Base.copy(phi.value)
  phi0[:,1,:]=0.5*(phi0[:,1,:]+phi0[:,2,:])
  phi0[:,end,:]=0.5*(phi0[:,end-1,:]+phi0[:,end,:])
  phi0[:,:,1]=0.5*(phi0[:,:,1]+phi0[:,:,1])
  phi0[:,:,end]=0.5*(phi0[:,:,end-1]+phi0[:,:,end])
  phi0[1,:,:]=0.5*(phi0[1,:,:]+phi0[2,:,:])
  phi0[end,:,:]=0.5*(phi0[end-1,:,:]+phi0[end,:,:])

  vmin = minimum(phi0)
  vmax = maximum(phi0)
  # 6 surfaces
  # surfaces 1,2 (x=x[1], x=x[end])
  mayavis.mesh(squeeze(X[floor(Integer,Nx/2.0),:,:],1),squeeze(Y[floor(Integer,Nx/2.0),:,:],1),squeeze(Z[floor(Integer,Nx/2.0),:,:],1),
    scalars=squeeze(phi0[floor(Integer,Nx/2.0)+1,:,:],1), vmin=vmin, vmax=vmax, opacity=0.8)
  mayavis.mesh(squeeze(X[Nx,:,:],1),squeeze(Y[Nx,:,:],1),squeeze(Z[Nx,:,:],1),
    scalars=squeeze(phi0[Nx+2,:,:],1), vmin=vmin, vmax=vmax, opacity=0.8)

  # surfaces 3,4 (y=y[1], y=y[end]
  mayavis.mesh(squeeze(X[:,floor(Integer,Ny/2.0),:],2),squeeze(Y[:,floor(Integer,Ny/2.0),:],2),squeeze(Z[:,floor(Integer,Ny/2.0),:],2),
    scalars=squeeze(phi0[:,floor(Integer,Ny/2.0)+1,:],2), vmin=vmin, vmax=vmax, opacity=0.8)
  mayavis.mesh(squeeze(X[:,Ny,:],2),squeeze(Y[:,Ny,:],2),squeeze(Z[:,Ny,:],2),
    scalars=squeeze(phi0[:,Ny+2,:],2), vmin=vmin, vmax=vmax, opacity=0.8)

  # surfaces 5,6 (z=z[1], z=z[end]
  mayavis.mesh(X[:,:,floor(Integer,Nz/2.0)],Y[:,:,floor(Integer,Nz/2.0)],Z[:,:,floor(Integer,Nz/2.0)],
    scalars=phi0[:,:,floor(Integer,Nz/2.0)+1], vmin=vmin, vmax=vmax, opacity=0.8)
  mayavis.mesh(X[:,:,Nz],Y[:,:,Nz],Z[:,:,Nz],
    scalars=phi0[:,:,Nz+1], vmin=vmin, vmax=vmax, opacity=0.8)
  mayavis.colorbar()
  mayavis.axes()
  mshot=mayavis.screenshot()
  mayavis.show()
  return mshot
end
end

########################## Visualize Vectors #####################
function visualizeCellVectors(phi::CellVector)
d=phi.domain.dimension
if d==1 || d==1.5
  println("Vector visualization works only in 2D and 3D")
elseif d==2 || d==2.5
  x = phi.domain.cellcenters.x
  y = phi.domain.cellcenters.y
  quiver(repmat(x, 1, length(y)), repmat(y', length(x), 1), phi.xvalue, phi.yvalue)
elseif d==2.8
  x = phi.domain.cellcenters.x
  y = phi.domain.cellcenters.y'
  subplot(111, polar="true")
  quiver(repmat(y, length(x), 1), repmat(x, 1, length(y)),
  phi.xvalue.*cos(y)-phi.yvalue.*sin(y), phi.xvalue.*sin(y)+phi.yvalue.*cos(y))
elseif d==3
  Nx = phi.domain.dims[1]
  Ny = phi.domain.dims[2]
  Nz = phi.domain.dims[3]
  x=phi.domain.cellcenters.x
  y=Array(Float64,1,Ny)
  y[:]=phi.domain.cellcenters.y
  z=Array(Float64,1,1,Nz)
  z[:]=phi.domain.cellcenters.z

  #vmin = minimum(phi.xvalue)
  #vmax = maximum(phi0)

  a=ones(Nx,Ny,Nz)
  X = x.*a
  Y = y.*a
  Z = z.*a

  mayavis.quiver3d(X,Y,Z, phi.xvalue, phi.yvalue, phi.zvalue)
  mayavis.outline()

#   # 6 surfaces
#   # surfaces 1,2 (x=x[1], x=x[end])
#   Y=repmat(y,1,Nz)
#   Z=repmat(z,1,Ny)
#   mayavis.mesh(x[1]*ones(Ny,Nz),Y,Z',scalars=squeeze(phi.value[2,2:end-1,2:end-1],1), vmin=vmin, vmax=vmax, opacity=0.8)
#   mayavis.mesh(x[end]*ones(Ny,Nz),Y,Z',scalars=squeeze(phi.value[end-1,2:end-1,2:end-1],1), vmin=vmin, vmax=vmax, opacity=0.8)
#
#   # surfaces 3,4 (y=y[1], y=y[end]
#   X = repmat(x,1,Nz)
#   Z = repmat(z,1,Nx)
#   mayavis.mesh(X,y[1]*ones(Nx,Nz),Z',scalars=squeeze(phi.value[2:end-1,2,2:end-1],2), vmin=vmin, vmax=vmax, opacity=0.8)
#   mayavis.mesh(X,y[end]*ones(Nx,Nz),Z',scalars=squeeze(phi.value[2:end-1,end-1,2:end-1],2), vmin=vmin, vmax=vmax, opacity=0.8)
#   mayavis.axes()
#
#   # surfaces 5,6 (z=z[1], z=z[end]
#   X = repmat(x,1,Ny)
#   Y = repmat(y,1,Nx)
#   mayavis.mesh(X,Y',z[1]*ones(Nx,Ny),scalars=phi.value[2:end-1,2:end-1,2], vmin=vmin, vmax=vmax, opacity=0.8)
#   mayavis.mesh(X,Y',z[end]*ones(Nx,Ny),scalars=phi.value[2:end-1,2:end-1,end-1], vmin=vmin, vmax=vmax, opacity=0.8)

  mayavis.colorbar()
  mayavis.axes()
  mshot= mayavis.screenshot()
  mayavis.show()
  return mshot

elseif d==3.2
  Nx = phi.domain.dims[1]
  Ny = phi.domain.dims[2]
  Nz = phi.domain.dims[3]
  r=phi.domain.cellcenters.x
  theta = Array(Float64,1,Ny)
  theta[:]=phi.domain.cellcenters.y
  z=Array(Float64,1,1,Nz)
  z[:]=phi.domain.cellcenters.z
  a=ones(Nx,Ny,Nz)
  R=r.*a
  TH = theta.*a
  Z = z.*a

  X=R.*cos(TH)
  Y=R.*sin(TH)

  #vmin = minimum(phi0)
  #vmax = maximum(phi0)
  # 6 surfaces
  # surfaces 1,2 (x=x[1], x=x[end])
  mayavis.quiver3d(X,Y,Z, phi.xvalue.*cos(TH)-phi.yvalue.*sin(TH),
  phi.xvalue.*sin(TH)+phi.yvalue.*cos(TH), phi.zvalue)
  mayavis.colorbar()
  mayavis.axes()
  mshot=mayavis.screenshot()
  mayavis.show()
  return mshot
end
end


text:
const RL = 3.09727
const BL = 1.9856
const FIELD_W = 220
const FIELD_H = 180
const SPEED_MAX = 125

function instruction(s::MyStrategy, field, id)
    self = field.self_robots[id]
    if id in [2,3,4,5]
        state = s.robot_state[id]
        return move(state, self, field.ball)
    end
    return (0, 0)
end

function drive(deg, maxspeed = SPEED_MAX)
    ω = deg / 180 * π
    v1 = maxspeed
    v2 = v1 - abs(ω) * RL
    return ω > 0 ? (v2, v1) : (v1, v2)
end

function angle(origin, target)
    (dx, dy) = target .- origin
    atand(dy, dx)
end

function anglediff(d1, d2)
    d=(d2 - d1 + 180) % 360 - 180
    return d < -180 ? d + 360 : d
end

function move(state, self, target)
    at = angle(self.position, target)
    as = self.rotation
    e = anglediff(as, at)
    #addsample(e)
    e_int = state[:e_int]
    e_last = state[:e_last]
    # You may need to adjust these parameters
    Kp = 5.5
    Ki = 0.25
    Kd = 2.7
    u = Kp * e+Ki * e_int + Kd * (e - e_last)
    state[:e_int] = e_int + e
    state[:e_last] = e
    drive(u, 50)
end

ext: jl

max_stars_repo_name: Runfa-Zhang/MeetUpMaterials

max_issues_repo_name: Runfa-Zhang/MeetUpMaterials

max_forks_repo_name: Runfa-Zhang/MeetUpMaterials


text:
using Test, ANSIColoredPrinters

function repr_color(printer::PlainTextPrinter)
    result = repr("text/plain", printer, context = :color => true)
    get(ENV, "JULIA_DEBUG", "") == "" || println(result)
    take!(printer.buf)
    return result
end

@testset "no color" begin
    buf = IOBuffer()
    printer = PlainTextPrinter(buf)
    print(buf, "\e[36m", " CyanFG ")
    print(buf, "\e[39m", " Normal ")
    result = repr("text/plain", printer)
    @test result == " CyanFG  Normal "
end

@testset "single modification" begin
    buf = IOBuffer()
    printer = PlainTextPrinter(buf)

    @testset "bold/faint" begin
        print(buf, "\e[0m", " Normal ")
        print(buf, "\e[1m", " Bold ")
        print(buf, "\e[2m", " Faint ")
        print(buf, "\e[22m", " Normal ")
        @test repr_color(printer) == " Normal \e[1m Bold \e[2m Faint \e[m Normal "
    end

    @testset "fg/bg" begin
        print(buf, "\e[0m", " Normal ")
        print(buf, "\e[36m", " CyanFG ")
        print(buf, "\e[39;44m", " BlueBG ")
        print(buf, "\e[49m", " Normal ")
        @test repr_color(printer) == " Normal \e[36m CyanFG \e[39;44m BlueBG \e[m Normal "
    end

    @testset "blink" begin
        print(buf, "\e[0m", " Normal ")
        print(buf, "\e[5m", " Blink ")
        print(buf, "\e[6m", " RapidBlink ")
        print(buf, "\e[25m", " Normal ")
        @test repr_color(printer) == " Normal \e[5m Blink \e[6m RapidBlink \e[m Normal "
    end
end

@testset "nested modification" begin
    buf = IOBuffer()
    printer = PlainTextPrinter(buf)

    @testset "bold/italic" begin
        print(buf, "\e[0m", " Normal ")
        print(buf, "\e[1m", " Bold ")
        print(buf, "\e[3m", " Bold-Italic ")
        print(buf, "\e[23m", " Bold ")
        print(buf, "\e[0m", " Normal ")
        result = repr_color(printer)
        @test result == " Normal \e[1m Bold \e[3m Bold-Italic \e[23m Bold \e[m Normal "

        print(buf, "\e[3m", " Italic ")
        print(buf, "\e[1m", " Bold-Italic ")
        print(buf, "\e[22m", " Italic ")
        print(buf, "\e[0m", " Normal ")
        result = repr_color(printer)
        @test result == "\e[3m Italic \e[1m Bold-Italic \e[22m Italic \e[m Normal "
    end

    @testset "bold/fg" begin
        print(buf, "\e[0m", " Normal ")
        print(buf, "\e[1m", " Bold ")
        print(buf, "\e[36m", " Bold-CyanFG ")
        print(buf, "\e[39m", " Bold ")
        print(buf, "\e[0m", " Normal ")
        result = repr_color(printer)
        @test result == " Normal \e[1m Bold \e[36m Bold-CyanFG \e[39m Bold \e[m Normal "
    end

    @testset "strike/blink" begin
        print(buf, "\e[0m", " Normal ")
        print(buf, "\e[9m", " Strike ")
        print(buf, "\e[5m", " Strike-Blink ")
        print(buf, "\e[25m", " Strike ")
        print(buf, "\e[29m", " Normal ")
        result = repr_color(printer)
        @test result == " Normal \e[9m Strike \e[5m Strike-Blink \e[25m Strike \e[m Normal "
    end

    @testset "fg/bg" begin
        print(buf, "\e[0m", " Normal ")
        print(buf, "\e[96m", " LightCyanFG ")
        print(buf, "\e[45m", " LightCyanFG-MagentaBG ")
        print(buf, "\e[49m", " LightCyanFG ")
        print(buf, "\e[39m", " Normal ")
        result = repr_color(printer)
        @test result == " Normal \e[96m LightCyanFG \e[45m LightCyanFG-MagentaBG " *
                        "\e[49m LightCyanFG \e[m Normal "
    end
end

@testset "overlapped modification" begin
    buf = IOBuffer()
    printer = PlainTextPrinter(buf)

    @testset "bold/fg" begin
        print(buf, "\e[0m", " Normal ")
        print(buf, "\e[1m", " Bold ")
        print(buf, "\e[36m", " Bold-CyanFG ")
        print(buf, "\e[22;36m", " CyanFG ")
        print(buf, "\e[39m", " Normal ")
        result = repr_color(printer)
        @test result == " Normal \e[1m Bold \e[36m Bold-CyanFG " *
                        "\e[22;36m CyanFG \e[m Normal "
    end

    @testset "fg/bg" begin
        print(buf, "\e[0m", " Normal ")
        print(buf, "\e[96m", " LightCyanFG ")
        print(buf, "\e[45m", " LightCyanFG-MagentaBG ")
        print(buf, "\e[39m", " MagentaBG ")
        print(buf, "\e[49m", " Normal ")
        result = repr_color(printer)
        @test result == " Normal \e[96m LightCyanFG \e[45m LightCyanFG-MagentaBG " *
                        "\e[39m MagentaBG \e[m Normal "
    end
end

@testset "force reset" begin
    buf = IOBuffer()
    printer = PlainTextPrinter(buf)

    print(buf, "\e[7m", " Invert ")
    print(buf, "\e[8m", " Conceal ")
    @test repr_color(printer) == "\e[7m Invert \e[8m Conceal \e[m"
end

@testset "256 colors" begin
    buf = IOBuffer()
    printer = PlainTextPrinter(buf)

    print(buf, "\e[38;5;4m", " Blue(FG) ")
    print(buf, "\e[48;5;16m", " #000(BG) ")
    print(buf, "\e[38;5;110m", " #87afd7(FG) ")
    print(buf, "\e[48;5;255m", " #eee(BG) ")

    result = repr_color(printer)
    @test result == "\e[34m Blue(FG) \e[48;5;16m #000(BG) " *
                    "\e[38;5;110m #87afd7(FG) \e[48;5;255m #eee(BG) \e[m"
end

@testset "24-bit colors" begin
    buf = IOBuffer()
    printer = PlainTextPrinter(buf)

    print(buf, "\e[38;2;0;128;255m", " #0080ff(FG) ")
    print(buf, "\e[48;2;238;238;238m", " #eee(BG) ")
    print(buf, "\e[38;2;170;187;204m", " #abc(FG) ")

    result = repr_color(printer)
    @test result == "\e[38;2;0;128;255m #0080ff(FG) \e[48;5;255m #eee(BG) " *
                    "\e[38;2;170;187;204m #abc(FG) \e[m"
end


text:
@testset "Testing IMT submodule" begin

fnMeas = joinpath(datadir,"imt","measurementIMT2D.h5")
fnCalib = joinpath(datadir,"imt","systemMatrixIMT2D.h5")

measIMT = MPIFile(fnMeas)
calibIMT = MPIFile(fnCalib)
@test typeof(measIMT) == IMTFileMeas
@test typeof(calibIMT) == IMTFileCalib

@test size( measData(measIMT) ) == (1632, 2, 1, 1)
@test size( measData(calibIMT) ) == (400, 817, 2, 1)

for imt in (measIMT, calibIMT)
  println("Test $imt")
  @test studyName(imt) == "n.a."
  @test studyNumber(imt) == 0
  @test studyDescription(imt) == "n.a."
  @test studyTime(imt) == nothing

  @test experimentName(imt) == "n.a."
  @test experimentNumber(imt) == 0
  @test experimentDescription(imt) == "n.a."
  @test experimentSubject(imt) == "n.a."
  @test experimentIsSimulation(imt) == true

  @test scannerFacility(imt) == "n.a."
  @test scannerOperator(imt) == "n.a."
  @test scannerManufacturer(imt) == "n.a."
  @test scannerName(imt) == "n.a."
  @test scannerTopology(imt) == "n.a."

  @test tracerName(imt) == ["n.a."]
  @test tracerBatch(imt) == ["n.a."]
  @test tracerVendor(imt) == ["n.a."]
  @test tracerVolume(imt) == [0.0]
  @test tracerConcentration(imt) == [0.0]
  @test tracerInjectionTime(imt) == [Dates.unix2datetime(0)]

  @test acqStartTime(imt) == Dates.unix2datetime(0)
  @test acqGradient(imt)[:,:,1,1] == [0 0 0; 0 0 0; 0 0 0]
  @test acqFramePeriod(imt) == 0.0006528
  @test acqNumPeriodsPerFrame(imt) == 1

  #@test acqOffsetFieldShift(imt)[:,1,1] == [NaN; NaN; NaN]
  @test acqNumAverages(imt) == 1
  @test acqNumFrames(imt) == 1
  #@test acqOffsetField(imt) == [0.0; 0.0; 0.0]
  @test acqNumPeriods(imt) == 1

  @test dfNumChannels(imt) == 2
  @test dfWaveform(imt) == "sine"
  @test dfStrength(imt)[:,:,1] == [0.0 0.0 0.0]
  @test dfPhase(imt)[:,:,1] == [0.0 0.0 0.0] #[1.5707963267948966 1.5707963267948966 1.5707963267948966]
  @test dfBaseFrequency(imt) == 2500000.0
  @test dfDivider(imt)[:,1] == [102; 96; 99]
  @test dfCycle(imt) == 0.0006528

  @test rxNumChannels(imt) == 2
  @test rxBandwidth(imt) == 1.25e6
  @test rxNumSamplingPoints(imt) == 1632
  @test rxUnit == MPIFiles.rxUnit
  @test rxDataConversionFactor(imt) == reshape(Float64[1.0 0.0 1.0 0.0], 2, 2)

  #@test size(getMeasurements(imt, numAverages=1,
  #            spectralLeakageCorrection=false)) == (53856,3,1,1)

  #@test size(getMeasurements(imt, numAverages=10,
  #            spectralLeakageCorrection=false)) == (1632,3,1,50)

  #@test size(getMeasurements(imt, numAverages=10, frames=1:500,
  #            spectralLeakageCorrection=true)) == (1632,3,1,50)

  #@test size(getMeasurementsFD(imt, numAverages=10, frames=1:500)) == (817,3,1,50)

  #@test size(getMeasurementsFD(imt, numAverages=10, frames=1:500, loadasreal=true)) == (1634,3,1,50)

  #@test size(getMeasurementsFD(imt,frequencies=1:10, numAverages=10)) == (10,1,50)

end

# Calibration File

  @test size( measData(calibIMT) ) == (400, 817, 2, 1)
  @test measIsFourierTransformed(calibIMT) == true
  @test measIsTFCorrected(calibIMT) == false
  @test measIsFastFrameAxis(calibIMT) == true
  @test measIsBGCorrected(calibIMT) == false

  @test size( calibSNR(calibIMT) ) == (817, 2, 1)
  @test calibFov(calibIMT) == [0.01, 0.01, 0.0001]
  @test calibFovCenter(calibIMT) == [0.0; 0.0; 0.0]
  @test calibSize(calibIMT) == [20; 20; 1]
  @test calibOrder(calibIMT) == "xyz"
  @test calibPositions(calibIMT) == nothing
  @test calibOffsetField(calibIMT) == nothing
  @test calibDeltaSampleSize(calibIMT) == [0.0; 0.0; 0.0]
  @test calibMethod(calibIMT) == "simulation"

#  @test size(filterFrequencies(sm, SNRThresh = 5)) == (147,)
#  #@test size(filterFrequencies(sm, numUsedFreqs = 100)) == (100,) # not working

#  @test size(getSystemMatrix(sm,1:10)) == (1936,10)
#  @test size(getSystemMatrix(sm,1:10,loadasreal=true)) == (1936,20)

end


text:
function my_retry(f::Function)
    delays = ExponentialBackOff(;
        n = 10,
        max_delay = 30.0,
    )
    return retry(f; delays)()
end


text:
DataGeneratingDistibutions = Union{Vector{UnivariateDistribution},Vector{Vector{UnivariateDistribution}}}
"""
    convert(::DataGeneratingDistibutions, v)

Converting DataGeneratingDistibutions to specific vector.
"""
function Base.convert(::Type{DataGeneratingDistibutions}, v::Vector{T}) where T <: UnivariateDistribution
    convert(Vector{UnivariateDistribution}, v)
end
function Base.convert(::Type{DataGeneratingDistibutions}, v::Vector{Vector{T}}) where T <: UnivariateDistribution
    convert(Vector{Vector{UnivariateDistribution}}, v)
end

"""
    Simulation(experiment, parameters, datagendists, maxsteps, onestepsizes, minsteps)

A simulation setup includes the experiment, data generating distributions, max number of steps
and minimum number of steps.
"""
struct Simulation
    experiment::Experiment
    parameters::Vector{Symbol}
    datagendists::DataGeneratingDistibutions
    maxsteps::Integer
    onestepsizes::Union{Vector{Int},Vector{Vector{Int}}}
    minsteps::Integer
    function Simulation(;experiment, datagendists,
        maxsteps, onestepsizes, minsteps=0, parameters::Union{Vector{Symbol},Nothing}=nothing)
        minsteps <= maxsteps || error("maxsteps should be equal or larger than minsteps")
        if parameters === nothing
            return new(experiment, defaultparams(experiment), datagendists, maxsteps, onestepsizes, minsteps)
        else
            return new(experiment, parameters, datagendists, maxsteps, onestepsizes, minsteps)
        end
    end
end

function Base.show(io::IO, simulation::Simulation)
    info = """
    $(simulation.experiment)
    maxsteps  : $(simulation.maxsteps)
    minstpes  : $(simulation.minsteps)
    """
    print(io, typeof(simulation), "\n", info)
end

function updateonce!(experiment::Experiment, datagendists, onestepsizes)
    length(experiment.models) == length(datagendists) ||
        error("Number of models should be equal to number of data gen distributions.")
    length(datagendists) == length(onestepsizes) ||
        error("Number of one step sizes should be equal to number of data gen distributions.")
    for (modelname, datagendist, onestepsize) in zip(experiment.modelnames, datagendists, onestepsizes)
        model = experiment.models[modelname]
        stats = samplestats(model, datagendist, onestepsize)
        update!(experiment, modelname, stats)
    end
    nothing
end

updateonce!(simulation::Simulation) = updateonce!(simulation.experiment, simulation.datagendists, simulation.onestepsizes)

function runonce(simulation::Simulation, numsamples::Integer)
    experiment = simulation.experiment
    runnum = 0
    metricvals = Vector{Vector{Float64}}()
    winner = nothing
    while experiment.winner === nothing && runnum < simulation.maxsteps
        runnum += 1
        updateonce!(simulation)
        _, metricval = metrics(experiment, simulation.parameters, numsamples=numsamples)
        push!(metricvals, metricval)
        if runnum < simulation.minsteps
            continue
        else
            winner, _ = decide!(experiment, simulation.parameters, numsamples=numsamples)
        end
    end
    metricvals, winner
end

function record!(winners, winner)
    if winner === nothing
        push!(winners, "nothing")
    else
        push!(winners, winner)
    end
    nothing
end

function runsequential(simulation::Simulation; numsamples=10_000, numsims=50)
    numsteps = Vector{Int}()
    winners = Vector{String}()
    metricvals = Vector{Vector{Vector{Float64}}}()
    for _ = 1:numsims
        simcopy = deepcopy(simulation)
        metrival, winner = runonce(simcopy, numsamples)
        push!(metricvals, metrival)
        push!(numsteps, length(metrival))
        record!(winners, winner)
    end
    numsteps, winners, metricvals
end

# TODO: speed up using StaticArrays and multiple threads.


text:
using Revise, Parameters, Logging, Test
using CUDA
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils
using ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils

##### iter data block
threads=(32,5);
blocks =17;
mainArrDims= (67,177,90);
dataBdim = (43,21,17);
mainArrCPU= falses(mainArrDims);
mainArrCPU[5,5,5]= true;
mainArrGPU = CuArray(mainArrCPU);
metaData = MetaDataUtils.allocateMetadata(mainArrDims,dataBdim);
metaDataDims= size(metaData);
inBlockLoopX,inBlockLoopY,inBlockLoopZ= (fld(dataBdim[1] ,threads[1]),fld(dataBdim[2] ,threads[2]),dataBdim[3]    );
#we are iterating here block by block sequentially
loopXMeta,loopYZMeta= (1,1)#(metaDataDims[1],fld(metaDataDims[2]*metaDataDims[3] ,blocks)  )

resShmem = CUDA.zeros(Bool,(dataBdim[1]+2,dataBdim[2]+2,dataBdim[3]+2 ));
sourceShmem = CUDA.zeros(Bool,(dataBdim));

function processDataKernel(mainArrGPU,sourceShmem,resShmem,mainArrDims,metaDataDims,loopXMeta,loopYZMeta,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ )
    @iter3dOuter(metaDataDims,loopXMeta,loopYZMeta,yTimesZmeta,
    begin
        # @ifXY 1 1    CUDA.@cuprint "  xMeta $(xMeta) yMeta $(yMeta)  zMeta $(zMeta) \n"

        @iterDataBlock(mainArrDims,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ,xMeta,yMeta,zMeta,
        begin
            #@ifXY 1 1    CUDA.@cuprint "  x $(x) y $(y)  z $(z) \n"

            maskBool=mainArrGPU[x,y,z]
            ProcessMainDataVerB.@processMaskData( maskBool)

    end)end)

    return
end

@cuda threads=threads blocks=blocks processDataKernel(mainArrGPU,sourceShmem,resShmem,mainArrDims,metaDataDims,loopXMeta,loopYZMeta,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ )
@test  resShmem[5+1,5+1,5+1]==false
@test  resShmem[7,6,6]==true
@test  resShmem[4+1,5+1,5+1]==true
@test  resShmem[5+1,4+1,5+1]==true
@test  resShmem[5+1,6+1,5+1]==true
@test  resShmem[5+1,5+1,4+1]==true
@test  resShmem[5+1,5+1,6+1]==true

########## loadMainValues
using Revise, Parameters, Logging, Test
using CUDA
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils
using ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils

threads=(32,5);
blocks =17;
mainArrDims= (67,177,90);
dataBdim = (43,21,17);
mainArrCPU= falses(mainArrDims);
mainArrCPU[5,5,5]= true;
mainArrCPU[33,10,7]= true;
mainArrGPU = CuArray(mainArrCPU);
metaData = MetaDataUtils.allocateMetadata(mainArrDims,dataBdim);
metaDataDims= size(metaData);
inBlockLoopX,inBlockLoopY,inBlockLoopZ= (fld(dataBdim[1] ,threads[1]),fld(dataBdim[2] ,threads[2]),dataBdim[3]    );
#we are iterating here block by block sequentially
loopXMeta,loopYZMeta= (1,1)#(metaDataDims[1],fld(metaDataDims[2]*metaDataDims[3] ,blocks)  )

resShmem = CUDA.zeros(Bool,(dataBdim[1]+2,dataBdim[2]+2,dataBdim[3]+2 ));
sourceShmem = CUDA.zeros(Bool,(dataBdim));

function processDataKernel(mainArrGPU,sourceShmem,resShmem,mainArrDims,metaDataDims,loopXMeta,loopYZMeta,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ )
    @iter3dOuter(metaDataDims,loopXMeta,loopYZMeta,yTimesZmeta,
    begin
        #@ifXY 1 1    CUDA.@cuprint "  xMeta $(xMeta) yMeta $(yMeta)  zMeta $(zMeta) \n"

        @iterDataBlock(mainArrDims,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ,xMeta,yMeta,zMeta,
        begin
            ProcessMainDataVerB.@loadMainValues( mainArrGPU,dataBdim)

    end)end)

    return
end

@cuda threads=threads blocks=blocks processDataKernel(mainArrGPU,sourceShmem,resShmem,mainArrDims,metaDataDims,loopXMeta,loopYZMeta,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ )
@test  resShmem[5+1,5+1,5+1]==false
@test  resShmem[7,6,6]==true
@test  resShmem[4+1,5+1,5+1]==true
@test  resShmem[5+1,4+1,5+1]==true
@test  resShmem[5+1,6+1,5+1]==true
@test  resShmem[5+1,5+1,4+1]==true
@test  resShmem[5+1,5+1,6+1]==true

@test sum(resShmem)==12

@test  resShmem[34,10,8]==true
@test  sourceShmem[5 ,5 ,5 ]==true
@test  sourceShmem[33 ,10 ,7 ]==true

filter(cart-> resShmem[cart], CartesianIndices(resShmem)  )

#################  paddingIter
using Revise, Parameters, Logging, Test
using CUDA
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils
using ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils,..ResultListUtils

threads=(32,5);
blocks =1;
mainArrDims= (67,177,90);
 #dataBdim = (7,3,1);
dataBdim = (43,21,17);
mainArrCPU= falses(mainArrDims);
mainArrCPU[5,5,5]= true;
mainArrCPU[33,10,7]= true;
mainArrGPU = CuArray(mainArrCPU);
refArrCPU = falses(mainArrDims);
refArrGPU= CuArray(refArrCPU);

refArrCPU[dataBdim[1]+1,dataBdim[2]+2,dataBdim[3]+2]

metaData = MetaDataUtils.allocateMetadata(mainArrDims,dataBdim);
metaDataDims= size(metaData);
totalFp,totalFn= 500,500
resList,resListIndicies,maxResListIndex= allocateResultLists(totalFp,totalFn)
iterNumb=3
isGold = 1
xMeta,yMeta,zMeta = 2,2,2
loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,inBlockLoopX,inBlockLoopY,inBlockLoopZ,metaDataLength,loopMeta,loopWarpMeta,clearIterResShmemLoop,clearIterSourceShmemLoop,clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength,sourceShmemTotalLength =calculateLoopsIter(dataBdim,threads[1],threads[2],metaDataDims,blocks)
singleVal = CUDA.zeros(Int32,1)

#for start we will get top padding
function paddingIterKernel(singleVal,metaDataDims,refArrGPU,mainArrDims,loopAZFixed,loopBZfixed,dataBdim,mainArrGPU,iterNumb,isGold,xMeta,yMeta,zMeta)
    # resShmem =  @cuDynamicSharedMem(Bool,(dataBdim[1]+2,dataBdim[2]+2,dataBdim[3]+2)) # we need this additional 33th an 34th spots
    # sync_threads()
    # @ifXY 1 1 resShmem[1,2,2]= true
    # @ifXY 1 2 resShmem[1,2,3]= true
    # @ifXY 1 3 resShmem[1,3,2]= true

    # sync_threads()

    # paddingIter(loopX,loopY,maxXdim, maxYdim,a,b,c , xMetaChange,yMetaChange,zMetaChange, mainArr,refArr, dir,iterNumb,queueNumber,xMeta,yMeta,zMeta)

    @paddingIter(loopAZFixed,loopBZfixed,dataBdim[1], dataBdim[2],begin
    CUDA.@atomic singleVal[1]+=1
end )

    return
end

@cuda threads=threads blocks=blocks paddingIterKernel(singleVal,metaDataDims,refArrGPU,mainArrDims,loopAZFixed,loopBZfixed,dataBdim,mainArrGPU,iterNumb,isGold,xMeta,yMeta,zMeta)
@test  sum(singleVal) == dataBdim[1]*dataBdim[2]

#################  paddingProcessCombined
using Revise, Parameters, Logging, Test
using CUDA
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils
using ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils,..ResultListUtils

threads=(32,5);
blocks =1;
mainArrDims= (67,177,90);
dataBdim = (43,21,17);
mainArrCPU= falses(mainArrDims);
mainArrCPU[5,5,5]= true;
mainArrCPU[33,10,7]= true;
mainArrGPU = CuArray(mainArrCPU);
refArrCPU = falses(mainArrDims);

metaData = MetaDataUtils.allocateMetadata(mainArrDims,dataBdim);
metaDataDims= size(metaData);
totalFp,totalFn= 500,500
resList,resListIndicies,maxResListIndex= allocateResultLists(totalFp,totalFn)
iterNumb=3
isGold = 1
queueNumber = 11
xMeta,yMeta,zMeta = 2,2,2
isToBeValidated = true
loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,inBlockLoopX,inBlockLoopY,inBlockLoopZ,metaDataLength,loopMeta,loopWarpMeta,clearIterResShmemLoop,clearIterSourceShmemLoop,clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength,sourceShmemTotalLength = calculateLoopsIter(dataBdim,threads[1],threads[2],metaDataDims,blocks)
resShmem = CUDA.zeros(Bool,(dataBdim[1]+2,dataBdim[2]+2,dataBdim[3]+2)) # we need this additional 33th an 34th spots
resShmem[2,2,1]=true
resShmem[3,2,1]=true
resShmem[4,3,1]=true #for start we will get top padding
refArrCPU[dataBdim[1]+2,dataBdim[2]+2,dataBdim[3]+1]= true
refArrGPU= CuArray(refArrCPU);

resShmem[queueNumber,dataBdim[2],1]=true

dataBdim[1]+2
dataBdim[2]+2
dataBdim[3]+1
metaData[xMeta,yMeta,zMeta-1,getResOffsetsBeg()+queueNumber]=7

function paddingIterKernel(resShmem,metaData,resList,resListIndicies,maxResListIndex,metaDataDims,refArrGPU,mainArrDims,loopAZFixed,loopBZfixed,dataBdim,mainArrGPU,iterNumb,isGold,xMeta,yMeta,zMeta)

    # paddingIter(loopX,loopY,maxXdim, maxYdim,a,b,c , xMetaChange,yMetaChange,zMetaChange, mainArr,refArr, dir,iterNumb,queueNumber,xMeta,yMeta,zMeta)

    @paddingProcessCombined(loopAZFixed,loopBZfixed,dataBdim[1], dataBdim[2], x,y,1 ,0,0,-1,  mainArrGPU,refArrGPU,5,iterNumb, 12-isGold,xMeta,yMeta,zMeta,isGold)

    return
end

@cuda threads=threads blocks=blocks paddingIterKernel(resShmem,metaData,resList,resListIndicies,maxResListIndex,metaDataDims,refArrGPU,mainArrDims,loopAZFixed,loopBZfixed,dataBdim,mainArrGPU,iterNumb,isGold,xMeta,yMeta,zMeta)
@test Int64(sum(resListIndicies))>0
@test  Int64.(Array(resList[8,:]))==[dataBdim[1]+2,dataBdim[2]+2,dataBdim[3]+1,1,5,3]

######### process padding
using Revise, Parameters, Logging, Test
using CUDA
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils
using ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils,..ResultListUtils

threads=(32,5);
blocks =1;
mainArrDims= (67,177,90);
dataBdim = (10,10,10)
mainArrCPU= falses(mainArrDims);
mainArrCPU[5,5,5]= true;
mainArrCPU[33,10,7]= true;
refArrCPU = falses(mainArrDims);

metaData = MetaDataUtils.allocateMetadata(mainArrDims,dataBdim);
metaDataDims= size(metaData);
totalFp,totalFn= 500,500
resList,resListIndicies,maxResListIndex= allocateResultLists(totalFp,totalFn)
iterNumb=3
isGold = 1
queueNumber = 11
xMeta,yMeta,zMeta = 1,1,1 # here by assumption that we have all meta indicies 0 based
isToBeValidated = true
loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,inBlockLoopX,inBlockLoopY,inBlockLoopZ,metaDataLength,loopMeta,loopWarpMeta,clearIterResShmemLoop,clearIterSourceShmemLoop,clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength,sourceShmemTotalLength = calculateLoopsIter(dataBdim,threads[1],threads[2],metaDataDims,blocks)
resShmem = CUDA.zeros(Bool,(dataBdim[1]+2,dataBdim[2]+2,dataBdim[3]+2)) # we need this additional 33th an 34th spots
# resShmem[2,2,1]=true
# resShmem[3,2,1]=true
# resShmem[4,3,1]=true #for start we will get top padding
refArrCPU[dataBdim[1]+2,dataBdim[2]+2,dataBdim[3]+1]= true
#setting some trues in both
begin

    # so we should get all mins as 2 and all maxes as 5

    # 1)   Left FP
    # 2)   Left FN
    mainArrCPU[11,13,13]=true
    refArrCPU[11,13,15]=true

   # 3)   Right FP
   # 4)   Right FN
    mainArrCPU[20,13,13]=true
    refArrCPU[20,13,15]=true

    mainArrCPU[20,14,13]=true
    refArrCPU[20,14,15]=true

 # 5)   Posterior FP

    # 6)   Posterior FN

    mainArrCPU[13,11,14]=true
    refArrCPU[14,11,13]=true

    mainArrCPU[15,11,13]=true
    refArrCPU[13,11,15]=true

    mainArrCPU[16,11,13]=true
    refArrCPU[13,11,16]=true

    # 7)   Anterior FP
    # 8)   Anterior FN
    mainArrCPU[14,20,13]=true
    refArrCPU[13,20,14]=true

    mainArrCPU[15,20,13]=true
    refArrCPU[14,20,16]=true

    mainArrCPU[17,20,13]=true
    refArrCPU[15,20,17]=true

    mainArrCPU[18,20,17]=true
    refArrCPU[16,20,18]=true

    # 9)   Top FP
    # 10)   Top FN
    mainArrCPU[14,15,11]=true
    refArrCPU[14,16,11]=true

    mainArrCPU[17,15,11]=true
    refArrCPU[17,16,11]=true

    mainArrCPU[15,15,11]=true
    refArrCPU[15,16,11]=true

    mainArrCPU[16,15,11]=true
    refArrCPU[16,16,11]=true

    mainArrCPU[18,15,11]=true
    refArrCPU[18,16,11]=true

    # 11)   Bottom FP
    # 12)   Bottom FN

    mainArrCPU[14,15,20]=true
    refArrCPU[14,16,20]=true

    mainArrCPU[17,15,20]=true
    refArrCPU[17,16,20]=true

    mainArrCPU[15,15,20]=true
    refArrCPU[15,16,20]=true

    mainArrCPU[16,15,20]=true
    refArrCPU[16,16,20]=true

    mainArrCPU[18,15,20]=true
    refArrCPU[18,16,20]=true

    mainArrCPU[19,15,20]=true
    refArrCPU[19,16,20]=true

    # 13)   main block Fp
    # 14)   main block Fn

    mainArrCPU[14,15,14]=true
    refArrCPU[14,16,16]=true

    mainArrCPU[17,15,17]=true
    refArrCPU[17,16,13]=true

    mainArrCPU[15,15,13]=true
    refArrCPU[15,16,13]=true

    mainArrCPU[16,15,13]=true
    refArrCPU[16,16,13]=true

    mainArrCPU[18,15,13]=true
    refArrCPU[18,16,13]=true

    mainArrCPU[19,15,13]=true
    refArrCPU[19,16,13]=true

    mainArrCPU[19,12,13]=true
    refArrCPU[19,12,13]=true

    mainArrCPU[12,12,17]=true
    refArrCPU[12,12,16]=true

    mainArrCPU[22,15,13]=true
    refArrCPU[18,22,13]=true

    mainArrCPU[19,22,13]=true
    refArrCPU[22,16,13]=true

    mainArrCPU[19,22,13]=true
    refArrCPU[19,12,22]=true

    mainArrCPU[12,12,17]=true
    refArrCPU[12,12,22]=true

    mainArrCPU[9,15,13]=true
    refArrCPU[18,9,13]=true

    mainArrCPU[19,22,9]=true
    refArrCPU[9,16,13]=true

    mainArrCPU[19,9,13]=true
    refArrCPU[19,12,9]=true

    mainArrCPU[9,12,17]=true
    refArrCPU[12,9,22]=true
end

mainArrGPU = CuArray(mainArrCPU);
## IMPORTANT in that way we set all points to agree so all should be aded to results !
# refArrGPU= CuArray(mainArrCPU);
refArrGPU= CUDA.ones(Bool,mainArrDims);
#setting offsets for all result queues
for qn in 1:14
    # resShmem[qn,1,1]=true
    #checking is it marked as to be validates
    # resShmem[qn,dataBdim[2],1]=true
    for a in [-1,0,1], b in [-1,0,1], c in [-1,0,1]
        metaData[xMeta+a+1,yMeta+b+1,zMeta+c+1,getResOffsetsBeg()+qn]=qn*50
    end
end

inBlockLoopX,inBlockLoopY,inBlockLoopZ= (fld(dataBdim[1] ,threads[1]),fld(dataBdim[2] ,threads[2]),dataBdim[3]    );
sourceShmem = CUDA.zeros(Bool,(dataBdim));
areToBeValidated= CUDA.ones(Bool,14)
isAnythingInPadding= CUDA.zeros(Bool,6)
function processPaddingKernel(isAnythingInPadding,areToBeValidated,sourceShmem,inBlockLoopX,inBlockLoopY,inBlockLoopZ,resShmem,metaData,resList,resListIndicies,maxResListIndex,metaDataDims,refArrGPU,mainArrDims,loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,dataBdim,mainArrGPU,iterNumb,isGold,xMeta,yMeta,zMeta,isToBeValidated)

    # paddingIter(loopX,loopY,maxXdim, maxYdim,a,b,c , xMetaChange,yMetaChange,zMetaChange, mainArr,refArr, dir,iterNumb,queueNumber,xMeta,yMeta,zMeta)
    @iterDataBlock(mainArrDims,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ,xMeta,yMeta,zMeta,
    begin
        ProcessMainDataVerB.@loadMainValues(mainArrGPU,dataBdim)
    end)

    sync_threads()

    @processPadding(isGold,xMeta,yMeta,zMeta,iterNumb,mainArrGPU,refArrGPU)

    return
end

@cuda threads=threads blocks=blocks processPaddingKernel(isAnythingInPadding,areToBeValidated,sourceShmem,inBlockLoopX,inBlockLoopY,inBlockLoopZ,resShmem,metaData,resList,resListIndicies,maxResListIndex,metaDataDims,refArrGPU,mainArrDims,loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,dataBdim,mainArrGPU,iterNumb,isGold,xMeta,yMeta,zMeta,isToBeValidated)

@test Int64(metaData[1,2,2,getIsToBeActivatedInGoldNumb()])==1
@test Int64(metaData[3,2,2,getIsToBeActivatedInGoldNumb()])==1
@test Int64(metaData[2,1,2,getIsToBeActivatedInGoldNumb()])==1
@test Int64(metaData[2,3,2,getIsToBeActivatedInGoldNumb()])==1
@test Int64(metaData[2,2,1,getIsToBeActivatedInGoldNumb()])==1
@test Int64(metaData[2,2,3,getIsToBeActivatedInGoldNumb()])==1

metaCorrX = 2
metaCorrY = 2
metaCorrZ = 2
### gold pass analysis
# 1)   Left FP  so this should be only to the right
@test Int64(metaData[metaCorrX+1,metaCorrY,metaCorrZ,getNewCountersBeg()+1])==2

# 3)   Right FP
@test Int64(metaData[metaCorrX-1,metaCorrY,metaCorrZ,getNewCountersBeg()+3])==1

# 5)   Posterior FP
@test Int64(metaData[metaCorrX,metaCorrY+1,metaCorrZ,getNewCountersBeg()+5])==4

# 7)   Anterior FP
@test Int64(metaData[metaCorrX,metaCorrY-1,metaCorrZ,getNewCountersBeg()+7])==3

# 9)   Top FP
@test Int64(metaData[metaCorrX,metaCorrY,metaCorrZ+1,getNewCountersBeg()+9])==6

# 11)   Bottom FP
@test Int64(metaData[metaCorrX,metaCorrY,metaCorrZ-1,getNewCountersBeg()+11])==5

# 13)   main block Fp
# @test Int64(metaData[metaCorrX,metaCorrY,metaCorrZ-1,getNewCountersBeg()+13])==8
@test length(filter(it->it>0,resListIndicies)) == (1+2+3+4+5+6)

@test sum(isAnythingInPadding)==6

######### process main part of data block
### data 1
using Revise, Parameters, Logging, Test
using CUDA
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils
using ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils,..ResultListUtils

threads=(32,5);
blocks =1;
mainArrDims= (67,177,90);
dataBdim = (10,10,10)
mainArrCPU= falses(mainArrDims);
mainArrCPU[5,5,5]= true;
mainArrCPU[33,10,7]= true;
refArrCPU = falses(mainArrDims);

metaData = MetaDataUtils.allocateMetadata(mainArrDims,dataBdim);
metaDataDims= size(metaData);
totalFp,totalFn= 500,500
resList,resListIndicies,maxResListIndex= allocateResultLists(totalFp,totalFn)
iterNumb=3
isGold = 1
queueNumber = 11
xMeta,yMeta,zMeta = 1,1,1 # here by assumption that we have all meta indicies 0 based
isToBeValidated = true
loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,inBlockLoopX,inBlockLoopY,inBlockLoopZ,metaDataLength,loopMeta,loopWarpMeta,clearIterResShmemLoop,clearIterSourceShmemLoop,clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength,sourceShmemTotalLength = calculateLoopsIter(dataBdim,threads[1],threads[2],metaDataDims,blocks)
resShmem = CUDA.zeros(Bool,(dataBdim[1]+2,dataBdim[2]+2,dataBdim[3]+2)) # we need this additional 33th an 34th spots
# resShmem[2,2,1]=true
# resShmem[3,2,1]=true
# resShmem[4,3,1]=true #for start we will get top padding
refArrCPU[dataBdim[1]+2,dataBdim[2]+2,dataBdim[3]+1]= true
#setting some trues in both
begin

    # so we should get all mins as 2 and all maxes as 5

    # 1)   Left FP
    # 2)   Left FN
    mainArrCPU[11,13,13]=true
    refArrCPU[11,13,15]=true

   # 3)   Right FP
   # 4)   Right FN
    mainArrCPU[20,13,13]=true
    refArrCPU[20,13,15]=true

    mainArrCPU[20,14,13]=true
    refArrCPU[20,14,15]=true

 # 5)   Posterior FP

    # 6)   Posterior FN

    mainArrCPU[13,11,14]=true
    refArrCPU[14,11,13]=true

    mainArrCPU[15,11,13]=true
    refArrCPU[13,11,15]=true

    mainArrCPU[16,11,13]=true
    refArrCPU[13,11,16]=true

    # 7)   Anterior FP
    # 8)   Anterior FN
    mainArrCPU[14,20,13]=true
    refArrCPU[13,20,14]=true

    mainArrCPU[15,20,13]=true
    refArrCPU[14,20,16]=true

    mainArrCPU[17,20,13]=true
    refArrCPU[15,20,17]=true

    mainArrCPU[18,20,17]=true
    refArrCPU[16,20,18]=true

    # 9)   Top FP
    # 10)   Top FN
    mainArrCPU[14,15,11]=true
    refArrCPU[14,16,11]=true

    mainArrCPU[17,15,11]=true
    refArrCPU[17,16,11]=true

    mainArrCPU[15,15,11]=true
    refArrCPU[15,16,11]=true

    mainArrCPU[16,15,11]=true
    refArrCPU[16,16,11]=true

    mainArrCPU[18,15,11]=true
    refArrCPU[18,16,11]=true

    # 11)   Bottom FP
    # 12)   Bottom FN

    mainArrCPU[14,15,20]=true
    refArrCPU[14,16,20]=true

    mainArrCPU[17,15,20]=true
    refArrCPU[17,16,20]=true

    mainArrCPU[15,15,20]=true
    refArrCPU[15,16,20]=true

    mainArrCPU[16,15,20]=true
    refArrCPU[16,16,20]=true

    mainArrCPU[18,15,20]=true
    refArrCPU[18,16,20]=true

    mainArrCPU[19,15,20]=true
    refArrCPU[19,16,20]=true

    # 13)   main block Fp
    # 14)   main block Fn

    mainArrCPU[14,15,14]=true
    refArrCPU[14,16,16]=true

    mainArrCPU[17,15,17]=true
    refArrCPU[17,16,13]=true

    mainArrCPU[15,15,13]=true
    refArrCPU[15,16,13]=true

    mainArrCPU[16,15,13]=true
    refArrCPU[16,16,13]=true

    mainArrCPU[18,15,13]=true
    refArrCPU[18,16,13]=true

    mainArrCPU[19,15,13]=true
    refArrCPU[19,16,13]=true

    mainArrCPU[19,12,13]=true
    refArrCPU[19,12,13]=true

    mainArrCPU[12,12,17]=true
    refArrCPU[12,12,16]=true

    mainArrCPU[22,15,13]=true
    refArrCPU[18,22,13]=true

    mainArrCPU[19,22,13]=true
    refArrCPU[22,16,13]=true

    mainArrCPU[19,22,13]=true
    refArrCPU[19,12,22]=true

    mainArrCPU[12,12,17]=true
    refArrCPU[12,12,22]=true

    mainArrCPU[9,15,13]=true
    refArrCPU[18,9,13]=true

    mainArrCPU[19,22,9]=true
    refArrCPU[9,16,13]=true

    mainArrCPU[19,9,13]=true
    refArrCPU[19,12,9]=true

    mainArrCPU[9,12,17]=true
    refArrCPU[12,9,22]=true
end

mainArrGPU = CuArray(mainArrCPU);
## IMPORTANT in that way we set all points to agree so all should be aded to results !
# refArrGPU= CuArray(mainArrCPU);
refArrGPU= CUDA.ones(Bool,mainArrDims);
#setting offsets for all result queues
for qn in 1:14
    # resShmem[qn,1,1]=true
    #checking is it marked as to be validates
    # resShmem[qn,dataBdim[2],1]=true
    for a in [-1,0,1], b in [-1,0,1], c in [-1,0,1]
        metaData[xMeta+a+1,yMeta+b+1,zMeta+c+1,getResOffsetsBeg()+qn]=qn*50
        metaData[xMeta+a+1,yMeta+b+1,zMeta+c+1,getIsToBeAnalyzedNumb()+qn]=1
    end
end

inBlockLoopX,inBlockLoopY,inBlockLoopZ= (fld(dataBdim[1] ,threads[1]),fld(dataBdim[2] ,threads[2]),dataBdim[3]    );
sourceShmem = CUDA.zeros(Bool,(dataBdim));
areToBeValidated= CUDA.ones(Bool,14)
isAnythingInPadding= CUDA.zeros(Bool,6)
function executeDataIterWithPaddingKernel(isAnythingInPadding,areToBeValidated,sourceShmem,inBlockLoopX,inBlockLoopY,inBlockLoopZ,resShmem,metaData,resList,resListIndicies,maxResListIndex,metaDataDims,refArrGPU,mainArrDims,loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,dataBdim,mainArrGPU,iterNumb,isGold,xMeta,yMeta,zMeta,isToBeValidated)
    isMaskFull= true
    @executeDataIterWithPadding(mainArrDims, inBlockLoopX,inBlockLoopY,inBlockLoopZ,mainArrGPU,refArrGPU,xMeta,yMeta,zMeta,isGold,iterNumb)

    return
end

@cuda threads=threads blocks=blocks executeDataIterWithPaddingKernel(isAnythingInPadding,areToBeValidated,sourceShmem,inBlockLoopX,inBlockLoopY,inBlockLoopZ,resShmem,metaData,resList,resListIndicies,maxResListIndex,metaDataDims,refArrGPU,mainArrDims,loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,dataBdim,mainArrGPU,iterNumb,isGold,xMeta,yMeta,zMeta,isToBeValidated)

@test Int64(metaData[1,2,2,getIsToBeActivatedInGoldNumb()])==1
@test Int64(metaData[3,2,2,getIsToBeActivatedInGoldNumb()])==1
@test Int64(metaData[2,1,2,getIsToBeActivatedInGoldNumb()])==1
@test Int64(metaData[2,3,2,getIsToBeActivatedInGoldNumb()])==1
@test Int64(metaData[2,2,1,getIsToBeActivatedInGoldNumb()])==1
@test Int64(metaData[2,2,3,getIsToBeActivatedInGoldNumb()])==1

metaCorrX = 2
metaCorrY = 2
metaCorrZ = 2
### gold pass analysis
# 1)   Left FP  so this should be only to the right
@test Int64(metaData[metaCorrX+1,metaCorrY,metaCorrZ,getNewCountersBeg()+1])==2

# 3)   Right FP
@test Int64(metaData[metaCorrX-1,metaCorrY,metaCorrZ,getNewCountersBeg()+3])==1

# 5)   Posterior FP
@test Int64(metaData[metaCorrX,metaCorrY+1,metaCorrZ,getNewCountersBeg()+5])==4

# 7)   Anterior FP
@test Int64(metaData[metaCorrX,metaCorrY-1,metaCorrZ,getNewCountersBeg()+7])==3

# 9)   Top FP
@test Int64(metaData[metaCorrX,metaCorrY,metaCorrZ+1,getNewCountersBeg()+9])==6

# 11)   Bottom FP
@test Int64(metaData[metaCorrX,metaCorrY,metaCorrZ-1,getNewCountersBeg()+11])==5

@test  metaData[metaCorrX,metaCorrY,metaCorrZ,getFullInSegmNumb()-1]==0

########## data 2

using Revise, Parameters, Logging, Test
using CUDA
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils
using ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils,..ResultListUtils

threads=(32,5);
blocks =1;
mainArrDims= (67,177,90);
dataBdim = (10,10,10)
mainArrCPU= falses(mainArrDims);

mainArrCPU[17,17,17]= true;

refArrCPU = falses(mainArrDims);

metaData = MetaDataUtils.allocateMetadata(mainArrDims,dataBdim);
metaDataDims= size(metaData);
totalFp,totalFn= 500,500
resList,resListIndicies,maxResListIndex= allocateResultLists(totalFp,totalFn)
iterNumb=3
isGold = 1
queueNumber = 11
xMeta,yMeta,zMeta = 1,1,1 # here by assumption that we have all meta indicies 0 based
isToBeValidated = true
loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,inBlockLoopX,inBlockLoopY,inBlockLoopZ,metaDataLength,loopMeta,loopWarpMeta,clearIterResShmemLoop,clearIterSourceShmemLoop,clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength,sourceShmemTotalLength = calculateLoopsIter(dataBdim,threads[1],threads[2],metaDataDims,blocks)
resShmem = CUDA.zeros(Bool,(dataBdim[1]+2,dataBdim[2]+2,dataBdim[3]+2)) # we need this additional 33th an 34th spots

mainArrGPU = CuArray(mainArrCPU);
## IMPORTANT in that way we set all points to agree so all should be aded to results !
# refArrGPU= CuArray(mainArrCPU);
refArrGPU= CUDA.ones(Bool,mainArrDims);
#setting offsets for all result queues
for qn in 1:14
    # resShmem[qn,1,1]=true
    #checking is it marked as to be validates
    # resShmem[qn,dataBdim[2],1]=true
    for a in [-1,0,1], b in [-1,0,1], c in [-1,0,1]
        metaData[xMeta+a+1,yMeta+b+1,zMeta+c+1,getResOffsetsBeg()+qn]=qn*50
        metaData[xMeta+a+1,yMeta+b+1,zMeta+c+1,getIsToBeAnalyzedNumb()+qn]=1
    end
end

inBlockLoopX,inBlockLoopY,inBlockLoopZ= (fld(dataBdim[1] ,threads[1]),fld(dataBdim[2] ,threads[2]),dataBdim[3]    );
sourceShmem = CUDA.zeros(Bool,(dataBdim));
areToBeValidated= CUDA.ones(Bool,14)
isAnythingInPadding= CUDA.zeros(Bool,6)
function executeDataIterWithPaddingKernel(isAnythingInPadding,areToBeValidated,sourceShmem,inBlockLoopX,inBlockLoopY,inBlockLoopZ,resShmem,metaData,resList,resListIndicies,maxResListIndex,metaDataDims,refArrGPU,mainArrDims,loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,dataBdim,mainArrGPU,iterNumb,isGold,xMeta,yMeta,zMeta,isToBeValidated)
    isMaskFull= true
    @executeDataIterWithPadding(mainArrDims, inBlockLoopX,inBlockLoopY,inBlockLoopZ,mainArrGPU,refArrGPU,xMeta,yMeta,zMeta,isGold,iterNumb)

    return
end

@cuda threads=threads blocks=blocks executeDataIterWithPaddingKernel(isAnythingInPadding,areToBeValidated,sourceShmem,inBlockLoopX,inBlockLoopY,inBlockLoopZ,resShmem,metaData,resList,resListIndicies,maxResListIndex,metaDataDims,refArrGPU,mainArrDims,loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,dataBdim,mainArrGPU,iterNumb,isGold,xMeta,yMeta,zMeta,isToBeValidated)

metaCorrX = 2
metaCorrY = 2
metaCorrZ = 2
@test  metaData[metaCorrX,metaCorrY,metaCorrZ,getFullInSegmNumb()-1]==0

# 13)   main block Fp
# @test Int64(metaData[metaCorrX,metaCorrY,metaCorrZ-1,getNewCountersBeg()+13])==8

### gold pass analysis
# 1)   Left FP  so this should be only to the right
@test Int64(metaData[metaCorrX,metaCorrY,metaCorrZ,getNewCountersBeg()+13])==6

length(filter(it->it>0,resListIndicies))

########## data 3

using Revise, Parameters, Logging, Test
using CUDA
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils
using ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils,..ResultListUtils

threads=(32,5);
blocks =1;
mainArrDims= (67,177,90);
dataBdim = (10,10,10)
mainArrCPU= trues(mainArrDims);

mainArrCPU[17,17,17]= true;

refArrCPU = falses(mainArrDims);

metaData = MetaDataUtils.allocateMetadata(mainArrDims,dataBdim);
metaDataDims= size(metaData);
totalFp,totalFn= 500,500
resList,resListIndicies,maxResListIndex= allocateResultLists(totalFp,totalFn)
iterNumb=3
isGold = 1
queueNumber = 11
xMeta,yMeta,zMeta = 1,1,1 # here by assumption that we have all meta indicies 0 based
isToBeValidated = true
loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,inBlockLoopX,inBlockLoopY,inBlockLoopZ,metaDataLength,loopMeta,loopWarpMeta,clearIterResShmemLoop,clearIterSourceShmemLoop,clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength,sourceShmemTotalLength = calculateLoopsIter(dataBdim,threads[1],threads[2],metaDataDims,blocks)
resShmem = CUDA.zeros(Bool,(dataBdim[1]+2,dataBdim[2]+2,dataBdim[3]+2)) # we need this additional 33th an 34th spots

mainArrGPU = CuArray(mainArrCPU);
## IMPORTANT in that way we set all points to agree so all should be aded to results !
# refArrGPU= CuArray(mainArrCPU);
refArrGPU= CUDA.ones(Bool,mainArrDims);
#setting offsets for all result queues
for qn in 1:14
    # resShmem[qn,1,1]=true
    #checking is it marked as to be validates
    # resShmem[qn,dataBdim[2],1]=true
    for a in [-1,0,1], b in [-1,0,1], c in [-1,0,1]
        metaData[xMeta+a+1,yMeta+b+1,zMeta+c+1,getResOffsetsBeg()+qn]=qn*50
        metaData[xMeta+a+1,yMeta+b+1,zMeta+c+1,getIsToBeAnalyzedNumb()+qn]=1
    end
end

inBlockLoopX,inBlockLoopY,inBlockLoopZ= (fld(dataBdim[1] ,threads[1]),fld(dataBdim[2] ,threads[2]),dataBdim[3]    );
sourceShmem = CUDA.zeros(Bool,(dataBdim));
areToBeValidated= CUDA.ones(Bool,14)
isAnythingInPadding= CUDA.zeros(Bool,6)
function executeDataIterWithPaddingKernel(isAnythingInPadding,areToBeValidated,sourceShmem,inBlockLoopX,inBlockLoopY,inBlockLoopZ,resShmem,metaData,resList,resListIndicies,maxResListIndex,metaDataDims,refArrGPU,mainArrDims,loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,dataBdim,mainArrGPU,iterNumb,isGold,xMeta,yMeta,zMeta,isToBeValidated)
    isMaskFull= true
    @executeDataIterWithPadding(mainArrDims, inBlockLoopX,inBlockLoopY,inBlockLoopZ,mainArrGPU,refArrGPU,xMeta,yMeta,zMeta,isGold,iterNumb)

    return
end

@cuda threads=threads blocks=blocks executeDataIterWithPaddingKernel(isAnythingInPadding,areToBeValidated,sourceShmem,inBlockLoopX,inBlockLoopY,inBlockLoopZ,resShmem,metaData,resList,resListIndicies,maxResListIndex,metaDataDims,refArrGPU,mainArrDims,loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,dataBdim,mainArrGPU,iterNumb,isGold,xMeta,yMeta,zMeta,isToBeValidated)

metaCorrX = 2
metaCorrY = 2
metaCorrZ = 2
@test  metaData[metaCorrX,metaCorrY,metaCorrZ,getFullInSegmNumb()-1]==1


text:
import SigDigits.count_digits
import SigDigits.count_decimalplaces
import SigDigits.count_nonzero_decimalplaces

@testset "Total digits" begin
    @test count_digits(0.0000007) == 8
    @test count_digits(5.034001) == 7
    @test count_digits(7000) == 4
    @test count_digits(1.000) == 1
    @test count_digits(6.983410) == 6
    @test count_digits(14080.093149) == 11
end

@testset "Decimal places" begin
    @test count_decimalplaces(0.0000007) == 7
    @test count_decimalplaces(5.034001) == 6
    @test count_decimalplaces(7000) == 0
    @test count_decimalplaces(1.000) == 0
    @test count_decimalplaces(6.983410) == 5
    @test count_decimalplaces(14080.093149) == 6
end

@testset "Nonzero decimal places" begin
    @test count_nonzero_decimalplaces(0.0000007) == 1
    @test count_nonzero_decimalplaces(5.034001) == 3
    @test count_nonzero_decimalplaces(7000) == 0
    @test count_nonzero_decimalplaces(1.000) == 0
    @test count_nonzero_decimalplaces(6.983410) == 5
    @test count_nonzero_decimalplaces(14080.093149) == 5
end


text:
const JD_PY_MIN_NAME = ".__py_tmp_minscript.py"

"""
$(SIGNATURES)

Does a full pass followed by a pre-rendering and minification step.

* `prerender=true`: whether to pre-render katex and highlight.js (requires `node.js`)
* `minify=true`:    whether to minify output (requires `python3` and `css_html_js_minify`)
* `sig=false`:      whether to return an integer indicating success (see [`publish`](@ref))
* `prepath=""`:     set this to something like "project-name" if it's a project page

Note: if the prerendering is set to `true`, the minification will take longer as the HTML files
will be larger (especially if you have lots of maths on pages).
"""
function optimize(; prerender::Bool=true, minify::Bool=true, sig::Bool=false,
                    prepath::String="", no_fail_prerender::Bool=true,
                    suppress_errors::Bool=true)::Union{Nothing,Bool}
    suppress_errors && (SUPPRESS_ERR[] = true)
    #
    # Prerendering
    #
    if prerender && !JD_CAN_PRERENDER
        @warn "I couldn't find node and so will not be able to pre-render javascript."
        prerender = false
    end
    if prerender && !JD_CAN_HIGHLIGHT
        @warn "I couldn't load 'highlight.js' so will not be able to pre-render code blocks. " *
              "You can install it with `npm install highlight.js`."
    end
    if !isempty(prepath)
        GLOBAL_PAGE_VARS["prepath"] = prepath => (String, )
    end
    # re-do a (silent) full pass
    start = time()
    fmsg = "\r→ Full pass"
    print(fmsg)
    withpre = fmsg * ifelse(prerender, rpad(" (with pre-rendering)", 24), rpad(" (no pre-rendering)", 24))
    print(withpre)
    succ = (serve(single=true, prerender=prerender, nomess=true, isoptim=true, no_fail_prerender=no_fail_prerender) === nothing)
    print_final(withpre, start)

    #
    # Minification
    #
    if minify && (succ || no_fail_prerender)
        if JD_CAN_MINIFY
            start = time()
            mmsg = rpad("→ Minifying *.[html|css] files...", 35)
            print(mmsg)
            # copy the script to the current dir
            cp(joinpath(dirname(pathof(JuDoc)), "scripts", "minify.py"), JD_PY_MIN_NAME; force=true)
            # run it
            succ = success(`$([e for e in split(PY)]) $JD_PY_MIN_NAME`)
            # remove the script file
            rm(JD_PY_MIN_NAME)
            print_final(mmsg, start)
        else
            @warn "I didn't find css_html_js_minify, you can install it via pip."*
                  "The output will not be minified."
        end
    end
    SUPPRESS_ERR[] = false
    return ifelse(sig, succ, nothing)
end

"""
$(SIGNATURES)

This is a simple wrapper doing a git commit and git push without much fanciness. It assumes the
current directory is a git folder.
This will work in most simple scenarios (e.g. there's only one person updating the website).
In other scenarios you should probably do this manually.

**Keyword arguments**

* `prerender=true`: prerender javascript before pushing see [`optimize`](@ref)
* `minify=true`:    minify output before pushing see [`optimize`](@ref)
* `nopass=false`:   set this to true if you have already run `optimize` manually.
* `prepath=""`:     set this to something like "project-name" if it's a project page
* `message="jd-update"`: add commit message.
"""
function publish(; prerender::Bool=true, minify::Bool=true, nopass::Bool=false,
                   prepath::String="", message::String="jd-update")::Nothing
    succ = true
    if !isempty(prepath) || !nopass
        succ = optimize(prerender=prerender, minify=minify, sig=true, prepath=prepath)
    end
    if succ
        start = time()
        pubmsg = rpad("→ Pushing updates with git...", 35)
        print(pubmsg)
        try
            run(`git add -A `)
            run(`git commit -m "$message" --quiet`)
            run(`git push --quiet`)
            print_final(pubmsg, start)
        catch e
            println("✘ Could not push updates, verify your connection and try manually.\n")
            @show e
        end
    else
        println("✘ Something went wrong in the optimisation step. Not pushing updates.")
    end
    return nothing
end

"""
$(SIGNATURES)

Cleanpull allows you to pull from your remote git repository after having removed the local
output directory. This will help avoid merge clashes.
"""
function cleanpull()::Nothing
    FOLDER_PATH[] = pwd()
    set_paths!()
    if isdir(PATHS[:pub])
        rmmsg = rpad("→ Removing local output dir...", 35)
        print(rmmsg)
        rm(PATHS[:pub], force=true, recursive=true)
        println("\r" * rmmsg * " [done ✔ ]")
    end
    try
        pmsg = rpad("→ Retrieving updates from the repository...", 35)
        print(pmsg)
        run(`git pull --quiet`)
        println("\r" * pmsg * " [done ✔ ]")
    catch e
        println("✘ Could not pull updates, verify your connection and try manually.\n")
        @show e
    end
    return nothing
end


text:
# If not integer...
function solve_ip!(m::JuMP.Model, d::StudentHousingData)
    for i = 1:nhouses
        setcategory(m[:λ][i], :Bin)
    end
    for λnew in λ_generated
        setcategory(λnew, :Bin)
    end
    @assert solve(m) == :Optimal
    m
end

function recover_integer_soln(m::JuMP.Model, d::StudentHousingData)
    solve_ip!(m, d)

    # In case one of our initial columns was in the optimal basis
    results = zeros(nhouses, length(d.patterns_allow))
    for k = 1:nhouses
        if getvalue(m[:λ][k]) ≈ 1.0
            results[k, StudentHousing.house_allowedby(k, d)[1]] = 1.0
        end
    end
    # Check all the other columns we added
    for k = 1:length(V_generated)
        if getvalue(λ_generated[k]) ≈ 1.0
            # Find the assignment that we picked for our house
            patterns = find(V_generated[k] .≈ 1.0)
            # Check the house index
            h = house_choice[k]
            # Mark the assignment
            results[h, patterns] .= 1.0
        end
    end
    results
end

# If not integer...
function recover_integer_soln_p(m::JuMP.Model, d::StudentHousingData)
    solve_ip!(m, d)
    npatterns = length(d.patterns_allow)

    results = zeros(nhouses, length(d.patterns_allow))
    for k = 1:npatterns
        if getvalue(m[:λ][k]) ≈ 1.0
            for i = 1:nhouses
                if house_fits_pattern(i, k, d)
                    results[i, k] .= 1.0
                    break
                end
            end
        end
    end
    # Check all the other columns we added
    for k = 1:length(U_generated)
        if getvalue(λ_generated[k]) ≈ 1.0
            # Find the assignment that we picked for our house
            houses = find(U_generated[k] .≈ 1.0)
            # Check the pattern index
            p = pattern_choice[k]
            # Mark the assignment
            results[houses, p] .= 1.0
        end
    end
    results
end


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

function write_reserve_margin_w(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
	T = inputs["T"]     # Number of time steps (hours)
	#dfResMar dataframe with weights included for calculations
	dfResMar_w = DataFrame(Constraint = [Symbol("t$t") for t in 1:T])
	temp_ResMar_w = transpose(dual.(EP[:cCapacityResMargin]))./inputs["omega"]
	if setup["ParameterScale"] == 1
		temp_ResMar_w = temp_ResMar_w * ModelScalingFactor # Convert from MillionUS$/GWh to US$/MWh
	end
	dfResMar_w = hcat(dfResMar_w, DataFrame(temp_ResMar_w, :auto))
	auxNew_Names_res=[Symbol("Constraint"); [Symbol("CapRes_$i") for i in 1:inputs["NCapacityReserveMargin"]]]
	rename!(dfResMar_w,auxNew_Names_res)
	CSV.write(joinpath(path, "ReserveMargin_w.csv"), dfResMar_w)
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
# JuLIP.jl master file.

module JuLIP

using Reexport
@reexport using NeighbourLists

# quickly switch between Matrices and Vectors of SVectors, etc
include("arrayconversions.jl")

# File IO
include("FIO.jl")
@reexport using JuLIP.FIO

# define types and abstractions of generic functions
include("abstractions.jl")

include("chemistry.jl")
@reexport using JuLIP.Chemistry

# the main atoms type
include("atoms.jl")
include("dofmanagement.jl")

# how to build some simple domains
include("build.jl")
@reexport using JuLIP.Build

# a few auxiliary routines
include("utils.jl")
@reexport using JuLIP.Utils

# interatomic potentials prototypes and some example implementations
include("Potentials.jl")
@reexport using JuLIP.Potentials

# basic preconditioning capabilities
include("preconditioners.jl")
@reexport using JuLIP.Preconditioners

# some solvers
include("Solve.jl")
@reexport using JuLIP.Solve

# experimental features
include("Experimental.jl")
@reexport using JuLIP.Experimental

# the following are some sub-modules that are primarily used
# to create further abstractions to be shared across several
# modules in the JuLIP-verse.
include("mlips.jl")
include("nbody.jl")

# codes to facilitate testing
include("Testing.jl")

end # module


text:
# probit
# https://stats.stackexchange.com/questions/155239/marginal-effects-from-bayesian-probit

using Turing
using Distributions
using Random
using LinearAlgebra
using StatsPlots
using MCMCChains
using DataFrames
using Printf
# using Plots # reesport
# using Statistics # reesport

using StatsFuns: normcdf, norminvcdf, normpdf

# cdf of normal distribution # Φ(x'β)
plot(normcdf.(Vector(-10:0.1:10)))
# probit function, inverse-cdf of normal distribution # Φ^{-1}(y)
plot(norminvcdf.(Vector(0:0.01:1)))

# geberate moc data from logit model
n = 1_000
num_var = 2
Random.seed!(99)
X = rand(Uniform(-sqrt(3.0), 2.0*sqrt(3.0)), n, num_var) # case-specific variables Z_i Matrix
X = [ones(n) X]
β = Vector([0.0, 0.5, -0.5]) # true coefficients
# V = X * β # representative utility, # utility (latent variable) y_star = U = Xβ + ϵ = V + ϵ
p = normcdf.(X * β) # probability (する確率)
y = rand.(Bernoulli.(p)) # observable choice 0 or 1

# set prior, β ~ MvNormal(β0, A0), no σ²
n, k = size(X)
β0 = zeros(k)
A0 = 0.01 * I

@model function probit_model(X, y)
    # Set the priors only for β.
    β ~ MvNormal(β0, inv(A0))

    # Write the likelihood, can't vectorize? no multivariate distribution
    n = size(X, 1)
    for i = 1:n
        y[i] ~ Bernoulli(normcdf(X[i, :]' * β))
    end
end

# Draw random number from posterior.
model = probit_model(X, y)

n_draws = 5_000
# n_chains = 4
n_chains = 1
n_tune = 1_000
chn = sample(model, NUTS(), MCMCThreads(), n_draws, n_chains)
# chn = sample(model, NUTS(), MCMCThreads(), n_draws, n_chains)[n_tune:n_draws, :, :]
β_chn = chn[n_tune+1:n_draws, :, :] # delete initial and exclude odd (hamiltonian bluh bluh) chain
β_chn = group(β_chn, :β)

summarize(β_chn)
quantile(β_chn)
plot(β_chn)

# Marginal Effect
"""
Marginal Effect of probit regression for continuous variable

∂Pr(Y=1|X)/∂X_j = ∂Pr(Φ(Xβ))/∂X_j
                = ϕ(Xβ)β_j, where ϕ() is standard normal pdf

The marginal effect on X_j depends on β_j, X_j and other β and X. Therefore,

    β_chn:          random drawed β
    β_interest:     variables of interest, :Symbol
                    one of a symbol in `chain` argument.
    x:              independent variable vector to calculate marginal
                    effect (individual, at mean, etc...).
    Output
    ------
    dydxs:   marginal effects regarding to each random drawed β, ((n_draws - n_tune) * n_chains) matrix
"""
function marginal_effect_probit(β_chn, β_interest, x)
    β_rand = get_params(β_chn).β # get random value
    β_rand_interest = β_chn[β_interest].data
    dydxs = []
    for which_chain = chains(β_chn)
        β_rand_single_chain = reduce(hcat, [β_rand[j].data[:, which_chain] for j=1:length(β_rand)])
        β_rand_interest_single_chain = β_rand_interest[:, which_chain]
        dydx = normpdf.(β_rand_single_chain * x) .* β_rand_interest_single_chain
        push!(dydxs, dydx)
    end
    dydxs = reduce(hcat, dydxs)
    return dydxs
end

rand_margins = marginal_effect_probit(β_chn, Symbol("β[2]"), X[100, :])
density(rand_margins)

"""
Average Marginal Effect of all sample.
"""
function marginal_effect_probit_mae(β_chn, β_interest, X)
    # https://discourse.julialang.org/t/converting-a-matrix-into-an-array-of-arrays/17038/7
    # https://discourse.julialang.org/t/how-to-broadcast-over-only-certain-function-arguments/19274
    return mean(marginal_effect_probit.((β_chn,), (β_interest,), [r for r in eachrow(X)]))
end

rand_mae = marginal_effect_probit_mae(β_chn, Symbol("β[2]"), X)
density(rand_mae)

# Plot Marginal Effect regarding to different value of X_j in interest
"""
Margins plot

Plot Marginal Effect regarding to different value of X_j in interest
# https://docs.julialang.org/en/v1/manual/functions/#Optional-Arguments
# https://docs.julialang.org/en/v1/manual/functions/#Keyword-Arguments
"""
function plot_marginal_effect(β_chn::Chains, β_interest::Symbol, x_interest_theoritical::Vector; at_mean::Bool=false, ci_p::Float64=0.95)
    interest_dim = findall(x -> x == β_interest, names(β_chn))
    other_dim = findall(x -> x != β_interest, names(β_chn))
    x_interest = X[:, interest_dim]
    # https://stackoverflow.com/questions/37661221/julia-select-all-but-one-element-in-array-matrix
    x_other = X[:, 1:end .!= interest_dim]
    # x_interest_theoritical = Vector(minimum(x_interest):step:maximum(x_interest))

    if at_mean
        # at mean other x
        X_theoritical_atmean = Matrix{Float64}(undef, size(x_interest_theoritical, 1), length(names(β_chn)))
        X_theoritical_atmean[:, interest_dim] = x_interest_theoritical
        for other in other_dim
            X_theoritical_atmean[:, other] .= mean(X, dims=1)[other]
        end
        rand_margins = marginal_effect_probit.((β_chn,), (β_interest,), [r for r in eachrow(X_theoritical_atmean)])
        ylabel = @sprintf("Marginal Effect at Means, CI(%d)", ci_p*100)
    else
        # Average other x
        X_theoritical_ames = []
        for x in x_interest_theoritical
            X_theoritical_ame = copy(X)
            X_theoritical_ame[:, interest_dim] .= x
            push!(X_theoritical_ames, X_theoritical_ame)
        end
        rand_margins = marginal_effect_probit_mae.((β_chn,), (β_interest,), X_theoritical_ames)
        ylabel = @sprintf("Average Marginal Effect, CI(%d)", ci_p*100)
    end

    # margins plot
    # https://discourse.julialang.org/t/asymmetric-error-bars-and-box-plots-in-julia/73647/2
    rand_margins_median = median.(rand_margins)
    ci_lower = quantile.(rand_margins, (1-ci_p)/2)
    ci_upper = quantile.(rand_margins, ci_p-(1-ci_p)/2)
    plot(x_interest_theoritical, rand_margins_median, color="black", xlabel=@sprintf("%s's x", β_interest), ylabel=ylabel,  legend=false)
    display(scatter!(x_interest_theoritical, rand_margins_median, yerror=(ci_lower, ci_upper), color="cornflowerblue"))

    return rand_margins
end

x_hypo = Vector(-2.0:0.5:3.5)

plot_marginal_effect(β_chn, Symbol("β[2]"), x_hypo, at_mean=false)
plot_marginal_effect(β_chn, Symbol("β[2]"), x_hypo, at_mean=true)
plot_marginal_effect(β_chn, Symbol("β[2]"), x_hypo, at_mean=true, ci_p=0.5)

plot_marginal_effect(β_chn, Symbol("β[3]"), x_hypo, at_mean=false)
plot_marginal_effect(β_chn, Symbol("β[3]"), x_hypo, at_mean=true)


text:
using EquivalentCircuits, DelimitedFiles, Test
using EquivalentCircuits: subcircuits, flatten, readablecircuit, simplifycircuit!, replace_redundant_cpes!
using EquivalentCircuits: isterminal, isoperation, generatekarva, karva_parameters,karva_to_tree, Circuit
using EquivalentCircuits: tree_to_function, circuitfunction

measurements = [5919.90 - 15.79im, 5919.58 - 32.68im, 5918.18 - 67.58im, 5912.24 - 139.49im,
 5887.12 - 285.74im, 5785.04 - 566.88im, 5428.94 - 997.19im, 4640.21 - 1257.83im, 3871.84 - 978.97im,
  3537.68 - 564.96im, 3442.94 - 315.40im, 3418.14 - 219.69im, 3405.51 - 242.57im, 3373.90 - 396.07im,
   3249.67 - 742.03im, 2808.42 - 1305.92im, 1779.41 - 1698.97im, 701.96 - 1361.47im, 208.29 - 777.65im, 65.93 - 392.51im];

frequencies = [0.10, 0.21, 0.43, 0.89, 1.83, 3.79, 7.85, 16.24, 33.60, 69.52, 143.84, 297.64, 615.85, 1274.27, 2636.65,
 5455.59, 11288.38, 23357.21, 48329.30, 100000.00];

library = [Circuit("+-R+C-RRCLCRCLLLR", Any[0, 0, 22.377192750081925, 0, 4.042529368461352e-9, 0, 3399.0384392519263, 2498.7384040467605, 4.0061261087766165e-6, 0.7948425599214597, 6.126315997167497e-5, 1373.5802658510875, 3.6191586501947424e-5, 0.17338854108986923, 0.46734655676131776, 0.24812830993871104, 1500.925039002129], 8.293095935567338e-6)
Circuit("+++--+RCRPRLRRCLL", Any[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 16.952525176348484, 4.019951476143508e-9, 3379.008162497164, [209604.83569111855, 0.9610007184206273], 2549.666570971039, 3.2866976510279723e-7, 2.418161729014113, 1654.6779402286252, 7.295488055251058e-5, 0.2994773175317633, 0.2905678390917996], 1.7509954992998087e-5)
Circuit("+-RRPCCRLLPPLLRCP", Any[0, 0, 0.00041729614400979383, 4595.3298212413, [6.1253600051612616e7, 0.8859174267613648], 8.481342404807135e-6, 8.005423544288912e-5, 1537.5259065364694, 0.7866046048091446, 0.4825771386985225, [3075.672528367832, 0.768918132091958], [421.20526643818135, 0.10530131660954534], 0.23435266803982535, 0.4782133397334183, 657.0140358658155, 5.755479554270154e-5, [524.214144301653, 0.13105353607541326]], 0.02705946970454961)];

library_fit = "R1-[C2,R3-[C4,R5]]";

head = 8;
encoding = generatekarva(head);
parameters = karva_parameters(encoding);
example_encoding = "+-+RRR-PCLLLCRPPP";
example_parameters = karva_parameters(encoding);
exam_params = [0.0, 0.0, 0.0, 1543.907617339479, 770.7386241387493, 1273.4835727121672, 0.0, [3499.829595663477, 0.8749573989158692],
 1.6880836678204815e-5, 0.04940667507884311, 0.6677953592773065, 0.42043991655286717, 9.192012270021564e-5, 3694.724527430516,
  [618.972522236354, 0.1547431305590885], [764.360425865072, 0.191090106466268], [3817.8583766637335, 0.9544645941659333]];
exam_params_redundant_CPE = [0.0, 0.0, 0.0, 1543.907617339479, 770.7386241387493, 1273.4835727121672, 0.0, [3499.829595663477, 1],
 1.6880836678204815e-5, 0.04940667507884311, 0.6677953592773065, 0.42043991655286717, 9.192012270021564e-5, 3694.724527430516,
  [618.972522236354, 0.1547431305590885], [764.360425865072, 0.191090106466268], [3817.8583766637335, 0.9544645941659333]];
example_tree = karva_to_tree(example_encoding,example_parameters);
example_circuit = Circuit(example_encoding,example_parameters,nothing);
example_circuit_redundant_CPE = Circuit(example_encoding,exam_params_redundant_CPE,nothing);
example_function = tree_to_function(example_tree);
Readable_example_circuit = readablecircuit(example_circuit);
example_function2 = circuitfunction(Readable_example_circuit);
circuit = Circuit(encoding,parameters,nothing);
readable_circuit = readablecircuit(circuit);
tree = karva_to_tree(encoding,parameters);

@testset "EquivalentCircuits.jl" begin
    # Check the circuitevolution function.
    @test library_fit == circuitevolution(measurements,frequencies,initial_population = library).Circuit
    # Evaluate properties of the generated circuit encoding.
    @test isoperation(encoding[1]) == true
    # Only terminals in the encoding's tail.
    @test all(isterminal.(collect(encoding[head+1:end]))) == true
    # conversion of circuit object to user-readable circuit.
    @test readablecircuit(example_circuit) == "[C1,P2]-R3-[R4,R5]"
    # Check the output of circuit simulation.
    @test example_function(exam_params,100) == 1789.4845259396134 - 10.86703970004648im
    # Replace redundant CPEs and simplify
    replace_redundant_cpes!(example_circuit_redundant_CPE) # converts the redundant CPE to an equivalent capacitor.
    simplifycircuit!(example_circuit_redundant_CPE) # simplifies all parralelly or serially connected similar components.
    # Circuit simplification should lead to a single resistor element R3.
    @test readablecircuit(example_circuit_redundant_CPE) == "C1-R2"
    # Parameteroptimisation : basic checks of solution lengths and bounds.
    optparams = parameteroptimisation("[C1,P2]-R3",measurements,frequencies)
    @test length(optparams) == 4
    C1,P2_1,P2_2,R3 = optparams
    @test 0<C1<10
    @test 0<P2_1<1.0e9
    @test 0<P2_2<1
    @test 0<R3<1.0e9
    # Subtrees length checking. There ought to be as many circuits as coding terminal elements.
    @test length(subcircuits(example_circuit)) == 5

end


text:
using LinearAlgebra
using StaticArrays

# node of a tree structured stack
mutable struct TSS{V <: AbstractVector{Int}}
	id::Int
	length::Int
	parent_id::Int
	parent_idx::Int # index of the parent at which stack joins to it
	chars::V
end

function TSS(id::Integer, parent_id::Integer, parent_idx::Integer, firstchars::Union{AbstractVector{T}, NTuple{N,T} where N} where {T<:Integer}, max_length::Integer)
	len = length(firstchars)
	chars = Vector{Int}(undef, max_length)
	chars[1:len] = collect(firstchars)
	return TSS(id, len, parent_id, parent_idx, chars)
end

TSS(id::Integer, firstchars::Union{AbstractVector{T}, NTuple{N,T} where N} where {T<:Integer}, max_length::Integer) =
TSS(id, 0, 0, firstchars, max_length)

function Base.push!(beam::TSS, char)
	beam.chars[beam.length += 1] = char
	return nothing
end

function Base.replace!(beam::TSS, parent_id::Integer, parent_idx::Integer, char::Int)
	beam.parent_id = parent_id
	beam.parent_idx = parent_idx
	beam.chars[1] = char
	beam.length = 1
	return nothing
end

@generated function getsequences(ps::StaticVector{N,<:StaticVector{D,T}}) where {T<:Real,D,N}
	head = quote
		scores = Vector{T}(undef, D^N)
		sequences = similar(scores, NTuple{N,Int})
		n = 1
	end

	indentation = "   "
	N_indentations = indentation^N
	strloop = string(
		"@inbounds ", # optional optimization
		ntuple(k -> indentation^(k-1) * "for (i$k, p$k) in enumerate(ps[$k])\n", N)...,
		N_indentations, "sequences[n] = (", ntuple(k -> "i$k, ", N)..., ")\n",
		N_indentations, "scores[n] = ", mapreduce(k -> "p$k", (ls, rs) -> ls * " * " * rs, 1:N), "\n",
		N_indentations, "n += 1\n",
		ntuple(k -> indentation^(N - k) * "end\n", N)...)

	# loop = Meta.quot(strloop)
	loop = Meta.parse(strloop)

	tail = :(return scores, sequences)
	expr = Expr(:block, head, loop, tail)
	return expr
end

function beam_search(ŷs::AbstractVector{<:AbstractVector}, width::Integer)::Vector{Vector{Int}}
	# alphabet length
	D = length(first(ŷs))
	# get integer initlen such that D^initlen >= beam width
	initlen = max(ceil(Int, log(D, width)), 1)
	# compute D^initlen number of candidate sequences and their corresponding scores
	scores, sequences = getsequences(SVector(ntuple(i -> Size(D)(ŷs[i]), initlen)))
	widthrng = 1:width
	# find the width number of top scoring candidate sequences
	topidxs = partialsortperm(scores, widthrng; rev=true)
	topscores = scores[topidxs]
	topsequences = sequences[topidxs]
	max_length = length(ŷs)
	# initialize beams with top candidate sequences
	beams = [ TSS(id, firstchars, max_length) for (id, (score, firstchars)) ∈ enumerate(zip(topscores, topsequences)) ]
	# preallocate arrays before main loop
	scores = Matrix{Float64}(undef, D, width)
	vecscores = vec(scores)
	idxs = collect(eachindex(vecscores))
	beamcounts = similar(widthrng)
	sortedbeamids = similar(beamcounts)
	@views ŷs = ŷs[(initlen+1):end]
	availablepairs = Vector{NTuple{2,Int}}(undef, width)

	for ŷ ∈ ŷs
		# compute the scores for all candidate sequences
		mul!(scores, ŷ, topscores')
		# pick the top scoring ones
		topidxs = partialsortperm!(idxs, vecscores, widthrng; rev=true)
		topscores = vecscores[topidxs]
		# each element is a (charidx, beamidx) tuple
		topsequences = CartesianIndices(scores)[topidxs]
		# how many times each beam is scored in the top
		fill!(beamcounts, 0)
		# flag indicating if the beam is unassigned
		unassigned = trues(width)
		i = 0
		for ci ∈ topsequences
			char, beamid = Tuple(ci)
			beamcounts[beamid] += 1
			if unassigned[beamid]
				# if the beam is not assigned, then assign it
				push!(beams[beamid], char)
				# and flag it as such
				unassigned[beamid] = false
			else # if assigned,
				# then add the pair to the list of available pairs
				availablepairs[i += 1] = (beamid, char)
			end
		end
		@views unassignedbeams = beams[unassigned]
		for (beam, (beamid, char)) ∈ zip(unassignedbeams, availablepairs)
			parentbeam = beams[beamid]
			parent_id = parentbeam.id
			parent_idx = parentbeam.length - 1
			replace!(beam, parent_id, parent_idx, char)
		end
	end

	predictions = Vector{Vector{Int}}(undef, width)
	for (i, beam) ∈ enumerate(beams)
		prediction = Vector{Int}(undef, max_length)
		Δ = max_length
		lastidx = beam.length
		while true
			rng = 1:lastidx
			Δ -= lastidx
			prediction[rng .+ Δ] = @view beam.chars[rng]
			parentbeam_id = beam.parent_id
			(parentbeam_id == 0) && break
			lastidx = beam.parent_idx
			beam = beams[parentbeam_id]
		end
		predictions[i] = prediction
	end

	return predictions
end

ŷs = [[0.1, 0.2, 0.3, 0.4, 0.5],
      [0.5, 0.4, 0.3, 0.2, 0.1],
      [0.1, 0.2, 0.3, 0.4, 0.5],
      [0.5, 0.4, 0.3, 0.2, 0.1],
      [0.1, 0.2, 0.3, 0.4, 0.5],
      [0.5, 0.4, 0.3, 0.2, 0.1],
      [0.1, 0.2, 0.3, 0.4, 0.5],
      [0.5, 0.4, 0.3, 0.2, 0.1],
      [0.1, 0.2, 0.3, 0.4, 0.5],
      [0.5, 0.4, 0.3, 0.2, 0.1]]
width = 3

predictions = beam_search(ŷs, width)

wan = ones(Int, length(first(predictions)))

paths = hcat((predictions .- (wan,))...)

truepaths = hcat([4, 0, 4, 0, 4, 0, 4, 0, 4, 0], [4, 0, 4, 0, 4, 0, 4, 0, 4, 1], [4, 0, 4, 0, 4, 0, 4, 0, 3, 0])


text:
module day14

struct Part1 end
struct Part2 end

include("inputs.jl")

ex1 = inputs.ex1_d14()
ex2 = inputs.ex2_d14()
input = inputs.input_d14()

function parse_input(input)
    lines = split(input, "\n")
    memrx = r"mem\[(\d+)\]"
    instructions = []
    for line in lines
        first, last= split(line, " = ")
        if startswith(line, "mem")
            idx = parse(Int, match(memrx, first)[1])
            val = parse(Int, last)
            msk = ""
        elseif startswith(line, "mask")
            msk = lpad(last, 64, '.') # make mask 64 bits for easy comparison
            idx, val = missing, missing
        end
        push!(instructions, (idx=idx, val=val, msk=msk))
    end
    return instructions
end

test_mask = "............................XXXXXXXXXXXXXXXXXXXXXXXXXXXXX1XXXX0X"
@assert parse_input(ex1)[1].msk == test_mask

# Part 1
function apply_mask(p::Part1, mask, val)
    result = Array{Char,1}()
    valbit = bitstring(UInt64(val))
    for bit = 29:64
        if mask[bit] == 'X'
            push!(result, valbit[bit])
        elseif mask[bit] == '1'
            push!(result, '1')
        elseif mask[bit] == '0'
            push!(result, '0')
        end
    end
    parse(Int, String(result), base=2)
end

@assert apply_mask(Part1(), test_mask, 11) == 73
@assert apply_mask(Part1(), test_mask, 101) == 101
@assert apply_mask(Part1(), test_mask, 0) == 64

write!(p::Part1, memory, mask, i) = memory[i.idx] = apply_mask(p, mask, i.val)

function initialize_program(p::Union{Part1,Part2}, input)
    insts = parse_input(input)
    memory = Dict{Int, Int}()
    mask = String(['X' for _ in 1:64])

    for i in insts
        if !isempty(i.msk)
            mask = i.msk
        else
            write!(p, memory, mask, i)
        end
    end
    return sum(values(memory))
end

part1(input) = initialize_program(Part1(), input)
@assert part1(ex1) == 165
@info "Day 14, Part 1 answer: $(part1(input))"

# Part 2
function apply_mask(p::Part2, mask, addr)
    result = Array{Char,1}()
    addrbit = bitstring(UInt64(addr))
    for bit = 29:64
        if mask[bit] == 'X'
            push!(result, 'X')
        elseif mask[bit] == '1'
            push!(result, '1')
        elseif mask[bit] == '0'
            push!(result, addrbit[bit])
        end
    end

    addrs = [result]
    while 'X' ∈ addrs[1]
        x = popfirst!(addrs)
        push!(addrs, replace(x, 'X'=>'0', count=1))
        push!(addrs, replace(x, 'X'=>'1', count=1))
    end
    [parse(Int, String(addr), base=2) for addr in addrs]
end

test_mask2 = lpad("000000000000000000000000000000X1001X", 64, '.')
@assert apply_mask(Part2(), test_mask2, 42) == [26, 27, 58, 59]

function write!(p::Part2, memory, mask, i)
    addrs = apply_mask(p, mask, i.idx)
    [memory[addr] = i.val for addr in addrs]
end

part2(input) = initialize_program(Part2(), input)
@assert part2(ex2) == 208
@info "Day 14, Part 2 answer: $(part2(input))"

end


text:
"""
    InfiniteOpt.add_point_variable(model::JuMP.Model,
                                   var::InfiniteOpt.PointVariable,
                                   key::Val{:TransData}
                                   )::InfiniteOpt.GeneralVariableRef

Make a `PointVariableRef` and map it to the appropriate transcription variable
and return the `GeneralVariableRef`. This is an extension of
[`add_point_variable`](@ref InfiniteOpt.add_point_variable(::JuMP.Model,::Any,::Any, ::Any))
for `TranscriptionOpt`.
"""
function InfiniteOpt.add_point_variable(
    model::JuMP.Model,
    ivref::InfiniteOpt.GeneralVariableRef,
    support::Vector{Float64},
    ::Val{:TransData}
    )::InfiniteOpt.GeneralVariableRef
    # check if an internal variable was already created
    data = transcription_data(model)
    internal_vref = get(data.point_lookup, (ivref, support), nothing)
    if !isnothing(internal_vref)
        return internal_vref
    end
    # check if whether the infinite model already has one, else create one
    inf_model = JuMP.owner_model(ivref)
    inf_model_index = get(inf_model.point_lookup, (ivref, support), nothing)
    if !isnothing(inf_model_index)
        return InfiniteOpt._make_variable_ref(inf_model, inf_model_index)
    else
        # make negative index to not conflict with the InfiniteModel
        raw_index = data.last_point_index -= 1
        # make the reference and map it to a transcription variable
        pvref = InfiniteOpt.GeneralVariableRef(JuMP.owner_model(ivref), raw_index,
                                               InfiniteOpt.PointVariableIndex)
        trans_var = lookup_by_support(model, ivref, support)
        data.finvar_mappings[pvref] = trans_var
        data.point_lookup[(ivref, support)] = pvref
        return pvref
    end
end

"""
    InfiniteOpt.add_semi_infinite_variable(model::JuMP.Model,
                                     var::InfiniteOpt.SemiInfiniteVariable,
                                     key::Val{:TransData}
                                     )::InfiniteOpt.GeneralVariableRef

Make a `SemiInfiniteVariableRef` and add `var` to the transcription data
and return the `GeneralVariableRef`. This is an extension of
[`add_semi_infinite_variable`](@ref InfiniteOpt.add_semi_infinite_variable(::JuMP.Model,::Any,::Any))
for `TranscriptionOpt`. Note that `internal_semi_infinite_variable` is also
extended to be able to access the `var`.
"""
function InfiniteOpt.add_semi_infinite_variable(
    model::JuMP.Model,
    var::InfiniteOpt.SemiInfiniteVariable,
    ::Val{:TransData}
    )::InfiniteOpt.GeneralVariableRef
    # check if an internal variable was already created
    ivref = var.infinite_variable_ref
    eval_supps = var.eval_supports
    data = transcription_data(model)
    internal_vref = get(data.semi_lookup, (ivref, eval_supps), nothing)
    if !isnothing(internal_vref)
        return internal_vref
    end
    # check if whether the infinite model already has one, else create one
    inf_model = JuMP.owner_model(ivref)
    inf_model_index = get(inf_model.semi_lookup, (ivref, eval_supps), nothing)
    if !isnothing(inf_model_index)
        return InfiniteOpt._make_variable_ref(inf_model, inf_model_index)
    else
        # make negative index to not conflict with the InfiniteModel
        semi_infinite_vars = transcription_data(model).semi_infinite_vars
        raw_index = -1 * (length(semi_infinite_vars) + 1)
        # make the reference and map it to a transcription variable
        rvref = InfiniteOpt.GeneralVariableRef(JuMP.owner_model(ivref), raw_index,
                                               InfiniteOpt.SemiInfiniteVariableIndex)
        push!(semi_infinite_vars, var)
        _set_semi_infinite_variable_mapping(model, var, rvref, InfiniteOpt._index_type(ivref))
        data.semi_lookup[(ivref, eval_supps)] = rvref
        return rvref
    end
end