PSQ_SearchForSpikes: Make it faster

[t-b]

From #2653 (comment):

Total time: 170.379, Time in Function code: 74.9018 (44%)
Top function percentages:
Function MIES_AnalysisFunctions_PatchSeq.ipf PSQ_SearchForSpikes: 44%
Function MIES_Cache.ipf CA_GetCacheIndex: 42%

Annotated Top Functions:

*******************************************************************************************
Function: MIES_AnalysisFunctions_PatchSeq.ipf PSQ_SearchForSpikes; Percent total 44%
*******************************************************************************************
[00]*         	|static Function/WAVE PSQ_SearchForSpikes(string device, variable type, WAVE sweepWave, variable headstage, variable offset, variable level, [variable searchEnd, variable numberOfSpikesReq, WAVE spikePositions, variable &numberOfSpikesFound])
[00]          	|
[00]          	|	variable first, last, overrideValue, rangeSearchLevel
[00]          	|	variable minVal, maxVal, numSpikesFoundOverride
[00]          	|	string msg
[00]          	|
[00]*         	|	WAVE spikeDetection = LBN_GetNumericWave()
[00]*         	|	spikeDetection = ((p == headstage) ? 0 : NaN)
[00]          	|
[00]*         	|	WAVE config = AFH_GetConfigWave(device, sweepWave)
[00]          	|
[00]*         	|	if(ParamIsDefault(searchEnd))
[00]*         	|		searchEnd = Inf
[00]*         	|	endif
[00]          	|
[00]*         	|	if(ParamIsDefault(numberOfSpikesReq))
[00]          	|		numberOfSpikesReq = 1
[00]          	|	else
[00]*         	|		ASSERT(numberOfSpikesReq > 0, "Invalid numberOfSpikesFound")
[00]*         	|		numberOfSpikesReq = trunc(numberOfSpikesReq)
[00]*         	|	endif
[00]          	|
[00]*         	|	if(!ParamIsDefault(numberOfSpikesFound))
[00]          	|		numberOfSpikesFound = NaN
[00]*         	|	endif
[00]          	|
[00]*         	|	sprintf msg, "Type %d, headstage %d, offset %g, numberOfSpikesReq %d", type, headstage, offset, numberOfSpikesReq
[00]*         	|	DEBUGPRINT(msg)
[00]          	|
[00]*         	|	if(type == PSQ_CHIRP || type == PSQ_TRUE_REST_VM)
[00]          	|		first = offset
[00]          	|		last  = searchEnd
[00]          	|	else
[00]          	|		// search pulse in DA and use the pulse as search region
[00]*         	|		WAVE singleDA = AFH_ExtractOneDimDataFromSweep(device, sweepWave, headstage, XOP_CHANNEL_TYPE_DAC, config = config)
[01]*         	|		[minVal, maxVal] = WaveMinAndMax(singleDA, offset, Inf)
[00]          	|
[00]*         	|		if(minVal == 0 && maxVal == 0)
[00]          	|			if(type == PSQ_SQUARE_PULSE)
[00]          	|				first = 0
[00]          	|				last  = searchEnd
[00]          	|			else
[00]          	|				return spikeDetection
[00]*         	|			endif
[00]          	|		else
[00]*         	|			rangeSearchLevel = minVal + GetMachineEpsilon(WaveType(singleDA))
[00]          	|
[00]*         	|			Make/FREE/D levels
[17]**        	|			FindLevels/R=(offset, Inf)/Q/N=2/DEST=levels singleDA, rangeSearchLevel
[00]*         	|			ASSERT(V_LevelsFound == 2, "Could not find two levels")
[00]*         	|			first = levels[0]
[00]          	|
[00]*         	|			if(type == PSQ_DA_SCALE)
[00]*         	|				last = levels[1]
[00]          	|			else
[00]*         	|				last = searchEnd
[00]          	|			endif
[00]          	|		endif
[00]*         	|	endif
[00]          	|
[00]*         	|	if(TestOverrideActive())
[00]          	|		WAVE overrideResults = GetOverrideResults()
[00]          	|		NVAR count           = $GetCount(device)
[00]          	|
[00]          	|		switch(type)
[00]          	|			case PSQ_TRUE_REST_VM:
[00]          	|				overrideValue          = overrideResults[0][count][1]
[00]          	|				numSpikesFoundOverride = overrideValue
[00]          	|				break
[00]          	|			case PSQ_CHIRP:
[00]          	|				overrideValue          = !overrideResults[0][count][3]
[00]          	|				numSpikesFoundOverride = overrideValue > 0
[00]          	|				break
[00]          	|			case PSQ_RHEOBASE:
[00]          	|				overrideValue          = overrideResults[0][count][1]
[00]          	|				numSpikesFoundOverride = overrideValue > 0
[00]          	|				break
[00]          	|			case PSQ_SQUARE_PULSE:
[00]          	|				overrideValue          = overrideResults[0][count][0]
[00]          	|				numSpikesFoundOverride = overrideValue > 0
[00]          	|				break
[00]          	|			case PSQ_RAMP:
[00]          	|				overrideValue          = overrideResults[0][count][1]
[00]          	|				numSpikesFoundOverride = overrideValue > 0
[00]          	|				break
[00]          	|			case PSQ_DA_SCALE:
[00]          	|				overrideValue = overrideResults[0][count][1]
[00]          	|				if(overrideValue > 0)
[00]          	|					numSpikesFoundOverride = overrideResults[0][count][2]
[00]          	|				else
[00]          	|					numSpikesFoundOverride = 0
[00]          	|				endif
[00]          	|				break
[00]          	|			default:
[00]          	|				FATAL_ERROR("unsupported type")
[00]          	|		endswitch
[00]          	|
[00]          	|		if(overrideValue == 0 || overrideValue == 1)
[00]          	|			spikeDetection[headstage] = overrideValue
[00]          	|		else
[00]          	|			spikeDetection[headstage] = overrideValue >= first && overrideValue <= last
[00]          	|		endif
[00]          	|
[00]          	|		if(!ParamIsDefault(spikePositions))
[00]          	|			ASSERT(WaveExists(spikePositions), "Wave spikePositions must exist")
[00]          	|			Redimension/D/N=(numSpikesFoundOverride) spikePositions
[00]          	|			spikePositions[] = overrideValue
[00]          	|		endif
[00]          	|
[00]          	|		if(!ParamIsDefault(numberOfSpikesFound))
[00]          	|			numberOfSpikesFound = numSpikesFoundOverride
[00]*         	|		endif
[00]          	|	else
[00]          	|
[00]*         	|		WAVE singleAD = AFH_ExtractOneDimDataFromSweep(device, sweepWave, headstage, XOP_CHANNEL_TYPE_ADC, config = config)
[00]*         	|		ASSERT(!cmpstr(WaveUnits(singleAD, -1), "mV"), "Unexpected AD Unit")
[00]          	|
[00]*         	|		if(type == PSQ_CHIRP)
[00]          	|			// use PA plot logic
[00]          	|			Duplicate/R=(first, last)/FREE singleAD, chirpChunk
[00]          	|
[00]          	|			ASSERT(numberOfSpikesReq == Inf, "Unexpected value of numberOfSpikesReq")
[00]          	|			WAVE spikePositionsResult = PA_SpikePositionsForNonVC(chirpChunk, level)
[00]          	|
[00]          	|			spikeDetection[headstage] = DimSize(spikePositionsResult, ROWS) > 0
[00]          	|
[00]          	|			if(!ParamIsDefault(numberOfSpikesFound))
[00]          	|				numberOfSpikesFound = DimSize(spikePositionsResult, ROWS)
[00]*         	|			endif
[00]          	|		else
[00]*         	|			if(numberOfSpikesReq == 1)
[00]          	|				// search the spike from the start till the end
[00]          	|				FindLevel/Q/R=(first, last)/B=3 singleAD, level
[00]          	|				spikeDetection[headstage] = !V_flag
[00]          	|
[00]          	|				if(!ParamIsDefault(spikePositions))
[00]          	|					ASSERT(WaveExists(spikePositions), "Wave spikePositions must exist")
[00]          	|					Redimension/D/N=(numberOfSpikesReq) spikePositions
[00]          	|					spikePositions[0] = V_LevelX
[00]          	|				endif
[00]          	|
[00]          	|				if(!ParamIsDefault(numberOfSpikesFound))
[00]          	|					numberOfSpikesFound = spikeDetection[headstage]
[00]*         	|				endif
[00]*         	|			elseif(numberOfSpikesReq > 1)
[00]*         	|				Make/D/FREE/N=0 crossings
[26]***       	|				FindLevels/Q/R=(first, last)/N=(numberOfSpikesReq)/DEST=crossings/EDGE=(FINDLEVEL_EDGE_INCREASING)/B=3 singleAD, level
[00]*         	|				spikeDetection[headstage] = IsFinite(numberOfSpikesReq) ? (!V_flag) : (V_LevelsFound > 0)
[00]          	|
[00]*         	|				if(!ParamIsDefault(spikePositions))
[00]*         	|					ASSERT(WaveExists(spikePositions), "Wave spikePositions must exist")
[00]*         	|					Redimension/D/N=(V_LevelsFound) spikePositions
[00]          	|
[00]*         	|					if(spikeDetection[headstage])
[00]*         	|						spikePositions[] = crossings[p]
[00]          	|					endif
[00]          	|				endif
[00]          	|
[00]*         	|				if(!ParamIsDefault(numberOfSpikesFound))
[00]          	|					numberOfSpikesFound = V_LevelsFound
[00]*         	|				endif
[00]          	|			else
[00]          	|				FATAL_ERROR("Invalid number of spikes value")
[00]*         	|			endif
[00]          	|		endif
[00]*         	|	endif
[00]          	|
[00]*         	|	ASSERT(IsFinite(spikeDetection[headstage]), "Expected finite result")
[00]          	|
[00]*         	|	return DEBUGPRINTw(spikeDetection)
[00]          	|End

