text:
# molecular dynamics code (first attempt) - version 1
# Vinicius Piccoli - 29/04/2020

using ProgressMeter

  # pbc

  function pbcposition(s,L)
    if s > L
      s = s - L
    elseif s < 0
      s = s + L
    end
    return s
  end

 function  pbcseparation(ds,L)
   if ds > 0.5*L
     ds = ds - L
   elseif ds < -0.5*L
     ds = ds + L
   end
   return ds
 end

  #velocities

  function  setVelocities(N,intEner)
    vxSum = 0.0
    vySum = 0.0

    # arryas

    vx = zeros(N)
    vy = zeros(N)

    for i in 1:N
      vx[i] = rand() - 0.5 # vx
      vy[i] = rand() - 0.5 # vy
      vxSum = vxSum + vx[i]
      vySum = vySum + vy[i]
    end

    # zero center of mass momentum

    vxcm = vxSum / N
    vycm = vySum / N

    for i in 1:N
      vx[i] = vx[i] - vxcm
      vy[i] = vy[i] - vycm
    end

    #rescale velocities to a desired initial kinetic energy

    v2sum = 0

    for i in 1:N
      v2sum = v2sum + vx[i]^2 + vy[i]^2
    end
    KinectEnergyPerParticle = 0.5*v2sum/N
    rescale = sqrt(intEner/KinectEnergyPerParticle)
    for i in 1:N
      vx[i] = vx[i] * rescale
      vy[i] = vy[i] * rescale
    end
    return vx,vy
  end

# set ramdom variables

# não sei ainda fazer isso funcionar

#code to set random positions restraints
# If the system is a dilute gas, we can choose the initial positions of particles by placing them at random
# , making sure that no two particles are too close to another.

  function setRandomPos(N, Lx, Ly)

    rMinimumSquared  = 2.0 ^ (1.0 / 3.0)
    lx = zeros(N)
    ly = zeros(N)

    # first positions

    lx[1] = Lx*rand()
    ly[1] = Ly*rand()

    j = 1

    for i in 2:N

      lx[i] = Lx*rand() # x
      ly[i] = Ly*rand() # y

      overlap=true

      while j < i && overlap==true

       dx = pbcseparation(lx[i] - lx[j], Lx)
       dy = pbcseparation(ly[i] - ly[j], Ly)
       if dx*dx + dy*dy <  rMinimumSquared
         while dx*dx + dy*dy <  rMinimumSquared

           lx[i] = Lx*rand()
           ly[i] = Ly*rand()
           dx = pbcseparation(lx[i] - lx[j], Lx)
           dy = pbcseparation(ly[i] - ly[j], Ly)
           if dx*dx + dy*dy > rMinimumSquared
             overlap=false
           end

         end
       end

       j = j + 1

      end
    end

    return lx, ly

  end

# Placement of particles on rectangular lattice
# We must specify the number of particles per row nx and the number per collum ny.
# The linear dimensions Lx and Ly are adjustable parameters and can be varied after
# initialization to be as close as possible to their desired values.

  function setRectangular(Lx, Ly, nx, ny)

    dx = Lx / nx # distance between collumns
    dy = Ly / ny # distance between rows
    numx = zeros(nx*ny)
    numy = zeros(ny*nx)

    i = 0
    for ix in 1:nx    # loop through particles in an row
      for iy in 1:ny  # loop through rows

        i = i + 1
        numx[i] = dx * (ix + 0.5)
        numy[i] = dy * (iy + 0.5)

      end
    end
   return numx, numy
  end

  # Computing acceleration

  function accel(N, numx, numy, Lx, Ly)

    ax = zeros(N)
    ay = zeros(N)

    oneOverR2 = 0
    oneOverR6 = 0

    tpea = 0
    virialAc = 0

    #lennard-jones parameters -

    e = 0.00868
    sig = 3.4

    for i in 1:N
      for j in i+1:N

        dx = pbcseparation(numx[i] - numx[j], Lx)
        dy = pbcseparation(numy[i] - numy[j], Ly)
        r2 = dx*dx + dy*dy

       # fOver =  24 * (e/r2) * (2*(sig/r2)^12 - (sig/r2)^6)
       # we will truncate the force in r2, thefore particles that have
       # distance between each other larger than r2 will not have force computed.

        oneOverR2 = 1.0 / r2
        oneOverR6 = oneOverR2 * oneOverR2 * oneOverR2
        fOver = 48 * oneOverR6 * (oneOverR6 - 0.5) * oneOverR2

        fx = fOver * dx
        fy = fOver * dy

        ax[i] += fx    # a = f due to reduced units (mass is equal to unity)
        ay[i] += fy
        ax[j] -= fx
        ay[j] -= fy

        tpea += 4.0*(oneOverR6 * oneOverR6 - oneOverR6) # total potential energy accumulated
        virialAc += dx*fx + dy*fy # related to mean pressure

      end
    end
    return ax, ay, tpea, virialAc

  end

# COLOCAR EM UM SCRIPT SEPARADO DE FORMA QUE EU POSSA LER UM ARQUIVO EXTERNO
# CONTENDO OS PARAÂMETROS NECESSÁRIOS PARA A REALIZAÇÃO DA SIMULAÇÃO
# initializing molecular dynamics

let

  num = 2500

# timestep
  dt = 0.01
  dt2 = dt*dt

  time = zeros(num)
  time[1] = 0.

# box length
  Lx = 50
  Ly = sqrt(3) * Lx / 2

# number of molecules
  nx = 15
  ny = 15
  N  = nx * ny # number of molecules

# Initial kinect energy
  iner = 6

# Initial configuration
  posx, posy = setRectangular(Lx,Ly,nx,ny) # pos(100, N, Lx, Ly) # setRandomPos(N, Lx, Ly)

# Initial velocities
  vx, vy = setVelocities(N,iner)
  totalKinectenergy = 0.

# calculation of kinect energy
  for j in 1:length(vx)
    totalKinectenergy += (vx[j]*vx[j] + vy[j]*vy[j])
  end

# Vectos for potential and kinect energies and for the pressure

 ener = open("energies.dat","w")

 println(ener,"#time    potential    pressure    kinect   square_kin")

 PE = zeros(num)
 vir = zeros(num)
 PK =  zeros(num)
 PK_sqr = zeros(num) # square of PK

 PK[1] = totalKinectenergy * 0.5
 PK_sqr[1] = PK[1]*PK[1]

# Accelerations

 ax, ay, PE[1], vir[1]  = accel(N, posx, posy, Lx, Ly)
 println(ener,round(time[1],digits=3),"    ",round(PE[1],digits=3),"    ",round(vir[1],digits=3),"    ",round(PK[1],digits=3),"    ",round(PK_sqr[1],digits=3))

# file contain just the atoms positions
  file2 = open("pos_traj.dat","w")
  println(file2,"# posx    posy")
  println(file2,"# frame 1")

  for i in 1:length(posx)
    println(file2," ", posx[i],"    ",posy[i])
  end

# file -
  file = open("traj.xyz","w")

  println(file,"$(nx*ny)")
  println(file,"Frame 1")
  for i in 1:length(posx)
   println(file,"H","   ", posx[i],"   ",posy[i] ,"   " ,0.0)
 end

# verlet integrations - num steps

 q = 1

 @showprogress 1 "Running- MD of a 2D lennard-jones Gas" for k in 2:num
   q = q + 1

   time[q] = time[q-1] + dt

   # new positions with old accelerations

   totalKinectenergy = 0.

   for i in 1:N

     # POSITIONS

     posx[i] = posx[i] + vx[i]*dt + 0.5*ax[i]*dt2
     posy[i] = posy[i] + vy[i]*dt + 0.5*ay[i]*dt2

     # periodic boundary conditions

     posx[i] = pbcposition(posx[i],Lx)
     posy[i] = pbcposition(posy[i],Ly)

     # VELOCITIES

     vx[i]   = vx[i] + 0.5*ax[i]*dt
     vy[i]   = vy[i] + 0.5*ay[i]*dt

     totalKinectenergy += (vx[i]*vx[i] + vy[i]*vy[i])

  end

  PK[k] = 0.5*totalKinectenergy
  PK_sqr[k] = PK[k]*PK[k]

   # new accelerations with new positions and new velocities

  ax, ay, PE[k],vir[k] = accel(N, posx, posy, Lx, Ly)
  println(ener,round(time[k],digits=3),"    ",round(PE[k],digits=3),"    ",round(vir[k],digits=3),"    ",round(PK[k],digits=3),"    ",round(PK_sqr[k],digits=3))

   # New velocities

  for i in 1:N
    vx[i] = vx[i] + 0.5*ax[i]*dt
    vy[i] = vy[i] + 0.5*ay[i]*dt
  end

  # file for the calculation of g(r)
  println(file2,"# frame $k")
  for i in 1:length(posx)
    println(file2," ", posx[i],"    ",posy[i])
  end

   # vmd file
   println(file,"$(nx*ny)")
   println(file,"Frame $k")
   for i in 1:length(posx)
     println(file,"H","   ", posx[i],"   ",posy[i] ,"   " ,0.0)
   end

 end
 close(ener)
 close(file)
 close(file2)
end


text:
@testset "transform" begin
    @testset "add self-loops" begin
        A = [1  1  0  0
             0  0  1  0
             0  0  0  1
             1  0  0  0]
        A2 =   [2  1  0  0
                0  1  1  0
                0  0  1  1
                1  0  0  1]

        g = GNNGraph(A; graph_type=GRAPH_T)
        fg2 = add_self_loops(g)
        @test adjacency_matrix(g) == A
        @test g.num_edges == sum(A)
        @test adjacency_matrix(fg2) == A2
        @test fg2.num_edges == sum(A2)
    end

    @testset "batch"  begin
        g1 = GNNGraph(random_regular_graph(10,2), ndata=rand(16,10), graph_type=GRAPH_T)
        g2 = GNNGraph(random_regular_graph(4,2), ndata=rand(16,4), graph_type=GRAPH_T)
        g3 = GNNGraph(random_regular_graph(7,2), ndata=rand(16,7), graph_type=GRAPH_T)

        g12 = Flux.batch([g1, g2])
        g12b = blockdiag(g1, g2)
        @test g12 == g12b

        g123 = Flux.batch([g1, g2, g3])
        @test g123.graph_indicator == [fill(1, 10); fill(2, 4); fill(3, 7)]

        s, t = edge_index(g123)
        @test s == [edge_index(g1)[1]; 10 .+ edge_index(g2)[1]; 14 .+ edge_index(g3)[1]]
        @test t == [edge_index(g1)[2]; 10 .+ edge_index(g2)[2]; 14 .+ edge_index(g3)[2]]
        @test node_features(g123)[:,11:14] ≈ node_features(g2)

        # scalar graph features
        g1 = GNNGraph(g1, gdata=rand())
        g2 = GNNGraph(g2, gdata=rand())
        g3 = GNNGraph(g3, gdata=rand())
        g123 = Flux.batch([g1, g2, g3])
        @test g123.gdata.u == [g1.gdata.u, g2.gdata.u, g3.gdata.u]

        # Batch of batches
        g123123 = Flux.batch([g123, g123])
        @test g123123.graph_indicator == [fill(1, 10); fill(2, 4); fill(3, 7); fill(4, 10); fill(5, 4); fill(6, 7)]
        @test g123123.num_graphs == 6
    end

    @testset "unbatch" begin
        g1 = rand_graph(10, 20, graph_type=GRAPH_T)
        g2 = rand_graph(5, 10, graph_type=GRAPH_T)
        g12 = Flux.batch([g1, g2])
        gs = Flux.unbatch([g1,g2])
        @test length(gs) == 2
        @test gs[1].num_nodes == 10
        @test gs[1].num_edges == 20
        @test gs[1].num_graphs == 1
        @test gs[2].num_nodes == 5
        @test gs[2].num_edges == 10
        @test gs[2].num_graphs == 1
    end

    @testset "getgraph"  begin
        g1 = GNNGraph(random_regular_graph(10,2), ndata=rand(16,10), graph_type=GRAPH_T)
        g2 = GNNGraph(random_regular_graph(4,2), ndata=rand(16,4), graph_type=GRAPH_T)
        g3 = GNNGraph(random_regular_graph(7,2), ndata=rand(16,7), graph_type=GRAPH_T)
        g = Flux.batch([g1, g2, g3])

        g2b, nodemap = getgraph(g, 2, nmap=true)
        s, t = edge_index(g2b)
        @test s == edge_index(g2)[1]
        @test t == edge_index(g2)[2]
        @test node_features(g2b) ≈ node_features(g2)

        g2c = getgraph(g, 2)
        @test g2c isa GNNGraph{typeof(g.graph)}

        g1b, nodemap = getgraph(g1, 1, nmap=true)
        @test g1b === g1
        @test nodemap == 1:g1.num_nodes
    end

    @testset "add_edges" begin
        if GRAPH_T == :coo
            s = [1,1,2,3]
            t = [2,3,4,5]
            g = GNNGraph(s, t, graph_type=GRAPH_T)
            snew = [1]
            tnew = [4]
            gnew = add_edges(g, snew, tnew)
            @test gnew.num_edges == 5
            @test sort(inneighbors(gnew, 4)) == [1, 2]

            g = GNNGraph(s, t, edata=(e1=rand(2,4), e2=rand(3,4)), graph_type=GRAPH_T)
            # @test_throws ErrorException add_edges(g, snew, tnew)
            gnew = add_edges(g, snew, tnew, edata=(e1=ones(2,1), e2=zeros(3,1)))
            @test all(gnew.edata.e1[:,5] .== 1)
            @test all(gnew.edata.e2[:,5] .== 0)
        end
    end

    @testset "add_nodes" begin
        if GRAPH_T == :coo
            g = rand_graph(6, 4, ndata=rand(2, 6), graph_type=GRAPH_T)
            gnew = add_nodes(g, 5, ndata=ones(2, 5))
            @test gnew.num_nodes == g.num_nodes + 5
            @test gnew.num_edges == g.num_edges
            @test gnew.num_graphs == g.num_graphs
            @test all(gnew.ndata.x[:,7:11] .== 1)
        end
    end

    @testset "remove_self_loops" begin
        if GRAPH_T == :coo # add_edges and set_edge_weight only implemented for coo
            g = rand_graph(10, 20, graph_type=GRAPH_T)
            g1 = add_edges(g, [1:5;], [1:5;])
            @test g1.num_edges == g.num_edges + 5
            g2 = remove_self_loops(g1)
            @test g2.num_edges == g.num_edges
            @test sort_edge_index(edge_index(g2)) == sort_edge_index(edge_index(g))

            # with edge features and weights
            g1 = GNNGraph(g1, edata=(e1=ones(3,g1.num_edges), e2=2*ones(g1.num_edges)))
            g1 = set_edge_weight(g1, 3*ones(g1.num_edges))
            g2 = remove_self_loops(g1)
            @test g2.num_edges == g.num_edges
            @test sort_edge_index(edge_index(g2)) == sort_edge_index(edge_index(g))
            @test size(get_edge_weight(g2)) == (g2.num_edges,)
            @test size(g2.edata.e1) == (3, g2.num_edges)
            @test size(g2.edata.e2) == (g2.num_edges,)

        end
    end

    @testset "remove_multi_edges" begin
        if GRAPH_T == :coo
            g = rand_graph(10, 20, graph_type=GRAPH_T)
            s, t = edge_index(g)
            g1 = add_edges(g, s[1:5], t[1:5])
            @test g1.num_edges == g.num_edges + 5
            g2 = remove_multi_edges(g1, aggr=+)
            @test g2.num_edges == g.num_edges
            @test sort_edge_index(edge_index(g2)) == sort_edge_index(edge_index(g))

            # Default aggregation is +
            g1 = GNNGraph(g1, edata=(e1=ones(3,g1.num_edges), e2=2*ones(g1.num_edges)))
            g1 = set_edge_weight(g1, 3*ones(g1.num_edges))
            g2 = remove_multi_edges(g1)
            @test g2.num_edges == g.num_edges
            @test sort_edge_index(edge_index(g2)) == sort_edge_index(edge_index(g))
            @test count(g2.edata.e1[:,i] == 2*ones(3) for i in 1:g2.num_edges) == 5
            @test count(g2.edata.e2[i] == 4 for i in 1:g2.num_edges) == 5
            w2 = get_edge_weight(g2)
            @test count(w2[i] == 6 for i in 1:g2.num_edges) == 5
        end
    end

    @testset "negative_sample" begin
        if GRAPH_T == :coo
            n, m = 10, 30
            g = rand_graph(n, m, bidirected=true, graph_type=GRAPH_T)

            # check bidirected=is_bidirected(g) default
            gneg = negative_sample(g, num_neg_edges=20)
            @test gneg.num_nodes == g.num_nodes
            @test gneg.num_edges == 20
            @test is_bidirected(gneg)
            @test intersect(g, gneg).num_edges == 0
        end
    end

    @testset "rand_edge_split" begin
        if GRAPH_T == :coo
            n, m = 100,300

            g = rand_graph(n, m, bidirected=true, graph_type=GRAPH_T)
            # check bidirected=is_bidirected(g) default
            g1, g2 = rand_edge_split(g, 0.9)
            @test is_bidirected(g1)
            @test is_bidirected(g2)
            @test intersect(g1, g2).num_edges == 0
            @test g1.num_edges + g2.num_edges == g.num_edges
            @test g2.num_edges < 50

            g = rand_graph(n, m, bidirected=false, graph_type=GRAPH_T)
            # check bidirected=is_bidirected(g) default
            g1, g2 = rand_edge_split(g, 0.9)
            @test !is_bidirected(g1)
            @test !is_bidirected(g2)
            @test intersect(g1, g2).num_edges == 0
            @test g1.num_edges + g2.num_edges == g.num_edges
            @test g2.num_edges < 50

            g1, g2 = rand_edge_split(g, 0.9, bidirected=false)
            @test !is_bidirected(g1)
            @test !is_bidirected(g2)
            @test intersect(g1, g2).num_edges == 0
            @test g1.num_edges + g2.num_edges == g.num_edges
            @test g2.num_edges < 50
        end
    end

    @testset "set_edge_weight" begin
        g = rand_graph(10, 20, graph_type=GRAPH_T)
        w = rand(20)

        gw = set_edge_weight(g, w)
        @test get_edge_weight(gw) == w

        # now from weighted graph
        s, t = edge_index(g)
        g2 = GNNGraph(s, t, rand(20), graph_type=GRAPH_T)
        gw2 = set_edge_weight(g2, w)
        @test get_edge_weight(gw2) == w
    end

    @testset "to_bidirected" begin
        if GRAPH_T == :coo
            s, t = [1, 2, 3, 3, 4], [2, 3, 4, 4, 4]
            w = [1.0, 2.0, 3.0, 4.0, 5.0]
            e = [10.0, 20.0, 30.0, 40.0, 50.0]
            g = GNNGraph(s, t, w, edata = e)

            g2 = to_bidirected(g)
            @test g2.num_nodes == g.num_nodes
            @test g2.num_edges == 7
            s2, t2 = edge_index(g2)
            @test is_bidirected(g2)
            @test !has_multi_edges(g2)

            s2, t2 = edge_index(g2)
            w2 = get_edge_weight(g2)
            @test s2 == [1, 2, 2, 3, 3, 4, 4]
            @test t2 == [2, 1, 3, 2, 4, 3, 4]
            @test w2 == [1, 1, 2, 2, 3.5, 3.5, 5]
            @test g2.edata.e == [10.0, 10.0, 20.0, 20.0, 35.0, 35.0, 50.0]
        end
    end
