Home > m > timetools > statespacefunctions > ltpda_ss_modify.m

ltpda_ss_modify

PURPOSE ^

paramSubstitute enables to modify parameters and retrieve numerical

SYNOPSIS ^

function varargout = ltpda_ss_modify(varargin)

DESCRIPTION ^

 paramSubstitute 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)
 [subystem] = paramSubstitute(subystem, operations)
 [Subystem_jac, subystem] = paramSubstitute(subystem, operations)
 [Subystem_hess, Subystem_jac, subystem] = paramSubstitute(subystem, operations)
 
 INPUTS: subystem - plist describing a subsystem (refer to ltpda_ss_check for details)
 operations - plist of elements for numerical values the of parameters

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Tue 26-Feb-2008 10:52:52 by m2html © 2003