diff --git a/QuickCSF/QuickCSF.py b/QuickCSF/QuickCSF.py index 58a23fc..8141957 100644 --- a/QuickCSF/QuickCSF.py +++ b/QuickCSF/QuickCSF.py @@ -59,14 +59,62 @@ def makeFrequencySpace(min=.2, max=36, count=20): return frequencySpace -def csf_unmapped(parameters, frequency): - '''The truncated log-parabola model for human contrast sensitivity +def makeLogLinearSpace(min, max, count): + '''Creates values at log-linear equally spaced intervals''' + + expRange = numpy.log10(max)-numpy.log10(min) + expMin = numpy.log10(min) + + logLinearSpace = numpy.array([0.] * count) + for i in range(count): + logLinearSpace[i] = ( + (i/(count-1)) * expRange + expMin + ) + return logLinearSpace + +def makePeakSensitivitySpace(min=2, max=1000, count=28): + '''Creates peak sensitivity values for searching at log-linear equally spaced intervals''' + + logger.debug('Making peak sensitivity space: ' + str(locals())) + + peakSensitivitySpace = makeLogLinearSpace(min, max, count) + + return peakSensitivitySpace + +def makePeakFrequencySpace(min=.2, max=20, count=21): + '''Creates peak frequency values for searching at log-linear equally spaced intervals''' + + logger.debug('Making peak frequency space: ' + str(locals())) + + peakFrequencySpace = makeLogLinearSpace(min, max, count) + + return peakFrequencySpace + +def makeBandwidthSpace(min=1, max=10, count=21): + '''Creates bandwidth values for searching at log-linear equally spaced intervals''' + + logger.debug('Making bandwidth space: ' + str(locals())) + + peakBandwidthSpace = makeLogLinearSpace(min, max, count) + return peakBandwidthSpace + +def makeLogDeltaSpace(min=.02, max=2, count=21): + '''Creates log delta values for searching at log-linear equally spaced intervals''' + + logger.debug('Making log delta space: ' + str(locals())) + + peakLogDeltaSpace = makeLogLinearSpace(min, max, count) + + return peakLogDeltaSpace + +def csf_unmapped(parameterIndices, parameterSpace, frequency, continuous=False): + '''The truncated log-parabola model for human contrast sensitivity Expects UNMAPPED parameters Param order = peak sensitivity, peak frequency, bandwidth, log delta ''' # Get everything into log-units - [peakSensitivity, peakFrequency, logBandwidth, delta] = mapCSFParams(parameters) + [peakSensitivity, peakFrequency, logBandwidth, delta] = mapCSFParams(parameterIndices, parameterSpace, False, continuous) return csf(peakSensitivity, peakFrequency, logBandwidth, delta, frequency) @@ -123,7 +171,14 @@ def myCSF(frequency): return area -def mapCSFParams(params, exponify=False): +def mapParams(paramIndices, paramSpace): + expMin = numpy.min(paramSpace) + expMax = numpy.max(paramSpace) + expRange = expMax-expMin + count = len(paramSpace) + return paramIndices/(count-1) * expRange + expMin + +def mapCSFParams(paramIndices, paramSpace, exponify=False, continuous=False): ''' Maps parameter indices to log values @@ -133,15 +188,22 @@ def mapCSFParams(params, exponify=False): Bandwidth: octaves Delta: 1/contrast (Difference between Peak Sensitivity and the truncation) ''' - peakSensitivity = 0.1*params[:,0] + 0.3 - peakFrequency = -0.7 + 0.1*params[:,1] - bandwidth = 0.05 * params[:,2] - logDelta = -1.7 + 0.1 * params[:,3] + # log-linear normalize + if continuous: + peakSensitivity = mapParams(paramIndices[:,0], paramSpace[0]) + peakFrequency = mapParams(paramIndices[:,1], paramSpace[1]) + bandwidth = mapParams(paramIndices[:,2], paramSpace[2]) + logDelta = mapParams(paramIndices[:,3], paramSpace[3]) + else: + peakSensitivity = paramSpace[0][paramIndices[:,0]] + peakFrequency = paramSpace[1][paramIndices[:,1]] + bandwidth = paramSpace[2][paramIndices[:,2]] + logDelta = paramSpace[3][paramIndices[:,3]] + delta = numpy.power(10, logDelta) if exponify: deltaDiff = numpy.power(10, peakSensitivity-delta) - peakSensitivity = numpy.power(10, peakSensitivity) peakFrequency = numpy.power(10, peakFrequency) bandwidth = numpy.power(10, bandwidth) @@ -153,26 +215,28 @@ def entropy(p): return numpy.multiply(-p, numpy.log(p)) - numpy.multiply(1-p, numpy.log(1-p)) class QuickCSFEstimator(): - def __init__(self, stimulusSpace=None): + def __init__(self, parameterSpace=None, stimulusSpace=None): '''Create a new QuickCSF estimator with the specified input/output spaces Args: + parameterSpace: 4,x numpy array of attributes to be used for four parameters searching + numpy.array([peakSensitivity, peakFrequency, bandwidth, logdelta]) stimulusSpace: 2,x numpy array of attributes to be used for stimulus generation numpy.array([contrasts, frequencies]) ''' if stimulusSpace is None: stimulusSpace = numpy.array([ - makeContrastSpace(.0001, .05), - makeFrequencySpace() + makeContrastSpace(.001, 1, 24), + makeFrequencySpace(.2, 36, 20) ]) - - parameterSpace = numpy.array([ - numpy.arange(0, 28), # Peak sensitivity - numpy.arange(0, 21), # Peak frequency - numpy.arange(0, 21), # Log bandwidth - numpy.arange(0, 21) # Low frequency truncation (log delta) - ]) - + if parameterSpace is None: + parameterSpace = numpy.array([ + makePeakSensitivitySpace(2, 1000, 28), # Peak sensitivity + makePeakFrequencySpace(0.2, 20, 21), # Peak frequency + makeBandwidthSpace(1, 10, 21), # Log bandwidth + makeLogDeltaSpace(0.02, 2, 21) # Low frequency truncation (log delta) + ]) + logger.info('Initializing QuickCSFEStimator') logger.debug('Initializing QuickCSFEstimator stimSpace='+str(stimulusSpace).replace('\n','')+', paramSpace='+str(parameterSpace).replace('\n','')) @@ -203,7 +267,7 @@ def next(self): # more probable stim params have higher weight of being sampled randomSampleCount = 100 - paramIndicies = numpy.random.choice( + paramIndices = numpy.random.choice( numpy.arange(self.paramComboCount), randomSampleCount, p=self.probabilities[:,0] @@ -211,8 +275,8 @@ def next(self): # calculate probabilities for all stimuli with all samples of parameters # @TODO: parallelize this - stimIndicies = numpy.arange(self.stimComboCount).reshape(-1,1) - p = self._pmeas(paramIndicies, stimIndicies) + stimIndices = numpy.arange(self.stimComboCount).reshape(-1,1) + p = self._pmeas(paramIndices, stimIndices) # Determine amount of information to be gained pbar = sum(p)/randomSampleCount @@ -269,7 +333,7 @@ def _pmeas(self, parameterIndex, stimulusIndex=None): stimulusIndices = self.inflateStimulusIndex(stimulusIndex) frequencies = self.stimulusSpace[1][stimulusIndices[:,1]].reshape(1,-1) - csfValues = csf_unmapped(parameters, frequencies) + csfValues = csf_unmapped(parameters, self.parameterSpace, frequencies) # Make vector of sensitivities contrast = self.stimulusSpace[0][stimulusIndices[:,0]] @@ -334,7 +398,7 @@ def getResults(self, leaveAsIndices=False): '''Calculate an estimate of all 4 parameters based on their probabilities Args: - leaveAsIndicies: if False, will output real-world, linear-scale values + leaveAsIndices: if False, will output real-world, linear-scale values if True, will output indices, which can be converted with `mapCSFParams()` ''' @@ -349,11 +413,11 @@ def getResults(self, leaveAsIndices=False): results = estimatedParamMeans.reshape(1,len(self.parameterRanges)) - if not leaveAsIndices: - results = mapCSFParams(results, True).T - + if leaveAsIndices: + return results + + results = mapCSFParams(results, self.parameterSpace, True, True).T results = results.reshape(4).tolist() - return { 'peakSensitivity': results[0], 'peakFrequency': results[1], diff --git a/QuickCSF/StimulusGenerators.py b/QuickCSF/StimulusGenerators.py index 498dd78..6114a6b 100644 --- a/QuickCSF/StimulusGenerators.py +++ b/QuickCSF/StimulusGenerators.py @@ -30,12 +30,22 @@ def __init__(self, size=100, orientation=None, minContrast=.01, maxContrast=1.0, contrastResolution=24, minFrequency=0.2, maxFrequency=36.0, frequencyResolution=20, + minPeakSensitivity=2, maxPeakSensitivity=1000, peakSensitivityResolution=28, + minPeakFrequency=0.2, maxPeakFrequency=20, peakFrequencyResolution=21, + minBandwidth=1, maxBandwidth=10, bandwidthResolution=21, + minLogDelta=0.02, maxLogDelta=2, logDeltaResolution=21, degreesToPixels=None ): super().__init__( stimulusSpace = numpy.array([ QuickCSF.makeContrastSpace(minContrast, maxContrast, contrastResolution), QuickCSF.makeFrequencySpace(minFrequency, maxFrequency, frequencyResolution) + ]), + parameterSpace = numpy.array([ + QuickCSF.makePeakSensitivitySpace(minPeakSensitivity, maxPeakSensitivity, peakSensitivityResolution), + QuickCSF.makePeakFrequencySpace(minPeakFrequency, maxPeakFrequency, peakFrequencyResolution), + QuickCSF.makeBandwidthSpace(minBandwidth, maxBandwidth, bandwidthResolution), + QuickCSF.makeLogDeltaSpace(minLogDelta, maxLogDelta, logDeltaResolution) ]) ) diff --git a/QuickCSF/app.py b/QuickCSF/app.py index 1933f64..01ccf5e 100644 --- a/QuickCSF/app.py +++ b/QuickCSF/app.py @@ -69,7 +69,7 @@ def onStateTransition(state, data): degreesToPixels = functools.partial(screens.degreesToPixels, distance_mm=settings['distance_mm']) - stimGenerator = StimulusGenerators.QuickCSFGenerator(degreesToPixels=degreesToPixels, **settings['Stimuli']) + stimGenerator = StimulusGenerators.QuickCSFGenerator(degreesToPixels=degreesToPixels, **settings['Stimuli'], **settings['Parameters']) controller = CSFController.Controller_2AFC(stimGenerator, **settings['Controller']) mainWindow.participantReady.connect(controller.onParticipantReady) @@ -114,7 +114,7 @@ def getSettings(): controllerSettings.add_argument('--waitForReady', default=False, action='store_true', help='Wait for the participant to indicate they are ready for the next trial') stimulusSettings = parser.add_argument_group('Stimuli') - stimulusSettings.add_argument('-minc', '--minContrast', type=float, default=.01, help='The lowest contrast value to measure (0.0-1.0)') + stimulusSettings.add_argument('-minc', '--minContrast', type=float, default=.001, help='The lowest contrast value to measure (0.0-1.0)') stimulusSettings.add_argument('-maxc', '--maxContrast', type=float, default=1.0, help='The highest contrast value to measure (0.0-1.0)') stimulusSettings.add_argument('-cr', '--contrastResolution', type=int, default=24, help='The number of contrast steps') @@ -125,6 +125,23 @@ def getSettings(): stimulusSettings.add_argument('--size', type=int, default=3, help='Gabor patch size in (degrees)') stimulusSettings.add_argument('--orientation', type=float, help='Orientation of gabor patch (degrees). If unspecified, each trial will be random') + parameterSettings = parser.add_argument_group('Parameters') + parameterSettings.add_argument('-minps', '--minPeakSensitivity', type=float, default=2.0, help='The lower bound of peak sensitivity value (>1.0)') + parameterSettings.add_argument('-maxps', '--maxPeakSensitivity', type=float, default=1000.0, help='The upper bound of peak sensitivity value') + parameterSettings.add_argument('-psr', '--peakSensitivityResolution', type=int, default=28, help='The number of peak sensitivity steps') + + parameterSettings.add_argument('-minpf', '--minPeakFrequency', type=float, default=.2, help='The lower bound of peak frequency value (>0)') + parameterSettings.add_argument('-maxpf', '--maxPeakFrequency', type=float, default=20.0, help='The upper bound of peak frequency value') + parameterSettings.add_argument('-pfr', '--peakFrequencyResolution', type=int, default=21, help='The number of peak frequency steps') + + parameterSettings.add_argument('-minb', '--minBandwidth', type=float, default=1.0, help='The lower bound of bandwidth value') + parameterSettings.add_argument('-maxb', '--maxBandwidth', type=float, default=10.0, help='The upper bound of bandwidth value') + parameterSettings.add_argument('-br', '--bandwidthResolution', type=int, default=21, help='The number of bandwidth steps') + + parameterSettings.add_argument('-mind', '--minLogDelta', type=float, default=.02, help='The lower bound of logdelta value') + parameterSettings.add_argument('-maxd', '--maxLogDelta', type=float, default=2.0, help='The upper bound of logdelta value') + parameterSettings.add_argument('-dr', '--logDeltaResolution', type=int, default=21, help='The number of logdelta steps') + settings = argparseqt.groupingTools.parseIntoGroups(parser) if None in [settings['sessionID'], settings['distance_mm']]: settings = ui.getSettings(parser, settings, ['sessionID', 'distance_mm']) @@ -136,6 +153,12 @@ def main(): settings = getSettings() if not settings is None: + # parameter validation: peak sensitivty space and peak frequency space should be subspaces of stimulus spaces (1/contrast, frequency) + if not (settings['Parameters']['minPeakSensitivity'] >= 1/settings['Stimuli']['maxContrast'] and settings['Parameters']['maxPeakSensitivity'] <= 1/settings['Stimuli']['minContrast']): + raise ValueError("Please increase the range of contrast space or decrease the range of peak sensitivity space.") + if not (settings['Parameters']['minPeakFrequency'] >= settings['Stimuli']['minFrequency'] and settings['Parameters']['maxPeakFrequency'] <= settings['Stimuli']['maxFrequency']): + raise ValueError("Please increase the range of frequency space or decrease the range of peak frequency space.") + logPath = pathlib.Path(settings['outputFile']).parent log.startLog(settings['sessionID'], logPath) run(settings) diff --git a/QuickCSF/simulate.py b/QuickCSF/simulate.py index b8a66c6..6da2f9a 100644 --- a/QuickCSF/simulate.py +++ b/QuickCSF/simulate.py @@ -31,7 +31,7 @@ def plot(qCSFEstimator, graph=None, unmappedTrueParams=None, showNumbers=True): frequencyDomain = QuickCSF.makeFrequencySpace(.005, 80, 50).reshape(-1,1) if unmappedTrueParams is not None: - truthData = QuickCSF.csf_unmapped(unmappedTrueParams.reshape(1, -1), frequencyDomain) + truthData = QuickCSF.csf_unmapped(unmappedTrueParams.reshape(1, -1), qCSFEstimator.parameterSpace, frequencyDomain) truthData = numpy.power(10, truthData) truthLine = graph.fill_between( frequencyDomain.reshape(-1), @@ -42,13 +42,7 @@ def plot(qCSFEstimator, graph=None, unmappedTrueParams=None, showNumbers=True): truthData = None estimatedParamMeans = qCSFEstimator.getResults(leaveAsIndices=True) - estimatedParamMeans = numpy.array([[ - estimatedParamMeans['peakSensitivity'], - estimatedParamMeans['peakFrequency'], - estimatedParamMeans['bandwidth'], - estimatedParamMeans['delta'], - ]]) - estimatedData = QuickCSF.csf_unmapped(estimatedParamMeans.reshape(1, -1), frequencyDomain) + estimatedData = QuickCSF.csf_unmapped(estimatedParamMeans.reshape(1, -1), qCSFEstimator.parameterSpace, frequencyDomain, True) estimatedData = numpy.power(10, estimatedData) estimatedLine = graph.fill_between( @@ -84,13 +78,13 @@ def plot(qCSFEstimator, graph=None, unmappedTrueParams=None, showNumbers=True): graph.grid() if showNumbers: - estimatedParamMeans = QuickCSF.mapCSFParams(estimatedParamMeans, exponify=True) + estimatedParamMeans = QuickCSF.mapCSFParams(estimatedParamMeans, qCSFEstimator.parameterSpace, exponify=True, continuous=True) estimatedParamMeans = estimatedParamMeans.reshape(1,-1).tolist()[0] paramEstimates = '%03.2f, %.4f, %.4f, %.4f' % tuple(estimatedParamMeans) estimatedLine.set_label(f'Estim: {paramEstimates}') if truthData is not None: - trueParams = QuickCSF.mapCSFParams(unmappedTrueParams, True).T.tolist()[0] + trueParams = QuickCSF.mapCSFParams(unmappedTrueParams, qCSFEstimator.parameterSpace, True).T.tolist()[0] trueParams = '%03.2f, %.4f, %.4f, %.4f' % tuple(trueParams) truthLine.set_label(f'Truth: {trueParams}') @@ -105,12 +99,16 @@ def runSimulation( imagePath=None, usePerfectResponses=False, stimuli={ - 'minContrast':0.01, 'maxContrast':1, 'contrastResolution':24, + 'minContrast':0.001, 'maxContrast':1, 'contrastResolution':24, 'minFrequency':.2, 'maxFrequency':36, 'frequencyResolution':20, }, parameters={ 'truePeakSensitivity':18, 'truePeakFrequency':11, 'trueBandwidth':12, 'trueDelta':11, + 'minPeakSensitivity':2, 'maxPeakSensitivity':1000, 'peakSensitivityResolution':28, + 'minPeakFrequency':0.2, 'maxPeakFrequency':20, 'peakFrequencyResolution': 21, + 'minBandwidth':1, 'maxBandwidth':10, 'bandwidthResolution': 21, + 'minLogDelta':0.02, 'maxLogDelta':2, 'logDeltaResolution': 21, }, ): logger.info('Starting simulation') @@ -124,6 +122,12 @@ def runSimulation( QuickCSF.makeContrastSpace(stimuli['minContrast'], stimuli['maxContrast'], stimuli['contrastResolution']), QuickCSF.makeFrequencySpace(stimuli['minFrequency'], stimuli['maxFrequency'], stimuli['frequencyResolution']) ]) + parameterSpace = numpy.array([ + QuickCSF.makePeakSensitivitySpace(parameters['minPeakSensitivity'], parameters['maxPeakSensitivity'], parameters['peakSensitivityResolution']), # Peak sensitivity + QuickCSF.makePeakFrequencySpace(parameters['minPeakFrequency'], parameters['maxPeakFrequency'], parameters['peakFrequencyResolution']), # Peak frequency + QuickCSF.makeBandwidthSpace(parameters['minBandwidth'], parameters['maxBandwidth'], parameters['bandwidthResolution']), # Log bandwidth + QuickCSF.makeLogDeltaSpace(parameters['minLogDelta'], parameters['maxLogDelta'], parameters['logDeltaResolution']) # Low frequency truncation (log delta) + ]) unmappedTrueParams = numpy.array([[ parameters['truePeakSensitivity'], @@ -131,7 +135,7 @@ def runSimulation( parameters['trueBandwidth'], parameters['trueDelta'], ]]) - qcsf = QuickCSF.QuickCSFEstimator(stimulusSpace) + qcsf = QuickCSF.QuickCSFEstimator(parameterSpace, stimulusSpace) graph = plot(qcsf, unmappedTrueParams=unmappedTrueParams) @@ -145,7 +149,7 @@ def runSimulation( if usePerfectResponses: logger.debug('Simulating perfect response') frequency = newStimValues[:,1] - trueSens = numpy.power(10, QuickCSF.csf_unmapped(unmappedTrueParams, numpy.array([frequency]))) + trueSens = numpy.power(10, QuickCSF.csf_unmapped(unmappedTrueParams, qcsf.parameterSpace, numpy.array([frequency]))) testContrast = newStimValues[:,0] testSens = 1 / testContrast @@ -175,7 +179,7 @@ def runSimulation( paramEstimates = qcsf.getResults() logger.info('Results: ' + str(paramEstimates)) - trueParams = QuickCSF.mapCSFParams(unmappedTrueParams, True).T + trueParams = QuickCSF.mapCSFParams(unmappedTrueParams, qcsf.parameterSpace, True).T print('******* Results *******') print(f'\tEstimates = {paramEstimates}') print(f'\tActuals = {trueParams}') @@ -215,7 +219,7 @@ def entropyPlot(qcsf): parser.add_argument('-perfect', '--usePerfectResponses', default=False, action='store_true', help='Whether to simulate perfect responses, rather than probablistic ones') stimuliSettings = parser.add_argument_group('stimuli') - stimuliSettings.add_argument('-minc', '--minContrast', type=float, default=.01, help='The lowest contrast value to measure (0.0-1.0)') + stimuliSettings.add_argument('-minc', '--minContrast', type=float, default=.001, help='The lowest contrast value to measure (0.0-1.0)') stimuliSettings.add_argument('-maxc', '--maxContrast', type=float, default=1.0, help='The highest contrast value to measure (0.0-1.0)') stimuliSettings.add_argument('-cr', '--contrastResolution', type=int, default=24, help='The number of contrast steps') @@ -229,6 +233,22 @@ def entropyPlot(qcsf): parameterSettings.add_argument('-b', '--trueBandwidth', type=int, default=12, help='True bandwidth (index)') parameterSettings.add_argument('-d', '--trueDelta', type=int, default=11, help='True delta truncation (index)') + parameterSettings.add_argument('-minps', '--minPeakSensitivity', type=float, default=2.0, help='The lower bound of peak sensitivity value (>1.0)') + parameterSettings.add_argument('-maxps', '--maxPeakSensitivity', type=float, default=1000.0, help='The upper bound of peak sensitivity value') + parameterSettings.add_argument('-psr', '--peakSensitivityResolution', type=int, default=28, help='The number of peak sensitivity steps') + + parameterSettings.add_argument('-minpf', '--minPeakFrequency', type=float, default=.2, help='The lower bound of peak frequency value (>0)') + parameterSettings.add_argument('-maxpf', '--maxPeakFrequency', type=float, default=20.0, help='The upper bound of peak frequency value') + parameterSettings.add_argument('-pfr', '--peakFrequencyResolution', type=int, default=21, help='The number of peak frequency steps') + + parameterSettings.add_argument('-minb', '--minBandwidth', type=float, default=1.0, help='The lower bound of bandwidth value') + parameterSettings.add_argument('-maxb', '--maxBandwidth', type=float, default=10.0, help='The upper bound of bandwidth value') + parameterSettings.add_argument('-br', '--bandwidthResolution', type=int, default=21, help='The number of bandwidth steps') + + parameterSettings.add_argument('-mind', '--minLogDelta', type=float, default=.02, help='The lower bound of logdelta value') + parameterSettings.add_argument('-maxd', '--maxLogDelta', type=float, default=2.0, help='The upper bound of logdelta value') + parameterSettings.add_argument('-dr', '--logDeltaResolution', type=int, default=21, help='The number of logdelta steps') + settings = argparseqt.groupingTools.parseIntoGroups(parser) if settings['trials'] is None: from qtpy import QtWidgets @@ -237,4 +257,10 @@ def entropyPlot(qcsf): settings = ui.getSettings(parser, settings, ['trials']) if settings is not None: + # parameter validation: peak sensitivty space and peak frequency space should be subspaces of stimulus spaces (1/contrast, frequency) + if not (settings['parameters']['minPeakSensitivity'] >= 1/settings['stimuli']['maxContrast'] and settings['parameters']['maxPeakSensitivity'] <= 1/settings['stimuli']['minContrast']): + raise ValueError("Please increase the range of contrast space or decrease the range of peak sensitivity space.") + if not (settings['parameters']['minPeakFrequency'] >= settings['stimuli']['minFrequency'] and settings['parameters']['maxPeakFrequency'] <= settings['stimuli']['maxFrequency']): + raise ValueError("Please increase the range of frequency space or decrease the range of peak frequency space.") + runSimulation(**settings)