diff m-toolbox/classes/@matrix/linlsqsvd.m @ 0:f0afece42f48

Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@matrix/linlsqsvd.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,196 @@
+% LINLSQSVD Linear least squares with singular value decomposition
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% DESCRIPTION: Linear least square problem with singular value
+% decomposition
+%
+% ALGORITHM: % It solves the problem
+%
+%        Y = HX
+%
+% where X are the parameters, Y the measurements, and H the linear
+% equations relating the two.
+% It is able to perform linear identification of the parameters of a
+% multichannel systems. The results of different experiments on the same
+% system can be passed as input. The algorithm, thanks to the singular
+% value decomposition, extract the maximum amount of information from each
+% single channel and for each experiment. Total information is then
+% combined to get the final result.
+%            
+% CALL:                   pars = linfitsvd(H1,...,HN,Y,pl);
+% 
+% If the experiment is 1 then H1,...,HN and Y are aos.
+% If the experiments are M, then H1,...,HN and Y are Mx1 matrix objects
+% with the aos relating to the given experiment in the proper position.
+% 
+% INPUT:
+%               - Hi represent the columns of H
+%               - Y represent the measurement set
+% 
+% OUTPUT:
+%               - pars: a pest object containing parameter estimation
+% 
+% 09-11-2010 L Ferraioli
+%       CREATION
+%
+% <a href="matlab:utils.helper.displayMethodInfo('matrix', 'linfitsvd')">Parameters Description</a>
+%
+% VERSION:     $Id: linlsqsvd.m,v 1.9 2011/05/02 14:18:19 luigi Exp $
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function varargout = linlsqsvd(varargin)
+  
+  %%% LTPDA stufs and get data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  
+  % Check if this is a call for parameters
+  if utils.helper.isinfocall(varargin{:})
+    varargout{1} = getInfo(varargin{3});
+    return
+  end
+  
+  import utils.const.*
+  utils.helper.msg(msg.OMNAME, 'running %s/%s', mfilename('class'), mfilename);
+  
+  % Collect input variable names
+  in_names = cell(size(varargin));
+  for ii = 1:nargin,in_names{ii} = inputname(ii);end
+  
+  % Collect all ltpdauoh objects
+  [mtxs, mtxs_invars] = utils.helper.collect_objects(varargin(:), 'matrix', in_names);
+  [pl, invars] = utils.helper.collect_objects(varargin(:), 'plist');
+  
+  
+  inhists = [mtxs(:).hist];
+
+  
+  %%% combine plists
+  pl = parse(pl, getDefaultPlist());
+  
+  %%% collect inputs names
+  argsname = mtxs(end).name;
+  
+  %%% get input params
+  kwnpars    = find(pl,'KnownParams');
+  sThreshold = find(pl,'sThreshold');
+  
+  
+  %%% do fit
+%   if ~isempty(kwnpars) && isfield(kwnpars,'pos')
+%     [a,Ca,Corra,Vu,bu,Cbu,Fbu,mse,dof,ppm] = utils.math.linlsqsvd(mtxs,kwnpars);
+%   else
+%     [a,Ca,Corra,Vu,bu,Cbu,Fbu,mse,dof,ppm] = utils.math.linlsqsvd(mtxs);
+%   end
+  
+  if ~isempty(kwnpars) && isfield(kwnpars,'pos')
+    if ~isempty(sThreshold)
+      [a,Ca,Corra,Vu,bu,Cbu,Fbu,mse,dof,ppm] = utils.math.linlsqsvd(mtxs,sThreshold,kwnpars);
+    else
+      [a,Ca,Corra,Vu,bu,Cbu,Fbu,mse,dof,ppm] = utils.math.linlsqsvd(mtxs,kwnpars);
+    end
+  else
+    if ~isempty(sThreshold)
+      [a,Ca,Corra,Vu,bu,Cbu,Fbu,mse,dof,ppm] = utils.math.linlsqsvd(mtxs,sThreshold);
+    else
+      [a,Ca,Corra,Vu,bu,Cbu,Fbu,mse,dof,ppm] = utils.math.linlsqsvd(mtxs);
+    end
+  end
+  
+  fitparams = cell(1,numel(a));
+  nmstr = '';
+  for kk=1:numel(a)
+    fitparams{kk} = sprintf('a%s',num2str(kk));
+    units{kk} = mtxs(end).objs(1).yunits / mtxs(kk).objs(1).yunits;
+    units{kk}.simplify;
+    if isempty(nmstr)
+      nmstr = sprintf('%s*%s',fitparams{kk},mtxs(kk).name);
+    else
+      nmstr = [nmstr '+' sprintf('%s*%s',fitparams{kk},mtxs(kk).name)];
+    end
+  end
+  
+  pe = pest();
+  pe.setY(a);
+  pe.setDy(sqrt(diag(Ca)));
+  pe.setCov(Ca);
+  pe.setChi2(mse);
+  pe.setNames(fitparams);
+  pe.setDof(dof);
+  pe.setYunits(units{:});
+  pe.name = nmstr;
+  pe.setModels(mtxs(1:end-1));
+  
+  % set History
+  pe.addHistory(getInfo('None'), pl, [mtxs_invars(:)], [inhists(:)]);
+
+  
+  if nargout == 1
+    varargout{1} = pe;
+  elseif nargout == 11
+    varargout{1} = pe;
+    varargout{2} = a;
+    varargout{3} = Ca;
+    varargout{4} = Corra;
+    varargout{5} = Vu;
+    varargout{6} = bu;
+    varargout{7} = Cbu;
+    varargout{8} = Fbu;
+    varargout{9} = mse;
+    varargout{10} = dof;
+    varargout{11} = ppm;
+  else
+    error('invalid number of outputs!')
+  end
+  
+end
+    
+
+%--------------------------------------------------------------------------
+% Get Info Object
+%--------------------------------------------------------------------------
+function ii = getInfo(varargin)
+  if nargin == 1 && strcmpi(varargin{1}, 'None')
+    sets = {};
+    pl   = [];
+  else
+    sets = {'Default'};
+    pl   = getDefaultPlist;
+  end
+  % Build info object
+  ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: linlsqsvd.m,v 1.9 2011/05/02 14:18:19 luigi Exp $', sets, pl);
+  ii.setArgsmin(2);
+  ii.setOutmin(1);
+%   ii.setOutmax(1);
+  ii.setModifier(false);
+end
+
+%--------------------------------------------------------------------------
+% Get Default Plist
+%--------------------------------------------------------------------------
+function plout = getDefaultPlist()
+  persistent pl;  
+  if exist('pl', 'var')==0 || isempty(pl)
+    pl = buildplist();
+  end
+  plout = pl;  
+end
+
+function pl = buildplist()
+  
+  pl = plist();
+  
+  p = param({'KnownParams', ['Known Parameters. A struct array with the fields:<ul>'...
+    '<li> pos - a number indicating the corresponding position of the parameter (corresponding column of H)</li>'...
+    '<li> value - the value for the parameter</li>'...
+    '<li> err - the uncertainty associated to the parameter</li>'...
+    '</ul>']}, paramValue.EMPTY_CELL);
+  pl.append(p);
+  
+  p = param({'sThreshold',['Fix upper treshold for singular values.'...
+    'Singular values larger than the value will be ignored.'...
+    'This correspon to consider only parameters combinations with error lower then the value']},...
+    paramValue.EMPTY_DOUBLE);
+  pl.append(p);
+  
+  
+end