Mercurial > hg > ltpda
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