Skip to content

Commit

Permalink
PSQ_DAScale: Fix DAScale estimation for PSQ_DS_ADAPT
Browse files Browse the repository at this point in the history 
When we have filled in AP frequencies via PSQ_DS_GatherOvershootCorrection
and don't have any futureDAScales left to measure, we need to calculate a
new DAscale value via PSQ_DS_CalculateDAScale.

But before this commit we did pass in the fI slope/offset, apfreq, DAScale from the just
filled in sweep, which is wrong. We need to pass in the data from the
sweep with the highest AP frequency to continue our search. To
complicate things even more we need to search backwards and use the most
recent sweep.

Bug present since feaf86f (PSQ_DaScale: Add new operation mode Adaptive
Suprathreshold, 2023-10-20).
  • Loading branch information
@t-b
t-b committed Aug 21, 2024
1 parent d94938e commit 11ce187
Show file tree
Hide file tree
Showing 2 changed files with 63 additions and 42 deletions.
64 changes: 43 additions & 21 deletions Packages/MIES/MIES_AnalysisFunctions_PatchSeq.ipf
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2604,14 +2604,10 @@ End
/// - Check if the resulting fit slope is smaller than `slopePercentage`
///
/// @retval futureDAScales future DAScale values including the historic ones, @see PSQ_DS_GatherFutureDAScalesAndFrequency
/// @retval fitOffset offset of the linear curve fit
/// @retval fitSlope slope of the linear curve fit
/// @retval DAScale current DAScale value
/// @retval apfreq current AP frequency value
static Function [WAVE futureDAScales, variable fitOffset, variable fitSlope, variable DAScale, variable apfreq] PSQ_DS_EvaluateAdaptiveThresholdSweep(string device, variable sweepNo, variable headstage, STRUCT PSQ_DS_DAScaleParams &cdp)
static Function [WAVE futureDAScales] PSQ_DS_EvaluateAdaptiveThresholdSweep(string device, variable sweepNo, variable headstage, STRUCT PSQ_DS_DAScaleParams &cdp)

string key, errMsg
variable maxSlope, validFit
variable maxSlope, validFit, fitOffset, fitSlope

