# Attempting to customise broadcasting for BioSequences

**URL:** <https://discourse.julialang.org/t/attempting-to-customise-broadcasting-for-biosequences/64892>\
**Category:** General Usage\
**Tags:** question, broadcast, broadcasting\
**Created:** [July 19, 2021, 1:16am UTC](https://discourse.julialang.org/t/attempting-to-customise-broadcasting-for-biosequences/64892 "2021-07-19T01:16:12Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ward9250](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ward9250/32/42768_2.png) [@Ward9250](https://discourse.julialang.org/u/Ward9250)\
**Post date:** [July 19, 2021, 1:16am UTC](https://discourse.julialang.org/t/attempting-to-customise-broadcasting-for-biosequences/64892/1 "2021-07-19T01:16:12Z")

</div>

Hello! I’m currently trying to customise and improve how Broadcasting is done with BioSequences.jl - types. Specifically `LongSequence` at this moment, for a v3 release. And I could do with some help and advice doing it properly.

So, BioSequence types are not strictly arrays in that they don’t inherit from any abstract array type. But they implement and satisfy indexing, iteration etc. As if they were vectors - they have one dimension, they have the `IndexStyle` of `Base.IndexLinear`, they have `length`, `size` will return `(length(x),)` and `axes` will return `(OneTo(length(x)), )`

We want to customise broadcasting such that - in some cases, the output type is a LongSequence - but tbh most of the time it’s probably a vector/array.

So I want to start with a simple example: `a .== DNA_G`, where a is a `LongSequence`.

```julia
julia> using BioSequences
[Info: Precompiling BioSequences [7e6ae17a-c86d-528c-b3b9-7f778a29fe59]

julia> a = dna"ATCGATCGATCGATCGATCG"
20nt DNA Sequence:
ATCGATCGATCGATCGATCG

julia> a .== DNA_G
20-element BitVector:
 0
 0
 0
 1
 0
 0
 0
 1
 0
 0
 0
 1
 0
 0
 0
 1
 0
 0
 0
 1

julia> Base.broadcastable(a)
20-element Vector{DNA}:
 DNA_A
 DNA_T
 DNA_C
 DNA_G
 DNA_A
 DNA_T
 DNA_C
 DNA_G
 DNA_A
 DNA_T
 DNA_C
 DNA_G
 DNA_A
 DNA_T
 DNA_C
 DNA_G
 DNA_A
 DNA_T
 DNA_C
 DNA_G

```

So for the simple example of `a .== DNA_G`, we expect an array (well a bitarray) out, and we get it. BUT we can see broadcasting spends the time and memory to build a vector of `DNA` (as per `Base.broadcastable`), rather than simply index into `a` itself.

So reading the customising broadcasting docs, I can see I can do something about this by implementing `Base.broadcastable`, so I do: `Base.broadcastable(x::LongSequence) = x`.

At this point trying again, I get an error:

```julia
julia> a .== DNA_G
ERROR: MethodError: no method matching getindex(::LongSequence{DNAAlphabet{4}}, ::CartesianIndex{1})
Closest candidates are:
  getindex(::LongSequence, ::UnitRange{var"#s30"} where var"#s30"<:Integer) at /Users/bward/repos/github/BioJulia/BioSequences/src/longsequences/indexing.jl:39
  getindex(::BioSequence, ::Colon) at /Users/bward/repos/github/BioJulia/BioSequences/src/biosequence/indexing.jl:66
  getindex(::BioSequence, ::AbstractVector{Bool}) at /Users/bward/repos/github/BioJulia/BioSequences/src/biosequence/indexing.jl:43

```

Ok, this is a bit weird to me, because my sequence type implements linear indexing, not cartesian, and satisfies the traits and methods for linear indexing. Indeed if I check the broadcast style:

```julia
julia> Base.BroadcastStyle(typeof(a))
Base.Broadcast.DefaultArrayStyle{1}()

```

So it know’s LongSequence is a 1 dimensional thing.

Ok so out of interest I do this: `Base.getindex(x::LongSequence, i::Base.CartesianIndex{1}) = x[first(i.I)]`

Such that trying to index with a 1d cartesian index just passes straight to linear indexing, and try again…

It looks like it worked…

Before

```julia
julia> using BenchmarkTools

julia> using BioSequences
[Info: Precompiling BioSequences [7e6ae17a-c86d-528c-b3b9-7f778a29fe59]

julia> a = dna"ATCGATCGATCGATCGATCG"
20nt DNA Sequence:
ATCGATCGATCGATCGATCG

julia> @benchmark a .== DNA_G
BenchmarkTools.Trial: 10000 samples with 175 evaluations.
 Range (min … max): 610.503 ns … 11.402 μs ┊ GC (min … max): 0.00% … 94.11%
 Time (median): 620.303 ns ┊ GC (median): 0.00%
 Time (mean ± σ): 636.760 ns ± 335.831 ns ┊ GC (mean ± σ): 1.81% ± 3.24%

   ▁▆▆█▃                                                         
  ▃█████▅▃▃▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂ ▃
  611 ns Histogram: frequency by time 753 ns <

 Memory estimate: 336 bytes, allocs estimate: 6.

```

After

```julia
julia> using BioSequences
[Info: Precompiling BioSequences [7e6ae17a-c86d-528c-b3b9-7f778a29fe59]

julia> a = dna"ATCGATCGATCGATCGATCG"
20nt DNA Sequence:
ATCGATCGATCGATCGATCG

julia> using BenchmarkTools

julia> @benchmark a .== DNA_G
BenchmarkTools.Trial: 10000 samples with 186 evaluations.
 Range (min … max): 554.548 ns … 9.028 μs ┊ GC (min … max): 0.00% … 93.01%
 Time (median): 565.371 ns ┊ GC (median): 0.00%
 Time (mean ± σ): 579.382 ns ± 237.237 ns ┊ GC (mean ± σ): 1.21% ± 2.78%

  ▃▇██▇▆▅▄▃▃▂▁▁▁▁ ▂
  ████████████████████▇▇██▆▇▆▆▆▇▇▅▄▅▇▆▆▆▆▅▅▃▆▇▇▇▆▅▄▅▄▄▃▄▃▃▅▁▅▁▅ █
  555 ns Histogram: log(frequency) by time 738 ns <

 Memory estimate: 224 bytes, allocs estimate: 5.

```

So with broadcasting using the sequence itself, there’s -1 allocation occurring and it’s going slightly faster.

So at this point, I should ask - was what I did the right thing to do, or have I borked something that kinda works but not really? It still feels weird to me to be indexing a linear thing with the cartesian.

I suppose to move onto cases like `a .= DNA_A` or say `complement.(a)` which should return a sequence and not an array, I’ll be going more into the machinery and implementing a style, and possibly implementing the `copy` and `copyto!`, methods, and figuring out precedence for styles. Has anyone implemented broadcasting for a type that behaves like an array but isn’t exactly `<: AbstractArray{1}`, like BioSequences, that might provide me a good example?

Many thanks,  
Ben.