# Function speed up, reduce allocations

**URL:** <https://discourse.julialang.org/t/function-speed-up-reduce-allocations/71647>\
**Category:** Performance\
**Created:** [November 17, 2021, 10:51am UTC](https://discourse.julialang.org/t/function-speed-up-reduce-allocations/71647 "2021-11-17T10:51:01Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![david.brust](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david.brust/32/30891_2.png) [@david.brust](https://discourse.julialang.org/u/david.brust)\
**Post date:** [November 17, 2021, 10:51am UTC](https://discourse.julialang.org/t/function-speed-up-reduce-allocations/71647/1 "2021-11-17T10:51:01Z")

</div>

Hi,

I am working on a combinatorial problem called Block Design:

[Block design - Wikipedia](https://en.wikipedia.org/wiki/Block_design) ( specifically Pairwise balanced uniform designs (2-designs or BIBDs) )

I try to group N=25 samples (e.g. integers 1:25) into “blocks” (aka pools) of 5 samples each, where each pair of samples only occurs in one block together. The blocks should be chosen so that 5 of these blocks form a partition of the set of all 25 samples. A partition of the set of all samples into 5 blocks as described above call a “layer”, I try to find as many unique layers as possible. (in the above example there are 6 layers)

I implemented a backtracking algorithm to find solutions to the problem, which currently finds fesible solutions for N=9, N=16,& N=25, but takes prohibitively long for N=36.

Timing for N=25:

> @time BlockDesign.run()  
> 1.606897 seconds (2.06 M allocations: 307.621 MiB, 0.59% gc time)

I think the backtracking algorithm itself can be improved, but looking at the large number of allocations and the Profiler output pointing to the **feasible** (PD, M, M\_) function, this could be address first.  
(Is there a way to attach the complete profiler output as .txt?):

Function call highlighted by profiler:  
` fs = myfind(feasible(PD, M, M_))`

In short: how can we make the function **feasible** (PD, M, M\_) faster, either with the current implementation using the Int64 matrices M & M\_ or in another way.

The functions + helper functions needed to run above function call:

```julia
struct PoolDesign{T<:Integer}
    N::T
    m::T
    q::T
    lmax::T
end

PD = PoolDesign(
    25, # number samples
    5, # samples per pool
    5, # pools per layer
    6, # layers in pooling design
)

function myfind(c)
    a = similar(c, Int)
    count = 1
    for i in eachindex(c)
        a[count] = i
        count += (c[i] != zero(eltype(c)))
    end
    return resize!(a, count-1)
end

 function pool_idx(M, PD)
    for (i, row) in enumerate(eachrow(M))
        if sum(row) < PD.m
            return i
        end
    end
    size(M,1)
end

layer_idx(M, PD) = ((pool_idx(M, PD) - 1)÷PD.q) + 1

layer(PD, l) = @view M[((l-1)*PD.q+1):(l * PD.q), :]

current_layer(M, PD) = layer(PD, layer_idx(M, PD))

function feasible(PD, M, M_)
    feas = collect(1:PD.N)
    cl = current_layer(M, PD)
    # if current layer is empty, all samples are feasible
    if !iszero(cl) 
        for i in 1:PD.N
            for q in 1:PD.q
                if cl[q, i] == 1
                    feas[i] = 0
                end
            end
        end
        # if current pool is empty,
        # all remaining unplaced samples in layer are feasible
        if !iszero(current_pool(M, PD))
            for k in myfind(current_pool(M, PD))
                for l in 1:PD.N
                    feas[l] = feas[l] * M_[k, l]
                end
                feas[k] = 0
            end
        end
    end
    return feas
end

dims = (PD.lmax * PD.m, PD.N)
M = zeros(Int64, dims)
M_ = ones(Int64, PD.N, PD.N)

fs = myfind(feasible(PD, M, M_))

```

The matrix M stores the Block design, assigning the samples (cols) to the different blocks (rows), where M(j,i) = 1 means that sample i is block j.  
The matrix M\_ is symmetric and stores the information about which pairs of samples already occured together in one block.

Ex. solution for N=9:

> 1 1 1 0 0 0 0 0 0 → (1,2,3) |  
> 0 0 0 1 1 1 0 0 0 → (4,5,6) | Layer 1  
> 0 0 0 0 0 0 1 1 1 → (7,8,9) |  
> 1 0 0 1 0 0 1 0 0 → (1,4,7) ]  
> 0 1 0 0 1 0 0 1 0 → (2,5,8) ] Layer 2  
> 0 0 1 0 0 1 0 0 1 → (3,6,9) ]  
> 1 0 0 0 1 0 0 0 1 → (1,5,9) }  
> 0 1 0 0 0 1 1 0 0 → (2,6,7) } Layer 3  
> 0 0 1 1 0 0 0 1 0 → (3,4,8) }  
> 1 0 0 0 0 1 0 1 0 → (1,6,8) ||  
> 0 1 0 1 0 0 0 0 1 → (2,4,9) || Layer 4  
> 0 0 1 0 1 0 1 0 0 → (3,5,7) ||