We use PSQ_SearchForSpikes in the following analysis functions:

• DAScale Supra
• Square Pulse
• Rheobase
• Ramp
• Chirp
• TrueRestingMembranePotential

[t-b]

From https://aiephys.slack.com/archives/C07LLQEDKU0/p1775585919730959?thread_ts=1774986572.840799&cid=C07LLQEDKU0 we know that the function profiling was run on PSQ_Ramp sweeps and there we use it during midsweep.

[t-b]

Sweep 59 covers numerical values LBN row range 3320 to 3356 and the difference of the first and lasttimestamps is ~58s, the stimset length (acquired) in comparison is 12.7s. For PSQ_Ramp we call PSQ_SearchForSpikes during mid sweep event.
On my machine we call PSQ_SearchForSpikes in the PS_RA2 test RunWithOpts(testcase = "PS_RA2") around 300 times which means every 55ms.

Ideas:

• We could cache the search in the DAC data, but that needs also to work when we change the DA data as we do for PSQ_Ramp
• Or maybe even switch to use PSQ_GetPulseDurations(device, PSQ_RAMP, s.sweepNo, totalOnsetDelay) for the search range instead of AFH_ExtractOneDimDataFromSweep and FindLevels
• We can have a WaveMinMax check before using FindLevel, as if the maximum is below level there is also no level crossing
• We can also not call PSQ_SearchForSpikes that often but only every XXX ms No we don't want that
• We can also try to multithread find levels with overlapping windows

