HolyLab / ImagineInterface.jl

Read and write Imagine analog and digital recordings and commands
Other
2 stars 2 forks source link

Slow append!() when input is large #30

Closed HaoyangRong closed 7 years ago

HaoyangRong commented 7 years ago

Hi Cody, append!() seems to be pretty slow when the input vector is large as we talked. The time consumption scales near linearly with the size of the input vector. I speculate compress!() might be the main cause. Here is my code:

using ImagineInterface
import Unitful: μm, s, V

#  Convert valve ID into a (stim_len x sample_rate)x5 boolean matrix
function ValveToSeq(valve)
  return repmat(ValveToStim(valve),stim_len*sample_rate,1)
end

# Convert valve ID into a 5-digit binary number
function ValveToStim(valve)
  if valve in 1:16
    EN="1"
    As=bin(valve-1,4)
    commd_s=string(As,EN)
    commd_b=map(x -> x=='1',collect(commd_s))
  else
    commd_b=falses(n_chn)
  end
  return commd_b'
end

########################Define parameters here#########################
sample_rate = 5000s^-1

ValveIDs=[1 2 3 5 9 11][:] # function shuffle requires a flat array
n_stimID=length(ValveIDs) # total number of stimuli IDs

stim_len= 1s # stimulus duration in each trial

ITI=60s # inter-trial-interval: gap between the 1st stim offset to 2nd stim onset

n_blck=2 # total number of blocks
IBI=120s # inter-block-interval

prestim=10s # time before 1st stimulus
poststim=10s # time after last stimulus

n_chn=5 # number of stimulus channels

##
############# Generate an array which holds stim sequence for the entire experiment ################

stim_pro=[shuffle(ValveIDs) for itr=1:n_blck]
stim_pro=hcat(stim_pro...)

######################Construct the stimuli matrix#########################
##
tic()

t_blck=n_stim*stim_len+(n_stim-1)*ITI # length of a block
t_total=prestim+t_blck*n_blck+IBI*(n_blck-1)+poststim # total time of the experiment
n_pts=t_total*sample_rate # total number of time points
# AllStim=falses(n_chn,n_pts)
AllStim=falses(n_pts,n_chn) # preallocate a matrix for all stimuli

cur_idx=prestim*sample_rate
next_idx=cur_idx+stim_len*sample_rate
for blck_id=1:n_blck
  for stim_itr=1:n_stim
    AllStim[cur_idx+1:next_idx,:]=ValveToSeq(stim_pro[stim_itr,blck_id])
    if stim_itr<n_stim # if the current stim is not the last one in this block
      cur_idx=next_idx+ITI*sample_rate # then cur_idx will be incremented by ITI
    else
      cur_idx=next_idx+IBI*sample_rate
    end
    next_idx=cur_idx+stim_len*sample_rate
  end
end
toc()

################## Append stimuli sequences to corresponding ImagineCommands#######
tic()
rig = "ocpi-2"
ocpi2 = rigtemplate(rig; sample_rate = sample_rate)
stims=getstimuli(ocpi2)
A0=stims[1]
A1=stims[2]
A2=stims[3]
A3=stims[4]
EN=stims[5]
# @profile append!(A0,"A0",AllStim[:,1])
# ProfileView.view()
append!(A0,"A0",AllStim[:,1])
append!(A1,"A1",AllStim[:,2])
append!(A2,"A2",AllStim[:,3])
append!(A3,"A3",AllStim[:,4])
append!(EN,"EN",AllStim[:,5])
toc()
Cody-G commented 7 years ago

Thanks for the report! The main problem was that this function had a type instability. It produces the conversion function that gets executed on every sample, so it caused a major slowdown. After fixing that I get a 10x improvement in speed.

I also realized that for digital signals we were doing two more mapping steps than necessary because the raw units (UInt8 true/false) are the same units as the vector that you're feeding in. After fixing that I got another 2x, for a total of 20x speed increase. Let me know how that works for you. If it's not fast enough there is still more that we can do.

See PR #31

HaoyangRong commented 7 years ago

With sample_rate=50000s^-1, n_blck=2, n_stimID=6 the original version took \~90s to append all stimuli, but it took only \~7s after your update. I noticed a significant increase in speed! When it comes to writing stimulus profile for a real experiment n_blck ranges from 4 to 6 and n_stimID\~=15, it took about 35s or more. Considering the full waveform file would include up to 5 additional channels for dual-color imaging (piezo, 2 cameras, 2 lasers), I hope the speed can be further increased if possible. The stimulus waveforms are pretty sparse, e.g. they only have a small portion of trues here and there. Maybe that's something to be utilized? I'm looking into the generation of other control commands. I'll keep you updated.

Cody-G commented 7 years ago

So my thought was that as a next step we can create a version of append! that accepts vectors that are already compressed, i.e. instances of RLEVector. That means it would be the responsibility of the user to do the compression beforehand. I don't think this would complicate things much for the user. An RLEVector is just a vector with elements that are RepeatedValues Here's an example of how you could construct a stim on, stim off sequence:


nsamps_on = 100
nsamps_off = 200
stim_vec = ImagineInterface.RepeatedValue{Bool}[]
stim_on = ImagineInterface.RepeatedValue(nsamps_on, true)
stim_off = ImagineInterface.RepeatedValue(nsamps_off, false)
push!(stim_vec, stim_on)
push!(stim_vec, stim_off)

#currently this isn't implemented, but I will add it
#(assume stim_com1 is an ImagineSignal for a stimulus command)
append!(stim_com1, "on_off", stim_vec)

