File Coverage

/root/.cpan/build/PDL-CCS-1.23.11-0/blib/lib/PDL/CCS/Nd.pm
Criterion Covered Total %
statement 51 915 5.5
branch 2 432 0.4
condition 0 120 0.0
subroutine 20 164 12.2
pod 61 111 54.9
total 134 1742 7.6


line stmt bran cond sub pod time code
1             ## File: PDL::CCS::Nd.pm
2             ## Author: Bryan Jurish
3             ## Description: N-dimensional CCS-encoded pseudo-PDL
4              
5             package PDL::CCS::Nd;
6 4     4   1805 use PDL::Lite qw();
  4         2703  
  4         97  
7 4     4   28 use PDL::VectorValued;
  4         9  
  4         26  
8 4     4   776 use PDL::CCS::Config qw(ccs_indx);
  4         8  
  4         217  
9 4     4   25 use PDL::CCS::Functions qw(ccs_decode ccs_pointerlen ccs_qsort);
  4         9  
  4         28  
10 4     4   403 use PDL::CCS::Utils qw(ccs_encode_pointers ccs_decode_pointer);
  4         10  
  4         27  
11 4     4   348 use PDL::CCS::Ufunc;
  4         8  
  4         27  
12 4     4   595 use PDL::CCS::Ops;
  4         11  
  4         22  
13 4     4   361 use PDL::CCS::MatrixOps;
  4         8  
  4         22  
14 4     4   477 use Carp;
  4         9  
  4         229  
15 4     4   25 use strict;
  4         26  
  4         219  
