Mercurial > forge
annotate extra/NaN/inst/nantest.m @ 12656:ddbd9137910a octave-forge
fix nantest to avoid crashing Octave 4.0.0 on Windows
| author | schloegl |
|---|---|
| date | Thu, 02 Jul 2015 11:14:16 +0000 |
| parents | 9a277c75a100 |
| children | aca981120490 |
| rev | line source |
|---|---|
| 2414 | 1 % NANTEST checks several mathematical operations and a few |
| 2 % statistical functions for their correctness related to NaN's. | |
| 3 % e.g. it checks norminv, normcdf, normpdf, sort, matrix division and multiplication. | |
| 4 % | |
| 5 % | |
| 6 % see also: NANINSTTEST | |
| 7 % | |
| 8 % REFERENCE(S): | |
| 9 % [1] W. Kahan (1996) Lecture notes on the Status of "IEEE Standard 754 for | |
| 10 % Binary Floating-point Arithmetic. | |
| 11 % | |
| 12 | |
| 13 | |
| 14 % This program is free software; you can redistribute it and/or modify | |
| 15 % it under the terms of the GNU General Public License as published by | |
| 16 % the Free Software Foundation; either version 2 of the License, or | |
| 17 % (at your option) any later version. | |
| 18 % | |
| 19 % This program is distributed in the hope that it will be useful, | |
| 20 % but WITHOUT ANY WARRANTY; without even the implied warranty of | |
| 21 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
| 22 % GNU General Public License for more details. | |
| 23 % | |
| 24 % You should have received a copy of the GNU General Public License | |
| 4404 | 25 % along with this program; If not, see <http://www.gnu.org/licenses/>. |
| 2414 | 26 |
| 27 % $Id$ | |
| 8037 | 28 % Copyright (C) 2000-2004,2009 by Alois Schloegl <alois.schloegl@gmail.com> |
| 2414 | 29 % This script is part of the NaN-toolbox |
| 7889 | 30 % http://pub.ist.ac.at/~schloegl/matlab/NaN/ |
| 2414 | 31 |
| 4352 | 32 %FLAG_WARNING = warning; |
| 33 %warning('off'); | |
| 2414 | 34 |
| 6534 | 35 try |
| 2414 | 36 x = randn([3,4,5]); |
| 37 x(~isnan(x)) = 0; | |
| 6534 | 38 catch |
| 2414 | 39 fprintf(1,'WARNING: NANTEST fails for 3-DIM matrices. \n'); |
| 40 end; | |
| 6534 | 41 try |
| 2414 | 42 [s,n] = sumskipnan([nan,1,4,5]); |
| 6534 | 43 catch |
| 2414 | 44 fprintf(1,'WARNING: SUMSKIPNAN is not avaible. \n'); |
| 45 end; | |
| 46 | |
| 47 % check NORMPDF, NORMCDF, NORMINV | |
| 48 x = [-inf,-2,-1,-.5,0,.5,1,2,3,inf,nan]'; | |
| 6534 | 49 if exist('normpdf','file')==2, |
| 2414 | 50 q(1) = sum(isnan(normpdf(x,2,0)))>sum(isnan(x)); |
| 51 if q(1), | |
| 52 fprintf(1,'NORMPDF cannot handle v=0.\n'); | |
| 53 fprintf(1,'-> NORMPDF should be replaced\n'); | |
| 54 end; | |
| 55 end; | |
| 56 | |
| 6534 | 57 if exist('normcdf','file')==2, |
| 2414 | 58 q(2) = sum(isnan(normcdf(x,2,0)))>sum(isnan(x)); |
| 59 if q(2), | |
| 60 fprintf(1,'NORMCDF cannot handle v=0.\n'); | |
| 61 fprintf(1,'-> NORMCDF should be replaced\n'); | |
| 62 end; | |
| 63 end; | |
| 64 | |
|
12548
9a277c75a100
fix compatibility issue with octave 3.8.2 in debian/jessie
schloegl
parents:
12515
diff
changeset
|
65 if ~(any(exist('erfinv') == [2,5])) |
| 12389 | 66 fprintf(1,'ERFINV is not available\n'); |
| 67 | |
| 68 elseif exist('norminv','file')==2, | |
| 2414 | 69 p = [-inf,-.2,0,.2,.5,1,2,inf,nan]; |
| 70 q(3) = sum(~isnan(norminv(p,2,0)))<4; | |
| 71 if q(3), | |
| 72 fprintf(1,'NORMINV cannot handle correctly v=0.\n'); | |
| 73 fprintf(1,'-> NORMINV should be replaced\n'); | |
| 74 end; | |
| 75 q(4) = ~isnan(norminv(0,NaN,0)); | |
| 6535 | 76 q(5) = any(norminv(0.5,[1 2 3],0)~=(1:3)); |
| 2414 | 77 end; |
| 78 | |
| 6534 | 79 if exist('tpdf','file')==2, |
| 2414 | 80 q(6) = ~isnan(tpdf(nan,4)); |
| 81 if q(6), | |
| 82 fprintf(1,'TPDF(NaN,4) does not return NaN\n'); | |
| 83 fprintf(1,'-> TPDF should be replaced\n'); | |
| 84 end; | |
| 85 end; | |
| 86 | |
| 6534 | 87 if exist('tcdf','file')==2, |
| 88 try | |
| 2414 | 89 q(7) = ~isnan(tcdf(nan,4)); |
| 6534 | 90 catch |
| 2414 | 91 q(7) = 1; |
| 92 end; | |
| 93 if q(7), | |
| 94 fprintf(1,'TCDF(NaN,4) does not return NaN\n'); | |
| 95 fprintf(1,'-> TCDF should be replaced\n'); | |
| 96 end; | |
| 97 end; | |
| 98 | |
| 6534 | 99 if exist('tinv','file')==2, |
| 100 try | |
| 2414 | 101 q(8) = ~isnan(tinv(nan,4)); |
| 6534 | 102 catch |
| 2414 | 103 q(8) = 1; |
| 104 end; | |
| 105 if q(8), | |
| 106 fprintf(1,'TINV(NaN,4) does not return NaN\n'); | |
| 107 fprintf(1,'-> TINV should be replaced\n'); | |
| 108 end; | |
| 109 end; | |
| 110 | |
| 6534 | 111 q(9) = isreal(double(2+3i)); |
|
5741
f0955d62224a
faster exist() improves performance (especially in covm.m)
schloegl
parents:
5425
diff
changeset
|
112 if q(9) |
|
f0955d62224a
faster exist() improves performance (especially in covm.m)
schloegl
parents:
5425
diff
changeset
|
113 printf('DOUBLE rejects imaginary part\n-> this can affect SUMSKIPNAN\n'); |
|
f0955d62224a
faster exist() improves performance (especially in covm.m)
schloegl
parents:
5425
diff
changeset
|
114 end; |
|
f0955d62224a
faster exist() improves performance (especially in covm.m)
schloegl
parents:
5425
diff
changeset
|
115 |
| 6519 | 116 try |
| 117 x = reshape(1:6,3,2); | |
| 118 [cc,nn] = covm(x+i*x,'e'); | |
| 119 q(10) = 0; | |
| 120 catch | |
| 121 q(10) = 1; | |
| 122 end; | |
|
5741
f0955d62224a
faster exist() improves performance (especially in covm.m)
schloegl
parents:
5425
diff
changeset
|
123 |
|
f0955d62224a
faster exist() improves performance (especially in covm.m)
schloegl
parents:
5425
diff
changeset
|
124 if 0, |
| 2414 | 125 %%%%% MOD |
| 126 if exist('mod')>1, | |
| 127 if (mod(5,0))~=0, | |
| 128 fprintf(1,'WARNING: MOD(x,0) does not return 0.\n'); | |
| 129 end; | |
| 130 if isnan(mod(5,0)), | |
| 131 fprintf(1,'WARNING: MOD(x,0) returns NaN.\n'); | |
| 132 end; | |
| 133 if isnan(mod(5,inf)), | |
| 134 fprintf(1,'WARNING: MOD(x,INF) returns NaN.\n'); | |
| 135 end; | |
| 136 end; | |
| 137 %%%%% REM | |
| 138 if exist('rem')>1, | |
| 139 if (rem(5,0))~=0, | |
| 140 fprintf(1,'WARNING: REM(x,0) does not return 0.\n'); | |
| 141 end; | |
| 142 if isnan(rem(5,0)), | |
| 143 fprintf(1,'WARNING: REM(x,0) returns NaN.\n'); | |
| 144 end; | |
| 145 if isnan(mod(5,inf)), | |
| 146 fprintf(1,'WARNING: REM(x,INF) returns NaN.\n'); | |
| 147 end; | |
| 148 end; | |
|
5741
f0955d62224a
faster exist() improves performance (especially in covm.m)
schloegl
parents:
5425
diff
changeset
|
149 end; |
| 2414 | 150 |
| 151 %%%%% NANSUM(NAN) - this test addresses a problem in Matlab 5.3, 6.1 & 6.5 | |
| 6534 | 152 if exist('nansum','file'), |
| 2414 | 153 if isnan(nansum(nan)), |
| 154 fprintf(1,'Warning: NANSUM(NaN) returns NaN instead of 0\n'); | |
| 155 fprintf(1,'-> NANSUM should be replaced\n'); | |
| 156 end; | |
| 157 end; | |
| 158 %%%%% NANSUM(NAN) - this test addresses a problem in Matlab 5.3, 6.1 & 6.5 | |
| 6534 | 159 if exist('nanstd','file'), |
| 2414 | 160 if ~isnan(nanstd(0)), |
| 161 fprintf(1,'Warning: NANSTD(x) with isscalar(x) returns 0 instead of NaN\n'); | |
| 162 fprintf(1,'-> NANSTD should be replaced\n'); | |
| 163 end; | |
| 164 end; | |
| 165 %%%%% GEOMEAN - this test addresses a problem in Octave | |
| 6534 | 166 if exist('geomean','file'), |
| 167 if isnan(geomean((0:3)')), | |
| 2414 | 168 fprintf(1,'Warning: GEOMEAN([0,1,2,3]) NaN instead of 0\n'); |
| 169 fprintf(1,'-> GEOMEAN should be replaced\n'); | |
| 170 end; | |
| 171 end; | |
| 172 %%%%% HARMMEAN - this test addresses a problem in Octave | |
| 6534 | 173 if exist('harmmean','file'), |
| 2414 | 174 if isnan(harmmean(0:3)), |
| 175 fprintf(1,'Warning: HARMMEAN([0,1,2,3]) NaN instead of 0\n'); | |
| 176 fprintf(1,'-> HARMMEAN should be replaced\n'); | |
| 177 end; | |
| 178 end; | |
| 179 %%%%% BITAND - this test addresses a problem in Octave | |
| 180 if exist('bitand')>1, | |
| 181 if isnan(bitand(2^33-1,13)), | |
| 182 fprintf(1,'BITAND can return NaN. \n'); | |
| 183 end; | |
| 184 end; | |
| 185 %%%%% BITSHIFT - this test addresses a problem in Octave | |
| 6534 | 186 if exist('bitshift','file'), |
| 2414 | 187 if isnan(bitshift(5,30,32)), |
| 188 fprintf(1,'BITSHIFT can return NaN.\n'); | |
| 189 end; | |
| 190 end; | |
| 4352 | 191 %%%%% ALL - this test addresses a problem in some old Octave and FreeMat v3.5 |
| 2414 | 192 if any(NaN)==1, |
| 4352 | 193 fprintf(1,'WARNING: ANY(NaN) returns 1 instead of 0\n'); |
| 2414 | 194 end; |
| 195 if any([])==1, | |
| 4352 | 196 fprintf(1,'WARNING: ANY([]) returns 1 instead of 0\n'); |
| 2414 | 197 end; |
| 4352 | 198 %%%%% ALL - this test addresses a problem in some old Octave and FreeMat v3.5 |
| 2414 | 199 if all(NaN)==0, |
| 4352 | 200 fprintf(1,'WARNING: ALL(NaN) returns 0 instead of 1\n'); |
| 2414 | 201 end; |
| 202 if all([])==0, | |
| 4352 | 203 fprintf(1,'WARNING: ALL([]) returns 0 instead of 1\n'); |
| 2414 | 204 end; |
| 205 | |
| 12389 | 206 %%%%% SORT - this was once a problem in Octave Version < 2.1.36, and still is in FreeMat 4.0 %%%% |
| 2414 | 207 if ~all(isnan(sort([3,4,NaN,3,4,NaN]))==[0,0,0,0,1,1]), |
| 12389 | 208 warning('Warning: SORT does not properly handle NaN.'); |
| 2414 | 209 end; |
| 210 | |
| 211 %%%%% commutativity of 0*NaN %%% This test adresses a problem in Octave | |
| 212 x=[-2:2;4:8]'; | |
| 213 y=x;y(2,1)=nan;y(4,2)=nan; | |
| 214 B=[1,0,2;0,3,1]; | |
| 215 if ~all(all(isnan(y*B)==isnan(B'*y')')), | |
| 216 fprintf(2,'WARNING: 0*NaN within matrix multiplication is not commutative\n'); | |
| 217 end; | |
| 218 | |
| 219 % from Kahan (1996) | |
| 220 tmp = (0-3*i)/inf; | |
| 221 if isnan(tmp) | |
| 222 fprintf(2,'WARNING: (0-3*i)/inf results in NaN instead of 0.\n'); | |
| 223 end; | |
| 224 | |
| 225 %(roots([5,0,0])-[0;0]) | |
| 226 %(roots([2,-10,12])-[3;2]) | |
| 227 %(roots([2e-37,-2,2])-[1e37;1]) | |
| 228 %%%%% check nan/nan %% this test addresses a problem in Matlab 5.3, 6.1 & 6.5 | |
| 229 p = 4; | |
|
12656
ddbd9137910a
fix nantest to avoid crashing Octave 4.0.0 on Windows
schloegl
parents:
12548
diff
changeset
|
230 tmp1 = repmat(nan, 4); |
|
ddbd9137910a
fix nantest to avoid crashing Octave 4.0.0 on Windows
schloegl
parents:
12548
diff
changeset
|
231 tmp2 = repmat(nan, 4); |
|
ddbd9137910a
fix nantest to avoid crashing Octave 4.0.0 on Windows
schloegl
parents:
12548
diff
changeset
|
232 if ~ispc |
|
12548
9a277c75a100
fix compatibility issue with octave 3.8.2 in debian/jessie
schloegl
parents:
12515
diff
changeset
|
233 try |
|
9a277c75a100
fix compatibility issue with octave 3.8.2 in debian/jessie
schloegl
parents:
12515
diff
changeset
|
234 tmp1 = repmat(nan,p)/repmat(nan,p); |
|
9a277c75a100
fix compatibility issue with octave 3.8.2 in debian/jessie
schloegl
parents:
12515
diff
changeset
|
235 catch % exception error in Octave 3.8.2 of debian wheezy |
|
12656
ddbd9137910a
fix nantest to avoid crashing Octave 4.0.0 on Windows
schloegl
parents:
12548
diff
changeset
|
236 fprintf(2,'repmat(nan,4)/repmat(nan,4) fails with an exception\n'); |
|
12548
9a277c75a100
fix compatibility issue with octave 3.8.2 in debian/jessie
schloegl
parents:
12515
diff
changeset
|
237 end; |
|
9a277c75a100
fix compatibility issue with octave 3.8.2 in debian/jessie
schloegl
parents:
12515
diff
changeset
|
238 try |
|
9a277c75a100
fix compatibility issue with octave 3.8.2 in debian/jessie
schloegl
parents:
12515
diff
changeset
|
239 tmp2 = repmat(nan,p)\repmat(nan,p); |
|
9a277c75a100
fix compatibility issue with octave 3.8.2 in debian/jessie
schloegl
parents:
12515
diff
changeset
|
240 catch % exception error in Octave 3.8.2 of debian wheezy |
|
12656
ddbd9137910a
fix nantest to avoid crashing Octave 4.0.0 on Windows
schloegl
parents:
12548
diff
changeset
|
241 fprintf(2,'repmat(nan,4)\repmat(nan,4) fails with an exception\n'); |
|
12548
9a277c75a100
fix compatibility issue with octave 3.8.2 in debian/jessie
schloegl
parents:
12515
diff
changeset
|
242 end |
|
12656
ddbd9137910a
fix nantest to avoid crashing Octave 4.0.0 on Windows
schloegl
parents:
12548
diff
changeset
|
243 end; |
| 2414 | 244 tmp3 = repmat(0,p)/repmat(0,p); |
| 245 tmp4 = repmat(0,p)\repmat(0,p); | |
| 246 tmp5 = repmat(0,p)*repmat(inf,p); | |
| 247 tmp6 = repmat(inf,p)*repmat(0,p); | |
| 248 x = randn(100,1)*ones(1,p); y=x'*x; | |
| 249 tmp7 = y/y; | |
| 250 tmp8 = y\y; | |
| 251 | |
| 252 if ~all(isnan(tmp1(:))), | |
| 4352 | 253 fprintf(1,'WARNING: matrix division NaN/NaN does not result in NaN\n'); |
| 2414 | 254 end; |
| 255 if ~all(isnan(tmp2(:))), | |
| 4352 | 256 fprintf(1,'WARNING: matrix division NaN\\NaN does not result in NaN\n'); |
| 2414 | 257 end; |
| 258 if ~all(isnan(tmp3(:))), | |
| 259 fprintf(2,'WARNING: matrix division 0/0 does not result in NaN\n'); | |
| 260 end; | |
| 261 if ~all(isnan(tmp4(:))), | |
| 262 fprintf(2,'WARNING: matrix division 0\\0 does not result in NaN\n'); | |
| 263 end; | |
| 264 if ~all(isnan(tmp5(:))), | |
| 265 fprintf(2,'WARNING: matrix multiplication 0*inf does not result in NaN\n'); | |
| 266 end; | |
| 267 if ~all(isnan(tmp6(:))), | |
| 268 fprintf(2,'WARNING: matrix multiplication inf*0 does not result in NaN\n'); | |
| 269 end; | |
| 270 if any(any(tmp7==inf)); | |
| 271 fprintf(2,'WARNING: right division of two singulare matrices return INF\n'); | |
| 272 end; | |
| 273 if any(any(tmp8==inf)); | |
| 274 fprintf(2,'WARNING: left division of two singulare matrices return INF\n'); | |
| 275 end; | |
| 276 | |
| 277 tmp = [tmp1;tmp2;tmp3;tmp4;tmp5;tmp6;tmp7;tmp8]; | |
| 278 | |
|
5741
f0955d62224a
faster exist() improves performance (especially in covm.m)
schloegl
parents:
5425
diff
changeset
|
279 |
|
f0955d62224a
faster exist() improves performance (especially in covm.m)
schloegl
parents:
5425
diff
changeset
|
280 |
| 4352 | 281 %warning(FLAG_WARNING); |
| 3761 | 282 |
| 283 | |
| 284 %%%%% QUANTILE TEST | |
|
12656
ddbd9137910a
fix nantest to avoid crashing Octave 4.0.0 on Windows
schloegl
parents:
12548
diff
changeset
|
285 d = [1 1 2 2 4 4 10 700]'; |
| 3761 | 286 q = [-1,0,.05,.1,.25,.49,.5,.51,.75,.8, .999999,1,2]; |
| 287 r = [ NaN, 1, 1, 1, 1.5, 2, 3, 4, 7, 10, 700, 700, NaN]; | |
|
12656
ddbd9137910a
fix nantest to avoid crashing Octave 4.0.0 on Windows
schloegl
parents:
12548
diff
changeset
|
288 if any( quantile(d, q)' - r>0) |
| 3761 | 289 fprintf(1,'Quantile(1): failed\n'); |
| 290 else | |
| 291 fprintf(1,'Quantile(1): OK\n'); | |
| 292 end; | |
| 293 if exist('histo3','file') | |
|
12656
ddbd9137910a
fix nantest to avoid crashing Octave 4.0.0 on Windows
schloegl
parents:
12548
diff
changeset
|
294 H = histo3(d); |
| 3761 | 295 else |
| 296 H.X = [1;2;4;10;700]; | |
| 297 H.H = [2;2;2;1;1]; | |
| 5425 | 298 H.datatype = 'HISTOGRAM'; |
| 3761 | 299 end; |
| 300 if any( quantile(H, q)' - r>0) | |
| 301 fprintf(1,'Quantile(2): failed\n'); | |
| 302 else | |
| 303 fprintf(1,'Quantile(2): OK\n'); | |
| 304 end; | |
| 305 |