end

ext: jl


text:
function modeorthonorm(basis,N,err = 1e-9,ω=1.)

x,w,T = nfieldtrans(N,2,ω=ω)
Normtest = T'*(T.*(w.*ones(N)))

#check orthonormality is close to identity matrix for the basis
maxerr = maximum(abs.(Normtest-eye(Normtest)))
maxerr < err
end


text:
# For low level cudnn functions that require a pointer to a number
cptr(x,a::DenseCuArray{Float64})=Float64[x]
cptr(x,a::DenseCuArray{Float32})=Float32[x]
cptr(x,a::DenseCuArray{Float16})=Float32[x]

# Conversion between Julia and CUDNN datatypes
cudnnDataType(::Type{Float16})=CUDNN_DATA_HALF
cudnnDataType(::Type{Float32})=CUDNN_DATA_FLOAT
cudnnDataType(::Type{Float64})=CUDNN_DATA_DOUBLE
cudnnDataType(::Type{Int8}) = CUDNN_DATA_INT8
cudnnDataType(::Type{UInt8}) = CUDNN_DATA_UINT8
cudnnDataType(::Type{Int32}) = CUDNN_DATA_INT32
# The following are 32-bit elements each composed of 4 8-bit integers, only supported with CUDNN_TENSOR_NCHW_VECT_C
# CUDNN_DATA_INT8x4,
# CUDNN_DATA_UINT8x4,
# CUDNN_DATA_INT8x32,
juliaDataType(a)=(a==CUDNN_DATA_HALF ? Float16 :
                  a==CUDNN_DATA_FLOAT ? Float32 :
                  a==CUDNN_DATA_DOUBLE ? Float64 :
                  a==CUDNN_DATA_INT8 ? Int8 :
                  a==CUDNN_DATA_UINT8 ? UInt8 :
                  a==CUDNN_DATA_INT32 ? Int32 : error())

tuple_strides(A::Tuple) = _strides((1,), A)
_strides(out::Tuple{Int}, A::Tuple{}) = ()
_strides(out::NTuple{N,Int}, A::NTuple{N}) where {N} = out
@inline function _strides(out::NTuple{M,Int}, A::Tuple) where M
    _strides((out..., out[M]*A[M]), A)
end

# The storage data types for alpha and beta are:
#     float for HALF and FLOAT tensors, and
#     double for DOUBLE tensors.
scalingParameter(T, val) = error("Unknown tensor type $T")
scalingParameter(::Type{Float16}, val) = Ref{Float32}(val)
scalingParameter(::Type{Float32}, val) = Ref{Float32}(val)
scalingParameter(::Type{Float64}, val) = Ref{Float64}(val)

# Create temporary reserveSpace. Use 128 to avoid alignment issues.
function cudnnTempSpace(nbytes)
    nbytes == 0 ? nothing : CuArray{Int128}(undef, (nbytes-1)÷sizeof(Int128)+1)
end


text:
include("LiftingDictionaries.jl")
include("LiftingProgressions.jl")
include("LiftingProgrammes.jl")

max_stars_repo_name: dcelisgarza/Lifting


text:
decode_id(c::IO) =
  let id = decode_int(c)
    if id == -1
      backend_error(c)
    else
      id
    end
  end

encode_BIMLevel = encode_int
decode_BIMLevel = decode_int_or_error
encode_FloorFamily = encode_int
decode_FloorFamily = decode_int_or_error

function request_operation(conn::IO, name)
  write(conn, Int32(0))
  encode_String(conn, name)
  op = read(conn, Int32)
  if op == -1
    error(name * " is not available")
  else
    op
  end
end

interrupt_processing(conn::IO) = write(conn, Int32(-1))
#=
create_op(name::ASCIIString, argtypes::Array{DataType}, rettype::DataType) = (
  op = request_operation(name);
  (args...) -> (conn = current_connection();
                write(conn, Int32(op));
                @assert length(args) == length(argtypes) "Incorrect number of args";
                for actual_arg in [argtype(arg) for (arg, argtype) in zip(args, argtypes)]
                  write(conn, actual_arg)
                end;
                read(conn, rettype)))

macro def_op(name, argtypes, rettype)
  :($(esc(name)) = create_op($(string(name)), $argtypes, $rettype))
end

@def_op(Circle, [XYZ, XYZ, Float64], Int32)
@def_op(Sphere, [XYZ, Float64], Int32)

circle = (op = request_operation("Circle");
          conn = current_connection();
          (c, n, r)-> (send_data(conn, Int8(op));
                       send_data(conn, c);
                       send_data(conn, n);
                       send_data(conn, Float64(r);
                       read(conn, Int32))))
=#

function parse_c_signature(sig)
  m = match(r"^ *(public|) *(\w+) *([\[\]]*) +((?:\w|:|<|>)+) *\( *(.*) *\)", sig)
  ret = Symbol(m.captures[2])
  array_level = count(c -> c=='[', something(m.captures[3], ""))
  name = Symbol(m.captures[4])
  params = split(m.captures[5], r" *, *", keepempty=false)
  function parse_c_decl(decl)
    m = match(r"^ *((?:\w|:|<|>)+) *([\[\]]*) *(\w+)$", decl)
    (Symbol(m.captures[1]), count(c -> c=='[', something(m.captures[2], "")), Symbol(m.captures[3]))
  end
  (name, [parse_c_decl(decl) for decl in params], (ret, array_level))
end

#=
const julia_type_for_c_type = Dict(
  :byte => :Int8,
  :double => :Float64,
  :float => :Float64,
  :int => :Int,
  :bool => :Bool,
  :Point3d => :XYZ,
  :Point2d => :XYZ,
  :Vector3d => :VXYZ,
  :string => :ASCIIString,
  :ObjectId => :Int32,
  :Entity => :Int32,
  :BIMLevel => :Int32,
  :FloorFamily => :Int32,
    ) #Either Int32 or Int64, depending on the architecture

julia_type(ctype, is_array) = is_array ? :(Vector{$(julia_type(ctype, false))}) : julia_type_for_c_type[ctype]
=#
export show_rpc, step_rpc
const show_rpc = Parameter(false)
const step_rpc = Parameter(false)
function initiate_rpc_call(conn, opcode, name)
    if step_rpc()
        print(stderr, "About to call $(name) [press ENTER]")
        readline()
    end
    if show_rpc()
        print(stderr, name)
    end
end
function complete_rpc_call(conn, opcode, result)
    if show_rpc()
        println(stderr, result == nothing ? "-> nothing" : "-> $(result)")
    end
    result
end

#=
This does not play well with reset_backend because the opcode are in global
scope and completely separated from the backend and, thus, the operation
still tries to use the same opcode even after a reset. Must think about
storing the generated functions on the backend itself, so that we can drop
it and restart.

Another option would be to have the table of generated methods separated
from the channel on the backend side, but I'm afraid that will prevent multiple
clients (not sure they are useful, though).
=#

function my_symbol(prefix, name)
  Symbol(prefix * replace(string(name), ":" => "_"))
end

function translate(sym)
  Symbol(replace(string(sym), r"[:<>]" => "_", ))
end

function rpc(prefix, str)
    name, params, ret = parse_c_signature(str)
    func_name = my_symbol(prefix, name)
    #Expr(:(=), Expr(:call, esc(name), [p[1] for p in params]...), Expr(:call, :+, params[1][1], 2))
    esc(quote
        $func_name =
          let opcode = -1
              buf = IOBuffer()
              (conn, $([:($(p[3])) for p in params]...)) -> begin
                opcode = Int32(request_operation(conn, $(string(name))))
                global $func_name = (conn, $([:($(p[3])) for p in params]...)) -> begin
                  initiate_rpc_call(conn, opcode, $(string(name)))
                  take!(buf) # Reset the buffer just in case there was an encoding error on a previous call
                  write(buf, opcode)
                  $([:($(Symbol("encode_", translate(p[1]), "_array"^p[2]))(buf, $(p[3]))) for p in params]...)
                  write(conn, take!(buf))
                  complete_rpc_call(conn, opcode, $(Symbol("decode_", translate(ret[1]), "_array"^ret[2]))(conn))
                end
                $func_name(conn, $([:($(p[3])) for p in params]...))
          end
        end
        #export $func_name
      end)
end

#=
We parameterize signatures to support different programming languages, e.g.
Signature{:CPP} is for C++ functions, while Signature{:CS} is for C# static methods
=#

struct Signature{T}
  description::AbstractString
end

# C++
parse_signature(sig::Signature{:CPP}) =
  let func_name(name) = replace(name, ":" => "_"),
      type_name(name) = replace(name, r"[:<>]" => "_"),
      m = match(r"^ *(public|) *(\w+) *([\[\]]*) +((?:\w|:|<|>)+) *\( *(.*) *\)", sig.description),
      ret = type_name(m.captures[2]),
      array_level = count(c -> c=='[', something(m.captures[3], "")),
      name = m.captures[4],
      params = split(m.captures[5], r" *, *", keepempty=false),
      parse_c_decl(decl) =
        let m = match(r"^ *((?:\w|:|<|>)+) *([\[\]]*) *(\w+)$", decl)
          (type_name(m.captures[1]), count(c -> c=='[', something(m.captures[2], "")), Symbol(m.captures[3]))
        end
    (func_name(name), name, [parse_c_decl(decl) for decl in params], (ret, array_level))
  end

# C#
parse_signature(sig::Signature{:CS}) =
  let m = match(r"^ *(public|) *(\w+) *([\[\]]*) +(\w+) *\( *(.*) *\)", sig.description),
      ret = m.captures[2],
      array_level = count(c -> c=='[', something(m.captures[3], "")),
      name = m.captures[4],
      params = split(m.captures[5], r" *, *", keepempty=false),
      parse_c_decl(decl) =
        let m = match(r"^ *(\w+) *([\[\]]*) *(\w+)$", decl)
          (m.captures[1], count(c -> c=='[', something(m.captures[2], "")), Symbol(m.captures[3]))
    end
    (name, name, [parse_c_decl(decl) for decl in params], (ret, array_level))
  end

parse_signature(sig::Signature{:Python}) =
  (:yeah, :whatever)

#=
A remote function encapsulates the information needed for communicating with remote
applications, including the opcode that represents the remote function and that
is generated by the remote application from the remote_name upon request.
=#

mutable struct RemoteFunction <: Function
  signature::Signature
  local_name::AbstractString
  remote_name::AbstractString
  opcode::Int32
  encoder::Function
  buffer::IOBuffer
end