If this post is in the wrong category, please feel free to move it appropriately.

Thanks  
David

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [November 17, 2021, 11:46am UTC](https://discourse.julialang.org/t/function-speed-up-reduce-allocations/71647/2 "2021-11-17T11:46:10Z")

</div>

It would simplify things considerably if you could post the MWE in one markup block. Then I wouldn’t need to manually compose it.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [November 17, 2021, 11:52am UTC](https://discourse.julialang.org/t/function-speed-up-reduce-allocations/71647/3 "2021-11-17T11:52:15Z")

</div>

> [@david.brust](#):
>
> ` for k in myfind(current_pool(M, PD))`

my guess just looking at the code is that you should probably expand this `myfind` function into this loop, to get rid of the array that is created by `myfind`, and just iterate over the elements that are necessary directly. That will make the code a little less modular, but at first sight it seems to be the cause of the allocations you are seeing.

That is, do something like:

```julia
for val in current_pool(M, PD)
    if val != zero(...)
        # ... do what you need to do
    end
end

```

---

<div class="post-metadata">

**Author:** ![david.brust](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david.brust/32/30891_2.png) [@david.brust](https://discourse.julialang.org/u/david.brust)\
**Post date:** [November 17, 2021, 12:17pm UTC](https://discourse.julialang.org/t/function-speed-up-reduce-allocations/71647/4 "2021-11-17T12:17:01Z")

</div>

Thanks for your answer. I also use the `myfind()` function in other places of the code to get the non-zero elements from the passed arrays.  
When ran with `myfind()` I get:

```julia
@time myfind(feasible(PD, M, M_))
  0.000012 seconds (7 allocations: 768 bytes)

```

Compared to:

```julia
@time feasible(PD, M, M_)
  0.000010 seconds (6 allocations: 480 bytes)

```

There are already allocations occuring in `feasible(PD, M, M_)`, where I tried to set the elements zero, and filter the zero elements out later in `myfind()`.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [November 17, 2021, 12:19pm UTC](https://discourse.julialang.org/t/function-speed-up-reduce-allocations/71647/5 "2021-11-17T12:19:51Z")

</div>

`myfind` is called from within `feasible`, isn’t it?

(by the way, you can write that myfind with `findall(!isequal(0),x)`, as far as I could understand it, but that won’t be any better)

---

<div class="post-metadata">

**Author:** ![david.brust](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david.brust/32/30891_2.png) [@david.brust](https://discourse.julialang.org/u/david.brust)\
**Post date:** [November 17, 2021, 12:30pm UTC](https://discourse.julialang.org/t/function-speed-up-reduce-allocations/71647/6 "2021-11-17T12:30:05Z")

</div>

Yes it is, sorry for the confusion.

I got now:

```julia
function feasible(PD, M, M_)
    feas = collect(1:PD.N)

    cl = current_layer(M, PD)
    # if current layer is empty, all samples are feasible
    if !iszero(cl)
                
        for i in 1:PD.N
            for q in 1:PD.q
                if cl[q, i] == 1
                    feas[i] = 0
                end
            end
            
        end

        # if current pool is empty,
        # all remaining unplaced samples in layer are feasible
        if !iszero(current_pool(M, PD))

            for k in current_pool(M, PD)
                if k != zero(k)
                    for l in 1:PD.N
                        feas[l] = feas[l] * M_[k, l]
                    end
                    feas[k] = 0
                end
            end
        end
    end
    return feas
end

@time feasible(PD, M, M_)
0.000011 seconds (6 allocations: 480 bytes)

```

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [November 17, 2021, 12:33pm UTC](https://discourse.julialang.org/t/function-speed-up-reduce-allocations/71647/7 "2021-11-17T12:33:36Z")

</div>

> [@david.brust](#):
>
> `layer(PD, l) = @view M[((l-1)*PD.q+1):(l * PD.q), :]`

I noticed now that `M` is not an input of `layer`. Probably it should be.

---

<div class="post-metadata">

**Author:** ![david.brust](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david.brust/32/30891_2.png) [@david.brust](https://discourse.julialang.org/u/david.brust)\
**Post date:** [November 17, 2021, 12:43pm UTC](https://discourse.julialang.org/t/function-speed-up-reduce-allocations/71647/8 "2021-11-17T12:43:16Z")

</div>

Good eye! This did the trick…  
It is down to

```julia
@time feasible(PD, M, M_)
  0.000004 seconds (1 allocation: 160 bytes)

```

Is the last allocation due to

```julia
 feas = collect(1:PD.N)

```

?

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [November 17, 2021, 12:43pm UTC](https://discourse.julialang.org/t/function-speed-up-reduce-allocations/71647/9 "2021-11-17T12:43:40Z")

</div>

That certainly is one allocation.

---

<div class="post-metadata">

**Author:** ![david.brust](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david.brust/32/30891_2.png) [@david.brust](https://discourse.julialang.org/u/david.brust)\
**Post date:** [November 17, 2021, 12:52pm UTC](https://discourse.julialang.org/t/function-speed-up-reduce-allocations/71647/10 "2021-11-17T12:52:37Z")

</div>

Since `feasible(PD, M, M_)` is called frequently, I will make a preallocated array `feas` that will get passed as a parameter and modified in the function. Thanks for your help!

---

<div class="post-metadata">

**Author:** ![nhavt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nhavt/32/16356_2.png) [@nhavt](https://discourse.julialang.org/u/nhavt)\
**Post date:** [November 17, 2021, 2:47pm UTC](https://discourse.julialang.org/t/function-speed-up-reduce-allocations/71647/11 "2021-11-17T14:47:35Z")

</div>

Hi, if you want to construct large balanced incomplete block designs (BIBDs), it might be good to have look at this work ([Memetic collaborative approaches for finding balanced incomplete block designs - ScienceDirect](https://www.sciencedirect.com/science/article/abs/pii/S0305054819302461)), where they proposed several algorithms to construct BIBDs. In your example it seems that you want to construct a resolvable balanced incomplete block design with v = 25, k =5, b = 5, r = 6?

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [November 17, 2021, 4:03pm UTC](https://discourse.julialang.org/t/function-speed-up-reduce-allocations/71647/12 "2021-11-17T16:03:16Z")

</div>

> [@lmiq](#):
>
> findall(!isequal(0),x)

Or, preferably, `findall(!iszero, x)`

---

<div class="post-metadata">

**Author:** ![david.brust](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david.brust/32/30891_2.png) [@david.brust](https://discourse.julialang.org/u/david.brust)\
**Post date:** [November 17, 2021, 8:27pm UTC](https://discourse.julialang.org/t/function-speed-up-reduce-allocations/71647/13 "2021-11-17T20:27:53Z")

</div>

Thanks for suggesting this paper, it is very relevant!

> you want toconstruct a resolvable balanced incomplete block design with v = 25, k =5, b = 5, r = 6?

Yes, exactly.