This file is indexed.

/usr/share/psychtoolbox-3/PsychCal/CalibrateFitGamma.m is in psychtoolbox-3-common 3.0.11.20140816.dfsg1-1.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
function cal = CalibrateFitGamma(cal,nInputLevels)
% cal = CalibrateFitGamma(cal,[nInputLevels])
%
% Fit the gamma function to the calibration measurements.  Options for field
% cal.describe.gamma.fitType are:
%    simplePower
%    crtLinear
%    crtPolyLinear
%    crtGamma
%    crtSumPow
%    betacdf
%    sigmoid
%    weibull
%
% Underlying fit routine is FitGamma for functional forms originally supported,
% and these rely on the optimization toolbox.
%
% Newer functions (e.g, crtSumPow, betacdf) use the curvefit toolbox and that's just
% done locally in this routine.  Much less cumbersome.
%
% NOTE (5/27/10, dhb): crtSumPow does not currently appear to normalize the
% measurements to unity, while the older methods do (in FitDeviceGamma).
% This may be a bug, but since we're not currently using crtSumPow I'm not
% going to look into it in detail right now.
%
% See also PsychGamma.
%
% 3/26/02  dhb  Pulled out of CalibrateMonDrvr.
% 11/14/06 dhb  Define nInputLevels and pass to underlying fit routine.
% 07/22/07 dhb  Add simplePower fitType.
% 08/02/07 dhb  Optional pass of nInputLevels.
%          dhb  Don't allow a long string of zeros at the start.
%          dhb  Reduce redundant code for higher order terms by pulling out of switch
% 08/03/07 dhb  Debug.  Add call to MakeMonotonic for first three components.
% 11/19/09 dhb  Added crtSumPow option, coded to [0-1] world and using curve fit toolbox.
% 3/07/10  dhb  Cosmetic to make m-lint happier, including some "|" -> "||"
% 3/07/10  dhb  Added crtLinear option.
%          dhb  contrasthThresh and fitBreakThresh values only set if not already in struct.
%          dhb  Call MakeGammaMonotonic rather than MakeMonotonic where appropriate.
%          dhb  Use linear interpolation for higher order linear model weights, rather than
%               a polynomial.  I now think that ringing is worse than not smoothing enough.
% 3/08/10  dhb  Update list of options in comment above.
% 5/26/10  dhb  Allow gamma input values to be either a single column or a matrix with same number of columns as devices.
% 6/5/10   dhb  Extend fix above to higher order terms in the gamma fit.
%          dhb  Fix or supress MATLAB lint warnings.
%          dhb  Add betacdf fit option, which seems to provide a flexible sigmoidally shaped fit.
% 6/8/10   dhb, ar Make sure to set cal.gammaInput in options that use curvefit toolbox method.
%               Add a call to MakeGammaMonotonic around input values for higher order linmod fit.
% 6/1010   dhb  Fix higher order fit in case where there are multiple gamma input columns.  Blew this the other day.
% 6/11/10  dhb  Allow passing of weighting parameter as part of cal.describe.gamma structure.  Change functional form of betacdf
%               to include wrapped power functions.
% 4/12/11  dhb  For simplePower option, return vector of exponents in cal.describe.exponents.

% Set nInputLevels
if (nargin < 2 || isempty(nInputLevels))
    nInputLevels = 1024;
end

