Mapping along a dimension with findfirst, sometimes failing

General Usage
#1

I have a 3d array of records from a model (looking at the spread of disease resistance), the first two dimensions are for some parameter values, the third is for time so X[1,1,:] would be how the prevalence of disease resistance changes over time. I’d like a concise way to find the first time that the prevalence passes some critical value (note: it might not reach that value).

I currently have
t50 = mapslices(t -> findfirst(x -> x > 0.5, t), X, dims=3)
which does findfirst for the index of the first time the prevalence passes 0.5.

However, sometimes the prevalence doesn’t reach 0.5, and so I get an error:
ERROR: MethodError: Cannot `convert` an object of type Int64 to an object of type Nothing
which I think is because findfirst is returning nothing. It doesn’t return an array of type Union{Float64, Nothing}, which would be what I’d expect.

What I’d ideally like it to do is to specify a default value (e.g. the time of the final index), and just substitute in that, but I’m not sure how to handle the nothings. Right now I’m doing it manually with 3 nested messy for loops, but a more elegant solution would be nice, if anyone has any ideas?

#2 
mapslices((t -> something(findfirst(x -> x > 0.5, t), default_value), X, dims=3)
2 Likes
#3

This seems to be a bug in mapslices: https://github.com/JuliaLang/julia/issues/26232. In the future the solution (other than substituting a default) will probably be a multidimensional eachslice (see https://github.com/JuliaLang/julia/pull/31217#issuecomment-468782768).

In the meantime you can substitute a default, or if you do happen to want the nothings you could try

map(t->findfirst(x->x>0.9, t), (view(X,i,j,:) for i=1:size(X,1), j=1:size(X,2)))
1 Like
#4

mapslices unfortunately doesn’t support unions like that. Instead, you can use map over eachslice:

t50 = map(t -> findfirst(x -> x > 0.5, t), eachslice(X, dims=3))
1 Like
#5

@Chris_Foster @mbauman Thanks, both of those look like pretty decent solutions! I didn’t know about eachslice, that’s handy to have.

#6

It looks like my final solution is:

t50 = map(t -> findfirst(x -> x>0.5, t),
		  (view(X, i,j,:) for i in 1:size(X,1), j in 1:size(X,2))) .|>
	t -> times[something(t), end]

which seems pretty readable, and half the length of my loopy based code. eachslice would be preferable, but it doesn’t yet support eachslice(X, dims=(1,2)), maybe that will be a later version!

1 Like
#7

I think @Tamas_Papp’s solution might work best (be most efficient, and shortest) for your immediate problem where you know what the nothing will be turned into.

2 Likes
#8

Ah, I didn’t realise that something was an actual function! Yes, that is a better solution to what I wrote, thank you @Tamas_Papp!

1 Like
#9

Haha! Yes, it’s somehow oddly named, even though the name makes complete sense once you know what it does.

2 Likes