[t-b]

According to https://aiephys.slack.com/archives/C07LLQEDKU0/p1775588521703379?thread_ts=1774986572.840799&cid=C07LLQEDKU0 this issue is new and was not present in the old version Release_2.9_20250502.

We did not change anything drastic, see

$ git nicelog Release_2.9_20250502.. MIES_AnalysisFunctions_PatchSeq.ipf
fafe7f046c (PSQ_DeterminePulseDurationFinder: Adapt the search method for PSQ_RAMP, 2026-03-03)
be13cc2561 (treewide: Prefer constants for FindLevel EDGE flag, 2026-03-03)
34fe0085e9 (Fix: Fix naming of HeadstageContigency -> HeadstageContingency, 2025-04-10)
01ce91d3da (treewide: Add fallthrough statements, 2025-06-24)
08a86ed084 (treewide: Remove bogus fallthrough statements, 2025-06-25)
cf744860b1 (Packages/MIES: Unify case comment for no-break to fallthrough, 2025-06-24)
a85bdb3abc (treewide: Introduce FATAL_ERROR and SFH_FATAL_ERROR, 2025-06-09)
e444f3cd39 (Reformat with newer ipt version, 2025-06-05)
591b396bd0 (HW: Factor out readback of ITC fifoPos to HW_ITC_ReadFifoPos, 2025-05-02)

[t-b]

Investigate how long a plain FindLevels call takes on sweep 59 AD data.

[t-b]

@timjarsky Executing FindLevels/Q/N=(5)/DEST=crossings/EDGE=(FINDLEVEL_EDGE_INCREASING)/B=3 root:MIES:HardwareDevices:ITC18USB:Device0:Data:X_59:AD_0, PSQ_SPIKE_LEVEL is basically instantaneous here.

[MichaelHuth]

Just a thought: It might be related to the /R flag used.

[t-b]

Just a thought: It might be related to the /R flag used.

Here this is the same speedwise:

Function Dostuff()

	variable ref

	WAVE wv = root:MIES:HardwareDevices:ITC18USB:Device0:Data:X_59:AD_0

	ref = StopMSTimer(-2)
	FindLevels/Q/N=(5)/DEST=crossings/EDGE=(FINDLEVEL_EDGE_INCREASING)/B=3 wv, PSQ_SPIKE_LEVEL
	printf "%g\r", (StopMSTimer(-2) - ref) * MICRO_TO_ONE

	ref = StopMSTimer(-2)
	FindLevels/R=(0, inf)/Q/N=(5)/DEST=crossings/EDGE=(FINDLEVEL_EDGE_INCREASING)/B=3 wv, PSQ_SPIKE_LEVEL
	printf "%g\r", (StopMSTimer(-2) - ref) * MICRO_TO_ONE
End

gives

•dostuff()
  0.0074736
  0.0074717