WAVE textualValues = GetLBTextualValues(device)
WAVE numericalValues = GetLBNumericalValues(device)
Expand Down Expand Up @@ -2643,23 +2639,11 @@ static Function [WAVE futureDAScales, variable fitOffset, variable fitSlope, var
fitSlope = NaN
endif

if(WaveExists(apfreqs))
apfreq = apfreqs[Inf]
else
apfreq = NaN
endif

if(WaveExists(DAScales))
DAScale = DAScales[Inf]
else
DAScale = NaN
endif

[maxSlope, WAVE fitSlopesAll, WAVE DAScalesAll] = PSQ_DS_CalculateMaxSlopeAndWriteToLabnotebook(device, sweepNo, headstage, fitSlope)

PSQ_DS_CalculateReachedFinalSlopeAndWriteToLabnotebook(device, sweepNo, fitSlopesAll, DAScalesAll, fitSlope, maxSlope, cdp.slopePercentage, validFit)

return [futureDAScales, fitOffset, fitSlope, DAScale, apfreq]
return [futureDAScales]
End

/// @brief Determine from the POST_SWEEP_EVENT if the set is already finished
Expand Down Expand Up @@ -3222,6 +3206,43 @@ static Function [variable daScaleStepMinNorm, variable daScaleStepMaxNorm] PSQ_D
return [daScaleMinStepWidth, daScaleMaxStepWidth]
End

static Function [variable fitOffset, variable fitSlope, variable DAScale, variable apfreq] PSQ_DS_GetValuesOfLargestAPFreq(string device, variable sweepNo, variable headstage)

variable emptySCI, offset, i, numEntries, maxValue, maxLoc

[WAVE apfreqAll, emptySCI] = PSQ_DS_GetLabnotebookData(device, sweepNo, headstage, PSQ_DS_APFREQ)
ASSERT(!emptySCI, "Unexpected emptySCI")
[WAVE fitSlopeAll, emptySCI] = PSQ_DS_GetLabnotebookData(device, sweepNo, headstage, PSQ_DS_FI_SLOPE)
ASSERT(!emptySCI, "Unexpected emptySCI")
[WAVE fitOffsetAll, emptySCI] = PSQ_DS_GetLabnotebookData(device, sweepNo, headstage, PSQ_DS_FI_OFFSET)
ASSERT(!emptySCI, "Unexpected emptySCI")
[WAVE DAScaleAll, emptySCI] = PSQ_DS_GetLabnotebookData(device, sweepNo, headstage, PSQ_DS_DASCALE)
ASSERT(!emptySCI, "Unexpected emptySCI")

// get the largest value with the highest index
maxValue = -Inf
numEntries = DimSize(apfreqAll, ROWS)
for(i = numEntries - 1; i >= 0; i -= 1)
if(apfreqAll[i] > maxValue)
maxValue = apfreqAll[i]
maxLoc = i
endif
endfor

// highest AP frequency is at V_maxloc
apfreq = apfreqAll[maxLoc]
DAScale = DAScaleAll[maxLoc]

// find the rest by indexing relative from the end as we have fewer fitOffset/fitSlope entries in total
// but a 1:1 correspondence in the SCI data which is at the end
offset = DimSize(apfreqAll, ROWS) - maxLoc

fitOffset = fitOffsetAll[DimSize(fitOffsetAll, ROWS) - offset]
fitSlope = fitOffsetAll[DimSize(fitSlopeAll, ROWS) - offset]

return [fitOffset, fitSlope, DAScale, apfreq]
End

static Structure PSQ_DS_DAScaleParams
variable absFrequencyMinDistance, maxFrequencyChangePercent
variable daScaleStepMinNorm, daScaleStepMaxNorm, slopePercentage
Expand Down Expand Up @@ -3679,7 +3700,7 @@ Function PSQ_DAScale(device, s)
variable sweepPassed, setPassed, length, minLength, reachedFinalSlope, fitOffset, fitSlope, apfreq, enoughFIPointsPassedQC
variable minimumSpikeCount, maximumSpikeCount, daScaleModifierParam, measuredAllFutureDAScales, fallbackDAScaleRangeFac
variable sweepsInSet, passesInSet, acquiredSweepsInSet, multiplier, asyncAlarmPassed, supraStimsetCycle
variable daScaleStepMinNorm, daScaleStepMaxNorm, maxSlope, validFit
variable daScaleStepMinNorm, daScaleStepMaxNorm, maxSlope, validFit, emptySCI
string msg, stimset, key, opMode, offsetOp, textboxString, str, errMsg
variable daScaleOffset
variable finalSlopePercent = NaN
Expand Down Expand Up @@ -3955,7 +3976,7 @@ Function PSQ_DAScale(device, s)
cdp.daScaleStepMinNorm = daScaleStepMinNorm
cdp.daScaleStepMaxNorm = daScaleStepMaxNorm

[WAVE futureDAScales, fitOffset, fitSlope, DAScale, apfreq] = PSQ_DS_EvaluateAdaptiveThresholdSweep(device, s.sweepNo, s.headstage, cdp)
[WAVE futureDAScales] = PSQ_DS_EvaluateAdaptiveThresholdSweep(device, s.sweepNo, s.headstage, cdp)

// prepare wave for setting DAScale; DAScalesIndex is incremented when required below
WAVE/ZZ DAScales = futureDAScales
Expand Down Expand Up @@ -4158,6 +4179,7 @@ Function PSQ_DAScale(device, s)

if(sweepPassed)
if(measuredAllFutureDAScales)
[fitOffset, fitSlope, DAScale, apfreq] = PSQ_DS_GetValuesOfLargestAPFreq(device, s.sweepNo, s.headstage)

dascale = PSQ_DS_CalculateDAScale(cdp, fitOffset, fitSlope, DAScale, apfreq)
Make/FREE/D DAScaleNew = {dascale}
Expand Down
41 changes: 20 additions & 21 deletions Packages/tests/HardwareAnalysisFunctions/UTF_PatchSeqDAScale_Adapt.ipf
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1512,7 +1512,7 @@ static Function PS_DS_AD10_REENTRY([string str])
CHECK_EQUAL_WAVES(entries[%samplingPass], {1, 1, 1}, mode = WAVE_DATA)

CHECK_EQUAL_WAVES(entries[%futureDAScalesPass], {0, 0, 0}, mode = WAVE_DATA)
CHECK_EQUAL_WAVES(entries[%fiSlopeReachedPass], {0, 0, 0}, mode = WAVE_DATA)
CHECK_EQUAL_WAVES(entries[%fiSlopeReachedPass], {0, 1, 0}, mode = WAVE_DATA)
CHECK_EQUAL_WAVES(entries[%enoughFIPointsPass], {1, 1, 1}, mode = WAVE_DATA)
CHECK_EQUAL_WAVES(entries[%validSlopePass], {1, 1, 1}, mode = WAVE_DATA)
CHECK_EQUAL_WAVES(entries[%initialValidSlopePass], {1, NaN, NaN}, mode = WAVE_DATA)
Expand All @@ -1531,14 +1531,14 @@ static Function PS_DS_AD10_REENTRY([string str])
CHECK_EQUAL_WAVES(entries[%dascaleFromRhSuAd], DAScalesFromRhSuAd, mode = WAVE_DATA)
CHECK_EQUAL_WAVES(entries[%sweepPassFromRhSuAd], sweepPassedFRomRhSuAd, mode = WAVE_DATA)

Make/FREE/D maxSlopeRef = {9.374999999999962e-10, 3.124999999999997e-09, 6.510416666666681e-09}
Make/FREE/D fiSlopeRef = {9.374999999999962e-10, 3.124999999999997e-09, 6.510416666666681e-09}
Make/FREE/D fiOffsetRef = {-21.49999999999984, -129.9999999999999, -314.1666666666805}
Make/FREE/T futureDAScalesRef = {"4.96;5.44;", "4.96;5.44;5.6704;", "4.96;5.44;5.6704;5.82246399999999;"}
Make/FREE/D maxSlopeRef = {9.374999999999962e-10, 9.374999999999962e-10, 9.374999999999962e-10}
Make/FREE/D fiSlopeRef = {9.374999999999962e-10, -8.272058823073271e-10, 3.124999999999996e-10}
Make/FREE/D fiOffsetRef = {-21.49999999999984, 66.02941176244343, 30.16666666697915}
Make/FREE/T futureDAScalesRef = {"4.96;3.14666666656668;", "4.96;3.14666666656668;7.94666666656668;", "4.96;3.14666666656668;7.94666666656668;14.5466666665667;"}

Make/FREE/D fiSlopesFromRhSuAdRef = {1e-10, 2e-10, 3e-10}
Make/FREE/D fiOffsetsFromRhSuAdRef = {9, 7, 4}
Make/FREE/D DAScalesRef = {4.96, 5.440000000000003, 5.670400000000001}
Make/FREE/D DAScalesRef = {4.96, 3.146666666566679, 7.946666666566679}

CHECK_EQUAL_WAVES(entries[%maxSlope], maxSlopeRef, mode = WAVE_DATA, tol = 1e-24)
CHECK_EQUAL_WAVES(entries[%fiSlope], fiSlopeRef, mode = WAVE_DATA, tol = 1e-24)
Expand Down Expand Up @@ -1646,16 +1646,16 @@ static Function PS_DS_AD11_REENTRY([string str])
CHECK_EQUAL_WAVES(entries[%sweepPassFromRhSuAd], sweepPassedFRomRhSuAd, mode = WAVE_DATA)

Make/FREE/D maxSlopeRef = {3e-10, 3e-10}
Make/FREE/D fiSlopeRef = {1.744186046511628e-10, 1.014061654948618e-10}
Make/FREE/D fiOffsetRef = {9.023255813953488, 13.6181719848567}
Make/FREE/T futureDAScalesRef = {"6.29333333333333;11.224;", \
"6.29333333333333;11.224;"}
Make/FREE/D fiSlopeRef = {1.744186046511628e-10, 9.722222222222222e-11}
Make/FREE/D fiOffsetRef = {9.023255813953488, 13.88148148148148}
Make/FREE/T futureDAScalesRef = {"6.29333333333333;11.4361904761905;", \
"6.29333333333333;11.4361904761905;"}

// we do have three pairs in apFrequenciesFromRhSuAd/DAScalesFromRhSuAd but a neighboring duplicate
// so only two valid slopes and offsets
Make/FREE/D fiSlopesFromRhSuAdRef = {PSQ_DS_SKIPPED_FI_SLOPE, 1.5e-10, 3e-10}
Make/FREE/D fiOffsetsFromRhSuAdRef = {PSQ_DS_SKIPPED_FI_SLOPE, 8.5, 4}
Make/FREE/D DAScalesRef = {6.293333333333333, 11.224}
Make/FREE/D DAScalesRef = {6.293333333333333, 11.43619047619048}

CHECK_EQUAL_WAVES(entries[%maxSlope], maxSlopeRef, mode = WAVE_DATA, tol = 1e-24)
CHECK_EQUAL_WAVES(entries[%fiSlope], fiSlopeRef, mode = WAVE_DATA, tol = 1e-24)
Expand Down Expand Up @@ -1713,10 +1713,10 @@ static Function PS_DS_AD12([string str])

WAVE wv = PSQ_CreateOverrideResults(str, PSQ_TEST_HEADSTAGE, PSQ_DA_SCALE, opMode = PSQ_DS_ADAPT)

wv[][0][%APFrequency] = 25
wv[][1][%APFrequency] = 21 // future DAScale (8.4)
wv[][2][%APFrequency] = 21 // redoing future DAScale (8.4), passed f-I slope QC (1.) and got another future DAScale 5.95
wv[][3][%APFrequency] = 21.1 // passed f-I slope QC (2.)
wv[][0][%APFrequency] = 27
wv[][1][%APFrequency] = 20 // future DAScale (8.22)
wv[][2][%APFrequency] = 20 // redoing future DAScale (8.22), passed f-I slope QC (1.) and got another future DAScale 5.86
wv[][3][%APFrequency] = 20.1 // passed f-I slope QC (2.)

wv[][][%AsyncQC] = 1
wv[][][%BaselineQC] = 1
Expand Down Expand Up @@ -1761,14 +1761,13 @@ static Function PS_DS_AD12_REENTRY([string str])
CHECK_EQUAL_WAVES(entries[%dascaleFromRhSuAd], DAScalesFromRhSuAd, mode = WAVE_DATA)
CHECK_EQUAL_WAVES(entries[%sweepPassFromRhSuAd], sweepPassedFRomRhSuAd, mode = WAVE_DATA)

Make/FREE/D maxSlopeRef = {3.416149068322985e-10, 3.416149068322985e-10, 3.416149068322985e-10, 3.416149068322985e-10}
Make/FREE/D fiSlopeRef = {3.416149068322985e-10, 3.416149068322985e-10, -2.376451525789897e-10, -4.078983962177433e-12}
Make/FREE/D fiOffsetRef = {2.043478260869549, 2.043478260869549, 40.9697542533081, 21.34276443867624}
Make/FREE/T futureDAScalesRef = {"6.72;8.40318181818182;", "6.72;8.40318181818182;", "6.72;8.40318181818182;5.95159090909091;", "6.72;8.40318181818182;5.95159090909091;"}

Make/FREE/D maxSlopeRef = {4.037267080745342e-10, 4.037267080745342e-10, 4.037267080745342e-10, 4.037267080745342e-10}
Make/FREE/D fiSlopeRef = {4.037267080745342e-10, 4.037267080745342e-10, -4.648501752436268e-10, -4.232032580229025e-12}
Make/FREE/D fiOffsetRef = {-0.1304347826086989, -0.1304347826086989, 58.23793177637171, 20.348121140308}
Make/FREE/T futureDAScalesRef = {"6.72;8.22586153836154;", "6.72;8.22586153836154;", "6.72;8.22586153836154;5.86293076918077;", "6.72;8.22586153836154;5.86293076918077;"}
Make/FREE/D fiSlopesFromRhSuAdRef = {2e-10, 2e-10, 1e-10}
Make/FREE/D fiOffsetsFromRhSuAdRef = {8, 8, 10.5}
Make/FREE/D DAScalesRef = {6.719999999999999, 8.403181818181814, 8.403181818181814, 5.951590909090905}
Make/FREE/D DAScalesRef = {6.719999999999999, 8.225861538361537, 8.225861538361537, 5.86293076918077}

CHECK_EQUAL_WAVES(entries[%maxSlope], maxSlopeRef, mode = WAVE_DATA, tol = 1e-24)
CHECK_EQUAL_WAVES(entries[%fiSlope], fiSlopeRef, mode = WAVE_DATA, tol = 1e-24)
Expand Down

0 comments on commit 11ce187

Please sign in to comment.