remote_function(sig::Signature, local_name::AbstractString, remote_name::AbstractString, encoder::Function) =
  RemoteFunction(sig, local_name, remote_name, -1, encoder, IOBuffer())

ensure_opcode(f::RemoteFunction, conn) =
  f.opcode == -1 ?
    f.opcode = Int32(request_operation(conn, f.remote_name)) :
    f.opcode

reset_opcode(f::RemoteFunction) =
  f.opcode = -1

call_remote(f::RemoteFunction, conn, args...) =
  f.encoder(ensure_opcode(f, conn), conn, f.buffer, args...)

(f::RemoteFunction)(conn, args...) = call_remote(f, conn, args...)

#=
The most important part is the lang_rpc function. It parses the string describing the
signature of the remote function, generates a function to encode arguments and
retrieve results and, finally, creates the remote_function object.
=#

remote_function_meta_program(sig, local_name, remote_name, params, ret) =
    esc(:(remote_function(
           $(sig),
           $(local_name),
           $(remote_name),
           (opcode, conn, buf, $([:($(p[3])) for p in params]...)) -> begin
              initiate_rpc_call(conn, opcode, $(remote_name))
              take!(buf) # Reset the buffer just in case there was an encoding error on a previous call
              write(buf, opcode)
              $([:($(Symbol("encode_", p[1], "_array"^p[2]))(buf, $(p[3]))) for p in params]...)
              write(conn, take!(buf))
              complete_rpc_call(conn, opcode, $(Symbol("decode_", ret[1], "_array"^ret[2]))(conn))
            end)))

lang_rpc(lang, sigs) =
  [let sig = Signature{lang}(str),
      (local_name, remote_name, params, ret) = parse_signature(sig)
    (Symbol(local_name), remote_function_meta_program(sig, local_name, remote_name, params, ret))
   end
   for str in split(sigs, "\n") if str != ""]

#=
Given that remote apps might fail, it might be necessary to reset a connection.
This can be done either by resetting the opcodes of all remote functions or by
simply recreating all of them. This second alternative looks better because it
also allows to access multiple instances of the same remote app (as long as it)
support multiple separate connections.

=#

#=
The idea is that a particular remote application will store all of its functions
in a struct, as follows:

"""

@remote_functions :CS """
  int add(int a, int b)
  int sub(int a, int b)
"""
=#

export remote_functions
macro remote_functions(lang, str)
  let remotes = lang_rpc(lang.value, str)
    Expr(:tuple, [Expr(:(=), remote...) for remote in remotes]...)
  end
end


text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2017 Varnerlab
# Robert Frederick Smith School of Chemical and Biomolecular Engineering
# Cornell University, Ithaca NY 14850
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #
function FluxDriver(data_dictionary)
	# ----------------------------------------------------------------------------------- #
	# FluxDriver.jl was generated using the Kwatee code generation system.
	# FluxDriver: Solves the flux balance analysis problem from TSTART to TSTOP given the model encoded in data_dictionary.
	# Username: jeffreyvarner
	# Type: NFBA-JULIA
	# Version: 1.0
	# Generation timestamp: 03-03-2016 21:13:51
	#
	# Input arguments:
	# data_dictionary  - Data dictionary instance (holds species and flux models etc)
	#
	# Return arguments:
	# objective_value	 - value of the objective returned by GLPK
	# calculated_flux_array	 - optimal metabolic flux array (number_of_fluxes x 1)
	# uptake_array	 - stoichiometrix_matrix*flux_array (number_of_species x 1)
	# exit_flag	 - Exit flag returned by GLPK
	# ----------------------------------------------------------------------------------- #

	# Get the stoichiometric_matrix from data_dictionary -
	stoichiometric_matrix = data_dictionary["stoichiometric_matrix"]
	(number_of_species,number_of_fluxes) = size(stoichiometric_matrix);

	# Setup the GLPK problem -
	lp_problem = GLPK.Prob();
	GLPK.set_prob_name(lp_problem, "sample");
	GLPK.set_obj_name(lp_problem, "objective")

	# Set solver parameters
	solver_parameters = GLPK.SimplexParam();
	solver_parameters.msg_lev = GLPK.MSG_OFF;
	solver_parameters.presolve = GLPK.OFF;

	# Are we doing min -or- max?
	min_flag = data_dictionary["is_minimum_flag"];
	if min_flag == true
		GLPK.set_obj_dir(lp_problem, GLPK.MIN);
	else
		GLPK.set_obj_dir(lp_problem, GLPK.MAX);
	end

	# Set the number of constraints and fluxes -
	GLPK.add_rows(lp_problem, number_of_species);
	GLPK.add_cols(lp_problem, number_of_fluxes);

	# Setup flux bounds, and objective function -
	default_bounds_array = data_dictionary["default_flux_bounds_array"]
	(number_of_fluxes,number_of_bounds) = size(default_bounds_array)
	for flux_index = 1:number_of_fluxes

		flux_lower_bound = default_bounds_array[flux_index,1]
		flux_upper_bound = default_bounds_array[flux_index,2]

		# Check bounds type ... default is DB -
		if (flux_upper_bound == flux_lower_bound)
			flux_constraint_type = GLPK.FX
		else
			flux_constraint_type = GLPK.DB
		end

		# flux symbol? (later use name - for now, fake it)
		flux_symbol = "R_"*string(flux_index)

		# Set the bounds in GLPK -
		GLPK.set_col_name(lp_problem, flux_index, flux_symbol);
		GLPK.set_col_bnds(lp_problem, flux_index, flux_constraint_type, flux_lower_bound, flux_upper_bound);
	end

	# Setup objective function -
	objective_coefficient_array = data_dictionary["objective_coefficient_array"]
	for (flux_index,obj_coeff) in enumerate(objective_coefficient_array)

		# Set the objective function value in GLPK -
		GLPK.set_obj_coef(lp_problem, flux_index, obj_coeff);
	end

	# Setup problem constraints for the metabolites -
	species_bounds_array = data_dictionary["species_bounds_array"]
	for species_index = 1:number_of_species

		species_lower_bound = species_bounds_array[species_index,1]
		species_upper_bound = species_bounds_array[species_index,2]

		# defualt
		species_constraint_type = GLPK.FX
		if (species_lower_bound != species_upper_bound)
			species_constraint_type = GLPK.DB
		end

		# set the symbol -
		species_symbol = "x_"*string(species_index)

		# Set the species bounds in GLPK -
		GLPK.set_row_name(lp_problem, species_index, species_symbol);
		GLPK.set_row_bnds(lp_problem, species_index, species_constraint_type, species_lower_bound, species_upper_bound);

	end

	# Setup the stoichiometric array -
	counter = 1;
	row_index_array = zeros(Int,number_of_species*number_of_fluxes);
	col_index_array = zeros(Int,number_of_species*number_of_fluxes);
	species_index_vector = collect(1:number_of_species);
	flux_index_vector = collect(1:number_of_fluxes);
	flat_stoichiometric_array = zeros(Float64,number_of_species*number_of_fluxes);
	for species_index in species_index_vector
		for flux_index in flux_index_vector
			row_index_array[counter] = species_index;
			col_index_array[counter] = flux_index;
			flat_stoichiometric_array[counter] = stoichiometric_matrix[species_index,flux_index];
			counter+=1;
		end
	end
	GLPK.load_matrix(lp_problem, number_of_species*number_of_fluxes, row_index_array, col_index_array, flat_stoichiometric_array);

	# Call the solver -
	exit_flag = GLPK.simplex(lp_problem, solver_parameters);

	# Get the objective function value -
	objective_value = GLPK.get_obj_val(lp_problem);

	# Get the calculated flux values from GLPK -
	calculated_flux_array = zeros(Float64,number_of_fluxes);
	for flux_index in flux_index_vector
		calculated_flux_array[flux_index] = GLPK.get_col_prim(lp_problem, flux_index);
	end

	# Get the dual values -
	dual_value_array = zeros(Float64,number_of_fluxes);
	for flux_index in flux_index_vector
		dual_value_array[flux_index] = GLPK.get_col_dual(lp_problem, flux_index);
	end

	# is this solution optimal?
	status_flag = GLPK.get_status(lp_problem)

	# Calculate the uptake array -
	uptake_array = stoichiometric_matrix*calculated_flux_array;

	# Formulate the return tuple -
	return (objective_value, calculated_flux_array, dual_value_array, uptake_array, exit_flag,status_flag);
end


text:
config_path(args...) = joinpath(abspath(first(DEPOT_PATH), "config", "DFControl"), args...)

# All the folder
mkpath(config_path())
mkpath(config_path("servers"))
mkpath(config_path("execs"))
mkpath(config_path("jobs"))
mkpath(config_path("logs"))
mkpath(config_path("logs/jobs"))

# All the job info
touch(config_path("jobs", "pending.txt"))
touch(config_path("jobs", "archived.txt"))
touch(config_path("jobs", "active.txt"))
touch(config_path("jobs", "running.txt"))


text:
"""
histogram(array; bins = 100)
# Description
return arrays for plotting histogram
"""
function histogram(
    array;
    bins = minimum([100, length(array)]),
    normalize = true,
)
    tmp = zeros(bins)
    down, up = extrema(array)
    down, up = down == up ? (down - 1, up + 1) : (down, up) # edge case
    bucket = collect(range(down, up, length = bins + 1))
    normalization = normalize ? length(array) : 1
    for i in eachindex(array)
        # normalize then multiply by bins
        val = (array[i] - down) / (up - down) * bins
        ind = ceil(Int, val)
        # handle edge cases
        ind = maximum([ind, 1])
        ind = minimum([ind, bins])
        tmp[ind] += 1 / normalization
    end
    return (bucket[2:end] + bucket[1:(end - 1)]) .* 0.5, tmp
end

function histogram_plot(array; bins = 100)
    vxs, vys = histogram(array, bins = bins)

    pdf = barplot(
        vxs,
        vys,
        color = :red,
        strokecolor = :red,
        strokewidth = 1,
    )
end


text:
module HotFuzz

using MacroTools

include("hooks.jl")
include("trace.jl")

end # module


text:
"""

Julia  の HypothesisTests には，母相関係数の検定がないので，以下の関数を定義する。

スピアマンの順位相関係数，ケンドールの順位相関係数においては，同順位がある場合には近似検定になる。

ケンドールの順位相関係数においては，信頼区間は求めない。
"""

using Statistics, StatsBase, Distributions
roundn(x, digits) = round(x, digits=digits)
formatp(p) = p < 0.00001 ? "< 0.00001" : "= " * string(roundn(p, 5))
index(method) = indexin(method[1], ['p', 's', 'k'])[1]
function CorTest(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}; level=0.95, method="pearson")
    func = [cor, corspearman, corkendall][index(method)]
    CorTest(length(x), func(x, y); level, method)
end
function CorTest(n::Int, r::Float64; level=0.95, method="pearson")
    if method != "kendall"
        t = r * sqrt((n-2)/(1-r^2))
        df = n-2
        p = 2ccdf(TDist(df), abs(t))
        q = quantile(Normal(), 0.5 - 0.5level)
        conf = tanh.(atanh(r) .+ [q, -q] ./ sqrt(n-3))
        println("n = $n, r = $(roundn(r, 5)), t = $(roundn(t, 5)), df = $df, p-value $(formatp(p))")
        println("$(Int(100level)) percent confidence interval = $(roundn.(conf, 5))")
    else
        z = r / sqrt((4n + 10) / (9n * (n - 1)))
        p = 2ccdf(Normal(), abs(z))
        println("n = $n, r = $(roundn(r, 5)), z = $(roundn(z, 5)), p-value $(formatp(p))")
    end
end;

"""
使用法

julia> x = [58.3, 43.6, 49.9, 64.1, 51.6, 36.6, 58.7, 66.7, 50.7, 53.1,
            56.1, 42.0, 45.5, 30.9, 42.3];

julia> y = [46.3, 51.5, 64.5, 55.6, 63.9, 35.6, 61.2, 56.9, 44.2, 57.3,
            42.7, 47.1, 32.4, 39.7, 51.0];

julia> CorTest(x, y)
n = 15, r = 0.52225, t = 2.20803, df = 13, p-value = 0.04582
95 percent confidence interval = [0.01363, 0.81616]

julia> CorTest(x, y, method="spearman")
n = 15, r = 0.47143, t = 1.92737, df = 13, p-value = 0.07607
95 percent confidence interval = [-0.05384, 0.79234]

julia> CorTest(x, y, method="kendall")
n = 15, r = 0.31429, z = 1.63308, p-value = 0.10245

julia> CorTest(x, y, level=0.99) # 信頼率の設定は，level= で指定する
n = 15, r = 0.52225, t = 2.20803, df = 13, p-value = 0.04582
99 percent confidence interval = [-0.16269, 0.86753]

julia> CorTest(24, 0.476, level=0.9) # サンプルサイズと相関係数を与える場合
n = 24, r = 0.476, t = 2.53869, df = 22, p-value = 0.01871
90 percent confidence interval = [0.15754, 0.70478]

"""

max_stars_repo_name: r-de-r/stats

max_issues_repo_name: r-de-r/stats

max_forks_repo_name: r-de-r/stats


text:
const ERRMSG_HIDE_COMMUNICATION = "@hide_communication: the code block must start with exactly one @parallel computation call and be followed by code to set boundary conditions and to perform communication. Type `?@hide_communication` for more information."
const ERRMSG_INVALID_STREAM = "@hide_communication: invalid keyword argument in @parallel: stream. Specifying the stream is not allowed in a @parallel call inside the code block argument of @hide_communication."
const ERRMSG_INVALID_COMP = "@hide_communication: invalid positional arguments in @parallel: no positional argument is allowed in the @parallel call to perform computations inside the code block argument of @hide_communication."
const ERRMSG_INVALID_BC_COMM = "@hide_communication: invalid positional arguments in @parallel: exactly one positional argument argument (ranges) must be passed in @parallel calls to set boundary conditions and to perform communication inside the code block argument of @hide_communication."

##
const HIDE_COMMUNICATION_DOC = """
    @hide_communication boundary_width block

!!! note "Advanced"
        @hide_communication ranges_outer ranges_inner block

Hide the communication behind the computation within the code `block`.

# Arguments
- `boundary_width::Tuple{Integer,Integer,Integer} | Tuple{Integer,Integer} | Tuple{Integer}`: width of the boundaries in each dimension. The boundaries must include (at least) all the data that is accessed in the communcation performed.
- `block`: code block wich starts with exactly one [`@parallel`](@ref) call to perform computations, followed by code to set boundary conditions and to perform communication (as e.g. `update_halo!` from the package `ImplicitGlobalGrid`). The [`@parallel`](@ref) call to perform computations cannot contain any positional arguments (ranges, nblocks or nthreads) nor the stream keyword argument (stream=...). The code to set boundary conditions and to perform communication must only access the elements in the boundary ranges of the fields modified in the [`@parallel`](@ref) call; all elements can be acccessed from other fields. Moreover, this code must not include statements in array broadcasting notation, because they are always run on the default CUDA stream (for CUDA.jl < v2.0), which makes CUDA stream overlapping impossible. Instead, boundary region elements can, e.g., be accessed with [`@parallel`](@ref) calls passing a ranges argument that ensures that no threads mapping to elements outside of `ranges_outer` are launched. Note that these [`@parallel`](@ref) `ranges` calls cannot contain any other positional arguments (nblocks or nthreads) nor the stream keyword argument (stream=...).

