function [fit_out,x,fitComment] = ...
FitGamma(values_in,measurements,values_out,fitType)
% [fit_out,x,fitComment] = ...
% FitGamma(values_in,measurements,values_out,[fitType])
%
% Fit a gamma function. This essentially a driver function.
% It has two main purposes.
%
% First it tries several different
% underlying parametric forms for the fit and chooses the best
% for a particular data set.
%
% Second, it does the bookkeeping for fitting each column of
% input measurements. (Each of the underlying fit functions expects
% only vector input.)
%
% To a large extent, the interface to the underlying fit functions
% (e.g. FitGammaPow, FitGammaSig, ...) is uniform. However, this routine
% does have to know a little bit about initial value dimension and choice.
% We have tried to localize this information in the initialization routines
% (e.g. InitialXPow, InitialXSig, ...) as much as possible, but some
% caution is advised.
%
% Optional argument fitType allows you to force the return of a particular
% parametric form. Currently:
% fitType == 1: Power function
% fitType == 2: Extended power function
% fitType == 3: Sigmoid
% fitType == 4: Weibull
% fitType == 5: Modified polynomial
% fitType == 6: Linear interpolation
% fitType == 7: Cubic spline
%
% All fit types are in a form such that the fit is forced through the
% origin for 0 input. This is because our convention is that gamma
% correction happens after subtraction of the ambient light.
%
% NOTE: FitGammaPow (and perhaps other subroutines) uses CONSTR, which is part of the
% Mathworks Optimization Toolbox.
%
% Also see FitGammaDemo.
% 10/3/93 dhb Removed polynomial fit from list tried with fitType == 0.
% Added Weibull function fit
% 3/15/94 dhb, jms Added linear interpolation.
% 7/18/94 dhb Added cubic spline interpolation.
% 8/7/00 dhb Fix bug. Spline was calling linear interpolation. Thanks to
% Chien-Chung Chen for notifying us of this bug.
% 11/14/06 dhb Modify how default type is set. Handle passed empty matrix.
% 3/07/10 dhb Cosmetic to make m-lint happier, including some "|" -> "||"
% 5/26/10 dhb Allow values_in to be either a single column or a matrix with same number of columns as values_out.
% 6/5/10 dhb Fix a lot of MATLAB lint warnings.
% dhb Fix error reporting to actually take mean across devices.
% dhb Rewrite how mean is taken for evaluation of best fit. I think this was done right.
% 11/07/10 dhb Print out fit exponents when gamma fit by a simple power function.
% Get sizes
[nDevices] = size(measurements,2);
[nOut] = size(values_out,1);
% If input comes as a single column, then replicate it to
% match number of devices
if (size(values_in,2) == 1)
values_in = repmat(values_in,1,nDevices);
end
% Set up number of fit types
nFitTypes = 5;
errv = zeros(nFitTypes,nDevices);
% Handle force fittting
if (nargin < 4 || isempty(fitType))
fitType = 0;
end
if ~ismember(fitType, 0:7)
error('Unsupported fitType requested. Only 0 to 7 are supported.');
end
if ismember(fitType, [0, 1, 2, 3, 4]) && IsOctave
v = version;
if str2num(v(1)) < 6
error('For use with Octave, you need at least Octave version 6.');
end
try
% Try loading the optim package with the optimization functions:
pkg load optim;
catch
error('For use with Octave, you must install the ''optim'' package from Octave Forge. See ''help pkg''.');
end
% Got optim package loaded. Does it support fmincon()?
if ~exist('fmincon')
error('For use with Octave, you need at least version 1.6.0 of the ''optim'' package from Octave Forge.');
end
end
% Fit with simple power function through origin
if (fitType == 0 || fitType == 1 || fitType == 2)
disp('Fitting with simple power function');
fit_out1 = zeros(nOut,nDevices);
[nParams] = size(InitialXPow,1);
x1 = zeros(nParams,nDevices);
for i = 1:nDevices
x0 = InitialXPow;
[fit_out1(:,i),x1(:,i),errv(1,i)] = ...
FitGammaPow(values_in(:,i),measurements(:,i),values_out,x0);
fprintf('Exponent for device %d is %g\n',i,x1(:,i));
end
fprintf('Simple power function fit, RMSE: %g\n',mean(errv(1,i)));
end
% Fit with extended power function. Use power function
% fit to drive parameters. InitialXExtP can take a two
% vector as input. This defines the parameters of a good fitting
% simple power function.
if (fitType == 0 || fitType == 2)
disp('Fitting with extended power function');
fit_out2 = zeros(nOut,nDevices);
[nParams] = size(InitialXExtP,1);
x2 = zeros(nParams,nDevices);
for i = 1:nDevices
x0 = InitialXExtP(x1(:,i));
[fit_out2(:,i),x2(:,i),errv(2,i)] = ...
FitGammaExtP(values_in(:,i),measurements(:,i),values_out,x0);
end
fprintf('Extended power function fit, RMSE: %g\n',mean(errv(2,i)));
end
% Fit with a sigmoidal shape. This works well for
% the dimmer packs controlling lights. InitialXSig can take
% a two vector as input. This defines roughly the input for
% half maximum and the maximum output value.
if (fitType == 0 || fitType == 3)
disp('Fitting with sigmoidal function');
fit_out3 = zeros(nOut,nDevices);
[nParams] = size(InitialXSig,1);
x3 = zeros(nParams,nDevices);
for i = 1:nDevices
maxVals = max(values_in(:,i));
x0 = InitialXSig(maxVals'/2);
[fit_out3(:,i),x3(:,i),errv(3,i)] = ...
FitGammaSig(values_in(:,i),measurements(:,i),values_out,x0);
end
fprintf('Sigmoidal fit, RMSE: %g\n',mean(errv(3,i)));
end
% Fit with Weibull
if (fitType == 0 || fitType == 4)
disp('Fitting with Weibull function');
fit_out4 = zeros(nOut,nDevices);
[nParams] = size(InitialXWeib(values_in(:,1),measurements(:,1)),1);
x4 = zeros(nParams,nDevices);
for i = 1:nDevices
x0 = InitialXWeib(values_in(:,i),measurements(:,i));
[fit_out4(:,i),x4(:,i),errv(4,i)] = ...
FitGammaWeib(values_in(:,i),measurements(:,i),values_out,x0);
end
fprintf('Weibull function fit, RMSE: %g\n',mean(errv(4,i)));
end
% Fit with polynomial. InitalXPoly is used mostly for consistency
% with other calling forms, since FitGammaPoly computes an analytic
% fit to start the search. But it serves to implicitly defines the
% order of the polynomial.
if (fitType == 0 || fitType == 5)
disp('Fitting with polynomial');
fit_out5 = zeros(nOut,nDevices);
[order5] = size(InitialXPoly,1);
x5 = zeros(order5,nDevices);
for i = 1:nDevices
[fit_out5(:,i),x5(:,i),errv(5,i)] = ...
FitGammaPoly(values_in(:,i),measurements(:,i),values_out);
end
fprintf('Polynomial fit, order %g, RMSE: %g\n',order5,mean(errv(5,i)));
end
% Linear interpolation. Variable x is bogus here, but
% we fill it in to keep the accountants upstream happy.
if (fitType == 6)
disp('Fitting with linear interpolation');
fit_out6 = zeros(nOut,nDevices);
for i = 1:nDevices
[fit_out6(:,i)] = ...
FitGammaLin(values_in(:,i),measurements(:,i),values_out);
end
x6 = [];
end
% Cubic spline. Variable x is bogus here, but
% we fill it in to keep the accountants upstream happy.
if (fitType == 7)
disp('Fitting with cubic spline');
fit_out7 = zeros(nOut,nDevices);
for i = 1:nDevices
[fit_out7(:,i)] = ...
FitGammaSpline(values_in(:,i),measurements(:,i),values_out);
end
x7 = [];
end
% If we are not forcing a fit type, find best fit.
% Don't check linear interpolation, as it has zero error always.
% Currently we take the minimum mean error over all devices.
% In principle, could use best fit type for each device. But
% that would make the interface tricky.
if (fitType == 0)
meanErr = mean(errv,2);
[minErr,bestFit] = min(meanErr);
fitType = bestFit;
end
if (fitType == 1)
fit_out = fit_out1;
x = x1;
fitComment = (sprintf('Simple power function fit, RMSE: %g',...
mean(errv(1,:))));
elseif (fitType == 2)
fit_out = fit_out2;
x = x2;
fitComment = (sprintf('Extended power function fit, RMSE: %g',...
mean(errv(2,:))));
elseif (fitType == 3)
fit_out = fit_out3;
x = x3;
fitComment = (sprintf('Sigmoidal fit, RMSE: %g',...
mean(errv(3,:))));
elseif (fitType == 4)
fit_out = fit_out4;
x = x4;
fitComment = (sprintf('Weibull fit, RMSE: %g',...
mean(errv(4,:))));
elseif (fitType == 5)
fit_out = fit_out5;
x = x5;
fitComment = (sprintf('Polynomial fit, RMSE: %g',...
mean(errv(5,:))));
elseif (fitType == 6)
fit_out = fit_out6;
x = x6;
fitComment = (sprintf('Linear interpolation fit'));
elseif (fitType == 7)
fit_out = fit_out7;
x = x7;
fitComment = (sprintf('Cubic spline fit'));
end
% Check that fit is non-decreasing
if (CheckMonotonic(fit_out) == 0)
disp('Warning, fit is not non-decreasing');
end