00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #if defined (__GNUG__) && defined (USE_PRAGMA_INTERFACE_IMPLEMENTATION)
00025 #pragma implementation
00026 #endif
00027
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031
00032 #include <cfloat>
00033 #include <vector>
00034
00035 #include "Array-util.h"
00036 #include "dNDArray.h"
00037 #include "mx-base.h"
00038 #include "f77-fcn.h"
00039 #include "lo-error.h"
00040 #include "lo-ieee.h"
00041 #include "lo-mappers.h"
00042
00043 #if defined (HAVE_FFTW3)
00044 #include "oct-fftw.h"
00045
00046 ComplexNDArray
00047 NDArray::fourier (int dim) const
00048 {
00049 dim_vector dv = dims ();
00050
00051 if (dim > dv.length () || dim < 0)
00052 return ComplexNDArray ();
00053
00054 int stride = 1;
00055 int n = dv(dim);
00056
00057 for (int i = 0; i < dim; i++)
00058 stride *= dv(i);
00059
00060 int howmany = numel () / dv (dim);
00061 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
00062 int nloop = (stride == 1 ? 1 : numel () / dv (dim) / stride);
00063 int dist = (stride == 1 ? n : 1);
00064
00065 const double *in (fortran_vec ());
00066 ComplexNDArray retval (dv);
00067 Complex *out (retval.fortran_vec ());
00068
00069
00070 for (int k = 0; k < nloop; k++)
00071 octave_fftw::fft (in + k * stride * n, out + k * stride * n,
00072 n, howmany, stride, dist);
00073
00074 return retval;
00075 }
00076
00077 ComplexNDArray
00078 NDArray::ifourier (int dim) const
00079 {
00080 dim_vector dv = dims ();
00081
00082 if (dim > dv.length () || dim < 0)
00083 return ComplexNDArray ();
00084
00085 int stride = 1;
00086 int n = dv(dim);
00087
00088 for (int i = 0; i < dim; i++)
00089 stride *= dv(i);
00090
00091 int howmany = numel () / dv (dim);
00092 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
00093 int nloop = (stride == 1 ? 1 : numel () / dv (dim) / stride);
00094 int dist = (stride == 1 ? n : 1);
00095
00096 ComplexNDArray retval (*this);
00097 Complex *out (retval.fortran_vec ());
00098
00099
00100 for (int k = 0; k < nloop; k++)
00101 octave_fftw::ifft (out + k * stride * n, out + k * stride * n,
00102 n, howmany, stride, dist);
00103
00104 return retval;
00105 }
00106
00107 ComplexNDArray
00108 NDArray::fourier2d (void) const
00109 {
00110 dim_vector dv = dims();
00111 if (dv.length () < 2)
00112 return ComplexNDArray ();
00113
00114 dim_vector dv2(dv(0), dv(1));
00115 const double *in = fortran_vec ();
00116 ComplexNDArray retval (dv);
00117 Complex *out = retval.fortran_vec ();
00118 int howmany = numel() / dv(0) / dv(1);
00119 int dist = dv(0) * dv(1);
00120
00121 for (int i=0; i < howmany; i++)
00122 octave_fftw::fftNd (in + i*dist, out + i*dist, 2, dv2);
00123
00124 return retval;
00125 }
00126
00127 ComplexNDArray
00128 NDArray::ifourier2d (void) const
00129 {
00130 dim_vector dv = dims();
00131 if (dv.length () < 2)
00132 return ComplexNDArray ();
00133
00134 dim_vector dv2(dv(0), dv(1));
00135 ComplexNDArray retval (*this);
00136 Complex *out = retval.fortran_vec ();
00137 int howmany = numel() / dv(0) / dv(1);
00138 int dist = dv(0) * dv(1);
00139
00140 for (int i=0; i < howmany; i++)
00141 octave_fftw::ifftNd (out + i*dist, out + i*dist, 2, dv2);
00142
00143 return retval;
00144 }
00145
00146 ComplexNDArray
00147 NDArray::fourierNd (void) const
00148 {
00149 dim_vector dv = dims ();
00150 int rank = dv.length ();
00151
00152 const double *in (fortran_vec ());
00153 ComplexNDArray retval (dv);
00154 Complex *out (retval.fortran_vec ());
00155
00156 octave_fftw::fftNd (in, out, rank, dv);
00157
00158 return retval;
00159 }
00160
00161 ComplexNDArray
00162 NDArray::ifourierNd (void) const
00163 {
00164 dim_vector dv = dims ();
00165 int rank = dv.length ();
00166
00167 ComplexNDArray tmp (*this);
00168 Complex *in (tmp.fortran_vec ());
00169 ComplexNDArray retval (dv);
00170 Complex *out (retval.fortran_vec ());
00171
00172 octave_fftw::ifftNd (in, out, rank, dv);
00173
00174 return retval;
00175 }
00176
00177 #else
00178
00179 extern "C"
00180 {
00181
00182
00183
00184
00185
00186 F77_RET_T
00187 F77_FUNC (cffti, CFFTI) (const int&, Complex*);
00188
00189 F77_RET_T
00190 F77_FUNC (cfftf, CFFTF) (const int&, Complex*, Complex*);
00191
00192 F77_RET_T
00193 F77_FUNC (cfftb, CFFTB) (const int&, Complex*, Complex*);
00194 }
00195
00196 ComplexNDArray
00197 NDArray::fourier (int dim) const
00198 {
00199 dim_vector dv = dims ();
00200
00201 if (dim > dv.length () || dim < 0)
00202 return ComplexNDArray ();
00203
00204 ComplexNDArray retval (dv);
00205 int npts = dv(dim);
00206 int nn = 4*npts+15;
00207 Array<Complex> wsave (nn);
00208 Complex *pwsave = wsave.fortran_vec ();
00209
00210 OCTAVE_LOCAL_BUFFER (Complex, tmp, npts);
00211
00212 int stride = 1;
00213
00214 for (int i = 0; i < dim; i++)
00215 stride *= dv(i);
00216
00217 int howmany = numel () / npts;
00218 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
00219 int nloop = (stride == 1 ? 1 : numel () / npts / stride);
00220 int dist = (stride == 1 ? npts : 1);
00221
00222 F77_FUNC (cffti, CFFTI) (npts, pwsave);
00223
00224 for (int k = 0; k < nloop; k++)
00225 {
00226 for (int j = 0; j < howmany; j++)
00227 {
00228 OCTAVE_QUIT;
00229
00230 for (int i = 0; i < npts; i++)
00231 tmp[i] = elem((i + k*npts)*stride + j*dist);
00232
00233 F77_FUNC (cfftf, CFFTF) (npts, tmp, pwsave);
00234
00235 for (int i = 0; i < npts; i++)
00236 retval ((i + k*npts)*stride + j*dist) = tmp[i];
00237 }
00238 }
00239
00240 return retval;
00241 }
00242
00243 ComplexNDArray
00244 NDArray::ifourier (int dim) const
00245 {
00246 dim_vector dv = dims ();
00247
00248 if (dim > dv.length () || dim < 0)
00249 return ComplexNDArray ();
00250
00251 ComplexNDArray retval (dv);
00252 int npts = dv(dim);
00253 int nn = 4*npts+15;
00254 Array<Complex> wsave (nn);
00255 Complex *pwsave = wsave.fortran_vec ();
00256
00257 OCTAVE_LOCAL_BUFFER (Complex, tmp, npts);
00258
00259 int stride = 1;
00260
00261 for (int i = 0; i < dim; i++)
00262 stride *= dv(i);
00263
00264 int howmany = numel () / npts;
00265 howmany = (stride == 1 ? howmany : (howmany > stride ? stride : howmany));
00266 int nloop = (stride == 1 ? 1 : numel () / npts / stride);
00267 int dist = (stride == 1 ? npts : 1);
00268
00269 F77_FUNC (cffti, CFFTI) (npts, pwsave);
00270
00271 for (int k = 0; k < nloop; k++)
00272 {
00273 for (int j = 0; j < howmany; j++)
00274 {
00275 OCTAVE_QUIT;
00276
00277 for (int i = 0; i < npts; i++)
00278 tmp[i] = elem((i + k*npts)*stride + j*dist);
00279
00280 F77_FUNC (cfftb, CFFTB) (npts, tmp, pwsave);
00281
00282 for (int i = 0; i < npts; i++)
00283 retval ((i + k*npts)*stride + j*dist) = tmp[i] /
00284 static_cast<double> (npts);
00285 }
00286 }
00287
00288 return retval;
00289 }
00290
00291 ComplexNDArray
00292 NDArray::fourier2d (void) const
00293 {
00294 dim_vector dv = dims();
00295 dim_vector dv2 (dv(0), dv(1));
00296 int rank = 2;
00297 ComplexNDArray retval (*this);
00298 int stride = 1;
00299
00300 for (int i = 0; i < rank; i++)
00301 {
00302 int npts = dv2(i);
00303 int nn = 4*npts+15;
00304 Array<Complex> wsave (nn);
00305 Complex *pwsave = wsave.fortran_vec ();
00306 Array<Complex> row (npts);
00307 Complex *prow = row.fortran_vec ();
00308
00309 int howmany = numel () / npts;
00310 howmany = (stride == 1 ? howmany :
00311 (howmany > stride ? stride : howmany));
00312 int nloop = (stride == 1 ? 1 : numel () / npts / stride);
00313 int dist = (stride == 1 ? npts : 1);
00314
00315 F77_FUNC (cffti, CFFTI) (npts, pwsave);
00316
00317 for (int k = 0; k < nloop; k++)
00318 {
00319 for (int j = 0; j < howmany; j++)
00320 {
00321 OCTAVE_QUIT;
00322
00323 for (int l = 0; l < npts; l++)
00324 prow[l] = retval ((l + k*npts)*stride + j*dist);
00325
00326 F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave);
00327
00328 for (int l = 0; l < npts; l++)
00329 retval ((l + k*npts)*stride + j*dist) = prow[l];
00330 }
00331 }
00332
00333 stride *= dv2(i);
00334 }
00335
00336 return retval;
00337 }
00338
00339 ComplexNDArray
00340 NDArray::ifourier2d (void) const
00341 {
00342 dim_vector dv = dims();
00343 dim_vector dv2 (dv(0), dv(1));
00344 int rank = 2;
00345 ComplexNDArray retval (*this);
00346 int stride = 1;
00347
00348 for (int i = 0; i < rank; i++)
00349 {
00350 int npts = dv2(i);
00351 int nn = 4*npts+15;
00352 Array<Complex> wsave (nn);
00353 Complex *pwsave = wsave.fortran_vec ();
00354 Array<Complex> row (npts);
00355 Complex *prow = row.fortran_vec ();
00356
00357 int howmany = numel () / npts;
00358 howmany = (stride == 1 ? howmany :
00359 (howmany > stride ? stride : howmany));
00360 int nloop = (stride == 1 ? 1 : numel () / npts / stride);
00361 int dist = (stride == 1 ? npts : 1);
00362
00363 F77_FUNC (cffti, CFFTI) (npts, pwsave);
00364
00365 for (int k = 0; k < nloop; k++)
00366 {
00367 for (int j = 0; j < howmany; j++)
00368 {
00369 OCTAVE_QUIT;
00370
00371 for (int l = 0; l < npts; l++)
00372 prow[l] = retval ((l + k*npts)*stride + j*dist);
00373
00374 F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave);
00375
00376 for (int l = 0; l < npts; l++)
00377 retval ((l + k*npts)*stride + j*dist) = prow[l] /
00378 static_cast<double> (npts);
00379 }
00380 }
00381
00382 stride *= dv2(i);
00383 }
00384
00385 return retval;
00386 }
00387
00388 ComplexNDArray
00389 NDArray::fourierNd (void) const
00390 {
00391 dim_vector dv = dims ();
00392 int rank = dv.length ();
00393 ComplexNDArray retval (*this);
00394 int stride = 1;
00395
00396 for (int i = 0; i < rank; i++)
00397 {
00398 int npts = dv(i);
00399 int nn = 4*npts+15;
00400 Array<Complex> wsave (nn);
00401 Complex *pwsave = wsave.fortran_vec ();
00402 Array<Complex> row (npts);
00403 Complex *prow = row.fortran_vec ();
00404
00405 int howmany = numel () / npts;
00406 howmany = (stride == 1 ? howmany :
00407 (howmany > stride ? stride : howmany));
00408 int nloop = (stride == 1 ? 1 : numel () / npts / stride);
00409 int dist = (stride == 1 ? npts : 1);
00410
00411 F77_FUNC (cffti, CFFTI) (npts, pwsave);
00412
00413 for (int k = 0; k < nloop; k++)
00414 {
00415 for (int j = 0; j < howmany; j++)
00416 {
00417 OCTAVE_QUIT;
00418
00419 for (int l = 0; l < npts; l++)
00420 prow[l] = retval ((l + k*npts)*stride + j*dist);
00421
00422 F77_FUNC (cfftf, CFFTF) (npts, prow, pwsave);
00423
00424 for (int l = 0; l < npts; l++)
00425 retval ((l + k*npts)*stride + j*dist) = prow[l];
00426 }
00427 }
00428
00429 stride *= dv(i);
00430 }
00431
00432 return retval;
00433 }
00434
00435 ComplexNDArray
00436 NDArray::ifourierNd (void) const
00437 {
00438 dim_vector dv = dims ();
00439 int rank = dv.length ();
00440 ComplexNDArray retval (*this);
00441 int stride = 1;
00442
00443 for (int i = 0; i < rank; i++)
00444 {
00445 int npts = dv(i);
00446 int nn = 4*npts+15;
00447 Array<Complex> wsave (nn);
00448 Complex *pwsave = wsave.fortran_vec ();
00449 Array<Complex> row (npts);
00450 Complex *prow = row.fortran_vec ();
00451
00452 int howmany = numel () / npts;
00453 howmany = (stride == 1 ? howmany :
00454 (howmany > stride ? stride : howmany));
00455 int nloop = (stride == 1 ? 1 : numel () / npts / stride);
00456 int dist = (stride == 1 ? npts : 1);
00457
00458 F77_FUNC (cffti, CFFTI) (npts, pwsave);
00459
00460 for (int k = 0; k < nloop; k++)
00461 {
00462 for (int j = 0; j < howmany; j++)
00463 {
00464 OCTAVE_QUIT;
00465
00466 for (int l = 0; l < npts; l++)
00467 prow[l] = retval ((l + k*npts)*stride + j*dist);
00468
00469 F77_FUNC (cfftb, CFFTB) (npts, prow, pwsave);
00470
00471 for (int l = 0; l < npts; l++)
00472 retval ((l + k*npts)*stride + j*dist) = prow[l] /
00473 static_cast<double> (npts);
00474 }
00475 }
00476
00477 stride *= dv(i);
00478 }
00479
00480 return retval;
00481 }
00482
00483 #endif
00484
00485
00486
00487 boolNDArray
00488 NDArray::operator ! (void) const
00489 {
00490 boolNDArray b (dims ());
00491
00492 for (int i = 0; i < length (); i++)
00493 b.elem (i) = ! elem (i);
00494
00495 return b;
00496 }
00497
00498 bool
00499 NDArray::any_element_is_negative (bool neg_zero) const
00500 {
00501 int nel = nelem ();
00502
00503 if (neg_zero)
00504 {
00505 for (int i = 0; i < nel; i++)
00506 if (lo_ieee_signbit (elem (i)))
00507 return true;
00508 }
00509 else
00510 {
00511 for (int i = 0; i < nel; i++)
00512 if (elem (i) < 0)
00513 return true;
00514 }
00515
00516 return false;
00517 }
00518
00519
00520 bool
00521 NDArray::any_element_is_inf_or_nan (void) const
00522 {
00523 int nel = nelem ();
00524
00525 for (int i = 0; i < nel; i++)
00526 {
00527 double val = elem (i);
00528 if (xisinf (val) || xisnan (val))
00529 return true;
00530 }
00531
00532 return false;
00533 }
00534
00535 bool
00536 NDArray::all_elements_are_int_or_inf_or_nan (void) const
00537 {
00538 int nel = nelem ();
00539
00540 for (int i = 0; i < nel; i++)
00541 {
00542 double val = elem (i);
00543 if (xisnan (val) || D_NINT (val) == val)
00544 continue;
00545 else
00546 return false;
00547 }
00548
00549 return true;
00550 }
00551
00552
00553
00554
00555 bool
00556 NDArray::all_integers (double& max_val, double& min_val) const
00557 {
00558 int nel = nelem ();
00559
00560 if (nel > 0)
00561 {
00562 max_val = elem (0);
00563 min_val = elem (0);
00564 }
00565 else
00566 return false;
00567
00568 for (int i = 0; i < nel; i++)
00569 {
00570 double val = elem (i);
00571
00572 if (val > max_val)
00573 max_val = val;
00574
00575 if (val < min_val)
00576 min_val = val;
00577
00578 if (D_NINT (val) != val)
00579 return false;
00580 }
00581
00582 return true;
00583 }
00584
00585 bool
00586 NDArray::too_large_for_float (void) const
00587 {
00588 int nel = nelem ();
00589
00590 for (int i = 0; i < nel; i++)
00591 {
00592 double val = elem (i);
00593
00594 if (val > FLT_MAX || val < FLT_MIN)
00595 return true;
00596 }
00597
00598 return false;
00599 }
00600
00601
00602
00603 boolNDArray
00604 NDArray::all (int dim) const
00605 {
00606 MX_ND_ANY_ALL_REDUCTION (MX_ND_ALL_EVAL (MX_ND_ALL_EXPR), true);
00607 }
00608
00609 boolNDArray
00610 NDArray::any (int dim) const
00611 {
00612 MX_ND_ANY_ALL_REDUCTION (MX_ND_ANY_EVAL (MX_ND_ANY_EXPR), false);
00613 }
00614
00615 NDArray
00616 NDArray::cumprod (int dim) const
00617 {
00618 MX_ND_CUMULATIVE_OP (NDArray, double, 1, *);
00619 }
00620
00621 NDArray
00622 NDArray::cumsum (int dim) const
00623 {
00624 MX_ND_CUMULATIVE_OP (NDArray, double, 0, +);
00625 }
00626
00627 NDArray
00628 NDArray::prod (int dim) const
00629 {
00630 MX_ND_REAL_OP_REDUCTION (*= elem (iter_idx), 1);
00631 }
00632
00633 NDArray
00634 NDArray::sumsq (int dim) const
00635 {
00636 MX_ND_REAL_OP_REDUCTION (+= std::pow (elem (iter_idx), 2), 0);
00637 }
00638
00639 NDArray
00640 NDArray::sum (int dim) const
00641 {
00642 MX_ND_REAL_OP_REDUCTION (+= elem (iter_idx), 0);
00643 }
00644
00645 NDArray
00646 NDArray::max (int dim) const
00647 {
00648 ArrayN<int> dummy_idx;
00649 return max (dummy_idx, dim);
00650 }
00651
00652 NDArray
00653 NDArray::max (ArrayN<int>& idx_arg, int dim) const
00654 {
00655 dim_vector dv = dims ();
00656 dim_vector dr = dims ();
00657
00658 if (dv.numel () == 0 || dim > dv.length () || dim < 0)
00659 return NDArray ();
00660
00661 dr(dim) = 1;
00662
00663 NDArray result (dr);
00664 idx_arg.resize (dr);
00665
00666 int x_stride = 1;
00667 int x_len = dv(dim);
00668 for (int i = 0; i < dim; i++)
00669 x_stride *= dv(i);
00670
00671 for (int i = 0; i < dr.numel (); i++)
00672 {
00673 int x_offset;
00674 if (x_stride == 1)
00675 x_offset = i * x_len;
00676 else
00677 {
00678 int x_offset2 = 0;
00679 x_offset = i;
00680 while (x_offset >= x_stride)
00681 {
00682 x_offset -= x_stride;
00683 x_offset2++;
00684 }
00685 x_offset += x_offset2 * x_stride * x_len;
00686 }
00687
00688 int idx_j;
00689
00690 double tmp_max = octave_NaN;
00691
00692 for (idx_j = 0; idx_j < x_len; idx_j++)
00693 {
00694 tmp_max = elem (idx_j * x_stride + x_offset);
00695
00696 if (! octave_is_NaN_or_NA (tmp_max))
00697 break;
00698 }
00699
00700 for (int j = idx_j+1; j < x_len; j++)
00701 {
00702 double tmp = elem (j * x_stride + x_offset);
00703
00704 if (octave_is_NaN_or_NA (tmp))
00705 continue;
00706 else if (tmp > tmp_max)
00707 {
00708 idx_j = j;
00709 tmp_max = tmp;
00710 }
00711 }
00712
00713 result.elem (i) = tmp_max;
00714 idx_arg.elem (i) = octave_is_NaN_or_NA (tmp_max) ? 0 : idx_j;
00715 }
00716
00717 return result;
00718 }
00719
00720 NDArray
00721 NDArray::min (int dim) const
00722 {
00723 ArrayN<int> dummy_idx;
00724 return min (dummy_idx, dim);
00725 }
00726
00727 NDArray
00728 NDArray::min (ArrayN<int>& idx_arg, int dim) const
00729 {
00730 dim_vector dv = dims ();
00731 dim_vector dr = dims ();
00732
00733 if (dv.numel () == 0 || dim > dv.length () || dim < 0)
00734 return NDArray ();
00735
00736 dr(dim) = 1;
00737
00738 NDArray result (dr);
00739 idx_arg.resize (dr);
00740
00741 int x_stride = 1;
00742 int x_len = dv(dim);
00743 for (int i = 0; i < dim; i++)
00744 x_stride *= dv(i);
00745
00746 for (int i = 0; i < dr.numel (); i++)
00747 {
00748 int x_offset;
00749 if (x_stride == 1)
00750 x_offset = i * x_len;
00751 else
00752 {
00753 int x_offset2 = 0;
00754 x_offset = i;
00755 while (x_offset >= x_stride)
00756 {
00757 x_offset -= x_stride;
00758 x_offset2++;
00759 }
00760 x_offset += x_offset2 * x_stride * x_len;
00761 }
00762
00763 int idx_j;
00764
00765 double tmp_min = octave_NaN;
00766
00767 for (idx_j = 0; idx_j < x_len; idx_j++)
00768 {
00769 tmp_min = elem (idx_j * x_stride + x_offset);
00770
00771 if (! octave_is_NaN_or_NA (tmp_min))
00772 break;
00773 }
00774
00775 for (int j = idx_j+1; j < x_len; j++)
00776 {
00777 double tmp = elem (j * x_stride + x_offset);
00778
00779 if (octave_is_NaN_or_NA (tmp))
00780 continue;
00781 else if (tmp < tmp_min)
00782 {
00783 idx_j = j;
00784 tmp_min = tmp;
00785 }
00786 }
00787
00788 result.elem (i) = tmp_min;
00789 idx_arg.elem (i) = octave_is_NaN_or_NA (tmp_min) ? 0 : idx_j;
00790 }
00791
00792 return result;
00793 }
00794
00795 NDArray
00796 NDArray::concat (const NDArray& rb, const Array<int>& ra_idx)
00797 {
00798 if (rb.numel () > 0)
00799 insert (rb, ra_idx);
00800 return *this;
00801 }
00802
00803 ComplexNDArray
00804 NDArray::concat (const ComplexNDArray& rb, const Array<int>& ra_idx)
00805 {
00806 ComplexNDArray retval (*this);
00807 if (rb.numel () > 0)
00808 retval.insert (rb, ra_idx);
00809 return retval;
00810 }
00811
00812 charNDArray
00813 NDArray::concat (const charNDArray& rb, const Array<int>& ra_idx)
00814 {
00815 charNDArray retval (dims ());
00816 int nel = numel ();
00817
00818 for (int i = 0; i < nel; i++)
00819 {
00820 double d = elem (i);
00821
00822 if (xisnan (d))
00823 {
00824 (*current_liboctave_error_handler)
00825 ("invalid conversion from NaN to character");
00826 return retval;
00827 }
00828 else
00829 {
00830 int ival = NINT (d);
00831
00832 if (ival < 0 || ival > UCHAR_MAX)
00833
00834
00835 ival = 0;
00836
00837 retval.elem (i) = static_cast<char>(ival);
00838 }
00839 }
00840
00841 if (rb.numel () == 0)
00842 return retval;
00843
00844 retval.insert (rb, ra_idx);
00845 return retval;
00846 }
00847
00848 NDArray
00849 real (const ComplexNDArray& a)
00850 {
00851 int a_len = a.length ();
00852 NDArray retval;
00853 if (a_len > 0)
00854 retval = NDArray (mx_inline_real_dup (a.data (), a_len), a.dims ());
00855 return retval;
00856 }
00857
00858 NDArray
00859 imag (const ComplexNDArray& a)
00860 {
00861 int a_len = a.length ();
00862 NDArray retval;
00863 if (a_len > 0)
00864 retval = NDArray (mx_inline_imag_dup (a.data (), a_len), a.dims ());
00865 return retval;
00866 }
00867
00868 NDArray&
00869 NDArray::insert (const NDArray& a, int r, int c)
00870 {
00871 Array<double>::insert (a, r, c);
00872 return *this;
00873 }
00874
00875 NDArray&
00876 NDArray::insert (const NDArray& a, const Array<int>& ra_idx)
00877 {
00878 Array<double>::insert (a, ra_idx);
00879 return *this;
00880 }
00881
00882 NDArray
00883 NDArray::abs (void) const
00884 {
00885 NDArray retval (dims ());
00886
00887 int nel = nelem ();
00888
00889 for (int i = 0; i < nel; i++)
00890 retval(i) = fabs (elem (i));
00891
00892 return retval;
00893 }
00894
00895 Matrix
00896 NDArray::matrix_value (void) const
00897 {
00898 Matrix retval;
00899
00900 int nd = ndims ();
00901
00902 switch (nd)
00903 {
00904 case 1:
00905 retval = Matrix (Array2<double> (*this, dimensions(0), 1));
00906 break;
00907
00908 case 2:
00909 retval = Matrix (Array2<double> (*this, dimensions(0), dimensions(1)));
00910 break;
00911
00912 default:
00913 (*current_liboctave_error_handler)
00914 ("invalid conversion of NDArray to Matrix");
00915 break;
00916 }
00917
00918 return retval;
00919 }
00920
00921 void
00922 NDArray::increment_index (Array<int>& ra_idx,
00923 const dim_vector& dimensions,
00924 int start_dimension)
00925 {
00926 ::increment_index (ra_idx, dimensions, start_dimension);
00927 }
00928
00929 int
00930 NDArray::compute_index (Array<int>& ra_idx,
00931 const dim_vector& dimensions)
00932 {
00933 return ::compute_index (ra_idx, dimensions);
00934 }
00935
00936
00937 std::ostream&
00938 operator << (std::ostream& os, const NDArray& a)
00939 {
00940 int nel = a.nelem ();
00941
00942 for (int i = 0; i < nel; i++)
00943 {
00944 os << " ";
00945 octave_write_double (os, a.elem (i));
00946 os << "\n";
00947 }
00948 return os;
00949 }
00950
00951 std::istream&
00952 operator >> (std::istream& is, NDArray& a)
00953 {
00954 int nel = a.nelem ();
00955
00956 if (nel < 1 )
00957 is.clear (std::ios::badbit);
00958 else
00959 {
00960 double tmp;
00961 for (int i = 0; i < nel; i++)
00962 {
00963 tmp = octave_read_double (is);
00964 if (is)
00965 a.elem (i) = tmp;
00966 else
00967 goto done;
00968 }
00969 }
00970
00971 done:
00972
00973 return is;
00974 }
00975
00976
00977
00978
00979 #define EMPTY_RETURN_CHECK(T) \
00980 if (nel == 0) \
00981 return T (dv);
00982
00983 NDArray
00984 min (double d, const NDArray& m)
00985 {
00986 dim_vector dv = m.dims ();
00987 int nel = dv.numel ();
00988
00989 EMPTY_RETURN_CHECK (NDArray);
00990
00991 NDArray result (dv);
00992
00993 for (int i = 0; i < nel; i++)
00994 {
00995 OCTAVE_QUIT;
00996 result (i) = xmin (d, m (i));
00997 }
00998
00999 return result;
01000 }
01001
01002 NDArray
01003 min (const NDArray& m, double d)
01004 {
01005 dim_vector dv = m.dims ();
01006 int nel = dv.numel ();
01007
01008 EMPTY_RETURN_CHECK (NDArray);
01009
01010 NDArray result (dv);
01011
01012 for (int i = 0; i < nel; i++)
01013 {
01014 OCTAVE_QUIT;
01015 result (i) = xmin (d, m (i));
01016 }
01017
01018 return result;
01019 }
01020
01021 NDArray
01022 min (const NDArray& a, const NDArray& b)
01023 {
01024 dim_vector dv = a.dims ();
01025 int nel = dv.numel ();
01026
01027 if (dv != b.dims ())
01028 {
01029 (*current_liboctave_error_handler)
01030 ("two-arg min expecting args of same size");
01031 return NDArray ();
01032 }
01033
01034 EMPTY_RETURN_CHECK (NDArray);
01035
01036 NDArray result (dv);
01037
01038 for (int i = 0; i < nel; i++)
01039 {
01040 OCTAVE_QUIT;
01041 result (i) = xmin (a (i), b (i));
01042 }
01043
01044 return result;
01045 }
01046
01047 NDArray
01048 max (double d, const NDArray& m)
01049 {
01050 dim_vector dv = m.dims ();
01051 int nel = dv.numel ();
01052
01053 EMPTY_RETURN_CHECK (NDArray);
01054
01055 NDArray result (dv);
01056
01057 for (int i = 0; i < nel; i++)
01058 {
01059 OCTAVE_QUIT;
01060 result (i) = xmax (d, m (i));
01061 }
01062
01063 return result;
01064 }
01065
01066 NDArray
01067 max (const NDArray& m, double d)
01068 {
01069 dim_vector dv = m.dims ();
01070 int nel = dv.numel ();
01071
01072 EMPTY_RETURN_CHECK (NDArray);
01073
01074 NDArray result (dv);
01075
01076 for (int i = 0; i < nel; i++)
01077 {
01078 OCTAVE_QUIT;
01079 result (i) = xmax (d, m (i));
01080 }
01081
01082 return result;
01083 }
01084
01085 NDArray
01086 max (const NDArray& a, const NDArray& b)
01087 {
01088 dim_vector dv = a.dims ();
01089 int nel = dv.numel ();
01090
01091 if (dv != b.dims ())
01092 {
01093 (*current_liboctave_error_handler)
01094 ("two-arg max expecting args of same size");
01095 return NDArray ();
01096 }
01097
01098 EMPTY_RETURN_CHECK (NDArray);
01099
01100 NDArray result (dv);
01101
01102 for (int i = 0; i < nel; i++)
01103 {
01104 OCTAVE_QUIT;
01105 result (i) = xmax (a (i), b (i));
01106 }
01107
01108 return result;
01109 }
01110
01111 NDS_CMP_OPS(NDArray, , double, )
01112 NDS_BOOL_OPS(NDArray, double, 0.0)
01113
01114 SND_CMP_OPS(double, , NDArray, )
01115 SND_BOOL_OPS(double, NDArray, 0.0)
01116
01117 NDND_CMP_OPS(NDArray, , NDArray, )
01118 NDND_BOOL_OPS(NDArray, NDArray, 0.0)
01119
01120
01121
01122
01123
01124