Home > m > timetools > statespacefunctions > ltpda_ss_modify.m

ltpda_ss_modify

PURPOSE ^

ltpda_ss_modify enables to modify parameters and retrieve numerical

SYNOPSIS ^

function varargout = ltpda_ss_modify(varargin)

DESCRIPTION ^

 ltpda_ss_modify enables to modify parameters and retrieve numerical
 values

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 DESCRIPTION: ltpda_ss_modify enables to modify parameters and retrieve
 system data in numrical or parametric form, eventually take the
 jacobian or hessian of its fields.

 CALL: varargout = ltpda_ss_modify(varargin)
 [subystem] = paramSubstitute(subystem, operations)
 
 INPUTS: subystem - (array) of plist describing a subsystem (refer to
 ltpda_ss_check for details) operations - plist of elements for numerical
 values the of parameters 
 'NOT_USED_PARAM_NAMES' cell array of parameter discarted in this
 analysis, substituted by their standard numerical value
 'NOT_USED_PARAM_ORDER' array of 1 / 0 depending whether parameter are
 effectively discarted or not
 
 'SET_PARAM_NAMES' Cell array of parameter names for numerical substitutions with user defined values
 'SET_PARAM_VALUES' Array of parameter means for numerical substitutions
 'SET_PARAM_SIGMAS' Array of parameter standard deviation for numerical substitutions (unit may be none!!!)
 'SET_PARAM_VALUES_ORDER' Array of 1 / 0 depending whether parameter means are effectively substituted or not
 'SET_PARAM_SIGMAS_ORDER' Array of 1 / 0 depending whether parameter means are effectively substituted or not
 
 'JAC_PARAM_NAMES' same for parameters used for evaluating a jacobian
 'JAC_PARAM_ORDER'
 'HESS_PARAM_NAMES' same for parameters used for evaluating a hessian
 'HESS_PARAM_ORDER' 
 'SYMBOLIC_OUTPUT_ALLOWED' if 0 checks output tensors are double and not symbollic arrays
 
 OUTPUTS: subystem - (array of) plist describing the subsystem (refer to
 test_subsys for details) with substitutions / modifications done 
 Here jacobians and hessians are stored within the AMAT,BMAT... arrays
 that are augmented by one dimension

 in both cases differenciated array are the same size for the first two
 dimensions, the third (resp. 4rth) one being for the direction of the
 differenciation
 The output array are of size Nsys*Ndiff

 VERSION: $Id: ltpda_ss_modify.m,v 1.5 2008/03/11 16:52:56 adrien Exp $

 HISTORY: 26-02-2008 A Grynagier
 29-01-2008 A Grynagier
 Creation

 % TO DO : 
 add reshpae for BMATS and DMATS at the beggining
 update help header (orders 'ALL', 'NO')
 display in window parametric substitutions and differenciation proceeded
 Define user scenarii to see ensure output is satisfying
 Various numerical checks

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function varargout = ltpda_ss_modify(varargin)
0002 % ltpda_ss_modify enables to modify parameters and retrieve numerical
0003 % values
0004 %
0005 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0006 %
0007 % DESCRIPTION: ltpda_ss_modify enables to modify parameters and retrieve
0008 % system data in numrical or parametric form, eventually take the
0009 % jacobian or hessian of its fields.
0010 %
0011 % CALL: varargout = ltpda_ss_modify(varargin)
0012 % [subystem] = paramSubstitute(subystem, operations)
0013 %
0014 % INPUTS: subystem - (array) of plist describing a subsystem (refer to
0015 % ltpda_ss_check for details) operations - plist of elements for numerical
0016 % values the of parameters
0017 % 'NOT_USED_PARAM_NAMES' cell array of parameter discarted in this
0018 % analysis, substituted by their standard numerical value
0019 % 'NOT_USED_PARAM_ORDER' array of 1 / 0 depending whether parameter are
0020 % effectively discarted or not
0021 %
0022 % 'SET_PARAM_NAMES' Cell array of parameter names for numerical substitutions with user defined values
0023 % 'SET_PARAM_VALUES' Array of parameter means for numerical substitutions
0024 % 'SET_PARAM_SIGMAS' Array of parameter standard deviation for numerical substitutions (unit may be none!!!)
0025 % 'SET_PARAM_VALUES_ORDER' Array of 1 / 0 depending whether parameter means are effectively substituted or not
0026 % 'SET_PARAM_SIGMAS_ORDER' Array of 1 / 0 depending whether parameter means are effectively substituted or not
0027 %
0028 % 'JAC_PARAM_NAMES' same for parameters used for evaluating a jacobian
0029 % 'JAC_PARAM_ORDER'
0030 % 'HESS_PARAM_NAMES' same for parameters used for evaluating a hessian
0031 % 'HESS_PARAM_ORDER'
0032 % 'SYMBOLIC_OUTPUT_ALLOWED' if 0 checks output tensors are double and not symbollic arrays
0033 %
0034 % OUTPUTS: subystem - (array of) plist describing the subsystem (refer to
0035 % test_subsys for details) with substitutions / modifications done
0036 % Here jacobians and hessians are stored within the AMAT,BMAT... arrays
0037 % that are augmented by one dimension
0038 %
0039 % in both cases differenciated array are the same size for the first two
0040 % dimensions, the third (resp. 4rth) one being for the direction of the
0041 % differenciation
0042 % The output array are of size Nsys*Ndiff
0043 %
0044 % VERSION: $Id: ltpda_ss_modify.m,v 1.5 2008/03/11 16:52:56 adrien Exp $
0045 %
0046 % HISTORY: 26-02-2008 A Grynagier
0047 % 29-01-2008 A Grynagier
0048 % Creation
0049 %
0050 % % TO DO :
0051 % add reshpae for BMATS and DMATS at the beggining
0052 % update help header (orders 'ALL', 'NO')
0053 % display in window parametric substitutions and differenciation proceeded
0054 % Define user scenarii to see ensure output is satisfying
0055 % Various numerical checks
0056 
0057 %% standard calls for LTPDA function data
0058 
0059 ALGONAME = mfilename;
0060 VERSION = '$Id: ltpda_ss_modify.m,v 1.5 2008/03/11 16:52:56 adrien Exp $';
0061 CATEGORY = 'STATESPACE';
0062 display(['starting ' ALGONAME]);
0063 
0064 if not(isempty(varargin))
0065     if isequal( varargin{1}, 'Version')
0066         varargout = VERSION;
0067         return;
0068     elseif isequal(varargin{1}, 'Params')
0069         varargout = plist( 'NOT_USED_PARAM_NAMES', {'param1', 'param2'}, 'NOT_USED_PARAM_ORDER', [1 0 ],...
0070             'SET_PARAM_NAMES', {'param2', 'param3'}, 'SET_PARAM_VALUES', [3 6], 'SET_PARAM_SIGMAS', [0.2 0.1],...
0071             'SET_PARAM_VALUES_ORDER', 'ALL', 'SET_PARAM_SIGMAS_ORDER', 'NO', ...
0072             'JAC_PARAM_NAMES', {'param2', 'param4'}, 'JAC_PARAM_ORDER', [1 1],...
0073             'HESS_PARAM_NAMES', {'param2', 'param3'}, 'HESS_PARAM_ORDER' , [1 0],...
0074             'SYMBOLIC_OUTPUT_ALLOWED', 1);
0075         return;
0076     elseif isequal(varargin{1}, 'Category')
0077         varargout = CATEGORY;
0078         return;
0079     end
0080 end
0081 
0082 if length(varargin)==2
0083     Systems = varargin{1};
0084     modifs = varargin{2};
0085     System_out = plist();
0086 else
0087     error('wrong number of inputs!')
0088 end
0089 
0090 for i_system = 1:length(Systems)
0091     for i_modif = 1 :length(modifs)
0092         System = Systems(i_system);
0093         modif = modifs(i_modif);
0094 
0095 %% Retriveing option fields for operations
0096 
0097 
0098         % fields from user options
0099         NotUsedParamNamesIni = find(modif,'NOT_USED_PARAM_NAMES') ;
0100         NotUsedParamOrder = find(modif,'NOT_USED_PARAM_ORDER') ;
0101         SetParamNamesIni = find(modif,'SET_PARAM_NAMES') ;
0102         SetParamValuesIni = find(modif,'SET_PARAM_VALUES') ;
0103         SetParamSigmasIni = find(modif,'SET_PARAM_SIGMAS') ;
0104         SetParamValuesOrder = find(modif,'SET_PARAM_VALUES_ORDER') ;
0105         SetParamSigmasOrder = find(modif,'SET_PARAM_SIGMAS_ORDER') ;
0106         JacParamNamesIni = find(modif,'JAC_PARAM_NAMES') ;
0107         JacParamOrder = find(modif,'JAC_PARAM_ORDER') ;
0108         HessParamNamesIni = find(modif,'HESS_PARAM_NAMES') ;
0109         HessParamOrder = find(modif,'HESS_PARAM_ORDER') ;
0110         SymbolicOutputAllowed = find(modif,'SYMBOLIC_OUTPUT_ALLOWED') ;
0111 
0112         % subsystem data matrices
0113         SystemAMat = find(System,'AMAT');
0114         SystemBMats = find(System,'BMATS');
0115         SystemCMat = find(System,'CMAT');
0116         SystemDMats = find(System,'DMATS');
0117         % subsystem parametric data
0118         SystemParamNames = find(System,'PARAMNAMES');
0119         SystemParamValues = find(System,'PARAMVALUES');
0120         SystemParamSigmas = find(System,'PARAMSIGMAS');
0121         SystemNbParams = length(SystemParamNames);
0122 
0123 %% building updated list for NotUsedParam...
0124         if isequal(NotUsedParamOrder,'ALL')
0125             NotUsedParamNames = NotUsedParamNamesIni;
0126         elseif isequal(NotUsedParamOrder,'NO')
0127             NotUsedParamNames = {};
0128         else
0129             NotUsedParamNames = {};
0130             j=1;
0131             for i=1:length(NotUsedParamNamesIni)
0132                 if NotUsedParamOrder(i)
0133                     NotUsedParamNames{j} = NotUsedParamNamesIni{i}; %#ok<AGROW>
0134                     j = j+1;
0135                 end
0136             end
0137         end
0138         NotUsedParamValues = [];
0139         for i=1:length(NotUsedParamNames)
0140             for j=1:length(SystemParamNames)
0141                 if isequal(NotUsedParamNames{i},SystemParamNames{j})
0142                     NotUsedParamValues(i)=SystemParamValues(j); %#ok<AGROW>
0143                 end
0144             end
0145         end
0146 %% building updated list for SetParam...
0147         if (isequal(SetParamValuesOrder,'ALL'))
0148             SetParamValuesOrder = ones(length(SetParamNamesIni),1);
0149         elseif (isequal(SetParamValuesOrder,'NO'))
0150             SetParamValuesOrder = zeros(length(SetParamNamesIni),1);
0151         end
0152         if (isequal(SetParamSigmasOrder,'ALL'))
0153             SetParamSigmasOrder = ones(length(SetParamNamesIni),1);
0154         elseif (isequal(SetParamSigmasOrder,'NO'))
0155             SetParamSigmasOrder = zeros(length(SetParamNamesIni),1);
0156         end
0157         SetParamNames = {};
0158         SetParamValues = [];
0159         SetParamSigmas = [];
0160         j =1;
0161         for i=1:length(SetParamNamesIni)
0162             if SetParamValuesOrder(i) || SetParamSigmasOrder(i)
0163                 SetParamNames{j} = SetParamNamesIni{i}; %#ok<AGROW>
0164                 if SetParamValuesOrder(i)
0165                     SetParamValues(j) = SetParamValuesIni(i); %#ok<AGROW>
0166                 else
0167                     for k=1:SystemNbParams
0168                         if isequal(SystemParamNames{k}, SetParamNamesIni{i})
0169                             SetParamValues(j) = SystemParamValues(k); %#ok<AGROW>
0170                         end
0171                     end
0172                 end
0173                 if   SetParamSigmasOrder(i)
0174                     SetParamSigmas(j) = SetParamSigmasIni(i); %#ok<AGROW>
0175                 else
0176                     for k=1:SystemNbParams
0177                         if isequal(SystemParamNames{k}, SetParamNamesIni{i})
0178                             SetParamSigmas(j) = SystemParamSigmas(k); %#ok<AGROW>
0179                         end
0180                     end
0181                 end
0182                 j = j+1;
0183             end
0184         end
0185 %% building updated list for JacParam...
0186         if isequal(JacParamOrder,'ALL')
0187             JacParamNames = JacParamNamesIni;
0188         elseif isequal(JacParamOrder,'NO')
0189             JacParamNames = {};
0190         else
0191             JacParamNames = {};
0192             j=1;
0193             for i=1:length(JacParamNamesIni)
0194                 if JacParamOrder(i)
0195                     JacParamNames{j} = JacParamNamesIni{i}; %#ok<AGROW>
0196                     j = j+1;
0197                 end
0198             end
0199         end
0200 %% building updated list for HessParam...
0201         if isequal(HessParamOrder,'ALL')
0202             HessParamNames = HessParamNamesIni;
0203         elseif isequal(HessParamOrder,'NO')
0204             HessParamNames = {};
0205         else
0206             HessParamNames = {};
0207             j=1;
0208             for i=1:length(HessParamNamesIni)
0209                 if HessParamOrder(i)
0210                     HessParamNames{j} = HessParamNamesIni{i}; %#ok<AGROW>
0211                     j = j+1;
0212                 end
0213             end
0214         end
0215         nNotUsed = length(NotUsedParamNames);
0216         nSet = length(SetParamNames);
0217         nJac = length(JacParamNames);
0218         nHess = length(HessParamNames);
0219 
0220         %% building param position tables
0221         InSet = zeros(SystemNbParams,1);
0222         InJac = zeros(SystemNbParams,1);
0223         InHess = zeros(SystemNbParams,1);
0224         InNotUsed = zeros(SystemNbParams,1);
0225         for i=1:SystemNbParams
0226             for j=1:nNotUsed
0227                 if isequal(SystemParamNames{i},NotUsedParamNames{j})
0228                     InNotUsed(i) = j;
0229                 end
0230             end
0231             for j=1:nSet
0232                 if isequal(SystemParamNames{i},SetParamNames{j})
0233                     InSet(i) = j;
0234                 end
0235             end
0236             for j=1:nJac
0237                 if isequal(SystemParamNames{i},JacParamNames{j})
0238                     InJac(i) = j;
0239                 end
0240             end
0241             for j=1:nHess
0242                 if isequal(SystemParamNames{i},HessParamNames{j})
0243                     InHess(i) = j;
0244                 end
0245             end
0246         end
0247 
0248 %% looking for non exisitng but modified parameter names
0249         if not(length(NotUsedParamNames)==sum(InNotUsed>0))
0250             msg = 'error with notusedparameter list because one does not exist in system.';
0251             error(msg);
0252         elseif not(length(SetParamNames)==sum(InSet>0))
0253             msg = 'error with setparameter list because one does not exist in system.';
0254             error(msg);
0255         elseif not(length(JacParamNames)==sum(InJac>0))
0256             msg = 'error with jacparameter list because one does not exist in system.';
0257             error(msg);
0258         elseif not(length(HessParamNames)==sum(InHess>0))
0259             msg = 'error with hessparameter list because one does not exist in system.';
0260             error(msg);
0261         end
0262 
0263 %% looking of unused but used parameters
0264         for i=1:nNotUsed
0265             for j =1:nSet
0266                 if isequal(NotUsedParamNames{i},SetParamNames{j})
0267                     msg = 'error because parameter set was supposed to be unused';
0268                     error(msg);
0269                 end
0270             end
0271             for j =1:nJac
0272                 if isequal(NotUsedParamNames{i},JacParamNames{j})
0273                     msg = 'error because parameter for jacobian was supposed to be unused';
0274                     error(msg);
0275                 end
0276             end
0277             for j =1:nHess
0278                 if isequal(NotUsedParamNames{i},HessParamNames{j})
0279                     msg = 'error because parameter for hessian was supposed to be unused';
0280                     error(msg);
0281                 end
0282             end
0283         end
0284 
0285 %% Substituting unused parameters
0286         for i=1:nNotUsed
0287             if InNotUsed(i)>0
0288                 SystemAMat{1} = subs(SystemAMat{1}, NotUsedParamNames(InNotUsed(i)) , NotUsedParamValues(InNotUsed(i)));
0289                 SystemCMat{1} = subs(SystemCMat{1}, NotUsedParamNames(InNotUsed(i)) , NotUsedParamValues(InNotUsed(i)));
0290                 for j=1:length(SystemBMats)
0291                     SystemBMats{j} = subs(SystemBMats{j}, NotUsedParamNames(InNotUsed(i)) , NotUsedParamValues(InNotUsed(i)));
0292                     SystemDMats{j} = subs(SystemDMats{j}, NotUsedParamNames(InNotUsed(i)) , NotUsedParamValues(InNotUsed(i)));
0293                 end
0294             end
0295         end
0296 
0297 %% modifying parameter tables in System (setparam list)
0298         for i=1:SystemNbParams
0299             if InSet(i)>0
0300                 SystemParamValues(i) = SetParamValues(InSet(i));
0301                 SystemParamSigmas(i) = SetParamSigmas(InSet(i));
0302             end
0303         end
0304 %% Choosing parameter lists for substitution
0305         if not(SymbolicOutputAllowed)
0306             display('Forbidden symbollic output case is not handled in the current version' );
0307             SubsNames = SystemParamNames;
0308             SubsValues = SystemParamValues;
0309         else
0310             SubsNames = SetParamNames;
0311             SubsValues = SetParamValues;
0312         end
0313 
0314 
0315 %% building jacobians - symbollic expression
0316         % creating symbollic variables vector
0317         if nJac>0
0318             JacParamVar = zeros(nJac,1);
0319             for i = 1:length(JacParamNames2)
0320                 nameVar = eval(['JacParamNames{', num2str(i),'}']);
0321                 cmd = [nameVar ,'=sym(', nameVar ,');']; %Declaration not needed !!
0322                 eval(cmd);
0323                 cmd2 = ['JacParamVar(',num2str(i),')=',nameVar];
0324                 eval(cmd2);
0325             end
0326             % differenciating
0327             SystemAMat{1,2} = subs( ltpda_jacobian(SystemAMat{1}, JacParamVar,1), SubsNames , SubsValues);
0328             SystemCMat{1,2} = subs( ltpda_jacobian(SystemCMat{1}, JacParamVar,1), SubsNames , SubsValues);
0329             for j=1:length(SystemBMats)
0330                 SystemBMats{j,2} = subs( ltpda_jacobian(SystemBMats{j}, JacParamVar,1), SubsNames , SubsValues);
0331                 SystemDMats{j,2} = subs( ltpda_jacobian(SystemDMats{j}, JacParamVar,1), SubsNames , SubsValues);
0332             end
0333         end
0334 
0335 %% building hessians - symbollic expression
0336 
0337         % creating symbollic variables vector
0338         if nHess>0
0339             HessParamVar = zeros(nHess,1);
0340             for j = 1:nHess
0341                 nameVar = eval(['HessParamNames2{', num2str(j),'}']);
0342                 cmd = [nameVar ,'=sym(', nameVar ,');']; %Declaration not needed !!
0343                 eval(cmd);
0344                 cmd2 = ['HessParamVar(',num2str(j),')=',nameVar];
0345                 eval(cmd2);
0346             end
0347             % differenciating
0348             SystemAMat{1,3} = subs( ltpda_jacobian(SystemAMat{1}, HessParamVar,2), SubsNames , SubsValues);
0349             SystemCMat{1,3} = subs( ltpda_jacobian(SystemCMat{1}, HessParamVar,2), SubsNames , SubsValues);
0350             for j=1:length(SystemBMats)
0351                 SystemBMats{j,3} = subs( ltpda_jacobian(SystemBMats{j}, HessParamVar,2), SubsNames , SubsValues);
0352                 SystemDMats{j,3} = subs( ltpda_jacobian(SystemDMats{j}, HessParamVar,2), SubsNames , SubsValues);
0353             end
0354         end
0355 
0356 %% Numerical value for system
0357 
0358         if nSet>0
0359             SystemAMat{1,1} = subs(SystemAMat{1}, SubsNames , SubsValues);
0360             SystemCMat{1,1} = subs(SystemCMat{1}, SubsNames , SubsValues);
0361             for j=1:length(SystemBMats)
0362                 SystemBMats{j,1} = subs(SystemBMats{j}, SubsNames , SubsValues);
0363                 SystemDMats{j,1} = subs(SystemDMats{j}, SubsNames , SubsValues);
0364             end
0365         end
0366 
0367 %% parsing into new system
0368         System_out(i_system, i_modif) = pset(System, 'AMAT',SystemAMat, 'BMATS', SystemBMats, 'CMAT',SystemCMat, 'DMATS', SystemDMats); %#ok<AGROW>
0369     end
0370 
0371 varargout = {System_out};
0372 % to be added  : modyfiying the parameter lists!
0373 
0374 end

Generated on Mon 31-Mar-2008 13:54:54 by m2html © 2003