"Fossies" - the Fresh Open Source Software Archive

Member "PDL-2.025/Basic/Core/Core.pm" (19 Nov 2020, 108598 Bytes) of package /linux/misc/PDL-2.025.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Perl source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "Core.pm" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 2.024_vs_2.025.

    1 package PDL::Core;
    2 
    3 # Core routines for PDL module
    4 
    5 use strict;
    6 use warnings;
    7 use PDL::Exporter;
    8 require PDL; # for $VERSION
    9 use DynaLoader;
   10 our @ISA    = qw( PDL::Exporter DynaLoader );
   11 our $VERSION = '2.025';
   12 bootstrap PDL::Core $VERSION;
   13 use PDL::Types ':All';
   14 use Config;
   15 
   16 our @EXPORT = qw( piddle pdl null barf ); # Only stuff always exported!
   17 my @convertfuncs = map PDL::Types::typefld($_,'convertfunc'), PDL::Types::typesrtkeys();
   18 my @exports_internal = qw(howbig threadids topdl);
   19 my @exports_normal   = (@EXPORT,
   20   @convertfuncs,
   21   qw(nelem dims shape null
   22       convert inplace zeroes zeros ones list listindices unpdl
   23       set at flows thread_define over reshape dog cat barf type diagonal
   24       dummy mslice approx flat sclr squeeze
   25       get_autopthread_targ set_autopthread_targ get_autopthread_actual
   26       get_autopthread_size set_autopthread_size) );
   27 our @EXPORT_OK = (@exports_internal, @exports_normal);
   28 our %EXPORT_TAGS = (
   29    Func     => [@exports_normal],
   30    Internal => [@exports_internal] );
   31 
   32 our ($level, @dims, $sep, $sep2, $match);
   33 
   34 # Important variables (place in PDL namespace)
   35 # (twice to eat "used only once" warning)
   36 
   37 $PDL::debug      =       # Debugging info
   38 $PDL::debug      = 0;
   39 $PDL::verbose      =         # Functions provide chatty information
   40 $PDL::verbose      = 0;
   41 $PDL::use_commas   = 0;        # Whether to insert commas when printing arrays
   42 $PDL::floatformat  = "%7g";    # Default print format for long numbers
   43 $PDL::doubleformat = "%10.8g";
   44 $PDL::indxformat   = "%12d";   # Default print format for PDL_Indx values
   45 $PDL::undefval     = 0;        # Value to use instead of undef when creating PDLs
   46 $PDL::toolongtoprint = 10000;  # maximum pdl size to stringify for printing
   47 
   48 ################ Exportable functions of the Core ######################
   49 
   50 # log10() is now defined in ops.pd
   51 
   52 *howbig       = \&PDL::howbig;    *unpdl    = \&PDL::unpdl;
   53 *nelem        = \&PDL::nelem;     *inplace  = \&PDL::inplace;
   54 *dims         = \&PDL::dims;      *list     = \&PDL::list;
   55 *threadids    = \&PDL::threadids; *listindices  = \&PDL::listindices;
   56 *null         = \&PDL::null;      *set      = \&PDL::set;
   57 *at     = \&PDL::at;      *flows    = \&PDL::flows;
   58 *sclr           = \&PDL::sclr;    *shape        = \&PDL::shape;
   59 
   60 for (map {
   61   [ PDL::Types::typefld($_,'convertfunc'), PDL::Types::typefld($_,'numval') ]
   62 } PDL::Types::typesrtkeys()) {
   63   my ($conv, $val) = @$_;
   64   no strict 'refs';
   65   *$conv = *{"PDL::$conv"} = sub {
   66     return bless [$val], "PDL::Type" unless @_;
   67     alltopdl('PDL', (scalar(@_)>1 ? [@_] : shift), PDL::Type->new($val));
   68   };
   69 }
   70 
   71 BEGIN {
   72     *thread_define = \&PDL::thread_define;
   73     *convert      = \&PDL::convert;   *over      = \&PDL::over;
   74     *dog          = \&PDL::dog;       *cat           = \&PDL::cat;
   75     *type         = \&PDL::type;      *approx        = \&PDL::approx;
   76     *diagonal     = \&PDL::diagonal;
   77     *dummy        = \&PDL::dummy;
   78     *mslice       = \&PDL::mslice;
   79     *isempty      = \&PDL::isempty;
   80     *string       = \&PDL::string;
   81 }
   82 
   83 =head1 NAME
   84 
   85 PDL::Core - fundamental PDL functionality and vectorization/threading
   86 
   87 =head1 DESCRIPTION
   88 
   89 Methods and functions for type conversions, PDL creation,
   90 type conversion, threading etc.
   91 
   92 =head1 SYNOPSIS
   93 
   94  use PDL::Core;             # Normal routines
   95  use PDL::Core ':Internal'; # Hairy routines
   96 
   97 =head1 VECTORIZATION/THREADING: METHOD AND NOMENCLATURE
   98 
   99 PDL provides vectorized operations via a built-in engine.
  100 Vectorization is called "threading" for historical reasons.
  101 The threading engine implements simple rules for each operation.
  102 
  103 Each PDL object has a "shape" that is a generalized N-dimensional
  104 rectangle defined by a "dim list" of sizes in an arbitrary
  105 set of dimensions.  A PDL with shape 2x3 has 6 elements and is
  106 said to be two-dimensional, or may be referred to as a 2x3-PDL.
  107 The dimensions are indexed numerically starting at 0, so a
  108 2x3-PDL has a dimension 0 (or "dim 0") with size 2 and a 1 dimension
  109 (or "dim 1") with size 3.
  110 
  111 PDL generalizes *all* mathematical operations with the notion of
  112 "active dims": each operator has zero or more active dims that are
  113 used in carrying out the operation.  Simple scalar operations like
  114 scalar multiplication ('*') have 0 active dims.  More complicated
  115 operators can have more active dims.  For example, matrix
  116 multiplication ('x') has 2 active dims.  Additional dims are
  117 automatically vectorized across -- e.g. multiplying a 2x5-PDL with a
  118 2x5-PDL requires 10 simple multiplication operations, and yields a
  119 2x5-PDL result.
  120 
  121 =head2 Threading rules
  122 
  123 In any PDL expression, the active dims appropriate for each operator
  124 are used starting at the 0 dim and working forward through the dim
  125 list of each object.  All additional dims after the active dims are
  126 "thread dims".  The thread dims do not have to agree exactly: they are
  127 coerced to agree according to simple rules:
  128 
  129 =over 3
  130 
  131 =item * Null PDLs match any dim list (see below).
  132 
  133 =item * Dims with sizes other than 1 must all agree in size.
  134 
  135 =item * Dims of size 1 are expanded as necessary.
  136 
  137 =item * Missing dims are expanded appropriately.
  138 
  139 =back
  140 
  141 The "size 1" rule implements "generalized scalar" operation, by
  142 analogy to scalar multiplication.  The "missing dims" rule
  143 acknowledges the ambiguity between a missing dim and a dim of size 1.
  144 
  145 =head2 Null PDLs
  146 
  147 PDLs on the left-hand side of assignment can have the special value
  148 "Null".  A null PDL has no dim list and no set size; its shape is
  149 determined by the computed shape of the expression being assigned to
  150 it.   Null PDLs contain no values and can only be assigned to.  When
  151 assigned to (e.g. via the C<.=> operator), they cease to be null PDLs.
  152 
  153 To create a null PDL, use C<PDL-E<gt>null()>.
  154 
  155 =head2 Empty PDLs
  156 
  157 PDLs can represent the empty set using "structured Empty" variables.
  158 An empty PDL is not a null PDL.
  159 
  160 Any dim of a PDL can be set explicitly to size 0.  If so, the PDL
  161 contains zero values (because the total number of values is the
  162 product of all the sizes in the PDL's shape or dimlist).
  163 
  164 Scalar PDLs are zero-dimensional and have no entries in the dim list,
  165 so they cannot be empty.  1-D and higher PDLs can be empty.  Empty
  166 PDLs are useful for set operations, and are most commonly encountered
  167 in the output from selection operators such as L<which|PDL::Primitive>
  168 and L<whichND|PDL::Primitive>.  Not all empty PDLs have the same
  169 threading properties -- e.g. a 2x0-PDL represents a collection of
  170 2-vectors that happens to contain no elements, while a simple 0-PDL
  171 represents a collection of scalar values (that also happens to contain
  172 no elements).
  173 
  174 Note that 0 dims are not adjustable via the threading rules -- a dim
  175 with size 0 can only match a corresponding dim of size 0 or 1.
  176 
  177 =head2 Thread rules and assignments
  178 
  179 Versions of PDL through 2.4.10 have some irregularity with threading and
  180 assignments.  Currently the threading engine performs a full expansion of
  181 both sides of the computed assignment operator C<.=> (which assigns values
  182 to a pre-existing PDL).  This leads to counter-intuitive behavior in
  183 some cases:
  184 
  185 =over 3
  186 
  187 =item * Generalized scalars and computed assignment
  188 
  189 If the PDL on the left-hand side of C<.=> has a dim of size 1, it can be
  190 treated as a generalized scalar, as in:
  191 
  192     $x = sequence(2,3);
  193     $y = zeroes(1,3);
  194     $y .= $x;
  195 
  196 In this case, C<$y> is automatically treated as a 2x3-PDL during the
  197 threading operation, but half of the values from C<$x> silently disappear.
  198 The output is, as Kernighan and Ritchie would say, "undefined".
  199 
  200 Further, if the value on the right of C<.=> is empty, then C<.=> becomes
  201 a silent no-op:
  202 
  203     $x = zeroes(0);
  204     $y = zeroes(1);
  205     $y .= $x+1;
  206     print $y;
  207 
  208 will print C<[0]>.  In this case, "$x+1" is empty, and "$y" is a generalized
  209 scalar that is adjusted to be empty, so the assignment is carried out for
  210 zero elements (a no-op).
  211 
  212 Both of these behaviors are considered harmful and should not be relied upon:
  213 they may be patched away in a future version of PDL.
  214 
  215 =item * Empty PDLs and generalized scalars
  216 
  217 Generalized scalars (PDLs with a dim of size 1) can match any size in the
  218 corresponding dim, including 0.  Thus,
  219 
  220     $x = ones(2,0);
  221     $y = sequence(2,1);
  222     $c = $x * $y;
  223     print $c;
  224 
  225 prints C<Empty[2,0]>.
  226 
  227 This behavior is counterintuitive but desirable, and will be preserved
  228 in future versions of PDL.
  229 
  230 =back
  231 
  232 =head1 VARIABLES
  233 
  234 These are important variables of B<global> scope and are placed
  235 in the PDL namespace.
  236 
  237 =head3 C<$PDL::debug>
  238 
  239 =over 4
  240 
  241 When true, PDL debugging information is printed.
  242 
  243 =back
  244 
  245 =head3 C<$PDL::verbose>
  246 
  247 =over 4
  248 
  249 When true, PDL functions provide chatty information.
  250 
  251 =back
  252 
  253 =head3 C<$PDL::use_commas>
  254 
  255 =over 4
  256 
  257 Whether to insert commas when printing pdls
  258 
  259 =back
  260 
  261 =head3 C<$PDL::floatformat>, C<$PDL::doubleformat>, C<$PDL::indxformat>
  262 
  263 =over 4
  264 
  265 The default print format for floats, doubles, and indx values,
  266 respectively.  The default default values are:
  267 
  268   $PDL::floatformat  = "%7g";
  269   $PDL::doubleformat = "%10.8g";
  270   $PDL::indxformat   = "%12d";
  271 
  272 =back
  273 
  274 =head3 C<$PDL::undefval>
  275 
  276 =over 4
  277 
  278 The value to use instead of C<undef> when creating pdls.
  279 
  280 =back
  281 
  282 =head3 C<$PDL::toolongtoprint>
  283 
  284 =over 4
  285 
  286 The maximal size pdls to print (defaults to 10000 elements)
  287 
  288 =back
  289 
  290 =head1 FUNCTIONS
  291 
  292 
  293 =head2 barf
  294 
  295 =for ref
  296 
  297 Standard error reporting routine for PDL.
  298 
  299 C<barf()> is the routine PDL modules should call to report errors. This
  300 is because C<barf()> will report the error as coming from the correct
  301 line in the module user's script rather than in the PDL module.
  302 
  303 For now, barf just calls Carp::confess()
  304 
  305 Remember C<barf()> is your friend. *Use* it!
  306 
  307 =for example
  308 
  309 At the perl level:
  310 
  311  barf("User has too low an IQ!");
  312 
  313 In C or XS code:
  314 
  315  barf("You have made %d errors", count);
  316 
  317 Note: this is one of the few functions ALWAYS exported
  318 by PDL::Core
  319 
  320 =cut
  321 
  322 use Carp;
  323 sub barf { goto &Carp::confess }
  324 sub cluck { goto &Carp::cluck }
  325 *PDL::barf  = \&barf;
  326 *PDL::cluck = \&cluck;
  327 
  328 ########## Set Auto-PThread Based On Environment Vars ############
  329 PDL::set_autopthread_targ( $ENV{PDL_AUTOPTHREAD_TARG} ) if( defined ( $ENV{PDL_AUTOPTHREAD_TARG} ) );
  330 PDL::set_autopthread_size( $ENV{PDL_AUTOPTHREAD_SIZE} ) if( defined ( $ENV{PDL_AUTOPTHREAD_SIZE} ) );
  331 ##################################################################
  332 
  333 =head2 pdl
  334 
  335 =for ref
  336 
  337 PDL constructor - creates new piddle from perl scalars/arrays, piddles, and strings
  338 
  339 =for usage
  340 
  341  $double_pdl = pdl(SCALAR|ARRAY REFERENCE|ARRAY|STRING);  # default type
  342  $type_pdl   = pdl(PDL::Type,SCALAR|ARRAY REFERENCE|ARRAY|STRING);
  343 
  344 =for example
  345 
  346  $x = pdl [1..10];                    # 1D array
  347  $x = pdl ([1..10]);                  # 1D array
  348  $x = pdl (1,2,3,4);                  # Ditto
  349  $y = pdl [[1,2,3],[4,5,6]];          # 2D 3x2 array
  350  $y = pdl "[[1,2,3],[4,5,6]]";        # Ditto (slower)
  351  $y = pdl "[1 2 3; 4 5 6]";           # Ditto
  352  $y = pdl q[1 2 3; 4 5 6];            # Ditto, using the q quote operator
  353  $y = pdl "1 2 3; 4 5 6";             # Ditto, less obvious, but still works
  354  $y = pdl 42                          # 0-dimensional scalar
  355  $c = pdl $x;                         # Make a new copy
  356 
  357  $u = pdl ushort(), 42                # 0-dimensional ushort scalar
  358  $y = pdl(byte(),[[1,2,3],[4,5,6]]);  # 2D byte piddle
  359 
  360  $n = pdl indx(), [1..5];             # 1D array of indx values
  361  $n = pdl indx, [1..5];               # ... can leave off parens
  362  $n = indx( [1..5] );                 # ... still the same!
  363 
  364  $x = pdl([1,2,3],[4,5,6]);           # 2D
  365  $x = pdl([1,2,3],[4,5,6]);           # 2D
  366 
  367 Note the last two are equivalent - a list is automatically
  368 converted to a list reference for syntactic convenience. i.e. you
  369 can omit the outer C<[]>
  370 
  371 You can mix and match arrays, array refs, and PDLs in your argument
  372 list, and C<pdl> will sort them out.  You get back a PDL whose last
  373 (slowest running) dim runs across the top level of the list you hand
  374 in, and whose first (fastest running) dim runs across the deepest
  375 level that you supply.
  376 
  377 At the moment, you cannot mix and match those arguments with string
  378 arguments, though we can't imagine a situation in which you would
  379 really want to do that.
  380 
  381 The string version of pdl also allows you to use the strings C<bad>, C<inf>,
  382 and C<nan>, and it will insert the values that you mean (and set the bad flag
  383 if you use C<bad>). You can mix and match case, though you shouldn't. Here are
  384 some examples:
  385 
  386  $bad = pdl q[1 2 3 bad 5 6];  # Set fourth element to the bad value
  387  $bad = pdl q[1 2 3 BAD 5 6];  # ditto
  388  $bad = pdl q[1 2 inf bad 5];  # now third element is IEEE infinite value
  389  $bad = pdl q[nan 2 inf -inf]; # first value is IEEE nan value
  390 
  391 The default constructor uses IEEE double-precision floating point numbers. You
  392 can use other types, but you will get a warning if you try to use C<nan> with
  393 integer types (it will be replaced with the C<bad> value) and you will get a
  394 fatal error if you try to use C<inf>.
  395 
  396 Throwing a PDL into the mix has the same effect as throwing in a list ref:
  397 
  398   pdl(pdl(1,2),[3,4])
  399 
  400 is the same as
  401 
  402   pdl([1,2],[3,4]).
  403 
  404 All of the dimensions in the list are "padded-out" with undefval to
  405 meet the widest dim in the list, so (e.g.)
  406 
  407   $x = pdl([[1,2,3],[2]])
  408 
  409 gives you the same answer as
  410 
  411   $x = pdl([[1,2,3],[2,undef,undef]]);
  412 
  413 If your PDL module has bad values compiled into it (see L<PDL::Bad>), 
  414 you can pass BAD values into the constructor within pre-existing PDLs.
  415 The BAD values are automatically kept BAD and propagated correctly.
  416 
  417 C<pdl()> is a functional synonym for the 'new' constructor,
  418 e.g.:
  419 
  420  $x = new PDL [1..10];
  421 
  422 In order to control how undefs are handled in converting from perl lists to
  423 PDLs, one can set the variable C<$PDL::undefval>.
  424 For example:
  425 
  426  $foo = [[1,2,undef],[undef,3,4]];
  427  $PDL::undefval = -999;
  428  $f = pdl $foo;
  429  print $f
  430  [
  431   [   1    2 -999]
  432   [-999    3    4]
  433  ]
  434 
  435 C<$PDL::undefval> defaults to zero.
  436 
  437 As a final note, if you include an Empty PDL in the list of objects to
  438 construct into a PDL, it is kept as a placeholder pane -- so if you feed
  439 in (say) 7 objects, you get a size of 7 in the 0th dim of the output PDL.
  440 The placeholder panes are completely padded out.  But if you feed in only
  441 a single Empty PDL, you get back the Empty PDL (no padding).
  442 
  443 =cut
  444 
  445 sub pdl {PDL->pdl(@_)}
  446 
  447 sub piddle {PDL->pdl(@_)}
  448 
  449 =head2 null
  450 
  451 =for ref
  452 
  453 Returns a 'null' piddle.
  454 
  455 =for usage
  456 
  457  $x = null;
  458 
  459 C<null()> has a special meaning to L<PDL::PP|PDL::PP>. It is used to
  460 flag a special kind of empty piddle, which can grow to
  461 appropriate dimensions to store a result (as opposed to
  462 storing a result in an existing piddle).
  463 
  464 =for example
  465 
  466  pdl> sumover sequence(10,10), $ans=null;p $ans
  467  [45 145 245 345 445 545 645 745 845 945]
  468 
  469 =cut
  470 
  471 sub PDL::null{
  472     my $class = scalar(@_) ? shift : undef; # if this sub called with no
  473                         #  class ( i.e. like 'null()', instead
  474                         #  of '$obj->null' or 'CLASS->null', setup
  475 
  476     if( defined($class) ){
  477         $class = ref($class) || $class;  # get the class name
  478     }
  479     else{
  480         $class = 'PDL';  # set class to the current package name if null called
  481                     # with no arguments
  482     }
  483 
  484     return $class->initialize();
  485 }
  486 
  487 =head2 nullcreate
  488 
  489 =for ref
  490 
  491 Returns a 'null' piddle.
  492 
  493 =for usage
  494 
  495  $x = PDL->nullcreate($arg)
  496 
  497 This is an routine used by many of the threading primitives
  498 (i.e. L<sumover|PDL::Ufunc/sumover>,
  499 L<minimum|PDL::Ufunc/minimum>, etc.) to generate a null piddle for the
  500 function's output that will behave properly for derived (or
  501 subclassed) PDL objects.
  502 
  503 For the above usage:
  504 If C<$arg> is a PDL, or a derived PDL, then C<$arg-E<gt>null> is returned.
  505 If C<$arg> is a scalar (i.e. a zero-dimensional PDL) then C<PDL-E<gt>null>
  506 is returned.
  507 
  508 =for example
  509 
  510  PDL::Derived->nullcreate(10)
  511    returns PDL::Derived->null.
  512  PDL->nullcreate($pdlderived)
  513    returns $pdlderived->null.
  514 
  515 =cut
  516 
  517 sub PDL::nullcreate{
  518     my ($type,$arg) = @_;
  519         return ref($arg) ? $arg->null : $type->null ;
  520 }
  521 
  522 =head2 nelem
  523 
  524 =for ref
  525 
  526 Return the number of elements in a piddle
  527 
  528 =for usage
  529 
  530  $n = nelem($piddle); $n = $piddle->nelem;
  531 
  532 =for example
  533 
  534  $mean = sum($data)/nelem($data);
  535 
  536 =head2 dims
  537 
  538 =for ref
  539 
  540 Return piddle dimensions as a perl list
  541 
  542 =for usage
  543 
  544  @dims = $piddle->dims;  @dims = dims($piddle);
  545 
  546 =for example
  547 
  548  pdl> p @tmp = dims zeroes 10,3,22
  549  10 3 22
  550 
  551 See also L<shape|shape> which returns a piddle instead.
  552 
  553 =head2 shape
  554 
  555 =for ref
  556 
  557 Return piddle dimensions as a piddle
  558 
  559 =for usage
  560 
  561  $shape = $piddle->shape;  $shape = shape($piddle);
  562 
  563 =for example
  564 
  565  pdl> p $shape = shape zeroes 10,3,22
  566  [10 3 22]
  567 
  568 See also L<dims|dims> which returns a perl list.
  569 
  570 =head2 ndims
  571 
  572 =for ref
  573 
  574 Returns the number of dimensions in a piddle. Alias
  575 for L<getndims|PDL::Core/getndims>.
  576 
  577 =head2 getndims
  578 
  579 =for ref
  580 
  581 Returns the number of dimensions in a piddle
  582 
  583 =for usage
  584 
  585  $ndims = $piddle->getndims;
  586 
  587 =for example
  588 
  589  pdl> p zeroes(10,3,22)->getndims
  590  3
  591 
  592 =head2 dim
  593 
  594 =for ref
  595 
  596 Returns the size of the given dimension of a piddle. Alias
  597 for L<getdim|PDL::Core/getdim>.
  598 
  599 =head2 getdim
  600 
  601 =for ref
  602 
  603 Returns the size of the given dimension.
  604 
  605 =for usage
  606 
  607  $dim0 = $piddle->getdim(0);
  608 
  609 =for example
  610 
  611  pdl> p zeroes(10,3,22)->getdim(1)
  612  3
  613 
  614 Negative indices count from the end of the dims array.
  615 Indices beyond the end will return a size of 1. This
  616 reflects the idea that any pdl is equivalent to an
  617 infinitely dimensional array in which only a finite number of
  618 dimensions have a size different from one. For example, in that sense a
  619 3D piddle of shape [3,5,2] is equivalent to a [3,5,2,1,1,1,1,1,....]
  620 piddle. Accordingly,
  621 
  622   print $x->getdim(10000);
  623 
  624 will print 1 for most practically encountered piddles.
  625 
  626 =head2 topdl
  627 
  628 =for ref
  629 
  630 alternate piddle constructor - ensures arg is a piddle
  631 
  632 =for usage
  633 
  634  $x = topdl(SCALAR|ARRAY REFERENCE|ARRAY);
  635 
  636 The difference between L<pdl()|/pdl> and C<topdl()> is that the
  637 latter will just 'fall through' if the argument is
  638 already a piddle. It will return a reference and I<NOT>
  639 a new copy.
  640 
  641 This is particularly useful if you are writing a function
  642 which is doing some fiddling with internals and assumes
  643 a piddle argument (e.g. for method calls). Using C<topdl()>
  644 will ensure nothing breaks if passed with '2'.
  645 
  646 Note that C<topdl()> is not exported by default (see example
  647 below for usage).
  648 
  649 =for example
  650 
  651  use PDL::Core ':Internal'; # use the internal routines of
  652                             # the Core module
  653 
  654  $x = topdl 43;             # $x is piddle with value '43'
  655  $y = topdl $piddle;        # fall through
  656  $x = topdl (1,2,3,4);      # Convert 1D array
  657 
  658 =head2 get_datatype
  659 
  660 =for ref
  661 
  662 Internal: Return the numeric value identifying the piddle datatype
  663 
  664 =for usage
  665 
  666  $x = $piddle->get_datatype;
  667 
  668 Mainly used for internal routines.
  669 
  670 NOTE: get_datatype returns 'just a number' not any special
  671 type object, unlike L<type|/type>.
  672 
  673 =head2 howbig
  674 
  675 =for ref
  676 
  677 Returns the sizeof a piddle datatype in bytes.
  678 
  679 Note that C<howbig()> is not exported by default (see example
  680 below for usage).
  681 
  682 =for usage
  683 
  684  use PDL::Core ':Internal'; # use the internal routines of
  685                             # the Core module
  686 
  687  $size = howbig($piddle->get_datatype);
  688 
  689 Mainly used for internal routines.
  690 
  691 NOTE: NOT a method! This is because get_datatype returns
  692 'just a number' not any special object.
  693 
  694 =for example
  695 
  696  pdl> p howbig(ushort([1..10])->get_datatype)
  697  2
  698 
  699 
  700 =head2 get_dataref
  701 
  702 =for ref
  703 
  704 Return the internal data for a piddle, as a perl SCALAR ref.
  705 
  706 Most piddles hold their internal data in a packed perl string, to take
  707 advantage of perl's memory management.  This gives you direct access
  708 to the string, which is handy when you need to manipulate the binary
  709 data directly (e.g. for file I/O).  If you modify the string, you'll
  710 need to call L<upd_data|upd_data> afterward, to make sure that the
  711 piddle points to the new location of the underlying perl variable.
  712 
  713 Calling C<get_dataref> automatically physicalizes your piddle (see
  714 L<make_physical|/PDL::make_physical>).  You definitely
  715 don't want to do anything to the SV to truncate or deallocate the
  716 string, unless you correspondingly call L<reshape|/reshape> to make the
  717 PDL match its new data dimension.
  718 
  719 You definitely don't want to use get_dataref unless you know what you
  720 are doing (or are trying to find out): you can end up scrozzling
  721 memory if you shrink or eliminate the string representation of the
  722 variable.  Here be dragons.
  723 
  724 =head2 upd_data
  725 
  726 =for ref
  727 
  728 Update the data pointer in a piddle to match its perl SV.
  729 
  730 This is useful if you've been monkeying with the packed string
  731 representation of the PDL, which you probably shouldn't be doing
  732 anyway.  (see L<get_dataref|get_dataref>.)
  733 
  734 =cut
  735 
  736 sub topdl {PDL->topdl(@_)}
  737 
  738 ####################### Overloaded operators #######################
  739 
  740 # This is to used warn if an operand is non-numeric or non-PDL.
  741 sub warn_non_numeric_op_wrapper {
  742     my ($cb, $op_name) = @_;
  743     return sub {
  744         my ($op1, $op2) = @_;
  745         unless( Scalar::Util::looks_like_number($op2)
  746             || ( Scalar::Util::blessed($op2) && $op2->isa('PDL') )
  747             ) {
  748             warn "'$op2' is not numeric nor a PDL in operator $op_name";
  749         };
  750         $cb->(@_);
  751     }
  752 }
  753 
  754 { package PDL;
  755   # use UNIVERSAL 'isa'; # need that later in info function
  756   use Carp;
  757 
  758   use overload (
  759         "+"     => \&PDL::plus,     # in1, in2
  760         "*"     => \&PDL::mult, # in1, in2
  761         "-"     => \&PDL::minus,    # in1, in2, swap if true
  762         "/"     => \&PDL::divide,   # in1, in2, swap if true
  763 
  764         "+="    => sub { PDL::plus     ($_[0], $_[1], $_[0], 0); $_[0]; }, # in1, in2, out, swap if true
  765         "*="    => sub { PDL::mult ($_[0], $_[1], $_[0], 0); $_[0]; }, # in1, in2, out, swap if true
  766         "-="    => sub { PDL::minus    ($_[0], $_[1], $_[0], 0); $_[0]; }, # in1, in2, out, swap if true
  767         "/="    => sub { PDL::divide   ($_[0], $_[1], $_[0], 0); $_[0]; }, # in1, in2, out, swap if true
  768 
  769         ">"     => \&PDL::gt,       # in1, in2, swap if true
  770         "<"     => \&PDL::lt,       # in1, in2, swap if true
  771         "<="    => \&PDL::le,       # in1, in2, swap if true
  772         ">="    => \&PDL::ge,       # in1, in2, swap if true
  773         "=="    => \&PDL::eq,       # in1, in2
  774         "eq"    => PDL::Core::warn_non_numeric_op_wrapper(\&PDL::eq, 'eq'),
  775                                     # in1, in2
  776         "!="    => \&PDL::ne,       # in1, in2
  777 
  778         "<<"    => \&PDL::shiftleft,  # in1, in2, swap if true
  779         ">>"    => \&PDL::shiftright, # in1, in2, swap if true
  780         "|"     => \&PDL::or2,        # in1, in2
  781         "&"     => \&PDL::and2,       # in1, in2
  782         "^"     => \&PDL::xor,        # in1, in2
  783 
  784         "<<="   => sub { PDL::shiftleft ($_[0], $_[1], $_[0], 0); $_[0]; }, # in1, in2, out, swap if true
  785         ">>="   => sub { PDL::shiftright($_[0], $_[1], $_[0], 0); $_[0]; }, # in1, in2, out, swap if true
  786         "|="    => sub { PDL::or2      ($_[0], $_[1], $_[0], 0); $_[0]; }, # in1, in2, out, swap if true
  787         "&="    => sub { PDL::and2     ($_[0], $_[1], $_[0], 0); $_[0]; }, # in1, in2, out, swap if true
  788         "^="    => sub { PDL::xor       ($_[0], $_[1], $_[0], 0); $_[0]; }, # in1, in2, out, swap if true
  789             "**="   => sub { PDL::power     ($_[0], $_[1], $_[0], 0); $_[0]; }, # in1, in2, out, swap if true
  790             "%="    => sub { PDL::modulo    ($_[0], $_[1], $_[0], 0); $_[0]; }, # in1, in2, out, swap if true
  791 
  792         "sqrt"  => sub { PDL::sqrt ($_[0]); },
  793         "abs"   => sub { PDL::abs  ($_[0]); },
  794         "sin"   => sub { PDL::sin  ($_[0]); },
  795         "cos"   => sub { PDL::cos  ($_[0]); },
  796 
  797         "!"     => sub { PDL::not  ($_[0]); },
  798         "~"     => sub { PDL::bitnot ($_[0]); },
  799 
  800         "log"   => sub { PDL::log   ($_[0]); },
  801         "exp"   => sub { PDL::exp   ($_[0]); },
  802 
  803             "**"    => \&PDL::power,          # in1, in2, swap if true
  804 
  805             "atan2" => \&PDL::atan2,          # in1, in2, swap if true
  806             "%"     => \&PDL::modulo,         # in1, in2, swap if true
  807 
  808             "<=>"   => \&PDL::spaceship,      # in1, in2, swap if true
  809 
  810         "="     =>  sub {$_[0]},          # Don't deep copy, just copy reference
  811 
  812         ".="    => sub {
  813                         my @args = reverse &PDL::Core::rswap;
  814                         PDL::Ops::assgn(@args);
  815                         return $args[1];
  816                     },
  817 
  818         'x'     =>  sub{my $foo = $_[0]->null();
  819                   PDL::Primitive::matmult(@_[0,1],$foo); $foo;},
  820 
  821         'bool'  => sub { return 0 if $_[0]->isnull;
  822                  croak("multielement piddle in conditional expression (see PDL::FAQ questions 6-10 and 6-11)")
  823                      unless $_[0]->nelem == 1;
  824                  $_[0]->clump(-1)->at(0); },
  825         "\"\""  =>  \&PDL::Core::string   );
  826 }
  827 
  828 sub rswap { if($_[2]) { return @_[1,0]; } else { return @_[0,1]; } }
  829 
  830 ##################### Data type/conversion stuff ########################
  831 
  832 
  833 # XXX Optimize!
  834 
  835 sub PDL::dims {  # Return dimensions as @list
  836    my $pdl = PDL->topdl (shift);
  837    my @dims = ();
  838    for(0..$pdl->getndims()-1) {push @dims,($pdl->getdim($_))}
  839    return @dims;
  840 }
  841 
  842 sub PDL::shape {  # Return dimensions as a pdl
  843    my $pdl = PDL->topdl (shift);
  844    my @dims = ();
  845    for(0..$pdl->getndims()-1) {push @dims,($pdl->getdim($_))}
  846    return indx(\@dims);
  847 }
  848 
  849 sub PDL::howbig {
  850     my $t = shift;
  851     if("PDL::Type" eq ref $t) {$t = $t->[0]}
  852     PDL::howbig_c($t);
  853 }
  854 
  855 =head2 threadids
  856 
  857 =for ref
  858 
  859 Returns the piddle thread IDs as a perl list
  860 
  861 Note that C<threadids()> is not exported by default (see example
  862 below for usage).
  863 
  864 =for usage
  865 
  866  use PDL::Core ':Internal'; # use the internal routines of
  867                             # the Core module
  868 
  869  @ids = threadids $piddle;
  870 
  871 =cut
  872 
  873 sub PDL::threadids {  # Return dimensions as @list
  874    my $pdl = PDL->topdl (shift);
  875    my @dims = ();
  876    for(0..$pdl->getnthreadids()) {push @dims,($pdl->getthreadid($_))}
  877    return @dims;
  878 }
  879 
  880 ################# Creation/copying functions #######################
  881 
  882 
  883 sub PDL::pdl { my $x = shift; return $x->new(@_) }
  884 
  885 =head2 doflow
  886 
  887 =for ref
  888 
  889 Turn on/off dataflow
  890 
  891 =for usage
  892 
  893  $x->doflow;  doflow($x);
  894 
  895 =cut
  896 
  897 sub PDL::doflow {
  898     my $this = shift;
  899     $this->set_dataflow_f(1);
  900     $this->set_dataflow_b(1);
  901 }
  902 
  903 =head2 flows
  904 
  905 =for ref
  906 
  907 Whether or not a piddle is indulging in dataflow
  908 
  909 =for usage
  910 
  911  something if $x->flows; $hmm = flows($x);
  912 
  913 =cut
  914 
  915 sub PDL::flows {
  916     my $this = shift;
  917          return ($this->fflows || $this->bflows);
  918 }
  919 
  920 =head2 new
  921 
  922 =for ref
  923 
  924 new piddle constructor method
  925 
  926 =for usage
  927 
  928  $x = PDL->new(SCALAR|ARRAY|ARRAY REF|STRING);
  929 
  930 =for example
  931 
  932  $x = PDL->new(42);             # new from a Perl scalar
  933  $x = new PDL 42;               # ditto
  934  $y = PDL->new(@list_of_vals);  # new from Perl list
  935  $y = new PDL @list_of_vals;    # ditto
  936  $z = PDL->new(\@list_of_vals); # new from Perl list reference
  937  $w = PDL->new("[1 2 3]");      # new from Perl string, using
  938                                 # Matlab constructor syntax
  939 
  940 Constructs piddle from perl numbers and lists
  941 and strings with Matlab/Octave style constructor
  942 syntax.
  943 
  944 The string input is fairly versatile though not
  945 performance optimized. The goal is to make it
  946 easy to copy and paste code from PDL output and
  947 to offer a familiar Matlab syntax for piddle
  948 construction. As of May, 2010, it is a new
  949 feature, so feel free to report bugs or suggest
  950 new features.  See documentation for L<pdl> for
  951 more examples of usage.
  952 
  953 
  954 =cut
  955 
  956 use Scalar::Util;       # for looks_like_number test
  957 use Carp 'carp';        # for carping (warnings in caller's context)
  958 
  959 # This is the code that handles string arguments. It has now gotten quite large,
  960 # so here's the basic explanation. I want to allow expressions like 2, 1e3, +4,
  961 # bad, nan, inf, and more. Checking this can get tricky. This croaks when it
  962 # finds:
  963 # 1) strings of e or E that are longer than 1 character long (like eeee)
  964 # 2) non-supported characters or strings
  965 # 3) expressions that are syntactically erroneous, like '1 2 3 ]', which has an
  966 #    extra bracket
  967 # 4) use of inf when the data type does not support inf (i.e. the integers)
  968 
  969 sub PDL::Core::new_pdl_from_string {
  970    my ($new, $original_value, $this, $type) = @_;
  971    my $value = $original_value;
  972 
  973    # Check for input that would generate empty piddles as output:
  974    my @types = PDL::Types::types;
  975    return zeroes($types[$type], 1)->where(zeroes(1) < 0)
  976       if ($value eq '' or $value eq '[]');
  977 
  978    # I check for invalid characters later, but arbitrary strings of e will
  979    # pass that check, so I'll check for that here, first.
  980 #   croak("PDL::Core::new_pdl_from_string: I found consecutive copies of e but\n"
  981 #      . "  I'm not sure what you mean. You gave me $original_value")
  982 #      if ($value =~ /ee/i);
  983    croak("PDL::Core::new_pdl_from_string: found 'e' as part of a larger word in $original_value")
  984       if $value =~ /e\p{IsAlpha}/ or $value =~ /\p{IsAlpha}e/;
  985 
  986    # Only a few characters are allowed in the expression, but we want to allow
  987    # expressions like 'inf' and 'bad'. As such, convert those values to internal
  988    # representations that will pass the invalid-character check. We'll replace
  989    # them with Perl-evalute-able strings in a little bit. Here, I represent
  990    #  bad => EE
  991    #  nan => ee
  992    #  inf => Ee
  993    #  pi  => eE
  994    # --( Bad )--
  995    croak("PDL::Core::new_pdl_from_string: found 'bad' as part of a larger word in $original_value")
  996       if $value =~ /bad\B/ or $value =~ /\Bbad/;
  997    my ($has_bad) = ($value =~ s/\bbad\b/EE/gi);
  998    # --( nan )--
  999    my ($has_nan) = 0;
 1000    croak("PDL::Core::new_pdl_from_string: found 'nan' as part of a larger word in $original_value")
 1001       if $value =~ /\Bnan/ or $value =~ /nan\B/;
 1002    $has_nan++ if ($value =~ s/\bnan\b/ee/gi);
 1003    # Strawberry Perl compatibility:
 1004    croak("PDL::Core::new_pdl_from_string: found '1.#IND' as part of a larger word in $original_value")
 1005       if $value =~ /IND\B/i;
 1006    $has_nan++ if ($value =~ s/1\.\#IND/ee/gi);
 1007    # --( inf )--
 1008    my ($has_inf) = 0;
 1009    # Strawberry Perl compatibility:
 1010    croak("PDL::Core::new_pdl_from_string: found '1.#INF' as part of a larger word in $original_value")
 1011       if $value =~ /INF\B/i;
 1012    $has_inf++ if ($value =~ s/1\.\#INF/Ee/gi);
 1013    # Other platforms:
 1014    croak("PDL::Core::new_pdl_from_string: found 'inf' as part of a larger word in $original_value")
 1015       if $value =~ /inf\B/ or $value =~ /\Binf/;
 1016    $has_inf++ if ($value =~ s/\binf\b/Ee/gi);
 1017    # --( pi )--
 1018    croak("PDL::Core::new_pdl_from_string: found 'pi' as part of a larger word in $original_value")
 1019       if $value =~ /pi\B/ or $value =~ /\Bpi/;
 1020    $value =~ s/\bpi\b/eE/gi;
 1021 
 1022    # Some data types do not support nan and inf, so check for and warn or croak,
 1023    # as appropriate:
 1024    if ($has_nan and not $types[$type]->usenan) {
 1025       carp("PDL::Core::new_pdl_from_string: no nan for type $types[$type]; converting to bad value");
 1026       $value =~ s/ee/EE/g;
 1027       $has_bad += $has_nan;
 1028       $has_nan = 0;
 1029    }
 1030    croak("PDL::Core::new_pdl_from_string: type $types[$type] does not support inf")
 1031       if ($has_inf and not $types[$type]->usenan);
 1032 
 1033    # Make the white-space uniform and see if any not-allowed characters are
 1034    # present:
 1035    $value =~ s/\s+/ /g;
 1036    if (my ($disallowed) = ($value =~ /([^\[\]\+\-0-9;,.eE ]+)/)) {
 1037       croak("PDL::Core::new_pdl_from_string: found disallowed character(s) '$disallowed' in $original_value");
 1038    }
 1039 
 1040    # Wrap the string in brackets [], so that the following works:
 1041    # $x = new PDL q[1 2 3];
 1042    # We'll have to check for dimensions of size one after we've parsed
 1043    # the string and built a PDL from the resulting array.
 1044    $value = '[' . $value . ']';
 1045 
 1046    # Make sure that each closing bracket followed by an opening bracket
 1047    # has a comma in between them:
 1048    $value =~ s/\]\s*\[/],[/g;
 1049 
 1050    # Semicolons indicate 'start a new row' and require special handling:
 1051    if ($value =~ /;/) {
 1052       $value =~ s/(\[[^\]]+;[^\]]+\])/[$1]/g;
 1053       $value =~ s/;/],[/g;
 1054    }
 1055 
 1056    # Remove ending decimal points and insert zeroes in front of starting
 1057    # decimal points. This makes the white-space-to-comma replacement
 1058    # in the next few lines much simpler.
 1059    $value =~ s/(\d\.)(z|[^\d])/${1}0$2/g;
 1060    $value =~ s/(\A|[^\d])\./${1}0./g;
 1061 
 1062    # Remove whitspace between signs and the numbers that follow them:
 1063    $value =~ s/([+\-])\s+/$1/g;
 1064 
 1065 #   # make unambiguous addition/subtraction (white-space on both sides
 1066 #   # of operator) by removing white-space from both sides
 1067 #   $value =~ s/([\dEe])\s+([+\-])\s+(?=[Ee\d])/$1$2/g;
 1068 
 1069    # Replace white-space separators with commas:
 1070    $value =~ s/([.\deE])\s+(?=[+\-eE\d])/$1,/g;
 1071 
 1072    # Remove all other white space:
 1073    $value =~ s/\s+//g;
 1074 
 1075    # Croak on operations with bad values. It might be nice to simply replace
 1076    # these with bad values, but that is more difficult that I like, so I'm just
 1077    # going to disallow that here:
 1078    croak("PDL::Core::new_pdl_from_string: Operations with bad values are not supported")
 1079       if($value =~ /EE[+\-]/ or $value =~ /[+\-]EE/);
 1080 
 1081    # Check for things that will evaluate as functions and croak if found
 1082    if (my ($disallowed) = ($value =~ /((\D+|\A)[eE]\d+)/)) {
 1083       croak("PDL::Core::new_pdl_from_string: syntax error, looks like an improper exponentiation: $disallowed\n"
 1084          . "You originally gave me $original_value\n");
 1085    }
 1086 
 1087    # Replace the place-holder strings with strings that will evaluate to their
 1088    # correct numerical values when we run the eval:
 1089    $value =~ s/\bEE\b/bad/g;
 1090    my $bad = $types[$type]->badvalue;
 1091    $value =~ s/\bee\b/nan/g;
 1092    my $inf = -pdl(0)->log;
 1093    $value =~ s/\bEe\b/inf/g;
 1094    my $nnan = $inf - $inf;
 1095    my $nan= $this->initialize();
 1096    $nan->set_datatype($nnan->get_datatype);
 1097    $nan->setdims([]);
 1098 
 1099    # pack("d*", "nan") will work here only on perls that numify the string "nan" to a NaN.
 1100    # pack( "d*", (-1.0) ** 0.5 ) will hopefully work in more places, though it seems both
 1101    # pack("d*", "nan") and pack( "d*", (-1.0) ** 0.5 ) fail on *old* MS Compilers (MSVC++ 6.0 and earlier).
 1102    # sisyphus 4 Jan 2013.
 1103    ${$nan->get_dataref}     = pack( "d*", (-1.0) ** 0.5 );
 1104 
 1105    $nan->upd_data();
 1106    $value =~ s/\beE\b/pi/g;
 1107 
 1108    my $val = eval {
 1109       # Install the warnings handler:
 1110       my $old_warn_handler = $SIG{__WARN__};
 1111       local $SIG{__WARN__} = sub {
 1112          if ($_[0] =~ /(Argument ".*" isn't numeric)/) {
 1113             # Send the error through die. This is *always* get caught, so keep
 1114             # it simple.
 1115             die "Incorrectly formatted input: $1\n";
 1116          }
 1117          elsif ($old_warn_handler) {
 1118             $old_warn_handler->(@_);
 1119          }
 1120          else {
 1121             warn @_;
 1122          }
 1123       };
 1124 
 1125       # Let's see if we can parse it as an array-of-arrays:
 1126       local $_ = $value;
 1127       return PDL::Core::parse_basic_string ($inf, $nan, $nnan, $bad);
 1128    };
 1129 
 1130    # Respect BADVAL_USENAN
 1131    require PDL::Config;
 1132    $has_bad += $has_inf + $has_nan if $PDL::Config{BADVAL_USENAN};
 1133 
 1134    if (ref $val eq 'ARRAY') {
 1135       my $to_return = PDL::Core::pdl_avref($val,$this,$type);
 1136       if( $to_return->dim(-1) == 1 ) {
 1137           if( $to_return->dims > 1 ) {
 1138               # remove potentially spurious last dimension
 1139               $to_return = $to_return->mv(-1,1)->clump(2)->sever;
 1140           } elsif( $to_return->dims == 1 ) {
 1141               # fix scalar values
 1142               $to_return->setdims([]);
 1143           }
 1144       }
 1145       # Mark bad if appropriate
 1146       $to_return->badflag($has_bad > 0);
 1147       return $to_return;
 1148    }
 1149    else {
 1150       my @message = ("PDL::Core::new_pdl_from_string: string input='$original_value', string output='$value'" );
 1151       if ($@) {
 1152          push @message, $@;
 1153       } else {
 1154          push @message, "Internal error: unexpected output type ->$val<- is not ARRAY ref";
 1155       }
 1156       croak join("\n  ", @message);
 1157    }
 1158 }
 1159 
 1160 sub PDL::Core::parse_basic_string {
 1161     # Assumes $_ holds the string of interest, and modifies that value
 1162     # in-place.
 1163 
 1164     use warnings;
 1165 
 1166     # Takes a string with proper bracketing, etc, and returns an array-of-arrays
 1167     # filled with numbers, suitable for use with pdl_avref. It uses recursive
 1168     # descent to handle the nested nature of the data. The string should have
 1169     # no whitespace and should be something that would evaluate into a Perl
 1170     # array-of-arrays (except that strings like 'inf', etc, are allowed).
 1171 
 1172     my ($inf, $nan, $nnan, $bad) = @_;
 1173 
 1174     # First character should be a bracket:
 1175     die "Internal error: input string -->$_<-- did not start with an opening bracket\n"
 1176         unless s/^\[//;
 1177 
 1178     my @to_return;
 1179     # Loop until we run into our closing bracket:
 1180     my $sign = 1;
 1181     my $expects_number = 0;
 1182     SYMBOL: until (s/^\]//) {
 1183         # If we have a bracket, then go recursive:
 1184         if (/^\[/) {
 1185             die "Expected a number but found a bracket at ... ", substr ($_, 0, 10), "...\n"
 1186                 if $expects_number;
 1187             push @to_return, PDL::Core::parse_basic_string(@_);
 1188             next SYMBOL;
 1189         }
 1190         elsif (s/^\+//) {
 1191             die "Expected number but found a plus sign at ... ", substr ($_, 0, 10), "...\n"
 1192                 if $expects_number;
 1193             $expects_number = 1;
 1194             redo SYMBOL;
 1195         }
 1196         elsif (s/^\-//) {
 1197             die "Expected number but found a minus sign at ... ", substr ($_, 0, 10), "...\n"
 1198                 if $expects_number;
 1199             $sign = -1;
 1200             $expects_number = 1;
 1201             redo SYMBOL;
 1202         }
 1203         elsif (s/^bad//i) {
 1204             push @to_return, $bad;
 1205         }
 1206         elsif (s/^inf//i or s/1\.\#INF//i) {
 1207             push @to_return, $sign * $inf;
 1208         }
 1209         elsif (s/^nan//i or s/^1\.\#IND//i) {
 1210                         if ($sign == -1) {
 1211                           push @to_return, $nnan;
 1212                         } else {
 1213                           push @to_return, $nan;
 1214                         }
 1215         }
 1216         elsif (s/^pi//i) {
 1217             push @to_return, $sign * 4 * atan2(1, 1);
 1218         }
 1219         elsif (s/^e//i) {
 1220             push @to_return, $sign * exp(1);
 1221         }
 1222         elsif (s/^([\d+\-e.]+)//i) {
 1223             # Note that improper numbers are handled by the warning signal
 1224             # handler
 1225                         my $val = $1;
 1226                         my $nval = $val + 0x0;
 1227                         push @to_return, ($sign>0x0) ? $nval : -$nval;
 1228         }
 1229         else {
 1230             die "Incorrectly formatted input at:\n  ", substr ($_, 0, 10), "...\n";
 1231         }
 1232     }
 1233     # Strip off any commas
 1234     continue {
 1235         $sign = 1;
 1236         $expects_number = 0;
 1237         s/^,//;
 1238     }
 1239 
 1240     return \@to_return;
 1241 }
 1242 
 1243 sub PDL::new {
 1244    # print "in PDL::new\n";
 1245    my $this = shift;
 1246    return $this->copy if ref($this);
 1247    my $type = ref($_[0]) eq 'PDL::Type' ? ${shift @_}[0]  : $PDL_D;
 1248    my $value = (@_ >1 ? [@_] : shift);  # ref thyself
 1249 
 1250    unless(defined $value) {
 1251        if($PDL::debug && $PDL::undefval) {
 1252        print STDERR "Warning: PDL::new converted undef to $PDL::undefval ($PDL::undefval)\n";
 1253        }
 1254        $value = $PDL::undefval+0
 1255    }
 1256 
 1257    return pdl_avref($value,$this,$type) if ref($value) eq "ARRAY";
 1258    my $new = $this->initialize();
 1259    $new->set_datatype($type);
 1260 
 1261 
 1262    if (ref(\$value) eq "SCALAR") {
 1263       # The string processing is extremely slow. Benchmarks indicated that it
 1264       # takes 10x longer to process a scalar number compared with normal Perl
 1265       # conversion of a string to a number. So, only use the string processing
 1266       # if the input looks like a real string, i.e. it doesn't look like a plain
 1267       # number. Note that for our purposes, looks_like_number incorrectly
 1268       # handles the strings 'inf' and 'nan' on Windows machines. We want to send
 1269       # those to the string processing, so this checks for them in a way that
 1270       # short-circuits the looks_like_number check.
 1271       if (PDL::Core::is_scalar_SvPOK($value)
 1272             and ($value =~ /inf/i or $value =~ /nan/i
 1273                or !Scalar::Util::looks_like_number($value))) {
 1274          # new was passed a string argument that doesn't look like a number
 1275          # so we can process as a Matlab-style data entry format.
 1276         return PDL::Core::new_pdl_from_string($new,$value,$this,$type);
 1277       } elsif ($Config{ivsize} < 8 && $pack[$new->get_datatype] eq 'q*') {
 1278          # special case when running on a perl without 64bit int support
 1279          # we have to avoid pack("q", ...) in this case
 1280          # because it dies with error: "Invalid type 'q' in pack"
 1281          $new->setdims([]);
 1282          set_c($new, [0], $value);
 1283       } else {
 1284          $new->setdims([]);
 1285          ${$new->get_dataref}     = pack( $pack[$new->get_datatype], $value );
 1286          $new->upd_data();
 1287       }
 1288    }
 1289    elsif (blessed($value)) { # Object
 1290        $new = $value->copy;
 1291    }
 1292    else {
 1293        barf("Can not interpret argument $value of type ".ref($value) );
 1294    }
 1295    return $new;
 1296 }
 1297 
 1298 
 1299 =head2 copy
 1300 
 1301 =for ref
 1302 
 1303 Make a physical copy of a piddle
 1304 
 1305 =for usage
 1306 
 1307  $new = $old->copy;
 1308 
 1309 Since C<$new = $old> just makes a new reference, the
 1310 C<copy> method is provided to allow real independent
 1311 copies to be made.
 1312 
 1313 =cut
 1314 
 1315 # Inheritable copy method
 1316 #
 1317 # XXX Must be fixed
 1318 # Inplace is handled by the op currently.
 1319 
 1320 sub PDL::copy {
 1321     my $value = shift;
 1322     barf("Argument is an ".ref($value)." not an object") unless blessed($value);
 1323     my $option  = shift;
 1324     $option = "" if !defined $option;
 1325     if ($value->is_inplace) {   # Copy protection
 1326        $value->set_inplace(0);
 1327        return $value;
 1328     }
 1329     # threadI(-1,[]) is just an identity vafftrans with threadId copying ;)
 1330     my $new = $value->threadI(-1,[])->sever;
 1331     return $new;
 1332 }
 1333 
 1334 =head2 hdr_copy
 1335 
 1336 =for ref
 1337 
 1338 Return an explicit copy of the header of a PDL.
 1339 
 1340 hdr_copy is just a wrapper for the internal routine _hdr_copy, which
 1341 takes the hash ref itself.  That is the routine which is used to make
 1342 copies of the header during normal operations if the hdrcpy() flag of
 1343 a PDL is set.
 1344 
 1345 General-purpose deep copies are expensive in perl, so some simple
 1346 optimization happens:
 1347 
 1348 If the header is a tied array or a blessed hash ref with an associated
 1349 method called C<copy>, then that ->copy method is called.  Otherwise, all
 1350 elements of the hash are explicitly copied.  References are recursively
 1351 deep copied.
 1352 
 1353 This routine seems to leak memory.
 1354 
 1355 =cut
 1356 
 1357 sub PDL::hdr_copy {
 1358   my $pdl = shift;
 1359   my $hdr = $pdl->gethdr;
 1360   return PDL::_hdr_copy($hdr);
 1361 }
 1362 
 1363 # Same as hdr_copy but takes a hash ref instead of a PDL.
 1364 sub PDL::_hdr_copy {
 1365   my $hdr = shift;
 1366   my $tobj;
 1367 
 1368   print "called _hdr_copy\n" if($PDL::debug);
 1369 
 1370   unless( (ref $hdr)=~m/HASH/ ) {
 1371     print"returning undef\n" if($PDL::debug);
 1372     return undef ;
 1373   }
 1374 
 1375   if($tobj = tied %$hdr) { #
 1376     print "tied..."if($PDL::debug);
 1377     if(UNIVERSAL::can($tobj,"copy")) {
 1378       my %rhdr;
 1379       tie(%rhdr, ref $tobj, $tobj->copy);
 1380       print "returning\n" if($PDL::debug);
 1381       return \%rhdr;
 1382     }
 1383 
 1384     # Astro::FITS::Header is special for now -- no copy method yet
 1385     # but it is recognized.  Once it gets a copy method this will become
 1386     # vestigial:
 1387 
 1388     if(UNIVERSAL::isa($tobj,"Astro::FITS::Header")) {
 1389       print "Astro::FITS::Header..." if($PDL::debug);
 1390       my @cards = $tobj->cards;
 1391       my %rhdr;
 1392       tie(%rhdr,"Astro::FITS::Header", new Astro::FITS::Header(Cards=>\@cards));
 1393       print "returning\n" if($PDL::debug);
 1394       return \%rhdr;
 1395     }
 1396   }
 1397   elsif(UNIVERSAL::can($hdr,"copy")) {
 1398     print "found a copy method\n" if($PDL::debug);
 1399     return $hdr->copy;
 1400   }
 1401 
 1402   # We got here if it's an unrecognized tie or if it's a vanilla hash.
 1403   print "Making a hash copy..." if($PDL::debug);
 1404 
 1405   return PDL::_deep_hdr_copy($hdr);
 1406 
 1407 }
 1408 
 1409 #
 1410 # Sleazy deep-copier that gets most cases
 1411 # --CED 14-April-2003
 1412 #
 1413 
 1414 sub PDL::_deep_hdr_copy {
 1415   my $val = shift;
 1416 
 1417   if(ref $val eq 'HASH') {
 1418     my (%a,$key);
 1419     for $key(keys %$val) {
 1420       my $value = $val->{$key};
 1421       $a{$key} = (ref $value) ? PDL::_deep_hdr_copy($value) : $value;
 1422     }
 1423     return \%a;
 1424   }
 1425 
 1426   if(ref $val eq 'ARRAY') {
 1427     my (@a,$z);
 1428     for $z(@$val) {
 1429       push(@a,(ref $z) ? PDL::_deep_hdr_copy($z) : $z);
 1430     }
 1431     return \@a;
 1432   }
 1433 
 1434   if(ref $val eq 'SCALAR') {
 1435     my $x = $$val;
 1436     return \$x;
 1437   }
 1438 
 1439   if(ref $val eq 'REF') {
 1440     my $x = PDL::_deep_hdr_copy($$val);
 1441     return \$x;
 1442   }
 1443 
 1444   # Special case for PDLs avoids potential nasty header recursion...
 1445   if(UNIVERSAL::isa($val,'PDL')) {
 1446     my $h;
 1447     $val->hdrcpy(0) if($h = $val->hdrcpy); # assignment
 1448     my $out = $val->copy;
 1449     $val->hdrcpy($h) if($h);
 1450     return $out;
 1451   }
 1452 
 1453   if(UNIVERSAL::can($val,'copy')) {
 1454     return $val->copy;
 1455   }
 1456 
 1457   $val;
 1458 }
 1459 
 1460 
 1461 =head2 unwind
 1462 
 1463 =for ref
 1464 
 1465 Return a piddle which is the same as the argument except
 1466 that all threadids have been removed.
 1467 
 1468 =for usage
 1469 
 1470  $y = $x->unwind;
 1471 
 1472 =head2 make_physical
 1473 
 1474 =for ref
 1475 
 1476 Make sure the data portion of a piddle can be accessed from XS code.
 1477 
 1478 =for example
 1479 
 1480  $x->make_physical;
 1481  $x->call_my_xs_method;
 1482 
 1483 Ensures that a piddle gets its own allocated copy of data. This obviously
 1484 implies that there are certain piddles which do not have their own data.
 1485 These are so called I<virtual> piddles that make use of the I<vaffine>
 1486 optimisation (see L<PDL::Indexing|PDL::Indexing>).
 1487 They do not have their own copy of
 1488 data but instead store only access information to some (or all) of another
 1489 piddle's data.
 1490 
 1491 Note: this function should not be used unless absolutely necessary
 1492 since otherwise memory requirements might be severely increased. Instead
 1493 of writing your own XS code with the need to call C<make_physical> you
 1494 might want to consider using the PDL preprocessor
 1495 (see L<PDL::PP|PDL::PP>)
 1496 which can be used to transparently access virtual piddles without the
 1497 need to physicalise them (though there are exceptions).
 1498 
 1499 =cut
 1500 
 1501 sub PDL::unwind {
 1502     my $value = shift;
 1503     my $foo = $value->null();
 1504     $foo .= $value->unthread();
 1505     return $foo;
 1506 }
 1507 
 1508 =head2 dummy
 1509 
 1510 =for ref
 1511 
 1512 Insert a 'dummy dimension' of given length (defaults to 1)
 1513 
 1514 No relation to the 'Dungeon Dimensions' in Discworld!
 1515 
 1516 Negative positions specify relative to last dimension,
 1517 i.e. C<dummy(-1)> appends one dimension at end,
 1518 C<dummy(-2)> inserts a dummy dimension in front of the
 1519 last dim, etc.
 1520 
 1521 If you specify a dimension position larger than the existing
 1522 dimension list of your PDL, the PDL gets automagically padded with extra
 1523 dummy dimensions so that you get the dim you asked for, in the slot you
 1524 asked for.  This could cause you trouble if, for example,
 1525 you ask for $x->dummy(5000,1) because $x will get 5,000 dimensions,
 1526 each of rank 1.
 1527 
 1528 Because padding at the beginning of the dimension list moves existing
 1529 dimensions from slot to slot, it's considered unsafe, so automagic
 1530 padding doesn't work for large negative indices -- only for large
 1531 positive indices.
 1532 
 1533 =for usage
 1534 
 1535  $y = $x->dummy($position[,$dimsize]);
 1536 
 1537 =for example
 1538 
 1539  pdl> p sequence(3)->dummy(0,3)
 1540  [
 1541   [0 0 0]
 1542   [1 1 1]
 1543   [2 2 2]
 1544  ]
 1545 
 1546  pdl> p sequence(3)->dummy(3,2)
 1547  [
 1548   [
 1549    [0 1 2]
 1550   ]
 1551   [
 1552    [0 1 2]
 1553   ]
 1554  ]
 1555 
 1556  pdl> p sequence(3)->dummy(-3,2)
 1557  Runtime error: PDL: For safety, <pos> < -(dims+1) forbidden in dummy.  min=-2, pos=-3
 1558 
 1559 =cut
 1560 
 1561 sub PDL::dummy($$;$) {
 1562    my ($pdl,$dim,$size) = @_;
 1563    barf("Missing position argument to dummy()") unless defined $dim;  # required argument
 1564    $dim = $pdl->getndims+1+$dim if $dim < 0;
 1565    $size = defined($size) ? (1 * $size) : 1;  # make $size a number (sf feature # 3479009)
 1566 
 1567    barf("For safety, <pos> < -(dims+1) forbidden in dummy.  min="
 1568      . -($pdl->getndims+1).", pos=". ($dim-1-$pdl->getndims) ) if($dim<0);
 1569 
 1570    # Avoid negative repeat count warning that came with 5.21 and later.
 1571    my $dim_diff = $dim - $pdl->getndims;
 1572    my($s) = ',' x ( $dim_diff > 0 ? $pdl->getndims : $dim );
 1573    $s .= '*1,'  x ( $dim_diff > 0 ? $dim_diff : 0 );
 1574    $s .= "*$size";
 1575 
 1576    $pdl->slice($s);
 1577 }
 1578 
 1579 
 1580 ## Cheesy, slow way
 1581 #   while ($dim>$pdl->getndims){
 1582 #     print STDERR "."; flush STDERR;
 1583 #     $pdl = $pdl->dummy($pdl->getndims,1);
 1584 #   }
 1585 #
 1586 #   barf ("too high/low dimension in call to dummy, allowed min/max=0/"
 1587 #    . $_[0]->getndims)
 1588 #     if $dim>$pdl->getndims || $dim < 0;
 1589 #
 1590 #   $_[2] = 1 if ($#_ < 2);
 1591 #   $pdl->slice((','x$dim)."*$_[2]");
 1592 
 1593 =head2 clump
 1594 
 1595 =for ref
 1596 
 1597 "clumps" several dimensions into one large dimension
 1598 
 1599 If called with one argument C<$n> clumps the first C<$n>
 1600 dimensions into one. For example, if C<$x> has dimensions
 1601 C<(5,3,4)> then after
 1602 
 1603 =for example
 1604 
 1605  $y = $x->clump(2);   # Clump 2 first dimensions
 1606 
 1607 the variable C<$y> will have dimensions C<(15,4)>
 1608 and the element C<$y-E<gt>at(7,3)> refers to the element
 1609 C<$x-E<gt>at(1,2,3)>.
 1610 
 1611 Use C<clump(-1)> to flatten a piddle. The method L<flat|PDL::Core/flat>
 1612 is provided as a convenient alias.
 1613 
 1614 Clumping with a negative dimension in general leaves that many
 1615 dimensions behind -- e.g. clump(-2) clumps all of the first few
 1616 dimensions into a single one, leaving a 2-D piddle.
 1617 
 1618 If C<clump> is called with an index list with more than one element
 1619 it is treated as a list of dimensions that should be clumped together
 1620 into one. The resulting
 1621 clumped dim is placed at the position of the lowest index in the list.
 1622 This convention ensures that C<clump> does the expected thing in
 1623 the usual cases. The following example demonstrates typical usage:
 1624 
 1625   $x = sequence 2,3,3,3,5; # 5D piddle
 1626   $c = $x->clump(1..3);    # clump all the dims 1 to 3 into one
 1627   print $c->info;          # resulting 3D piddle has clumped dim at pos 1
 1628  PDL: Double D [2,27,5]
 1629 
 1630 =cut
 1631 
 1632 sub PDL::clump {
 1633   my $ndims = $_[0]->getndims;
 1634   if ($#_ < 2) {
 1635     return &PDL::_clump_int(@_);
 1636   } else {
 1637     my ($this,@dims) = @_;
 1638     my $targd = $ndims-1;
 1639     my @dimmark = (0..$ndims-1);
 1640     barf "too many dimensions" if @dims > $ndims;
 1641     for my $dim (@dims) {
 1642       barf "dimension index $dim larger than greatest dimension"
 1643     if $dim > $ndims-1 ;
 1644       $targd = $dim if $targd > $dim;
 1645       barf "duplicate dimension $dim" if $dimmark[$dim]++ > $dim;
 1646     }
 1647     my $clumped = $this->thread(@dims)->unthread(0)->clump(scalar @dims);
 1648     $clumped = $clumped->mv(0,$targd) if $targd > 0;
 1649     return $clumped;
 1650   }
 1651 }
 1652 
 1653 =head2 thread_define
 1654 
 1655 =for ref
 1656 
 1657 define functions that support threading at the perl level
 1658 
 1659 =for example
 1660 
 1661  thread_define 'tline(a(n);b(n))', over {
 1662   line $_[0], $_[1]; # make line compliant with threading
 1663  };
 1664 
 1665 
 1666 C<thread_define> provides some support for threading (see
 1667 L<PDL::Indexing>) at the perl level. It allows you to do things for
 1668 which you normally would have resorted to PDL::PP (see L<PDL::PP>);
 1669 however, it is most useful to wrap existing perl functions so that the
 1670 new routine supports PDL threading.
 1671 
 1672 C<thread_define> is used to define new I<threading aware>
 1673 functions. Its first argument is a symbolic repesentation of the new
 1674 function to be defined. The string is composed of the name of the new
 1675 function followed by its signature (see L<PDL::Indexing> and L<PDL::PP>)
 1676 in parentheses. The second argument is a subroutine that will be
 1677 called with the slices of the actual runtime arguments as specified by
 1678 its signature. Correct dimension sizes and minimal number of
 1679 dimensions for all arguments will be checked (assuming the rules of
 1680 PDL threading, see L<PDL::Indexing>).
 1681 
 1682 The actual work is done by the C<signature> class which parses the signature
 1683 string, does runtime dimension checks and the routine C<threadover> that
 1684 generates the loop over all appropriate slices of pdl arguments and creates
 1685 pdls as needed.
 1686 
 1687 Similar to C<pp_def> and its C<OtherPars> option it is possible to
 1688 define the new function so that it accepts normal perl args as well as
 1689 piddles. You do this by using the C<NOtherPars> parameter in the
 1690 signature. The number of C<NOtherPars> specified will be passed
 1691 unaltered into the subroutine given as the second argument of
 1692 C<thread_define>. Let's illustrate this with an example:
 1693 
 1694  PDL::thread_define 'triangles(inda();indb();indc()), NOtherPars => 2',
 1695   PDL::over {
 1696     ${$_[3]} .= $_[4].join(',',map {$_->at} @_[0..2]).",-1,\n";
 1697   };
 1698 
 1699 This defines a function C<triangles> that takes 3 piddles as input
 1700 plus 2 arguments which are passed into the routine unaltered. This routine
 1701 is used to collect lists of indices into a perl scalar that is passed by
 1702 reference. Each line is preceded by a prefix passed as C<$_[4]>. Here is
 1703 typical usage:
 1704 
 1705  $txt = '';
 1706  triangles(pdl(1,2,3),pdl(1),pdl(0),\$txt," "x10);
 1707  print $txt;
 1708 
 1709 resulting in the following output
 1710 
 1711  1,1,0,-1,
 1712  2,1,0,-1,
 1713  3,1,0,-1,
 1714 
 1715 which is used in
 1716 L<PDL::Graphics::TriD::VRML|PDL::Graphics::TriD::VRML>
 1717 to generate VRML output.
 1718 
 1719 Currently, this is probably not much more than a POP (proof of principle)
 1720 but is hoped to be useful enough for some real life work.
 1721 
 1722 Check L<PDL::PP|PDL::PP> for the format of the signature. Currently, the
 1723 C<[t]> qualifier and all type qualifiers are ignored.
 1724 
 1725 =cut
 1726 
 1727 sub PDL::over (&) { $_[0] }
 1728 sub PDL::thread_define ($$) {
 1729   require PDL::PP::Signature;
 1730   my ($str,$sub) = @_;
 1731   my $others = 0;
 1732   if ($str =~ s/[,]*\s*NOtherPars\s*=>\s*([0-9]+)\s*[,]*//) {$others = $1}
 1733   barf "invalid string $str" unless $str =~ /\s*([^(]+)\((.+)\)\s*$/x;
 1734   my ($name,$sigstr) = ($1,$2);
 1735   print "defining '$name' with signature '$sigstr' and $others extra args\n"
 1736                           if $PDL::debug;
 1737   my $sig = new PDL::PP::Signature($sigstr);
 1738   my $args = @{$sig->names}; # number of piddle arguments
 1739   barf "no piddle args" if $args == 0;
 1740   $args--;
 1741   # TODO: $sig->dimcheck(@_) + proper creating generation
 1742   my $def = "\@_[0..$args] = map {PDL::Core::topdl(\$_)} \@_[0..$args];\n".
 1743             '$sig->checkdims(@_);
 1744          PDL::threadover($others,@_,$sig->realdims,$sig->creating,$sub)';
 1745   my $package = caller;
 1746   local $^W = 0; # supress the 'not shared' warnings
 1747   print "defining...\nsub $name { $def }\n" if $PDL::debug;
 1748   eval ("package $package; sub $name { $def }");
 1749   barf "error defining $name: $@\n" if $@;
 1750 }
 1751 
 1752 =head2 thread
 1753 
 1754 =for ref
 1755 
 1756 Use explicit threading over specified dimensions (see also L<PDL::Indexing>)
 1757 
 1758 =for usage
 1759 
 1760  $y = $x->thread($dim,[$dim1,...])
 1761 
 1762 =for example
 1763 
 1764  $x = zeroes 3,4,5;
 1765  $y = $x->thread(2,0);
 1766 
 1767 Same as L<PDL::thread1|/PDL::thread1>, i.e. uses thread id 1.
 1768 
 1769 =cut
 1770 
 1771 sub PDL::thread {
 1772     my $var = shift;
 1773     $var->threadI(1,\@_);
 1774 }
 1775 
 1776 =head2 diagonal
 1777 
 1778 =for ref
 1779 
 1780 Returns the multidimensional diagonal over the specified dimensions.
 1781 
 1782 =for usage
 1783 
 1784  $d = $x->diagonal(dim1, dim2,...)
 1785 
 1786 =for example
 1787 
 1788  pdl> $x = zeroes(3,3,3);
 1789  pdl> ($y = $x->diagonal(0,1))++;
 1790  pdl> p $x
 1791  [
 1792   [
 1793    [1 0 0]
 1794    [0 1 0]
 1795    [0 0 1]
 1796   ]
 1797   [
 1798    [1 0 0]
 1799    [0 1 0]
 1800    [0 0 1]
 1801   ]
 1802   [
 1803    [1 0 0]
 1804    [0 1 0]
 1805    [0 0 1]
 1806   ]
 1807  ]
 1808 
 1809 =cut
 1810 
 1811 sub PDL::diagonal {
 1812     my $var = shift;
 1813     $var->diagonalI(\@_);
 1814 }
 1815 
 1816 =head2 thread1
 1817 
 1818 =for ref
 1819 
 1820 Explicit threading over specified dims using thread id 1.
 1821 
 1822 =for usage
 1823 
 1824  $xx = $x->thread1(3,1)
 1825 
 1826 =for example
 1827 
 1828  Wibble
 1829 
 1830 Convenience function interfacing to
 1831 L<PDL::Slices::threadI|PDL::Slices/threadI>.
 1832 
 1833 =cut
 1834 
 1835 sub PDL::thread1 {
 1836     my $var = shift;
 1837     $var->threadI(1,\@_);
 1838 }
 1839 
 1840 =head2 thread2
 1841 
 1842 =for ref
 1843 
 1844 Explicit threading over specified dims using thread id 2.
 1845 
 1846 =for usage
 1847 
 1848  $xx = $x->thread2(3,1)
 1849 
 1850 =for example
 1851 
 1852  Wibble
 1853 
 1854 Convenience function interfacing to
 1855 L<PDL::Slices::threadI|PDL::Slices/threadI>.
 1856 
 1857 =cut
 1858 
 1859 sub PDL::thread2 {
 1860     my $var = shift;
 1861     $var->threadI(2,\@_);
 1862 }
 1863 
 1864 =head2 thread3
 1865 
 1866 =for ref
 1867 
 1868 Explicit threading over specified dims using thread id 3.
 1869 
 1870 =for usage
 1871 
 1872  $xx = $x->thread3(3,1)
 1873 
 1874 =for example
 1875 
 1876  Wibble
 1877 
 1878 Convenience function interfacing to
 1879 L<PDL::Slices::threadI|PDL::Slices/threadI>.
 1880 
 1881 =cut
 1882 
 1883 sub PDL::thread3 {
 1884     my $var = shift;
 1885     $var->threadI(3,\@_);
 1886 }
 1887 
 1888 my %info = (
 1889         D => {
 1890           Name => 'Dimension',
 1891           Sub => \&PDL::Core::dimstr,
 1892          },
 1893         T => {
 1894           Name => 'Type',
 1895           Sub => sub { return $_[0]->type->shortctype; },
 1896          },
 1897         S => {
 1898           Name => 'State',
 1899           Sub => sub { my $state = '';
 1900                    $state .= 'P' if $_[0]->allocated;
 1901                    $state .= 'V' if $_[0]->vaffine &&
 1902                  !$_[0]->allocated; # apparently can be both?
 1903                    $state .= '-' if $state eq '';   # lazy eval
 1904                    $state .= 'C' if $_[0]->anychgd;
 1905                    $state .= 'B' if $_[0]->badflag;
 1906                    $state;
 1907                  },
 1908          },
 1909         F => {
 1910           Name => 'Flow',
 1911           Sub => sub { my $flows = '';
 1912                    $flows = ($_[0]->bflows ? 'b':'') .
 1913                  '~' . ($_[0]->fflows ? 'f':'')
 1914                    if ($_[0]->flows);
 1915                    $flows;
 1916                  },
 1917          },
 1918         M => {
 1919           Name => 'Mem',
 1920           Sub => sub { my ($size,$unit) = ($_[0]->allocated ?
 1921                            $_[0]->nelem*
 1922                       PDL::howbig($_[0]->get_datatype)/1024 : 0, 'KB');
 1923                    if ($size > 0.01*1024) { $size /= 1024;
 1924                             $unit = 'MB' };
 1925                    return sprintf "%6.2f%s",$size,$unit;
 1926                  },
 1927          },
 1928         C => {
 1929           Name => 'Class',
 1930           Sub => sub { ref $_[0] }
 1931          },
 1932         A => {
 1933           Name => 'Address',
 1934           Sub => sub { use Config;
 1935                                my $ivdformat = $Config{ivdformat};
 1936                                $ivdformat =~ s/"//g;
 1937                                sprintf "%$ivdformat", $_[0]->address }
 1938          },
 1939        );
 1940 
 1941 my $allowed = join '',keys %info;
 1942 
 1943 # print the dimension information about a pdl in some appropriate form
 1944 sub dimstr {
 1945   my $this = shift;
 1946 
 1947   my @dims = $this->dims;
 1948   my @ids  = $this->threadids;
 1949   my ($nids,$i) = ($#ids - 1,0);
 1950   my $dstr = 'D ['. join(',',@dims[0..($ids[0]-1)]) .']';
 1951   if ($nids > 0) {
 1952     for $i (1..$nids) {
 1953       $dstr .= " T$i [". join(',',@dims[$ids[$i]..$ids[$i+1]-1]) .']';
 1954     }
 1955   }
 1956   return $dstr;
 1957 }
 1958 
 1959 =head2 sever
 1960 
 1961 =for ref
 1962 
 1963 sever any links of this piddle to parent piddles
 1964 
 1965 In PDL it is possible for a piddle to be just another
 1966 view into another piddle's data. In that case we call
 1967 this piddle a I<virtual piddle> and the original piddle owning
 1968 the data its parent. In other languages these alternate views
 1969 sometimes run by names such as I<alias> or I<smart reference>.
 1970 
 1971 Typical functions that return such piddles are C<slice>, C<xchg>,
 1972 C<index>, etc. Sometimes, however, you would like to separate the
 1973 I<virtual piddle> from its parent's data and just give it a life of
 1974 its own (so that manipulation of its data doesn't change the parent).
 1975 This is simply achieved by using C<sever>. For example,
 1976 
 1977 =for example
 1978 
 1979    $x = $pdl->index(pdl(0,3,7))->sever;
 1980    $x++;       # important: $pdl is not modified!
 1981 
 1982 In many (but not all) circumstances it acts therefore similar to
 1983 L<copy|PDL::Core/copy>.
 1984 However, in general performance is better with C<sever> and secondly,
 1985 C<sever> doesn't lead to futile copying when used on piddles that
 1986 already have their own data. On the other hand, if you really want to make
 1987 sure to work on a copy of a piddle use L<copy|PDL::Core/copy>.
 1988 
 1989    $x = zeroes(20);
 1990    $x->sever;   # NOOP since $x is already its own boss!
 1991 
 1992 Again note: C<sever> I<is not> the same as L<copy|PDL::Core/copy>!
 1993 For example,
 1994 
 1995    $x = zeroes(1); # $x does not have a parent, i.e. it is not a slice etc
 1996    $y = $x->sever; # $y is now pointing to the same piddle as $x
 1997    $y++;
 1998    print $x;
 1999  [1]
 2000 
 2001 but
 2002 
 2003    $x = zeroes(1);
 2004    $y = $x->copy; # $y is now pointing to a new piddle
 2005    $y++;
 2006    print $x;
 2007  [0]
 2008 
 2009 
 2010 =head2 info
 2011 
 2012 =for ref
 2013 
 2014 Return formatted information about a piddle.
 2015 
 2016 =for usage
 2017 
 2018  $x->info($format_string);
 2019 
 2020 =for example
 2021 
 2022  print $x->info("Type: %T Dim: %-15D State: %S");
 2023 
 2024 Returns a string with info about a piddle. Takes an optional
 2025 argument to specify the format of information a la sprintf.
 2026 Format specifiers are in the form C<%E<lt>widthE<gt>E<lt>letterE<gt>>
 2027 where the width is optional and the letter is one of
 2028 
 2029 =over 7
 2030 
 2031 =item T
 2032 
 2033 Type
 2034 
 2035 =item D
 2036 
 2037 Formatted Dimensions
 2038 
 2039 =item F
 2040 
 2041 Dataflow status
 2042 
 2043 =item S
 2044 
 2045 Some internal flags (P=physical,V=Vaffine,C=changed,B=may contain bad data)
 2046 
 2047 =item C
 2048 
 2049 Class of this piddle, i.e. C<ref $pdl>
 2050 
 2051 =item A
 2052 
 2053 Address of the piddle struct as a unique identifier
 2054 
 2055 =item M
 2056 
 2057 Calculated memory consumption of this piddle's data area
 2058 
 2059 =back
 2060 
 2061 =cut
 2062 
 2063 sub PDL::info {
 2064     my ($this,$str) = @_;
 2065     $str = "%C: %T %D" unless defined $str;
 2066     return ref($this)."->null" if $this->isnull;
 2067     my @hash = split /(%[-,0-9]*[.]?[0-9]*\w)/, $str;
 2068     my @args = ();
 2069     my $nstr = '';
 2070     for my $form (@hash) {
 2071     if ($form =~ s/^%([-,0-9]*[.]?[0-9]*)(\w)$/%$1s/) {
 2072         barf "unknown format specifier $2" unless defined $info{$2};
 2073         push @args, &{$info{$2}->{Sub}}($this);
 2074     }
 2075     $nstr .= $form;
 2076     }
 2077     return sprintf $nstr, @args;
 2078 }
 2079 
 2080 =head2 approx
 2081 
 2082 =for ref
 2083 
 2084 test for approximately equal values (relaxed C<==>)
 2085 
 2086 =for example
 2087 
 2088   # ok if all corresponding values in
 2089   # piddles are within 1e-8 of each other
 2090   print "ok\n" if all approx $x, $y, 1e-8;
 2091 
 2092 C<approx> is a relaxed form of the C<==> operator and
 2093 often more appropriate for floating point types (C<float>
 2094 and C<double>).
 2095 
 2096 Usage:
 2097 
 2098 =for usage
 2099 
 2100   $res = approx $x, $y [, $eps]
 2101 
 2102 The optional parameter C<$eps> is remembered across invocations
 2103 and initially set to 1e-6, e.g.
 2104 
 2105   approx $x, $y;         # last $eps used (1e-6 initially)
 2106   approx $x, $y, 1e-10;  # 1e-10
 2107   approx $x, $y;         # also 1e-10
 2108 
 2109 =cut
 2110 
 2111 my $approx = 1e-6;  # a reasonable init value
 2112 sub PDL::approx {
 2113   my ($x,$y,$eps) = @_;
 2114   $eps = $approx unless defined $eps;  # the default eps
 2115   $approx = $eps;    # remember last eps
 2116   # NOTE: ($x-$y)->abs breaks for non-piddle inputs
 2117   return abs($x-$y) < $eps;
 2118 }
 2119 
 2120 =head2 mslice
 2121 
 2122 =for ref
 2123 
 2124 Convenience interface to L<slice|PDL::Slices/slice>,
 2125 allowing easier inclusion of dimensions in perl code.
 2126 
 2127 =for usage
 2128 
 2129  $w = $x->mslice(...);
 2130 
 2131 =for example
 2132 
 2133  # below is the same as $x->slice("5:7,:,3:4:2")
 2134  $w = $x->mslice([5,7],X,[3,4,2]);
 2135 
 2136 =cut
 2137 
 2138 # called for colon-less args
 2139 # preserves parens if present
 2140 sub intpars { $_[0] =~ /\(.*\)/ ? '('.int($_[0]).')' : int $_[0] }
 2141 
 2142 sub PDL::mslice {
 2143         my($pdl) = shift;
 2144         return $pdl->slice(join ',',(map {
 2145                         !ref $_ && $_ eq "X" ? ":" :
 2146                ref $_ eq "ARRAY" ? $#$_ > 1 && @$_[2] == 0 ?
 2147                "(".int(@$_[0]).")" : join ':', map {int $_} @$_ :
 2148                         !ref $_ ? intpars $_ :
 2149                         die "INVALID SLICE DEF $_"
 2150                 } @_));
 2151 }
 2152 
 2153 =head2 nslice_if_pdl
 2154 
 2155 =for ref
 2156 
 2157 If C<$self> is a PDL, then calls C<slice> with all but the last
 2158 argument, otherwise $self->($_[-1]) is called where $_[-1} is the
 2159 original argument string found during PDL::NiceSlice filtering.
 2160 
 2161 DEVELOPER'S NOTE: this routine is found in Core.pm.PL but would be
 2162 better placed in Slices/slices.pd.  It is likely to be moved there
 2163 and/or changed to "slice_if_pdl" for PDL 3.0.
 2164 
 2165 =for usage
 2166 
 2167  $w = $x->nslice_if_pdl(...,'(args)');
 2168 
 2169 =cut
 2170 
 2171 sub PDL::nslice_if_pdl {
 2172    my ($pdl) = shift;
 2173    my ($orig_args) = pop;
 2174 
 2175    # warn "PDL::nslice_if_pdl called with (@_) args, originally ($orig_args)\n";
 2176 
 2177    if (ref($pdl) eq 'CODE') {
 2178       # barf('PDL::nslice_if_pdl tried to process a sub ref, please use &$subref() syntax')
 2179       @_ = eval $orig_args;
 2180       goto &$pdl;
 2181    }
 2182 
 2183    unshift @_, $pdl;
 2184    goto &PDL::slice;
 2185 }
 2186 
 2187 =head2 nslice
 2188 
 2189 =for ref
 2190 
 2191 C<nslice> was an internally used interface for L<PDL::NiceSlice|PDL::NiceSlice>,
 2192 but is now merely a springboard to L<PDL::Slice|PDL::Slice>.  It is deprecated
 2193 and likely to disappear in PDL 3.0.
 2194 
 2195 =cut
 2196 sub PDL::nslice {
 2197     unless($PDL::nslice_warning_issued) {
 2198     $PDL::nslice_warning_issued = 1;
 2199     warn "WARNING: deprecated call to PDL::nslice detected.  Use PDL::slice instead.\n (Warning will be issued only once per session)\n";
 2200     }
 2201     goto &PDL::slice;
 2202 }
 2203 
 2204 sub blessed {
 2205     my $ref = ref(shift);
 2206     return $ref =~ /^(REF|SCALAR|ARRAY|HASH|CODE|GLOB||)$/ ? 0 : 1;
 2207 }
 2208 
 2209 # Convert numbers to PDL if not already
 2210 
 2211 sub PDL::topdl {
 2212     return $_[0]->new(@_[1..$#_]) if($#_ > 1); # PDLify an ARRAY
 2213     return $_[1] if blessed($_[1]); # Fall through
 2214     return $_[0]->new($_[1]) if ref(\$_[1]) eq  'SCALAR' or
 2215            ref($_[1]) eq 'ARRAY';
 2216     barf("Can not convert a ".ref($_[1])." to a ".$_[0]);
 2217 0;}
 2218 
 2219 # Convert everything to PDL if not blessed
 2220 
 2221 sub alltopdl {
 2222     if (ref $_[2] eq 'PDL::Type') {
 2223       return convert($_[1], $_[2]) if blessed($_[1]);
 2224       return $_[0]->new($_[2], $_[1]) if $_[0] eq 'PDL';
 2225     }
 2226     return $_[1] if blessed($_[1]); # Fall through
 2227     return $_[0]->new($_[1]);
 2228 0;}
 2229 
 2230 
 2231 =head2 inplace
 2232 
 2233 =for ref
 2234 
 2235 Flag a piddle so that the next operation is done 'in place'
 2236 
 2237 =for usage
 2238 
 2239  somefunc($x->inplace); somefunc(inplace $x);
 2240 
 2241 In most cases one likes to use the syntax C<$y = f($x)>, however
 2242 in many case the operation C<f()> can be done correctly
 2243 'in place', i.e. without making a new copy of the data for
 2244 output. To make it easy to use this, we write C<f()> in such
 2245 a way that it operates in-place, and use C<inplace> to hint
 2246 that a new copy should be disabled. This also makes for
 2247 clear syntax.
 2248 
 2249 Obviously this will not work for all functions, and if in
 2250 doubt see the function's documentation. However one
 2251 can assume this is
 2252 true for all elemental functions (i.e. those which just
 2253 operate array element by array element like C<log10>).
 2254 
 2255 =for example
 2256 
 2257  pdl> $x = xvals zeroes 10;
 2258  pdl> log10(inplace $x)
 2259  pdl> p $x
 2260  [-inf 0    0.30103 0.47712125 0.60205999    0.69897 0.77815125 0.84509804 0.90308999 0.95424251]
 2261 
 2262 =cut
 2263 
 2264 # Flag pdl for in-place operations
 2265 
 2266 sub PDL::inplace {
 2267     my $pdl = PDL->topdl(shift); $pdl->set_inplace(1); return $pdl;
 2268 }
 2269 
 2270 # Copy if not inplace
 2271 
 2272 
 2273 =head2 is_inplace
 2274 
 2275 =for ref
 2276 
 2277 Test the in-place flag on a piddle
 2278 
 2279 =for usage
 2280 
 2281   $out = ($in->is_inplace) ? $in : zeroes($in);
 2282   $in->set_inplace(0)
 2283 
 2284 Provides access to the L<inplace|/inplace> hint flag, within the perl millieu.
 2285 That way functions you write can be inplace aware... If given an
 2286 argument the inplace flag will be set or unset depending on the value
 2287 at the same time. Can be used for shortcut tests that delete the
 2288 inplace flag while testing:
 2289 
 2290   $out = ($in->is_inplace(0)) ? $in : zeroes($in); # test & unset!
 2291 
 2292 =head2 set_inplace
 2293 
 2294 =for ref
 2295 
 2296 Set the in-place flag on a piddle
 2297 
 2298 =for usage
 2299 
 2300   $out = ($in->is_inplace) ? $in : zeroes($in);
 2301   $in->set_inplace(0);
 2302 
 2303 Provides access to the L<inplace|/inplace> hint flag, within the perl millieu.
 2304 Useful mainly for turning it OFF, as L<inplace|/inplace> turns it ON more
 2305 conveniently.
 2306 
 2307 =head2 new_or_inplace
 2308 
 2309 =for usage
 2310 
 2311     $w = new_or_inplace(shift());
 2312     $w = new_or_inplace(shift(),$preferred_type);
 2313 
 2314 =for ref
 2315 
 2316 Return back either the argument pdl or a copy of it depending on whether
 2317 it be flagged in-place or no.  Handy for building inplace-aware functions.
 2318 
 2319 If you specify a preferred type (must be one of the usual PDL type strings,
 2320 a list ref containing several of them, or a string containing several of them),
 2321 then the copy is coerced into the first preferred type listed if it is not
 2322 already one of the preferred types.
 2323 
 2324 Note that if the inplace flag is set, no coersion happens even if you specify
 2325 a preferred type.
 2326 
 2327 =cut
 2328 
 2329 sub new_or_inplace {
 2330     my $pdl = shift;
 2331     my $preferred = shift;
 2332     my $force = shift;
 2333     if($pdl->is_inplace) {
 2334         $pdl->set_inplace(0);
 2335         return $pdl;
 2336     } else {
 2337         unless(defined($preferred)) {
 2338         return $pdl->copy;
 2339         } else {
 2340         $preferred = join(",",@$preferred) if(ref($preferred) eq 'ARRAY');
 2341         my $s = "".$pdl->type;
 2342         if($preferred =~ m/(^|\,)$s(\,|$)/i) {
 2343             # Got a match - the PDL is one of the preferred types.
 2344             return $pdl->copy();
 2345         } else {
 2346             # No match - promote it to the first in the list.
 2347             $preferred =~ s/\,.*//;
 2348             my $out = PDL::new_from_specification('PDL',new PDL::Type($preferred),$pdl->dims);
 2349             $out .= $pdl;
 2350             return $out;
 2351         }
 2352         }
 2353     }
 2354     barf "PDL::Core::new_or_inplace - This can never happen!";
 2355 }
 2356 *PDL::new_or_inplace = \&new_or_inplace;
 2357 
 2358 # Allow specifications like zeroes(10,10) or zeroes($x)
 2359 # or zeroes(inplace $x) or zeroes(float,4,3)
 2360 
 2361 =head2 new_from_specification
 2362 
 2363 =for ref
 2364 
 2365 Internal method: create piddle by specification
 2366 
 2367 This is the argument processing method called by L<zeroes|/zeroes>
 2368 and some other functions
 2369 which constructs piddles from argument lists of the form:
 2370 
 2371  [type], $nx, $ny, $nz,...
 2372 
 2373 For C<$nx>, C<$ny>, etc. 0 and 1D piddles are allowed.
 2374 Giving those has the same effect as if saying C<$arg-E<gt>list>,
 2375 e.g.
 2376 
 2377    1, pdl(5,2), 4
 2378 
 2379 is equivalent to
 2380 
 2381    1, 5, 2, 4
 2382 
 2383 Note, however, that in all functions using C<new_from_specification>
 2384 calling C<func $piddle> will probably not do what you want. So to play safe
 2385 use (e.g. with zeroes)
 2386 
 2387   $pdl = zeroes $dimpdl->list;
 2388 
 2389 Calling
 2390 
 2391   $pdl = zeroes $dimpdl;
 2392 
 2393 will rather be equivalent to
 2394 
 2395   $pdl = zeroes $dimpdl->dims;
 2396 
 2397 However,
 2398 
 2399   $pdl = zeroes ushort, $dimpdl;
 2400 
 2401 will again do what you intended since it is interpreted
 2402 as if you had said
 2403 
 2404   $pdl = zeroes ushort, $dimpdl->list;
 2405 
 2406 This is unfortunate and confusing but no good solution seems
 2407 obvious that would not break existing scripts.
 2408 
 2409 =cut
 2410 
 2411 sub PDL::new_from_specification{
 2412     my $class = shift;
 2413     my $type = ref($_[0]) eq 'PDL::Type' ? ${shift @_}[0]  : $PDL_D;
 2414     my $nelems = 1; my @dims;
 2415     for (@_) {
 2416        if (ref $_) {
 2417          barf "Trying to use non-piddle as dimensions?" unless $_->isa('PDL');
 2418          barf "Trying to use multi-dim piddle as dimensions?"
 2419               if $_->getndims > 1;
 2420          warn "creating > 10 dim piddle (piddle arg)!"
 2421               if $_->nelem > 10;
 2422          for my $dim ($_->list) {$nelems *= $dim; push @dims, $dim}
 2423        } else {
 2424           if ($_) {  # quiet warnings when $_ is the empty string
 2425              barf "Dimensions must be non-negative" if $_<0;
 2426              $nelems *= $_; push @dims, $_
 2427           } else {
 2428              $nelems *= 0; push @dims, 0;
 2429           }
 2430        }
 2431     }
 2432     my $pdl = $class->initialize();
 2433     $pdl->set_datatype($type);
 2434     $pdl->setdims([@dims]);
 2435     print "Dims: ",(join ',',@dims)," DLen: ",(length $ {$pdl->get_dataref}),"\n" if $PDL::debug;
 2436     return $pdl;
 2437 }
 2438 
 2439 =head2 isnull
 2440 
 2441 =for ref
 2442 
 2443 Test whether a piddle is null
 2444 
 2445 =for usage
 2446 
 2447  croak("Input piddle mustn't be null!")
 2448      if $input_piddle->isnull;
 2449 
 2450 This function returns 1 if the piddle is null, zero if it is not. The purpose
 2451 of null piddles is to "tell" any PDL::PP methods to allocate new memory for
 2452 an output piddle, but only when that PDL::PP method is called in full-arg
 2453 form. Of course, there's no reason you couldn't commandeer the special value
 2454 for your own purposes, for which this test function would prove most helpful.
 2455 But in general, you shouldn't need to test for a piddle's nullness.
 2456 
 2457 See L</Null PDLs> for more information.
 2458 
 2459 =head2 isempty
 2460 
 2461 =for ref
 2462 
 2463 Test whether a piddle is empty
 2464 
 2465 =for usage
 2466 
 2467  print "The piddle has zero dimension\n" if $pdl->isempty;
 2468 
 2469 This function returns 1 if the piddle has zero elements. This is
 2470 useful in particular when using the indexing function which. In the
 2471 case of no match to a specified criterion, the returned piddle has
 2472 zero dimension.
 2473 
 2474  pdl> $w=sequence(10)
 2475  pdl> $i=which($w < -1)
 2476  pdl> print "I found no matches!\n" if ($i->isempty);
 2477  I found no matches!
 2478 
 2479 Note that having zero elements is rather different from the concept
 2480 of being a null piddle, see the L<PDL::FAQ|PDL::FAQ> and
 2481 L<PDL::Indexing|PDL::Indexing>
 2482 manpages for discussions of this.
 2483 
 2484 =cut
 2485 
 2486 sub PDL::isempty {
 2487     my $pdl=shift;
 2488     return ($pdl->nelem == 0);
 2489 }
 2490 
 2491 =head2 zeroes
 2492 
 2493 =for ref
 2494 
 2495 construct a zero filled piddle from dimension list or template piddle.
 2496 
 2497 Various forms of usage,
 2498 
 2499 (i) by specification or (ii) by template piddle:
 2500 
 2501 =for usage
 2502 
 2503  # usage type (i):
 2504  $w = zeroes([type], $nx, $ny, $nz,...);
 2505  $w = PDL->zeroes([type], $nx, $ny, $nz,...);
 2506  $w = $pdl->zeroes([type], $nx, $ny, $nz,...);
 2507  # usage type (ii):
 2508  $w = zeroes $y;
 2509  $w = $y->zeroes
 2510  zeroes inplace $w;     # Equivalent to   $w .= 0;
 2511  $w->inplace->zeroes;   #  ""
 2512 
 2513 =for example
 2514 
 2515  pdl> $z = zeroes 4,3
 2516  pdl> p $z
 2517  [
 2518   [0 0 0 0]
 2519   [0 0 0 0]
 2520   [0 0 0 0]
 2521  ]
 2522  pdl> $z = zeroes ushort, 3,2 # Create ushort array
 2523  [ushort() etc. with no arg returns a PDL::Types token]
 2524 
 2525 See also L<new_from_specification|/PDL::new_from_specification>
 2526 for details on using piddles in the dimensions list.
 2527 
 2528 =cut
 2529 
 2530 sub zeroes { ref($_[0]) && ref($_[0]) ne 'PDL::Type' ? PDL::zeroes($_[0]) : PDL->zeroes(@_) }
 2531 sub PDL::zeroes {
 2532     my $class = shift;
 2533     my $pdl = scalar(@_)? $class->new_from_specification(@_) : $class->new_or_inplace;
 2534     $pdl.=0;
 2535     return $pdl;
 2536 }
 2537 
 2538 # Create convenience aliases for zeroes
 2539 
 2540 =head2 zeros
 2541 
 2542 =for ref
 2543 
 2544 construct a zero filled piddle (see zeroes for usage)
 2545 
 2546 =cut
 2547 
 2548 *zeros = \&zeroes;
 2549 *PDL::zeros = \&PDL::zeroes;
 2550 
 2551 =head2 ones
 2552 
 2553 =for ref
 2554 
 2555 construct a one filled piddle
 2556 
 2557 =for usage
 2558 
 2559  $w = ones([type], $nx, $ny, $nz,...);
 2560  etc. (see 'zeroes')
 2561 
 2562 =for example
 2563 
 2564  see zeroes() and add one
 2565 
 2566 See also L<new_from_specification|/PDL::new_from_specification>
 2567 for details on using piddles in the dimensions list.
 2568 
 2569 =cut
 2570 
 2571 sub ones { ref($_[0]) && ref($_[0]) ne 'PDL::Type' ? PDL::ones($_[0]) : PDL->ones(@_) }
 2572 sub PDL::ones {
 2573     my $class = shift;
 2574     my $pdl = scalar(@_)? $class->new_from_specification(@_) : $class->new_or_inplace;
 2575     $pdl.=1;
 2576     return $pdl;
 2577 }
 2578 
 2579 =head2 reshape
 2580 
 2581 =for ref
 2582 
 2583 Change the shape (i.e. dimensions) of a piddle, preserving contents.
 2584 
 2585 =for usage
 2586 
 2587  $x->reshape(NEWDIMS); reshape($x, NEWDIMS);
 2588 
 2589 The data elements are preserved, obviously they will wrap
 2590 differently and get truncated if the new array is shorter.
 2591 If the new array is longer it will be zero-padded.
 2592 
 2593 ***Potential incompatibility with earlier versions of PDL****
 2594 If the list of C<NEWDIMS> is empty C<reshape> will just drop
 2595 all dimensions of size 1 (preserving the number of elements):
 2596 
 2597   $w = sequence(3,4,5);
 2598   $y = $w(1,3);
 2599   $y->reshape();
 2600   print $y->info;
 2601  PDL: Double D [5]
 2602 
 2603 Dimensions of size 1 will also be dropped if C<reshape> is
 2604 invoked with the argument -1:
 2605 
 2606   $y = $w->reshape(-1);
 2607 
 2608 As opposed to C<reshape> without arguments, C<reshape(-1)>
 2609 preserves dataflow:
 2610 
 2611   $w = ones(2,1,2);
 2612   $y = $w(0)->reshape(-1);
 2613   $y++;
 2614   print $w;
 2615  [
 2616   [
 2617    [2 1]
 2618   ]
 2619   [
 2620    [2 1]
 2621   ]
 2622  ]
 2623 
 2624 Important: Piddles are changed inplace!  
 2625 
 2626 Note: If C<$x> is connected to any other PDL (e.g. if it is a slice)
 2627 then the connection is first severed.
 2628 
 2629 =for example
 2630 
 2631  pdl> $x = sequence(10)
 2632  pdl> reshape $x,3,4; p $x
 2633  [
 2634   [0 1 2]
 2635   [3 4 5]
 2636   [6 7 8]
 2637   [9 0 0]
 2638  ]
 2639  pdl> reshape $x,5; p $x
 2640  [0 1 2 3 4]
 2641 
 2642 =cut
 2643 
 2644 *reshape = \&PDL::reshape;
 2645 sub PDL::reshape{
 2646     if (@_ == 2 && $_[1] == -1) {  # a slicing reshape that drops 1-dims
 2647     return $_[0]->slice( map { $_==1 ? [0,0,0] : [] } $_[0]->dims);
 2648     }
 2649     my $pdl = topdl($_[0]);
 2650     $pdl->sever;
 2651     my $nelem = $pdl->nelem;
 2652     my @dims = grep defined, @_[1..$#_];
 2653     for my $dim(@dims) { barf "reshape: invalid dim size '$dim'" if $dim < 0 }
 2654     @dims = grep($_ != 1, $pdl->dims) if @dims == 0; # get rid of dims of size 1
 2655     $pdl->setdims([@dims]);
 2656     $pdl->upd_data;
 2657     if ($pdl->nelem > $nelem) {
 2658     my $tmp=$pdl->clump(-1)->slice("$nelem:-1");
 2659     $tmp .= 0;
 2660     }
 2661     $_[0] = $pdl;
 2662     return $pdl;
 2663 }
 2664 
 2665 =head2 squeeze
 2666 
 2667 =for ref
 2668 
 2669 eliminate all singleton dimensions (dims of size 1)
 2670 
 2671 =for example
 2672 
 2673  $y = $w(0,0)->squeeze;
 2674 
 2675 Alias for C<reshape(-1)>. Removes all singleton dimensions
 2676 and preserves dataflow. A more concise interface is
 2677 provided by L<PDL::NiceSlice|PDL::NiceSlice> via modifiers:
 2678 
 2679  use PDL::NiceSlice;
 2680  $y = $w(0,0;-); # same as $w(0,0)->squeeze
 2681 
 2682 =cut
 2683 
 2684 *squeeze = \&PDL::squeeze;
 2685 sub PDL::squeeze { return $_[0]->reshape(-1) }
 2686 
 2687 =head2 flat
 2688 
 2689 =for ref
 2690 
 2691 flatten a piddle (alias for C<< $pdl->clump(-1) >>)
 2692 
 2693 =for example
 2694 
 2695   $srt = $pdl->flat->qsort;
 2696 
 2697 Useful method to make a 1D piddle from an
 2698 arbitrarily sized input piddle. Data flows
 2699 back and forth as usual with slicing routines.
 2700 Falls through if argument already E<lt>= 1D.
 2701 
 2702 =cut
 2703 
 2704 *flat = \&PDL::flat;
 2705 sub PDL::flat { # fall through if < 2D
 2706   return my $dummy = $_[0]->getndims != 1 ? $_[0]->clump(-1) : $_[0];
 2707 }
 2708 
 2709 =head2 convert
 2710 
 2711 =for ref
 2712 
 2713 Generic datatype conversion function
 2714 
 2715 =for usage
 2716 
 2717  $y = convert($x, $newtypenum);
 2718 
 2719 =for example
 2720 
 2721  $y = convert $x, long
 2722  $y = convert $x, ushort
 2723 
 2724 C<$newtype> is a type B<number>, for convenience they are
 2725 returned by C<long()> etc when called without arguments.
 2726 
 2727 =cut
 2728 
 2729 # type to type conversion functions (with automatic conversion to pdl vars)
 2730 
 2731 sub PDL::convert {
 2732   # we don't allow inplace conversion at the moment
 2733   # (not sure what needs to be changed)
 2734   barf 'Usage: $y = convert($x, $newtypenum)'."\n" if $#_!=1;
 2735   my ($pdl,$type)= @_;
 2736   $pdl = pdl($pdl) unless ref $pdl; # Allow normal numbers
 2737   $type = $type->enum if ref($type) eq 'PDL::Type';
 2738   barf 'Usage: $y = convert($x, $newtypenum)'."\n" unless Scalar::Util::looks_like_number($type);
 2739   return $pdl if $pdl->get_datatype == $type;
 2740   # make_physical-call: temporary stopgap to work around core bug
 2741   my $conv = $pdl->flowconvert($type)->make_physical->sever;
 2742   return $conv;
 2743 }
 2744 
 2745 =head2 Datatype_conversions
 2746 
 2747 =for ref
 2748 
 2749 byte|short|ushort|long|indx|longlong|float|double (shorthands to convert datatypes)
 2750 
 2751 =for usage
 2752 
 2753  $y = double $x; $y = ushort [1..10];
 2754  # all of the above listed shorthands behave similarly
 2755 
 2756 When called with a piddle argument, they convert to the specific
 2757 datatype.
 2758 
 2759 When called with a numeric, list, listref, or string argument they
 2760 construct a new piddle. This is a convenience to avoid having to be
 2761 long-winded and say C<$x = long(pdl(42))>
 2762 
 2763 Thus one can say:
 2764 
 2765  $w = float(1,2,3,4);           # 1D
 2766  $w = float q[1 2 3; 4 5 6];    # 2D
 2767  $w = float([1,2,3],[4,5,6]);   # 2D
 2768  $w = float([[1,2,3],[4,5,6]]); # 2D
 2769 
 2770 Note the last three give identical results, and the last two are exactly
 2771 equivalent - a list is automatically converted to a list reference for
 2772 syntactic convenience. i.e. you can omit the outer C<[]>
 2773 
 2774 When called with no arguments, these functions return a special type token.
 2775 This allows syntactical sugar like:
 2776 
 2777  $x = ones byte, 1000,1000;
 2778 
 2779 This example creates a large piddle directly as byte datatype in
 2780 order to save memory.
 2781 
 2782 In order to control how undefs are handled in converting from perl lists to
 2783 PDLs, one can set the variable C<$PDL::undefval>;
 2784 see the function L<pdl()|/pdl> for more details.
 2785 
 2786 =for example
 2787 
 2788  pdl> p $x=sqrt float [1..10]
 2789  [1 1.41421 1.73205 2 2.23607 2.44949 2.64575 2.82843 3 3.16228]
 2790  pdl> p byte $x
 2791  [1 1 1 2 2 2 2 2 3 3]
 2792 
 2793 =head2 byte
 2794 
 2795 Convert to byte datatype
 2796 
 2797 =head2 short
 2798 
 2799 Convert to short datatype
 2800 
 2801 =head2 ushort
 2802 
 2803 Convert to ushort datatype
 2804 
 2805 =head2 long
 2806 
 2807 Convert to long datatype
 2808 
 2809 =head2 indx
 2810 
 2811 Convert to indx datatype
 2812 
 2813 =head2 longlong
 2814 
 2815 Convert to longlong datatype
 2816 
 2817 =head2 float
 2818 
 2819 Convert to float datatype
 2820 
 2821 =head2 double
 2822 
 2823 Convert to double datatype
 2824 
 2825 =head2 type
 2826 
 2827 =for ref
 2828 
 2829 return the type of a piddle as a blessed type object
 2830 
 2831 A convenience function for use with the piddle constructors, e.g.
 2832 
 2833 =for example
 2834 
 2835  $y = PDL->zeroes($x->type,$x->dims,3);
 2836  die "must be float" unless $x->type == float;
 2837 
 2838 See also the discussion of the C<PDL::Type> class in L<PDL::Types>.
 2839 Note that the C<PDL::Type> objects have overloaded comparison and
 2840 stringify operators so that you can compare and print types:
 2841 
 2842  $x = $x->float if $x->type < float;
 2843  $t = $x->type; print "Type is $t\";
 2844 
 2845 =cut
 2846 
 2847 sub PDL::type { return PDL::Type->new($_[0]->get_datatype); }
 2848 
 2849 ##################### Printing ####################
 2850 
 2851 # New string routine
 2852 
 2853 $PDL::_STRINGIZING = 0;
 2854 
 2855 sub PDL::string {
 2856     my($self,$format)=@_;
 2857     my $to_return = eval {
 2858         if($PDL::_STRINGIZING) {
 2859             return "ALREADY_STRINGIZING_NO_LOOPS";
 2860         }
 2861         local $PDL::_STRINGIZING = 1;
 2862         my $ndims = $self->getndims;
 2863         if($self->nelem > $PDL::toolongtoprint) {
 2864             return "TOO LONG TO PRINT";
 2865         }
 2866         if ($ndims==0) {
 2867         if ( $self->badflag() and $self->isbad() ) {
 2868             return "BAD";
 2869         } else {
 2870             my @x = $self->at();
 2871             return ($format ? sprintf($format, $x[0]) : "$x[0]");
 2872         }
 2873         }
 2874         return "Null" if $self->isnull;
 2875         return "Empty[".join("x",$self->dims)."]" if $self->isempty; # Empty piddle
 2876         local $sep  = $PDL::use_commas ? "," : " ";
 2877         local $sep2 = $PDL::use_commas ? "," : "";
 2878         if ($ndims==1) {
 2879            return str1D($self,$format);
 2880         }
 2881         else{
 2882            return strND($self,$format,0);
 2883         }
 2884     };
 2885     if ($@) {
 2886         # Remove reference to this line:
 2887         $@ =~ s/\s*at .* line \d+\s*\.\n*/./;
 2888         PDL::Core::barf("Stringizing problem: $@");
 2889     }
 2890     return $to_return;
 2891 }
 2892 
 2893 ############## Section/subsection functions ###################
 2894 
 2895 =head2 list
 2896 
 2897 =for ref
 2898 
 2899 Convert piddle to perl list
 2900 
 2901 =for usage
 2902 
 2903  @tmp = list $x;
 2904 
 2905 Obviously this is grossly inefficient for the large datasets PDL is designed to
 2906 handle. This was provided as a get out while PDL matured. It  should now be mostly
 2907 superseded by superior constructs, such as PP/threading. However it is still
 2908 occasionally useful and is provied for backwards compatibility.
 2909 
 2910 =for example
 2911 
 2912  for (list $x) {
 2913    # Do something on each value...
 2914  }
 2915 
 2916 If you compile PDL with bad value support (the default), your machine's
 2917 docs will also say this:
 2918 
 2919 =for bad
 2920 
 2921 list converts any bad values into the string 'BAD'.
 2922 
 2923 =cut
 2924 
 2925 # No threading, just the ordinary dims.
 2926 sub PDL::list{ # pdl -> @list
 2927      barf 'Usage: list($pdl)' if $#_!=0;
 2928      my $pdl = PDL->topdl(shift);
 2929      return () if nelem($pdl)==0;
 2930      @{listref_c($pdl)};
 2931 }
 2932 
 2933 =head2 unpdl
 2934 
 2935 =for ref
 2936 
 2937 Convert piddle to nested Perl array references
 2938 
 2939 =for usage
 2940 
 2941  $arrayref = unpdl $x;
 2942 
 2943 This function returns a reference to a Perl list-of-lists structure
 2944 equivalent to the input piddle (within the limitation that while values
 2945 of elements should be preserved, the detailed datatypes will not as
 2946 perl itself basically has "number" data rather than byte, short, int...
 2947 E.g., C<< sum($x - pdl( $x->unpdl )) >> should equal 0.
 2948 
 2949 Obviously this is grossly inefficient in memory and processing for the
 2950 large datasets PDL is designed to handle. Sometimes, however, you really
 2951 want to move your data back to Perl, and with proper dimensionality,
 2952 unlike C<list>.
 2953 
 2954 =for example
 2955 
 2956  use JSON;
 2957  my $json = encode_json unpdl $pdl;
 2958 
 2959 If you compile PDL with bad value support (the default), your machine's
 2960 docs will also say this:
 2961 
 2962 =cut
 2963 
 2964 =for bad
 2965 
 2966 unpdl converts any bad values into the string 'BAD'.
 2967 
 2968 =cut
 2969 
 2970 sub PDL::unpdl {
 2971     barf 'Usage: unpdl($pdl)' if $#_ != 0;
 2972     my $pdl = PDL->topdl(shift);
 2973     return [] if $pdl->nelem == 0;
 2974     return _unpdl_int($pdl);
 2975 }
 2976 
 2977 sub _unpdl_int {
 2978     my $pdl = shift;
 2979     if ($pdl->ndims > 1) {
 2980         return [ map { _unpdl_int($_) } dog $pdl ];
 2981     } else {
 2982         return listref_c($pdl);
 2983     }
 2984 }
 2985 
 2986 =head2 listindices
 2987 
 2988 =for ref
 2989 
 2990 Convert piddle indices to perl list
 2991 
 2992 =for usage
 2993 
 2994  @tmp = listindices $x;
 2995 
 2996 C<@tmp> now contains the values C<0..nelem($x)>.
 2997 
 2998 Obviously this is grossly inefficient for the large datasets PDL is designed to
 2999 handle. This was provided as a get out while PDL matured. It  should now be mostly
 3000 superseded by superior constructs, such as PP/threading. However it is still
 3001 occasionally useful and is provied for backwards compatibility.
 3002 
 3003 =for example
 3004 
 3005  for $i (listindices $x) {
 3006    # Do something on each value...
 3007  }
 3008 
 3009 =cut
 3010 
 3011 sub PDL::listindices{ # Return list of index values for 1D pdl
 3012      barf 'Usage: list($pdl)' if $#_!=0;
 3013      my $pdl = shift;
 3014      return () if nelem($pdl)==0;
 3015      barf 'Not 1D' if scalar(dims($pdl)) != 1;
 3016      return (0..nelem($pdl)-1);
 3017 }
 3018 
 3019 =head2 set
 3020 
 3021 =for ref
 3022 
 3023 Set a single value inside a piddle
 3024 
 3025 =for usage
 3026 
 3027  set $piddle, @position, $value
 3028 
 3029 C<@position> is a coordinate list, of size equal to the
 3030 number of dimensions in the piddle. Occasionally useful,
 3031 mainly provided for backwards compatibility as superseded
 3032 by use of L<slice|PDL::Slices/slice> and assignment operator C<.=>.
 3033 
 3034 =for example
 3035 
 3036  pdl> $x = sequence 3,4
 3037  pdl> set $x, 2,1,99
 3038  pdl> p $x
 3039  [
 3040   [ 0  1  2]
 3041   [ 3  4 99]
 3042   [ 6  7  8]
 3043   [ 9 10 11]
 3044  ]
 3045 
 3046 =cut
 3047 
 3048 sub PDL::set{    # Sets a particular single value
 3049     barf 'Usage: set($pdl, $x, $y,.., $value)' if $#_<2;
 3050     my $self  = shift; my $value = pop @_;
 3051     set_c ($self, [@_], $value);
 3052     return $self;
 3053 }
 3054 
 3055 =head2 at
 3056 
 3057 =for ref
 3058 
 3059 Returns a single value inside a piddle as perl scalar.
 3060 
 3061 =for usage
 3062 
 3063  $z = at($piddle, @position); $z=$piddle->at(@position);
 3064 
 3065 C<@position> is a coordinate list, of size equal to the
 3066 number of dimensions in the piddle. Occasionally useful
 3067 in a general context, quite useful too inside PDL internals.
 3068 
 3069 =for example
 3070 
 3071  pdl> $x = sequence 3,4
 3072  pdl> p $x->at(1,2)
 3073  7
 3074 
 3075 If you compile PDL with bad value support (the default), your machine's
 3076 docs will also say this:
 3077 
 3078 =for bad
 3079 
 3080 at converts any bad values into the string 'BAD'.
 3081 
 3082 =cut
 3083 
 3084 sub PDL::at {     # Return value at ($x,$y,$z...)
 3085     barf 'Usage: at($pdl, $x, $y, ...)' if $#_<0;
 3086     my $self = shift;
 3087     at_bad_c ($self, [@_]);
 3088 }
 3089 
 3090 =head2 sclr
 3091 
 3092 =for ref
 3093 
 3094 return a single value from a piddle as a scalar
 3095 
 3096 =for example
 3097 
 3098   $val = $x(10)->sclr;
 3099   $val = sclr inner($x,$y);
 3100 
 3101 The C<sclr> method is useful to turn a piddle into a normal Perl
 3102 scalar. Its main advantage over using C<at> for this purpose is the fact
 3103 that you do not need to worry if the piddle is 0D, 1D or higher dimensional.
 3104 Using C<at> you have to supply the correct number of zeroes, e.g.
 3105 
 3106   $x = sequence(10);
 3107   $y = $x->slice('4');
 3108   print $y->sclr; # no problem
 3109   print $y->at(); # error: needs at least one zero
 3110 
 3111 C<sclr> is generally used when a Perl scalar is required instead
 3112 of a one-element piddle. If the input is a multielement piddle
 3113 the first value is returned as a Perl scalar. You can optionally
 3114 switch on checks to ensure that the input piddle has only one element:
 3115 
 3116   PDL->sclr({Check => 'warn'}); # carp if called with multi-el pdls
 3117   PDL->sclr({Check => 'barf'}); # croak if called with multi-el pdls
 3118 
 3119 are the commands to switch on warnings or raise an error if
 3120 a multielement piddle is passed as input. Note that these options
 3121 can only be set when C<sclr> is called as a class method (see
 3122 example above). Use
 3123 
 3124   PDL->sclr({Check=>0});
 3125 
 3126 to switch these checks off again (default setting);
 3127 When called as a class method the resulting check mode is returned
 3128 (0: no checking, 1: warn, 2: barf).
 3129 
 3130 =cut
 3131 
 3132 my $chkmode = 0; # default mode no checks
 3133 use PDL::Options;
 3134 sub PDL::sclr {
 3135   my $this = shift;
 3136   if (ref $this) { # instance method
 3137     carp "multielement piddle in 'sclr' call"
 3138       if ($chkmode == 1 && $this->nelem > 1);
 3139     croak "multielement piddle in 'sclr' call"
 3140       if ($chkmode == 2 && $this->nelem > 1);
 3141     return sclr_c($this);
 3142   } else {  # class method
 3143     my $check = (iparse({Check=>0},ifhref($_[0])))[1];
 3144     if (lc($check) eq 'warn') {$chkmode = 1}
 3145     elsif (lc($check) eq 'barf') {$chkmode = 2}
 3146     else {$chkmode = $check != 0 ? 1 : 0}
 3147     return $chkmode;
 3148   }
 3149 }
 3150 
 3151 =head2 cat
 3152 
 3153 =for ref
 3154 
 3155 concatenate piddles to N+1 dimensional piddle
 3156 
 3157 Takes a list of N piddles of same shape as argument,
 3158 returns a single piddle of dimension N+1.
 3159 
 3160 =for example
 3161 
 3162  pdl> $x = cat ones(3,3),zeroes(3,3),rvals(3,3); p $x
 3163  [
 3164   [
 3165    [1 1 1]
 3166    [1 1 1]
 3167    [1 1 1]
 3168   ]
 3169   [
 3170    [0 0 0]
 3171    [0 0 0]
 3172    [0 0 0]
 3173   ]
 3174   [
 3175    [1 1 1]
 3176    [1 0 1]
 3177    [1 1 1]
 3178   ]
 3179  ]
 3180 
 3181 If you compile PDL with bad value support (the default), your machine's
 3182 docs will also say this:
 3183 
 3184 =for bad
 3185 
 3186 The output piddle is set bad if any input piddles have their bad flag set.
 3187 
 3188 Similar functions include L<append|PDL::Primitive/append>, which
 3189 appends only two piddles along their first dimension, and
 3190 L<glue|PDL::Primitive/glue>, which can append more than two piddles
 3191 along an arbitrary dimension.
 3192 
 3193 Also consider the generic constructor L<pdl|pdl>, which can handle
 3194 piddles of different sizes (with zero-padding), and will return a
 3195 piddle of type 'double' by default, but may be considerably faster (up
 3196 to 10x) than cat.
 3197 
 3198 =cut
 3199 
 3200 sub PDL::cat {
 3201     my $res;
 3202     my $old_err = $@;
 3203     $@ = '';
 3204     eval {
 3205         $res = $_[0]->initialize;
 3206         $res->set_datatype((sort {$b<=>$a} map{$_->get_datatype} @_)[0] );
 3207 
 3208         my @resdims = $_[0]->dims;
 3209         for my $i(0..$#_){
 3210             my @d = $_[$i]->dims;
 3211             for my $j(0..$#d) {
 3212             $resdims[$j] = $d[$j] if( !defined($resdims[$j]) or $resdims[$j]==1 );
 3213             die "mismatched dims\n" if($d[$j] != 1 and $resdims[$j] != $d[$j]);
 3214             }
 3215         }
 3216         $res->setdims( [@resdims,scalar(@_) ]);
 3217         my ($i,$t); my $s = ":,"x@resdims;
 3218         for (@_) { $t = $res->slice($s."(".$i++.")"); $t .= $_}
 3219 
 3220         # propagate any bad flags
 3221         for (@_) { if ( $_->badflag() ) { $res->badflag(1); last; } }
 3222     };
 3223     if ($@ eq '') {
 3224         # Restore the old error and return
 3225         $@ = $old_err;
 3226         return $res;
 3227     }
 3228 
 3229     # If we've gotten here, then there's been an error, so check things
 3230     # and barf out a meaningful message.
 3231 
 3232     if  ($@ =~ /PDL::Ops::assgn|mismatched/
 3233       or $@ =~ /"badflag"/
 3234       or $@ =~ /"initialize"/) {
 3235         my (@mismatched_dims, @not_a_piddle);
 3236         my $i = 0;
 3237 
 3238         # non-piddles and/or dimension mismatch.  The first argument is
 3239         # ok unless we have the "initialize" error:
 3240         if ($@ =~ /"initialize"/) {
 3241             # Handle the special case that there are *no* args passed:
 3242             barf("Called PDL::cat without any arguments") unless @_;
 3243 
 3244             while ($i < @_ and not eval{ $_[$i]->isa('PDL')}) {
 3245                 push (@not_a_piddle, $i);
 3246                 $i++;
 3247             }
 3248         }
 3249 
 3250         # Get the dimensions of the first actual piddle in the argument
 3251         # list:
 3252         my $first_piddle_argument = $i;
 3253         my @dims = $_[$i]->dims if ref($_[$i]) =~ /PDL/;
 3254 
 3255         # Figure out all the ways that the caller screwed up:
 3256         while ($i < @_) {
 3257             my $arg = $_[$i];
 3258             # Check if not a piddle
 3259             if (not eval{$arg->isa('PDL')}) {
 3260                 push @not_a_piddle, $i;
 3261             }
 3262             # Check if different number of dimensions
 3263             elsif (@dims != $arg->ndims) {
 3264                 push @mismatched_dims, $i;
 3265             }
 3266             # Check if size of dimensions agree
 3267             else {
 3268                 DIMENSION: for (my $j = 0; $j < @dims; $j++) {
 3269                     if ($dims[$j] != $arg->dim($j)) {
 3270                         push @mismatched_dims, $i;
 3271                         last DIMENSION;
 3272                     }
 3273                 }
 3274             }
 3275             $i++;
 3276         }
 3277 
 3278         # Construct a message detailing the results
 3279         my $message = "bad arguments passed to function PDL::cat\n";
 3280         if (@mismatched_dims > 1) {
 3281             # Many dimension mismatches
 3282             $message .= "The dimensions of arguments "
 3283                         . join(', ', @mismatched_dims[0 .. $#mismatched_dims-1])
 3284                         . " and $mismatched_dims[-1] do not match the\n"
 3285                         . "   dimensions of the first piddle argument (argument $first_piddle_argument).\n";
 3286         }
 3287         elsif (@mismatched_dims) {
 3288             # One dimension mismatch
 3289             $message .= "The dimensions of argument $mismatched_dims[0] do not match the\n"
 3290                         . "   dimensions of the first piddle argument (argument $first_piddle_argument).\n";
 3291         }
 3292         if (@not_a_piddle > 1) {
 3293             # many non-piddles
 3294             $message .= "Arguments " . join(', ', @not_a_piddle[0 .. $#not_a_piddle-1])
 3295                         . " and $not_a_piddle[-1] are not piddles.\n";
 3296         }
 3297         elsif (@not_a_piddle) {
 3298             # one non-piddle
 3299             $message .= "Argument $not_a_piddle[0] is not a piddle.\n";
 3300         }
 3301 
 3302         # Handle the edge case that something else happened:
 3303         if (@not_a_piddle == 0 and @mismatched_dims == 0) {
 3304             barf("cat: unknown error from the internals:\n$@");
 3305         }
 3306 
 3307         $message .= "(Argument counting starts from zero.)";
 3308         croak($message);
 3309     }
 3310     else {
 3311         croak("cat: unknown error from the internals:\n$@");
 3312     }
 3313 }
 3314 
 3315 =head2 dog
 3316 
 3317 =for ref
 3318 
 3319 Opposite of 'cat' :). Split N dim piddle to list of N-1 dim piddles
 3320 
 3321 Takes a single N-dimensional piddle and splits it into a list of N-1 dimensional
 3322 piddles. The breakup is done along the last dimension.
 3323 Note the dataflown connection is still preserved by default,
 3324 e.g.:
 3325 
 3326 =for example
 3327 
 3328  pdl> $p = ones 3,3,3
 3329  pdl> ($x,$y,$c) = dog $p
 3330  pdl> $y++; p $p
 3331  [
 3332   [
 3333    [1 1 1]
 3334    [1 1 1]
 3335    [1 1 1]
 3336   ]
 3337   [
 3338    [2 2 2]
 3339    [2 2 2]
 3340    [2 2 2]
 3341   ]
 3342   [
 3343    [1 1 1]
 3344    [1 1 1]
 3345    [1 1 1]
 3346   ]
 3347  ]
 3348 
 3349 =for options
 3350 
 3351  Break => 1   Break dataflow connection (new copy)
 3352 
 3353 If you compile PDL with bad value support (the default), your machine's
 3354 docs will also say this:
 3355 
 3356 =for bad
 3357 
 3358 The output piddles are set bad if the original piddle has its bad flag set.
 3359 
 3360 =cut
 3361 
 3362 sub PDL::dog {
 3363   my $opt = pop @_ if ref($_[-1]) eq 'HASH';
 3364   my $p = shift;
 3365   my @res; my $s = ":,"x($p->getndims-1);
 3366   for my $i (0..$p->getdim($p->getndims-1)-1) {
 3367      $res[$i] = $p->slice($s."(".$i.")");
 3368      $res[$i] = $res[$i]->copy if $$opt{Break};
 3369      $i++;
 3370   }
 3371   return @res;
 3372 }
 3373 
 3374 ###################### Misc internal routines ####################
 3375 
 3376 # Recursively pack an N-D array ref in format [[1,1,2],[2,2,3],[2,2,2]] etc
 3377 # package vars $level and @dims must be initialised first.
 3378 
 3379 sub rpack {
 3380     my ($ptype,$x) = @_;  my ($ret,$type);
 3381 
 3382     $ret = "";
 3383     if (ref($x) eq "ARRAY") {
 3384 
 3385        if (defined($dims[$level])) {
 3386            barf 'Array is not rectangular' unless $dims[$level] == scalar(@$x);
 3387        }else{
 3388           $dims[$level] = scalar(@$x);
 3389        }
 3390 
 3391        $type = ref($$x[0]);
 3392         if ($type) {
 3393         $level++;
 3394         for(@$x) {
 3395             barf 'Array is not rectangular' unless $type eq ref($_); # Equal types
 3396             $ret .= rpack($ptype,$_);
 3397         }
 3398         $level--;
 3399         } else {
 3400         # These are leaf nodes
 3401         $ret = pack $ptype, map {defined($_) ? $_ : $PDL::undefval} @$x;
 3402       }
 3403     } elsif (ref($x) eq "PDL") {
 3404     barf 'Cannot make a new piddle from two or more piddles, try "cat"';
 3405     } else {
 3406         barf "Don't know how to make a PDL object from passed argument";
 3407     }
 3408     return $ret;
 3409 }
 3410 
 3411 sub rcopyitem {        # Return a deep copy of an item - recursively
 3412     my $x = shift;
 3413     my ($y, $key, $value);
 3414     if (ref(\$x) eq "SCALAR") {
 3415        return $x;
 3416     }elsif (ref($x) eq "SCALAR") {
 3417        $y = $$x; return \$y;
 3418     }elsif (ref($x) eq "ARRAY") {
 3419        $y = [];
 3420        for (@$x) {
 3421            push @$y, rcopyitem($_);
 3422        }
 3423        return $y;
 3424     }elsif (ref($x) eq "HASH") {
 3425        $y={};
 3426        while (($key,$value) = each %$x) {
 3427           $$y{$key} = rcopyitem($value);
 3428        }
 3429        return $y;
 3430     }elsif (blessed($x)) {
 3431        return $x->copy;
 3432     }else{
 3433        barf ('Deep copy of object failed - unknown component with type '.ref($x));
 3434     }
 3435 0;}
 3436 
 3437 # N-D array stringifier
 3438 
 3439 sub strND {
 3440     my($self,$format,$level)=@_;
 3441 #    $self->make_physical();
 3442     my @dims = $self->dims;
 3443     # print "STRND, $#dims\n";
 3444 
 3445     if ($#dims==1) { # Return 2D string
 3446        return str2D($self,$format,$level);
 3447     }
 3448     else { # Return list of (N-1)D strings
 3449        my $secbas = join '',map {":,"} @dims[0..$#dims-1];
 3450        my $ret="\n"." "x$level ."["; my $j;
 3451        for ($j=0; $j<$dims[$#dims]; $j++) {
 3452            my $sec = $secbas . "($j)";
 3453 #      print "SLICE: $sec\n";
 3454 
 3455            $ret .= strND($self->slice($sec),$format, $level+1);
 3456        chop $ret; $ret .= $sep2;
 3457        }
 3458        chop $ret if $PDL::use_commas;
 3459        $ret .= "\n" ." "x$level ."]\n";
 3460        return $ret;
 3461     }
 3462 }
 3463 
 3464 
 3465 # String 1D array in nice format
 3466 
 3467 sub str1D {
 3468     my($self,$format)=@_;
 3469     barf "Not 1D" if $self->getndims()!=1;
 3470     my $x = listref_c($self);
 3471     my ($ret,$dformat,$t);
 3472     $ret = "[";
 3473     my $dtype = $self->get_datatype();
 3474     $dformat = $PDL::floatformat  if $dtype == $PDL_F;
 3475     $dformat = $PDL::doubleformat if $dtype == $PDL_D;
 3476     $dformat = $PDL::indxformat if $dtype == $PDL_IND;
 3477 
 3478     my $badflag = $self->badflag();
 3479     for $t (@$x) {
 3480     if ( $badflag and $t eq "BAD" ) {
 3481         # do nothing
 3482         } elsif ($format) {
 3483       $t = sprintf $format,$t;
 3484     } else{ # Default
 3485         if ($dformat && length($t)>7) { # Try smaller
 3486         $t = sprintf $dformat,$t;
 3487         }
 3488     }
 3489        $ret .= $t.$sep;
 3490     }
 3491 
 3492     chop $ret; $ret.="]";
 3493     return $ret;
 3494 }
 3495 
 3496 # String 2D array in nice uniform format
 3497 
 3498 sub str2D{
 3499     my($self,$format,$level)=@_;
 3500 #    print "STR2D:\n"; $self->printdims();
 3501     my @dims = $self->dims();
 3502     barf "Not 2D" if scalar(@dims)!=2;
 3503     my $x = listref_c($self);
 3504     my ($i, $f, $t, $len, $ret);
 3505 
 3506     my $dtype = $self->get_datatype();
 3507     my $badflag = $self->badflag();
 3508 
 3509     my $findmax = 1;
 3510     if (!defined $format || $format eq "") {
 3511     # Format not given? - find max length of default
 3512     $len=0;
 3513 
 3514     if ( $badflag ) {
 3515         for (@$x) {
 3516         if ( $_ eq "BAD" ) { $i = 3; }
 3517         else               { $i = length($_); }
 3518         $len = $i>$len ? $i : $len;
 3519         }
 3520     } else {
 3521         for (@$x) {$i = length($_); $len = $i>$len ? $i : $len };
 3522     }
 3523 
 3524     $format = "%".$len."s";
 3525 
 3526     if ($len>7) { # Too long? - perhaps try smaller format
 3527         if ($dtype == $PDL_F) {
 3528         $format = $PDL::floatformat;
 3529         } elsif ($dtype == $PDL_D) {
 3530         $format = $PDL::doubleformat;
 3531         } elsif ($dtype == $PDL_IND) {
 3532         $format = $PDL::indxformat;
 3533         } else {
 3534         # Stick with default
 3535         $findmax = 0;
 3536         }
 3537     }
 3538     else {
 3539         # Default ok
 3540         $findmax = 0;
 3541     }
 3542     }
 3543 
 3544     if($findmax) {
 3545     # Find max length of strings in final format
 3546     $len=0;
 3547 
 3548     if ( $badflag ) {
 3549         for (@$x) {
 3550         if ( $_ eq "BAD" ) { $i = 3; }
 3551         else               { $i = length(sprintf $format,$_); }
 3552         $len = $i>$len ? $i : $len;
 3553         }
 3554     } else {
 3555         for (@$x) {
 3556         $i = length(sprintf $format,$_); $len = $i>$len ? $i : $len;
 3557         }
 3558     }
 3559     } # if: $findmax
 3560 
 3561     $ret = "\n" . " "x$level . "[\n";
 3562     {
 3563     my $level = $level+1;
 3564     $ret .= " "x$level ."[";
 3565     for ($i=0; $i<=$#$x; $i++) {
 3566 
 3567         if ( $badflag and $$x[$i] eq "BAD" ) {
 3568         $f = "BAD";
 3569         } else {
 3570         $f = sprintf $format,$$x[$i];
 3571         }
 3572 
 3573         $t = $len-length($f); $f = " "x$t .$f if $t>0;
 3574         $ret .= $f;
 3575         if (($i+1)%$dims[0]) {
 3576         $ret.=$sep;
 3577         }
 3578         else{ # End of output line
 3579         $ret.="]";
 3580         if ($i==$#$x) { # very last number
 3581             $ret.="\n";
 3582         }
 3583         else{
 3584             $ret.= $sep2."\n" . " "x$level ."[";
 3585         }
 3586         }
 3587     }
 3588     }
 3589     $ret .= " "x$level."]\n";
 3590     return $ret;
 3591 }
 3592 
 3593 #
 3594 # Sleazy hcpy saves me time typing
 3595 #
 3596 sub PDL::hcpy {
 3597   $_[0]->hdrcpy($_[1]);
 3598   $_[0];
 3599 }
 3600 
 3601 ########## Docs for functions in Core.xs ##################
 3602 # Pod docs for functions that are imported from Core.xs and are
 3603 #  not documented elsewhere. Currently this is not a complete
 3604 #  list. There are others.
 3605 
 3606 =head2 gethdr
 3607 
 3608 =for ref
 3609 
 3610 Retrieve header information from a piddle
 3611 
 3612 =for example
 3613 
 3614  $pdl=rfits('file.fits');
 3615  $h=$pdl->gethdr;
 3616  print "Number of pixels in the X-direction=$$h{NAXIS1}\n";
 3617 
 3618 The C<gethdr> function retrieves whatever header information is contained
 3619 within a piddle. The header can be set with L<sethdr|/sethdr> and is always a
 3620 hash reference or undef.
 3621 
 3622 C<gethdr> returns undef if the piddle has not yet had a header
 3623 defined; compare with C<hdr> and C<fhdr>, which are guaranteed to return a
 3624 defined value.
 3625 
 3626 Note that gethdr() works by B<reference>: you can modify the header
 3627 in-place once it has been retrieved:
 3628 
 3629   $x  = rfits($filename);
 3630   $xh = $x->gethdr();
 3631   $xh->{FILENAME} = $filename;
 3632 
 3633 It is also important to realise that in most cases the header is not
 3634 automatically copied when you copy the piddle.  See L<hdrcpy|/hdrcpy>
 3635 to enable automatic header copying.
 3636 
 3637 Here's another example: a wrapper around rcols that allows your piddle
 3638 to remember the file it was read from and the columns could be easily
 3639 written (here assuming that no regexp is needed, extensions are left
 3640 as an exercise for the reader)
 3641 
 3642  sub ext_rcols {
 3643     my ($file, @columns)=@_;
 3644     my $header={};
 3645     $$header{File}=$file;
 3646     $$header{Columns}=\@columns;
 3647 
 3648     @piddles=rcols $file, @columns;
 3649     foreach (@piddles) { $_->sethdr($header); }
 3650     return @piddles;
 3651  }
 3652 
 3653 =head2 hdr
 3654 
 3655 =for ref
 3656 
 3657 Retrieve or set header information from a piddle
 3658 
 3659 =for example
 3660 
 3661  $pdl->hdr->{CDELT1} = 1;
 3662 
 3663 The C<hdr> function allows convenient access to the header of a
 3664 piddle.  Unlike C<gethdr> it is guaranteed to return a defined value,
 3665 so you can use it in a hash dereference as in the example.  If the
 3666 header does not yet exist, it gets autogenerated as an empty hash.
 3667 
 3668 Note that this is usually -- but not always -- What You Want.  If you
 3669 want to use a tied L<Astro::FITS::Header|Astro::FITS::Header> hash,
 3670 for example, you should either construct it yourself and use C<sethdr>
 3671 to put it into the piddle, or use L<fhdr|fhdr> instead.  (Note that
 3672 you should be able to write out the FITS file successfully regardless
 3673 of whether your PDL has a tied FITS header object or a vanilla hash).
 3674 
 3675 =head2 fhdr
 3676 
 3677 =for ref
 3678 
 3679 Retrieve or set FITS header information from a piddle
 3680 
 3681 =for example
 3682 
 3683  $pdl->fhdr->{CDELT1} = 1;
 3684 
 3685 The C<fhdr> function allows convenient access to the header of a
 3686 piddle.  Unlike C<gethdr> it is guaranteed to return a defined value,
 3687 so you can use it in a hash dereference as in the example.  If the
 3688 header does not yet exist, it gets autogenerated as a tied
 3689 L<Astro::FITS::Header|Astro::FITS::Header> hash.
 3690 
 3691 Astro::FITS::Header tied hashes are better at matching the behavior of
 3692 FITS headers than are regular hashes.  In particular, the hash keys
 3693 are CAsE INsEnSItiVE, unlike normal hash keys.  See
 3694 L<Astro::FITS::Header> for details.
 3695 
 3696 If you do not have Astro::FITS::Header installed, you get back a
 3697 normal hash instead of a tied object.
 3698 
 3699 =head2 sethdr
 3700 
 3701 =for ref
 3702 
 3703 Set header information of a piddle
 3704 
 3705 =for example
 3706 
 3707  $pdl = zeroes(100,100);
 3708  $h = {NAXIS=>2, NAXIS1=>100, NAXIS=>100, COMMENT=>"Sample FITS-style header"};
 3709  # add a FILENAME field to the header
 3710  $$h{FILENAME} = 'file.fits';
 3711  $pdl->sethdr( $h );
 3712 
 3713 The C<sethdr> function sets the header information for a piddle.
 3714 You must feed in a hash ref or undef, and the header field of the PDL is
 3715 set to be a new ref to the same hash (or undefined).
 3716 
 3717 The hash ref requirement is a speed bump put in place since the normal
 3718 use of headers is to store fits header information and the like.  Of course,
 3719 if you want you can hang whatever ugly old data structure you want off
 3720 of the header, but that makes life more complex.
 3721 
 3722 Remember that the hash is not copied -- the header is made into a ref
 3723 that points to the same underlying data.  To get a real copy without
 3724 making any assumptions about the underlying data structure, you
 3725 can use one of the following:
 3726 
 3727   use PDL::IO::Dumper;
 3728   $pdl->sethdr( deep_copy($h) );
 3729 
 3730 (which is slow but general), or
 3731 
 3732   $pdl->sethdr( PDL::_hdr_copy($h) )
 3733 
 3734 (which uses the built-in sleazy deep copier), or (if you know that all
 3735 the elements happen to be scalars):
 3736 
 3737   { my %a = %$h;
 3738     $pdl->sethdr(\%a);
 3739   }
 3740 
 3741 which is considerably faster but just copies the top level.
 3742 
 3743 The C<sethdr> function must be given a hash reference or undef.  For
 3744 further information on the header, see L<gethdr|/gethdr>, L<hdr|/hdr>,
 3745 L<fhdr|/fhdr> and L<hdrcpy|/hdrcpy>.
 3746 
 3747 =head2 hdrcpy
 3748 
 3749 =for ref
 3750 
 3751 switch on/off/examine automatic header copying
 3752 
 3753 =for example
 3754 
 3755  print "hdrs will be copied" if $x->hdrcpy;
 3756  $x->hdrcpy(1);       # switch on automatic header copying
 3757  $y = $x->sumover;    # and $y will inherit $x's hdr
 3758  $x->hdrcpy(0);       # and now make $x non-infectious again
 3759 
 3760 C<hdrcpy> without an argument just returns the current setting of the
 3761 flag.  See also "hcpy" which returns its PDL argument (and so is useful
 3762 in method-call pipelines).
 3763 
 3764 Normally, the optional header of a piddle is not copied automatically
 3765 in pdl operations. Switching on the hdrcpy flag using the C<hdrcpy>
 3766 method will enable automatic hdr copying. Note that an actual deep
 3767 copy gets made, which is rather processor-inefficient -- so avoid
 3768 using header copying in tight loops!
 3769 
 3770 Most PDLs have the C<hdrcpy> flag cleared by default; however, some
 3771 routines (notably L<rfits|PDL::IO::FITS/rfits()>) set it by default
 3772 where that makes more sense.
 3773 
 3774 The C<hdrcpy> flag is viral: if you set it for a PDL, then derived
 3775 PDLs will get copies of the header and will also have their C<hdrcpy>
 3776 flags set.  For example:
 3777 
 3778   $x = xvals(50,50);
 3779   $x->hdrcpy(1);
 3780   $x->hdr->{FOO} = "bar";
 3781   $y = $x++;
 3782   $c = $y++;
 3783   print $y->hdr->{FOO}, " - ", $c->hdr->{FOO}, "\n";
 3784   $y->hdr->{FOO} = "baz";
 3785   print $x->hdr->{FOO}, " - ", $y->hdr->{FOO}, " - ", $c->hdr->{FOO}, "\n";
 3786 
 3787 will print:
 3788 
 3789   bar - bar
 3790   bar - baz - bar
 3791 
 3792 Performing an operation in which more than one PDL has its hdrcpy flag
 3793 causes the resulting PDL to take the header of the first PDL:
 3794 
 3795   ($x,$y) = sequence(5,2)->dog;
 3796   $x->hdrcpy(1); $y->hdrcpy(1);
 3797   $x->hdr->{foo} = 'a';
 3798   $y->hdr->{foo} = 'b';
 3799   print (($x+$y)->hdr->{foo} , ($y+$x)->hdr->{foo});
 3800 
 3801 will print:
 3802 
 3803   a b
 3804 
 3805 =head2 hcpy
 3806 
 3807 =for ref
 3808 
 3809 Switch on/off automatic header copying, with PDL pass-through
 3810 
 3811 =for example
 3812 
 3813   $x = rfits('foo.fits')->hcpy(0);
 3814   $x = rfits('foo.fits')->hcpy(1);
 3815 
 3816 C<hcpy> sets or clears the hdrcpy flag of a PDL, and returns the PDL
 3817 itself.  That makes it convenient for inline use in expressions.
 3818 
 3819 =head2 set_autopthread_targ
 3820 
 3821 =for ref
 3822 
 3823 Set the target number of processor threads (pthreads) for multi-threaded processing.
 3824 
 3825 =for usage
 3826 
 3827  set_autopthread_targ($num_pthreads);
 3828 
 3829 C<$num_pthreads> is the target number of pthreads the auto-pthread process will try to achieve.
 3830 
 3831 See L<PDL::ParallelCPU> for an overview of the auto-pthread process.
 3832 
 3833 =for example
 3834 
 3835   # Example turning on auto-pthreading for a target of 2 pthreads and for functions involving
 3836   #   PDLs with greater than 1M elements
 3837   set_autopthread_targ(2);
 3838   set_autopthread_size(1);
 3839 
 3840   # Execute a pdl function, processing will split into two pthreads as long as
 3841   #  one of the pdl-threaded dimensions is divisible by 2.
 3842   $x = minimum($y);
 3843 
 3844   # Get the actual number of pthreads that were run.
 3845   $actual_pthread = get_autopthread_actual();
 3846 
 3847 =cut
 3848 
 3849 *set_autopthread_targ       = \&PDL::set_autopthread_targ;
 3850 
 3851 =head2 get_autopthread_targ
 3852 
 3853 =for ref
 3854 
 3855 Get the current target number of processor threads (pthreads) for multi-threaded processing.
 3856 
 3857 =for usage
 3858 
 3859  $num_pthreads = get_autopthread_targ();
 3860 
 3861 C<$num_pthreads> is the target number of pthreads the auto-pthread process will try to achieve.
 3862 
 3863 See L<PDL::ParallelCPU> for an overview of the auto-pthread process.
 3864 
 3865 =cut
 3866 
 3867 *get_autopthread_targ       = \&PDL::get_autopthread_targ;
 3868 
 3869 =head2 get_autopthread_actual
 3870 
 3871 =for ref
 3872 
 3873 Get the actual number of pthreads executed for the last pdl processing function.
 3874 
 3875 =for usage
 3876 
 3877  $autopthread_actual = get_autopthread_actual();
 3878 
 3879 C<$autopthread_actual> is the actual number of pthreads executed for the last pdl processing function.
 3880 
 3881 See L<PDL::ParallelCPU> for an overview of the auto-pthread process.
 3882 
 3883 =cut
 3884 
 3885 *get_autopthread_actual      = \&PDL::get_autopthread_actual;
 3886 
 3887 =head2 set_autopthread_size
 3888 
 3889 =for ref
 3890 
 3891 Set the minimum size (in M-elements or 2^20 elements) of the largest PDL involved in a function where auto-pthreading will
 3892 be performed. For small PDLs, it probably isn't worth starting multiple pthreads, so this function
 3893 is used to define a minimum threshold where auto-pthreading won't be attempted.
 3894 
 3895 =for usage
 3896 
 3897  set_autopthread_size($size);
 3898 
 3899 C<$size> is the mimumum size, in M-elements or 2^20 elements (approx 1e6 elements) for the largest PDL involved in a function.
 3900 
 3901 See L<PDL::ParallelCPU> for an overview of the auto-pthread process.
 3902 
 3903 =for example
 3904 
 3905   # Example turning on auto-pthreading for a target of 2 pthreads and for functions involving
 3906   #   PDLs with greater than 1M elements
 3907   set_autopthread_targ(2);
 3908   set_autopthread_size(1);
 3909 
 3910   # Execute a pdl function, processing will split into two pthreads as long as
 3911   #  one of the pdl-threaded dimensions is divisible by 2.
 3912   $x = minimum($y);
 3913 
 3914   # Get the actual number of pthreads that were run.
 3915   $actual_pthread = get_autopthread_actual();
 3916 
 3917 =cut
 3918 
 3919 *set_autopthread_size       = \&PDL::set_autopthread_size;
 3920 
 3921 =head2 get_autopthread_size
 3922 
 3923 =for ref
 3924 
 3925 Get the current autopthread_size setting.
 3926 
 3927 =for usage
 3928 
 3929  $autopthread_size = get_autopthread_size();
 3930 
 3931 C<$autopthread_size> is the mimumum size limit for auto_pthreading to occur, in M-elements or 2^20 elements (approx 1e6 elements) for the largest PDL involved in a function
 3932 
 3933 See L<PDL::ParallelCPU> for an overview of the auto-pthread process.
 3934 
 3935 =cut
 3936 
 3937 *get_autopthread_size       = \&PDL::get_autopthread_size;
 3938 
 3939 =head1 AUTHOR
 3940 
 3941 Copyright (C) Karl Glazebrook (kgb@aaoepp.aao.gov.au),
 3942 Tuomas J. Lukka, (lukka@husc.harvard.edu) and Christian
 3943 Soeller (c.soeller@auckland.ac.nz) 1997.
 3944 Modified, Craig DeForest (deforest@boulder.swri.edu) 2002.
 3945 All rights reserved. There is no warranty. You are allowed
 3946 to redistribute this software / documentation under certain
 3947 conditions. For details, see the file COPYING in the PDL
 3948 distribution. If this file is separated from the PDL distribution,
 3949 the copyright notice should be included in the file.
 3950 
 3951 =cut
 3952 
 3953 #
 3954 # Easier to implement in perl than in XS...
 3955 #   -- CED
 3956 #
 3957 
 3958 sub PDL::fhdr {
 3959     my $pdl = shift;
 3960 
 3961     return $pdl->hdr
 3962     if( (defined $pdl->gethdr) ||
 3963     !defined $Astro::FITS::Header::VERSION
 3964         );
 3965 
 3966     # Avoid bug in 1.15 and earlier Astro::FITS::Header
 3967     my @hdr = ("SIMPLE  =                    T");
 3968     my $hdr = new Astro::FITS::Header(Cards=>\@hdr);
 3969     tie my %hdr, "Astro::FITS::Header", $hdr;
 3970     $pdl->sethdr(\%hdr);
 3971     return \%hdr;
 3972 }
 3973 
 3974 use Fcntl;
 3975 
 3976 BEGIN {
 3977    eval 'use File::Map 0.47 qw(:all)';
 3978    if ($@) {
 3979       carp "No File::Map found, using legacy mmap (if available)\n" if $PDL::verbose;
 3980       sub sys_map;
 3981       sub PROT_READ();
 3982       sub PROT_WRITE();
 3983       sub MAP_SHARED();
 3984       sub MAP_PRIVATE();
 3985    }
 3986 }
 3987 
 3988 # Implement File::Map::sys_map bug fix.  Also, might be possible
 3989 # to implement without so many external (non-Core perl) modules.
 3990 #
 3991 # sub pdl_do_sys_map {
 3992 #         my (undef, $length, $protection, $flags, $fh, $offset) = @_;
 3993 #         my $utf8 = File::Map::_check_layers($fh);
 3994 #         my $fd = ($flags & MAP_ANONYMOUS) ? (-1) : fileno($fh);
 3995 #         $offset ||= 0;
 3996 #         File::Map::_mmap_impl($_[0], $length, $protection, $flags, $fd, $offset, $utf8);
 3997 #         return;
 3998 # }
 3999 
 4000 sub PDL::set_data_by_file_map {
 4001    my ($pdl,$name,$len,$shared,$writable,$creat,$mode,$trunc) = @_;
 4002    my $pdl_dataref = $pdl->get_dataref();
 4003 
 4004    # Assume we have no data to free for now
 4005    # pdl_freedata($pdl);
 4006 
 4007    sysopen(my $fh, $name, ($writable && $shared ? O_RDWR : O_RDONLY) | ($creat ? O_CREAT : 0), $mode)
 4008       or die "Error opening file '$name'\n";
 4009 
 4010    binmode $fh;
 4011 
 4012    if ($trunc) {
 4013       truncate($fh,0) or die "set_data_by_mmap: truncate('$name',0) failed, $!";
 4014       truncate($fh,$len) or die "set_data_by_mmap: truncate('$name',$len) failed, $!";
 4015    }
 4016 
 4017    if ($len) {
 4018 
 4019       #eval {
 4020       # pdl_do_sys_map(  # will croak if the mapping fails
 4021       if ($PDL::debug) {
 4022          printf STDERR
 4023          "set_data_by_file_map: calling sys_map(%s,%d,%d,%d,%s,%d)\n",
 4024          $pdl_dataref,
 4025          $len,
 4026          PROT_READ | ($writable ?  PROT_WRITE : 0),
 4027          ($shared ? MAP_SHARED : MAP_PRIVATE),
 4028          $fh,
 4029          0;
 4030       }
 4031 
 4032       sys_map(  # will croak if the mapping fails
 4033          ${$pdl_dataref},
 4034          $len,
 4035          PROT_READ | ($writable ?  PROT_WRITE : 0),
 4036          ($shared ? MAP_SHARED : MAP_PRIVATE),
 4037          $fh,
 4038          0
 4039       );
 4040       #};
 4041 
 4042          #if ($@) {
 4043          #die("Error mmapping!, '$@'\n");
 4044          #}
 4045 
 4046       $pdl->upd_data;
 4047 
 4048       if ($PDL::debug) {
 4049          printf STDERR "set_data_by_file_map: length \${\$pdl_dataref} is %d.\n", length ${$pdl_dataref};
 4050       }
 4051       $pdl->set_state_and_add_deletedata_magic( length ${$pdl_dataref} );
 4052 
 4053    } else {
 4054 
 4055       #  Special case: zero-length file
 4056       $_[0] = undef;
 4057    }
 4058 
 4059    # PDLDEBUG_f(printf("PDL::MMap: mapped to %p\n",$pdl->data));
 4060    close $fh ;
 4061 }
 4062 
 4063 1;