!!! note "Advanced"
    - `ranges_outer::`Tuple with one or multiple `ranges` as required by the corresponding argument of [`@parallel`](@ref): the `ranges` must together span (at least) all the data that is accessed in the communcation and boundary conditions performed.
    - `ranges_inner::`Tuple with one or multiple `ranges` as required by the corresponding argument of [`@parallel`](@ref): the `ranges` must together span the data that is not included by `ranges_outer`.

# Examples
    @hide_communication (16, 2, 2) begin
        @parallel diffusion3D_step!(Te2, Te, Ci, lam, dt, dx, dy, dz);
        update_halo!(Te2);
    end

    @hide_communication (16, 2) begin
        @parallel diffusion2D_step!(Te2, Te, Ci, lam, dt, dx, dy);
        update_halo!(Te2);
    end

    @hide_communication ranges_outer ranges_inner begin
        @parallel diffusion3D_step!(Te2, Te, Ci, lam, dt, dx, dy, dz);
        update_halo!(Te2);
    end

    @parallel_indices (iy,iz) function bc_x(A)
        A[  1, iy,  iz] = A[    2,   iy,   iz]
        A[end, iy,  iz] = A[end-1,   iy,   iz]
        return
    end
    @parallel_indices (ix,iz) function bc_y(A)
        A[ ix,  1,  iz] = A[   ix,    2,   iz]
        A[ ix,end,  iz] = A[   ix,end-1,   iz]
        return
    end
    @parallel_indices (ix,iy) function bc_z(A)
        A[ ix,  iy,  1] = A[   ix,   iy,    2]
        A[ ix,  iy,end] = A[   ix,   iy,end-1]
        return
    end
    @hide_communication (16, 2, 2) begin
        @parallel diffusion3D_step!(Te2, Te, Ci, lam, dt, dx, dy, dz);
        @parallel (1:size(Te,2), 1:size(Te,3)) bc_x(Te);
        @parallel (1:size(Te,1), 1:size(Te,3)) bc_y(Te);
        @parallel (1:size(Te,1), 1:size(Te,2)) bc_z(Te);
        update_halo!(Te2);
    end

!!! note "Developers note"
    The communcation should not perform any blocking operations to enable a maximal overlap of communication with computation.

See also: [`@parallel`](@ref)
"""
@doc HIDE_COMMUNICATION_DOC
macro hide_communication(args...) check_initialized(); checkargs_hide_communication(args...); esc(hide_communication(args...)); end

## ARGUMENT CHECKS

function checkargs_hide_communication(args...)
    if (length(args) < 2 || length(args) > 3) @ArgumentError("wrong number of arguments.") end
    if !is_block(args[end]) @ArgumentError("the last argument must be a code block (obtained: $(args[end])).") end
end

## GATEWAY FUNCTIONS

function hide_communication(args::Union{Integer,Symbol,Expr}...; package::Symbol=get_package())
    if     (package == PKG_CUDA)    hide_communication_cuda(args...)
    elseif (package == PKG_THREADS) hide_communication_threads(args...)
    else                            @KeywordArgumentError("$ERRMSG_UNSUPPORTED_PACKAGE (obtained: $package).")
    end
end

## @HIDE_COMMUNICATION FUNCTIONS

function hide_communication_cuda(ranges_outer::Union{Symbol,Expr}, ranges_inner::Union{Symbol,Expr}, block::Expr)
    compcall, bc_and_commcalls = extract_calls(block)
    parallel_args = extract_args(compcall, Symbol("@parallel"))
    posargs, kwargs, compkernelcall = split_parallel_args(parallel_args)
    if any([x.args[1]==:stream for x in kwargs]) @ArgumentError(ERRMSG_INVALID_STREAM) end
    if (length(posargs) > 0) @ArgumentError(ERRMSG_INVALID_COMP) end
    bc_and_commcalls = process_bc_and_commcalls(bc_and_commcalls)
    quote
        for i in 1:length($ranges_outer)
            ParallelStencil.ParallelKernel.@parallel_async $ranges_outer[i] stream=ParallelStencil.ParallelKernel.get_priority_custream(i) $compkernelcall
        end
        for i in 1:length($ranges_inner)
            ParallelStencil.ParallelKernel.@parallel_async $ranges_inner[i] stream=ParallelStencil.ParallelKernel.get_custream(i) $compkernelcall
        end
        for i in 1:length($ranges_outer) ParallelStencil.ParallelKernel.@synchronize(ParallelStencil.ParallelKernel.get_priority_custream(i)); end
        $bc_and_commcalls
        for i in 1:length($ranges_inner) ParallelStencil.ParallelKernel.@synchronize(ParallelStencil.ParallelKernel.get_custream(i)); end
    end
end

function hide_communication_threads(ranges_outer::Union{Symbol,Expr}, ranges_inner::Union{Symbol,Expr}, block::Expr)
    return block #NOTE: This is currently not implemented, but enables correct execution nevertheless.
end

function hide_communication_cuda(boundary_width::Union{Integer,Symbol,Expr}, block::Expr)
    compcall, = extract_calls(block)
    parallel_args = extract_args(compcall, Symbol("@parallel"))
    posargs, kwargs, compkernelcall = split_parallel_args(parallel_args)
    if any([x.args[1]==:stream for x in kwargs]) @ArgumentError(ERRMSG_INVALID_STREAM) end
    if (length(posargs) > 0) @ArgumentError(ERRMSG_INVALID_COMP) end
    ranges = :(ParallelStencil.ParallelKernel.get_ranges($(compkernelcall.args[2:end]...)))
    ranges_outer = gensym_world("ranges_outer", @__MODULE__)
    ranges_inner = gensym_world("ranges_inner", @__MODULE__)
    quote
        $ranges_outer = ParallelStencil.ParallelKernel.get_ranges_outer($boundary_width, $ranges)
        $ranges_inner = ParallelStencil.ParallelKernel.get_ranges_inner($boundary_width, $ranges)
        ParallelStencil.ParallelKernel.@hide_communication $ranges_outer $ranges_inner $block  # NOTE: to enable macros like @hide_communication_cuda, one would need to use everywhere in hide_communication_cuda _cuda and _thread suffixed macros and enable their processing. This is more complex and generates less readable code though. => Design decision: only enable _cuda/_threads macros for macros that do not return package dependent macros.
    end
end

function hide_communication_threads(boundary_width::Union{Integer,Symbol,Expr}, block::Expr)
    return block #NOTE: This is currently not implemented, but enables correct execution nevertheless.
end

## FUNCTIONS TO EXTRACT AND PROCESS COMPUTATION AND BOUNDARY CONDITIONS CALLS / COMMUNICATION CALLS

function extract_calls(block::Expr)
    if (sum([isa(arg, Expr) for arg in block.args]) < 2) @ArgumentError(ERRMSG_HIDE_COMMUNICATION) end
    compcall         = block.args[2]
    bc_and_commcalls = quote $(block.args[3:end]...) end
    if !is_parallel_call(compcall) @ArgumentError(ERRMSG_HIDE_COMMUNICATION) end
    return compcall, bc_and_commcalls
end

function process_bc_and_commcalls(block::Expr)
    return postwalk(block) do x
        if !is_parallel_call(x) return x; end
        @capture(x, @parallel args__) || @ModuleInternalError("a @parallel call could not be pattern matched.")
        posargs, kwargs = split_parallel_args(args)
        if any([x.args[1]==:stream for x in kwargs]) @ArgumentError(ERRMSG_INVALID_STREAM) end
        if (length(posargs) != 1) @ArgumentError(ERRMSG_INVALID_BC_COMM) end
        return :(@parallel $(args[1:end-1]...) stream = ParallelStencil.ParallelKernel.get_priority_custream(1) $(args[end]))
    end
end

## FUNCTIONS TO GET INNER AND OUTER RANGES AND TO PROMOTE BOUNDARY_WIDTH TO 3D

promote_boundary_width(boundary_width::BOUNDARY_WIDTH_TYPE_1D)       = (boundary_width,    0, 0)
promote_boundary_width(boundary_width::BOUNDARY_WIDTH_TYPE_1D_TUPLE) = (boundary_width..., 0, 0)
promote_boundary_width(boundary_width::BOUNDARY_WIDTH_TYPE_2D)       = (boundary_width..., 0)
promote_boundary_width(boundary_width::BOUNDARY_WIDTH_TYPE)          = boundary_width
promote_boundary_width(boundary_width)                               = @ArgumentError("@hide_communication: boundary_width must be a Tuple of Integer of size 1, 2 or 3 (obtained: $boundary_width; its type is: $(typeof(boundary_width))).")

function get_ranges_outer(boundary_width, ranges::RANGES_TYPE) where T <:Integer
    boundary_width = promote_boundary_width(boundary_width)
    validate_ranges_args(boundary_width, ranges)
    ms = length.(ranges)
    bw = boundary_width
    if ms[3] > 1 # 3D
        ranges_outer = ((1:ms[1],             1:ms[2],             1:bw[3]),
                        (1:ms[1],             1:ms[2],             ms[3]-bw[3]+1:ms[3]),
                        (1:ms[1],             1:bw[2],             bw[3]+1:ms[3]-bw[3]),
                        (1:ms[1],             ms[2]-bw[2]+1:ms[2], bw[3]+1:ms[3]-bw[3]),
                        (1:bw[1],             bw[2]+1:ms[2]-bw[2], bw[3]+1:ms[3]-bw[3]),
                        (ms[1]-bw[1]+1:ms[1], bw[2]+1:ms[2]-bw[2], bw[3]+1:ms[3]-bw[3]),
                       )
    elseif ms[2] > 1 # 2D
        ranges_outer = ((1:ms[1],             1:bw[2],             1:1),
                        (1:ms[1],             ms[2]-bw[2]+1:ms[2], 1:1),
                        (1:bw[1],             bw[2]+1:ms[2]-bw[2], 1:1),
                        (ms[1]-bw[1]+1:ms[1], bw[2]+1:ms[2]-bw[2], 1:1),
                       )
    elseif ms[1] > 1  # 1D
        ranges_outer = ((1:bw[1],             1:1, 1:1),
                        (ms[1]-bw[1]+1:ms[1], 1:1, 1:1),
                       )
    else
        @ModuleInternalError("invalid argument 'ranges'.")
    end
    return Tuple([r for r in ranges_outer if all(length.(r) .!= 0)])
end

function get_ranges_inner(boundary_width, ranges::RANGES_TYPE) where T <:Integer
    boundary_width = promote_boundary_width(boundary_width)
    validate_ranges_args(boundary_width, ranges)
    ms = length.(ranges)
    bw = boundary_width
    if     (ms[3] > 1) return ( (bw[1]+1:ms[1]-bw[1], bw[2]+1:ms[2]-bw[2], bw[3]+1:ms[3]-bw[3]), ) # 3D
    elseif (ms[2] > 1) return ( (bw[1]+1:ms[1]-bw[1], bw[2]+1:ms[2]-bw[2], 1:1), ) # 2D
    elseif (ms[1] > 1) return ( (bw[1]+1:ms[1]-bw[1], 1:1, 1:1), ) # 1D
    else @ModuleInternalError("invalid argument 'ranges'.")
    end
end

function validate_ranges_args(boundary_width::Tuple{T,T,T}, ranges::RANGES_TYPE) where T <:Integer
    maxsize = length.(ranges)
    if any([r.start != 1 for r in ranges]) @ArgumentError("unsupported call to `@hide_communication boundary_width block`: automatic ranges determination for @hide_communication is not supported for ranges with start != 1. Use the explicit version `@hide_communication ranges_outer ranges_inner block` instead.") end
    if !all(maxsize .> 0) @IncoherentArgumentError("incoherent arguments in @hide_communication: the following must be true: all(length.(ranges) .> 0) ") end
    if !all(maxsize .> 2 .* boundary_width) @IncoherentArgumentError("incoherent arguments in @hide_communication: the following must be true: all(length.(ranges) .> 2 .* boundary_width) ") end # NOTE: this ensures among other things that if maxsize == 1, then boundary_width==0.
    if any((maxsize .> 1) .& (boundary_width .== 0)) @IncoherentArgumentError("incoherent arguments in @hide_communication: the following must be false: any((length.(ranges) .> 1) .& (boundary_width .== 0)) ") end
    if all(maxsize .== 1) @IncoherentArgumentError("incoherent arguments in @hide_communication: the following must be false: all(length.(ranges) .== 1)") end
end

## FUNCTIONS TO GET CREATE AND MANAGE CUDA STREAMS

@static if ENABLE_CUDA
    let
        global get_priority_custream, get_custream
        priority_custreams = Array{CuStream}(undef, 0)
        custreams          = Array{CuStream}(undef, 0)

        function get_priority_custream(id::Integer)
            while (id > length(priority_custreams)) push!(priority_custreams, CuStream(; flags=CUDA.STREAM_NON_BLOCKING, priority=CUDA.priority_range()[end])) end # CUDA.priority_range()[end] is max priority. # NOTE: priority_range cannot be called outside the function as only at runtime sure that CUDA is functional.
            return priority_custreams[id]
        end

        function get_custream(id::Integer)
            priority_min =
            while (id > length(custreams)) push!(custreams, CuStream(; flags=CUDA.STREAM_NON_BLOCKING, priority=CUDA.priority_range()[1])) end # CUDA.priority_range()[1] is min priority. # NOTE: priority_range cannot be called outside the function as only at runtime sure that CUDA is functional.
            return custreams[id]
        end
    end
end


text:
## This is based on the definition of a Beta(α,β) random variable
@inline function sampleBeta(α::Float64, β::Float64,
  rng::RNG = GLOBAL_RNG) where RNG <: AbstractRNG
  X::Float64 = sampleGamma(α, rng)
  Y::Float64 = sampleGamma(β, rng)
  return X/(X+Y)
end


text:
import Symata: unprotect, protect, @sjdoc

## These macro calls are also valid Symata. But, we do not include them in the
## autoloaded Symata code, because we want the text to be searchable before it
## is loadable. e.g. h"interval"

@sjdoc ExpToTrig """
    ExpToTrig(expr)

replace exponentials with trigonometric functions in `expr`.
"""

protect(:ExpToTrig)

@sjdoc Array """
    Array(f,n)

return a `List` of length `n` with elements `f(i)` for `i=1,...,n`

    Array(f,r,n)

return a `List` of length `n` with origin `r`.

    Array(f, n, [a,b])

return a `List` of length `n` of elements `f(v)` with `v` taking equally spaced values from `a` to `b`.

    Array(f, [n,m])

return a `List` of `List`s with elements `f(i,j)` for `i=1,...,n`, `j=1,...,m`
"""

protect(:Array)

@sjdoc Subdivide """
    Subdivide(n)

