ltpda: m-toolbox/test/test_ao_polynomfit.m annotate

author	Daniele Nicolodi <nicolodi@science.unitn.it>
date	Wed, 23 Nov 2011 19:22:13 +0100
parents
children

rev	line source
0 f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	1 % TEST_AO_POLYNOMFIT tests the polynomfit method of the AO class:
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	2 % - functionality
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	3 % - against lscov
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	4 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	5 % M Hueller and D Nicolodi 07-01-10
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	6 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	7 % $Id: test_ao_polynomfit.m,v 1.2 2010/04/08 04:30:02 mauro Exp $
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	8 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	9
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	10 function test_ao_polynomfit()
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	11 %% test polynomfit vs lscov
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	12 disp(' ************** ');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	13 disp('Example with combination of noise terms');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	14 disp(' ');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	15
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	16 fs = 1;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	17 nsecs = 10;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	18 n = [0 1 -2];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	19
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	20 X = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	21 N = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', 'm'));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	22 C = [ao(1, plist('yunits', 'm')) ao(4, plist('yunits', 'm/m')) ao(2, plist('yunits', 'm/m^(-2)'))];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	23 Y = C(1) * X.^0 + C(2) * X.^1 + C(3) * X.^(-2) + N;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	24 Y.simplifyYunits;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	25
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	26 % Fit with UTN double/polynom_fit
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	27 fprintf('UTN POLYNOM_FIT:\n');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	28 try
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	29 out = polynom_fit(X.y, Y.y, n, 'nopl');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	30 catch err
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	31 disp(err.message);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	32 disp('UTN method polynom_fit.m for double not available');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	33 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	34
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	35 % Fit with Matlab lscov
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	36 fprintf('LSCOV:\n');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	37 x = X.y;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	38 y = Y.y;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	39 m = zeros(length(x), length(n));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	40 for k = 1:length(n)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	41 m(:,k) = x .^ n(k);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	42 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	43
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	44 [p, stdx, mse] = lscov(m, y);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	45 p
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	46 stdx
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	47 mse
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	48
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	49 % Fit with LTPDA ao/polynomfit
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	50 fprintf('AO/POLYNOMFIT:\n');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	51 out = polynomfit(X, Y, plist('orders', n))
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	52
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	53
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	54 %% More tests
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	55 %% Make fake AO from polyval
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	56 nsecs = 5;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	57 fs = 10;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	58 n = [0 1 -2];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	59
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	60 unit_list = unit.supportedUnits;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	61 u1 = unit(cell2mat(utils.math.randelement(unit_list,1)));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	62
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	63 pl1 = plist('nsecs', nsecs, 'fs', fs, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	64 'tsfcn', sprintf('t.^%d + t.^%d + t.^%d + randn(size(t))', n), ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	65 'xunits', 's', 'yunits', u1);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	66 a1 = ao(pl1);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	67
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	68 C = [ao(1, plist('yunits', a1.yunits)) ao(4, plist('yunits', '')) ao(2, plist('yunits', simplify(a1.yunits/(a1.yunits.^-2))))];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	69 N = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', nsecs, 'yunits', u1));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	70 a2 = C(1) * a1.^n(1) + C(2) * a1.^n(2) + C(3) * a1.^n(3) + N;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	71 a2.simplifyYunits;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	72
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	73 %% Fit with a polynomial function
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	74
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	75 p11 = polynomfit(a1, plist('orders', n ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	76 ));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	77 p12 = polynomfit(a1, plist('orders', n, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	78 'dy', 0.1));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	79 p22 = polynomfit(a1, plist('orders', n, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	80 'dx', 0.1, 'dy', 0.1, 'P0', zeros(size(n))));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	81 p13 = polynomfit(a1, plist('orders', n, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	82 'dy', 0.1*ones(size(a1.y))));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	83 p23 = polynomfit(a1, plist('orders', n, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	84 'dx', 0.1ones(size(a1.x)), 'dy', 0.1ones(size(a1.y)), 'P0', zeros(size(n))));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	85 p14 = polynomfit(a1, plist('orders', n, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	86 'dy', ao(0.1, plist('yunits', a1.yunits))));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	87 p24 = polynomfit(a1, plist('orders', n, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	88 'dx', ao(0.1, plist('yunits', a1.xunits)), 'dy', ao(0.1, plist('yunits', a1.yunits)), 'P0', zeros(size(n))));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	89 p15 = polynomfit(a1, plist('orders', n, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	90 'dy', ao(0.1*ones(size(a1.y)), plist('yunits', a1.yunits))));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	91 p25 = polynomfit(a1, plist('orders', n, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	92 'dx', ao(0.1ones(size(a1.x)), plist('yunits', a1.xunits)), 'dy', ao(0.1ones(size(a1.y)), plist('yunits', a1.yunits)), 'P0', zeros(size(n))));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	93 p26 = polynomfit(a1, plist('orders', n, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	94 'dx', 0.1, 'dy', 0.1, 'P0', ao(zeros(size(n)))));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	95 p27 = polynomfit(a1, plist('orders', n, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	96 'dx', 0.1, 'dy', 0.1, 'P0', p11));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	97
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	98 %% Compute fit: evaluating pest
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	99
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	100 b11 = p11.eval(plist('type', 'tsdata', 'XData', a1));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	101 b11.removeVal(plist('value', Inf));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	102 b12 = p12.eval(plist('type', 'tsdata', 'XData', a1));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	103 b12.removeVal(plist('value', Inf));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	104 b22 = p22.eval(a1, plist('type', 'tsdata'));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	105 b22.removeVal(plist('value', Inf));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	106 b13 = p13.eval(plist('type', 'tsdata', 'XData', a1));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	107 b13.removeVal(plist('value', Inf));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	108 b23 = p23.eval(a1, plist('type', 'tsdata'));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	109 b23.removeVal(plist('value', Inf));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	110 b14 = p14.eval(plist('type', 'tsdata', 'XData', a1));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	111 b14.removeVal(plist('value', Inf));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	112 b24 = p24.eval(a1, plist('type', 'tsdata'));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	113 b24.removeVal(plist('value', Inf));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	114 b15 = p15.eval(a1, plist('type', 'tsdata'));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	115 b15.removeVal(plist('value', Inf));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	116 b25 = p25.eval(plist('type', 'tsdata', 'XData', a1));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	117 b25.removeVal(plist('value', Inf));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	118 b26 = p26.eval(plist('type', 'tsdata', 'XData', a1.x));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	119 b26.removeVal(plist('value', Inf));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	120 b27 = p27.eval(plist('type', 'tsdata', 'XData', a1.x));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	121 b27.removeVal(plist('value', Inf));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	122
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	123 %% Plot fit
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	124 iplot(a1, b11, b12, b13, b14, b15, b26, b27, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	125 plist('LineStyles', {'', '--'}));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	126
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	127 %% Remove linear trend
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	128 c = a1.removeVal(plist('value', Inf))-b11;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	129 iplot(c)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	130
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	131 plot(c.hist)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	132
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	133 %% Reproduce from history
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	134 disp('Try rebuilding')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	135 a_out = rebuild(c);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	136 iplot(a_out)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	137 plot(a_out.hist)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	138
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	139 %% Fit with a polynomial function
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	140
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	141 p11 = polynomfit(a1, a2, plist('orders', n ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	142 ));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	143 p12 = polynomfit(a1, a2, plist('orders', n, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	144 'dy', 0.1));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	145 p22 = polynomfit(a1, a2, plist('orders', n, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	146 'dx', 0.1, 'dy', 0.1, 'P0', zeros(size(n))));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	147 p13 = polynomfit(a1, a2, plist('orders', n, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	148 'dy', 0.1*ones(size(a1.x))));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	149 p23 = polynomfit(a1, a2, plist('orders', n, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	150 'dx', 0.1ones(size(a1.x)), 'dy', 0.1ones(size(a1.x)), 'P0', zeros(size(n))));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	151 p14 = polynomfit(a1, a2, plist('orders', n, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	152 'dy', ao(0.1, plist('yunits', a2.yunits))));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	153 p24 = polynomfit(a1, a2, plist('orders', n, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	154 'dx', ao(0.1, plist('yunits', a1.yunits)), 'dy', ao(0.1, plist('yunits', a2.yunits)), 'P0', zeros(size(n))));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	155 p15 = polynomfit(a1, a2, plist('orders', n, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	156 'dy', ao(0.1*ones(size(a2.y)), plist('yunits', a2.yunits))));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	157 p25 = polynomfit(a1, a2, plist('orders', n, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	158 'dx', ao(0.1ones(size(a1.y)), plist('yunits', a1.yunits)), 'dy', ao(0.1ones(size(a2.y)), plist('yunits', a2.yunits)), 'P0', zeros(size(n))));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	159 p26 = polynomfit(a1, a2, plist('orders', n, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	160 'dx', 0.1, 'dy', 0.1, 'P0', ao(zeros(size(n)))));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	161 p27 = polynomfit(a1, a2, plist('orders', n, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	162 'dx', 0.1, 'dy', 0.1, 'P0', p11));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	163
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	164 %% Compute fit: evaluating pest
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	165
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	166 b11 = p11.eval(plist('type', 'xydata', 'XData', a1.y));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	167 b12 = p12.eval(plist('type', 'xydata', 'XData', a1));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	168 b22 = p22.eval(a1, plist('type', 'xydata'));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	169 b13 = p13.eval(plist('type', 'xydata', 'XData', a1.y));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	170 b23 = p23.eval(a1, plist('type', 'xydata'));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	171 b14 = p14.eval(plist('type', 'xydata', 'XData', a1.y));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	172 b24 = p24.eval(a1, plist('type', 'xydata'));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	173 b15 = p15.eval(plist('type', 'xydata', 'XData', a1.y));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	174 b25 = p25.eval(plist('type', 'xydata', 'XData', a1.y));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	175 b26 = p26.eval(plist('type', 'xydata', 'XData', a1.y));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	176 b27 = p27.eval(plist('type', 'xydata', 'XData', a1.y));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	177
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	178 % Build reference object
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	179 a12 = ao(plist('xvals', a1.y, 'yvals', a2.y, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	180 'xunits', a1.yunits, 'yunits', a2.yunits));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	181 %% Plot fit
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	182 iplot(a12, b11, b12, b23, b14, b25, b26, b27, ...
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	183 plist('LineStyles', {'', '--'}));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	184
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	185 %% Remove linear trend
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	186 c = a12-b27;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	187 iplot(c)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	188
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	189 plot(c.hist)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	190
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	191 %% Reproduce from history
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	192 disp('Try rebuilding')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	193 a_out = rebuild(c);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	194 iplot(a_out)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	195 plot(a_out.hist)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	196
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	197 % end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	198 % END
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	199 end

0

f0afece42f48 Import.