Mercurial > hg > ltpda
diff m-toolbox/classes/@matrix/linfitsvd.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/linfitsvd.m Wed Nov 23 19:22:13 2011 +0100 @@ -0,0 +1,708 @@ +% LINFITSVD Linear fit with singular value decomposition +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% DESCRIPTION: Linear least square problem with singular value +% decomposition +% +% ALGORITHM: 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(os1,...,osn,pl); +% +% INPUT: +% - osi are vector of system output signals. They must be +% Nx1 matrix objects, where N is the output dimension of the +% system +% +% OUTPUT: +% - pars: a pest object containing parameter estimation +% +% <a href="matlab:utils.helper.displayMethodInfo('matrix', 'linfitsvd')">Parameters Description</a> +% +% VERSION: $Id: linfitsvd.m,v 1.38 2011/04/08 08:56:32 hewitson Exp $ +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function varargout = linfitsvd(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()); + + %%% get input parameters + % if the model is a matrix of smodels + lmod = find(pl, 'dmodel'); + % if the model is ssm + inNames = find(pl,'InNames'); + outNames = find(pl,'OutNames'); + % common parameters + mod = find(pl,'Model'); + fitparams = find(pl,'FitParams'); + inputdat = find(pl,'Input'); + WF = find(pl,'WhiteningFilter'); + nloops = find(pl,'Nloops'); + ncut = find(pl,'Ncut'); + npad = find(pl,'Npad'); + kwnpars = find(pl,'KnownParams'); + tol = find(pl,'tol'); + fastopt = find(pl,'fast'); + setalias = find(pl,'SetAlias'); + sThreshold = find(pl,'sThreshold'); + diffStep = find(pl,'diffStep'); + + boundedpars = find(pl,'BoundedParams'); + boudary = find(pl,'BoundVals'); + + % check if there are bounded parameters + if ~isempty(boundedpars) + boundparams = true; + else + boundparams = false; + end + + + % check the class of the model + if isa(mod,'ssm') + ssmmod = true; + else + ssmmod = false; + end + + %%% some sanity checks + if ~ssmmod + if numel(mtxs)~= numel(inputdat.objs) + error('Number of input data vectors must be the same of fit data vectors') + end + end + + + [fitparams,idx] = sort(fitparams); + fitparvals = fitparams; + if ssmmod + %%% get fit parameters + for ii=1:numel(fitparams) + fitparvals{ii} = mod.getParameters(plist('names',fitparams{ii})); + end + if isempty(diffStep) + sdiffStep = cell2mat(fitparvals).*0.01; + idz = sdiffStep == 0; + sdiffStep(idz) = 1e-7; + else + sdiffStep = diffStep(idx); + end + else + %%% get a single set of parameters + totparnames = {}; + totparvals = {}; + for ii=1:numel(mod.objs) + aa = mod.objs(ii).params; + cc = mod.objs(ii).values; + % get total parameter names + [bb,i1,i2]=union(totparnames,aa); + totparnames = [totparnames(i1),aa(i2)]; + % get total parameter values + totparvals = [totparvals(i1),cc(i2)]; + end + [totparnames,id] = sort(totparnames); + totparvals = totparvals(id); + + %%% get fit parameters + [nn,i1,i2] = intersect(totparnames,fitparams); + fitparams = totparnames(i1); + fitparvals = totparvals(i1); + [fitparams,id] = sort(fitparams); + fitparvals = fitparvals(id); + end + + if ~ssmmod + %%% linearize model with respect to fit parameters + if isempty(lmod) + lmod = linearize(mod,plist('Params',fitparams,'Sorting',false)); + end + end + + if isempty(WF) + wfdat = copy(mtxs,1); + elseif ~fastopt + %%% whitening fit data + wfdat = copy(mtxs,1); + for ii=1:numel(mtxs) + wfdat(ii) = filter(mtxs(ii),WF); + end + end + + % decide to pad in any case, assuming the objects have the same length + if isempty(npad) + npad = length(mtxs(1).objs(1).data.y) - 1; + end + + % set alias if there are + if setalias && (~ssmmod && ~isempty(mod.objs(1).aliasNames)) + nsecs = mtxs(1).objs(1).data.nsecs; + fs = mtxs(1).objs(1).data.fs; + + plalias = plist('nsecs',nsecs,'npad',npad,'fs',fs); + for ii=1:numel(mod.objs) + mod.objs(ii).assignalias(mod.objs(ii),plalias); + end + for jj=1:numel(lmod.objs) + for ii=1:numel(lmod.objs{jj}.objs) + lmod.objs{jj}.objs(ii).assignalias(lmod.objs{jj}.objs(ii),plalias); + end + end + end + + % do a copy to add at the output pest + outmod = copy(mod,1); + + % check if the fast option is active + if ~ssmmod && fastopt + % set length of fft (this should match the operation made in fftfilt_core) + nfft = length(mtxs(1).objs(1).data.y) + npad; + fs = mtxs(1).objs(1).data.fs; + % get fft freqs for current data. type option must match the one used + % in fftfilt_core for fft_core + fftfreq = utils.math.getfftfreq(nfft,fs,'one'); + % calculate freq resp of diagonal elements of WF + rWF = getWFresp(WF,fftfreq,fs); + % combine symbolic models with rWF + mod = joinmodelandfilter(mod,rWF); + lmod = joinmodelandfilter(lmod,rWF); + WF = []; + %%% whitening fit data + wfdat = copy(mtxs,1); + for ii=1:numel(mtxs) + for jj=1:numel(mtxs(ii).objs) + wfdat(ii).objs(jj) = fftfilt_core(wfdat(ii).objs(jj),rWF.objs(jj),npad); + end + end + clear rWF + end + + % init storage struct + loop_results = struct('a',cell(1),... + 'Ca',cell(1),... + 'Corra',cell(1),... + 'Vu',cell(1),... + 'bu',cell(1),... + 'Cbu',cell(1),... + 'Fbu',cell(1),... + 'mse',cell(1),... + 'params',cell(1),... + 'ppm',cell(1)); + + % init user interaction variable + reply = 'N'; + + %%% run fit loop + + % This causes problems on some machines so we remove it for now until we + % can investigate further. + % fftw('planner', 'exhaustive'); + + for kk=1:nloops + + % init index variable + xxx = 1; + + % init data struct + exps = struct(); + + %%% Set fit parameters into model + if ssmmod + fs = wfdat(1).objs(1).fs; + mod.doSetParameters(fitparams,cell2mat(fitparvals)); + lmod = parameterDiff(mod,plist('names',fitparams,'values',sdiffStep)); + lmod.modifyTimeStep(plist('newtimestep',1/fs)); + else + % fitparvals are updated at each fit loop + if fastopt + for ii = 1:numel(mod.objs) + mod.objs(ii).objs{2}.setParams(fitparams,fitparvals); + end + else + for ii = 1:numel(mod.objs) + mod.objs(ii).setParams(fitparams,fitparvals); + end + end + end + + %%% run over input data + + for ii=1:numel(inputdat.objs) + if ssmmod + %%% extract input + if isa(inputdat.objs{ii},'ao') + in = inputdat.objs{ii}; + elseif isa(inputdat.objs{ii},'matrix') + in = inputdat.objs{ii}.objs(:); + else + error('Unknown Input data type.') + end + + %%% calculates zero order response + plsym = plist('AOS VARIABLE NAMES',inNames{ii},... + 'RETURN OUTPUTS',outNames{ii},... + 'AOS',in); + %eo = simulate(lmod,plsym); + %zor = matrix(eo,plist('shape',[numel(eo) 1])); + zor = simulate(lmod,plsym); + % check dimensions + if size(zor.objs,1)<size(zor.objs,2) + % do transpose + zor = transpose(zor); + end + + %%% calculates first order response + for jj=1:numel(fitparams) + % get output ports names + [token, remain] = strtok(outNames{ii},'.'); + loutNames = token; + for zz=1:numel(token) + loutNames{zz} = sprintf('%s_DIFF_%s%s',token{zz},fitparams{jj},remain{zz}); + end + plsym = plist('AOS VARIABLE NAMES',inNames{ii},... + 'RETURN OUTPUTS',loutNames,... + 'AOS',in); + %eol = simulate(lmod,plsym); + %fstor(jj) = matrix(eol,plist('shape',[numel(eol) 1])); + fstor(jj) = simulate(lmod,plsym); + % check dimensions + if size(fstor(jj).objs,1)<size(fstor(jj).objs,2) + % do transpose + fstor(jj) = transpose(fstor(jj)); + end + fstor(jj).setName(fitparams{jj}); + end + else + %%% calculates zero order response + zor = fftfilt(inputdat.objs{ii},mod,plist('Npad',npad)); + + %%% calculates first order response + for jj=1:numel(lmod.objs) + fstor(jj) = fftfilt(inputdat.objs{ii},lmod.objs{jj},plist('Npad',npad)); + fstor(jj).setName(lmod.objs{jj}.name); + end + end + + if isempty(WF) + wzor = zor; + wfstor = fstor; + else + %%% whitening zor + wzor = filter(zor,WF); + + %%% whitening fstor + for jj=1:numel(fstor) + wfstor(jj) = filter(fstor(jj),WF); + wfstor(jj).setName(fstor(jj).name); + end + end + + %%% Collect object for the fit procedure + for jj=1:numel(wfdat(ii).objs) + % get difference between fit data and zero order response + tfdat = wfdat(ii).objs(jj) - wzor.objs(jj); + % remove whitening filter transient + tfdats = tfdat.split(plist('samples',[ncut+1 numel(tfdat.y)])); + % insert into exps struct + fitdata(xxx,1) = tfdats; + + % build fit basis + for gg=1:numel(fitparams) + for hh=1:numel(wfstor) + if strcmp(fitparams(gg),wfstor(hh).name) + bsel = wfstor(hh).objs(jj); + % remove whitening filter transient + bsels = bsel.split(plist('samples',[ncut+1 numel(tfdat.y)])); + end + end + bs(gg) = bsels; + end + % insert basis + fitbasis(xxx,:) = bs; + + % step up xxx + xxx = xxx + 1; + end %jj=1:numel(wfdat(ii).objs) + + end %ii=1:numel(inputdat.objs) + + %%% build input objects + [NN,MM] = size(fitbasis); + + for zz=1:MM + H(1,zz) = matrix(fitbasis(:,zz)); + end + + Y = matrix(fitdata); + + %%% Insert known parameters + if ~isempty(kwnpars) + kwnparmanes = kwnpars.names; + kwnparvals = kwnpars.y; + kwnparerrors = kwnpars.dy; + + % init struct + groundexps = struct; + + for ii=1:numel(kwnparmanes) + for jj=1:numel(fitparams) + if strcmp(kwnparmanes{ii},fitparams{jj}) + groundexps(ii).pos = jj; + groundexps(ii).value = kwnparvals(ii); + groundexps(ii).err = kwnparerrors(ii); + end + end + end + end + + + %%% do fit + if ~isempty(kwnpars) && isfield(groundexps,'pos') + plfit = plist('KnownParams',groundexps,'sThreshold',sThreshold); + [out,a,Ca,Corra,Vu,bu,Cbu,Fbu,mse,dof,ppm] = linlsqsvd(H,Y,plfit); + else + plfit = plist('sThreshold',sThreshold); + [out,a,Ca,Corra,Vu,bu,Cbu,Fbu,mse,dof,ppm] = linlsqsvd(H,Y,plfit); + end + + %%% update parameters values + for ii=1:numel(fitparams) + fitparvals{ii} = fitparvals{ii} + a(ii); + end + + %%% check for bouded params + if boundparams + for pp=1:numel(fitparams) + for qq=1:numel(boundedpars) + if strcmp(fitparams{pp},boundedpars{qq}); + % check boudaries + bd = boudary{qq}; + if fitparvals{pp}<bd(1) + fitparvals{pp} = bd(1); + elseif fitparvals{pp}>bd(2) + fitparvals{pp} = bd(2); + end + end + end + end + end + + %%% store intermediate results + loop_results(kk).a = a; + loop_results(kk).Ca = Ca; + loop_results(kk).Corra = Corra; + loop_results(kk).Vu = Vu; + loop_results(kk).bu = bu; + loop_results(kk).Cbu = Cbu; + loop_results(kk).Fbu = Fbu; + loop_results(kk).mse = mse; + loop_results(kk).params = fitparvals; + loop_results(kk).ppm = ppm; + + utils.helper.msg(msg.IMPORTANT, 'loop %d, mse %d\n',kk,mse); + + % check fit stability and accuracy + fitmsg = checkfit(Vu,a); + if ~isempty(fitmsg) && ~strcmpi(reply,'Y') + % display message + utils.helper.msg(msg.IMPORTANT, fitmsg); + % decide if stop for cycle + reply = input('Do you want to carry on with fit iteration? Y/N [Y]: ', 's'); + if isempty(reply) + reply = 'Y'; + end + if strcmpi(reply,'N') + break + end + end + + % check convergence + condvec = (abs(a).^2)./diag(Ca); % parameters a are going to zero during the fit iterations + + if all(condvec < tol) + condmsg = sprintf(['Fit parameters have reached convergence.\n'... + 'Fit loop was terminated at iteration %s.\n'],num2str(kk)); + utils.helper.msg(msg.IMPORTANT, condmsg); + break + end + + + end %for kk=1:nloops + + %%% output data + % get minimum mse + if ~isempty(fitmsg) + val = mse; + idx = kk; + else + mseprog = zeros(numel(loop_results),1); + for ii=1:numel(loop_results) + mseprog(ii) = loop_results(ii).mse; + end + [val,idx] = min(abs(mseprog-1)); % get the nearest to 1 value + utils.helper.msg(msg.PROC1, 'Output values at minimum mse, mse = %d\n',val); + end + + % output pest object + pe = pest(); + pe.setY(cell2mat(loop_results(idx).params)); + pe.setDy(sqrt(diag(loop_results(idx).Ca))); + pe.setCov(loop_results(idx).Ca); + pe.setChi2(loop_results(idx).mse); + pe.setNames(fitparams); + pe.setDof(dof); + pe.setModels(outmod); + pe.setName(sprintf('linfitsvd(%s)', mod.name)); + pe.setProcinfo(plist('loop_results',loop_results)); + + % set History + pe.addHistory(getInfo('None'), pl, [mtxs_invars(:)], [inhists(:)]); + varargout{1} = pe; + + + +end + +%-------------------------------------------------------------------------- +% check fit accuracy and stability +%-------------------------------------------------------------------------- +function msg = checkfit(V,aa) + + if size(V,1)<numel(aa) + % The number of parameters combinations is less than the number of fit + % parameters. Information cannot be recostructed fit results will be + % compromised + VV = abs(V).^2; + num = numel(aa)-size(V,1); + mVV = max(VV); + % try to identify non measured params + unmparams = []; + for jj = 1:num + [vl,idx] = min(mVV); + unmparams = [unmparams idx]; + mVV(idx) = []; + end + msg = sprintf(['!!! The number of parameters combinations is less than the number of fit parameters. \n' ... + 'Information cannot be recostructed and fit results will be compromised. \n'... + 'Try to remove parameters %s from the fit parameters list or add information with more experiments !!!\n'],num2str(unmparams)); + + else + + unmparams = []; + trh1 = 1e-4; + + % eigenvectors are normalized, therefore square of the rows of V are sum + % to one. Each column of V store the coefficients for a given parameter + % for the set of eigenvectors + for jj = 1:size(V,2) + cl = abs(V(:,jj)).^2; + if all(cl<trh1) + unmparams = [unmparams jj]; + end + end + if ~isempty(unmparams) + msg = sprintf(['!!! Parameter/s %s is/are not well measured. \n'... + 'Fit accuracy could be impaired. \n'... + 'Try to remove such parameters from the fit parameters list !!!\n'],num2str(unmparams)); + else + msg = ''; + end + end + + + + + +end + +%-------------------------------------------------------------------------- +% calculate frequency response of diagonal elements of the whitening filter +%-------------------------------------------------------------------------- +function rsp = getWFresp(wf,f,fs) + % run over wf elements + obj = wf.objs; + [rw,cl] = size(obj); + if rw~=cl + error('??? Matrix of whitening filter must be square'); + end + amdl = ao.initObjectWithSize(rw,1); + for jj=1:rw + % check filter type + switch lower(class(obj(jj,jj))) + case 'filterbank' + % get filter response on given frequencies + amdly = utils.math.mtxiirresp(obj(jj,jj).filters,f,fs,obj(jj,jj).type); + amdl(jj,1) = ao(fsdata(f, amdly, fs)); + case 'miir' + % get filter response on given frequencies + amdly = utils.math.mtxiirresp(obj(jj,jj),f,fs,[]); + amdl(jj,1) = ao(fsdata(f, amdly, fs)); + end + end + rsp = matrix(amdl,plist('shape',[rw 1])); + +end + +%-------------------------------------------------------------------------- +% Join Symbolic model and whitening filter for fast calculations +%-------------------------------------------------------------------------- +function jmod = joinmodelandfilter(smod,fil) + switch class(smod) + case 'matrix' + mobj = smod.objs; + [nn,mm] = size(mobj); + nmobj = collection.initObjectWithSize(nn,mm); + for ii=1:nn + for jj=1:mm + nmobj(ii,jj) = collection(fil.objs(ii,1),mobj(ii,jj)); + nmobj(ii,jj).setName(mobj(ii,jj).name); + end + end + jmod = matrix(nmobj, plist('shape',[nn,mm])); + jmod.setName(smod.name); + case 'collection' + matobj = matrix.initObjectWithSize(1,numel(smod.objs)); + for kk=1:numel(smod.objs) + mobj = smod.objs{kk}.objs; + [nn,mm] = size(mobj); + nmobj = collection.initObjectWithSize(nn,mm); + for ii=1:nn + for jj=1:mm + nmobj(ii,jj) = collection(fil.objs(ii,1),mobj(ii,jj)); + nmobj(ii,jj).setName(mobj(ii,jj).name); + end + end + matobj(kk) = matrix(nmobj); + matobj(kk).setName(smod.objs{kk}.name); + %smod.objs{kk} = matobj; + end + jmod = collection(matobj); + 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, 'matrix', 'ltpda', utils.const.categories.sigproc, '$Id: linfitsvd.m,v 1.38 2011/04/08 08:56:32 hewitson Exp $', sets, pl); + ii.setArgsmin(1); + 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() + + % General plist for moltichannel fits + pl = plist.MCH_FIT_PLIST; + + +% p = param({'FitParams','A cell array with the names of the fit parameters'}, {}); +% pl.append(p); + + p = param({'BoundedParams','A cell array with the names of the bounded fit parameters'}, {}); + pl.append(p); + + p = param({'BoundVals','A cell array with the boundaries values for the bounded fit parameters'}, {}); + pl.append(p); + +% p = param({'Model','System parametric model. A matrix of smodel objects or a ssm object'}, paramValue.EMPTY_DOUBLE); +% pl.append(p); + + p = param({'dModel','Partial derivatives of the system parametric model. A matrix of smodel objects'}, paramValue.EMPTY_DOUBLE); + pl.append(p); + +% p = param({'InNames','A cell array containing cell arrays of the input ports names for each experiment. Used only with ssm models.'}, {}); +% pl.append(p); +% +% p = param({'OutNames','A cell array containing cell arrays of the output ports names for each experiment. Used only with ssm models.'}, {}); +% pl.append(p); + +% p = param({'Input','Collection of input signals'},paramValue.EMPTY_DOUBLE); +% pl.append(p); + + p = param({'WhiteningFilter','The multichannel whitening filter. A matrix object of filters'},paramValue.EMPTY_DOUBLE); + pl.append(p); + + p = param({'Nloops', 'Number of desired iteration loops.'}, paramValue.DOUBLE_VALUE(1)); + pl.append(p); + + p = param({'Ncut', 'Number of bins to be discharged in order to cut whitening filter transients'}, paramValue.DOUBLE_VALUE(100)); + pl.append(p); + + % Number of points for zero padding + p = param({'Npad', 'Number of points for zero padding.'}, paramValue.EMPTY_DOUBLE); + pl.append(p); + + p = param({'KnownParams', 'Known Parameters. A pest object containing parameters values, names and errors'}, paramValue.EMPTY_DOUBLE); + pl.append(p); + + p = param({'tol','Convergence threshold for fit parameters'}, paramValue.DOUBLE_VALUE(1)); + pl.append(p); + + p = param({'fast',['Using fast option causes the whitening filter to be applied in frequency domain.'... + 'The filter matrix is considered diagonal. The method skip time domain filtering saving some process time'... + 'It works only when the imput model is a matrix of smodels']}, paramValue.FALSE_TRUE); + pl.append(p); + + p = param({'SetAlias','Set to true in order to aassign internally the values to the model alias'}, paramValue.FALSE_TRUE); + 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.DOUBLE_VALUE(1)); + pl.append(p); + + p = plist({'diffStep','Numerical differentiation step for ssm models'}, paramValue.EMPTY_DOUBLE); + pl.append(p); + + + +end