return a list that divides the unit interval into `n` intervals of equal length.

    Subdivide(xmax,n)

divide the interval `(0,xmax)` into `n` intervals

    Subdivide(xmin,xmax,n)

divide the interval `(xmin,xmax)` into `n` intervals.
"""

@sjdoc TakeDrop """
    TakeDrop(expr, seq)

returns `[Take(expr,seq), Drop(expr,seq)]`.
"""

protect(:TakeDrop)

@sjdoc UpTo """
    Take(expr, UpTo(n))

return `expr` keeping only the first `n` elements, or all elements
is there are less than `n` of them.
"""

protect(:UpTo)

@sjdoc ArrayDepth """
    ArrayDepth(m)

return the depth to which `m` is an array. If `m` is viewed as
a tensor, this is the tensor rank.

`ArrayDepth(m)` is equivalent to `Length(Dimensions(m))`.
"""

@sjdoc TensorRank """
    TensorRank(m)

return the depth to which `m` is an array. If `m` is viewed as
a tensor, this is the tensor rank.
"""

@sjdoc Divide """
    Divide(x,y)

is equivalent to `x * y^(-1)`.
"""


text:
"""
# module laminart

- Julia version: 1.4
- Author: niallcullinane
- Date: 2020-06-07

# Examples

```jldoctest
julia>
```
"""
module LaminartConv
include("./LamKernels.jl")

using NNlib, ImageFiltering, Images, OffsetArrays, CUDA
# , MEngProject.LamKernels

export I_u, fun_v_C, fun_equ

struct MyFunction{T} <: Function
  x_lgn::T
  C::T
  H_z::T
end

# function (ff::MyFunction)(du, u, p, t)
# # function f!(du, u, p, t)
#     @inbounds begin
#         x = CuArray(@view u[:, :, 1:p.K,:])
#         y = CuArray(@view u[:, :, p.K+1:2*p.K,:])
#         m = CuArray(@view u[:, :, 2*p.K+1:3*p.K,:])
#         z = CuArray(@view u[:, :, 3*p.K+1:4*p.K,:])
#         s = CuArray(@view u[:, :, 4*p.K+1:5*p.K,:])

#         v_p = CuArray(@view u[:, :, 5*p.K+1:5*p.K+1,:])
#         v_m = CuArray(@view u[:, :, 5*p.K+2:5*p.K+2,:])

#         dx = CuArray(@view du[:, :, 1:p.K,:])
#         dy = CuArray(@view du[:, :, p.K+1:2*p.K,:])
#         dm = CuArray(@view du[:, :, 2*p.K+1:3*p.K,:])
#         dz = CuArray(@view du[:, :, 3*p.K+1:4*p.K,:])
#         ds = CuArray(@view du[:, :, 4*p.K+1:5*p.K,:])

#         dv_p = CuArray(@view du[:, :, 5*p.K+1:5*p.K+1,:])
#         dv_m = CuArray(@view du[:, :, 5*p.K+2:5*p.K+2,:])

#         fun_x_lgn!(ff.x_lgn, x, p)
#         fun_v_C!(ff.C, v_p, v_m, p)
#         fun_H_z!(ff.H_z, z, p)

#         fun_dv!(dv_p, v_p, p.r, ff.x_lgn, p)
#         fun_dv!(dv_m, v_m, .-p.r, ff.x_lgn, p)
#         fun_dx_v1!(dx, x, ff.C, z, p.x_V2, p)
#         fun_dy!(dy, y, ff.C, x, m, p)
#         fun_dm!(dm, m, x, p)
#         fun_dz!(dz, z, y, ff.H_z, s, p)
#         fun_ds!(ds, s, ff.H_z, p)

#     end
#     return nothing
# end

function (ff::MyFunction)(du, u, p, t)
# function f!(du, u, p, t)
#     @inbounds begin
        x = @view u[:, :, 1:p.K,:]
        y = @view u[:, :, p.K+1:2*p.K,:]
        m = @view u[:, :, 2*p.K+1:3*p.K,:]
        z = @view u[:, :, 3*p.K+1:4*p.K,:]
        s = @view u[:, :, 4*p.K+1:5*p.K,:]

        v_p = @view u[:, :, 5*p.K+1:5*p.K+1,:]
        v_m = @view u[:, :, 5*p.K+2:5*p.K+2,:]

        dx = @view du[:, :, 1:p.K,:]
        dy = @view du[:, :, p.K+1:2*p.K,:]
        dm = @view du[:, :, 2*p.K+1:3*p.K,:]
        dz = @view du[:, :, 3*p.K+1:4*p.K,:]
        ds = @view du[:, :, 4*p.K+1:5*p.K,:]

        dv_p = @view du[:, :, 5*p.K+1:5*p.K+1,:]
        dv_m = @view du[:, :, 5*p.K+2:5*p.K+2,:]

        fun_x_lgn!(ff.x_lgn, x, p)
        fun_v_C!(ff.C, v_p, v_m, p)
        fun_H_z!(ff.H_z, z, p)

        fun_dv!(dv_p, v_p, p.r, ff.x_lgn, p)
        fun_dv!(dv_m, v_m, .-p.r, ff.x_lgn, p)
        fun_dx_v1!(dx, x, ff.C, z, p.x_V2, p)
#         fun_dy!(dy, y, ff.C, x, m, p)
#         fun_dm!(dm, m, x, p)
#         fun_dz!(dz, z, y, ff.H_z, s, p)
#         fun_ds!(ds, s, ff.H_z, p)

#     end
    return nothing
end

# function (ff::MyFunction)(du, u, p, t)
# # function f!(du, u, p, t)
#     @inbounds begin
#         x = @view u[:, :, 1:p.K,:]
#         y = @view u[:, :, p.K+1:2*p.K,:]
#         m = @view u[:, :, 2*p.K+1:3*p.K,:]
#         z = @view u[:, :, 3*p.K+1:4*p.K,:]
#         s = @view u[:, :, 4*p.K+1:5*p.K,:]

#         v_p = @view u[:, :, 5*p.K+1:5*p.K+1,:]
#         v_m = @view u[:, :, 5*p.K+2:5*p.K+2,:]

#         dx = @view du[:, :, 1:p.K,:]
#         dy = @view du[:, :, p.K+1:2*p.K,:]
#         dm = @view du[:, :, 2*p.K+1:3*p.K,:]
#         dz = @view du[:, :, 3*p.K+1:4*p.K,:]
#         ds = @view du[:, :, 4*p.K+1:5*p.K,:]

#         dv_p = @view du[:, :, 5*p.K+1:5*p.K+1,:]
#         dv_m = @view du[:, :, 5*p.K+2:5*p.K+2,:]

#         fun_x_lgn!(ff.x_lgn[:,:,:,:], x[:,:,:,:], p)
#         fun_v_C!(ff.C[:,:,:,:], v_p[:,:,:,:], v_m[:,:,:,:], p)
#         fun_H_z!(ff.H_z[:,:,:,:], z[:,:,:,:], p)

#         fun_dv!(dv_p[:,:,:,:], v_p[:,:,:,:], p.r[:,:,:,:], ff.x_lgn[:,:,:,:], p)
#         fun_dv!(dv_m[:,:,:,:], v_m[:,:,:,:], .-p.r[:,:,:,:], ff.x_lgn[:,:,:,:], p)
#         fun_dx_v1!(dx[:,:,:,:], x[:,:,:,:], ff.C[:,:,:,:], z[:,:,:,:], p.x_V2[:,:,:,:], p)
#         fun_dy!(dy[:,:,:,:], y[:,:,:,:], ff.C[:,:,:,:], x[:,:,:,:], m[:,:,:,:], p)
#         fun_dm!(dm[:,:,:,:], m[:,:,:,:], x[:,:,:,:], p)
#         fun_dz!(dz[:,:,:,:], z[:,:,:,:], y[:,:,:,:], ff.H_z[:,:,:,:], s[:,:,:,:], p)
#         fun_ds!(ds[:,:,:,:], s[:,:,:,:], ff.H_z[:,:,:,:], p)

#     end
#     return nothing
# end

function kernels(img::AbstractArray, p::NamedTuple)
       C_Q_temp = reshape(
        Array{eltype(img)}(undef, p.C_AB_l, p.C_AB_l * p.K),
        p.C_AB_l,
        p.C_AB_l,
        1,
		p.K
    )
C_P_temp = similar(C_Q_temp)
	    H_temp = reshape(
        zeros(eltype(img), p.H_l, p.H_l * p.K * p.K),
        p.H_l,
        p.H_l,
        p.K,
    p.K)
 T_temp = reshape(Array{eltype(img)}(undef, p.K * p.K), 1, 1, p.K, p.K)
 W_temp =
        reshape(Array{eltype(img)}(undef, p.W_l, p.W_l * p.K * p.K), p.W_l, p.W_l, p.K, p.K)
    for k ∈ 1:p.K
        θ = π * (k - 1.0f0) / p.K
        C_Q_temp[:, :, 1,k] = LamKernels.kern_d(p.σ_2, θ)
        C_P_temp[:, :, 1,k] = LamKernels.kern_b(p.σ_2, θ)
        H_temp[:, :, k,k] = p.H_fact .* LamKernels.gaussian_rot(p.H_σ_x, p.H_σ_y, θ, p.H_l)
# 		todo make T kernel more general for higher K
        T_temp[1, 1, k,1] = p.T_fact[k]
        T_temp[1, 1, 2,2] = p.T_fact[1]
        T_temp[1, 1, 1,2] = p.T_fact[2]
        #todo: generalise T and W for higher K
        #         T_temp[:,:,k] = KernelFactors.gaussian(p.T_σ, p.K)
        #         for l ∈ 1:p.K
        #             W_temp[:,:,l,k] =
        #         end
    end

    W_temp[:, :, 1, 1] =
        5f0 .* LamKernels.gaussian_rot(3f0, 0.8f0, 0f0, p.W_l) .+
        LamKernels.gaussian_rot(0.4f0, 1f0, 0f0, p.W_l)

	W_temp[:, :, 2, 2] =
        5f0 .* LamKernels.gaussian_rot(3f0, 0.8f0,  π / 2f0, p.W_l) .+
        LamKernels.gaussian_rot(0.4f0, 1f0, π / 2f0, p.W_l)

	W_temp[:, :, 1, 2] = relu.(
        0.2f0 .- LamKernels.gaussian_rot(2f0, 0.6f0, 0f0, p.W_l) .-
        LamKernels.gaussian_rot(0.3f0, 1.2f0, 0f0, p.W_l))

	W_temp[:, :, 2, 1] = relu.(
        0.2f0 .- LamKernels.gaussian_rot(2f0, 0.6f0, π / 2f0, p.W_l) .-
        LamKernels.gaussian_rot(0.3f0, 1.2f0, π / 2f0, p.W_l))

#     W_temp[:, :, 1, 1] =
#         5f0 .* LamKernels.gaussian_rot(3f0, 0.8f0, 0f0, p.W_l) .+
#         LamKernels.gaussian_rot(0.4f0, 1f0, 0f0, p.W_l)
#     W_temp[:, :, 2, 2] =
#         5f0 .* LamKernels.gaussian_rot(3f0, 0.8f0, 0f0, p.W_l) .+
#         LamKernels.gaussian_rot(0.4f0, 1f0, π / 2f0, p.W_l)
#     W_temp[:, :, 1, 2] = relu.(
#         0.2f0 .- LamKernels.gaussian_rot(2f0, 0.6f0, 0f0, p.W_l) .-
#         LamKernels.gaussian_rot(0.3f0, 1.2f0, 0f0, p.W_l))
#     W_temp[:, :, 2, 1] = relu.(
#         0.2f0 .- LamKernels.gaussian_rot(2f0, 0.6f0, 0f0, p.W_l) .-
#         LamKernels.gaussian_rot(0.3f0, 1.2f0, π / 2f0, p.W_l))

temp_out = (
        k_gauss_1 = reshape2d_4d(Kernel.gaussian(p.σ_1)),
        k_gauss_2 = reshape2d_4d(Kernel.gaussian(p.σ_2)),
        k_C_d = C_Q_temp,
        k_C_b = C_P_temp,

# 		todo use mean of x_lgn?
		k_x_lgn = reshape(ones(Float32,1,p.K),1,1,p.K,1),
# 		k_x_lgn = reshape(ones(Float32,1,p.K),1,1,p.K,1))./p.K,
        k_W_p = W_temp,
        k_W_m = W_temp,
        k_H = H_temp,
        k_T_p = T_temp,
        k_T_m = p.T_p_m .* T_temp,
        k_T_p_v2 = p.T_v2_fact .* T_temp,
        k_T_m_v2 = p.T_v2_fact .* p.T_p_m .* T_temp,
        dim_i = size(img)[1],
        dim_j = size(img)[2],
        x_V2 = reshape(zeros(Float32, size(img)[1], size(img)[2] * p.K), size(img)[1], size(img)[2],p.K,1),
ν_pw_n= p.ν^p.n, )
merge(p, temp_out)
end

function reshape2d_4d(img::AbstractArray)
    reshape(img, size(img)[1], size(img)[2], 1, 1)
end

# function reshape_ijk_ij1k(img::AbstractArray, p::NamedTuple)
#     reshape(img, size(img)[1], size(img)[2], 1, p.K)
# end

function conv!(out::AbstractArray, img::AbstractArray, kern::AbstractArray, p::NamedTuple)
	out_ = @view out[:,:,:,:]
    out_ .= NNlib.conv(img, kern, pad=(size(kern)[1]>>1, size(kern)[1]>>1, size(kern)[2]>>1, size(kern)[2]>>1), flipped=true)
# 	@. out_ = out
    return nothing
end

function conv_dw!(out::AbstractArray, img::AbstractArray, kern::AbstractArray, p::NamedTuple)
     out = NNlib.depthwiseconv(img, kern, pad=(size(kern)[1]>>1, size(kern)[2]>>1, size(kern)[3]>>1, size(kern)[4]>>1), flipped=true)
    return nothing
end

# function varables(I::AbstractArray, p::NamedTuple)
#    v_p = zeros(typeof(I[1,1]), size(I)[1], size(I)[2])
#    v_m = copy(v_p)
#    x_lgn = copy(v_p)
#    x = reshape(zeros(typeof(I[1,1]), size(I)[1], size(I)[2] * p.K), size(I)[1], size(I)[2], p.K)
#    y = copy(x)
#    m = copy(x)
#    z = copy(x)
#    s = copy(x)
#    C = copy(x)
#    H_z = copy(x)
#    x_V2 = copy(x)
#    y_V2 = copy(x)
#    m_V2 = copy(x)
#    z_V2 = copy(x)
#    s_V2 = copy(x)
#    H_z_V2 = copy(x)
#    return [v_p, v_m, x_lgn, x, y, m, z, s, C, H_z, x_V2, y_V2, m_V2, z_V2, s_V2, H_z_V2]
# end
#
#
# function variables_sep(I::AbstractArray, p::NamedTuple)
#    v_p = zeros(typeof(I[1,1]), size(I)[1], size(I)[2])
#    v_m = copy(v_p)
#    x_lgn = copy(v_p)
#    x = reshape(zeros(typeof(I[1,1]), size(I)[1], size(I)[2] * p.K), size(I)[1], size(I)[2], p.K)
#    y = copy(x)
#    m = copy(x)
#    z = copy(x)
#    s = copy(x)
#    C = copy(x)
#    H_z = copy(x)
#    x_V2 = copy(x)
#    y_V2 = copy(x)
#    m_V2 = copy(x)
#    z_V2 = copy(x)
#    s_V2 = copy(x)
#    H_z_V2 = copy(x)
#    return v_p, v_m, x_lgn, x, y, m, z, s, C, H_z, x_V2, y_V2, m_V2, z_V2, s_V2, H_z_V2
# end
#
# function variables_dict(I::AbstractArray, p::NamedTuple)
#    v_p = zeros(typeof(I[1,1]), size(I)[1], size(I)[2])
#    v_m = copy(v_p)
#    x_lgn = copy(v_p)
#    x = reshape(zeros(typeof(I[1,1]), size(I)[1], size(I)[2] * p.K), size(I)[1], size(I)[2], p.K)
#    y = copy(x)
#    m = copy(x)
#    z = copy(x)
#    s = copy(x)
#    C = copy(x)
#    H_z = copy(x)
#    x_V2 = copy(x)
#    y_V2 = copy(x)
#    m_V2 = copy(x)
#    z_V2 = copy(x)
#    s_V2 = copy(x)
#    H_z_V2 = copy(x)
#    return v_p, v_m, x_lgn, x, y, m, z, s, C, H_z, x_V2, y_V2, m_V2, z_V2, s_V2, H_z_V2
# end

# function varables(I::AbstractArray, p::NamedTuple)
#    ij = zeros(typeof(I[1,1]), size(I)[1], size(I)[2])
#    ijk = reshape(zeros(typeof(I[1,1]), size(I)[1], size(I)[2] * p.K), size(I)[1], size(I)[2], p.K)
#    temp_out = [v_p = copy(ij),
#    v_m = copy(ij),
#    x_lgn = copy(ij),
#    x = copy(ijk),
#    y = copy(ijk),
#    m = copy(ijk),
#    z = copy(ijk),
#    s = copy(ijk),
#    C = copy(ijk),
#    H_z = copy(ijk),
#    x_V2 = copy(ijk),
#    y_V2 = copy(ijk),
#    m_V2 = copy(ijk),
#    z_V2 = copy(ijk),
#    s_V2 = copy(ijk),
#    H_z_V2 = copy(ijk)]
#    return temp_out
# end

function add_I_u_p(I::AbstractArray, p::NamedTuple)
# 	todo fix
	I_4d = CuArray(reshape2d_4d(I))
	r = similar(I_4d)
	I_u!(r, I_4d, p)
    temp_out = (I = I_4d, r = r)
    return merge(p, temp_out)
end

# retina

function I_u!(r::AbstractArray, I::AbstractArray, p::NamedTuple)

	conv!(r, I, p.k_gauss_1, p)
    @. r = I - r
end

# lgn feedback

function fun_x_lgn!(x_lgn::AbstractArray, x::AbstractArray, p::NamedTuple)
	out_ = @view x_lgn[:,:,:,:]
	    out_ .= NNlib.conv(x[:,:,:,:], p.k_x_lgn, pad=0,flipped=true)
    return nothing
end

# function fun_F(value::Real, p::NamedTuple)
#     max.(value - p.Γ, 0)
# end
#
#
# # williomson uses differnt F, relu with threshold
# function fun_F_willimson(value::Real, p::NamedTuple)
#     value < p.Γ ? zero(value) : value
# end

function fun_f!(f_out::AbstractArray, x::AbstractArray, p::NamedTuple)
    @. f_out = (p.μ * x^p.n) / (p.ν_pw_n + x^p.n)
    return nothing
end

function func_filter_W!(
    W_out::AbstractArray,
    img::AbstractArray,
    W::AbstractArray,
    p::NamedTuple,
)
    temp_k = similar(W_out[:, :, 1])
    # @inbounds begin
        for k ∈ 1:p.K
            #     todo fix W
            img_k = @view img[:, :, k]
            out_k = @view W_out[:, :, k]
            imfilter!(out_k, img_k, centered(W[:, :, k, k]), p.filling)
            for l ∈ 1:p.K
                if l ≠ k
                    img_l = @view img[:, :, l]
                    imfilter!(
                        temp_k,
                        img_l,
                        centered(W[:, :, k, l]),
                        p.filling,
                    )
                    @. out_k += temp_k
                end
            end
        end
    # end
    return nothing
end

# LGN
function fun_dv!(
    dv::AbstractArray,
    v::AbstractArray,
    u::AbstractArray,
    x_lgn::AbstractArray,
    p::NamedTuple,
)
	conv!(dv, x_lgn, p.k_gauss_1, p)
    @. dv =
        p.δ_v * (
            -v + ((1f0 - v) * max(u, 0f0) * (1f0 + p.C_1 * x_lgn)) -
            ((1f0 + v) * p.C_2 * dv)
        )
    return nothing
end

# lgn to l6 and l4
function fun_v_C!(
    v_C::AbstractArray,
    v_p::AbstractArray,
    v_m::AbstractArray,
    p::NamedTuple,
)
    V = similar(v_p)
    temp = similar(v_p)

    @. temp = exp(-1.0f0 / 8.0f0) * (max(v_p, 0) - max(v_m, 0))
    conv!(V, temp, p.k_gauss_2, p)

    A = similar(v_C)
    #     allocate B to v_C

	conv!(A, V, p.k_C_d, p)
	conv!(v_C, V, p.k_C_b, p)

    @. v_C = p.γ * (max(A - abs(v_C), 0) + max(-A - abs(v_C), 0))
    return nothing
end

# L6
function fun_dx_v1!(
    dx::AbstractArray,
    x::AbstractArray,
    C::AbstractArray,
    z::AbstractArray,
    x_v2::AbstractArray,
    p::NamedTuple,
)
    @. dx =
        p.δ_c * (
            -x + (
                (1.0f0 - x) * (
                    (p.α * C) + (p.ϕ * max(z - p.Γ, 0)) .+ (p.v_21 * x_v2) +
                    p.att
                )
            )
        )
    return nothing
end

#     L4 excit
function fun_dy!(
    dy::AbstractArray,
    y::AbstractArray,
    C::AbstractArray,
    x::AbstractArray,
    m::AbstractArray,
    p::NamedTuple,
)
    conv!(dy, m, p.k_W_p, p)
    @. dy = m * dy
    fun_f!(dy, dy, p)
    @. dy = p.δ_c * (-y + ((1 - y) * (C + (p.η_p * x))) - ((1 + y) * dy))
    return nothing
end

#     l4 inhib
function fun_dm!(
    dm::AbstractArray,
    m::AbstractArray,
    x::AbstractArray,
    p::NamedTuple,
)
    conv!(dm, m, p.k_W_m, p)
    fun_f!(dm, dm, p)
    @. dm = p.δ_m * (-m + (p.η_m * x) - (m * dm))
    return nothing
end

#     L2/3 excit
function fun_dz!(
    dz::AbstractArray,
    z::AbstractArray,
    y::AbstractArray,
    H_z::AbstractArray,
    s::AbstractArray,
    p::NamedTuple,
)
    conv!(dz, s, p.k_T_p, p)
    @. dz =
        p.δ_z * (
            -z + ((1 - z) * ((p.λ * max(y, 0)) + H_z + (p.a_23_ex * p.att))) -
            ((z + p.ψ) * dz)
        )
    return nothing
end

#     L2/3 inhib
function fun_ds!(
    ds::AbstractArray,
    s::AbstractArray,
    H_z::AbstractArray,
    p::NamedTuple,
)
    conv!(ds, s, p.k_T_m, p)
    @. ds = p.δ_s * (-s + H_z + (p.a_23_in * p.att) - (s * ds))
    return nothing
end

function fun_H_z!(H_z::AbstractArray, z::AbstractArray, p::NamedTuple)
	conv!(H_z, max.(z .- p.Γ, 0f0), p.k_H, p)
    return nothing
end

#     V2 L6
function fun_dx_v2!(
    dx::AbstractArray,
    x_v2::AbstractArray,
    z_v2::AbstractArray,
    z::AbstractArray,
    p::NamedTuple,
)
    @. dx =
        p.δ_c * (
            -x_v2 + (
                (1 - x_v2) * (
                    (p.v12_6 * max(z - p.Γ, 0)) +
                    (p.ϕ * max(z_v2 - p.Γ, 0)) +
                    p.att
                )
            )
        )
    return nothing
end

# V2 L4 excit
function fun_dy_v2!(
    dy::AbstractArray,
    y_v2::AbstractArray,
    z::AbstractArray,
    x_v2::AbstractArray,
    m_v2::AbstractArray,
    p::NamedTuple,
)
    conv!(dy, m_v2, p.k_W_p, p)
    fun_f!(dy, dy, p)
    @. dy =
        p.δ_c * (
            -y_v2 +
            ((1.0f0 - y_v2) * ((p.v12_4 * max(z - p.Γ, 0)) + (p.η_p * x_v2))) -
            ((1.0f0 + y_v2) * dy)
        )
    return nothing
end

# # for equilabrium
# fun_equ(x) = x/(1+x)
#
#
# # lgn
# function fun_v_equ(u::AbstractArray, x::AbstractArray, p::NamedTuple)
#     return fun_equ.(relu.(u) .* (p.lgn_equ_u .+ (p.lgn_equ_A .* fun_A.(x))) .- (p.lgn_equ_B .* fun_B.(x)))
# end
#
#
# # lgn, no L6 feedback, light
# function fun_v_equ_noFb(u::AbstractArray, x::AbstractArray, p::NamedTuple)
#     return fun_equ.(relu.(u) .* (p.lgn_para_u))
# end
#
#
# # l6
# function fun_x_equ(C::AbstractArray, z::AbstractArray,  p::NamedTuple)
#     return fun_equ.((p.α .* C) .+ (p.ϕ .* fun_F.(z,p.Γ)) + (p.V_21 .* x_v2) .+ p.att)
# end
#
#
# # l6, no V2 feedback, light
# function fun_x_equ_noV2(C::AbstractArray, z::AbstractArray, p::NamedTuple)
#     return fun_equ.((p.α .* C) .+ (p.ϕ .* fun_F.(z,p.Γ)))
# end
#
#
# # l4 excit
# function fun_y_equ(C::AbstractArray, x::AbstractArray,  m::AbstractArray, p::NamedTuple)
#     return fun_equ.(C .+ (η_p .* x - fun_f.(imfilter(m, p.k_W_p), p.μ, p.ν, p.n)))
# end
#
# # l4 inhib, needs initial condition of itself
# function fun_m_equ(C::AbstractArray, x::AbstractArray, m_init::AbstractArray, p::NamedTuple)
#     return (p.η_m .* x ./ (1 .+  fun_f.(imfilter(m_init, p.k_W_m), p.μ, p.ν, p.n)))
# end
#
#
# # l4 inhib - no feedback kernel
# function fun_m_equ_noFb(C::AbstractArray, x::AbstractArray, p::NamedTuple)
#     return (p.η_m .* x)
# end
#
#
# #  l2/3 excit, needs initial condition of itself
# function fun_z_equ(y::AbstractArray, s::AbstractArray,  z_init::AbstractArray,  p::NamedTuple)
#     return (p.λ .* relu.(y)) .+ imfilter(fun_F.(z_init, p.Γ), k_H) .+ (p.a_23_ex .* p.att) .- (p.ϕ .* imfilter(s, k_T_p)) ./
#     (1 .+ (p.λ .* relu.(y)) .+ imfilter(fun_F.(z_init, p.Γ), k_H) .+ (p.a_23_ex .* p.att) .+ imfilter(s, k_T_p))
#     return
# end
#
#
# #  l2/3 excit - no feedback kernel
# function fun_z_equ_noFb(y::AbstractArray, s::AbstractArray, p::NamedTuple)
#     return (p.λ .* relu.(y)) .- (p.ϕ .* imfilter(s, k_T_p)) ./
#     (1 .+ (p.λ .* relu.(y)) .+ imfilter(s, k_T_p))
#     return
# end
#
#
# #  l2/3 inhib, needs initial condition of itself
# function fun_s_equ(z::AbstractArray, s::AbstractArray,  s_init::AbstractArray, p::NamedTuple)
#     return ((imfilter(fun_F.(z, p.Γ), k_H) + p.a_23_in .* p.att ) ./ (1 .+ imfilter(s_init, k_T_m)))           #????????? is T right?
# end
#
#
# #  l2/3 inhib - no feedback kernel
# function fun_s_equ_noFb(z::AbstractArray, s::AbstractArray, p::NamedTuple)
#     return (imfilter(fun_F.(z, p.Γ), k_H))
# end
#
#
# # V2 L6
# function fun_xV2_equ(x_v2::AbstractArray, z::AbstractArray, z_v2::AbstractArray, p::NamedTuple)
#     return fun_equ.((p.V_12_6 .* fun_F.(z, p.Γ) .+ (p.ϕ .* fun_F.(z_v2, p.Γ))))
# end
#
#
# #  V2 L4 excit
# function fun_yV2_equ(z::AbstractArray, x::AbstractArray,  m::AbstractArray, p::NamedTuple)
#     return (p.V_12_4 .* fun_F.(z,p.Γ) .+ (p.η_p .* x) - fun_f.(imfilter(m, k_W_p), p.μ, p.ν, p.n) ./
#    (1 .+ p.V_12_4 .* fun_F.(z,p.Γ) .+ (p.η_p .* x) .+ fun_f.(imfilter(m, k_W_p), p.μ, p.ν, p.n)))
# end

end


text:
# MLP Dynamics conditioned on the Time Step
struct TDChain{T<:Tuple}
    layers::T
    TDChain(xs...) = new{typeof(xs)}(xs)
end

@forward TDChain.layers Base.getindex,
Base.length,
Base.first,
Base.last,
Base.iterate,
Base.lastindex

Flux.functor(::Type{<:TDChain}, c) = c.layers, ls -> TDChain(ls...)

applytdchain(::Tuple{}, x, t) = x
applytdchain(fs::Tuple, x, t) = applytdchain(Base.tail(fs), first(fs)(vcat(x, t)), t)

function (c::TDChain)(x::AbstractMatrix, t)
    _t = similar(x, 1, size(x, 2))
    fill!(_t, t)
    return applytdchain(c.layers, x, _t)
end

function (c::TDChain)(x::AbstractMatrix, t::TrackedReal)
    _t = CUDA.ones(Float32, 1, size(x, 2)) .* t
    return applytdchain(c.layers, x, _t)
end

Base.getindex(c::TDChain, i::AbstractArray) = TDChain(c.layers[i]...)

Flux.testmode!(m::TDChain, mode = true) = (map(x -> testmode!(x, mode), m.layers); m)

function Base.show(io::IO, c::TDChain)
    print(io, "TimeDependentChain(")
    join(io, c.layers, ", ")
    print(io, ")")
end

# Some Common Network Layers

## Recognition RNN for mapping from time series data to a latent vector
struct RecognitionRNN{P,Q,R}
    i2h::Dense{P,Q,R}
    h2o::Dense{P,Q,R}
end

Flux.@functor RecognitionRNN

RecognitionRNN(; latent_dim::Int, nhidden::Int, obs_dim::Int) =
    RecognitionRNN(Dense(obs_dim + nhidden, nhidden), Dense(nhidden, latent_dim * 2))

function (rrnn::RecognitionRNN)(x, h)
    h = tanh.(rrnn.i2h(vcat(x, h)))
    return rrnn.h2o(h), h
end


text:
"""
	sphinterpolate(cmd0::String="", arg1=nothing, kwargs...)