% Fit gamma functions.
switch(cal.describe.gamma.fitType)
    
    case 'simplePower',
        mGammaMassaged = cal.rawdata.rawGammaTable(:,1:cal.nDevices);
        for i = 1:cal.nDevices
            mGammaMassaged(:,i) = MakeGammaMonotonic(HalfRect(mGammaMassaged(:,i)));
        end
        fitType = 1;
        [mGammaFit1a,cal.gammaInput,nil,theExponents] = FitDeviceGamma(...
            mGammaMassaged,cal.rawdata.rawGammaInput,fitType,nInputLevels);
        mGammaFit1 = mGammaFit1a;
        cal.describe.gamma.exponents = theExponents;
        
    case 'crtLinear'
        % Set to zero the raw data we believe to be below reliable measurement
        % threshold, and then fit the rest by linear interpolation.  Force answer
        % to be monotonic.
        if (~isfield(cal.describe.gamma,'contrastThresh'))
            cal.describe.gamma.contrastThresh = 0.001;
        end
        mGammaMassaged = cal.rawdata.rawGammaTable(:,1:cal.nDevices);
        massIndex = find(mGammaMassaged < cal.describe.gamma.contrastThresh);
        mGammaMassaged(massIndex) = zeros(length(massIndex),1);
        for i = 1:cal.nDevices
            mGammaMassaged(:,i) = MakeGammaMonotonic(HalfRect(mGammaMassaged(:,i)));
        end
        
        fitType = 6;
        [mGammaFit1,cal.gammaInput] = FitDeviceGamma(...
            mGammaMassaged,cal.rawdata.rawGammaInput,fitType,nInputLevels);
        
    case 'crtPolyLinear',
        % For fitting, we set to zero the raw data we
        % believe to be below reliable measurement threshold (contrastThresh).
        % Currently we are fitting both with polynomial and a linear interpolation,
        % using the latter for low measurement values.  The fit break point is
        % given by fitBreakThresh.   This technique was developed
        % through bitter experience and is not theoretically driven.
        if (~isfield(cal.describe.gamma,'contrastThresh'))
            cal.describe.gamma.contrastThresh = 0.001;
        end
        if (~isfield(cal.describe.gamma,'fitBreakThresh'))
            cal.describe.gamma.fitBreakThresh = 0.02;
        end
        mGammaMassaged = cal.rawdata.rawGammaTable(:,1:cal.nDevices);
        massIndex = find(mGammaMassaged < cal.describe.gamma.contrastThresh);
        mGammaMassaged(massIndex) = zeros(length(massIndex),1);
        for i = 1:cal.nDevices
            mGammaMassaged(:,i) = MakeGammaMonotonic(HalfRect(mGammaMassaged(:,i)));
        end
        fitType = 7;
        [mGammaFit1a,cal.gammaInput] = FitDeviceGamma(...
            mGammaMassaged,cal.rawdata.rawGammaInput,fitType,nInputLevels);
        fitType = 6;
        [mGammaFit1b,cal.gammaInput] = FitDeviceGamma(...
            mGammaMassaged,cal.rawdata.rawGammaInput,fitType,nInputLevels);
        mGammaFit1 = mGammaFit1a;
        for i = 1:cal.nDevices
            indexLin = find(mGammaMassaged(:,i) < cal.describe.gamma.fitBreakThresh);
            if (~isempty(indexLin))
                breakIndex = max(indexLin);
                breakInput = cal.rawdata.rawGammaInput(breakIndex);
                inputIndex = find(cal.gammaInput <= breakInput);
                if (~isempty(inputIndex))
                    mGammaFit1(inputIndex,i) = mGammaFit1b(inputIndex,i);
                end
            end
        end
        
    case 'crtGamma',
        mGammaMassaged = cal.rawdata.rawGammaTable(:,1:cal.nDevices);
        for i = 1:cal.nDevices
            mGammaMassaged(:,i) = MakeGammaMonotonic(HalfRect(mGammaMassaged(:,i)));
        end
        fitType = 2;
        [mGammaFit1a,cal.gammaInput] = FitDeviceGamma(...
            mGammaMassaged,cal.rawdata.rawGammaInput,fitType,nInputLevels);
        mGammaFit1 = mGammaFit1a;
        
    case 'crtSumPow',
        if (~exist('fit','file'))
            error('Fitting with the sum of exponentials requires the curve fitting toolbox\n');
        end
        if (max(cal.rawdata.rawGammaInput(:)) > 1)
            error('crtSumPower option assumes [0-1] specification of input\n');
        end
        mGammaMassaged = cal.rawdata.rawGammaTable(:,1:cal.nDevices);
        for i = 1:cal.nDevices
            mGammaMassaged(:,i) = MakeGammaMonotonic(HalfRect(mGammaMassaged(:,i)));
        end
        
        fitEqStr = 'a*x^b + (1-a)*x^c';
        a = 1;
        b = 2;
        c = 0;
        startPoint = [a b c];
        lowerBounds = [0 0.1 0.01];
        upperBounds = [1 10 10];
        
        % Fit and predictions
        fOptions = fitoptions('Method','NonlinearLeastSquares','Robust','on');
        fOptions1 = fitoptions(fOptions,'StartPoint',startPoint,'Lower',lowerBounds,'Upper',upperBounds);
        for i = 1:cal.nDevices
            if (size(cal.rawdata.rawGammaInput,2) == 1)
                fitstruct = fit(cal.rawdata.rawGammaInput,mGammaMassaged(:,i),fitEqStr,fOptions1);
            else
                fitstruct = fit(cal.rawdata.rawGammaInput(:,i),mGammaMassaged(:,i),fitEqStr,fOptions1);
            end
            mGammaFit1a(:,i) = feval(fitstruct,linspace(0,1,nInputLevels)); %#ok<*AGROW>
        end
        mGammaFit1 = mGammaFit1a;
        cal.gammaInput = linspace(0,1,nInputLevels)';
  
    case 'betacdf',
        if (~exist('fit','file'))
            error('Fitting with the betacdf requires the curve fitting toolbox\n');
        end
        if (~exist('betacdf','file'))
            error('Fitting with the betacdf requires the stats toolbox\n');
        end
        if (max(cal.rawdata.rawGammaInput(:)) > 1)
            error('betacdf option assumes [0-1] specification of input\n');
        end
        mGammaMassaged = cal.rawdata.rawGammaTable(:,1:cal.nDevices);
        for i = 1:cal.nDevices
            mGammaMassaged(:,i) = MakeGammaMonotonic(HalfRect(mGammaMassaged(:,i)));
        end
        
        fitEqStr = 'betacdf(betacdf(x.^f,a,b),c,d).^e';
        a = 1;
        b = 1;
        c = 1;
        d = 1;
        e = 1;
        f = 1;
        startPoint = [a b c d e f];
        lowerBounds = [1e-3 1e-3 1e-3 1e-3 1e-3 1e-3];
        upperBounds = [1e3 1e3 1e3 1e3 1e3 1e3];
        
        % Fit and predictions
        fOptions = fitoptions('Method','NonlinearLeastSquares','Robust','on','Display','off');
        fOptions1 = fitoptions(fOptions,'StartPoint',startPoint,'Lower',lowerBounds,'Upper',upperBounds,'MaxFunEvals',2000);
        for i = 1:cal.nDevices
            if (isfield(cal.describe.gamma,'useweight') && cal.describe.gamma.useweight >= 0)
                fOptionsUse = fitoptions(fOptions1,'Weights',1./(mGammaMassaged(:,i)+cal.describe.gamma.useweight));
            else
                fOptionsUse = fOptions1;
            end
            if (size(cal.rawdata.rawGammaInput,2) == 1)
                fitstruct = fit(cal.rawdata.rawGammaInput,mGammaMassaged(:,i),fitEqStr,fOptionsUse);
            else
                fitstruct = fit(cal.rawdata.rawGammaInput(:,i),mGammaMassaged(:,i),fitEqStr,fOptionsUse);
            end
            mGammaFit1a(:,i) = feval(fitstruct,linspace(0,1,nInputLevels)); %#ok<*AGROW>
        end
        mGammaFit1 = mGammaFit1a;
        cal.gammaInput = linspace(0,1,nInputLevels)';
		
    case 'sigmoid',
        mGammaMassaged = cal.rawdata.rawGammaTable(:,1:cal.nDevices);
        for i = 1:cal.nDevices
            mGammaMassaged(:,i) = MakeGammaMonotonic(HalfRect(mGammaMassaged(:,i)));
        end
        fitType = 3;
        [mGammaFit1a,cal.gammaInput] = FitDeviceGamma(...
            mGammaMassaged,cal.rawdata.rawGammaInput,fitType,nInputLevels);
        mGammaFit1 = mGammaFit1a;
        
    case 'weibull',
        mGammaMassaged = cal.rawdata.rawGammaTable(:,1:cal.nDevices);
        for i = 1:cal.nDevices
            mGammaMassaged(:,i) = MakeGammaMonotonic(HalfRect(mGammaMassaged(:,i)));
        end
        fitType = 4;
        [mGammaFit1a,cal.gammaInput] = FitDeviceGamma(...
            mGammaMassaged,cal.rawdata.rawGammaInput,fitType,nInputLevels);
        mGammaFit1 = mGammaFit1a;
        
    otherwise
        error('Unsupported gamma fit string passed');
        