16              
17             BEGIN {
18 4     4   18 *isa = \&UNIVERSAL::isa;
19 4         7139 *can = \&UNIVERSAL::can;
20             }
21              
22             our $VERSION = '1.23.11'; ##-- update with perl-reversion from Perl::Version module
23             our @ISA = qw();
24             our %EXPORT_TAGS =
25             (
26             ##-- respect PDL conventions (hopefully)
27             Func => [
28             ##-- Encoding/Decoding
29             qw(toccs todense),
30             ],
31             vars => [
32             qw($PDIMS $VDIMS $WHICH $VALS $PTRS $FLAGS $USER),
33             qw($BINOP_BLOCKSIZE_MIN $BINOP_BLOCKSIZE_MAX),
34             ],
35             flags => [
36             qw($CCSND_BAD_IS_MISSING $CCSND_NAN_IS_MISSING $CCSND_INPLACE $CCSND_FLAGS_DEFAULT),
37             ],
38             );
39             $EXPORT_TAGS{all} = [map {@$_} values(%EXPORT_TAGS)];
40             our @EXPORT = @{$EXPORT_TAGS{Func}};
41             our @EXPORT_OK = @{$EXPORT_TAGS{all}};
42              
43             ##--------------------------------------------------------------
44             ## Global variables for block-wise computation of binary operations
45              
46             ##-- some (hopefully sensible) defaults
47             #our $BINOP_BLOCKSIZE_MIN = 64;
48             #our $BINOP_BLOCKSIZE_MAX = undef; ##-- undef or zero: no maximum
49              
50             ##-- debug/devel defaults
51             our $BINOP_BLOCKSIZE_MIN = 1;
52             our $BINOP_BLOCKSIZE_MAX = 0;
53              
54             ##======================================================================
55             ## Globals
56              
57             our $PDIMS = 0;
58             our $VDIMS = 1;
59             our $WHICH = 2;
60             our $VALS = 3;
61             our $PTRS = 4;
62             our $FLAGS = 5;
63             our $USER = 6;
64              
65             ##-- flags
66             our $CCSND_BAD_IS_MISSING = 1;
67             our $CCSND_NAN_IS_MISSING = 2;
68             our $CCSND_INPLACE = 4;
69             our $CCSND_FLAGS_DEFAULT = 0; ##-- default flags
70              
71             ##-- pdl constants
72             our $P_BYTE = PDL::byte();
73             our $P_LONG = PDL::long();
74             our $P_INDX = ccs_indx();
75              
76 0 0   0   0 sub _min2 ($$) { $_[0]<$_[1] ? $_[0] : $_[1]; }
77 0 0   0   0 sub _max2 ($$) { $_[0]>$_[1] ? $_[0] : $_[1]; }
78              
79             ##======================================================================
80             ## Constructors etc.
81              
82             ## $obj = $class_or_obj->newFromDense($denseND);
83             ## $obj = $class_or_obj->newFromDense($denseND,$missing);
84             ## $obj = $class_or_obj->newFromDense($denseND,$missing,$flags);
85             ## + object structure: ARRAY
86             ## $PDIMS => $pdims, ##-- pdl(indx,$NPdims) : physical dimension sizes : $pdim_i => $dimSize_i
87             ## $VDIMS => $vdims, ##-- pdl(indx,$NVdims) : virtual dimension sizes
88             ## ## + $vdim_i => / -$vdimSize_i if $vdim_i is dummy
89             ## ## \ $pdim_i otherwise
90             ## ## + s.t. $whichND_logical_physical = $whichND->dice_axis(0,$vdims->where($vdims>=0));
91             ## $WHICH => $whichND, ##-- pdl(indx,$NPdims,$Nnz) ~ $dense_orig->whichND
92             ## ## + guaranteed to be sorted as for qsortvec() specs
93             ## ## + NOT changed by dimension-shuffling transformations
94             ## $VALS => $vals, ##-- pdl( ? ,$Nnz+1) ~ $dense->where($dense)->append($missing)
95             ## $PTRS => \@PTRS, ##-- array of ccsutils-pointers by physical dimension number
96             ## $FLAGS => $flags, ##-- integer holding some flags
97             ##
98             ## + each element of @PTRS is itself an array:
99             ## $PTRS[$i] => [ $PTR, $NZI ]
100             ##
101             sub newFromDense :lvalue {
102 0     0 1 0 my $that = shift;
103 0   0     0 return my $tmp=(bless [], ref($that)||$that)->fromDense(@_);
104             }
105              
106             ## $obj = $obj->fromDense($denseND,$missing,$flags)
107             sub fromDense :lvalue {
108 0     0 1 0 my ($obj,$p,$missing,$flags) = @_;
109 0         0 $p = PDL->topdl($p);
110 0 0       0 $p = $p->slice("*1") if (!$p->dims);
111 0 0       0 $missing = (defined($missing)
    0          
112             ? PDL->pdl($p->type,$missing)
113             : ($p->badflag
114             ? PDL->pdl($p->type,0)->setvaltobad(0)
115             : PDL->pdl($p->type,0)));
116 0 0       0 $flags = $CCSND_FLAGS_DEFAULT if (!defined($flags));
117 0 0       0 my $pwhichND = ($missing->isbad ? $p->isgood() : ($p != $missing))->whichND->vv_qsortvec;
118 0         0 my $pnz = $p->indexND($pwhichND)->append($missing);
119 0         0 $pnz->sever; ##-- always sever nzvals ?
120 0         0 my $pdims = PDL->pdl($P_INDX,[$p->dims]);
121 0         0 $obj->[$PDIMS] = $pdims;
122 0 0       0 $obj->[$VDIMS] = $pdims->isempty ? $pdims->pdl : $pdims->sequence;
123 0         0 $obj->[$WHICH] = $pwhichND;
124 0         0 $obj->[$VALS] = $pnz;
125 0         0 $obj->[$PTRS] = []; ##-- do we really need this ... yes
126 0         0 $obj->[$FLAGS] = $flags;
127 0         0 return $obj;
128             }
129              
130             ## $obj = $class_or_obj->newFromWhich($whichND,$nzvals,%options);
131             ## $obj = $class_or_obj->newFromWhich($whichND,$nzvals);
132             ## + %options: see $obj->fromWhich()
133             sub newFromWhich :lvalue {
134 0     0 1 0 my $that = shift;
135 0   0     0 return my $tmp=bless([],ref($that)||$that)->fromWhich(@_);
136             }
137              
138             ## $obj = $obj->fromWhich($whichND,$nzvals,%options);
139             ## $obj = $obj->fromWhich($whichND,$nzvals);
140             ## + %options:
141             ## sorted => $bool, ##-- if true, $whichND is assumed to be pre-sorted
142             ## steal => $bool, ##-- if true, $whichND and $nzvals are used literally (formerly implied 'sorted')
143             ## ## + in this case, $nzvals should really be: $nzvals->append($missing)
144             ## pdims => $pdims, ##-- physical dimension list; default guessed from $whichND (alias: 'dims')
145             ## missing => $missing, ##-- default: BAD if $nzvals->badflag, 0 otherwise
146             ## vdims => $vdims, ##-- virtual dims (default: sequence($nPhysDims)); alias: 'xdims'
147             ## flags => $flags, ##-- flags
148             sub fromWhich :lvalue {
149 0     0 1 0 my ($obj,$wnd,$nzvals,%opts) = @_;
150             my $missing = (defined($opts{missing})
151             ? PDL->pdl($nzvals->type,$opts{missing})
152 0 0       0 : ($nzvals->badflag
    0          
153             ? PDL->pdl($nzvals->type,0)->setvaltobad(0)
154             : PDL->pdl($nzvals->type,0)));
155              
156             ##-- get dims
157 0   0     0 my $pdims = $opts{pdims} // $opts{dims} // PDL->pdl($P_INDX, [($wnd->xchg(0,1)->maximum+1)->list]);
      0        
158 0 0       0 $pdims = PDL->pdl($P_INDX, $pdims) if (!UNIVERSAL::isa($pdims,'PDL'));
159              
160 0   0     0 my $vdims = $opts{vdims} // $opts{xdims} // $pdims->sequence;
      0        
161 0 0       0 $vdims = PDL->pdl($P_INDX, $vdims) if (!UNIVERSAL::isa($vdims,'PDL'));
162              
163             ##-- maybe sort & copy
164 0 0       0 if (!$opts{steal}) {
    0          
165             ##-- not stolen: copy or sever
166 0 0       0 if (!$opts{sorted}) {
167 0 0       0 my $wi = $wnd->isempty ? PDL->null->ccs_indx() : $wnd->vv_qsortveci;
168 0         0 $wnd = $wnd->dice_axis(1,$wi);
169 0         0 $nzvals = $nzvals->index($wi);
170             }
171 0         0 $wnd->sever; ##-- sever (~ copy)
172 0         0 $nzvals = $nzvals->append($missing); ##-- copy (b/c append)
173             }
174             elsif (!$opts{sorted}) {
175             ##-- "stolen" but un-sorted: we have "missing" value in $vals
176 0         0 my $wi = PDL->zeroes(ccs_indx, $wnd->dim(1)+1);
177 0         0 $wnd->vv_qsortveci($wi->slice("0:-2"));
178 0         0 $wi->set($wnd->dim(1) => $nzvals->nelem-1);
179 0         0 $wnd = $wnd->dice_axis(1,$wi->slice("0:-2"));
180 0         0 $nzvals = $nzvals->index($wi);
181             }
182              
183             ##-- setup and return
184 0         0 $obj->[$PDIMS] = $pdims;
185 0         0 $obj->[$VDIMS] = $vdims;
186 0         0 $obj->[$WHICH] = $wnd;
187 0         0 $obj->[$VALS] = $nzvals;
188 0         0 $obj->[$PTRS] = [];
189 0 0       0 $obj->[$FLAGS] = defined($opts{flags}) ? $opts{flags} : $CCSND_FLAGS_DEFAULT;
190 0         0 return $obj;
191             }
192              
193              
194             ## DESTROY : avoid PDL inheritance
195       0     sub DESTROY { ; }
196              
197             ## $ccs = $ccs->insertWhich($whichND,$whichVals)
198             ## + set or insert $whichND=>$whichVals
199             ## + implicitly calls make_physically_indexed
200             sub insertWhich :lvalue {
201 0     0 1 0 my ($ccs,$which,$vals) = @_;
202 0         0 $ccs->make_physically_indexed();
203              
204             ##-- sanity check
205 0 0       0 if ($which->dim(0) != $ccs->[$WHICH]->dim(0)) {
206 0         0 PDL::Lite::barf(ref($ccs)."::insertWhich(): wrong number of index dimensions in whichND argument:",
207             " is ", $which->dim(0), ", should be ", $ccs->[$WHICH]->dim(0));
208             }
209              
210             ##-- check for existing indices (potentially slow)
211 0         0 my $nzi = $ccs->indexNDi($which);
212 0         0 my ($nzi_new,$nzi_old) = ($nzi==$ccs->[$WHICH]->dim(1))->which_both;
213              
214             ##-- just set values for existing indices
215 0         0 $ccs->[$VALS]->index($nzi->index($nzi_old)) .= $vals->index($nzi_old);
216              
217             ##-- delegate insertion of new values to appendWhich()
218 0         0 my ($tmp);
219 0 0       0 return $tmp=$ccs->sortwhich if ($nzi_new->isempty);
220 0         0 return $tmp=$ccs->appendWhich($which->dice_axis(1,$nzi_new), $vals->index($nzi_new));
221             }
222              
223             ## $ccs = $ccs->appendWhich($whichND,$whichVals)
224             ## + inserts $whichND=>$whichVals into $ccs which are assumed NOT to be already present
225             ## + implicitly calls make_physically_indexed
226             sub appendWhich :lvalue {
227 0     0 1 0 my ($ccs,$which,$vals) = @_;
228 0         0 $ccs->make_physically_indexed();
229              
230             ##-- sanity check
231             #if ($which->dim(0) != $ccs->[$WHICH]->dim(0))
232 0 0       0 if ($which->dim(0) != $ccs->[$PDIMS]->nelem)
233             {
234 0         0 PDL::Lite::barf(ref($ccs)."::appendWhich(): wrong number of index dimensions in whichND argument:",
235             " is ", $which->dim(0), ", should be ", $ccs->[$PDIMS]->nelem);
236             }
237              
238             ##-- append: which
239 0 0       0 if (!$which->isempty) {
240 0         0 $ccs->[$WHICH] = $ccs->[$WHICH]->reshape($which->dim(0), $ccs->[$WHICH]->dim(1)+$which->dim(1));
241 0         0 $ccs->[$WHICH]->slice(",-".$which->dim(1).":-1") .= $which;
242             }
243              
244             ##-- append: vals
245 0 0       0 if (!$vals->isempty) {
246 0         0 my $missing = $ccs->missing;
247 0         0 $ccs->[$VALS] = $ccs->[$VALS]->reshape($ccs->[$VALS]->dim(0) + $vals->dim(0));
248 0         0 $ccs->[$VALS]->slice("-".($vals->dim(0)+1).":-2") .= $vals;
249 0         0 $ccs->[$VALS]->slice("-1") .= $missing;
250             }
251              
252 0         0 return $ccs->sortwhich();
253             }
254              
255             ## $ccs = $pdl->toccs()
256             ## $ccs = $pdl->toccs($missing)
257             ## $ccs = $pdl->toccs($missing,$flags)
258             *PDL::toccs = \&toccs;
259             sub toccs :lvalue {
260 0 0   0 1 0 return $_[0] if (isa($_[0],__PACKAGE__));
261 0         0 return my $tmp=__PACKAGE__->newFromDense(@_);
262             }
263              
264             ## $ccs = $ccs->copy()
265 4     4   10539 BEGIN { *clone = \© }
266             sub copy :lvalue {
267 0     0 1 0 my $ccs1 = shift;
268 0         0 my $ccs2 = bless [], ref($ccs1);
269 0         0 $ccs2->[$PDIMS] = $ccs1->[$PDIMS]->pdl;
270 0         0 $ccs2->[$VDIMS] = $ccs1->[$VDIMS]->pdl;
271 0         0 $ccs2->[$WHICH] = $ccs1->[$WHICH]->pdl;
272 0         0 $ccs2->[$VALS] = $ccs1->[$VALS]->pdl;
273 0 0       0 $ccs2->[$PTRS] = [ map {defined($_) ? [map {$_->pdl} @$_] : undef} @{$ccs1->[$PTRS]} ]; ##-- copy pointers?
  0         0  
  0         0  
  0         0  
274 0         0 $ccs2->[$FLAGS] = $ccs1->[$FLAGS];
275 0         0 return $ccs2;
276             }
277              
278             ## $ccs2 = $ccs->copyShallow()
279             ## + a very shallow version of copy()
280             ## + Copied : $PDIMS, @$PTRS, @{$PTRS->[*]}, $FLAGS
281             ## + Referenced: $VDIMS, $WHICH, $VALS, $PTRS->[*][*]
282             sub copyShallow :lvalue {
283 0     0 1 0 my $ccs = bless [@{$_[0]}], ref($_[0]);
  0         0  
284             ##
285             ##-- do copy some of it
286 0         0 $ccs->[$PDIMS] = $ccs->[$PDIMS]->pdl;
287             #$ccs->[$VDIMS] = $ccs->[$VDIMS]->pdl;
288 0 0       0 $ccs->[$PTRS] = [ map {defined($_) ? [@$_] : undef} @{$ccs->[$PTRS]} ];
  0         0  
  0         0  
289 0         0 $ccs;
290             }
291              
292             ## $ccs2 = $ccs->shadow(%args)
293             ## + args:
294             ## to => $ccs2, ##-- default: new
295             ## pdims => $pdims2, ##-- default: $pdims1->pdl (alias: 'dims')
296             ## vdims => $vdims2, ##-- default: $vdims1->pdl (alias: 'xdims')
297             ## ptrs => \@ptrs2, ##-- default: []
298             ## which => $which2, ##-- default: undef
299             ## vals => $vals2, ##-- default: undef ; if specified, should include final 'missing' element
300             ## flags => $flags, ##-- default: $flags1
301             sub shadow :lvalue {
302 0     0 1 0 my ($ccs,%args) = @_;
303 0 0 0     0 my $ccs2 = defined($args{to}) ? $args{to} : bless([], ref($ccs)||$ccs);
304 0 0       0 $ccs2->[$PDIMS] = (defined($args{pdims}) ? $args{pdims} : (defined($args{dims}) ? $args{dims} : $ccs->[$PDIMS]->pdl));
    0          
305 0 0       0 $ccs2->[$VDIMS] = (defined($args{vdims}) ? $args{vdims} : (defined($args{xdims}) ? $args{xdims} : $ccs->[$VDIMS]->pdl));
    0          
306 0 0       0 $ccs2->[$PTRS] = $args{ptrs} ? $args{ptrs} : [];
307 0         0 $ccs2->[$WHICH] = $args{which};
308 0         0 $ccs2->[$VALS] = $args{vals};
309 0 0       0 $ccs2->[$FLAGS] = defined($args{flags}) ? $args{flags} : $ccs->[$FLAGS];
310 0         0 return $ccs2;
311             }
312              
313              
314             ##--------------------------------------------------------------
315             ## Maintenance
316              
317             ## $ccs = $ccs->recode()
318             ## + recodes object, removing any missing values from $nzvals
319             sub recode :lvalue {
320 0     0 1 0 my $ccs = shift;
321 0         0 my $nz = $ccs->_nzvals;
322 0         0 my $z = $ccs->[$VALS]->slice("-1");
323              
324             ##-- get mask of "real" non-zero values
325 0         0 my $nzmask = PDL->zeroes($P_BYTE,$nz->nelem);
326 0         0 my ($nzmask1);
327 0 0       0 if ($z->isbad) {
328 0         0 $nz->isgood($nzmask);
329             }
330             else {
331 0         0 $nz->ne($z, $nzmask, 0);
332 0 0       0 if ($ccs->[$FLAGS] & $CCSND_BAD_IS_MISSING) {
333 0         0 $nzmask1 = $nzmask->pdl;
334 0         0 $nz->isgood($nzmask1);
335 0         0 $nzmask &= $nzmask1;
336             }
337             }
338 0 0       0 if ($ccs->[$FLAGS] & $CCSND_NAN_IS_MISSING) {
339 0 0       0 $nzmask1 = $nzmask->pdl if (!defined($nzmask1));
340 0         0 $nz->isfinite($nzmask1);
341 0         0 $nzmask &= $nzmask1;
342             }
343              
344             ##-- maybe recode
345 0 0       0 if (!$nzmask->all) {
346 0         0 my $nzi = $nzmask->which;
347 0         0 $ccs->[$WHICH] = $ccs->[$WHICH]->dice_axis(1,$nzi);
348 0         0 $ccs->[$VALS] = $ccs->[$VALS]->index($nzi)->append($z);
349 0         0 @{$ccs->[$PTRS]} = qw(); ##-- clear pointers
  0         0  
350             }
351              
352 0         0 return $ccs;
353             }
354              
355             ## $ccs = $ccs->sortwhich()
356             ## + sorts on $ccs->[$WHICH]
357             ## + may be DANGEROUS to indexing methods, b/c it alters $VALS
358             ## + clears pointers
359             sub sortwhich :lvalue {
360 0 0   0 1 0 return $_[0] if ($_[0][$WHICH]->isempty);
361 0         0 my $sorti = $_[0][$WHICH]->vv_qsortveci;
362 0         0 $_[0][$WHICH] = $_[0][$WHICH]->dice_axis(1,$sorti);
363 0         0 $_[0][$VALS] = $_[0][$VALS]->index($sorti->append($_[0][$WHICH]->dim(1)));
364             #
365             #-- DANGEROUS: pointer copy
366             # foreach (grep {defined($_)} @{$_[0][$PTRS]}) {
367             # $_->[1]->index($sorti) .= $_->[1];
368             # }
369             #--/DANGEROUS: pointer copy
370             #
371 0 0       0 @{$_[0][$PTRS]} = qw() if (! ($sorti==PDL->sequence($P_INDX,$sorti->dims))->all );
  0         0  
372 0         0 return $_[0];
373             }
374              
375              
376             ##--------------------------------------------------------------
377             ## Decoding
378              
379             ## $dense = $ccs->decode()
380             ## $dense = $ccs->decode($dense)
381             sub decode :lvalue {
382             ##-- decode physically stored index+value pairs
383 0     0 1 0 my $dense = ccs_decode($_[0][$WHICH],
384             $_[0]->_nzvals,
385             $_[0]->missing,
386             [ $_[0][$PDIMS] ],
387             );
388              
389             ##-- map physical dims with reorder()
390 0         0 my $porder = $_[0][$VDIMS]->where($_[0][$VDIMS]>=0);
391 0         0 $dense = $dense->reorder($porder->list); #if (($porder!=$_[0][$PDIMS]->sequence)->any);
392              
393             ##-- map virtual dims with dummy()
394 0         0 my @vdims = $_[0][$VDIMS]->list;
395 0         0 foreach (grep {$vdims[$_]<0} (0..$#vdims)) {
  0         0  
396 0         0 $dense = $dense->dummy($_, -$vdims[$_]);
397             }
398              
399             ##-- assign if $dense was specified by the user
400 0 0       0 if (defined($_[1])) {
401 0         0 $_[1] .= $dense;
402 0         0 return $_[1];
403             }
404              
405 0         0 return $dense;
406             }
407              
408             ## $dense = $ccs_or_dense->todense()
409             *PDL::todense = \&todense;
410 0 0   0 1 0 sub todense :lvalue { isa($_[0],__PACKAGE__) ? (my $tmp=$_[0]->decode(@_[1..$#_])) : $_[0]; }
411              
412             ##--------------------------------------------------------------
413             ## PDL API: Basic Properties
414              
415             ## $type = $obj->type()
416 0     0 0 0 sub type { $_[0][$VALS]->type; }
417              
418             ## $obj2 = $obj->convert($type)
419             ## + unlike PDL function, respects 'inplace' flag
420             sub convert :lvalue {
421 0 0   0 0 0 if ($_[0][$FLAGS] & $CCSND_INPLACE) {
422 0         0 $_[0][$VALS] = $_[0][$VALS]->convert($_[1]);
423 0         0 $_[0][$FLAGS] &= ~$CCSND_INPLACE;
424 0         0 return $_[0];
425             }
426 0         0 return my $tmp=$_[0]->shadow(which=>$_[0][$WHICH]->pdl, vals=>$_[0][$VALS]->convert($_[1]));
427             }
428              
429             ## byte,short,ushort,long,double,...
430             sub _pdltype_sub {
431 32     32   161 my $pdltype = shift;
432 32 0   0   308 return sub { return $pdltype if (!@_); convert(@_,$pdltype); };
  0         0  
  0         0  
433             }
434             foreach my $pdltype (map {$_->{convertfunc}} values %PDL::Types::typehash) {
435             #qw(byte short ushort long longlong indx float double)
436             eval "*${pdltype} = _pdltype_sub(PDL::${pdltype}());";
437             }
438              
439             ## $dimpdl = $obj->dimpdl()
440             ## + values in $dimpdl are negative for virtual dimensions
441             sub dimpdl :lvalue {
442 0     0 0 0 my $dims = $_[0][$VDIMS]->pdl;
443 0         0 my $physi = ($_[0][$VDIMS]>=0)->which;
444 0         0 (my $tmp=$dims->index($physi)) .= $_[0][$PDIMS]->index($_[0][$VDIMS]->index($physi));
445 0         0 return $dims;
446             }
447              
448             ## @dims = $obj->dims()
449 0     0 0 0 sub dims { $_[0]->dimpdl->abs->list; }
450              
451             ## $dim = $obj->dim($dimi)
452 0     0 0 0 sub dim { $_[0]->dimpdl->abs->at($_[1]); }
453             *getdim = \&dim;
454              
455             ## $ndims = $obj->ndims()
456 0     0 0 0 sub ndims { $_[0][$VDIMS]->nelem; }
457             *getndims = \&ndims;
458              
459             ## $nelem = $obj->nelem
460 0     0 0 0 sub nelem { $_[0]->dimpdl->abs->dprod; }
461              
462             ## $bool = $obj->isnull
463 0     0 0 0 sub isnull { $_[0][$VALS]->isnull; }
464              
465             ## $bool = $obj->isempty
466 0     0 0 0 sub isempty { $_[0]->nelem==0; }
467              
468             ##--------------------------------------------------------------
469             ## Low-level CCS access
470              
471             ## $bool = $ccs->is_physically_indexed()
472             ## + returns true iff only physical dimensions are present
473             sub is_physically_indexed {
474             (
475 0 0   0 1 0 $_[0][$VDIMS]->ndims==$_[0][$PDIMS]->ndims
476             &&
477             ($_[0][$VDIMS]==$_[0][$VDIMS]->sequence)->all
478             );
479             }
480              
481             ## $ccs2 = $ccs->to_physically_indexed()
482             ## + ensures that all non-missing elements are physically indexed
483             ## + just returns $ccs if all non-missing elements are already physically indexed
484             sub to_physically_indexed {
485 0 0   0 1 0 return $_[0] if ($_[0]->is_physically_indexed);
486 0         0 my $ccs = shift;
487 0         0 my $which = $ccs->whichND;
488 0         0 my $vals = $ccs->whichVals;
489 0         0 my $sorti = $which->vv_qsortveci;
490 0         0 return $ccs->shadow(
491             pdims=>$ccs->dimpdl->abs,
492             vdims=>$ccs->[$VDIMS]->sequence,
493             which=>$which->dice_axis(1,$sorti),
494             vals =>$vals->index($sorti)->append($ccs->missing),
495             )->sever;
496             }
497              
498             ## $ccs = $ccs->make_physically_indexed()
499             *make_physical = \&make_physically_indexed;
500             sub make_physically_indexed {
501 0 0   0 1 0 return $_[0] if ($_[0]->is_physically_indexed);
502 0         0 @{$_[0]} = @{$_[0]->to_physically_indexed};
  0         0  
  0         0  
503 0         0 return $_[0];
504             }
505              
506             ## $pdims = $obj->pdims()
507             ## $vdims = $obj->vdims()
508 0     0 1 0 sub pdims :lvalue { $_[0][$PDIMS]; }
509 0     0 1 0 sub vdims :lvalue { $_[0][$VDIMS]; }
510              
511              
512             ## $nelem_p = $obj->nelem_p : maximum number of physically addressable elements
513             ## $nelem_v = $obj->nelem_v : maximum number of virtually addressable elements
514 0     0 1 0 sub nelem_p { $_[0][$PDIMS]->dprod; }
515             *nelem_v = \&nelem;
516              
517             ## $v_per_p = $obj->_ccs_nvperp() : number of virtual elements per physical element
518 0     0   0 sub _ccs_nvperp { $_[0][$VDIMS]->where($_[0][$VDIMS]<0)->abs->dprod; }
519              
520             ## $nstored_p = $obj->nstored_p : actual number of physically stored elements
521             ## $nstored_v = $obj->nstored_v : actual number of physically+virtually stored elements
522 0     0 1 0 sub nstored_p { $_[0][$WHICH]->dim(1); }
523 0     0 1 0 sub nstored_v { $_[0][$WHICH]->dim(1) * $_[0]->_ccs_nvperp; }
524             *nstored = \&nstored_v;
525              
526              
527             ## $nnz = $obj->_nnz_p : returns actual $obj->[$VALS]->dim(0)-1
528             ## $nnz = $obj->_nnz_v : returns virtual $obj->[$VALS]->dim(0)-1
529 0     0   0 sub _nnz_p { $_[0][$VALS]->dim(0)-1; }
530 0     0   0 sub _nnz_v { ($_[0][$VALS]->dim(0)-1) * $_[0]->_ccs_nvperp; }
531             *_nnz = \&_nnz_v;
532              
533             ## $nmissing_p = $obj->nmissing_p()
534             ## $nmissing_v = $obj->nmissing_v()
535 0     0 1 0 sub nmissing_p { $_[0]->nelem_p - $_[0]->nstored_p; }
536 0     0 1 0 sub nmissing_v { $_[0]->nelem_v - $_[0]->nstored_v; }
537             *nmissing = \&nmissing_v;
538              
539             ## $bool = $obj->allmissing
540             ## + true if no non-missing values are stored
541 0     0 1 0 sub allmissing { $_[0][$VALS]->nelem <= 1; }
542              
543              
544             ## $missing = $obj->missing()
545             ## $missing = $obj->missing($missing)
546             sub missing {
547 0 0   0 1 0 $_[0][$VALS]->set(-1,$_[1]) if (@_>1);
548 0         0 $_[0][$VALS]->slice("-1");
549             }
550              
551             ## $obj = $obj->_missing($missingVal)
552             sub _missing :lvalue {
553 0 0   0   0 $_[0][$VALS]->set(-1,$_[1]) if (@_>1);
554 0         0 $_[0];
555             }
556              
557             ## $whichND_stored = $obj->_whichND()
558             ## $whichND_stored = $obj->_whichND($whichND)
559             sub _whichND :lvalue {
560 0 0   0   0 $_[0][$WHICH] = $_[1] if (@_>1);
561 0         0 $_[0][$WHICH];
562             }
563              
564             ## $_nzvals = $obj->_nzvals()
565             ## $_nzvals = $obj->_nzvals($nzvals)
566             ## + physical storage only
567 4     4   52865 BEGIN { *_whichVals = \&_nzvals; }
568             sub _nzvals :lvalue {
569 0     0   0 my ($tmp);
570 0 0       0 $_[0][$VALS]=$_[1]->append($_[0][$VALS]->slice("-1")) if (@_ > 1);
571 0 0       0 return $tmp=$_[0][$VALS]->index(PDL->null) if ($_[0][$VALS]->dim(0)<=1);
572 0         0 return $tmp=$_[0][$VALS]->slice("0:-2");
573             }
574              
575             ## $vals = $obj->_vals()
576             ## $vals = $obj->_vals($storedvals)
577             ## + physical storage only
578             sub _vals :lvalue {
579 0 0   0   0 $_[0][$VALS]=$_[1] if (@_ > 1);
580 0         0 $_[0][$VALS];
581             }
582              
583              
584             ## $ptr = $obj->ptr($dim_p); ##-- scalar context
585             ## ($ptr,$pi2nzi) = $obj->ptr($dim_p); ##-- list context
586             ## + returns cached value in $ccs->[$PTRS][$dim_p] if present
587             ## + caches value in $ccs->[$PTRS][$dim_p] otherwise
588             ## + $dim defaults to zero, for compatibility
589             ## + if $dim is zero, all($pi2nzi==sequence($obj->nstored))
590             ## + physical dimensions ONLY
591             sub ptr {
592 0     0 1 0 my ($ccs,$dim) = @_;
593 0 0       0 $dim = 0 if (!defined($dim));
594 0 0       0 $ccs->[$PTRS][$dim] = [$ccs->getptr($dim)] if (!$ccs->hasptr($dim));
595 0 0       0 return wantarray ? @{$ccs->[$PTRS][$dim]} : $ccs->[$PTRS][$dim][0];
  0         0  
596             }
597              
598             ## $bool = $obj->hasptr($dim_p)
599             ## + returns true iff $obj has a cached pointer for physical dim $dim_p
600             sub hasptr {
601 0     0 1 0 my ($ccs,$dim) = @_;
602 0 0       0 $dim = 0 if (!defined($dim));
603 0 0       0 return defined($ccs->[$PTRS][$dim]) ? scalar(@{$ccs->[$PTRS][$dim]}) : 0;
  0         0  
604             }
605              
606             ## ($ptr,$pi2nzi) = $obj->getptr($dim_p);
607             ## + as for ptr(), but does NOT cache anything, and does NOT check the cache
608             ## + physical dimensions ONLY
609 0     0 1 0 sub getptr { ccs_encode_pointers($_[0][$WHICH]->slice("($_[1]),"), $_[0][$PDIMS]->index($_[1])); }
610              
611             ## ($ptr,$pi2nzi) = $obj->setptr($dim_p, $ptr,$pi2nzi );
612             ## + low-level: set pointer for $dim_p
613             sub setptr {
614 0 0   0 0 0 if (UNIVERSAL::isa($_[2],'ARRAY')) {
615 0         0 $_[0][$PTRS][$_[1]] = $_[2];
616             } else {
617 0         0 $_[0][$PTRS][$_[1]] = [$_[2],$_[3]];
618             }
619 0         0 return $_[0]->ptr($_[1]);
620             }
621              
622             ## $obj = $obj->clearptrs()
623 0     0 1 0 sub clearptrs :lvalue { @{$_[0][$PTRS]}=qw(); return $_[0]; }
  0         0  
  0         0  
624              
625             ## $obj = $obj->clearptr($dim_p)
626             ## + low-level: clear pointer(s) for $dim_p
627             sub clearptr :lvalue {
628 0     0 1 0 my ($ccs,$dim) = @_;
629 0 0       0 return $ccs->clearptrs() if (!defined($dim));
630 0         0 $ccs->[$PTRS][$dim] = undef;
631 0         0 return $ccs;
632             }
633              
634             ## $flags = $obj->flags()
635             ## $flags = $obj->flags($flags)
636             ## + get local flags
637 0 0   0 1 0 sub flags { $_[0][$FLAGS] = $_[1] if (@_ > 1); $_[0][$FLAGS]; }
  0         0  
638              
639             ## $bool = $obj->bad_is_missing()
640             ## $bool = $obj->bad_is_missing($bool)
641             sub bad_is_missing {
642 0 0   0 1 0 if (@_ > 1) {
643 0 0       0 if ($_[1]) { $_[0][$FLAGS] |= $CCSND_BAD_IS_MISSING; }
  0         0  
644 0         0 else { $_[0][$FLAGS] &= ~$CCSND_BAD_IS_MISSING; }
645             }
646 0         0 $_[0][$FLAGS] & $CCSND_BAD_IS_MISSING;
647             }
648              
649             ## $obj = $obj->badmissing()
650 0     0 1 0 sub badmissing { $_[0][$FLAGS] |= $CCSND_BAD_IS_MISSING; $_[0]; }
  0         0  
651              
652             ## $bool = $obj->nan_is_missing()
653             ## $bool = $obj->nan_is_missing($bool)
654             sub nan_is_missing {
655 0 0   0 1 0 if (@_ > 1) {
656 0 0       0 if ($_[1]) { $_[0][$FLAGS] |= $CCSND_NAN_IS_MISSING; }
  0         0  
657 0         0 else { $_[0][$FLAGS] &= ~$CCSND_NAN_IS_MISSING; }
658             }
659 0         0 $_[0][$FLAGS] & $CCSND_NAN_IS_MISSING;
660             }
661              
662             ## $obj = $obj->nanmissing()
663 0     0 1 0 sub nanmissing { $_[0][$FLAGS] |= $CCSND_NAN_IS_MISSING; $_[0]; }
  0         0  
664              
665              
666             ## undef = $obj->set_inplace($bool)
667             ## + sets local inplace flag
668             sub set_inplace ($$) {
669 0 0   0 0 0 if ($_[1]) { $_[0][$FLAGS] |= $CCSND_INPLACE; }
  0         0  
670 0         0 else { $_[0][$FLAGS] &= ~$CCSND_INPLACE; }
671             }
672              
673             ## $bool = $obj->is_inplace()
674 0 0   0 0 0 sub is_inplace ($) { ($_[0][$FLAGS] & $CCSND_INPLACE) ? 1 : 0; }
675              
676             ## $obj = $obj->inplace()
677             ## + sets local inplace flag
678 0     0 0 0 sub inplace ($) { $_[0][$FLAGS] |= $CCSND_INPLACE; $_[0]; }
  0         0  
679              
680             ## $bool = $obj->badflag()
681             ## $bool = $obj->badflag($bool)
682             ## + wraps $obj->[$WHICH]->badflag, $obj->[$VALS]->badflag()
683             sub badflag {
684 0 0   0 0 0 if (@_ > 1) {
685 0         0 $_[0][$WHICH]->badflag($_[1]);
686 0         0 $_[0][$VALS]->badflag($_[1]);
687             }
688 0   0     0 return $_[0][$WHICH]->badflag || $_[0][$VALS]->badflag;
689             }
690              
691             ## $obj = $obj->sever()
692             ## + severs all sub-pdls
693             sub sever {
694 0     0 0 0 $_[0][$PDIMS]->sever;
695 0         0 $_[0][$VDIMS]->sever;
696 0         0 $_[0][$WHICH]->sever;
697 0         0 $_[0][$VALS]->sever;
698 0         0 foreach (grep {defined($_)} (@{$_[0][$PTRS]})) {
  0         0  
  0         0  
699 0         0 $_->[0]->sever;
700 0         0 $_->[1]->sever;
701             }
702 0         0 $_[0];
703             }
704              
705             ## \&code = _setbad_sub($pdlcode)
706             ## + returns a sub implementing setbadtoval(), setvaltobad(), etc.
707             sub _setbad_sub {
708 16     16   34 my $pdlsub = shift;
709             return sub {
710 0 0   0   0 if ($_[0]->is_inplace) {
711 0         0 $pdlsub->($_[0][$VALS]->inplace, @_[1..$#_]);
712 0         0 $_[0]->set_inplace(0);
713 0         0 return $_[0];
714             }
715 0         0 $_[0]->shadow(
716             which=>$_[0][$WHICH]->pdl,
717             vals=>$pdlsub->($_[0][$VALS],@_[1..$#_]),
718             );
719 16         219 };
720             }
721              
722             ## $obj = $obj->setnantobad()
723             foreach my $badsub (qw(setnantobad setbadtonan setbadtoval setvaltobad)) {
724             eval "*${badsub} = _setbad_sub(PDL->can('$badsub'));";
725             }
726              
727             ##--------------------------------------------------------------
728             ## Dimension Shuffling
729              
730             ## $ccs = $ccs->setdims_p(@dims)
731             ## + sets physical dimensions
732             *setdims = \&setdims_p;
733 0     0 1 0 sub setdims_p { $_[0][$PDIMS] = PDL->pdl($P_INDX,@_[1..$#_]); }
734              
735             ## $ccs2 = $ccs->dummy($vdim_index)
736             ## $ccs2 = $ccs->dummy($vdim_index, $vdim_size)
737             sub dummy :lvalue {
738 0     0 1 0 my ($ccs,$vdimi,$vdimsize) = @_;
739 0         0 my @vdims = $ccs->[$VDIMS]->list;
740 0 0       0 $vdimsize = 1 if (!defined($vdimsize));
741 0 0       0 $vdimi = 0 if (!defined($vdimi));
742 0 0       0 $vdimi = @vdims + $vdimi + 1 if ($vdimi < 0);
743 0 0       0 if ($vdimi < 0) {
744 0         0 PDL::Lite::barf(ref($ccs). "::dummy(): negative dimension number ", ($vdimi+@vdims), " exceeds number of dims ", scalar(@vdims));
745             }
746 0         0 splice(@vdims,$vdimi,0,-$vdimsize);
747 0         0 my $ccs2 = $ccs->copyShallow;
748 0         0 $ccs2->[$VDIMS] = PDL->pdl($P_INDX,\@vdims);
749 0         0 return $ccs2;
750             }
751              
752             ## $ccs2 = $ccs->reorder_pdl($vdim_index_pdl)
753             sub reorder_pdl :lvalue {
754 0     0 0 0 my $ccs2 = $_[0]->copyShallow;
755 0         0 $ccs2->[$VDIMS] = $ccs2->[$VDIMS]->index($_[1]);
756 0         0 $ccs2->[$VDIMS]->sever;
757 0         0 $ccs2;
758             }
759              
760             ## $ccs2 = $ccs->reorder(@vdim_list)
761 0     0 1 0 sub reorder :lvalue { $_[0]->reorder_pdl(PDL->pdl($P_INDX,@_[1..$#_])); }
762              
763             ## $ccs2 = $ccs->xchg($vdim1,$vdim2)
764             sub xchg :lvalue {
765 0     0 1 0 my $dimpdl = PDL->sequence($P_INDX,$_[0]->ndims);
766 0         0 my $tmp = $dimpdl->at($_[1]);
767 0         0 $dimpdl->set($_[1], $dimpdl->at($_[2]));
768 0         0 $dimpdl->set($_[2], $tmp);
769 0         0 return $tmp=$_[0]->reorder_pdl($dimpdl);
770             }
771              
772             ## $ccs2 = $ccs->mv($vDimFrom,$vDimTo)
773             sub mv :lvalue {
774 0     0 1 0 my ($d1,$d2) = @_[1,2];
775 0         0 my $ndims = $_[0]->ndims;
776 0 0       0 $d1 = $ndims+$d1 if ($d1 < 0);
777 0 0       0 $d2 = $ndims+$d2 if ($d2 < 0);
778 0 0       0 return my $tmp=$_[0]->reorder($d1 < $d2
779             ? ((0..($d1-1)), (($d1+1)..$d2), $d1, (($d2+1)..($ndims-1)))
780             : ((0..($d2-1)), $d1, ($d2..($d1-1)), (($d1+1)..($ndims-1)))
781             );
782             }
783              
784             ## $ccs2 = $ccs->transpose()
785             ## + always copies
786             sub transpose :lvalue {
787 0     0 1 0 my ($tmp);
788 0 0       0 if ($_[0]->ndims==1) {
789 0         0 return $tmp=$_[0]->dummy(0,1)->copy;
790             } else {
791 0         0 return $tmp=$_[0]->xchg(0,1)->copy;
792             }
793             }
794              
795              
796              
797             ##--------------------------------------------------------------
798             ## PDL API: Indexing
799              
800             sub slice { #:lvalue
801 0     0 0 0 PDL::Lite::barf(ref($_[0])."::slice() is not implemented yet (try dummy, dice_axis, indexND, etc.)");
802             }
803              
804             ## $nzi = $ccs->indexNDi($ndi)
805             ## + returns Nnz indices for virtual ND-index PDL $ndi
806             ## + index values in $ndi which are not present in $ccs are returned in $nzi as:
807             ## $ccs->[$WHICH]->dim(1) == $ccs->_nnz_p
808             sub indexNDi :lvalue {
809 0     0 1 0 my ($ccs,$ndi) = @_;
810             ##
811             ##-- get physical dims
812 0         0 my $dims = $ccs->[$VDIMS];
813 0         0 my $whichdimp = ($dims>=0)->which;
814 0         0 my $pdimi = $dims->index($whichdimp);
815             ##
816             #$ndi = $ndi->dice_axis(0,$whichdimp) ##-- BUG?!
817 0 0 0     0 $ndi = $ndi->dice_axis(0,$pdimi)
818             if ( $ndi->dim(0)!=$ccs->[$WHICH]->dim(0) || ($pdimi!=PDL->sequence($ccs->[$WHICH]->dim(0)))->any );
819             ##
820 0         0 my $foundi = $ndi->vsearchvec($ccs->[$WHICH]);
821 0         0 my $foundi_mask = ($ndi==$ccs->[$WHICH]->dice_axis(1,$foundi))->andover;
822 0         0 $foundi_mask->inplace->not;
823 0         0 (my $tmp=$foundi->where($foundi_mask)) .= $ccs->[$WHICH]->dim(1);
824 0         0 return $foundi;
825             }
826              
827             ## $vals = $ccs->indexND($ndi)
828 0     0 1 0 sub indexND :lvalue { my $tmp=$_[0][$VALS]->index($_[0]->indexNDi($_[1])); }
829              
830             ## $vals = $ccs->index2d($xi,$yi)
831 0     0 1 0 sub index2d :lvalue { my $tmp=$_[0]->indexND($_[1]->cat($_[2])->xchg(0,1)); }
832              
833             ## $nzi = $ccs->xindex1d($xi)
834             ## + nzi indices for dice_axis(0,$xi)
835             ## + physically indexed only
836             sub xindex1d :lvalue {
837 0     0 1 0 my ($ccs,$xi) = @_;
838 0         0 $ccs->make_physically_indexed;
839 0         0 my $nzi = $ccs->[$WHICH]->ccs_xindex1d($xi);
840 0         0 $nzi->sever;
841 0         0 return $nzi;
842             }
843              
844             ## $subset = $ccs->xsubset1d($xi)
845             ## + subset object like dice_axis(0,$xi) without $xi-renumbering
846             ## + returned object should participate in dataflow
847             ## + physically indexed only
848             sub xsubset1d :lvalue {
849 0     0 1 0 my ($ccs,$xi) = @_;
850 0         0 my $nzi = $ccs->xindex1d($xi);
851 0         0 return $ccs->shadow(which=>$ccs->[$WHICH]->dice_axis(1,$nzi),
852             vals =>$ccs->[$VALS]->index($nzi->append($ccs->_nnz)));
853             }
854              
855             ## $nzi = $ccs->pxindex1d($dimi,$xi)
856             ## + nzi indices for dice_axis($dimi,$xi), using ptr($dimi)
857             ## + physically indexed only
858             sub pxindex1d :lvalue {
859 0     0 1 0 my ($ccs,$dimi,$xi) = @_;
860 0         0 $ccs->make_physically_indexed();
861 0         0 my ($ptr,$pix) = $ccs->ptr($dimi);
862 0         0 my $xptr = $ptr->index($xi);
863 0         0 my $xlen = $ptr->index($xi+1) - $xptr;
864 0 0       0 my $nzi = defined($pix) ? $pix->index($xlen->rldseq($xptr))->qsort : $xlen->rldseq($xptr);
865 0         0 $nzi->sever;
866 0         0 return $nzi;
867             }
868              
869             ## $subset = $ccs->pxsubset1d($dimi,$xi)
870             ## + subset object like dice_axis($dimi,$xi) without $xi-renumbering, using ptr($dimi)
871             ## + returned object should participate in dataflow
872             ## + physically indexed only
873             sub pxsubset1d {
874 0     0 1 0 my ($ccs,$dimi,$xi) = @_;
875 0         0 my $nzi = $ccs->pxindex1d($dimi,$xi);
876 0         0 return $ccs->shadow(which=>$ccs->[$WHICH]->dice_axis(1,$nzi),
877             vals =>$ccs->[$VALS]->index($nzi->append($ccs->_nnz)));
878             }
879              
880             ## $nzi = $ccs->xindex2d($xi,$yi)
881             ## + returns nz-index piddle matching any index-pair in Cartesian product ($xi x $yi)
882             ## + caller object must be a ccs-encoded 2d matrix
883             ## + physically indexed only
884             sub xindex2d :lvalue {
885 0     0 1 0 my ($ccs,$xi,$yi) = @_;
886 0         0 $ccs->make_physically_indexed;
887 0         0 my $nzi = $ccs->[$WHICH]->ccs_xindex2d($xi,$yi);
888 0         0 $nzi->sever;
889 0         0 return $nzi;
890             }
891              
892             ## $subset = $ccs->xsubset2d($xi,$yi)
893             ## + returns a subset CCS object for all index-pairs in $xi,$yi
894             ## + caller object must be a ccs-encoded 2d matrix
895             ## + returned object should participate in dataflow
896             ## + physically indexed only
897             sub xsubset2d :lvalue {
898 0     0 1 0 my ($ccs,$xi,$yi) = @_;
899 0         0 my $nzi = $ccs->xindex2d($xi,$yi);
900 0         0 return $ccs->shadow(which=>$ccs->[$WHICH]->dice_axis(1,$nzi),
901             vals =>$ccs->[$VALS]->index($nzi->append($ccs->_nnz)));
902             }
903              
904             ## $vals = $ccs->index($flati)
905             sub index :lvalue {
906 0     0 1 0 my ($ccs,$i) = @_;
907 0         0 my $dummy = PDL->pdl(0)->slice(join(',', map {"*$_"} $ccs->dims));
  0         0  
908 0         0 my @coords = $dummy->one2nd($i);
909 0         0 my $ind = PDL->zeroes($P_INDX,$ccs->ndims,$i->dims);
910 0         0 my ($tmp);
911 0         0 ($tmp=$ind->slice("($_),")) .= $coords[$_] foreach (0..$#coords);
912 0         0 return $tmp=$ccs->indexND($ind);
913             }
914              
915             ## $ccs2 = $ccs->dice_axis($axis_v, $axisi)
916             ## + returns a new ccs object, should participate in dataflow
917             sub dice_axis :lvalue {
918 0     0 1 0 my ($ccs,$axis_v,$axisi) = @_;
919             ##
920             ##-- get
921 0         0 my $ndims = $ccs->ndims;
922 0 0       0 $axis_v = $ndims + $axis_v if ($axis_v < 0);
923 0 0 0     0 PDL::Lite::barf(ref($ccs)."::dice_axis(): axis ".($axis_v<0 ? ($axis_v+$ndims) : $axis_v)." out of range: should be 0<=dim<$ndims")
    0          
924             if ($axis_v < 0 || $axis_v >= $ndims);
925 0         0 my $axis = $ccs->[$VDIMS]->at($axis_v);
926 0 0       0 my $asize = $axis < 0 ? -$axis : $ccs->[$PDIMS]->at($axis);
927 0         0 $axisi = PDL->topdl($axisi);
928 0         0 my ($aimin,$aimax) = $axisi->minmax;
929 0 0       0 PDL::Lite::barf(ref($ccs)."::dice_axis(): invalid index $aimin (valid range 0..".($asize-1).")") if ($aimin < 0);
930 0 0       0 PDL::Lite::barf(ref($ccs)."::dice_axis(): invalid index $aimax (valid range 0..".($asize-1).")") if ($aimax >= $asize);
931             ##
932             ##-- check for virtual
933 0 0       0 if ($axis < 0) {
934             ##-- we're dicing a virtual axis: ok, but why?
935 0         0 my $naxisi = $axisi->nelem;
936 0         0 my $ccs2 = $ccs->copyShallow();
937 0         0 $ccs2->[$VDIMS] = $ccs->[$VDIMS]->pdl;
938 0         0 $ccs2->[$VDIMS]->set($axis_v, -$naxisi);
939 0         0 return $ccs2;
940             }
941             ##-- ok, we're dicing on a real axis
942 0         0 my ($ptr,$pi2nzi) = $ccs->ptr($axis);
943 0         0 my ($ptrix,$pi2nzix) = $ptr->ccs_decode_pointer($axisi);
944 0 0       0 my $nzix = defined($pi2nzi) ? $pi2nzi->index($pi2nzix) : $pi2nzix;
945 0         0 my $which = $ccs->[$WHICH]->dice_axis(1,$nzix);
946 0         0 $which->sever;
947 0 0       0 (my $tmp=$which->slice("($axis),")) .= $ptrix if (!$which->isempty); ##-- isempty() fix: v1.12
948 0         0 my $nzvals = $ccs->[$VALS]->index($nzix->append($ccs->[$WHICH]->dim(1)));
949             ##
950             ##-- construct output object
951 0         0 my $ccs2 = $ccs->shadow();
952 0         0 $ccs2->[$PDIMS]->set($axis, $axisi->nelem);
953 0         0 $ccs2->[$WHICH] = $which;
954 0         0 $ccs2->[$VALS] = $nzvals;
955             ##
956             ##-- sort output object (if not dicing on 0th dimension)
957 0 0       0 return $axis==0 ? $ccs2 : ($tmp=$ccs2->sortwhich());
958             }
959              
960             ## $onedi = $ccs->n2oned($ndi)
961             ## + returns a pseudo-index
962             sub n2oned :lvalue {
963 0     0 1 0 my $dimsizes = PDL->pdl($P_INDX,1)->append($_[0]->dimpdl->abs)->slice("0:-2")->cumuprodover;
964 0         0 return my $tmp=($_[1] * $dimsizes)->sumover;
965             }
966              
967             ## $whichND = $obj->whichND
968             ## + just returns the literal index PDL if possible: beware of dataflow!
969             ## + indices are NOT guaranteed to be returned in any surface-logical order,
970             ## although physically indexed dimensions should be sorted in physical-lexicographic order
971             sub whichND :lvalue {
972 0     0 1 0 my $vpi = ($_[0][$VDIMS]>=0)->which;
973 0         0 my ($wnd);
974 0 0       0 if ( $_[0][$VDIMS]->nelem==$_[0][$PDIMS]->nelem ) {
975 0 0       0 if (($_[0][$VDIMS]->index($vpi)==$_[0][$PDIMS]->sequence)->all) {
976             ##-- all literal & physically ordered
977 0         0 $wnd=$_[0][$WHICH];
978             }
979             else {
980             ##-- all physical, but shuffled
981 0         0 $wnd=$_[0][$WHICH]->dice_axis(0,$_[0][$VDIMS]->index($vpi));
982             }
983 0 0       0 return wantarray ? $wnd->xchg(0,1)->dog : $wnd;
984             }
985             ##-- virtual dims are in the game: construct output pdl
986 0         0 my $ccs = shift;
987 0         0 my $nvperp = $ccs->_ccs_nvperp;
988 0         0 my $nv = $ccs->nstored_v;
989 0         0 $wnd = PDL->zeroes($P_INDX, $ccs->ndims, $nv);
990 0         0 (my $tmp=$wnd->dice_axis(0,$vpi)->flat) .= $ccs->[$WHICH]->dummy(1,$nvperp)->flat;
991 0 0       0 if (!$wnd->isempty) {
992 0         0 my $nzi = PDL->sequence($P_INDX,$nv);
993 0         0 my @vdims = $ccs->[$VDIMS]->list;
994 0         0 my ($vdimi);
995 0         0 foreach (grep {$vdims[$#vdims-$_]<0} (0..$#vdims)) {
  0         0  
996 0         0 $vdimi = $#vdims-$_;
997 0         0 $nzi->modulo(-$vdims[$vdimi], $wnd->slice("($vdimi),"), 0);
998             }
999             }
1000 0 0       0 return wantarray ? $wnd->xchg(0,1)->dog : $wnd;
1001             }
1002              
1003             ## $whichVals = $ccs->whichVals()
1004             ## + returns $VALS corresponding to whichND() indices
1005             ## + beware of dataflow!
1006             sub whichVals :lvalue {
1007 0     0 1 0 my $vpi = ($_[0][$VDIMS]>=0)->which;
1008 0         0 my ($tmp);
1009 0 0       0 return $tmp=$_[0]->_nzvals() if ( $_[0][$VDIMS]->nelem==$_[0][$PDIMS]->nelem ); ##-- all physical
1010             ##
1011             ##-- virtual dims are in the game: construct output pdl
1012 0         0 return $tmp=$_[0]->_nzvals->slice("*".($_[0]->_ccs_nvperp))->flat;
1013             }
1014              
1015             ## $which = $obj->which()
1016             ## + not guaranteed to be returned in any meaningful order
1017 0     0 1 0 sub which :lvalue { my $tmp=$_[0]->n2oned(scalar $_[0]->whichND); }
1018              
1019             ## $val = $ccs->at(@index)
1020 0     0 1 0 sub at { $_[0]->indexND(PDL->pdl($P_INDX,@_[1..$#_]))->sclr; }
1021              
1022             ## $val = $ccs->set(@index,$value)
1023             sub set {
1024 0     0 1 0 my $foundi = $_[0]->indexNDi(PDL->pdl($P_INDX,@_[1..($#_-1)]));
1025 0 0       0 if ( ($foundi==$_[0][$WHICH]->dim(1))->any ) {
1026 0         0 carp(ref($_[0]).": cannot set() a missing value!")
1027             } else {
1028 0         0 (my $tmp=$_[0][$VALS]->index($foundi)) .= $_[$#_];
1029             }
1030 0         0 return $_[0];
1031             }
1032              
1033             ##--------------------------------------------------------------
1034             ## Mask Utilities
1035              
1036             ## $missing_mask = $ccs->ismissing()
1037             sub ismissing :lvalue {
1038 0     0 0 0 $_[0]->shadow(which=>$_[0][$WHICH]->pdl, vals=>$_[0]->_nzvals->zeroes->ccs_indx->append(1));
1039             }
1040              
1041             ## $nonmissing_mask = $ccs->ispresent()
1042             sub ispresent :lvalue {
1043 0     0 0 0 $_[0]->shadow(which=>$_[0][$WHICH]->pdl, vals=>$_[0]->_nzvals->ones->ccs_indx->append(0));
1044             }
1045              
1046             ##--------------------------------------------------------------
1047             ## Ufuncs
1048              
1049             ## $ufunc_sub = _ufuncsub($subname, \&ccs_accum_sub, $allow_bad_missing)
1050             sub _ufuncsub {
1051 56     56   133 my ($subname,$accumsub,$allow_bad_missing) = @_;
1052 56 50       144 PDL::Lite::barf(__PACKAGE__, "::_ufuncsub($subname): no underlying CCS accumulator func!") if (!defined($accumsub));
1053             return sub :lvalue {
1054 0     0   0 my $ccs = shift;
1055             ##
1056             ##-- preparation
1057 0         0 my $which = $ccs->whichND;
1058 0         0 my $vals = $ccs->whichVals;
1059 0         0 my $missing = $ccs->missing;
1060 0         0 my @dims = $ccs->dims;
1061 0         0 my ($which1,$vals1);
1062 0 0       0 if ($which->dim(0) <= 1) {
1063             ##-- flat sum
1064 0         0 $which1 = PDL->zeroes($P_INDX,1,$which->dim(1)); ##-- dummy
1065 0         0 $vals1 = $vals;
1066             } else {
1067 0         0 $which1 = $which->slice("1:-1,");
1068 0         0 my $sorti = $which1->vv_qsortveci;
1069 0         0 $which1 = $which1->dice_axis(1,$sorti);
1070 0         0 $vals1 = $vals->index($sorti);
1071             }
1072             ##
1073             ##-- guts
1074 0 0 0     0 my ($which2,$nzvals2) = $accumsub->($which1,$vals1,
1075             ($allow_bad_missing || $missing->isgood ? ($missing,$dims[0]) : (0,0))
1076             );
1077             ##
1078             ##-- get output pdl
1079 0         0 shift(@dims);
1080 0         0 my ($tmp);
1081 0 0       0 return $tmp=$nzvals2->squeeze if (!@dims); ##-- just a scalar: return a plain PDL
1082             ##
1083 0         0 my $newdims = PDL->pdl($P_INDX,\@dims);
1084 0         0 return $tmp=$ccs->shadow(
1085             pdims =>$newdims,
1086             vdims =>$newdims->sequence,
1087             which =>$which2,
1088             vals =>$nzvals2->append($missing->convert($nzvals2->type)),
1089             );
1090 56         569 };
1091             }
1092              
1093             foreach my $ufunc (
1094             qw(prod dprod sum dsum),
1095             qw(and or band bor),
1096             )
1097             {
1098             eval "*${ufunc}over = _ufuncsub('${ufunc}over', PDL::CCS::Ufunc->can('ccs_accum_${ufunc}'))";
1099             }
1100             foreach my $ufunc (qw(maximum minimum average))
1101             {
1102             eval "*${ufunc} = _ufuncsub('${ufunc}', PDL::CCS::Ufunc->can('ccs_accum_${ufunc}'))";
1103             }
1104              
1105             *nbadover = _ufuncsub('nbadover', PDL::CCS::Ufunc->can('ccs_accum_nbad'), 1);
1106             *ngoodover = _ufuncsub('ngoodover', PDL::CCS::Ufunc->can('ccs_accum_ngood'), 1);
1107             *nnz = _ufuncsub('nnz', PDL::CCS::Ufunc->can('ccs_accum_nnz'), 1);
1108              
1109             sub average_nz :lvalue {
1110 0     0 0 0 my $ccs = shift;
1111 0         0 return my $tmp=($ccs->sumover / $ccs->nnz);
1112             }
1113             #sub average {
1114             # my $ccs = shift;
1115             # my $missing = $ccs->missing;
1116             # return $ccs->sumover / $ccs->dim(0) if ($missing==0);
1117             # return ($ccs->sumover + (-$ccs->nnz+$ccs->dim(0))*$missing) / $ccs->dim(0);
1118             #}
1119              
1120 0 0   0 0 0 sub sum { my $z=$_[0]->missing; $_[0]->_nzvals->sum + ($z->isgood ? ($z->sclr * $_[0]->nmissing) : 0); }
  0         0  
1121 0 0   0 0 0 sub dsum { my $z=$_[0]->missing; $_[0]->_nzvals->dsum + ($z->isgood ? ($z->sclr * $_[0]->nmissing) : 0); }
  0         0  
1122 0 0   0 0 0 sub prod { my $z=$_[0]->missing; $_[0]->_nzvals->prod * ($z->isgood ? ($z->sclr ** $_[0]->nmissing) : 1); }
  0         0  
1123 0 0   0 0 0 sub dprod { my $z=$_[0]->missing; $_[0]->_nzvals->dprod * ($z->isgood ? ($z->sclr ** $_[0]->nmissing) : 1); }
  0         0  
1124 0     0 0 0 sub min { $_[0][$VALS]->min; }
1125 0     0 0 0 sub max { $_[0][$VALS]->max; }
1126 0     0 0 0 sub minmax { $_[0][$VALS]->minmax; }
1127              
1128 0 0   0 0 0 sub nbad { my $z=$_[0]->missing; $_[0]->_nzvals->nbad + ($z->isbad ? $_[0]->nmissing : 0); }
  0         0  
1129 0 0   0 0 0 sub ngood { my $z=$_[0]->missing; $_[0]->_nzvals->ngood + ($z->isgood ? $_[0]->nmissing : 0); }
  0         0  
1130              
1131 0     0 0 0 sub any { $_[0][$VALS]->any; }
1132 0     0 0 0 sub all { $_[0][$VALS]->all; }
1133              
1134              
1135             sub avg {
1136 0     0 0 0 my $z=$_[0]->missing;
1137 0         0 return ($_[0]->_nzvals->sum + ($_[0]->nelem-$_[0]->_nnz)*$z->sclr) / $_[0]->nelem;
1138             }
1139 0     0 0 0 sub avg_nz { $_[0]->_nzvals->avg; }
1140              
1141             sub isbad {
1142 0     0 0 0 my ($a,$out) = @_;
1143 0         0 return $a->shadow(which=>$a->[$WHICH]->pdl,vals=>$a->[$VALS]->isbad,to=>$out);
1144             }
1145             sub isgood {
1146 0     0 0 0 my ($a,$out) = @_;
1147 0         0 return $a->shadow(which=>$a->[$WHICH]->pdl,vals=>$a->[$VALS]->isgood,to=>$out);
1148             }
1149              
1150             ##--------------------------------------------------------------
1151             ## Index-Ufuncs
1152              
1153             sub _ufunc_ind_sub {
1154 8     8   20 my ($subname,$accumsub,$allow_bad_missing) = @_;
1155 8 50       28 PDL::Lite::barf(__PACKAGE__, "::_ufuncsub($subname): no underlying CCS accumulator func!") if (!defined($accumsub));
1156             return sub :lvalue {
1157 0     0   0 my $ccs = shift;
1158             ##
1159             ##-- preparation
1160 0         0 my $which = $ccs->whichND;
1161 0         0 my $vals = $ccs->whichVals;
1162 0         0 my $missing = $ccs->missing;
1163 0         0 my @dims = $ccs->dims;
1164 0         0 my ($which0,$which1,$vals1);
1165 0 0       0 if ($which->dim(0) <= 1) {
1166             ##-- flat X_ind
1167 0         0 $which0 = $which->slice("(0),");
1168 0         0 $which1 = PDL->zeroes($P_INDX,1,$which->dim(1)); ##-- dummy
1169 0         0 $vals1 = $vals;
1170             } else {
1171 0         0 my $sorti = $which->dice_axis(0, PDL->sequence($P_INDX,$which->dim(0))->rotate(-1))->vv_qsortveci;
1172 0         0 $which1 = $which->slice("1:-1,")->dice_axis(1,$sorti);
1173 0         0 $which0 = $which->slice("(0),")->index($sorti);
1174 0         0 $vals1 = $vals->index($sorti);
1175             }
1176             ##
1177             ##-- guts
1178 0 0 0     0 my ($which2,$nzvals2) = $accumsub->($which1,$vals1,
1179             ($allow_bad_missing || $missing->isgood ? ($missing,$dims[0]) : (0,0))
1180             );
1181             ##
1182             ##-- get output pdl
1183 0         0 shift(@dims);
1184 0         0 my $nzi2 = $nzvals2;
1185 0         0 my $nzi2_ok = ($nzvals2>=0);
1186 0         0 my ($tmp);
1187 0         0 ($tmp=$nzi2->where($nzi2_ok)) .= $which0->index($nzi2->where($nzi2_ok));
1188 0 0       0 return $tmp=$nzi2->squeeze if (!@dims); ##-- just a scalar: return a plain PDL
1189             ##
1190 0         0 my $newdims = PDL->pdl($P_INDX,\@dims);
1191 0         0 return $tmp=$ccs->shadow(
1192             pdims =>$newdims,
1193             vdims =>$newdims->sequence,
1194             which =>$which2,
1195             vals =>$nzi2->append(-1),
1196             );
1197 8         83 };
1198             }
1199              
1200             *maximum_ind = _ufunc_ind_sub('maximum_ind', PDL::CCS::Ufunc->can('ccs_accum_maximum_nz_ind'),1);
1201             *minimum_ind = _ufunc_ind_sub('minimum_ind', PDL::CCS::Ufunc->can('ccs_accum_minimum_nz_ind'),1);
1202              
1203             ##--------------------------------------------------------------
1204             ## Ufuncs: qsort (from CCS::Functions)
1205              
1206             ## ($which0,$nzVals0, $nzix,$nzenum, $whichOut) = $ccs->_qsort()
1207             ## ($which0,$nzVals0, $nzix,$nzenum, $whichOut) = $ccs->_qsort([o]nzix(NNz), [o]nzenum(Nnz))
1208             sub _qsort {
1209 0     0   0 my $ccs = shift;
1210 0         0 my $which0 = $ccs->whichND;
1211 0         0 my $nzvals0 = $ccs->whichVals;
1212 0         0 return ($which0,$nzvals0, ccs_qsort($which0->slice("1:-1,"),$nzvals0, $ccs->missing,$ccs->dim(0), @_));
1213             }
1214              
1215             ## $ccs_sorted = $ccs->qsort()
1216             ## $ccs_sorted = $ccs->qsort($ccs_sorted)
1217             sub qsort :lvalue {
1218 0     0 0 0 my $ccs = shift;
1219 0         0 my ($which0,$nzvals0,$nzix,$nzenum) = $ccs->_qsort();
1220 0         0 my $newdims = PDL->pdl($P_INDX,[$ccs->dims]);
1221 0         0 return my $tmp=$ccs->shadow(
1222             to => $_[0],
1223             pdims =>$newdims,
1224             vdims =>$newdims->sequence,
1225             which =>$nzenum->slice("*1,")->glue(0,$which0->slice("1:-1,")->dice_axis(1,$nzix)),
1226             vals =>$nzvals0->index($nzix)->append($ccs->missing),
1227             );
1228             }
1229              
1230             ## $ccs_sortedi = $ccs->qsorti()
1231             ## $ccs_sortedi = $ccs->qsorti($ccs_sortedi)
1232             sub qsorti :lvalue {
1233 0     0 0 0 my $ccs = shift;
1234 0         0 my ($which0,$nzvals0,$nzix,$nzenum) = $ccs->_qsort();
1235 0         0 my $newdims = PDL->pdl($P_INDX,[$ccs->dims]);
1236 0         0 return my $tmp=$ccs->shadow(
1237             to => $_[0],
1238             pdims =>$newdims,
1239             vdims =>$newdims->sequence,
1240             which =>$nzenum->slice("*1,")->glue(0,$which0->slice("1:-1,")->dice_axis(1,$nzix)),
1241             vals =>$which0->slice("(0),")->index($nzix)->append(-1),
1242             );
1243             }
1244              
1245             ##--------------------------------------------------------------
1246             ## Unary Operations
1247              
1248             ## $sub = _unary_op($opname,$pdlsub)
1249             sub _unary_op {
1250 36     36   82 my ($opname,$pdlsub) = @_;
1251             return sub :lvalue {
1252 0 0   0   0 if ($_[0]->is_inplace) {
1253 0         0 $pdlsub->($_[0][$VALS]->inplace);
1254 0         0 $_[0]->set_inplace(0);
1255 0         0 return $_[0];
1256             }
1257 0         0 return my $tmp=$_[0]->shadow(which=>$_[0][$WHICH]->pdl, vals=>$pdlsub->($_[0][$VALS]));
1258 36         377 };
1259             }
1260              
1261             foreach my $unop (qw(bitnot sqrt abs sin cos not exp log log10))
1262             {
1263             eval "*${unop} = _unary_op('${unop}',PDL->can('${unop}'));";
1264             }
1265              
1266             ##--------------------------------------------------------------
1267             ## OLD (but still used): Binary Operations: missing-is-annihilator
1268              
1269             ## ($rdpdl,$pdimsc,$vdimsc,$apcp,$bpcp) = _ccsnd_binop_align_dims($pdimsa,$vdimsa, $pdimsb,$vdimsb, $opname)
1270             # + returns:
1271             ## $rdpdl : (indx,2,$nrdims) : [ [$vdimai,$vdimbi], ...] s.t. $vdimai should align with $vdimbi
1272             ## $pdimsc : (indx,$ndimsc) : physical dim-size pdl for CCS output $c()
1273             ## $vdimsc : (indx,$ndimsc) : virtual dim-size pdl for CCS output $c()
1274             ## $apcp : (indx,2,$nac) : [ [$apdimi,$cpdimi], ... ] s.t. $cpdimi aligns 1-1 with $apdimi
1275             ## $bpcp : (indx,2,$nbc) : [ [$bpdimi,$cpdimi], ... ] s.t. $cpdimi aligns 1-1 with $bpdimi
1276             sub _ccsnd_binop_align_dims {
1277 0     0   0 my ($pdimsa,$vdimsa,$pdimsb,$vdimsb, $opname) = @_;
1278 0 0       0 $opname = '_ccsnd_binop_relevant_dims' if (!defined($opname));
1279              
1280             ##-- init
1281 0         0 my @pdimsa = $pdimsa->list;
1282 0         0 my @pdimsb = $pdimsb->list;
1283 0         0 my @vdimsa = $vdimsa->list;
1284 0         0 my @vdimsb = $vdimsb->list;
1285              
1286             ##-- get alignment-relevant dims
1287 0         0 my @rdims = qw();
1288 0         0 my ($vdima,$vdimb, $dimsza,$dimszb);
1289 0 0       0 foreach (0..($#vdimsa < $#vdimsb ? $#vdimsa : $#vdimsb)) {
1290 0         0 $vdima = $vdimsa[$_];
1291 0         0 $vdimb = $vdimsb[$_];
1292              
1293             ##-- get (virtual) dimension sizes
1294 0 0       0 $dimsza = $vdima>=0 ? $pdimsa[$vdima] : -$vdima;
1295 0 0       0 $dimszb = $vdimb>=0 ? $pdimsb[$vdimb] : -$vdimb;
1296              
1297             ##-- check for (virtual) size mismatch
1298 0 0 0     0 next if ($dimsza==1 || $dimszb==1); ##... ignoring (virtual) dims of size 1
1299 0 0       0 PDL::Lite::barf( __PACKAGE__ , "::$opname(): dimension size mismatch on dim($_): $dimsza != $dimszb")
1300             if ($dimsza != $dimszb);
1301              
1302             ##-- dims match: only align if both are physical
1303 0 0 0     0 push(@rdims, [$vdima,$vdimb]) if ($vdima>=0 && $vdimb>=0);
1304             }
1305 0         0 my $rdpdl = PDL->pdl($P_INDX,\@rdims);
1306              
1307             ##-- get output dimension sources
1308 0         0 my @_cdsrc = qw(); ##-- ( $a_or_b_for_dim0, ... )
1309 0 0       0 foreach (0..($#vdimsa > $#vdimsb ? $#vdimsa : $#vdimsb)) {
1310 0 0       0 push(@vdimsa, -1) if ($_ >= @vdimsa);
1311 0 0       0 push(@vdimsb, -1) if ($_ >= @vdimsb);
1312 0         0 $vdima = $vdimsa[$_];
1313 0         0 $vdimb = $vdimsb[$_];
1314 0 0       0 $dimsza = $vdima>=0 ? $pdimsa[$vdima] : -$vdima;
1315 0 0       0 $dimszb = $vdimb>=0 ? $pdimsb[$vdimb] : -$vdimb;
1316 0 0       0 if ($vdima>=0) {
    0          
1317 0 0       0 if ($vdimb>=0) { push(@_cdsrc, $dimsza>=$dimszb ? 0 : 1); } ##-- a:p, b:p --> c:p[max2(sz(a),sz(b))]
  0 0       0  
1318 0         0 else { push(@_cdsrc, 0); } ##-- a:p, b:v --> c:p[a]
1319             }
1320 0         0 elsif ($vdimb>=0) { push(@_cdsrc, 1); } ##-- a:v, b:p --> c:p[b]
1321 0 0       0 else { push(@_cdsrc, $dimsza>=$dimszb ? 0 : 1); } ##-- a:v, b:v --> c:v[max2(sz(a),sz(b))]
1322             }
1323 0         0 my $_cdsrcp = PDL->pdl($P_INDX,@_cdsrc);
1324              
1325             ##-- get c() dimension pdls
1326 0         0 my @pdimsc = qw();
1327 0         0 my @vdimsc = qw();
1328 0         0 my @apcp = qw(); ##-- ([$apdimi,$cpdimi], ...)
1329 0         0 my @bpcp = qw(); ##-- ([$bpdimi,$bpdimi], ...)
1330 0         0 foreach (0..$#_cdsrc) {
1331 0 0       0 if ($_cdsrc[$_]==0) { ##-- source(dim=$_) == a
1332 0 0       0 if ($vdimsa[$_]<0) { $vdimsc[$_]=$vdimsa[$_]; }
  0         0  
1333             else {
1334 0         0 $vdimsc[$_] = @pdimsc;
1335 0         0 push(@apcp, [$vdimsa[$_],scalar(@pdimsc)]);
1336 0         0 push(@pdimsc, $pdimsa[$vdimsa[$_]]);
1337             }
1338             } else { ##-- source(dim=$_) == b
1339 0 0       0 if ($vdimsb[$_]<0) { $vdimsc[$_]=$vdimsb[$_]; }
  0         0  
1340             else {
1341 0         0 $vdimsc[$_] = @pdimsc;
1342 0         0 push(@bpcp, [$vdimsb[$_],scalar(@pdimsc)]);
1343 0         0 push(@pdimsc, $pdimsb [$vdimsb[$_]]);
1344             }
1345             }
1346             }
1347 0         0 my $pdimsc = PDL->pdl($P_INDX,\@pdimsc);
1348 0         0 my $vdimsc = PDL->pdl($P_INDX,\@vdimsc);
1349 0         0 my $apcp = PDL->pdl($P_INDX,\@apcp);
1350 0         0 my $bpcp = PDL->pdl($P_INDX,\@bpcp);
1351              
1352 0         0 return ($rdpdl,$pdimsc,$vdimsc,$apcp,$bpcp);
1353             }
1354              
1355             ##-- OLD (but still used)
1356             ## \&code = _ccsnd_binary_op_mia($opName, \&pdlSub, $defType, $noSwap)
1357             ## + returns code for wrapping a builtin PDL binary operation \&pdlSub under the name "$opName"
1358             ## + $opName is just used for error reporting
1359             ## + $defType (if specified) is the default output type of the operation (e.g. PDL::long())
1360             sub _ccsnd_binary_op_mia {
1361 80     80   300 my ($opname,$pdlsub,$deftype,$noSwap) = @_;
1362              
1363             return sub :lvalue {
1364 0     0     my ($a,$b,$swap) = @_;
1365 0           my ($tmp);
1366 0 0         $swap=0 if (!defined($swap));
1367              
1368             ##-- check for & dispatch scalar operations
1369 0 0 0       if (!ref($b) || $b->nelem==1) {
1370 0 0         if ($a->is_inplace) {
1371 0 0         $pdlsub->($a->[$VALS]->inplace, todense($b), ($noSwap ? qw() : $swap));
1372 0           $a->set_inplace(0);
1373 0           return $tmp=$a->recode;
1374             }
1375 0 0         return $tmp=$a->shadow(
1376             which => $a->[$WHICH]->pdl,
1377             vals => $pdlsub->($a->[$VALS], todense($b), ($noSwap ? qw() : $swap))
1378             )->recode;
1379             }
1380              
1381             ##-- convert b to CCS
1382 0           $b = toccs($b);
1383              
1384             ##-- align dimensions & determine output sources
1385 0           my ($rdpdl,$pdimsc,$vdimsc,$apcp,$bpcp) = _ccsnd_binop_align_dims(@$a[$PDIMS,$VDIMS],
1386             @$b[$PDIMS,$VDIMS],
1387             $opname);
1388 0           my $nrdims = $rdpdl->dim(1);
1389              
1390             ##-- get & sort relevant indices, vals
1391 0           my $ixa = $a->[$WHICH];
1392 0           my $avals = $a->[$VALS];
1393 0           my $nixa = $ixa->dim(1);
1394 0           my $ra = $rdpdl->slice("(0)");
1395 0           my ($ixar,$avalsr);
1396 0 0         if ( $rdpdl->isempty ) {
    0          
1397             ##-- a: no relevant dims: align all pairs using a pseudo-dimension
1398 0           $ixar = PDL->zeroes($P_INDX, 1,$nixa);
1399 0           $avalsr = $avals;
1400             } elsif ( ($ra==PDL->sequence($P_INDX,$nrdims))->all ) {
1401             ##-- a: relevant dims are a prefix of physical dims, e.g. pre-sorted
1402 0 0         $ixar = $nrdims==$ixa->dim(0) ? $ixa : $ixa->slice("0:".($nrdims-1));
1403 0           $avalsr = $avals;
1404             } else {
1405 0           $ixar = $ixa->dice_axis(0,$ra);
1406 0 0         my $ixar_sorti = $ixar->isempty ? PDL->null->ccs_indx() : $ixar->vv_qsortveci;
1407 0           $ixa = $ixa->dice_axis(1,$ixar_sorti);
1408 0           $ixar = $ixar->dice_axis(1,$ixar_sorti);
1409 0           $avalsr = $avals->index($ixar_sorti);
1410             }
1411             ##
1412 0           my $ixb = $b->[$WHICH];
1413 0           my $bvals = $b->[$VALS];
1414 0           my $nixb = $ixb->dim(1);
1415 0           my $rb = $rdpdl->slice("(1)");
1416 0           my ($ixbr,$bvalsr);
1417 0 0         if ( $rdpdl->isempty ) {
    0          
1418             ##-- b: no relevant dims: align all pairs using a pseudo-dimension
1419 0           $ixbr = PDL->zeroes($P_INDX, 1,$nixb);
1420 0           $bvalsr = $bvals;
1421             } elsif ( ($rb==PDL->sequence($P_INDX,$nrdims))->all ) {
1422             ##-- b: relevant dims are a prefix of physical dims, e.g. pre-sorted
1423 0 0         $ixbr = $nrdims==$ixb->dim(0) ? $ixb : $ixb->slice("0:".($nrdims-1));
1424 0           $bvalsr = $bvals;
1425             } else {
1426 0           $ixbr = $ixb->dice_axis(0,$rb);
1427 0 0         my $ixbr_sorti = $ixbr->isempty ? PDL->null->ccs_indx() : $ixbr->vv_qsortveci;
1428 0           $ixb = $ixb->dice_axis(1,$ixbr_sorti);
1429 0           $ixbr = $ixbr->dice_axis(1,$ixbr_sorti);
1430 0           $bvalsr = $bvals->index($ixbr_sorti);
1431             }
1432              
1433              
1434             ##-- initialize: state vars
1435 0 0         my $blksz = $nixa > $nixb ? $nixa : $nixb;
1436 0 0 0       $blksz = $BINOP_BLOCKSIZE_MIN if ($BINOP_BLOCKSIZE_MIN && $blksz < $BINOP_BLOCKSIZE_MIN);
1437 0 0 0       $blksz = $BINOP_BLOCKSIZE_MAX if ($BINOP_BLOCKSIZE_MAX && $blksz > $BINOP_BLOCKSIZE_MAX);
1438 0           my $istate = PDL->zeroes($P_INDX,7); ##-- [ nnzai,nnzai_nxt, nnzbi,nnzbi_nxt, nnzci,nnzci_nxt, cmpval ]
1439 0           my $ostate = $istate->pdl;
1440              
1441             ##-- initialize: output vectors
1442 0           my $nzai = PDL->zeroes($P_INDX,$blksz);
1443 0           my $nzbi = PDL->zeroes($P_INDX,$blksz);
1444 0 0         my $nzc = PDL->zeroes((defined($deftype)
    0          
1445             ? $deftype
1446             : ($avals->type > $bvals->type
1447             ? $avals->type
1448             : $bvals->type)),
1449             $blksz);
1450 0           my $ixc = PDL->zeroes($P_INDX, $pdimsc->nelem, $blksz);
1451 0           my $nnzc = 0;
1452 0 0         my $zc = $pdlsub->($avals->slice("-1"), $bvals->slice("-1"), ($noSwap ? qw() : $swap))->convert($nzc->type);
1453 0           my $nanismissing = ($a->[$FLAGS]&$CCSND_NAN_IS_MISSING);
1454 0           my $badismissing = ($a->[$FLAGS]&$CCSND_BAD_IS_MISSING);
1455 0 0 0       $zc = $zc->setnantobad() if ($nanismissing && $badismissing);
1456 0 0         my $zc_isbad = $zc->isbad ? 1 : 0;
1457              
1458             ##-- block-wise variables
1459             ## + there are way too many of these...
1460 0           my ($nzai_prv,$nzai_pnx, $nzbi_prv,$nzbi_pnx, $nzci_prv,$nzci_pnx,$cmpval_prv);
1461 0           my ($nzai_cur,$nzai_nxt, $nzbi_cur,$nzbi_nxt, $nzci_cur,$nzci_nxt,$cmpval);
1462 0           my ($nzci_max, $blk_slice, $nnzc_blk,$nnzc_slice_blk);
1463 0           my ($nzai_blk,$nzbi_blk,$ixa_blk,$ixb_blk,$ixc_blk,$nzc_blk,$cimask_blk,$ciwhich_blk);
1464 0           my $nnzc_prev=0;
1465 0   0       do {
1466             ##-- align a block of data
1467 0           ccs_binop_align_block_mia($ixar,$ixbr,$istate, $nzai,$nzbi,$ostate);
1468              
1469             ##-- parse current alignment algorithm state
1470 0           ($nzai_prv,$nzai_pnx, $nzbi_prv,$nzbi_pnx, $nzci_prv,$nzci_pnx,$cmpval_prv) = $istate->list;
1471 0           ($nzai_cur,$nzai_nxt, $nzbi_cur,$nzbi_nxt, $nzci_cur,$nzci_nxt,$cmpval) = $ostate->list;
1472 0           $nzci_max = $nzci_cur-1;
1473              
1474 0 0         if ($nzci_max >= 0) {
1475             ##-- construct block output pdls: nzvals
1476 0           $blk_slice = "${nzci_prv}:${nzci_max}";
1477 0           $nzai_blk = $nzai->slice($blk_slice);
1478 0           $nzbi_blk = $nzbi->slice($blk_slice);
1479 0 0         $nzc_blk = $pdlsub->($avalsr->index($nzai_blk), $bvalsr->index($nzbi_blk), ($noSwap ? qw() : $swap));
1480              
1481             ##-- get indices of non-$missing c() values
1482 0 0 0       $cimask_blk = $zc_isbad || $nzc_blk->badflag ? $nzc_blk->isgood : ($nzc_blk!=$zc);
1483 0 0 0       $cimask_blk &= $nzc_blk->isgood if (!$zc_isbad && $badismissing);
1484 0 0         $cimask_blk &= $nzc_blk->isfinite if ($nanismissing);
1485 0 0         if ($cimask_blk->any) {
1486 0           $ciwhich_blk = $cimask_blk->which;
1487 0           $nzc_blk = $nzc_blk->index($ciwhich_blk);
1488              
1489 0           $nnzc_blk = $nzc_blk->nelem;
1490 0           $nnzc += $nnzc_blk;
1491 0           $nnzc_slice_blk = "${nnzc_prev}:".($nnzc-1);
1492              
1493             ##-- construct block output pdls: ixc
1494 0           $ixc_blk = $ixc->slice(",$nnzc_slice_blk");
1495 0 0         if (!$apcp->isempty) {
1496 0           $ixa_blk = $ixa->dice_axis(1,$nzai_blk->index($ciwhich_blk));
1497 0           ($tmp=$ixc_blk->dice_axis(0,$apcp->slice("(1),"))) .= $ixa_blk->dice_axis(0,$apcp->slice("(0),"));
1498             }
1499 0 0         if (!$bpcp->isempty) {
1500 0           $ixb_blk = $ixb->dice_axis(1,$nzbi_blk->index($ciwhich_blk));
1501 0           ($tmp=$ixc_blk->dice_axis(0,$bpcp->slice("(1),"))) .= $ixb_blk->dice_axis(0,$bpcp->slice("(0),"));
1502             }
1503              
1504             ##-- construct block output pdls: nzc
1505 0           ($tmp=$nzc->slice($nnzc_slice_blk)) .= $nzc_blk;
1506             }
1507             }
1508              
1509             ##-- possibly allocate for another block
1510 0 0 0       if ($nzai_cur < $nixa || $nzbi_cur < $nixb) {
1511 0           $nzci_nxt -= $nzci_cur;
1512 0           $nzci_cur = 0;
1513              
1514 0 0         if ($nzci_nxt+$blksz > $nzai->dim(0)) {
1515 0           $nzai = $nzai->reshape($nzci_nxt+$blksz);
1516 0           $nzbi = $nzbi->reshape($nzci_nxt+$blksz);
1517             }
1518 0           $ixc = $ixc->reshape($ixc->dim(0), $ixc->dim(1)+$nzai->dim(0));
1519 0           $nzc = $nzc->reshape($nzc->dim(0)+$nzai->dim(0));
1520              
1521 0           ($tmp=$istate) .= $ostate;
1522 0           $istate->set(4, $nzci_cur);
1523 0           $istate->set(5, $nzci_nxt);
1524             }
1525 0           $nnzc_prev = $nnzc;
1526              
1527             } while ($nzai_cur < $nixa || $nzbi_cur < $nixb);
1528              
1529             ##-- trim output pdls
1530 0 0         if ($nnzc > 0) {
1531             ##-- usual case: some values are non-missing
1532 0           $ixc = $ixc->slice(",0:".($nnzc-1));
1533 0           my $ixc_sorti = $ixc->vv_qsortveci;
1534 0           $nzc = $nzc->index($ixc_sorti)->append($zc->convert($nzc->type));
1535 0           $nzc->sever;
1536 0           $ixc = $ixc->dice_axis(1,$ixc_sorti);
1537 0           $ixc->sever;
1538             } else {
1539             ##-- pathological case: all values are "missing"
1540 0           $ixc = $ixc->dice_axis(1,PDL->pdl([]));
1541 0           $ixc->sever;
1542 0           $nzc = $zc->convert($zc->type);
1543             }
1544              
1545             ##-- set up final output object
1546 0           my $c = $a->shadow(
1547             pdims => $pdimsc,
1548             vdims => $vdimsc,
1549             which => $ixc,
1550             vals => $nzc,
1551             );
1552 0 0         if ($a->is_inplace) {
1553 0           @$a = @$c;
1554 0           $a->set_inplace(0);
1555 0           return $a;
1556             }
1557 0           return $c;
1558 80         1928 };
1559             }
1560              
1561             ##--------------------------------------------------------------
1562             ## NEW (but unused): Binary Operations: missing-is-annihilator: alignment
1563              
1564             ## \@parsed = _ccsnd_parse_signature($sig)
1565             ## \@parsed = _ccsnd_parse_signature($sig, $errorName)
1566             ## + parses a PDL-style signature
1567             ## + returned array has the form:
1568             ## ( $parsed_arg1, $parsed_arg2, ..., $parsed_argN )
1569             ## + where $parsed_arg$i =
1570             ## { name=>$argName, type=>$type, flags=>$flags, dims=>\@argDimNames, ... }
1571             ## + $flags is the string inside [] between type and arg name, if any
1572             sub _ccsnd_parse_signature {
1573 0     0     my ($sig,$errname) = @_;
1574 0 0         if ($sig =~ /^\s*\(/) {
1575             ##-- remove leading and trailing parentheses from signature
1576 0           $sig =~ s/^\s*\(\s*//;
1577 0           $sig =~ s/\s*\)\s*//;
1578             }
1579 0           my @args = ($sig =~ /[\s;]*([^\;]+)/g);
1580 0           my $parsed = [];
1581 0           my ($argName,$dimStr,$type,$flags,@dims);
1582 0           foreach (@args) {
1583 0           ($type,$flags) = ('','');
1584              
1585             ##-- check for type
1586 0 0         if ($_ =~ s/^\s*(byte|short|ushort|int|long|longlong|indx|float|double)\s*//) {
1587 0           $type = $1;
1588             }
1589              
1590             ##-- check for []-flags
1591 0 0         if ($_ =~ s/^\s*\[([^\]]*)\]\s*//g) {
1592 0           $flags = $1;
1593             }
1594              
1595             ##-- create output list: $argNumber=>{name=>$argName, dims=>[$dimNumber=>$dimName]}
1596 0 0         if ($_ =~ /^\s*(\S+)\s*\(([^\)]*)\)\s*$/) {
1597 0           ($argName,$dimStr) = ($1,$2);
1598 0 0         @dims = grep {defined($_) && $_ ne ''} split(/\,\s*/, $dimStr);
  0            
1599 0           push(@$parsed,{type=>$type,flags=>$flags,name=>$argName,dims=>[@dims]});
1600             } else {
1601 0 0         $errname = __PACKAGE__ . "::_ccsnd_parse_signature()" if (!defined($errname));
1602 0           die("${errname}: could not parse argument string '$_' for signature '$sig'");
1603             }
1604             }
1605 0           return $parsed;
1606             }
1607              
1608             ## \%dims = _ccsnd_align_dims(\@parsedSig, \@ccs_arg_pdls)
1609             ## \%dims = _ccsnd_align_dims(\@parsedSig, \@ccs_arg_pdls, $opName)
1610             ## + returns an dimension-alignment structure for @parsedSig with args @ccs_arg_pdls
1611             ## + returned %dims:
1612             ## ( $dimName => {size=>$dimSize, phys=>\@physical }, ... )
1613             ## - dim names "__thread_dim_${i}" are reserved
1614             ## - \@physical = [ [$argi,$pdimi_in_argi], ... ]
1615             sub _ccsnd_align_dims {
1616 0     0     my ($sig,$args,$opName) = @_;
1617 0 0         $opName = __PACKAGE__ . "::_ccsnd_align_dims()" if (!defined($opName));
1618              
1619             ##-- init: get virtual & physical dimension lists for arguments
1620 0           my @vdims = map { [$_->[$VDIMS]->list] } @$args;
  0            
1621 0           my @pdims = map { [$_->[$PDIMS]->list] } @$args;
  0            
1622              
1623             ##-- %dims = ($dimName => {size=>$dimSize, phys=>\@physical,... })
1624             ## + dim names "__thread_dim_${i}" are reserved
1625             ## + \@physical = [ [$argi,$pdimi], ... ]
1626 0           my %dims = map {($_=>undef)} map {@{$_->{dims}}} @$sig;
  0            
  0            
  0            
1627 0           my $nthreads = 0; ##-- number of threaded dims
1628              
1629             ##-- iterate over signature arguments, getting & checking dimension sizes
1630 0           my ($threadi, $argi,$arg_sig,$arg_ccs, $maxdim,$dimi,$pdimi,$dim_sig,$dim_ccs,$dimName, $dimsize,$isvdim);
1631 0           foreach $argi (0..$#$sig) {
1632 0           $arg_sig = $sig->[$argi];
1633 0           $arg_ccs = $args->[$argi];
1634              
1635             ##-- check for unspecified args
1636 0 0         if (!defined($arg_ccs)) {
1637 0 0         next if ($arg_sig->{flags} =~ /[ot]/); ##-- ... but not output or temporaries
1638 0           croak("$opName: argument '$arg_sig->{name}' not specified!");
1639             }
1640              
1641             ##-- reset thread counter
1642 0           $threadi=0;
1643              
1644             ##-- check dimension sizes
1645 0           $maxdim = _max2($#{$arg_sig->{dims}}, $#{$vdims[$argi]});
  0            
  0            
1646 0           foreach $dimi (0..$maxdim) {
1647 0 0         if (defined($dim_sig = $arg_sig->{dims}[$dimi])) {
1648             ##-- explicit dimension: name it
1649 0           $dimName = $dim_sig;
1650             } else {
1651 0           $dimName = "__thread_dim_".($threadi++);
1652             }
1653              
1654 0 0         if ($#{$vdims[$argi]} >= $dimi) {
  0            
1655 0           $pdimi = $vdims[$argi][$dimi];
1656 0 0         if ($pdimi >= 0) {
1657 0           $dimsize = $pdims[$argi][$pdimi];
1658 0           $isvdim = 0;
1659             } else {
1660 0           $dimsize = -$pdimi;
1661 0           $isvdim = 1;
1662             }
1663             } else {
1664 0           $dimsize = 1;
1665 0           $isvdim = 1;
1666             }
1667              
1668 0 0         if (!defined($dims{$dimName})) {
    0          
1669             ##-- new dimension
1670 0           $dims{$dimName} = { size=>$dimsize, phys=>[] };
1671             }
1672             elsif ($dims{$dimName}{size} != $dimsize) {
1673 0 0         if ($dims{$dimName}{size}==1) {
    0          
1674             ##-- ... we already had it, but as size=1 : override the stored size
1675 0           $dims{$dimName}{size} = $dimsize;
1676             }
1677             elsif ($dimsize != 1) {
1678             ##-- ... this is a non-trivial (size>1) dim which doesn't match: complain
1679 0           croak("$opName: size mismatch on dimension '$dimName' in argument '$arg_sig->{name}'",
1680             ": is $dimsize, should be $dims{$dimName}{size}");
1681             }
1682             }
1683 0 0         if (!$isvdim) {
1684             ##-- physical dim: add to alignment structure
1685 0           push(@{$dims{$dimName}{phys}}, [$argi,$pdimi]);
  0            
1686             }
1687             }
1688 0 0         $nthreads = $threadi if ($threadi > $nthreads);
1689             }
1690              
1691             ##-- check for undefined dims
1692 0           foreach (grep {!defined($dims{$_})} keys(%dims)) {
  0            
1693             #croak("$opName: cannot determine size for dimension '$_'");
1694             ##
1695             ##-- just set it to 1?
1696 0           $dims{$_} = {size=>1,phys=>[]};
1697             }
1698              
1699 0           return \%dims;
1700             }
1701              
1702             ##--------------------------------------------------------------
1703             ## Binary Operations: missing-is-annihilator: wrappers
1704              
1705             ##-- arithmetical & comparison operations
1706             foreach my $binop (
1707             qw(plus minus mult divide modulo power),
1708             qw(gt ge lt le eq ne spaceship),
1709             )
1710             {
1711             eval "*${binop} = *${binop}_mia = _ccsnd_binary_op_mia('${binop}',PDL->can('${binop}'));";
1712             die(__PACKAGE__, ": could not define binary operation $binop: $@") if ($@);
1713             }
1714              
1715             *pow = *pow_mia = _ccsnd_binary_op_mia('power',PDL->can('pow'),undef,1);
1716              
1717             ##-- integer-only operations
1718             foreach my $intop (
1719             qw(and2 or2 xor shiftleft shiftright),
1720             )
1721             {
1722             my $deftype = PDL->can($intop)->(PDL->pdl(0),PDL->pdl(0),0)->type->ioname;
1723             eval "*${intop} = *${intop}_mia = _ccsnd_binary_op_mia('${intop}',PDL->can('${intop}'),PDL::${deftype}());";
1724             die(__PACKAGE__, ": could not define integer operation $intop: $@") if ($@);
1725             }
1726              
1727             ## rassgn_mia($to,$from): binary assignment operation with missing-annihilator assumption
1728             ## + argument order is REVERSE of PDL 'assgn()' argument order
1729             *rassgn_mia = _ccsnd_binary_op_mia('rassgn', sub { PDL::assgn($_[1],$_[0]); $_[1]; });
1730              
1731             ## $to = $to->rassgn($from)
1732             ## + calls newFromDense() with $to flags if $from is dense
1733             ## + otherwise, copies $from to $to
1734             ## + argument order is REVERSED wrt PDL::assgn()
1735             sub rassgn :lvalue {
1736 0     0 0   my ($to,$from) = @_;
1737 0 0 0       if (!ref($from) || $from->nelem==1) {
1738             ##-- assignment from a scalar: treat the Nd object as a mask of available values
1739 0           (my $tmp=$to->[$VALS]) .= todense($from);
1740 0           return $to;
1741             }
1742 0 0         if (isa($from,__PACKAGE__)) {
1743             ##-- assignment from a CCS object: copy on a full dim match or an empty "$to"
1744 0           my $fromdimp = $from->dimpdl;
1745 0           my $todimp = $to->dimpdl;
1746 0 0 0       if ( $to->[$VALS]->dim(0)<=1 || $todimp->isempty || ($fromdimp==$todimp)->all ) {
      0        
1747 0           @$to = @{$from->copy};
  0            
1748 0           return $to;
1749             }
1750             }
1751             ##-- $from is something else: pass it on to 'rassgn_mia': effectively treat $to->[$WHICH] as a mask for $from
1752 0           $to->[$FLAGS] |= $CCSND_INPLACE;
1753 0           return my $tmp=$to->rassgn_mia($from);
1754             }
1755              
1756             ## $to = $from->assgn($to)
1757             ## + obeys PDL conventions
1758 0     0 0   sub assgn :lvalue { return my $tmp=$_[1]->rassgn($_[0]); }
1759              
1760              
1761             ##--------------------------------------------------------------
1762             ## CONTINUE HERE
1763              
1764             ## TODO:
1765             ## + virtual dimensions: clump
1766             ## + OPERATIONS:
1767             ## - accumulators: (some still missing: statistical, extrema-indices, atan2, ...)
1768              
1769             ##--------------------------------------------------------------
1770             ## Matrix operations
1771              
1772             ## $c = $a->inner($b)
1773             ## + inner product (may produce a large temporary)
1774 0     0 0   sub inner :lvalue { $_[0]->mult_mia($_[1],0)->sumover; }
1775              
1776             ## $c = $a->matmult($b)
1777             ## + mostly ganked from PDL::Primitive::matmult
1778             sub matmult :lvalue {
1779 0 0   0 0   PDL::Lite::barf("Invalid number of arguments for ", __PACKAGE__, "::matmult") if ($#_ < 1);
1780 0           my ($a,$b,$c) = @_; ##-- no $c!
1781 0 0 0       $c = undef if (!ref($c) && defined($c) && $c eq ''); ##-- strangeness: getting $c=''
      0        
1782              
1783 0           $b=toccs($b); ##-- ensure 2nd arg is a CCS object
1784              
1785             ##-- promote if necessary
1786 0           while ($a->getndims < 2) {$a = $a->dummy(-1)}
  0            
1787 0           while ($b->getndims < 2) {$b = $b->dummy(-1)}
  0            
1788              
1789             ##-- vector multiplication (easy)
1790 0 0 0       if ( ($a->dim(0)==1 && $a->dim(1)==1) || ($b->dim(0)==1 && $b->dim(1)==1) ) {
      0        
      0        
1791 0 0         if (defined($c)) { @$c = @{$a*$b}; return $c; }
  0            
  0            
  0            
1792 0           return $c=($a*$b);
1793             }
1794              
1795 0 0         if ($b->dim(1) != $a->dim(0)) {
1796 0           PDL::Lite::barf(sprintf("Dim mismatch in ", __PACKAGE__ , "::matmult of [%dx%d] x [%dx%d]: %d != %d",
1797             $a->dim(0),$a->dim(1),$b->dim(0),$b->dim(1),$a->dim(0),$b->dim(1)));
1798             }
1799              
1800 0           my $_c = $a->dummy(1)->inner($b->xchg(0,1)->dummy(2)); ##-- ye olde guttes
1801 0 0         if (defined($c)) { @$c = @$_c; return $c; }
  0            
  0            
1802              
1803 0           return $_c;
1804             }
1805              
1806             ## $c_dense = $a->matmult2d_sdd($b_dense)
1807             ## $c_dense = $a->matmult2d_sdd($b_dense, $zc)
1808             ## + signature as for PDL::Primitive::matmult()
1809             ## + should handle missing values correctly (except for BAD, inf, NaN, etc.)
1810             ## + see PDL::CCS::MatrixOps(3pm) for details
1811             sub matmult2d_sdd :lvalue {
1812 0     0 0   my ($a,$b,$c, $zc) = @_;
1813 0 0 0       $c = undef if (!ref($c) && defined($c) && $c eq ''); ##-- strangeness: getting $c=''
      0        
1814              
1815             ##-- promote if necessary
1816 0           while ($a->getndims < 2) {$a = $a->dummy(-1)}
  0            
1817 0           while ($b->getndims < 2) {$b = $b->dummy(-1)}
  0            
1818              
1819             ##-- vector multiplication (easy)
1820 0 0 0       if ( ($a->dim(0)==1 && $a->dim(1)==1) || ($b->dim(0)==1 && $b->dim(1)==1) ) {
      0        
      0        
1821 0 0         if (defined($c)) { @$c = @{$a*$b}; return $c; }
  0            
  0            
  0            
1822 0           return $c=($a*$b);
1823             }
1824              
1825             ##-- check dim sizes
1826 0 0         if ($b->dim(1) != $a->dim(0)) {
1827 0           PDL::Lite::barf(sprintf("Dim mismatch in ", __PACKAGE__, "::matmult2d [%dx%d] x [%dx%d] : %d != %d",
1828             $a->dims,$b->dims, $a->dim(0),$b->dim(1)));
1829             }
1830              
1831             ##-- ensure $b dense, $a physically indexed ccs
1832 0 0         $b = todense($b) if ($b->isa(__PACKAGE__));
1833 0           $a = $a->to_physically_indexed();
1834 0 0         if (!defined($c)) {
1835 0 0         my $ctype = $a->type > $b->type ? $a->type : $b->type;
1836 0           $c = PDL->zeroes($ctype, $b->dim(0),$a->dim(1));
1837             }
1838              
1839             ##-- compute $zc if required
1840 0 0         if (!defined($zc)) {
1841 0           $zc = (($a->missing + PDL->zeroes($a->type, $a->dim(0), 1)) x $b)->flat;
1842             }
1843              
1844 0           ccs_matmult2d_sdd($a->_whichND,$a->_nzvals,$a->missing, $b, $zc, $c);
1845              
1846 0           return $c;
1847             }
1848              
1849             ## $c_dense = $a->matmult2d_zdd($b_dense)
1850             ## + signature as for PDL::Primitive::matmult()
1851             ## + assumes $a->missing==0
1852             sub matmult2d_zdd :lvalue {
1853 0     0 0   my ($a,$b,$c) = @_;
1854 0 0 0       $c = undef if (!ref($c) && defined($c) && $c eq ''); ##-- strangeness: getting $c=''
      0        
1855              
1856             ##-- promote if necessary
1857 0           while ($a->getndims < 2) {$a = $a->dummy(-1)}
  0            
1858 0           while ($b->getndims < 2) {$b = $b->dummy(-1)}
  0            
1859              
1860             ##-- vector multiplication (easy)
1861 0 0 0       if ( ($a->dim(0)==1 && $a->dim(1)==1) || ($b->dim(0)==1 && $b->dim(1)==1) ) {
      0        
      0        
1862 0 0         if (defined($c)) { @$c = @{$a*$b}; return $c; }
  0            
  0            
  0            
1863 0           return $c=($a*$b);
1864             }
1865              
1866             ##-- check dim sizes
1867 0 0         if ($b->dim(1) != $a->dim(0)) {
1868 0           PDL::Lite::barf(sprintf("Dim mismatch in ", __PACKAGE__, "::matmult2d [%dx%d] x [%dx%d] : %d != %d",
1869             $a->dims,$b->dims, $a->dim(0),$b->dim(1)));
1870             }
1871              
1872             ##-- ensure $b dense, $a physically indexed ccs
1873 0 0         $b = todense($b) if ($b->isa(__PACKAGE__));
1874 0           $a = $a->to_physically_indexed();
1875 0 0         if (!defined($c)) {
1876 0 0         my $ctype = $a->type > $b->type ? $a->type : $b->type;
1877 0           $c = PDL->zeroes($ctype, $b->dim(0),$a->dim(1));
1878             }
1879 0           ccs_matmult2d_zdd($a->_whichND,$a->_nzvals, $b, $c);
1880              
1881 0           return $c;
1882             }
1883              
1884             ## $vnorm_dense = $a->vnorm($pdimi, ?$vnorm_dense)
1885             ## + assumes $a->missing==0
1886             sub vnorm {
1887 0     0 0   my ($a,$pdimi,$vnorm) = @_;
1888              
1889             ##-- ensure $a physically indexed ccs, $vnorm defined
1890 0           $a = $a->to_physically_indexed();
1891 0 0         $vnorm = PDL->zeroes($a->type, $a->dim($pdimi)) if (!defined($vnorm));
1892              
1893 0           ccs_vnorm($a->_whichND->slice("($pdimi),"), $a->_nzvals, $vnorm, $a->dim($pdimi));
1894 0           return $vnorm;
1895             }
1896              
1897              
1898             ## $vcos_dense = $a->vcos_zdd($b_dense, ?$vcos_dense, ?$norm_dense)
1899             ## + assumes $a->missing==0
1900             sub vcos_zdd {
1901 0     0 0   my $a = shift;
1902 0           my $b = shift;
1903              
1904             ##-- ensure $b dense, $a physically indexed ccs
1905 0 0         $b = todense($b) if (!UNIVERSAL::isa($b,__PACKAGE__));
1906 0           $a = $a->to_physically_indexed();
1907              
1908             ##-- guts
1909 0           return ccs_vcos_zdd($a->_whichND, $a->_nzvals, $b, $a->dim(0), @_);
1910             }
1911              
1912             ## $vcos_dense = $a->vcos_pzd($b_sparse, ?$norm_dense, ?$vcos_dense)
1913             ## + assumes $a->missing==0
1914             ## + uses $a->ptr(1)
1915             sub vcos_pzd {
1916 0     0 0   my $a = shift;
1917 0           my $b = shift;
1918              
1919             ##-- ensure $b dense, $a physically indexed ccs
1920 0 0         $b = toccs($b) if (!UNIVERSAL::isa($b,__PACKAGE__));
1921 0           $a = $a->to_physically_indexed();
1922 0           $b = $b->to_physically_indexed();
1923              
1924             ##-- get params
1925 0           my ($aptr,$aqsi) = $a->ptr(1);
1926 0           my $arows = $a->[$WHICH]->slice("(0),")->index($aqsi);
1927 0           my $avals = $a->[$VALS]->index($aqsi);
1928 0 0         my $anorm = @_ ? shift : $a->vnorm(0);
1929 0           my $brows = $b->[$WHICH]->slice("(0),");
1930 0           my $bvals = $b->_nzvals;
1931              
1932             ##-- guts
1933 0           return ccs_vcos_pzd($aptr,$arows,$avals, $brows,$bvals, $anorm, @_);
1934             }
1935              
1936              
1937             ##--------------------------------------------------------------
1938             ## Interpolation
1939              
1940             ## ($yi,$err) = $xi->interpolate($x,$y)
1941             ## + Signature: (xi(); x(n); y(n); [o] yi(); int [o] err())
1942             ## + routine for 1D linear interpolation
1943             ## + Given a set of points "($x,$y)", use linear interpolation to find the values $yi at a set of points $xi.
1944             ## + see PDL::Primitive::interpolate()
1945             sub interpolate {
1946 0     0 0   my ($xi,$x,$y, $yi,$err) = @_;
1947 0 0         $yi = $xi->clone if (!defined($yi));
1948 0 0         $err = $xi->clone if (!defined($err));
1949 0           $xi->[$VALS]->interpolate($x,$y, $yi->[$VALS], $err->[$VALS]);
1950 0 0         return wantarray ? ($yi,$err) : $yi;
1951             }
1952              
1953             ## $yi = $xi->interpolate($x,$y)
1954             ## + Signature: (xi(); x(n); y(n); [o] yi())
1955             ## + routine for 1D linear interpolation
1956             ## + see PDL::Primitive::interpol()
1957             sub interpol :lvalue {
1958 0     0 0   my ($xi,$x,$y, $yi) = @_;
1959 0 0         $yi = $xi->clone if (!defined($yi));
1960 0           $xi->[$VALS]->interpol($x,$y, $yi->[$VALS]);
1961 0           return $yi;
1962             }
1963              
1964              
1965             ##--------------------------------------------------------------
1966             ## General Information
1967              
1968             ## $density = $ccs->density()
1969             ## + returns PDL density as a scalar (lower is more sparse)
1970 0     0 1   sub density { $_[0][$WHICH]->dim(1) / $_[0]->nelem; }
1971              
1972             ## $compressionRate = $ccs->compressionRate()
1973             ## + higher is better
1974             ## + negative value indicates that dense storage would be more memory-efficient
1975             ## + pointers aren't included in the statistics: just which,nzvals,missing
1976             sub compressionRate {
1977 0     0 1   my $ccs = shift;
1978 0           my $dsize = PDL->pdl($ccs->nelem) * PDL::howbig($ccs->type);
1979 0           my $ccssize = (0
1980             + PDL->pdl($ccs->[$WHICH]->nelem) * PDL::howbig($ccs->[$WHICH]->type)
1981             + PDL->pdl($ccs->[$VALS]->nelem) * PDL::howbig($ccs->[$VALS]->type)
1982             + PDL->pdl($ccs->[$PDIMS]->nelem) * PDL::howbig($ccs->[$PDIMS]->type)
1983             + PDL->pdl($ccs->[$VDIMS]->nelem) * PDL::howbig($ccs->[$VDIMS]->type)
1984             );
1985 0           return (($dsize - $ccssize) / $dsize)->sclr;
1986             }
1987              
1988             ##--------------------------------------------------------------
1989             ## Stringification & Viewing
1990              
1991             ## $dimstr = _dimstr($pdl)
1992 0     0     sub _dimstr { return $_[0]->type.'('.join(',',$_[0]->dims).')'; }
1993 0     0     sub _pdlstr { return _dimstr($_[0]).'='.$_[0]; }
1994              
1995             ## $str = $obj->string()
1996             sub string {
1997 0     0 0   my ($pdims,$vdims,$which,$vals) = @{$_[0]}[$PDIMS,$VDIMS,$WHICH,$VALS];
  0            
1998 0 0         my $whichstr = ''.($which->isempty ? "Empty" : $which->xchg(0,1));
1999 0           $whichstr =~ s/^([^A-Z])/ $1/mg;
2000 0           chomp($whichstr);
2001             return
2002             (
2003 0           ''
2004             .ref($_[0]) . ':' . _dimstr($_[0]) ."\n"
2005             ." pdims:" . _pdlstr($pdims) ."\n"
2006             ." vdims:" . _pdlstr($vdims) ."\n"
2007             ." which:" . _dimstr($which)."^T=" . $whichstr . "\n"
2008             ." vals:" . _pdlstr($vals) ."\n"
2009             ." missing:" . _pdlstr($_[0]->missing) ."\n"
2010             );
2011             }
2012              
2013             ## $pstr = $obj->lstring()
2014             ## + literal perl-type string
2015 0     0 0   sub lstring { return overload::StrVal($_[0]); }
2016              
2017              
2018             ##======================================================================
2019             ## AUTOLOAD: pass to nonzero-PDL
2020             ## + doesn't seem to work well
2021             ##======================================================================
2022              
2023             #our $AUTOLOAD;
2024             #sub AUTOLOAD {
2025             # my $d = shift;
2026             # return undef if (!defined($d) || !defined($d->[$VALS]));
2027             # (my $name = $AUTOLOAD) =~ s/.*:://; ##-- strip qualification
2028             # my ($sub);
2029             # if (!($sub=UNIVERSAL::can($d->[$VALS],$name))) {
2030             # croak( ref($d) , "::$name() not defined for nzvals in ", __PACKAGE__ , "::AUTOLOAD.\n");
2031             # }
2032             # return $sub->($d->[$VALS],@_);
2033             #}
2034              
2035             ##--------------------------------------------------------------
2036             ## Operator overloading
2037              
2038             use overload (
2039             ##-- Binary ops: arithmetic
2040             "+" => \&plus_mia,
2041             "-" => \&minus_mia,
2042             "*" => \&mult_mia,
2043             "/" => \÷_mia,
2044             "%" => \&modulo_mia,
2045             "**" => \&power_mia,
2046 0     0   0 '+=' => sub { $_[0]->inplace->plus_mia(@_[1..$#_]); },
2047 0     0   0 '-=' => sub { $_[0]->inplace->minus_mia(@_[1..$#_]); },
2048 0     0   0 '*=' => sub { $_[0]->inplace->mult_mia(@_[1..$#_]); },
2049 0     0   0 '%=' => sub { $_[0]->inplace->divide_mia(@_[1..$#_]); },
2050 0     0   0 '**=' => sub { $_[0]->inplace->modulo_mia(@_[1..$#_]); },
2051              
2052             ##-- Binary ops: comparisons
2053             ">" => \>_mia,
2054             "<" => \<_mia,
2055             ">=" => \&ge_mia,
2056             "<=" => \&le_mia,
2057             "<=>" => \&spaceship_mia,
2058             "==" => \&eq_mia,
2059             "!=" => \&ne_mia,
2060             #"eq" => \&eq_mia
2061              
2062             ##-- Binary ops: bitwise & logic
2063             "|" => \&or2_mia,
2064             "&" => \&and2_mia,
2065             "^" => \&xor_mia,
2066             "<<" => \&shiftleft_mia,
2067             ">>" => \&shiftright_mia,
2068 0     0   0 '|=' => sub { $_[0]->inplace->or2_mia(@_[1..$#_]); },
2069 0     0   0 '&=' => sub { $_[0]->inplace->and2_mia(@_[1..$#_]); },
2070 0     0   0 '^=' => sub { $_[0]->inplace->xor_mia(@_[1..$#_]); },
2071 0     0   0 '<<=' => sub { $_[0]->inplace->shiftleft_mia(@_[1..$#_]); },
2072 0     0   0 '>>=' => sub { $_[0]->inplace->shiftright_mia(@_[1..$#_]); },
2073              
2074             ##-- Unary operations
2075             "!" => \¬,
2076             "~" => \&bitnot,
2077             "sqrt" => \&sqrt,
2078             "abs" => \&abs,
2079             "sin" => \&sin,
2080             "cos" => \&cos,
2081             "log" => \&log,
2082             "exp" => \&exp,
2083              
2084             ##-- assignment & assigning variants
2085             ".=" => \&rassgn,
2086              
2087             ##-- matrix operations
2088             'x' => \&matmult,
2089              
2090             ##-- Stringification & casts
2091             'bool' => sub {
2092 0     0   0 my $nelem = $_[0]->nelem;
2093 0 0       0 return 0 if ($nelem==0);
2094 0 0       0 croak("multielement ", __PACKAGE__, " pseudo-piddle in conditional expression") if ($nelem!=1);
2095 0         0 $_[0][$VALS]->at(0);
2096             },
2097 4         242 "\"\"" => \&string,
2098 4     4   55 );
  4         12  
2099              
2100              
2101             1; ##-- make perl happy
2102              
2103             ##======================================================================
2104             ## PODS: Header Administrivia
2105             ##======================================================================
2106             =pod
2107              
2108             =head1 NAME
2109              
2110             PDL::CCS::Nd - N-dimensional sparse pseudo-PDLs
2111              
2112             =head1 SYNOPSIS
2113              
2114             use PDL;
2115             use PDL::CCS::Nd;
2116              
2117             ##---------------------------------------------------------------------
2118             ## Example data
2119              
2120             $missing = 0; ##-- missing values
2121             $dense = random(@dims); ##-- densely encoded pdl
2122             $dense->where(random(@dims)<=0.95) .= $missing; ## ... made sparse
2123              
2124             $whichND = $dense->whichND; ##-- which values are present?
2125             $nzvals = $dense->indexND($whichND); ##-- ... and what are they?
2126              
2127              
2128             ##---------------------------------------------------------------------
2129             ## Constructors etc.
2130              
2131             $ccs = PDL::CCS:Nd->newFromDense($dense,%args); ##-- construct from dense matrix
2132             $ccs = PDL::CCS:Nd->newFromWhich($whichND,$nzvals,%args); ##-- construct from index+value pairs
2133              
2134             $ccs = $dense->toccs(); ##-- ensure PDL::CCS::Nd-hood
2135             $ccs = $ccs->toccs(); ##-- ... analogous to PDL::topdl()
2136             $ccs = $dense->toccs($missing,$flags); ##-- ... with optional arguments
2137              
2138             $ccs2 = $ccs->copy(); ##-- copy constructor
2139             $ccs2 = $ccs->copyShallow(); ##-- shallow copy, mainly for internal use
2140             $ccs2 = $ccs->shadow(%args); ##-- flexible copy method, for internal use
2141              
2142             ##---------------------------------------------------------------------
2143             ## Maintenance & Decoding
2144              
2145             $ccs = $ccs->recode(); ##-- remove missing values from stored VALS
2146             $ccs = $ccs->sortwhich(); ##-- internal use only
2147              
2148             $dense2 = $ccs->decode(); ##-- extract to a (new) dense matrix
2149              
2150             $dense2 = $ccs->todense(); ##-- ensure dense storage
2151             $dense2 = $dense2->todense(); ##-- ... analogous to PDL::topdl()
2152              
2153             ##---------------------------------------------------------------------
2154             ## PDL API: Basic Properties
2155              
2156             ##---------------------------------------
2157             ## Type conversion & Checking
2158             $ccs2 = $ccs->convert($type);
2159             $ccs2 = $ccs->byte;
2160             $ccs2 = $ccs->short;
2161             $ccs2 = $ccs->ushort;
2162             $ccs2 = $ccs->long;
2163             $ccs2 = $ccs->longlong;
2164             $ccs2 = $ccs->float;
2165             $ccs2 = $ccs->double;
2166              
2167             ##---------------------------------------
2168             ## Dimensions
2169             @dims = $ccs->dims();
2170             $ndims = $ccs->ndims();
2171             $dim = $ccs->dim($dimi);
2172             $nelem = $ccs->nelem;
2173             $bool = $ccs->isnull;
2174             $bool = $ccs->isempty;
2175              
2176             ##---------------------------------------
2177             ## Inplace & Dataflow
2178             $ccs = $ccs->inplace();
2179             $bool = $ccs->is_inplace;
2180             $bool = $ccs->set_inplace($bool);
2181             $ccs = $ccs->sever;
2182              
2183             ##---------------------------------------
2184             ## Bad Value Handling
2185              
2186             $bool = $ccs->bad_is_missing(); ##-- treat BAD values as missing?
2187             $bool = $ccs->bad_is_missing($bool);
2188             $ccs = $ccs->badmissing(); ##-- ... a la inplace()
2189              
2190             $bool = $ccs->nan_is_missing(); ##-- treat NaN values as missing?
2191             $bool = $ccs->nan_is_missing($bool);
2192             $ccs = $ccs->nanmissing(); ##-- ... a la inplace()
2193              
2194             $ccs2 = $ccs->setnantobad();
2195             $ccs2 = $ccs->setbadtonan();
2196             $ccs2 = $ccs->setbadtoval($val);
2197             $ccs2 = $ccs->setvaltobad($val);
2198              
2199             ##---------------------------------------------------------------------
2200             ## PDL API: Dimension Shuffling
2201              
2202             $ccs2 = $ccs->dummy($vdimi,$size);
2203             $ccs2 = $ccs->reorder(@vdims);
2204             $ccs2 = $ccs->xchg($vdim1,$vdim2);
2205             $ccs2 = $ccs->mv($vdimFrom,$vdimTo);
2206             $ccs2 = $ccs->transpose();
2207              
2208             ##---------------------------------------------------------------------
2209             ## PDL API: Indexing
2210              
2211             $nzi = $ccs->indexNDi($ndi); ##-- guts for indexing methods
2212             $ndi = $ccs->n2oned($ndi); ##-- returns 1d pseudo-index for $ccs
2213              
2214             $ivals = $ccs->indexND($ndi);
2215             $ivals = $ccs->index2d($xi,$yi);
2216             $ivals = $ccs->index($flati); ##-- buggy: no pseudo-threading!
2217             $ccs2 = $ccs->dice_axis($vaxis,$vaxis_ix);
2218              
2219             $nzi = $ccs->xindex1d($xi); ##-- nz-indices along 0th dimension
2220             $nzi = $ccs->pxindex1d($dimi,$xi); ##-- ... or any dimension, using ptr()
2221             $nzi = $ccs->xindex2d($xi,$yi); ##-- ... or for Cartesian product on 2d matrix
2222              
2223             $ccs2 = $ccs->xsubset1d($xi); ##-- subset along 0th dimension
2224             $ccs2 = $ccs->pxsubset1d($dimi,$xi); ##-- ... or any dimension, using ptr()
2225             $ccs2 = $ccs->xsubset2d($xi,$yi); ##-- ... or for Cartesian product on 2d matrix
2226              
2227             $whichND = $ccs->whichND();
2228             $vals = $ccs->whichVals(); ##-- like $ccs->indexND($ccs->whichND), but faster
2229             $which = $ccs->which()
2230              
2231             $value = $ccs->at(@index);
2232             $ccs = $ccs->set(@index,$value);
2233              
2234             ##---------------------------------------------------------------------
2235             ## PDL API: Ufuncs
2236              
2237             $ccs2 = $ccs->prodover;
2238             $ccs2 = $ccs->dprodover;
2239             $ccs2 = $ccs->sumover;
2240             $ccs2 = $ccs->dsumover;
2241             $ccs2 = $ccs->andover;
2242             $ccs2 = $ccs->orover;
2243             $ccs2 = $ccs->bandover;
2244             $ccs2 = $ccs->borover;
2245             $ccs2 = $ccs->maximum;
2246             $ccs2 = $ccs->minimum;
2247             $ccs2 = $ccs->maximum_ind; ##-- -1 indicates "missing" value is maximal
2248             $ccs2 = $ccs->minimum_ind; ##-- -1 indicates "missing" value is minimal
2249             $ccs2 = $ccs->nbadover;
2250             $ccs2 = $ccs->ngoodover;
2251             $ccs2 = $ccs->nnz;
2252              
2253             $sclr = $ccs->prod;
2254             $sclr = $ccs->dprod;
2255             $sclr = $ccs->sum;
2256             $sclr = $ccs->dsum;
2257             $sclr = $ccs->nbad;
2258             $sclr = $ccs->ngood;
2259             $sclr = $ccs->min;
2260             $sclr = $ccs->max;
2261             $bool = $ccs->any;
2262             $bool = $ccs->all;
2263              
2264             ##---------------------------------------------------------------------
2265             ## PDL API: Unary Operations (Overloaded)
2266              
2267             $ccs2 = $ccs->bitnot; $ccs2 = ~$ccs;
2268             $ccs2 = $ccs->not; $ccs2 = !$ccs;
2269             $ccs2 = $ccs->sqrt;
2270             $ccs2 = $ccs->abs;
2271             $ccs2 = $ccs->sin;
2272             $ccs2 = $ccs->cos;
2273             $ccs2 = $ccs->exp;
2274             $ccs2 = $ccs->log;
2275             $ccs2 = $ccs->log10;
2276              
2277             ##---------------------------------------------------------------------
2278             ## PDL API: Binary Operations (missing is annihilator)
2279             ## + $b may be a perl scalar, a dense PDL, or a PDL::CCS::Nd object
2280             ## + $c is always returned as a PDL::CCS::Nd ojbect
2281              
2282             ##---------------------------------------
2283             ## Arithmetic
2284             $c = $ccs->plus($b); $c = $ccs1 + $b;
2285             $c = $ccs->minus($b); $c = $ccs1 - $b;
2286             $c = $ccs->mult($b); $c = $ccs1 * $b;
2287             $c = $ccs->divide($b); $c = $ccs1 / $b;
2288             $c = $ccs->modulo($b); $c = $ccs1 % $b;
2289             $c = $ccs->power($b); $c = $ccs1 ** $b;
2290              
2291             ##---------------------------------------
2292             ## Comparisons
2293             $c = $ccs->gt($b); $c = ($ccs > $b);
2294             $c = $ccs->ge($b); $c = ($ccs >= $b);
2295             $c = $ccs->lt($b); $c = ($ccs < $b);
2296             $c = $ccs->le($b); $c = ($ccs <= $b);
2297             $c = $ccs->eq($b); $c = ($ccs == $b);
2298             $c = $ccs->ne($b); $c = ($ccs != $b);
2299             $c = $ccs->spaceship($b); $c = ($ccs <=> $b);
2300              
2301             ##---------------------------------------
2302             ## Bitwise Operations
2303             $c = $ccs->and2($b); $c = ($ccs & $b);
2304             $c = $ccs->or2($b); $c = ($ccs | $b);
2305             $c = $ccs->xor($b); $c = ($ccs ^ $b);
2306             $c = $ccs->shiftleft($b); $c = ($ccs << $b);
2307             $c = $ccs->shiftright($b); $c = ($ccs >> $b);
2308              
2309             ##---------------------------------------
2310             ## Matrix Operations
2311             $c = $ccs->inner($b);
2312             $c = $ccs->matmult($b); $c = $ccs x $b;
2313             $c_dense = $ccs->matmult2d_sdd($b_dense, $zc);
2314             $c_dense = $ccs->matmult2d_zdd($b_dense);
2315            
2316             $vnorm = $ccs->vnorm($pdimi);
2317             $vcos = $ccs->vcos_zdd($b_dense);
2318             $vcos = $ccs->vcos_pzd($b_ccs);
2319              
2320             ##---------------------------------------
2321             ## Other Operations
2322             $ccs->rassgn($b); $ccs .= $b;
2323             $str = $ccs->string(); $str = "$ccs";
2324              
2325             ##---------------------------------------------------------------------
2326             ## Indexing Utilities
2327              
2328              
2329             ##---------------------------------------------------------------------
2330             ## Low-Level Object Access
2331              
2332             $num_v_per_p = $ccs->_ccs_nvperp; ##-- num virtual / num physical
2333             $pdims = $ccs->pdims; $vdims = $ccs->vdims; ##-- physical|virtual dim pdl
2334             $nelem = $ccs->nelem_p; $nelem = $ccs->nelem_v; ##-- physical|virtual nelem
2335             $nstored = $ccs->nstored_p; $nstored = $ccs->nstored_v; ##-- physical|virtual Nnz+1
2336             $nmissing = $ccs->nmissing_p; $nmissing = $ccs->nmissing_v; ##-- physical|virtual nelem-Nnz
2337              
2338             $ccs = $ccs->make_physically_indexed(); ##-- ensure all dimensions are physically indexed
2339              
2340             $bool = $ccs->allmissing(); ##-- are all values missing?
2341              
2342             $missing_val = $ccs->missing; ##-- get missing value
2343             $missing_val = $ccs->missing($missing_val); ##-- set missing value
2344             $ccs = $ccs->_missing($missing_val); ##-- ... returning the object
2345              
2346             $whichND_phys = $ccs->_whichND(); ##-- get/set physical indices
2347             $whichND_phys = $ccs->_whichND($whichND_phys);
2348              
2349             $nzvals_phys = $ccs->_nzvals(); ##-- get/set physically indexed values
2350             $nzvals_phys = $ccs->_nzvals($vals_phys);
2351              
2352             $vals_phys = $ccs->_vals(); ##-- get/set physically indexed values
2353             $vals_phys = $ccs->_vals($vals_phys);
2354              
2355             $bool = $ccs->hasptr($pdimi); ##-- check for cached Harwell-Boeing pointer
2356             ($ptr,$ptrix) = $ccs->ptr($pdimi); ##-- ... get one, caching for later use
2357             ($ptr,$ptrix) = $ccs->getptr($pdimi); ##-- ... compute one, regardless of cache
2358             ($ptr,$ptrix) = $ccs->setptr($pdimi,$p,$pix); ##-- ... set a cached pointer
2359             $ccs->clearptr($pdimi); ##-- ... clear a cached pointer
2360             $ccs->clearptrs(); ##-- ... clear all cached pointers
2361              
2362             $flags = $ccs->flags(); ##-- get/set object-local flags
2363             $flags = $ccs->flags($flags);
2364              
2365             $density = $ccs->density; ##-- get object density
2366             $crate = $ccs->compressionRate; ##-- get compression rate
2367              
2368             =cut
2369              
2370             ##======================================================================
2371             ## Description
2372             ##======================================================================
2373             =pod
2374              
2375             =head1 DESCRIPTION
2376              
2377             PDL::CCS::Nd provides an object-oriented implementation of
2378             sparse N-dimensional vectors & matrices using a set of low-level
2379             PDLs to encode non-missing values.
2380             Currently, only a portion of the PDL API is implemented.
2381              
2382             =cut
2383              
2384             ##======================================================================
2385             ## Globals
2386             ##======================================================================
2387             =pod
2388              
2389             =head1 GLOBALS
2390              
2391             The following package-global variables are defined:
2392              
2393             =cut
2394              
2395             ##--------------------------------------------------------------
2396             ## Globals: Block Sizes
2397             =pod
2398              
2399             =head2 Block Size Constants
2400              
2401             $BINOP_BLOCKSIZE_MIN = 1;
2402             $BINOP_BLOCKSIZE_MAX = 0;
2403              
2404             Minimum (maximum) block size for block-wise incremental computation of binary operations.
2405             Zero or undef indicates no minimum (maximum).
2406              
2407             =cut
2408              
2409             ##--------------------------------------------------------------
2410             ## Globals: Object structure
2411             =pod
2412              
2413             =head2 Object Structure
2414              
2415             PDL::CCS::Nd object are implemented as perl ARRAY-references.
2416             For more intuitive access to object components, the following
2417             package-global variables can be used as array indices to access
2418             internal object structure:
2419              
2420             =over 4
2421              
2422             =item $PDIMS
2423              
2424             Indexes a pdl(long,$NPdims) of physically indexed dimension sizes:
2425              
2426             $ccs->[$PDIMS]->at($pdim_i) == $dimSize_i
2427              
2428             =item $VDIMS
2429              
2430             Indexes a pdl(long,$NVdims) of "virtual" dimension sizes:
2431              
2432             $ccs->[$VDIMS]->at($vdim_i) == / -$vdimSize_i if $vdim_i is a dummy dimension
2433             \ $pdim_i otherwise
2434              
2435             The $VDIMS piddle is used for dimension-shuffling transformations such as xchg()
2436             and reorder(), as well as for dummy().
2437              
2438             =item $WHICH
2439              
2440             Indexes a pdl(long,$NPdims,$Nnz) of the "physical indices" of all non-missing values
2441             in the non-dummy dimensions of the corresponding dense matrix.
2442             Vectors in $WHICH are guaranteed to be sorted in lexicographic order.
2443             If your $missing value is zero, and if your qsortvec() function works,
2444             it should be the case that:
2445              
2446             all( $ccs->[$WHICH] == $dense->whichND->qsortvec )
2447              
2448             A "physically indexed dimension" is just a dimension
2449             corresponding tp a single column of the $WHICH pdl, whereas a dummy dimension does
2450             not correspond to any physically indexed dimension.
2451              
2452             =item $VALS
2453              
2454             Indexes a vector pdl($valType, $Nnz+1) of all values in the sparse matrix,
2455             where $Nnz is the number of non-missing values in the sparse matrix. Non-final
2456             elements of the $VALS piddle are interpreted as the values of the corresponding
2457             indices in the $WHICH piddle:
2458              
2459             all( $ccs->[$VALS]->slice("0:-2") == $dense->indexND($ccs->[$WHICH]) )
2460              
2461             The final element of the $VALS piddle is referred to as "$missing", and
2462             represents the value of all elements of the dense physical matrix whose
2463             indices are not explicitly listed in the $WHICH piddle:
2464              
2465             all( $ccs->[$VALS]->slice("-1") == $dense->flat->index(which(!$dense)) )
2466              
2467             =item $PTRS
2468              
2469             Indexes an array of arrays containing Harwell-Boeing "pointer" piddle pairs
2470             for the corresponding physically indexed dimension.
2471             For a physically indexed dimension $d of size $N, $ccs-E[$PTRS][$d]
2472             (if it exists) is a pair [$ptr,$ptrix] as returned by
2473             PDL::CCS::Utils::ccs_encode_pointers($WHICH,$N), which are such that:
2474              
2475             =over 4
2476              
2477             =item $ptr
2478              
2479             $ptr is a pdl(long,$N+1) containing the offsets in $ptrix corresponding
2480             to the first non-missing value in the dimension $d.
2481             For all $i, 0 E= $i E $N, $ptr($i) contains the
2482             index of the first non-missing value (if any) from column $i of $dense(...,N,...)
2483             encoded in the $WHICH piddle. $ptr($N+1) contains the number of
2484             physically indexed cells in the $WHICH piddle.
2485              
2486             =item $ptrix
2487              
2488             Is an index piddle into dim(1) of $WHICH rsp. dim(0) of $VALS whose key
2489             positions correspond to the offsets listed in $ptr. The point here is
2490             that:
2491              
2492             $WHICH->dice_axis(1,$ptrix)
2493              
2494             is guaranteed to be primarily sorted along the pointer dimension $d, and
2495             stably sorted along all other dimensions, e.g. should be identical to:
2496              
2497             $WHICH->mv($d,0)->qsortvec->mv(0,$d)
2498              
2499             =back
2500              
2501              
2502             =item $FLAGS
2503              
2504             Indexes a perl scalar containing some object-local flags. See
2505             L<"Object Flags"> for details.
2506              
2507             =item $USER
2508              
2509             Indexes the first unused position in the object array.
2510             If you derive a class from PDL::CCS::Nd, you should use this
2511             position to place any new object-local data.
2512              
2513             =back
2514              
2515             =cut
2516              
2517              
2518             ##--------------------------------------------------------------
2519             ## Globals: Object Flags
2520             =pod
2521              
2522             =head2 Object Flags
2523              
2524             The following object-local constants are defined as bitmask flags:
2525              
2526             =over 4
2527              
2528             =item $CCSND_BAD_IS_MISSING
2529              
2530             Bitmask of the "bad-is-missing" flag. See the bad_is_missing() method.
2531              
2532             =item $CCSND_NAN_IS_MISSING
2533              
2534             Bitmask of the "NaN-is-missing" flag. See the nan_is_missing() method.
2535              
2536             =item $CCSND_INPLACE
2537              
2538             Bitmask of the "inplace" flag. See PDL::Core for details.
2539              
2540             =item $CCSND_FLAGS_DEFAULT
2541              
2542             Default flags for new objects.
2543              
2544             =back
2545              
2546             =cut
2547              
2548             ##======================================================================
2549             ## Methods
2550             ##======================================================================
2551             =pod
2552              
2553             =head1 METHODS
2554              
2555             =cut
2556              
2557             ##======================================================================
2558             ## Methods: Constructors etc.
2559             ##======================================================================
2560             =pod
2561              
2562             =head2 Constructors, etc.
2563              
2564             =over 4
2565              
2566             =item $class_or_obj-EnewFromDense($dense,$missing,$flags)
2567              
2568             =for sig
2569              
2570             Signature ($class_or_obj; dense(N1,...,NNdims); missing(); int flags)
2571              
2572             Class method. Create and return a new PDL::CCS::Nd object from a dense N-dimensional
2573             PDL $dense. If specified, $missing is used as the value for "missing" elements,
2574             and $flags are used to initialize the object-local flags.
2575              
2576             $missing defaults to BAD if the bad flag of $dense is set, otherwise
2577             $missing defaults to zero.
2578              
2579              
2580             =item $ccs-EfromDense($dense,$missing,$flags)
2581              
2582             =for sig
2583              
2584             Signature ($ccs; dense(N1,...,NNdims); missing(); int flags)
2585              
2586             Object method. Populate a sparse matrix object from a dense piddle $dense.
2587             See newFromDense().
2588              
2589              
2590             =item $class_or_obj-EnewFromWhich($whichND,$nzvals,%options)
2591              
2592             =for sig
2593              
2594             Signature ($class_or_obj; int whichND(Ndims,Nnz); nzvals(Nnz+1); %options)
2595              
2596             Class method. Create and return a new PDL::CCS::Nd object from a set
2597             of indices $whichND of non-missing elements in a (hypothetical) dense piddle
2598             and a vector $nzvals of the corresponding values. Known %options:
2599              
2600             sorted => $bool, ##-- if true, $whichND is assumed to be pre-sorted
2601             steal => $bool, ##-- if true, $whichND and $nzvals are used literally (formerly implied 'sorted')
2602             ## + in this case, $nzvals should really be: $nzvals->append($missing)
2603             pdims => $pdims, ##-- physical dimension list; default guessed from $whichND (alias: 'dims')
2604             vdims => $vdims, ##-- virtual dims (default: sequence($nPhysDims)); alias: 'xdims'
2605             missing => $missing, ##-- default: BAD if $nzvals->badflag, 0 otherwise
2606             flags => $flags ##-- flags
2607              
2608             =item $ccs-EfromWhich($whichND,$nzvals,%options)
2609              
2610             Object method. Guts for newFromWhich().
2611              
2612              
2613             =item $a-Etoccs($missing,$flags)
2614              
2615             Wrapper for newFromDense(). Return a PDL::CCS::Nd object for any piddle or
2616             perl scalar $a.
2617             If $a is already a PDL::CCS::Nd object, just returns $a.
2618             This method gets exported into the PDL namespace for ease of use.
2619              
2620              
2621             =item $ccs = $ccs-Ecopy()
2622              
2623             Full copy constructor.
2624              
2625             =item $ccs2 = $ccs-EcopyShallow()
2626              
2627             Shallow copy constructor, used e.g. by dimension-shuffling transformations.
2628             Copied components:
2629              
2630             $PDIMS, @$PTRS, @{$PTRS->[*]}, $FLAGS
2631              
2632             Referenced components:
2633              
2634             $VDIMS, $WHICH, $VALS, $PTRS->[*][*]
2635              
2636              
2637             =item $ccs2 = $ccs1-Eshadow(%args)
2638              
2639             Flexible constructor for computed PDL::CCS::Nd objects.
2640             Known %args:
2641              
2642             to => $ccs2, ##-- default: new
2643             pdims => $pdims2, ##-- default: $pdims1->pdl (alias: 'dims')
2644             vdims => $vdims2, ##-- default: $vdims1->pdl (alias: 'xdims')
2645             ptrs => \@ptrs2, ##-- default: []
2646             which => $which2, ##-- default: undef
2647             vals => $vals2, ##-- default: undef
2648             flags => $flags, ##-- default: $flags1
2649              
2650             =back
2651              
2652             =cut
2653              
2654              
2655             ##======================================================================
2656             ## Methods: Maintenance & Decoding
2657             ##======================================================================
2658             =pod
2659              
2660             =head2 Maintenance & Decoding
2661              
2662             =over 4
2663              
2664             =item $ccs = $ccs-Erecode()
2665              
2666             Recodes the PDL::CCS::Nd object, removing any missing values from its $VALS piddle.
2667              
2668             =item $ccs = $ccs-Esortwhich()
2669              
2670             Lexicographically sorts $ccs-E[$WHICH], altering $VALS accordingly.
2671             Clears $PTRS.
2672              
2673              
2674             =item $dense = $ccs-Edecode()
2675              
2676             =item $dense = $ccs-Edecode($dense)
2677              
2678             Decode a PDL::CCS::Nd object to a dense piddle.
2679             Dummy dimensions in $ccs should be created as dummy dimensions in $dense.
2680              
2681             =item $dense = $a-Etodense()
2682              
2683             Ensures that $a is not a PDL::CCS::Nd by wrapping decode().
2684             For PDLs or perl scalars, just returns $a.
2685              
2686             =back
2687              
2688             =cut
2689              
2690             ##======================================================================
2691             ## Methods: PDL API: Basic Properties
2692             ##======================================================================
2693             =pod
2694              
2695             =head2 PDL API: Basic Properties
2696              
2697             The following basic PDL API methods are implemented and/or wrapped
2698             for PDL::CCS::Nd objects:
2699              
2700             =over 4
2701              
2702             =item Type Checking & Conversion
2703              
2704             type, convert, byte, short, ushort, long, double
2705              
2706             Type-checking and conversion routines are passed on to the $VALS sub-piddle.
2707              
2708             =item Dimension Access
2709              
2710             dims, dim, getdim, ndims, getndims, nelem, isnull, isempty
2711              
2712             Note that nelem() returns the number of hypothetically addressable
2713             cells -- the number of cells in the corresponding dense matrix, rather
2714             than the number of non-missing elements actually stored.
2715              
2716             =item Inplace Operations
2717              
2718             set_inplace($bool), is_inplace(), inplace()
2719              
2720             =item Dataflow
2721              
2722             sever
2723              
2724             =item Bad Value Handling
2725              
2726             setnantobad, setbadtonan, setbadtoval, setvaltobad
2727              
2728             See also the bad_is_missing() and nan_is_missing() methods, below.
2729              
2730             =back
2731              
2732             =cut
2733              
2734             ##======================================================================
2735             ## Methods: PDL API: Dimension Shuffling
2736             ##======================================================================
2737             =pod
2738              
2739             =head2 PDL API: Dimension Shuffling
2740              
2741             The following dimension-shuffling methods are supported,
2742             and should be compatible to their PDL counterparts:
2743              
2744             =over 4
2745              
2746             =item dummy($vdimi)
2747              
2748             =item dummy($vdimi, $size)
2749              
2750             Insert a "virtual" dummy dimension of size $size at dimension index $vdimi.
2751              
2752              
2753             =item reorder(@vdim_list)
2754              
2755             Reorder dimensions according to @vdim_list.
2756              
2757             =item xchg($vdim1,$vdim2)
2758              
2759             Exchange two dimensions.
2760              
2761             =item mv($vdimFrom, $vdimTo)
2762              
2763             Move a dimension to another position, shoving remaining
2764             dimensions out of the way to make room.
2765              
2766             =item transpose()
2767              
2768             Always copies, unlike xchg(). Also unlike xchg(), works for 1d row-vectors.
2769              
2770             =back
2771              
2772             =cut
2773              
2774             ##======================================================================
2775             ## Methods: PDL API: Indexing
2776             ##======================================================================
2777             =pod
2778              
2779             =head2 PDL API: Indexing
2780              
2781             =over 4
2782              
2783             =item indexNDi($ndi)
2784              
2785             =for sig
2786              
2787             Signature: ($ccs; int ndi(NVdims,Nind); int [o]nzi(Nind))
2788              
2789             Guts for indexing methods. Given an N-dimensional index piddle $ndi, return
2790             a 1d index vector into $VALS for the corresponding values.
2791             Missing values are returned in $nzi as $Nnz == $ccs-E_nnz_p;
2792              
2793             Uses PDL::VectorValues::vsearchvec() internally, so expect O(Ndims * log(Nnz)) complexity.
2794             Although the theoretical complexity is tough to beat, this method could be
2795             made much faster in the usual (read "sparse") case by an intelligent use of $PTRS if
2796             and when available.
2797              
2798             =item indexND($ndi)
2799              
2800             =item index2d($xi,$yi)
2801              
2802             Should be mostly compatible to the PDL functions of the same names,
2803             but without any boundary handling.
2804              
2805             =item index($flati)
2806              
2807             Implicitly flattens the source pdl.
2808             This ought to be fixed.
2809              
2810             =item dice_axis($axis_v, $axisi)
2811              
2812             Should be compatible with the PDL function of the same name.
2813             Returns a new PDL::CCS::Nd object which should participate
2814             in dataflow.
2815              
2816             =item xindex1d($xi)
2817              
2818             Get non-missing indices for any element of $xi along 0th dimension;
2819             $xi must be sorted in ascending order.
2820              
2821             =item pxindex1d($dimi,$xi)
2822              
2823             Get non-missing indices for any element of $xi along physically indexed dimension $dimi,
2824             using L.
2825             $xi must be sorted in ascending order.
2826              
2827             =item xindex2d($xi,$yi)
2828              
2829             Get non-missing indices for any element in Cartesian product ($xi x $yi) for 2d sparse
2830             matrix.
2831             $xi and $yi must be sorted in ascending order.
2832              
2833             =item xsubset1d($xi)
2834              
2835             Returns a subset object similar to L,
2836             but without renumbering of indices along the diced dimension;
2837             $xi must be sorted in ascending order.
2838              
2839             =item pxsubset1d($dimi,$xi)
2840              
2841             Returns a subset object similar to L,
2842             but without renumbering of indices along the diced dimension;
2843             $xi must be sorted in ascending order.
2844              
2845             =item xsubset2d($xi,$yi)
2846              
2847             Returns a subset object similar to
2848             indexND( $xi-Eslice("*1,")-Ecat($yi)-Eclump(2)-Exchg(0,1) ),
2849             but without renumbering of indices;
2850             $xi and $yi must be sorted in ascending order.
2851              
2852              
2853             =item n2oned($ndi)
2854              
2855             Returns a 1d pseudo-index, used for implementation of which(), etc.
2856              
2857             =item whichND()
2858              
2859             Should behave mostly like the PDL function of the same name.
2860              
2861             Just returns the literal $WHICH piddle if possible: beware of dataflow!
2862             Indices are NOT guaranteed to be returned in any surface-logical order,
2863             although physically indexed dimensions should be sorted in physical-lexicographic
2864             order.
2865              
2866             =item whichVals()
2867              
2868             Returns $VALS indexed to correspond to the indices returned by whichND().
2869             The only reason to use whichND() and whichVals() rather than $WHICH and $VALS
2870             would be a need for physical representations of dummy dimension indices: try
2871             to avoid it if you can.
2872              
2873             =item which()
2874              
2875             As for the builtin PDL function.
2876              
2877              
2878             =item at(@index)
2879              
2880             Return a perl scalar corresponding to the Nd index @index.
2881              
2882             =item set(@index, $value)
2883              
2884             Set a non-missing value at index @index to $value.
2885             barf()s if @index points to a missing value.
2886              
2887             =back
2888              
2889             =cut
2890              
2891             ##======================================================================
2892             ## Methods: Operations: Ufuncs
2893             ##======================================================================
2894             =pod
2895              
2896             =head2 Ufuncs
2897              
2898             The following functions from PDL::Ufunc are implemented, and
2899             ought to handle missing values correctly (i.e. as their dense
2900             counterparts would):
2901              
2902             prodover
2903             prod
2904             dprodover
2905             dprod
2906             sumover
2907             sum
2908             dsumover
2909             dsum
2910             andover
2911             orover
2912             bandover
2913             borover
2914             maximum
2915             maximum_ind ##-- goofy if "missing" value is maximal
2916             max
2917             minimum
2918             minimum_ind ##-- goofy if "missing" value is minimal
2919             min
2920             nbadover
2921             nbad
2922             ngoodover
2923             ngood
2924             nnz
2925             any
2926             all
2927              
2928             Some Ufuncs are still unimplemented. see PDL::CCS::Ufunc for details.
2929              
2930             =cut
2931              
2932             ##======================================================================
2933             ## Methods: Operations: Unary
2934             ##======================================================================
2935             =pod
2936              
2937             =head2 Unary Operations
2938              
2939             The following unary operations are supported:
2940              
2941             FUNCTION OVERLOADS
2942             bitnot ~
2943             not !
2944             sqrt
2945             abs
2946             sin
2947             cos
2948             exp
2949             log
2950             log10
2951              
2952             Note that any pointwise unary operation can be performed directly on
2953             the $VALS piddle. You can wrap such an operation MY_UNARY_OP on piddles
2954             into a PDL::CCS::Nd method using the idiom:
2955              
2956             package PDL::CCS::Nd;
2957             *MY_UNARY_OP = _unary_op('MY_UNARY_OP', PDL->can('MY_UNARY_OP'));
2958              
2959             Note also that unary operations may change the "missing" value associated
2960             with the sparse matrix. This is easily seen to be the Right Way To Do It
2961             if you consider unary "not" over a very sparse (say 99% missing)
2962             binary-valued matrix: is is much easier and more efficient to alter only
2963             the 1% of physically stored (non-missing) values as well as the missing value
2964             than to generate a new matrix with 99% non-missing values, assuming $missing==0.
2965              
2966             =cut
2967              
2968             ##======================================================================
2969             ## Methods: Operations: Binary
2970             ##======================================================================
2971             =pod
2972              
2973             =head2 Binary Operations
2974              
2975             A number of basic binary operations on PDL::CCS::Nd operations are supported,
2976             which will produce correct results only under the assumption that "missing" values
2977             C<$missing> are annihilators for the operation in question. For example, if
2978             we want to compute:
2979              
2980             $c = OP($a,$b)
2981              
2982             for a binary operation OP on PDL::CCS::Nd objects C<$a> and C<$b>, the
2983             current implementation will produce the correct result for $c only if
2984             for all values C<$av> in C<$a> and C<$bv> in C<$b>:
2985              
2986             OP($av,$b->missing) == OP($a->missing,$b->missing) , and
2987             OP($a->missing,$bv) == OP($a->missing,$b->missing)
2988              
2989             This is true in general for OP==\&mult and $missing==0,
2990             but not e.g. for OP==\&plus and $missing==0.
2991             It should always hold for $missing==BAD (except in the case of assignment,
2992             which is a funny kind of operation anyways).
2993              
2994             Currently, the only way to ensure that all values are computed correctly
2995             in the general case is for $a and $b to contain exactly the same physically
2996             indexed values, which rather defeats the purposes of sparse storage,
2997             particularly if implicit pseudo-threading is involved (because then we would
2998             likely wind up instantiating -- or at least inspecting -- the entire dense
2999             matrix). Future implementations may relax these restrictions somewhat.
3000              
3001             The following binary operations are implemented:
3002              
3003             =over 4
3004              
3005             =item Arithmetic Operations
3006              
3007             FUNCTION OVERLOADS
3008             plus +
3009             minus -
3010             mult *
3011             divide /
3012             modulo %
3013             power **
3014              
3015             =item Comparisons
3016              
3017             FUNCTION OVERLOADS
3018             gt >
3019             ge >=
3020             lt <
3021             le <=
3022             eq ==
3023             ne !=
3024             spaceship <=>
3025              
3026             =item Bitwise Operations
3027              
3028             FUNCTION OVERLOADS
3029             and2 &
3030             or2 |
3031             xor ^
3032             shiftleft <<
3033             shiftright >>
3034              
3035             =item Matrix Operations
3036              
3037             FUNCTION OVERLOADS
3038             inner (none)
3039             matmult x
3040              
3041             =item Other Operations
3042              
3043             FUNCTION OVERLOADS
3044             rassgn .=
3045             string ""
3046              
3047             =back
3048              
3049             All supported binary operation functions obey the PDL input calling conventions
3050             (i.e. they all accept a third argument C<$swap>), and delegate computation
3051             to the underlying PDL functions. Note that the PDL::CCS::Nd methods currently
3052             do B support a third "output" argument.
3053             To wrap a new binary operation MY_BINOP into
3054             a PDL::CCS::Nd method, you can use the following idiom:
3055              
3056             package PDL::CCS::Nd;
3057             *MY_BINOP = _ccsnd_binary_op_mia('MY_BINOP', PDL->can('MY_BINOP'));
3058              
3059             The low-level alignment of physically indexed values
3060             for binary operations is performed by the
3061             function PDL::CCS::ccs_binop_align_block_mia().
3062             Computation is performed block-wise at the perl level to avoid
3063             over- rsp. underflow of the space requirements for the output PDL.
3064              
3065             =cut
3066              
3067              
3068             ##======================================================================
3069             ## Methods: Low-Level Object Access
3070             ##======================================================================
3071             =pod
3072              
3073             =head2 Low-Level Object Access
3074              
3075             The following methods provide low-level access to
3076             PDL::CCS::Nd object structure:
3077              
3078             =over 4
3079              
3080             =item insertWhich
3081              
3082             =for sig
3083              
3084             Signature: ($ccs; int whichND(Ndims,Nnz1); vals(Nnz1))
3085              
3086             Set or insert values in C<$ccs> for the indices in C<$whichND> to C<$vals>.
3087             C<$whichND> need not be sorted.
3088             Implicitly makes C<$ccs> physically indexed.
3089             Returns the (destructively altered) C<$ccs>.
3090              
3091              
3092             =item appendWhich
3093              
3094             =for sig
3095              
3096             Signature: ($ccs; int whichND(Ndims,Nnz1); vals(Nnz1))
3097              
3098             Like insertWhich(), but assumes that no values for any of the $whichND
3099             indices are already present in C<$ccs>. This is faster (because indexNDi
3100             need not be called), but less safe.
3101              
3102              
3103             =item is_physically_indexed()
3104              
3105             Returns true iff only physical dimensions are present.
3106              
3107             =item to_physically_indexed()
3108              
3109             Just returns the calling object if all non-missing elements are already physically indexed.
3110             Otherwise, returns a new PDL::CCS::Nd object identical to the caller
3111             except that all non-missing elements are physically indexed. This may gobble a large
3112             amount of memory if the calling element has large dummy dimensions.
3113             Also ensures that physical dimension order is identical to logical dimension order.
3114              
3115             =item make_physically_indexed
3116              
3117             Wrapper for to_physically_indexed() which eliminates dummy dimensions
3118             destructively in the calling object.
3119              
3120             Alias: make_physical().
3121              
3122              
3123             =item pdims()
3124              
3125             Returns the $PDIMS piddle. See L<"Object Structure">, above.
3126              
3127             =item vdims()
3128              
3129             Returns the $VDIMS piddle. See L<"Object Structure">, above.
3130              
3131              
3132             =item setdims_p(@dims)
3133              
3134             Sets $PDIMS piddle. See L<"Object Structure">, above.
3135             Returns the calling object.
3136             Alias: setdims().
3137              
3138             =item nelem_p()
3139              
3140             Returns the number of physically addressable elements.
3141              
3142             =item nelem_v()
3143              
3144             Returns the number of virtually addressable elements.
3145             Alias for nelem().
3146              
3147              
3148             =item _ccs_nvperp()
3149              
3150             Returns number of virtually addressable elements per physically
3151             adressable element, which should be a positive integer.
3152              
3153              
3154             =item nstored_p()
3155              
3156             Returns actual number of physically addressed stored elements
3157             (aka $Nnz aka $WHICH-Edim(1)).
3158              
3159             =item nstored_v()
3160              
3161             Returns actual number of physically+virtually addressed stored elements.
3162              
3163              
3164             =item nmissing_p()
3165              
3166             Returns number of physically addressable elements minus the number of
3167             physically stored elements.
3168              
3169             =item nmissing_v()
3170              
3171             Returns number of physically+virtually addressable elements minus the number of
3172             physically+virtually stored elements.
3173              
3174              
3175             =item allmissing()
3176              
3177             Returns true iff no non-missing values are stored.
3178              
3179              
3180             =item missing()
3181              
3182             =item missing($missing)
3183              
3184             Get/set the value to use for missing elements.
3185             Returns the (new) value for $missing.
3186              
3187              
3188             =item _whichND()
3189              
3190             =item _whichND($whichND)
3191              
3192             Get/set the underlying $WHICH piddle.
3193              
3194              
3195             =item _nzvals()
3196              
3197             =item _nzvals($storedvals)
3198              
3199             Get/set the slice of the underlying $VALS piddle correpsonding for non-missing values only.
3200             Alias: whichVals().
3201              
3202              
3203             =item _vals()
3204              
3205             =item _vals($storedvals)
3206              
3207             Get/set the underlying $VALS piddle.
3208              
3209             =item hasptr($pdimi)
3210              
3211             Returns true iff a pointer for physical dim $pdimi is cached.
3212              
3213             =item ptr($pdimi)
3214              
3215             Get a pointer pair for a physically indexed dimension $pdimi.
3216             Uses cached piddles in $PTRS if present, computes & caches otherwise.
3217              
3218             $pdimi defaults to zero. If $pdimi is zero, then it should hold that:
3219              
3220             all( $pi2nzi==sequence($ccs->nstored_p) )
3221              
3222             =item getptr($pdimi)
3223              
3224             Guts for ptr(). Does not check $PTRS and does not cache anything.
3225              
3226             =item clearptr($pdimi)
3227              
3228             Clears any cached Harwell-Boeing pointers for physically indexed dimension $pdimi.
3229              
3230             =item clearptrs()
3231              
3232             Clears any cached Harwell-Boeing pointers.
3233              
3234             =item flags()
3235              
3236             =item flags($flags)
3237              
3238             Get/set object-local $FLAGS.
3239              
3240              
3241             =item bad_is_missing()
3242              
3243             =item bad_is_missing($bool)
3244              
3245             Get/set the value of the object-local "bad-is-missing" flag.
3246             If this flag is set, BAD values in $VALS are considered "missing",
3247             regardless of the current value of $missing.
3248              
3249             =item badmissing()
3250              
3251             Sets the "bad-is-missing" flag and returns the calling object.
3252              
3253             =item nan_is_missing()
3254              
3255             =item nan_is_missing($bool)
3256              
3257             Get/set the value of the object-local "NaN-is-missing" flag.
3258             If this flag is set, NaN (and +inf, -inf) values in $VALS are considered "missing",
3259             regardless of the current value of $missing.
3260              
3261             =item nanmissing()
3262              
3263             Sets the "nan-is-missing" flag and returns the calling object.
3264              
3265             =back
3266              
3267             =cut
3268              
3269             ##======================================================================
3270             ## Methods: General Information
3271             ##======================================================================
3272             =pod
3273              
3274             =head2 General Information
3275              
3276             =over 4
3277              
3278             =item density()
3279              
3280             Returns the number of non-missing values divided by the number
3281             of indexable values in the sparse object as a perl scalar.
3282              
3283             =item compressionRate()
3284              
3285             Returns the compression rate of the PDL::CCS::Nd object
3286             compared to a dense piddle of the physically indexable dimensions.
3287             Higher values indicate better compression (e.g. lower density).
3288             Negative values indicate that dense storage would be more
3289             memory-efficient. Pointers are not included in the computation
3290             of the compression rate.
3291              
3292             =back
3293              
3294             =cut
3295              
3296             ##======================================================================
3297             ## Footer Administrivia
3298             ##======================================================================
3299              
3300             ##---------------------------------------------------------------------
3301             =pod
3302              
3303             =head1 ACKNOWLEDGEMENTS
3304              
3305             Perl by Larry Wall.
3306              
3307             PDL by Karl Glazebrook, Tuomas J. Lukka, Christian Soeller, and others.
3308              
3309             =cut
3310              
3311             ##----------------------------------------------------------------------
3312             =pod
3313              
3314             =head1 KNOWN BUGS
3315              
3316             Many.
3317              
3318             =cut
3319              
3320              
3321             ##---------------------------------------------------------------------
3322             =pod
3323              
3324             =head1 AUTHOR
3325              
3326             Bryan Jurish Emoocow@cpan.orgE
3327              
3328             =head2 Copyright Policy
3329              
3330             Copyright (C) 2007-2018, Bryan Jurish. All rights reserved.
3331              
3332             This package is free software, and entirely without warranty.
3333             You may redistribute it and/or modify it under the same terms
3334             as Perl itself.
3335              
3336             =head1 SEE ALSO
3337              
3338             perl(1),
3339             PDL(3perl),
3340             PDL::SVDLIBC(3perl),
3341             PDL::CCS::Nd(3perl),
3342              
3343             SVDLIBC: http://tedlab.mit.edu/~dr/SVDLIBC/
3344              
3345             SVDPACKC: http://www.netlib.org/svdpack/
3346              
3347             =cut