Spherical gridding in tension of data on a sphere

Full option list at [`sphinterpolate`]($(GMTdoc)sphinterpolate .html)

Parameters
----------

- **D** | **duplicates** :: [Type => Bool]

    Delete any duplicate points [Default assumes there are no duplicates].
    ($(GMTdoc)sphinterpolate.html#d)
- **G** | **save** | **outgrid** | **outfile** :: [Type => Str]

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = sphinterpolate(....) form.
    ($(GMTdoc)sphinterpolate.html#g)
- $(GMT.opt_I)
    ($(GMTdoc)sphinterpolate.html#i)
- **Q** | **tension** :: [Type => Number | Str]     ``Arg = mode[/options]``

    Specify one of four ways to calculate tension factors to preserve local shape properties or satisfy arc constraints.
    ($(GMTdoc)sphinterpolate.html#q)
- **T** | **var_tension** :: [Type => Bool | Str]

    Use variable tension (ignored with -Q0 [constant]
    ($(GMTdoc)sphinterpolate.html#t)
- **Z** | **scale** :: [Type => Bool | Str]

    Before interpolation, scale data by the maximum data range [no scaling].
    ($(GMTdoc)sphinterpolate.html#z)
- $(GMT.opt_R)
- $(GMT.opt_V)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_r)
- $(GMT.opt_swap_xy)
"""
function sphinterpolate(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && return monolitic("sphinterpolate ", cmd0, arg1)

	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode
	cmd, = parse_common_opts(d, "", [:G :I :R :V_params :bi :di :e :h :i :r :yx])
	cmd  = parse_these_opts(cmd, d, [[:D :duplicates], [:Q :tension], [:T :nodetable], [:Z :scale]])

	common_grd(d, cmd0, cmd, "sphinterpolate ", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
sphinterpolate(arg1, cmd0::String=""; kw...) = sphinterpolate(cmd0, arg1; kw...)


text:
@testset "MetisPartitioner.jl" begin
    mesh = buildSquare(0, 1, 0, 1, 2, 2)

    MOOSE.partition!(mesh, MOOSE.MetisPartitioner)
end


text:
using TikzPictures
using Test

svgBackends = [
    "testPic.pdf2svg.svg" => PdfToSvgBackend(),
    "testPic.poppler.svg" => PopplerBackend(),
    "testPic.dvisvgm.svg" => DVIBackend()
]
if VERSION < v"1.3"
    deleteat!(svgBackends, 2)
    @warn "Not testing PopplerBackend on Julia versions < 1.3"
end

# Pre-test cleanup (for repeated tests)
for file in ["testPic.pdf", "testPic.svg", "testDoc.pdf", "testDoc.tex", first.(svgBackends)...]
  if isfile(file)
    rm(file)
  end
end

# Run tests
data = "\\draw (0,0) -- (10,10);\n\\draw (10,0) -- (0,10);\n\\node at (5,5) {tikz \$\\sqrt{\\pi}\$};"
tp = TikzPicture(data, options="scale=0.25", preamble="")
td = TikzDocument()
push!(td, tp, caption="hello")

save(TEX("testPic"), tp)
@test isfile("testPic.tex")

# check that the TEX file contains the desired environments
function has_environment(content::String, environment::String)
    has_begin = occursin("\\begin{$environment}", content)
    has_end = occursin("\\end{$environment}", content)
    if has_begin && has_end
        return true # has both
    elseif !has_begin && !has_end
        return false # has neither
    else
        error("\\begin{$environment} and \\end{$environment} do not match")
    end
end

filecontent = join(readlines("testPic.tex", keep=true)) # read with line breaks
@test occursin(data, filecontent) # also check that the data is contained
@test has_environment(filecontent, "tikzpicture")
@test has_environment(filecontent, "document")

# same check for include_preamble=false and limit_to=:picture
save(TEX("testPic"; include_preamble=false), tp)
filecontent = join(readlines("testPic.tex", keep=true))
@test occursin(data, filecontent)
@test has_environment(filecontent, "tikzpicture") # must occur
@test !has_environment(filecontent, "document") # must not occur

# same check for limit_to=:data
save(TEX("testPic"; limit_to=:data), tp)
filecontent = join(readlines("testPic.tex", keep=true))
@test occursin(data, filecontent)
@test !has_environment(filecontent, "tikzpicture")
@test !has_environment(filecontent, "document")

save(TEX("testPic"), tp) # save again with limit_to=:all
if success(`lualatex -v`)
  save(PDF("testPic"), tp)
  @test isfile("testPic.pdf")

    save(SVG("testPic"), tp)
    @test isfile("testPic.svg") # default SVG backend

    @testset for (k, v) in svgBackends
        svgBackend(v)
        for usetectonic in [true, false]
            if usetectonic && v == DVIBackend()
                continue
            end
            tikzUseTectonic(usetectonic)
            save(SVG(k), tp)
            @test isfile(k)
            rm(k)
        end
    end

    save(PDF("testDoc"), td)
    @test isfile("testDoc.pdf")
else
    @warn "lualatex is missing; can not test compilation"
end

# Test tikz-cd

data = "A\\arrow{rd}\\arrow{r} & B \\\\& C"
tp = TikzPicture(data, options="scale=0.25", environment="tikzcd", preamble="\\usepackage{tikz-cd}")
td = TikzDocument()
push!(td, tp, caption="hello")

save(TEX("testCD"), tp)
@test isfile("testCD.tex")

filecontent = join(readlines("testCD.tex", keep=true)) # read with line breaks
@test occursin(data, filecontent) # also check that the data is contained
@test has_environment(filecontent, "tikzcd")
@test has_environment(filecontent, "document")

if success(`lualatex -v`)
    save(PDF("testCD"), tp)
    @test isfile("testCD.pdf")

    save(PDF("testCDDoc"), td)
    @test isfile("testCDDoc.pdf")
else
    @warn "lualatex is missing; can not test compilation"
end

# Test adjustbox width/height
adjustbox_width = TikzPicture(L"\node (sigmoid) [circle, draw=black, label=left:{$\sigma(z)$}] {};", width="2cm")
save(TEX("adjustbox_width"), adjustbox_width)
@test isfile("adjustbox_width.tex")

adjustbox_height = TikzPicture(L"\node (sigmoid) [circle, draw=black, label=left:{$\sigma(z)$}] {};", height="10cm")
save(TEX("adjustbox_height"), adjustbox_height)
@test isfile("adjustbox_height.tex")

adjustbox_aspect = TikzPicture(L"\node (sigmoid) [circle, draw=black, label=left:{$\sigma(z)$}] {};", height="10cm", keepAspectRatio=false)
save(TEX("adjustbox_aspect"), adjustbox_aspect)
@test isfile("adjustbox_aspect.tex")


text:
using ProgressBars, Printf, LinearAlgebra, Random

include("src/ISSNeuralODE.jl")
include("src/netlist.jl")

# set random seed
Random.seed!(0);

""" Generate synthetic training data """
datapath = joinpath("data", "demo")
if ~ispath(datapath)
    mkpath(datapath)
end

n, m = 3, 2

function sample_params(n, m)
    function _sample()
        W = randn(n, n)
        A = randn(n, n)
        B = randn(n, m)
        return W, A, B
    end
    W, A, B = _sample()
    while maximum(eigvals(W * A - I(n) + A' * W' - I(n))) >= 0
        W, A, B = _sample()
    end
    return W, A, B
end

W, A, B = sample_params(n, m)

# specify dynamics
function f(x, u)
    return -x + W * tanh.(A * x + B * u)
end

# data parameters
t_max = 50.0
num_knots = 20
num_samples = 500
δₜ = 0.1

# generate training data
# size (data dimension, time horizon, num_samples)
times = collect(0.0:δₜ:t_max)
udata = zeros(m, length(times), num_samples)
ydata = zeros(n, length(times), num_samples)

@printf "Generating training data..."
for i in ProgressBar(1:num_samples)
    # generate input signal
    tspan = (0.0, t_max)
    knots = vcat(0.0, sort(t_max * rand(num_knots - 2)), t_max)

    peaks = zeros(2, num_knots)
    for j = 1:num_knots
        peaks[:, j] = 2.0 * (rand(2) .- 0.5)
    end

    # linear interpolation for input
    u = (t) -> linear_interpolate(t, knots, peaks)

    # set initial condition to be equilibrium point
    ssprob = SteadyStateProblem((x, p, t) -> f(x, u(0.0)), zeros(3))
    x₀ = Array(solve(ssprob, SSRootfind()))

    # solve controlled ODE
    prob = ODEProblem((x, p, t) -> f(x, u(t)), x₀, tspan)
    y = Array(solve(prob, saveat = times, solver = Tsit5(), reltol = 1e-6, absol = 1e-9))

    for (j, s) in enumerate(times)
        udata[:, j, i] = u(s)
    end

    ydata[:, :, i] = y
end

""" preprocess data """
# dilate time
t_scale = 10.0
times *= t_scale

# scale inputs and outputs
umin = minimum(udata, dims = [2, 3])
umax = maximum(udata, dims = [2, 3])
ymin = minimum(ydata, dims = [2, 3])
ymax = maximum(ydata, dims = [2, 3])

uscale = 0.5 * (umax - umin)
yscale = 0.5 * (ymax - ymin)

ubias = -umin ./ uscale .- 1
ybias = -ymin ./ yscale .- 1

udata = udata ./ uscale .+ ubias
ydata = ydata ./ yscale .+ ybias

dataset = Dict([("t", times), ("u", udata), ("y", ydata)])

metadata = Dict([
    ("dt", δₜ),
    ("u_scale", uscale),
    ("u_bias", ubias),
    ("y_scale", yscale),
    ("y_bias", ybias),
    ("t_scale_factor", 1.0 / t_scale),
])

serialize("data/demo/demo_data.jls", dataset)
serialize("data/demo/demo_metadata.jls", metadata)

""" Learn model """#

include("configs/config_demo.jl")
run_config(demo_config)

max_stars_repo_name: syanga/ISSNeuralODE

max_issues_repo_name: syanga/ISSNeuralODE

max_forks_repo_name: syanga/ISSNeuralODE


text:
struct Bernoulli <: Distribution{Bool} end

"""
    bernoulli(prob_true::Real)

Samples a `Bool` value which is true with given probability
"""
const bernoulli = Bernoulli()

function logpdf(::Bernoulli, x::Bool, prob::Real)
    x ? log(prob) : log(1. - prob)
end

function logpdf_grad(::Bernoulli, x::Bool, prob::Real)
    prob_grad = x ? 1. / prob : -1. / (1-prob)
    (nothing, prob_grad)
end

random(::Bernoulli, prob::Real) = rand() < prob

(::Bernoulli)(prob) = random(Bernoulli(), prob)

has_output_grad(::Bernoulli) = false
has_argument_grads(::Bernoulli) = (true,)

export bernoulli

max_stars_repo_name: dendisuhubdy/Gen

max_issues_repo_name: kumsh/Gen

max_forks_repo_name: kumsh/Gen


text:
const _πA(mod) = _πACGT(mod)
const _πC(mod) = _πACGT(mod)
const _πG(mod) = _πACGT(mod)
const _πT(mod) = _πACGT(mod)
const _πACGU(mod) = _πACGT(mod)

_πACGT(mod::NASM) = 0.25

_π(mod::NASM) = SVector(_πA(mod), _πC(mod), _πG(mod), _πT(mod))

_πR(mod::NASM) = _πA(mod) + _πG(mod)

_πY(mod::NASM) = _πT(mod) + _πC(mod)

"""
Generate a Q matrix for a `NucleicAcidSubstitutionModel`, of the form:

\$Q = \\begin{bmatrix}
Q_{A, A} & Q_{A, C} & Q_{A, G} & Q_{A, T} \\\\
Q_{C, A} & Q_{C, C} & Q_{C, G} & Q_{C, T} \\\\
Q_{G, A} & Q_{G, C} & Q_{G, G} & Q_{G, T} \\\\
Q_{T, A} & Q_{T, C} & Q_{T, G} & Q_{T, T} \\end{bmatrix}\$
"""
@inline function Q(mod::NASM)
    α = _α(mod)
    β = _β(mod)
    γ = _γ(mod)
    δ = _δ(mod)
    ϵ = _ϵ(mod)
    η = _η(mod)
    πA = _πA(mod)
    πC = _πC(mod)
    πG = _πG(mod)
    πT = _πT(mod)

    return Qmatrix(-(δ * πC + η * πG + β * πT), δ * πA, η * πA, β * πA,
                   δ * πC, -(δ * πA + ϵ * πG + α * πT), ϵ * πC, α * πC,
                   η * πG, ϵ * πG, -(η * πA + ϵ * πC + γ * πT), γ * πG,
                   β * πT, α * πT, γ * πT, -(β * πA + α * πC + γ * πG))
end

@inline function P_generic(mod::NASM, t::Float64)
    if t < 0.0
        error("t must be positive")
    end
    return exp(Q(mod) * t)
end

function P_generic(mod::NASM, t::Array{Float64})
    if any(t .< 0.0)
        error("t must be positive")
    end
    try
        eig_vals, eig_vecs = eig(Q(mod))
        return [eig_vecs * exp(diagm(eig_vals)*i) * eig_vecs' for i in t]
    catch
        eig_vals, eig_vecs = eig(Array(Q(mod)))
        return [SMatrix(eig_vecs * exp(diagm(eig_vals)*i) * eig_vecs') for i in t]
    end
end

"""
Generate a P matrix for a `NucleicAcidSubstitutionModel`, of the form:

\$P = \\begin{bmatrix}
P_{A, A} & P_{A, C} & P_{A, G} & P_{A, T} \\\\
P_{C, A} & P_{C, C} & P_{C, G} & P_{C, T} \\\\
P_{G, A} & P_{G, C} & P_{G, G} & P_{G, T} \\\\
P_{T, A} & P_{T, C} & P_{T, G} & P_{T, T} \\end{bmatrix}\$

for specified time
"""
P(mod::NASM, t) = P_generic(mod, t)


text:
module Scenarios

export Scenario

immutable Scenario{AgentState, Action}
    AGENTS_NO::Int64
    REWARDERS_NO::Int64
    init::Function
    perceive::Function
    validActions::Array{Action, 1}
    doActions!::Function
    neighbours::Function
end

end

max_stars_repo_name: tudor-berariu/phd

max_issues_repo_name: tudor-berariu/phd

max_forks_repo_name: tudor-berariu/phd


text:
using RayTracer
using LinearAlgebra
using Test

@testset "creating_4_x_4_matrix" begin
    m = [1 2 3 4; 5.5 6.5 7.5 8.5; 9 10 11 12; 13.5 14.5 15.5 16.5]
    @test m[1,1] ≈ 1
    @test m[1,4] ≈ 4
    @test m[2,1] ≈ 5.5
    @test m[2,3] ≈ 7.5
    @test m[3,3] ≈ 11
    @test m[4,1] ≈ 13.5
    @test m[4,3] ≈ 15.5
end

@testset "equal_matrices" begin
    m = [1 2 3 4; 5.5 6.5 7.5 8.5; 9 10 11 12; 13.5 14.5 15.5 16.5]
    n = [1 2 3 4; 5.5 6.5 7.5 8.5; 9 10 11 12; 13.5 14.5 15.5 16.5]
    @test m ≈ n
end

@testset "different_matrices" begin
    m = [1 2 3 4; 5 6 7 8; 9 10 11 12; 13 14 15 16]
    n = [2 3 4 5; 6 7 8 9; 10 11 12 13; 14 15 16 1]
    @test m ≠ n
end

@testset "matrix_tuple_multiplication" begin
    a = [1 2 3 4; 2 4 4 2; 8 6 4 1; 0 0 0 1]
    b = RayTracer.tuple(1, 2, 3, 1)
    @test a*b ≈ RayTracer.tuple(18,24,33,1)
end

@testset "creating_submatix" begin
    a = [1 5 0; -3 2 7; 0 6 -3]
    b = RayTracer.submatrix(a,1,3)
    @test b ≈ [-3 2; 0 6]
end

@testset "computing minor" begin
    a = [3 5 0; 2 -1 -7; 6 -1 5]
    b = RayTracer.submatrix(a,2,1)
    @test det(b) ≈ 25
    @test RayTracer.minor(a,2,1) ≈ 25
end

@testset "computing_cofactor" begin
    a = [3 5 0; 2 -1 -7; 6 -1 5]
    @test RayTracer.minor(a,1,1) ≈ -12
    @test RayTracer.cofactor(a,1,1) ≈ -12
    @test RayTracer.minor(a,2,1) ≈ 25
    @test RayTracer.cofactor(a,2,1) ≈ -25
end

@testset "multiplying_by_translation_matrix" begin
    transform = RayTracer.translation(5,-3,2)
    p = RayTracer.point(-3,4,5)
    @test transform * p ≈ RayTracer.point(2,1,7)
end

@testset "multiplying_by_inverse_translation" begin
    transform = RayTracer.translation(5,-3,2)
    inv_transform = inv(transform)
    p = RayTracer.point(-3,4,5)
    @test inv_transform * p ≈ RayTracer.point(-8,7,3)
end

@testset "translation_does_not_affect_vectors" begin
    transform = RayTracer.translation(5,-3,2)
    v = RayTracer.vector(-3,4,5)
    @test transform*v ≈ v
end

@testset "scaling_applied_to_point" begin
    transform = RayTracer.scaling(2,3,4)
    p = RayTracer.point(-4,6,8)
    @test transform*p ≈ RayTracer.point(-8,18,32)
end

@testset "scaling_applied_to_vector" begin
    transform = RayTracer.scaling(2,3,4)
    v = RayTracer.vector(-4,6,8)
    @test transform*v ≈ RayTracer.vector(-8,18,32)
end

@testset "multiplying_inverse_scaling" begin
    transform = RayTracer.scaling(2,3,4)
    inv_transform = inv(transform)
    v = RayTracer.vector(-4,6,8)
    @test inv_transform*v ≈ RayTracer.vector(-2,2,2)
end

@testset "reflection_is_scaling_negative" begin
    transform = RayTracer.scaling(-1,1,1)
    p = RayTracer.point(2,3,4)
    @test transform*p ≈ RayTracer.point(-2,3,4)
end

@testset "rotation_around_x_axis" begin
    p = RayTracer.point(0,1,0)
    eith_turn = RayTracer.rotation_x(π/4)
    quarter_turn = RayTracer.rotation_x(π/2)
    @test eith_turn*p ≈ RayTracer.point(0,√2/2,√2/2)
    @test ≈(quarter_turn*p,RayTracer.point(0,0,1),1e-10)
end

@testset "inverse_rotation_x_axis" begin
    p = RayTracer.point(0,1,0)
    eith_turn = RayTracer.rotation_x(π/4)
    inv_transform = inv(eith_turn)
    @test inv_transform*p ≈ RayTracer.point(0,√2/2,-√2/2)
end

@testset "rotation_around_y_axis" begin
    p = RayTracer.point(0,0,1)
    eith_turn = RayTracer.rotation_y(π/4)
    quarter_turn = RayTracer.rotation_y(π/2)
    @test eith_turn*p ≈ RayTracer.point(√2/2,0,√2/2)
    @test ≈(quarter_turn*p,RayTracer.point(1,0,0),1e-10)
end

@testset "rotation_around_z_axis" begin
    p = RayTracer.point(0,1,0)
    eith_turn = RayTracer.rotation_z(π/4)
    quarter_turn = RayTracer.rotation_z(π/2)
    @test eith_turn*p ≈ RayTracer.point(-√2/2,√2/2,0)
    @test ≈(quarter_turn*p,RayTracer.point(-1,0,0),1e-10)
end

@testset "shearing_x_wrt_y" begin
    transform = RayTracer.shearing(1,0,0,0,0,0)
    p = RayTracer.point(2,3,4)
    @test transform*p ≈ RayTracer.point(5,3,4)
end

@testset "shearing_x_wrt_z" begin
    transform = RayTracer.shearing(0,1,0,0,0,0)
    p = RayTracer.point(2,3,4)
    @test transform*p ≈ RayTracer.point(6,3,4)
end

@testset "shearing_y_wrt_x" begin
    transform = RayTracer.shearing(0,0,1,0,0,0)
    p = RayTracer.point(2,3,4)
    @test transform*p ≈ RayTracer.point(2,5,4)
end

@testset "shearing_y_wrt_z" begin
    transform = RayTracer.shearing(0,0,0,1,0,0)
    p = RayTracer.point(2,3,4)
    @test transform*p ≈ RayTracer.point(2,7,4)
end

@testset "shearing_z_wrt_x" begin
    transform = RayTracer.shearing(0,0,0,0,1,0)
    p = RayTracer.point(2,3,4)
    @test transform*p ≈ RayTracer.point(2,3,6)
end

@testset "shearing_z_wrt_y" begin
    transform = RayTracer.shearing(0,0,0,0,0,1)
    p = RayTracer.point(2,3,4)
    @test transform*p ≈ RayTracer.point(2,3,7)
end

@testset "individual_transformations_sequential" begin
    p = RayTracer.point(1,0,1)
    A = RayTracer.rotation_x(π/2)
    B = RayTracer.scaling(5,5,5)
    C = RayTracer.translation(10,5,7)
    p2 = A*p
    @test ≈(p2,RayTracer.point(1,-1,0),1e-10)
    p3 = B*p2
    @test ≈(p3,RayTracer.point(5,-5,0),1e-10)
    p4 = C*p3
    @test ≈(p4,RayTracer.point(15,0,7),1e-10)
end

@testset "chained_transformations" begin
    p = RayTracer.point(1,0,1)
    A = RayTracer.rotation_x(π/2)
    B = RayTracer.scaling(5,5,5)
    C = RayTracer.translation(10,5,7)
    D = C*B*A
    @test ≈(D*p,RayTracer.point(15,0,7),1e-10)
end


text:
module PuzzleMaker

export
    MotifStitcher,
    RandomMotif,
    add_motif!,
    add_motifs!,
    add_motif_arrangement!,
    generate_random_motif_sequences!,
    make_puzzle,
    SamplerWeighted

using BioSequences
using FASTX
using GenomeGraphs
import Random

const MUTATIONS = Dict(
    DNA_A => (DNA_C, DNA_G, DNA_T),
    DNA_C => (DNA_A, DNA_G, DNA_T),
    DNA_G => (DNA_A, DNA_C, DNA_T),
    DNA_T => (DNA_A, DNA_C, DNA_G),
)
const DEFAULT_SAMPLER = SamplerWeighted(dna"ATCG", fill(0.25, 3))
const SEQ_T = LongSequence{DNAAlphabet{2}}

struct RandomMotif
    len::Int
    sampler::SamplerWeighted{DNA}
end

struct SiblingMotif
    source_motif::Int
    poly::Float64
end

mutable struct MotifStitcher
    n_motifs::Int
    random_motifs::Dict{Int,RandomMotif}
    fixed_motifs::Dict{Int,SEQ_T}
    sibling_motifs::Dict{Int,SiblingMotif}
    motif_order::Vector{Vector{Int}}
end

function MotifStitcher()
    return MotifStitcher(
        0,
        Dict{Int,RandomMotif}(),
        Dict{Int,SEQ_T}(),
        Dict{Int,SiblingMotif}(),
        Vector{Vector{Int}}(),
    )
end

@inline function add_motif!(ms::MotifStitcher, rm::RandomMotif)
    ms.n_motifs += 1
    ms.random_motifs[ms.n_motifs] = rm
    return ms
end

@inline function add_motif!(ms::MotifStitcher, len::Int)
    ms.n_motifs += 1
    ms.random_motifs[ms.n_motifs] = RandomMotif(len, DEFAULT_SAMPLER)
    return ms
end

@inline function add_motif!(ms::MotifStitcher, fm::SEQ_T)
    ms.n_motifs += 1
    ms.fixed_motifs[ms.n_motifs] = fm
    return ms
end

@inline function add_motif!(ms::MotifStitcher, fm::LongDNASeq)
    return add_motif!(ms, LongSequence{DNAAlphabet{2}}(fm))
end

#=
function fix_motif_sequence!(ms::MotifStitcher, motifid::Int, seq::SEQ_T)
    delete!(ms.random_motifs, motifid)
    ms.fixed_motifs[motifid] = seq
    ms.motif_lengths[motifid] = length(seq)
    return ms
end
=#

@inline function add_motif!(ms::MotifStitcher, sm::SiblingMotif)
    ms.n_motifs += 1
    ms.sibling_motifs[ms.n_motifs] = sm
    return ms
end

@inline function add_motif!(ms::MotifStitcher, p::Pair{Int,Float64})
    return add_motif!(ms, SiblingMotif(first(p), last(p)))
end

function add_motifs!(ms::MotifStitcher, args...)
    for arg in args
        add_motif!(ms, arg)
    end
    return ms
end

function add_motif_arrangement!(ms::MotifStitcher, arr::Vector{Int})
    arr′ = abs.(arr)
    @assert all(1 ≤ x ≤ ms.n_motifs for x in arr′)
    push!(ms.motif_order, copy(arr))
    return ms
end

#=
function _delete_and_rekey_dictionary!(d::Dict{Int,T}, key::Int) where {T}
    rekeyed = Dict(k-1 => v for (k, v) in d if k > key)
    for k in keys(d)
        if k >= key
            delete!(d, k)
        end
    end
    merge!(d, rekeyed)
    return d
end

function delete_motif!(ms::MotifStitcher, motifid::Int)
    deleteat!(ms.motif_lens, motifid)
    _delete_and_rekey_dictionary!(ms.random_motifs, motifid)
    _delete_and_rekey_dictionary!(ms.fixed_motifs, motifid)
    _delete_and_rekey_dictionary!(ms.motif_similarities, motifid)
    # TODO: Check the motif similarities, if any of them are based off the deleted motif,
    # issue a warning and convert it to a random sequence.
    return ms
end
=#

#=
function set_motif_lengths!(pm::MotifStitcher, lens::Vector{Int})
    @assert all(lens .> 0)
    resize!(pm.motif_lens, length(lens))
    copyto!(pm.motif_lens, lens)
    return pm
end
=#

function generate_motif_sequences(ms::MotifStitcher)
    sequences = Dict{Int, SEQ_T}()
    for (k, v) in ms.random_motifs
        sequences[k] = randseq(DNAAlphabet{2}(), v.sampler, v.len)
    end
    merge!(sequences, ms.fixed_motifs)
    # TODO: Calls rand a lot.
    #for (k, v) in keys(ms.sibling_motifs
    for k in sort([keys(ms.sibling_motifs)...])
        v = ms.sibling_motifs[k]
        s = copy(sequences[v.source_motif])
        npoly = Int(ceil(v.poly * length(s)))
        sites = view(Random.randperm(length(s)), 1:npoly)
        @inbounds for si in sites
            s[si] = rand(MUTATIONS[s[si]])
        end
        sequences[k] = s
    end
    return sequences
end

function make_puzzle(ms::MotifStitcher)
    sequences = Vector{SEQ_T}(undef, length(ms.motif_order))
    motif_seqs = generate_motif_sequences(ms)
    for i in eachindex(ms.motif_order)
        s = SEQ_T()
        mo = ms.motif_order[i]
        for m in mo
            if m < 0
                append!(s, reverse_complement(motif_seqs[abs(m)]))
            else
                append!(s, motif_seqs[m])
            end
        end
        sequences[i] = s
    end
    return sequences
end

function make_puzzle(ms::MotifStitcher, fastafile::String)
    seqs = make_puzzle(ms)
    open(FASTA.Writer, fastafile) do wtr
        for (i, s) in enumerate(seqs)
            rec = FASTA.Record(string("haplotype_", i), s)
            write(wtr, rec)
        end
    end
    return seqs
end

end # module


text:
# Each new term in the Fibonacci sequence is generated by adding the previous
# two terms. By starting with 1 and 2, the first 10 terms will be:
#
# 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, ...
#
# By considering the terms in the Fibonacci sequence whose values do not exceed
# four million, find the sum of the even-valued terms.

using ProjectEulerSolutions

# Brute force solution, cycle through Fibonacci numbers up to 4 mil and sum the
# even ones.
function p002solution_running_sum(n::Integer=3)::Integer
    fibs = fibonacci_bound(BigInt(n))
    running_sum = 0
    @simd for val in fibs
        if val % 2 == 0
            running_sum += val
        end
    end
    return running_sum
end

# With a little thought, we see that it is only the 3*k+2 Fibonacci numbers that
# are even, so just sum them.
function p002solution_sum(n::Integer=3)::Integer
    fibs = fibonacci_bound(BigInt(n))
    return sum(fibs[2:3:length(fibs)])
end

p002 = Problems.Problem(Dict("Running sum" => p002solution_running_sum, "Sum" => p002solution_sum))

Problems.benchmark(p002, 4_000_000 - 1)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
import Pkg # stdlib
import PredictMDExtra
import Test # stdlib

package_list = PredictMDExtra._package_list()

for p in package_list
    try
        eval(
            Base.Meta.parse(
                string(
                    "import ",
                    p,
                    ),
                ),
            )
    catch e
        @debug("Ignoring exception: ", e,)
    end
end