This will be much faster because it doesn't require iterating through the sample vector for compression. Now I will work on a PR that implements the above version of append! and also exports the RepeatedValue and RLEVector types so that you don't have to type ImagineInterface. to use them.

Note that this works for digital signals, but it's less useful for analog signals like the piezo that are almost always changing. But since the piezo signal is usually repeated every stack it should not cause a performance problem as long as you append one stack and then duplicate it.

Cody-G commented 7 years ago

So this is done, see #32

Here's an updated example script:

using ImagineInterface, Unitful
import Unitful:s

ocpi2 = rigtemplate("ocpi-2"; sample_rate = 50000s^-1)
stim_com1 = getstimuli(ocpi2)[1]

nsamps_on = 100000
nsamps_off = 200000
stim_vec = RepeatedValue{rawtype(stim_com1)}[]
stim_on = RepeatedValue(nsamps_on, true)
stim_off = RepeatedValue(nsamps_off, false)
push!(stim_vec, stim_on)
push!(stim_vec, stim_off)

append!(stim_com1, "on_off", stim_vec)

Using this approach stimuli signals can be generated instantly. I still recommend that for camera and piezo commands you reuse the signal for a single stack instead of doing it this way

HaoyangRong commented 7 years ago

That's awesome! I'll try it out.

HaoyangRong commented 7 years ago

The change resulted in a huge improvement in speed! Thanks!

Cody-G commented 7 years ago

Glad it helped!

On Wed, Aug 16, 2017 at 1:42 PM, HaoyangRong notifications@github.com wrote:

Closed #30 https://github.com/HolyLab/ImagineInterface/issues/30.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/HolyLab/ImagineInterface/issues/30#event-1208851006, or mute the thread https://github.com/notifications/unsubscribe-auth/AE78hJm_MYERnTKq-ArusyWENPOlFIXFks5sYzgXgaJpZM4OqZqn .

HaoyangRong commented 7 years ago

My complete stimulus sequences consists of individual pre-defined sequences arranged in an order. Sometimes, such an individual sequence can have a zero length, since its corresponding time length might be 0 due to the protocol.

append! has trouble working with zero-length sequence, and it seems to be sequence value dependent. Running code:

using ImagineInterface, Unitful
import Unitful:s

ocpi3 = rigtemplate("ocpi-2"; sample_rate = 5s^-1)
stim_com1 = getstimuli(ocpi3)[1]

nsamps_on = 10
nsamps_off = 20
stim_vec = RepeatedValue{rawtype(stim_com1)}[]
stim_on = RepeatedValue(nsamps_on, true)
stim_off = RepeatedValue(nsamps_off, false)
stim_0=RepeatedValue(0, false)
push!(stim_vec, stim_on)
push!(stim_vec,stim_0)
push!(stim_vec, stim_off)

append!(stim_com1, "on_off", stim_vec)

using Plots, ImaginePlots
plot(ocpi3)

I got image , which looks pretty good.

If I switch push!(stim_vec, stim_on) and push!(stim_vec,stim_0) (stim_0 is a zero length sequence):

using ImagineInterface, Unitful
import Unitful:s

ocpi3 = rigtemplate("ocpi-2"; sample_rate = 5s^-1)
stim_com1 = getstimuli(ocpi3)[1]

nsamps_on = 10
nsamps_off = 20
stim_vec = RepeatedValue{rawtype(stim_com1)}[]
stim_on = RepeatedValue(nsamps_on, true)
stim_off = RepeatedValue(nsamps_off, false)
stim_0=RepeatedValue(0, false)
push!(stim_vec,stim_0)
push!(stim_vec, stim_on)
push!(stim_vec, stim_off)

append!(stim_com1, "on_off", stim_vec)

using Plots, ImaginePlots
plot(ocpi3)

I got image It failed to append any true value to the stimulus channel.

If I keep the original order of push!(stim_vec, stim_on) and push!(stim_vec,stim_0), and instead I set stim_0=RepeatedValue(0, true)

using ImagineInterface, Unitful
import Unitful:s

ocpi3 = rigtemplate("ocpi-2"; sample_rate = 5s^-1)
stim_com1 = getstimuli(ocpi3)[1]

nsamps_on = 10
nsamps_off = 20
stim_vec = RepeatedValue{rawtype(stim_com1)}[]
stim_on = RepeatedValue(nsamps_on, true)
stim_off = RepeatedValue(nsamps_off, false)
stim_0=RepeatedValue(0, true)
push!(stim_vec, stim_on)
push!(stim_vec,stim_0)
push!(stim_vec, stim_off)

append!(stim_com1, "on_off", stim_vec)

using Plots, ImaginePlots
plot(ocpi3)

The same problem also happens: image

Cody-G commented 7 years ago

Is this still a problem?

Oops, didn't get the message. Reading it now.

On Sep 27, 2017 11:52 AM, "HaoyangRong" notifications@github.com wrote:

Reopened #30 https://github.com/HolyLab/ImagineInterface/issues/30.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/HolyLab/ImagineInterface/issues/30#event-1268194017, or mute the thread https://github.com/notifications/unsubscribe-auth/AE78hJvB0VhqLhpY2zIrCuQ8iTMnyMHrks5smn06gaJpZM4OqZqn .

Cody-G commented 7 years ago

Thanks for the detail, as always! This was actually a problem with get_samples. It had an assumption that each RLEVector has length >=1. I have a fix, will submit the PR now.

Cody-G commented 7 years ago

I'll merge #48 first since I made the fix on top of that.

Cody-G commented 7 years ago

Should be fixed, but please reopen if not.