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 committed Sep 3, 2024
1 parent cc451f8 commit 2606254
Show file tree
Hide file tree
Showing 4 changed files with 65 additions and 44 deletions.
64 changes: 43 additions & 21 deletions Packages/MIES/MIES_AnalysisFunctions_PatchSeq.ipf
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
2 changes: 1 addition & 1 deletion Packages/MIES/MIES_Constants.ipf
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ Constant RESULTS_VERSION = 3
Constant PSQ_PIPETTE_BATH_VERSION = 4
Constant PSQ_ACC_RES_SMOKE_VERSION = 2
Constant PSQ_CHIRP_VERSION = 13
Constant PSQ_DA_SCALE_VERSION = 8
Constant PSQ_DA_SCALE_VERSION = 9
Constant PSQ_RAMP_VERSION = 6
Constant PSQ_RHEOBASE_VERSION = 5
Constant PSQ_SQUARE_PULSE_VERSION = 4
Expand Down
2 changes: 1 addition & 1 deletion Packages/doc/dot/patch-seq-dascale.dot
Original file line number Diff line number Diff line change
Expand Up @@ -209,7 +209,7 @@ QC and f-I slope QC passing?"];
n109 -> n100 [label=Yes];
n103 [label="Measured all future DAScale values?"];
n100 -> n103 [label=Yes];
n106 [label="Calculate new DAScale value by\n extrapolating fit slope and offset\n and using MaximumChangePercent - 2 as frequency distance"];
n106 [label="Calculate new DAScale value by\n extrapolating fit slope and offset\n and using MaximumChangePercent - 2 as frequency distance.\n Uses the fit slopes and offsets from the\n passing sweep with the highest AP frequency searching from the back."];
n103 -> n106 [label=Yes];
n108 [label="Use the next DAScale value"];
n103 -> n108 [label=No];
Expand Down
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 2606254

Please sign in to comment.