end

% Fix contingous zeros at start problem
mGammaFit1 = FixZerosAtStart(mGammaFit1);
for j = 1:size(mGammaFit1,2)
    mGammaFit1(:,j) = MakeGammaMonotonic(mGammaFit1(:,j));
end

% Handle higher order terms, which are just fit with a polynomial
if (cal.nPrimaryBases > 1)
    m = size(mGammaFit1,1);
    mGammaFit2 = zeros(m,cal.nDevices*(cal.nPrimaryBases-1));
    
    % OLDFIT path does not contain option of handling data with independent input values
    % for measurements for each device primary.
    OLDFIT = 0;
    if (OLDFIT)
        for j = 1:cal.nDevices*(cal.nPrimaryBases-1)
            mGammaFit2(:,j) = ...
                FitGammaPolyR(cal.rawdata.rawGammaInput,cal.rawdata.rawGammaTable(:,cal.nDevices+j), ...
                cal.gammaInput);
        end
        
    % This is the code we're currently using.  It works for the case where different input levels are specified for
    % the measurments for each primary.
    else
        k = 1;
        for j = 1:cal.nDevices*(cal.nPrimaryBases-1)
            if (size(cal.rawdata.rawGammaInput,2) > 1)
                mGammaFit2(:,j) = interp1(MakeGammaMonotonic([0 ; cal.rawdata.rawGammaInput(:,k)]),[0 ; cal.rawdata.rawGammaTable(:,cal.nDevices+j)],cal.gammaInput,'linear');
            else
                mGammaFit2(:,j) = interp1(MakeGammaMonotonic([0 ; cal.rawdata.rawGammaInput]),[0 ; cal.rawdata.rawGammaTable(:,cal.nDevices+j)],cal.gammaInput,'linear');
            end
            k = k+1;
            if (k == cal.nDevices+1)
                k = 1;
            end
        end
    end
    
    mGammaFit = [mGammaFit1 , mGammaFit2];
else
    mGammaFit = mGammaFit1;
end

% Save information in form for calibration routines.
cal.gammaFormat = 0;
cal.gammaTable = mGammaFit;

return

% output = FixZerosAtStart(input)
%
% The OS/X routines need the fit gamma function to be monotonically
% increasing.  One way that sometimes fails is when a whole bunch of
% entries at the start are zero.  This routine fixes that up.
function output = FixZerosAtStart(input)

output = input;
for j = 1:size(input,2)
    for i = 1:size(input,1)
        if (input(i,j) > 0)
            break;
        end
    end
    if (i == size(input,1))
        error('Entire passed gamma function is zero');
    end
    output(1:i,j) = linspace(0,min([0.0001 input(i+1,j)/2]),i)';
end

return