Path: utzoo!utgpu!news-server.csri.toronto.edu!cs.utexas.edu!uunet!sparky!kent From: arg@ccvr1.cc.ncsu.edu (Ron Gallant) Newsgroups: comp.sources.misc Subject: v16i066: nlmdl - Estimate nonlinear statistical models, Part04/06 Message-ID: <1991Jan12.050650.8742@sparky.IMD.Sterling.COM> Date: 12 Jan 91 05:06:50 GMT Sender: kent@sparky.IMD.Sterling.COM (Kent Landfield) Organization: Sterling Software, IMD Lines: 1355 Approved: kent@sparky.imd.sterling.com X-Checksum-Snefru: 3a73dba5 e051be8c 6dc80daa cf231341 Submitted-by: arg@ccvr1.cc.ncsu.edu (Ron Gallant) Posting-number: Volume 16, Issue 66 Archive-name: nlmdl/part04 # This is a shell archive. Remove anything before this line, then unpack # it by saving it into a file and typing "sh file". To overwrite existing # files, type "sh file -c". You can also feed this as standard input via # unshar, or by typing "sh 'ch5eg1/electa.dat' <<'END_OF_FILE' X 1 1 0.056731 0.280382 0.662888 0.46931 X 2 1 0.103444 0.252128 0.644427 0.79539 X 3 1 0.158353 0.270089 0.571558 0.45756 X 4 1 0.108075 0.305072 0.586853 0.94713 X 5 1 0.083921 0.211656 0.704423 1.22054 X 6 1 0.112165 0.290532 0.597302 0.93181 X 7 1 0.071274 0.240518 0.688208 1.79152 X 8 1 0.076510 0.210503 0.712987 0.51442 X 9 1 0.066173 0.202999 0.730828 0.78407 X 10 1 0.094836 0.270281 0.634883 1.01354 X 11 1 0.078501 0.293953 0.627546 0.83854 X 12 1 0.059530 0.228752 0.711718 1.53957 X 13 1 0.208982 0.328053 0.462965 1.06694 X 14 1 0.083702 0.297272 0.619027 0.82437 X 15 1 0.138705 0.358329 0.502966 0.80712 X 16 1 0.111378 0.322564 0.566058 0.53169 X 17 1 0.092919 0.259633 0.647448 0.85439 X 18 1 0.039353 0.158205 0.802442 1.93326 X 19 1 0.066577 0.247454 0.685970 1.37160 X 20 2 0.102844 0.244335 0.652821 0.92766 X 21 2 0.125485 0.230305 0.644210 1.80934 X 22 2 0.154316 0.235135 0.610549 2.41501 X 23 2 0.165714 0.276980 0.557305 0.84658 X 24 2 0.145370 0.173112 0.681518 1.60788 X 25 2 0.184467 0.268865 0.546668 0.73838 X 26 2 0.162269 0.280939 0.556792 0.81116 X 27 2 0.112016 0.220850 0.667133 2.01503 X 28 2 0.226863 0.257833 0.515304 2.32035 X 29 2 0.118028 0.219830 0.662142 2.40172 X 30 2 0.137761 0.345117 0.517122 0.57141 X 31 2 0.079115 0.257319 0.663566 0.94474 X 32 2 0.185022 0.265051 0.549928 1.63778 X 33 2 0.144524 0.276133 0.579343 0.75816 X 34 2 0.201734 0.241966 0.556300 1.00136 X 35 2 0.094890 0.227651 0.677459 1.11384 X 36 2 0.102843 0.264515 0.632642 1.07185 X 37 2 0.107760 0.214232 0.678009 1.53659 X 38 2 0.156552 0.236422 0.607026 0.24099 X 39 2 0.088431 0.222746 0.688822 0.58066 X 40 2 0.146236 0.301884 0.551880 2.52983 X 41 3 0.080802 0.199005 0.720192 1.14741 X 42 3 0.100711 0.387758 0.511531 0.97934 X 43 3 0.073483 0.335280 0.591237 1.09361 X 44 3 0.059455 0.259823 0.680722 2.19468 X 45 3 0.076195 0.378371 0.545434 1.98221 X 46 3 0.076926 0.325032 0.598042 1.78194 X 47 3 0.086052 0.339653 0.574295 3.24274 X 48 3 0.069359 0.278369 0.652272 0.47593 X 49 3 0.071265 0.273866 0.654869 1.38369 X 50 3 0.100562 0.306247 0.593191 1.57831 X 51 3 0.050203 0.294285 0.655513 2.16900 X 52 3 0.059627 0.311932 0.628442 2.11575 X 53 3 0.081433 0.328604 0.589962 0.35681 X 54 3 0.075762 0.285972 0.638265 1.55275 X 55 3 0.042910 0.372337 0.584754 1.06305 X 56 3 0.086846 0.340184 0.572970 4.02013 X 57 3 0.102537 0.335535 0.561928 0.60712 X 58 3 0.068766 0.310782 0.620452 1.15334 X 59 3 0.058405 0.307111 0.634485 2.43797 X 60 4 0.055227 0.300839 0.643934 0.10082 X 61 4 0.107435 0.273937 0.618628 0.69302 X 62 4 0.105958 0.291205 0.602837 1.12592 X 63 4 0.132278 0.279429 0.588293 1.84425 X 64 4 0.094195 0.328866 0.576940 1.57972 X 65 4 0.115259 0.401079 0.483663 1.27034 X 66 4 0.150229 0.317866 0.531905 0.56330 X 67 4 0.168780 0.307669 0.523551 3.43139 X 68 4 0.118222 0.318080 0.563698 1.00979 X 69 4 0.103394 0.307671 0.588936 2.08458 X 70 4 0.124007 0.362115 0.513879 1.30410 X 71 4 0.197987 0.280130 0.521884 3.48146 X 72 4 0.108083 0.337004 0.554913 0.53206 X 73 5 0.088798 0.232568 0.678634 3.28987 X 74 5 0.100508 0.272139 0.627353 0.32678 X 75 5 0.127303 0.298519 0.574178 0.52452 X 76 5 0.109718 0.228172 0.662109 0.36622 X 77 5 0.130080 0.231037 0.638883 0.63788 X 78 5 0.148562 0.323579 0.527859 1.42239 X 79 5 0.106306 0.252137 0.641556 0.93535 X 80 5 0.080877 0.214172 0.704951 1.26243 X 81 5 0.081810 0.135665 0.782525 1.51472 X 82 5 0.131749 0.278338 0.589913 2.07858 X 83 5 0.059180 0.254533 0.686287 1.60681 X 84 5 0.078620 0.267252 0.654128 1.54706 X 85 5 0.090220 0.293831 0.615949 2.61162 X 86 5 0.086916 0.193967 0.719117 2.96418 X 87 5 0.132383 0.230489 0.637127 0.26912 X 88 5 0.085560 0.252321 0.662120 0.42554 X 89 5 0.071368 0.276238 0.652393 1.01926 X 90 5 0.061196 0.245025 0.693780 1.53807 X 91 5 0.086608 0.233981 0.679411 0.75711 X 92 5 0.105628 0.305471 0.588901 0.83647 X 93 5 0.078158 0.202536 0.719307 1.92096 X 94 5 0.048632 0.216807 0.734560 1.57795 X 95 5 0.094527 0.224344 0.681128 0.83216 X 96 5 0.092809 0.209154 0.698037 1.39364 X 97 5 0.035751 0.166231 0.798018 1.72697 X 98 5 0.065205 0.205058 0.729736 2.04120 X 99 5 0.092561 0.193848 0.713591 2.04708 X 100 5 0.063119 0.234114 0.702767 3.43969 X 101 5 0.091186 0.224488 0.684326 2.66918 X 102 5 0.047291 0.262623 0.690086 2.71072 X 103 5 0.081575 0.206400 0.712025 3.36803 X 104 5 0.108165 0.243650 0.648185 0.65682 X 105 5 0.079534 0.320450 0.600017 0.95523 X 106 5 0.084828 0.247189 0.667984 0.61441 X 107 5 0.063747 0.210343 0.725910 1.85034 X 108 5 0.081108 0.249960 0.668932 2.11274 X 109 5 0.089942 0.206601 0.703457 1.54120 X 110 5 0.046717 0.224784 0.728499 3.54351 X 111 5 0.114925 0.272279 0.612796 2.61769 X 112 5 0.115055 0.264415 0.620530 3.00236 X 113 5 0.081511 0.223870 0.694618 1.74166 X 114 5 0.109658 0.343593 0.546750 1.17640 X 115 5 0.114263 0.304761 0.580976 0.74566 X 116 5 0.115089 0.226412 0.658499 1.30392 X 117 5 0.040622 0.198986 0.760392 2.13339 X 118 5 0.073245 0.238522 0.688234 2.83039 X 119 5 0.087954 0.287450 0.624596 1.62179 X 120 5 0.091967 0.206131 0.701902 2.18534 X 121 5 0.142746 0.302939 0.554315 0.26503 X 122 5 0.117972 0.253811 0.628217 0.05082 X 123 5 0.071573 0.248324 0.680103 0.42740 X 124 5 0.073628 0.290586 0.635786 0.47979 X 125 5 0.121075 0.350781 0.528145 0.59551 X 126 5 0.077335 0.339358 0.583307 0.47506 X 127 5 0.074766 0.167202 0.758032 2.11867 X 128 5 0.208580 0.331363 0.460058 1.13621 X 129 5 0.080195 0.210619 0.709185 2.61204 X 130 5 0.066156 0.204118 0.729726 1.45227 X 131 5 0.112282 0.252638 0.635080 0.79071 X 132 5 0.041310 0.093106 0.865584 1.30697 X 133 5 0.102675 0.297009 0.600316 0.93691 X 134 5 0.102902 0.270832 0.626266 0.98718 X 135 5 0.118932 0.250104 0.630964 1.40085 X 136 5 0.139760 0.322394 0.537846 1.78710 X 137 5 0.121616 0.214626 0.663758 8.46237 X 138 5 0.065701 0.263818 0.670481 1.58663 X 139 5 0.034029 0.175181 0.790790 2.62535 X 140 5 0.074476 0.194744 0.730780 4.29430 X 141 5 0.059568 0.229705 0.710727 0.65404 X 142 5 0.088128 0.295546 0.616326 0.41292 X 143 5 0.075522 0.213622 0.710856 2.02370 X 144 5 0.057089 0.195720 0.747190 1.76998 X 145 5 0.096331 0.301692 0.601977 0.99891 X 146 5 0.120824 0.250280 0.628896 0.27942 X 147 6 0.034529 0.193456 0.772015 0.91673 X 148 6 0.026971 0.180848 0.792181 1.15617 X 149 6 0.045271 0.141894 0.812835 1.57107 X 150 6 0.067708 0.219302 0.712990 1.24515 X 151 6 0.079335 0.230693 0.689972 1.70748 X 152 6 0.022703 0.178896 0.798401 1.79959 X 153 6 0.043053 0.157142 0.799805 4.61665 X 154 6 0.057157 0.245931 0.696912 0.59504 X 155 6 0.063229 0.136192 0.800579 1.42499 X 156 6 0.076873 0.214209 0.708918 1.34371 X 157 6 0.027353 0.124894 0.847753 2.74908 X 158 6 0.067823 0.146994 0.785183 1.84628 X 159 6 0.056388 0.189185 0.754428 3.82472 X 160 6 0.036841 0.194994 0.768165 1.18199 X 161 6 0.059160 0.138681 0.802158 2.07338 X 162 6 0.051980 0.215700 0.732320 0.80376 X 163 6 0.027300 0.145072 0.827628 1.52316 X 164 6 0.014790 0.179619 0.805591 3.17526 X 165 6 0.047865 0.167561 0.784574 3.30794 X 166 6 0.115629 0.231381 0.652990 0.72456 X 167 7 0.104970 0.147525 0.747505 0.50274 X 168 7 0.119254 0.187409 0.693337 1.22571 X 169 7 0.042564 0.112839 0.844596 2.13534 X 170 7 0.096756 0.150178 0.753066 5.56011 X 171 7 0.063013 0.168422 0.768565 3.11725 X 172 7 0.080060 0.143934 0.776006 0.99796 X 173 7 0.097493 0.173391 0.729116 0.67859 X 174 7 0.102526 0.220954 0.676520 0.79027 X 175 7 0.085538 0.195686 0.718776 2.24498 X 176 7 0.068733 0.166248 0.765019 2.01993 X 177 7 0.094915 0.140119 0.764966 4.07330 X 178 7 0.076163 0.132046 0.791792 3.66432 X 179 7 0.099943 0.176885 0.723172 0.40768 X 180 7 0.081494 0.175082 0.743425 1.09065 X 181 7 0.196026 0.299348 0.504626 1.35008 X 182 7 0.093173 0.235816 0.671011 1.06138 X 183 7 0.172293 0.173032 0.654675 0.99219 X 184 7 0.067736 0.159600 0.772663 3.69199 X 185 7 0.102033 0.171697 0.726271 2.36676 X 186 7 0.067977 0.151109 0.780914 1.84563 X 187 8 0.071073 0.238985 0.689942 0.18316 X 188 8 0.049453 0.286788 0.663759 2.23986 X 189 8 0.062748 0.255129 0.682123 3.48084 X 190 8 0.032376 0.154905 0.812719 7.26135 X 191 8 0.055055 0.225296 0.719648 1.68814 X 192 8 0.037829 0.179051 0.783120 1.13804 X 193 8 0.020102 0.172396 0.807502 1.40894 X 194 8 0.021917 0.149092 0.828992 3.47472 X 195 8 0.047590 0.174735 0.777675 3.37689 X 196 8 0.063446 0.235823 0.700731 3.14810 X 197 8 0.034719 0.159398 0.805883 3.21710 X 198 8 0.055428 0.200488 0.744084 1.13941 X 199 8 0.058074 0.254823 0.687103 2.55414 X 200 8 0.060719 0.209763 0.729518 0.29071 X 201 8 0.045681 0.206177 0.748142 1.21336 X 202 8 0.040151 0.263161 0.696688 1.02370 X 203 8 0.072230 0.281460 0.646310 1.40580 X 204 8 0.064366 0.269816 0.665819 0.97704 X 205 8 0.035993 0.191422 0.772585 2.09909 X 206 9 0.091638 0.215290 0.693073 1.03679 X 207 9 0.072171 0.236658 0.691171 2.36788 X 208 9 0.056187 0.195345 0.748468 3.45908 X 209 9 0.095888 0.229586 0.674526 3.63796 X 210 9 0.069809 0.219558 0.710633 2.56887 X 211 9 0.142920 0.223801 0.633279 2.00319 X 212 9 0.087323 0.196401 0.716276 2.40644 X 213 9 0.064517 0.218711 0.716772 2.58552 X 214 9 0.086882 0.194778 0.718341 8.94023 X 215 9 0.067463 0.219228 0.713309 3.75275 X 216 9 0.105610 0.230661 0.663730 0.34082 X 217 9 0.138992 0.283123 0.577885 1.62649 X 218 9 0.081364 0.186967 0.731670 2.31678 X 219 9 0.114535 0.221751 0.663714 1.77709 X 220 9 0.069940 0.280622 0.649438 1.38765 X 221 9 0.073137 0.143219 0.783643 3.46442 X 222 9 0.096326 0.243241 0.660434 1.74696 X 223 9 0.083284 0.202951 0.713765 1.28613 X 224 9 0.179133 0.299403 0.521465 1.15897 END_OF_FILE if test 16800 -ne `wc -c <'ch5eg1/electa.dat'`; then echo shar: \"'ch5eg1/electa.dat'\" unpacked with wrong size! fi # end of 'ch5eg1/electa.dat' fi if test -f 'realmat.cc' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'realmat.cc'\" else echo shar: Extracting \"'realmat.cc'\" \(13906 characters\) sed "s/^X//" >'realmat.cc' <<'END_OF_FILE' X/* ---------------------------------------------------------------------------- X Xrealmat: realmat.cc X Xrealmat is a C++ matrix class. The header file realmat.h describes its use. X XCopyright (C) 1990. X XA. Ronald Gallant XP.O. Box 5513 XRaleigh NC 27650-5513 XUSA X XPermission to use, copy, modify, and distribute this software and its Xdocumentation for any purpose and without fee is hereby granted, provided Xthat the above copyright notice appear in all copies and that both that Xcopyright notice and this permission notice appear in supporting Xdocumentation. X XThis software is provided "as is" without any expressed or implied warranty. X X---------------------------------------------------------------------------- */ X#include "realmat.h" X Xrealmat::realmat(INTEGER r, INTEGER c, REAL fill_value) X{ X if (r<=0) error("Error, realmat, realmat, Number of rows not positive"); X if (c<=0) error("Error, realmat, realmat, Number of columns not positive"); X rows=r; X cols=c; X len=r*c; X x = new REAL[len]; X if (x == 0) error("Error, realmat, realmat, Operator new failed"); X REAL* top = &(x[len]); X REAL* t = x; X while (t < top) *t++ = fill_value; X} X Xvoid realmat::resize(INTEGER r, INTEGER c, REAL fill_value) X{ X if (r<=0) error("Error, realmat, resize, Number of rows not positive"); X if (c<=0) error("Error, realmat, resize, Number of columns not positive"); X delete x; X rows=r; X cols=c; X len=r*c; X x = new REAL[len]; X if (x == 0) error("Error, realmat, resize, Operator new failed"); X REAL* top = &(x[len]); X REAL* t = x; X while (t < top) *t++ = fill_value; X} X Xrealmat::realmat(realmat& a) X{ X rows=a.rows; X cols=a.cols; X len=a.rows*a.cols; X if (len > 0) { X x = new REAL[len]; X if (x == 0) error("Error, realmat, realmat, Operator new failed"); X REAL* top = &(x[len]); X REAL* t = x; X REAL* u = a.x; X while (t < top) *t++ = *u++; X } X else { X rows = 0; cols = 0; len = 0; x = 0; X } X} X Xrealmat& realmat::operator=(realmat& a) X{ X rows=a.rows; X cols=a.cols; X len=a.len; X if (this != &a) { X delete x; X if (len > 0) { X x = new REAL[len]; X if (x == 0) error("Error, realmat, realmat, Operator new failed"); X REAL* top = &(x[len]); X REAL* t = x; X REAL* u = a.x; X while (t < top) *t++ = *u++; X } X else { X rows = 0; cols = 0; len = 0; x = 0; X } X } X return *this; X} X Xrealmat& realmat::operator+=(realmat& a) X{ X if ( cols != a.cols || rows != a.rows ) X error("Error, realmat, operator+=, Matrices not conformable."); X REAL* top = &(x[len]); X REAL* t = x; X REAL* u = a.x; X while (t < top) *t++ += *u++; X return *this; X} X Xrealmat& realmat::operator-=(realmat& a) X{ X if ( cols != a.cols || rows != a.rows ) X error("Error, realmat, operator-=, Matrices not conformable."); X REAL* top = &(x[len]); X REAL* t = x; X REAL* u = a.x; X while (t < top) *t++ -= *u++; X return *this; X} X Xrealmat& realmat::operator++() X{ X if (len == 0) X error("Error, realmat, operator++, Can't increment a null matrix"); X REAL* top = &(x[len]); X REAL* t = x; X while (t < top) ++(*t++); X return *this; X} X Xrealmat& realmat::operator--() X{ X if (len == 0) X error("Error, realmat, operator--, Can't decrement a null matrix"); X REAL* top = &(x[len]); X REAL* t = x; X while (t < top) --(*t++); X return *this; X} X Xint operator==(realmat& a, realmat& b) X{ X if ( a.cols != b.cols || a.rows != b.rows ) return 0; X if ( a.len == 0 ) return 1; X REAL* top = &(a.x[a.len]); X REAL* t = a.x; X REAL* u = b.x; X while (t < top) if(*t++ != *u++) return 0; X return 1; X} X Xint operator<(realmat& a, realmat& b) X{ X if ( a.cols != b.cols || a.rows != b.rows ) return 0; X if ( a.len == 0 ) return 0; X REAL* top = &(a.x[a.len]); X REAL* t = a.x; X REAL* u = b.x; X while (t < top) if(*t++ >= *u++) return 0; X return 1; X} X Xint operator<=(realmat& a, realmat& b) X{ X if ( a.cols != b.cols || a.rows != b.rows ) return 0; X if ( a.len == 0 ) return 1; X REAL* top = &(a.x[a.len]); X REAL* t = a.x; X REAL* u = b.x; X while (t < top) if(*t++ > *u++) return 0; X return 1; X} X Xint operator>(realmat& a, realmat& b) X{ X if ( a.cols != b.cols || a.rows != b.rows ) return 0; X if ( a.len == 0 ) return 0; X REAL* top = &(a.x[a.len]); X REAL* t = a.x; X REAL* u = b.x; X while (t < top) if(*t++ <= *u++) return 0; X return 1; X} X Xint operator>=(realmat& a, realmat& b) X{ X if ( a.cols != b.cols || a.rows != b.rows ) return 0; X if ( a.len == 0 ) return 1; X REAL* top = &(a.x[a.len]); X REAL* t = a.x; X REAL* u = b.x; X while (t < top) if(*t++ < *u++) return 0; X return 1; X} X XREAL& realmat::check1(INTEGER i) X{ X if ( i < 1 || len < i ) X error("Error, realmat, check1, Index out of range"); X return x[i-1]; X} X XREAL& realmat::check2(INTEGER i, INTEGER j) X{ X if ( i < 1 || rows < i || j < 1 || cols < j ) X error("Error, realmat, check2, Index out of range"); X return x[i + rows*j - rows - 1]; // return x[rows*(j-1)+i-1] X} X X#ifdef GNU_GPP_COMPILER Xrealmat T(realmat& a) return d; { //This prevents a constructor call on return. X#else Xrealmat T(realmat& a) { realmat d; X#endif X INTEGER newrows = a.cols; X INTEGER newcols = a.rows; X INTEGER newlen = newrows*newcols; X REAL* newx = new REAL[newlen]; X if (newx == 0) a.error("Error, realmat, T, Operator new failed"); X X /* X The intention: X for (INTEGER j = 0; j < a.cols; j++) X for (INTEGER i = 0; i < a.rows; i++) X newx[j + newrows*i] = a.x[i + a.rows*j]; X X The optimized version: X */ X REAL* t; X REAL* t_j = newx; X REAL* t_top = &(newx[newlen]); X REAL* u = a.x; X REAL* u_top = &(a.x[a.len]); X X while (u < u_top) { X t = t_j++; X while (t < t_top) { X *t = *u++; X t += newrows; X } X } X X d.resize(newrows,newcols,newx); X return d; X} X X#ifdef GNU_GPP_COMPILER Xrealmat operator+(realmat& a, realmat& b) return d; { X#else Xrealmat operator+(realmat& a, realmat& b) { realmat d; X#endif X if ( a.cols != b.cols || a.rows != b.rows ) X a.error("Error, realmat, operator+, Matrices not conformable."); X REAL* newx = new REAL[a.len]; X if (newx == 0) a.error("Error, realmat, operator+, Operator new failed"); X REAL* top = &(newx[a.len]); X REAL* t = newx; X REAL* u = a.x; X REAL* v = b.x; X while (t < top) *t++ = *u++ + *v++; X d.resize(a.rows,a.cols,newx); X return d; X} X X#ifdef GNU_GPP_COMPILER Xrealmat operator+(realmat& a) return d; { X#else Xrealmat operator+(realmat& a) { realmat d; X#endif X REAL* newx = new REAL[a.len]; X if (newx == 0) a.error("Error, realmat, operator+, Operator new failed"); X REAL* top = &(newx[a.len]); X REAL* t = newx; X REAL* u = a.x; X while (t < top) *t++ = *u++; X d.resize(a.rows,a.cols,newx); X return d; X} X X#ifdef GNU_GPP_COMPILER Xrealmat operator-(realmat& a, realmat& b) return d; { X#else Xrealmat operator-(realmat& a, realmat& b) { realmat d; X#endif X if ( a.cols != b.cols || a.rows != b.rows ) X a.error("Error, realmat, operator-, Matrices not conformable."); X REAL* newx = new REAL[a.len]; X if (newx == 0) a.error("Error, realmat, operator-, Operator new failed"); X REAL* top = &(newx[a.len]); X REAL* t = newx; X REAL* u = a.x; X REAL* v = b.x; X while (t < top) *t++ = *u++ - *v++; X d.resize(a.rows,a.cols,newx); X return d; X} X X#ifdef GNU_GPP_COMPILER Xrealmat operator-(realmat& a) return d; { X#else Xrealmat operator-(realmat& a) { realmat d; X#endif X REAL* newx = new REAL[a.len]; X if (newx == 0) a.error("Error, realmat, operator-, Operator new failed"); X REAL* top = &(newx[a.len]); X REAL* t = newx; X REAL* u = a.x; X while (t < top) *t++ = - *u++; X d.resize(a.rows,a.cols,newx); X return d; X} X X#ifdef GNU_GPP_COMPILER Xrealmat operator*(realmat& a, realmat& b) return d; { X#else Xrealmat operator*(realmat& a, realmat& b) { realmat d; X#endif X X INTEGER acols = a.cols; X INTEGER arows = a.rows; X INTEGER brows = b.rows; X INTEGER bcols = b.cols; X REAL* ax = a.x; X REAL* bx = b.x; X X if (acols != brows) X a.error("Error, realmat, operator*, Matrices not conformable."); X X INTEGER newlen = arows*bcols; X REAL* newx = new REAL[newlen]; X if (newx == 0) a.error("Error, realmat, operator*, Operator new failed"); X X if (newlen == (INTEGER) 1 ) { //Special case of an inner product. X REAL sum = REAL_ZERO; X REAL* stop = ax + acols; X while (ax < stop) sum += (*ax++) * (*bx++); X *newx = sum; X } X else { //General case. X REAL* top = &(newx[newlen]); X REAL* t = newx; X while (t < top) *t++ = REAL_ZERO; X X /* X The intention: X INTEGER i,j,k; X for (j = 0; j < b.cols; j++) X for (k = 0; k < a.cols; k++) X for (i = 0; i < a.rows; i++) X newx[i+a.rows*j] += a.x[i+a.rows*k] * b.x[k+b.rows*j]; X X The optimized version: X */ X REAL* aik; X REAL* bkj; X REAL* cij; X REAL* cij_sav; X REAL* ai_top; X REAL* bj_top = b.x; X REAL* cj_top = newx; X X while (cj_top < top) { X ai_top = ax; X bkj = bj_top; X bj_top += brows; X cij_sav = cj_top; X cj_top += arows; X X while (bkj < bj_top) { X aik = ai_top; X ai_top += arows; X cij = cij_sav; X X while (cij < cj_top) { X *cij++ += (*bkj)*(*aik++); X } X bkj++; X } X } X } X d.resize(arows,bcols,newx); X return d; X} X X#ifdef GNU_GPP_COMPILER Xrealmat operator*(REAL& a, realmat& b) return d; { X#else Xrealmat operator*(REAL& a, realmat& b) { realmat d; X#endif X REAL* newx = new REAL[b.len]; X if (newx == 0) b.error("Error, realmat, operator*, Operator new failed"); X REAL* top = &(newx[b.len]); X REAL* t = newx; X REAL* u = b.x; X while (t < top) *t++ = a * *u++; X d.resize(b.rows,b.cols,newx); X return d; X} X X#ifdef GNU_GPP_COMPILER Xrealmat operator*(INTEGER& a, realmat& b) return d; { X#else Xrealmat operator*(INTEGER& a, realmat& b) { realmat d; X#endif X REAL* newx = new REAL[b.len]; X if (newx == 0) b.error("Error, realmat, operator*, Operator new failed"); X REAL f = a; X REAL* top = &(newx[b.len]); X REAL* t = newx; X REAL* u = b.x; X while (t < top) *t++ = f * *u++; X d.resize(b.rows,b.cols,newx); X return d; X} X Xvoid default_realmat_error_handler(const char* msg) X{ X cerr << msg << "\n"; X exit(1); X} X XONE_ARG_ERROR_HANDLER_T realmat_error_handler = default_realmat_error_handler; X XONE_ARG_ERROR_HANDLER_T set_realmat_error_handler(ONE_ARG_ERROR_HANDLER_T f) X{ X ONE_ARG_ERROR_HANDLER_T old = realmat_error_handler; X realmat_error_handler = f; X return old; X} X Xvoid realmat::error(const char* msg) X{ X (*realmat_error_handler)(msg); X} X Xostream& operator<<(ostream& stream, realmat& a) X{ X char line[256], eol[2] = {'\n','\0'}, bl[3] = {' ','\n','\0'}; X char *next,*save,*fcode; X REAL f,af; X INTEGER i,j,linesize,maxcol,pad,start,stop,r,c; X X r=rows(a); X c=cols(a); X X X ((LINESIZE<72)||(LINESIZE>133)) ? (linesize=133) : (linesize=LINESIZE); X X if ( r <= 0 || c <= 0 ) { X pad = linesize/2 - 5; X memset(line,' ',pad); line[0]='\n'; line[pad]='\0'; X stream << line << "Null matrix\n"; X return stream; X } X X maxcol = (linesize - 8)/12; X if (c < maxcol) maxcol = c; X pad = (linesize - 8 - 12*maxcol)/2 + 1; X memset(line,' ',pad); X save = line + pad; X X X start=1; X do { X stop = start - 1 + maxcol; X if (stop > c) stop = c; X X stream << bl; X stream << bl; X X next = save; X memset(next,' ',6); X next = next + 6; X X for (j = start; j <= stop; j++) { X if (j < 1000) {sprintf(next," Col%3i",j); next = next+12;} X else if (j <32760) {sprintf(next," C%5i" ,j); next = next+12;} X else {sprintf(next," TooBig" ); next = next+12;} X } X X memcpy(next,eol,2); X X stream << line; X stream << bl; X X for (i = 1; i <= r; i++) { X next = save; X if (i < 1000) {sprintf(next,"Row%3i",i); next = next+6;} X else if (i <32760) {sprintf(next,"R%5i" ,i); next = next+6;} X else {sprintf(next,"TooBig" ); next = next+6;} X for (j = start; j <= stop; j++) { X fcode = "%12.3e"; X f = a.elem(i,j); X af = f; if (f < REAL_ZERO) af = -f ; X if (af < 1.E+8 ) fcode = "%12.0f"; X if (af < 1.E+5 ) fcode = "%12.1f"; X if (af < 1.E+4 ) fcode = "%12.2f"; X if (af < 1.E+3 ) fcode = "%12.3f"; X if (af < 1.E+2 ) fcode = "%12.4f"; X if (af < 1.E+1 ) fcode = "%12.5f"; X if (af < 1.E+0 ) fcode = "%12.6f"; X if (af < 1.E-1 ) fcode = "%12.7f"; X if (af < 1.E-2 ) fcode = "%12.8f"; X if (af < 1.E-4 ) fcode = "%12.3e"; X if (af < 1.E-30) fcode = "%12.1f"; X sprintf(next,fcode,f); next = next+12; X } X memcpy(next,eol,2); X stream << line; X } X start = stop + 1; X } while (stop < c); X X return stream; X} X X#ifdef GNU_GPP_COMPILER Xrealmat invpsd(realmat& a, REAL eps) return d; { X#else Xrealmat invpsd(realmat& a, REAL eps) { realmat d; X#endif X X if (a.rows <= 0 || a.cols <= 0) X a.error("Error, realmat, invpsd, Null matrix."); X X if (a.cols != a.rows) X a.error("Error, realmat, invpsd, Matrix not square."); X X /* X The call to sqrt in dcond will catch this: X for (INTEGER i=0; i'status.cc' <<'END_OF_FILE' X/* ---------------------------------------------------------------------------- X Xnlmdl: status.cc X Xnlmdl is a C++ implementation of the statistical methods in A. Ronald XGallant, "Nonlinear Statistical Models," New York: John Wiley and Sons, X1987, ISBN 0-471-80260-3, using a matrix class realmat that is distributed Xwith it. The header files nlmdl.h and realmat.h describe the use of the Xprogram and matrix class, respectively. X XCopyright (C) 1990. X XA. Ronald Gallant XP.O. Box 5513 XRaleigh NC 27650-5513 XUSA X XPermission to use, copy, modify, and distribute this software and its Xdocumentation for any purpose and without fee is hereby granted, provided Xthat the above copyright notice appear in all copies and that both that Xcopyright notice and this permission notice appear in supporting Xdocumentation. X XThis software is provided "as is" without any expressed or implied warranty. X X---------------------------------------------------------------------------- */ X/* status is a class used by nlmdl */ X X#include "status.h" X Xstatus::status() X{ X switches[0]='\0'; X method[0]='\0'; X n=0; X M=0; X K=0; X p=0; X itheta=0; X ivar=0; X vartype[0]='\0'; X MA=0; X weights[0]='\0'; X tol=0; X eps=0; X detail[0]='\0'; X rank=0; X df[0]='\0'; X obj=0; X starting=STARTING_FILENAME; X ending=ENDING_FILENAME; X} X Xstatus::~status() { } X X Xint status::from(char* filename) X{ X X#ifdef GNU_GPP_COMPILER X X#ifdef USE_ATT_STYLE_IO_WITH_GNU X filebuf ib; X if( ib.open(filename, input) == 0 ){ X cerr << "Error, status::from, Cannot open " << filename << "\n"; X exit(1); X } X istream ifrom(&ib); X#endif X X#ifdef USE_GNU_STYLE_IO_WITH_GNU X istream ifrom(filename,io_readonly,a_useonly); X if ( !ifrom ) { X cerr << "Error, status::from, Cannot open " << filename << "\n"; X exit(1); X } X#endif X X#endif X X#ifdef TURBO_CPP_COMPILER X ifstream ifrom(filename); X if ( !ifrom ) { X cerr << "Error, status::from, Cannot open " << filename << "\n"; X exit(1); X } X#endif X X char junk[MAX_STATUS_LINE]; X char newline; X INTEGER i,len; X double x = 0; X X if (!ifrom.get(switches,MAX_STATUS_LINE,'\n')) { X cerr << "Error, status::from, Bad data (switches) in " << filename << "\n"; X exit(1); X } X X if (!(ifrom>>method)) return 0; X if ( strcmp(method,"SUR") != 0 X && strcmp(method,"TSLS") != 0 X && strcmp(method,"GMM") != 0 ) { X cerr << "Error, status::from, Bad data (method) in " << filename << "\n"; X exit(1); X } X if (!ifrom.get(junk,MAX_STATUS_LINE,'\n')) return 0; X X if (!(ifrom>>n)) return 0; X if (!ifrom.get(junk,MAX_STATUS_LINE,'\n')) return 0; X X if (!(ifrom>>M)) return 0; X if (!ifrom.get(junk,MAX_STATUS_LINE,'\n')) return 0; X X if (!(ifrom>>K)) return 0; X if (!ifrom.get(junk,MAX_STATUS_LINE,'\n')) return 0; X X if (!(ifrom>>p)) return 0; X if (!ifrom.get(junk,MAX_STATUS_LINE,'\n')) return 0; X X if (!(ifrom>>itheta)) return 0; X if (!ifrom.get(junk,MAX_STATUS_LINE,'\n')) return 0; X X if (!(ifrom>>ivar)) return 0; X if (!ifrom.get(junk,MAX_STATUS_LINE,'\n')) return 0; X X if (!(ifrom>>vartype)) return 0; X if ( strcmp(vartype,"homoskedastic") != 0 X && strcmp(vartype,"heteroskedastic") != 0 ) { X cerr << "Error, status::from, Bad data (vartype) in " << filename << "\n"; X exit(1); X } X if (!ifrom.get(junk,MAX_STATUS_LINE,'\n')) return 0; X X if (!(ifrom>>MA)) return 0; X if ( MA<0 || MA>=n ) { X cerr << "Error, status::from, Bad data (MA) in " << filename << "\n"; X exit(1); X } X if (!ifrom.get(junk,MAX_STATUS_LINE,'\n')) return 0; X X if (!(ifrom>>weights)) return 0; X if ( strcmp(weights,"none") != 0 && strcmp(weights,"Parzen") != 0 ) { X cerr << "Error, status::from, Bad data (weights) in " << filename << "\n"; X exit(1); X } X if (!ifrom.get(junk,MAX_STATUS_LINE,'\n')) return 0; X X if (!(ifrom>>x)) X return 0; X else X tol = (REAL)x; X if (!ifrom.get(junk,MAX_STATUS_LINE,'\n')) return 0; X X if (!(ifrom>>x)) X return 0; X else X eps = (REAL)x; X if (!ifrom.get(junk,MAX_STATUS_LINE,'\n')) return 0; X X if (!(ifrom>>detail)) return 0; X if ( strcmp(detail,"none") != 0 X && strcmp(detail,"minimal") != 0 X && strcmp(detail,"full") != 0 ) { X cerr << "Error, status::from, Bad data (detail) in " << filename << "\n"; X exit(1); X } X if (!ifrom.get(junk,MAX_STATUS_LINE,'\n')) return 0; X X if (!ifrom.get(newline)) return 0; X X for (i=1; i<=5; i++) { X if (!ifrom.get(junk,MAX_STATUS_LINE,'\n')) return 0; X if (!ifrom.get(newline)) return 0; X } X X#ifdef GNU_GPP_COMPILER X if (!(ifrom>>x)) return 0; X if (!ifrom.get(junk,MAX_STATUS_LINE,'\n')) { X cerr << "Error, status::from, Bad data (theta) in " << filename << "\n"; X exit(1); X } X#endif X X#ifdef TURBO_CPP_COMPILER X double y; X if (!ifrom.get(junk,MAX_STATUS_LINE,'\n')) return 0; X if (sscanf(junk, "%le", &y) <= 0) return 0; X if (!ifrom.get(junk,MAX_STATUS_LINE,'\n')) { X cerr << "Error, status::from, Bad data (theta) in " << filename << "\n"; X exit(1); X } X x = y; X#endif X X theta.resize(p,(INTEGER)1,REAL_ZERO); X theta[1] = (REAL)x; X X if (p==1) return 0; X for (i=2; i<=p; i++) { X if(!(ifrom>>x)) { X cerr << "Error, status::from, Bad data (theta) in " << filename << "\n"; X exit(1); X } X else { X theta[i] = (REAL)x; X } X if (!ifrom.get(junk,MAX_STATUS_LINE,'\n')) { X cerr << "Error, status::from, Bad data (theta) in " << filename << "\n"; X exit(1); X } X } X X if (!ifrom.get(newline)) return 0; X X#ifdef GNU_GPP_COMPILER X if (!(ifrom>>x)) return 0; X if (!ifrom.get(junk,MAX_STATUS_LINE,'\n')) { X cerr << "Error, status::from, Bad data (var) in " << filename << "\n"; X exit(1); X } X#endif X X#ifdef TURBO_CPP_COMPILER X if (!ifrom.get(junk,MAX_STATUS_LINE,'\n')) return 0; X if (sscanf(junk, "%le", &y) <= 0) return 0; X if (!ifrom.get(junk,MAX_STATUS_LINE,'\n')) { X cerr << "Error, status::from, Bad data (var) in " << filename << "\n"; X exit(1); X } X x = y; X#endif X X if ( strcmp(method,"SUR") == 0 || strcmp(method,"TSLS") == 0) { X var.resize(M,M); X var[1] = (REAL)x; X if (M==1) return 0; X for (i=2; i<=M*M; i++) { X if(!(ifrom>>x)) { X cerr << "Error, status::from, Bad data (var) in " << filename << "\n"; X exit(1); X } X else { X var[i] = (REAL)x; X } X if (!ifrom.get(junk,MAX_STATUS_LINE,'\n')) { X cerr << "Error, status::from, Bad data (var) in " << filename << "\n"; X exit(1); X } X } X } X else if ( strcmp(method,"GMM") == 0 ) { X len = M*K; X var.resize(len,len); X var[1] = (REAL)x; X if (len==1) return 0; X for (i=2; i<=len*len; i++) { X if(!(ifrom>>x)) { X cerr << "Error, status::from, Bad data (var) in " << filename << "\n"; X exit(1); X } X else { X var[i] = (REAL)x; X } X if (!ifrom.get(junk,MAX_STATUS_LINE,'\n')) { X cerr << "Error, status::from, Bad data (var) in " << filename << "\n"; X exit(1); X } X } X } X else { X cerr << "Error, status::from, Bad data (var) in " << filename << "\n"; X exit(1); X } X X ifrom.close(); X return 0; X X} X Xint status::to(char* filename) X{ X X#ifdef GNU_GPP_COMPILER X X#ifdef USE_ATT_STYLE_IO_WITH_GNU X filebuf ob; X if( ob.open(filename, output) == 0 ) { X cerr << "Error, status::to, Cannot open " << filename << "\n"; X exit(1); X } X ostream oto(&ob); X#endif X X#ifdef USE_GNU_STYLE_IO_WITH_GNU X ostream oto(filename,io_writeonly,a_create); X if ( !oto ) { X cerr << "Error, status::to, Cannot open " << filename << "\n"; X exit(1); X } X#endif X X#endif X X#ifdef TURBO_CPP_COMPILER X ofstream oto(filename); X if ( !oto ) { X cerr << "Error, status::to, Cannot open " << filename << "\n"; X exit(1); X } X#endif X X INTEGER i,j,len; X char temp[MAX_STATUS_LINE]; X X if (strlen(switches) == 0) { X if (!(oto << "\tThis line is passed to class model as a string.\n")) { X cerr << "Error, status::to, Error writing to " << filename << "\n"; X exit(1); X } X } X else { X if (!(oto << switches << "\n")) { X cerr << "Error, status::to, Error writing to " << filename << "\n"; X exit(1); X } X } X X sprintf(temp,"%-19s", method); X if (!(oto<0 is unwise.\n")) { X cerr << "Error, status::to, Error writing to " << filename << "\n"; X exit(1); X } X X sprintf(temp,"%-19.6e",tol); X if (!(oto<