--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/test_matrix_linearize.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,75 @@
+% test for matrix/linearize
+% 
+% 15-01-2010 L Ferraioli
+%       CREATION
+% 
+% 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% $Id: mdc3_exp3_loop_v3.m,v 1.2 2009/09/24 09:48:12 luigi Exp $
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+%% Loading data
+
+fprintf('===== loading data... =====\n')
+
+% laod parnames and values
+load C:\Users\Luigi\ltp_data_analysis\MDCs\MDC3\lf_analysis\Models\parnames_10perc.mat
+load C:\Users\Luigi\ltp_data_analysis\MDCs\MDC3\lf_analysis\Models\exp3_2_10perc_nomvalues.mat
+load C:\Users\Luigi\ltp_data_analysis\MDCs\MDC3\lf_analysis\Models\exp3_2_10perc_truevalues.mat
+% load C:\Users\Luigi\ltp_data_analysis\MDCs\MDC3\lf_analysis\Models\usedparams_10perc.mat
+
+
+% set ordered used parameters
+usedparams = {'dH','dsH','dS11','dS1D','dSD1','dSDD',...
+  'dh2','dsh2','dx1','dx2'};
+
+[nonfitparnms,ia,ib] = setxor(parnames,usedparams);
+nonfitnomvals = exp3_nomvalues(ia);
+nonfittruevals = exp3_truevalues(ia);
+
+[fitparnms,ia,ib] = setxor(parnames,nonfitparnms);
+fitnomvals = exp3_nomvalues(ia);
+fittruevals = exp3_truevalues(ia);
+
+%% get non-linear response model
+% the model is non-linear in the parameter dependence
+fprintf('===== Get TF Model... =====\n')
+
+H = matrix(plist('built-in','mdc3_ifo2ifo_v2'));
+
+% substitute for non fit parameters
+% subs for non fit parameters
+for ii = 1:numel(H.objs)
+  H.objs(ii).subs(plist('Params', nonfitparnms, 'Values',nonfittruevals));
+end
+
+
+%% Linearize model
+
+Hl = linearize(H);
+
+%% Random test
+
+% get random index
+i1 = round(numel(Hl.objs)*rand(1)); % choose the parameter
+i2 = round(numel(H.objs)*rand(1)); % choose the transfer function
+
+% extract element to be tested
+dH = Hl.objs{i1}.objs(i2);
+
+% get the parameter name
+pars = sort(usedparams);
+prm = pars{i1};
+
+% differentiate
+tdH = diff(H.objs(i2),prm);
+
+% evaluate response
+dHr = eval(dH);
+tdHr = eval(tdH);
+
+% plot
+iplot(dHr./tdHr)
+
+