## File: PDL::CCS::Compat.pm ## Author: Bryan Jurish ## Description: backwards-compatibility hacks for PDL::CCS package PDL::CCS::Compat; use PDL; use PDL::VectorValued; use PDL::CCS::Config qw(ccs_indx); use PDL::CCS::Functions; use PDL::CCS::Utils; use PDL::CCS::Ufunc; use PDL::CCS::Ops; use strict; our $VERSION = '1.24.1'; ##-- update with perl-reversion from Perl::Version module our @ISA = ('PDL::Exporter'); our @ccs_binops = (qw(plus minus mult divide modulo power), qw(gt ge lt le eq ne spaceship), qw(and2 or2 xor shiftleft shiftright), ); our @EXPORT_OK = ( ## ##-- Encoding qw(ccs_encode_compat), qw(ccsencode ccsencode_nz ccsencodefull ccsencodefull_nz), qw(ccsencodea ccsencode_naz ccsencodefulla ccsencodefull_naz), qw(ccsencodeg ccsencode_g ccsencodefullg ccsencodefull_g), qw(ccsencodei ccsencode_i ccsencodefulli ccsencodefull_i), qw(ccsencodei2d ccsencode_i2d ccsencodefulli2d ccsencodefull_i2d), ## ##-- Decoding qw(_ccsdecodecols ccsdecodecols), qw(ccsdecode ccsdecodefull), qw(ccsdecode_g ccsdecodeg ccsdecodefull_g ccsdecodefullg), ## ##-- Indexing qw(ccsiNDtonzi ccsi2dtonzi ccsitonzi), qw(ccswhichND ccswhich2d ccswhichfull ccswhich), qw(ccstranspose ccstransposefull), ## ##-- Lookup qw(ccsget ccsget2d), ## ##-- Operations (map {("ccs${_}_cv","ccs${_}_rv")} (@ccs_binops,qw(add diff))), ## ##-- Ufuncs (map {("ccs${_}","ccs${_}t")} qw(sumover prodover)), ); our %EXPORT_TAGS = ( Func => [@EXPORT_OK], ##-- respect PDL conventions (hopefully) ); ##====================================================================== ## pod: headers =pod =head1 NAME PDL::CCS::Compat - Backwards-compatibility module for PDL::CCS =head1 SYNOPSIS use PDL; use PDL::CCS::Compat; ##-- source pdl $a = random($N=8,$M=7); ##--------------------------------------------------------------------- ## Non-missing value counts $nnz = $a->flat->nnz; ##-- "missing" == 0 $nnaz = $a->flat->nnza(1e-6); ##-- "missing" ~= 0 #$ngood = $a->ngood; ##-- "missing" == BAD (see PDL::Bad) ##--------------------------------------------------------------------- ## CCS Encoding ($ptr,$rowids,$vals) = ccsencode_nz ($a); # missing == 0 ($ptr,$rowids,$vals) = ccsencode_naz($a,$eps); # missing ~= 0 ($ptr,$rowids,$vals) = ccsencode_g ($a); # missing == BAD ($ptr,$rowids,$vals) = ccsencode_i ($i,$ivals,$N); # generic flat ($ptr,$rowids,$vals) = ccsencode_i2d($xi,$yi,$ivals); # generic 2d ##--------------------------------------------------------------------- ## CCS Decoding $cols = ccsdecodecols($ptr,$rowids,$nzvals, $xvals $a2 = ccsdecode ($ptr,$rowids,$vals); # missing == 0 $a2 = ccsdecode_g($ptr,$rowids,$vals); # missing == BAD ##--------------------------------------------------------------------- ## CCS Index Conversion $nzi = ccsitonzi ($ptr,$rowids, $ix, $missing); # ix => nzi $nzi = ccsi2dtonzi($ptr,$rowids, $xi,$yi, $missing); # 2d => nzi $ix = ccswhich ($ptr,$rowids,$vals); # CCS => ix ($xi,$yi) = ccswhichND($ptr,$rowids,$vals); # CCS => 2d $xyi = ccswhichND($ptr,$rowids,$vals); # ...as scalar ##--------------------------------------------------------------------- ## CCS Lookup $ixvals = ccsget ($ptr,$rowids,$vals, $ix,$missing); # ix => values $ixvals = ccsget2d($ptr,$rowids,$vals, $xi,$yi,$missing); # 2d => values ##--------------------------------------------------------------------- ## CCS Operations ($ptrT,$rowidsT,$valsT) = ccstranspose($ptr,$rowids,$vals); # CCS<->CRS ##--------------------------------------------------------------------- ## Vector Operations, by column $nzvals_out = ccsadd_cv ($ptr,$rowids,$nzvals, $colvec); $nzvals_out = ccsdiff_cv($ptr,$rowids,$nzvals, $colvec); $nzvals_out = ccsmult_cv($ptr,$rowids,$nzvals, $colvec); $nzvals_out = ccsdiv_cv ($ptr,$rowids,$nzvals, $colvec); ##--------------------------------------------------------------------- ## Vector Operations, by row $nzvals_out = ccsadd_rv ($ptr,$rowids,$nzvals, $rowvec); $nzvals_out = ccsdiff_rv($ptr,$rowids,$nzvals, $rowvec); $nzvals_out = ccsmult_rv($ptr,$rowids,$nzvals, $rowvec); $nzvals_out = ccsdiv_rv ($ptr,$rowids,$nzvals, $rowvec); ##--------------------------------------------------------------------- ## Scalar Operations $nzvals_out = $nzvals * 42; # ... or whatever ##--------------------------------------------------------------------- ## Accumulators $rowsumover = ccssumover ($ptr,$rowids,$nzvals); ##-- like $a->sumover() $colsumovert = ccssumovert($ptr,$rowids,$nzvals); ##-- like $a->xchg(0,1)->sumover =cut ##====================================================================== ## Encoding =pod =head1 Encoding =cut ##--------------------------------------------------------------- ## Encoding: generic =pod =head2 ccs_encode_compat =for sig Signature: (indx awhich(2,Nnz); avals(Nnz); indx $N; indx $M; indx [o]ptr(N); indx [o]rowids(Nnz); [o]nzvals(Nnz)) Generic wrapper for backwards-compatible ccsencode() variants. =cut *ccs_encode_compat = \&PDL::ccs_encode_compat; sub PDL::ccs_encode_compat { my ($aw,$avals,$N,$M,$ptr,$rowids,$nzvals) = @_; $N = $aw->slice("(0),")->max+1 if (!defined($N)); $M = $aw->slice("(1),")->max+1 if (!defined($M)); my ($ptr1,$awi) = ccs_encode_pointers($aw->slice("(0),"), $N); if (defined($ptr)) { $ptr .= $ptr1->slice("0:-2"); } else { $ptr = $ptr1->slice("0:-2"); $ptr->sever; } if (defined($rowids)) { $rowids .= $aw->slice("(1),")->index($awi); } else { $rowids = $aw->slice("(1),")->index($awi); $rowids->sever; } if (defined($nzvals)) { $nzvals .= $avals->index($awi); } else { $nzvals = $avals->index($awi); $nzvals->sever; } return ($ptr,$rowids,$nzvals); } ##--------------------------------------------------------------- ## Encoding: MISSING=zero =pod =head2 ccsencode =head2 ccsencode_nz =for sig Signature: (a(N,M); indx [o]ptr(N); indx [o]rowids(Nnz); [o]nzvals(Nnz)) Encodes matrix $a() in compressed column format, interpreting zeroes as "missing" values. Allocates output vectors if required. =cut *ccsencode = *ccsencodefull = *ccsencodefull_nz = *PDL::ccsencode = *PDL::ccsencode_nz = *PDL::ccsencodefull = *PDL::ccsencodefull_nz = \&ccsencode_nz; sub ccsencode_nz { #my ($a,$ptr,$rowids,$nzvals) = @_; my $a = shift; $a = $a->clump(1+$a->ndims-2); ##-- clump(-2) broken in PDL-2.0.14 my $aw = $a->whichND; return ccs_encode_compat($aw, $a->indexND($aw), $a->dims, @_); } ##--------------------------------------------------------------- ## Encoding: MISSING=ZERO (approx) =pod =head2 ccsencodea =head2 ccsencode_naz =for sig Signature: (a(N,M); eps(); indx [o]ptr(N); indx [o]rowids(Nnz); [o]nzvals(Nnz)) Encodes matrix $a() in CCS format interpreting approximate zeroes as "missing" values. This function is just like ccsencode_nz(), but uses the tolerance parameter $eps() to determine which elements are to be treated as zeroes. Allocates output vectors if required. =cut *ccsencodea = *ccsencodefulla = *ccsencodefull_naz = *PDL::ccsencodea = *PDL::ccsencode_naz = *PDL::ccsencodefulla = *PDL::ccsencodefull_naz = \&ccsencode_naz; sub ccsencode_naz { #my ($a,$eps,$ptr,$rowids,$nzvals) = @_; my $a = shift; my $eps = shift; $a = $a->clump(1+$a->ndims-2); ##-- clump(-2) is broken in PDL-2.014 my $aw = $a->approx(0,$eps)->inplace->not->whichND; ##-- FIXME: optimize return ccs_encode_compat($aw, $a->indexND($aw), $a->dims, @_); } ##--------------------------------------------------------------- ## Encoding: MISSING=BAD =pod =head2 ccsencodeg =head2 ccsencode_g =for sig Signature: (a(N,M); indx [o]ptr(N); indx [o]rowids(Nnz); [o]nzvals(Nnz)) Encodes matrix $a() in CCS format interpreting BAD values as "missing". Requires bad-value support built into PDL. Allocates output vectors if required. =cut *ccsencodeg = *ccsencodefullg = *ccsencodefull_g = *PDL::ccsencodeg = *PDL::ccsencode_g = *PDL::ccsencodefullg = *PDL::ccsencodefull_g = \&ccsencode_g; sub ccsencode_g { #my ($a,$ptr,$rowids,$nzvals) = @_; my $a = shift; $a = $a->clump(1+$a->ndims-2); ##-- clump(-2) is broken in PDL-v2.014 my $amask = zeroes(byte,$a->dims); $a->isgood($amask); my $aw = $amask->whichND; return ccs_encode_compat($aw, $a->indexND($aw), $a->dims, @_); } ##--------------------------------------------------------------- ## Encoding: from flat index =pod =head2 ccsencode_i =for sig Signature: (indx ix(Nnz); nzvals(Nnz); indx $N; int [o]ptr(N); indx [o]rowids(Nnz); [o]nzvals_enc(Nnz)) General-purpose CCS encoding method for flat indices. Encodes values $nzvals() from flat-index locations $ix() into a CCS matrix ($ptr(), $rowids(), $nzvals_enc()). Allocates output vectors if required. $N (~ $a-Edim(0)) must be specified. =cut *ccsencodei = *ccsencodefulli = *ccsencodefull_i = *PDL::ccsencodei = *PDL::ccsencode_i = *PDL::ccsencodefulli = *PDL::ccsencodefull_i = \&ccsencode_i; sub ccsencode_i { #my ($iflat,$avals,$N_optional,$ptr,$rowids,$nzvals) = @_; my ($iflat,$avals) = splice(@_,0,2); my $N = defined($_[0]) && (!ref($_[0]) || $_[0]->nelem==1) ? shift : $_[0]->nelem; my $aw = ($iflat % $N)->cat($iflat/$N)->xchg(0,1); return ccs_encode_compat($aw, $avals, $N, undef, @_); } ##--------------------------------------------------------------- ## Encoding: from 2d index =pod =head2 ccsencode_i2d =for sig Signature: ( indx xvals(Nnz) ; indx yvals(Nnz) ; nzvals(Nnz) ; indx $N ; ##-- optional indx [o]ptr(N) ; indx [o]rowids(Nnz) ; [o]nzvals_enc(Nnz); ) General-purpose encoding method. Encodes values $nzvals() from 2d-index locations ($xvals(), $yvals()) in an $N-by-(whatever) PDL into a CCS matrix $ptr(), $rowids(), $nzvals_enc(). Allocates output vectors if required. If $N is omitted, it defaults to the maximum column index given in $xvals(). =cut *ccsencodei2d = *ccsencodefulli2d = *ccsencodefull_i2d = *PDL::ccsencodei2d = *PDL::ccsencode_i2d = *PDL::ccsencodefulli2d = *PDL::ccsencodefull_i2d = \&ccsencode_i2d; sub ccsencode_i2d { #my ($whichx,$whichy,$avals,$N_optional,$ptr,$rowids,$nzvals) = @_; my ($whichx,$whichy,$avals) = splice(@_, 0, 3); my $aw = $whichx->cat($whichy)->xchg(0,1); my $N = defined($_[0]) && (!ref($_[0]) || $_[0]->nelem==1) ? shift : ($whichx->max+1); return ccs_encode_compat($aw, $avals, $N, undef, @_); } ##====================================================================== ## Decoding =pod =head1 Decoding =cut ##--------------------------------------------------------------- ## Decoding: column-wise =pod =head2 ccsdecodecols =for sig Signature: ( indx ptr (N) ; indx rowids (Nnz); nzvals (Nnz); indx xvals (I) ; # default=sequence($N) missing() ; # default=0 M () ; # default=rowids->max+1 [o]cols (I,M); # default=new ) Extract dense columns from a CCS-encoded matrix (no dataflow). Allocates output matrix if required. If $a(N,M) was the dense source matrix for the CCS-encoding, and if missing values are zeros, then the following two calls are equivalent (modulo data flow): $cols = $a->dice_axis(1,$col_ix); $cols = ccsdecodecols($ptr,$rowids,$nzvals, $col_ix,0); =cut *PDL::_ccsdecodecols = \&_ccsdecodecols; #Pars => 'indx ptr(N); indx rowids(Nnz); nzvals(Nnz); indx col_ix(I); missing(); [o]cols(I,M);', sub _ccsdecodecols { ccsdecodecols(@_[0,1,2], $_[3],$_[4], undef, $_[5]); } *PDL::ccsdecodecols = \&ccsdecodecols; sub ccsdecodecols { my ($ptr,$rowids,$nzvals, $coli,$missing,$M, $cols) = @_; $coli = sequence(ccs_indx,$ptr->dim(0)) if (!defined($coli)); $coli = pdl(ccs_indx,$coli) if (!ref($coli)); my $ptr1 = zeroes(ccs_indx,$ptr->nelem+1); $ptr1->slice("0:-2") .= $ptr; $ptr1->set(-1 => $nzvals->nelem); $M = $rowids->max+1 if (!defined($M)); my ($ptrix,$nzix) = ccs_decode_pointer($ptr1,$coli); my $which = $ptrix->cat($rowids->index($nzix))->xchg(0,1); if (!defined($cols)) { $cols = ccs_decode($which, $nzvals->index($nzix), $missing, [$coli->nelem,$M]); $cols->sever; ##-- compat } else { ccs_decode($which, $nzvals->index($nzix), $missing, [$coli->nelem,$M], $cols); } return $cols; } ##--------------------------------------------------------------- ## Decoding: MISSING=0 =pod =head2 ccsdecode =for sig Signature: (indx ptr(N); indx rowids(Nnz); nzvals(Nnz); $M; [o]dense(N,M)) Decodes compressed column format vectors $ptr(), $rowids(), and $nzvals() into dense output matrix $a(). Allocates the output matrix if required. Note that if the original matrix (pre-encoding) contained trailing rows with no nonzero elements, such rows will not be allocated by this method (unless you specify either $M or $dense). In such cases, you might prefer to call ccsdecodecols() directly. =cut *PDL::ccsdecodefull = \&ccsdecodefull; ##-- (indx ptr(N); indx rowids(Nnz); nzvals(Nnz); [o]dense(N,M)) sub ccsdecodefull { ccsdecodecols(@_[0,1,2], undef,0,undef, @_[3..$#_]); } *PDL::ccsdecode = \&ccsdecode; sub ccsdecode { my ($ptr,$rowids,$nzvals, $M, $dense)=@_; if (!defined($dense)) { ##-- check for old calling convention (is $M a multi-dim PDL?) if (ref($M) && UNIVERSAL::isa($M, 'PDL') && $M->dim(0)==$ptr->dim(0)) { $dense = $M; } else { $M = $rowids->max+1 if (!defined($M)); $dense = zeroes($nzvals->type,$ptr->dim(0),$M); } } ccsdecodecols($ptr,$rowids,$nzvals, undef,0,$M, $dense); return $dense; } ##--------------------------------------------------------------- ## Decoding: MISSING=BAD =pod =head2 ccsdecode_g =for sig Signature: (indx ptr(N); indx rowids(Nnz); nzvals(Nnz); $M; [o]dense(N,M)) Convenience method. Like ccsdecode() but sets "missing" values to BAD. =cut *ccsdecodefullg = *PDL::ccsdecodefullg = *PDL::ccsdecodefull_g = \&ccsdecodefull_g; sub ccsdecodefull_g { my $badval = pdl($_[2]->type,0)->setvaltobad(0); ccsdecodecols(@_[0,1,2], undef,$badval,undef,undef, @_[3..$#_]); } *ccsdecodeg = *PDL::ccsdecodeg = *PDL::ccsdecode_g = \&ccsdecode_g; sub ccsdecode_g { my ($ptr,$rowids,$nzvals, $M, $dense)=@_; if (!defined($dense)) { ##-- check for old calling convention (is $M a multi-dim PDL?) if (ref($M) && UNIVERSAL::isa($M, 'PDL') && $M->dim(0)==$ptr->dim(0)) { $dense = $M; } else { $M = $rowids->max+1 if (!defined($M)); $dense = zeroes($nzvals->type,$ptr->dim(0),$M); } } my $badval = pdl($nzvals->type,0)->setvaltobad(0); ccsdecodecols($ptr,$rowids,$nzvals, undef,$badval,$M, $dense); return $dense; } ##====================================================================== ## Index Conversion ##====================================================================== =pod =head1 Index Conversion =cut ##------------------------------------------------------ ## ccsiNDtonzi() : index conversion: N-dimensional =pod =for sig Signature: (indx ptr(N); indx rowids(Nnz); indx ind(2,I); indx missing(); indx [o]nzix(I)) =head2 ccsiNDtonzi Convert N-dimensional index values $ind() appropriate for a dense matrix (N,M) into indices $nzix() appropriate for the $rowids() and/or $nzvals() components of the CCS-encoded matrix ($ptr(),$rowids(),$nzvals()). Missing values are returned in $nzix() as $missing(). =cut *PDL::ccsiNDtonzi = \&ccsiNDtonzi; sub ccsiNDtonzi { my ($ptr,$rowids,$ind, $missing, $nzix) = @_; my ($ptri,$ptrnzi) = ccs_decode_pointer($ptr->append($rowids->nelem)); my $ccswnd = $ptri->cat($rowids->index($ptrnzi))->xchg(0,1)->vv_qsortvec; $nzix = $ind->vsearchvec($ccswnd); my $nzix_mask = ($ind==$ccswnd->dice_axis(1,$nzix))->andover; $nzix_mask->inplace->not; #(my $tmp = $nzix->where($nzix_mask)) .= $missing; ##-- fix "Can't modify non-lvalue subroutine call" in 5.15.x (perl bug #107366) $nzix->where($nzix_mask) .= $missing; return $nzix; } ##------------------------------------------------------ ## ccsi2dtonzi() : index conversion: 2d =pod =head2 ccsi2dtonzi =for sig Signaure: (indx ptr(N); indx rowids(Nnz); indx col_ix(I); indx row_ix(I); indx missing(); indx [o]nzix(I)) Convert 2d index values $col_ix() and $row_ix() appropriate for a dense matrix (N,M) into indices $nzix() appropriate for the $rowids() and/or $nzvals() components of the CCS-encoded matrix ($ptr(),$rowids(),$nzvals()). Missing values are returned in $nzix() as $missing(). =cut *PDL::ccsi2dtonzi = \&ccsi2dtonzi; sub ccsi2dtonzi { my ($ptr,$rowids,$xi,$yi, $missing, $nzix) = @_; return ccsiNDtonzi($ptr,$rowids, $xi->cat($yi)->xchg(0,1), $missing,$nzix); } ##------------------------------------------------------ ## ccsitonzi() : index conversion: flat =pod =for sig Signature: (indx ptr(N); indx rowids(Nnz); indx ix(I); indx missing(); indx [o]nzix(I)) =head2 ccsitonzi Convert flat index values $ix() appropriate for a dense matrix (N,M) into indices $nzix() appropriate for the $rowids() and/or $nzvals() components of the CCS-encoded matrix ($ptr(),$rowids(),$nzvals()). Missing values are returned in $nzix() as $missing(). =cut *PDL::ccsitonzi = \&ccsitonzi; sub ccsitonzi { my ($ptr,$rowids,$ix, $missing, $nzix) = @_; my $dummy = pdl(byte,0)->slice("*".($ptr->dim(0)).",*".($rowids->max+1)); my ($xi,$yi) = $dummy->one2nd($ix); return ccsiNDtonzi($ptr,$rowids, $xi->cat($yi)->xchg(0,1), $missing,$nzix); } ##------------------------------------------------------ ## ccswhichND: get indices (N-dimensional) =pod =head2 ccswhichND =head2 ccswhich2d =head2 ccswhichfull =for sig Signature: (indx ptr(N); indx rowids(Nnz); nzvals(Nnz); indx [o]which_cols(Nnz); indx [o]which_rows(Nnz)', In scalar context, returns concatenation of $which_cols() and $which_rows(), similar to the builtin whichND(). Note however that ccswhichND() may return its index PDLs sorted in a different order than the builtin whichND() method for dense matrices. Use the qsort() or qsorti() methods if you need sorted index PDLs. =cut *ccswhich2d = *PDL::which2d = *PDL::ccswhichND = *ccswhichfull = *PDL::ccswhichfull = \&ccswhichND; sub ccswhichND { my ($ptr,$rowids,$nzvals, $which_cols,$which_rows) = @_; my ($ptrnzi); ($which_cols,$ptrnzi) = ccs_decode_pointer($ptr->append($rowids->nelem), sequence(ccs_indx, $ptr->nelem), $which_cols ); $which_rows = zeroes(ccs_indx, $rowids->nelem) if (!defined($which_rows)); $which_rows .= $rowids->index($ptrnzi); return wantarray ? ($which_cols,$which_rows) : $which_cols->cat($which_rows)->xchg(0,1); } ##------------------------------------------------------ ## ccswhich(): get indices (flat) =head2 ccswhich =for sig Signature: (indx ptr(N); indx rowids(Nnz); nzvals(Nnz); indx [o]which(Nnz); indx [t]wcols(Nnz)', Convenience method. Calls ccswhichfull(), and scales the output PDLs to correspond to a flat enumeration. The output PDL $which() is B guaranteed to be sorted in any meaningful order. Use the qsort() method if you need sorted output. =cut *PDL::ccswhich = \&ccswhich; sub ccswhich { my ($ptr,$rowids,$nzvals, $which, $wcols) = @_; my $nnz = $rowids->dim(0); $which = zeroes(ccs_indx,$nnz) if (!defined($which)); $wcols = zeroes(ccs_indx,$nnz) if (!defined($wcols)); ccswhichfull($ptr,$rowids,$nzvals, $wcols,$which); $which *= $ptr->dim(0); $which += $wcols; return $which; } ##------------------------------------------------------ ## ccstranspose() : transposition (convenience) =pod =head2 ccstranspose =for sig Signature: (indx ptr(N); indx rowids(Nnz); nzvals(Nnz); indx [o]ptrT(M); indx [o]rowidsT(Nnz); [o]nzvalsT(Nnz)', Transpose a compressed matrix. =cut *ccstransposefull = *PDL::ccstransposefull = *PDL::ccstranspose = \&ccstranspose; sub ccstranspose { my ($ptr,$rowids,$nzvals, $ptrT,$rowidsT,$nzvalsT)=@_; my $N = $ptr->dim(0); my $M = defined($ptrT) ? $ptrT->dim(0) : $rowids->max+1; my $wnd = ccswhichND($ptr,$rowids,$nzvals)->slice("1:0,"); return ccs_encode_compat($wnd,$nzvals,$M,$N, $ptrT,$rowidsT,$nzvalsT); } ##====================================================================== ## Lookup ##====================================================================== =pod =head1 Lookup =cut ##------------------------------------------------------ ## ccsget2d() : lookup: 2d =pod =head2 ccsget2d =for sig Signature: (indx ptr(N); indx rowids(Nnz); nzvals(Nnz); indx xvals(I); indx yvals(I); missing(); [o]ixvals(I)) Lookup values in a CCS-encoded PDL by 2d source index (no dataflow). Pretty much like ccsi2dtonzi(), but returns values instead of indices. If you know that your index PDLs $xvals() and $yvals() do not refer to any missing values in the CCS-encoded matrix, then the following two calls are equivalent (modulo dataflow): $ixvals = ccsget2d ($ptr,$rowids,$nzvals, $xvals,$yvals,0); $ixvals = index($nzvals, ccsi2dtonzi($ptr,$rowids, $xvals,$yvals,0)); The difference is that only the second incantation will cause subsequent changes to $ixvals to be propagated back into $nzvals. =cut *PDL::ccsget2d = \&ccsget2d; sub ccsget2d { my ($ptr,$rowids,$nzvals, $xi,$yi, $missing, $ixnzvals) = @_; my $nzi = ccsi2dtonzi($ptr,$rowids, $xi,$yi, -1); my $nzi_isgood = ($nzi != -1); $ixnzvals = zeroes($nzvals->type, $xi->nelem) if (!defined($ixnzvals)); if (!all($nzi_isgood)) { my $tmp; ($tmp=$ixnzvals->where( $nzi_isgood)) .= $nzvals->index($nzi->where($nzi_isgood)); ($tmp=$ixnzvals->where(!$nzi_isgood)) .= $missing; $ixnzvals->badflag(1) if (PDL->topdl($missing)->badflag); } return $ixnzvals; } ##------------------------------------------------------ ## ccsget() : lookup: flat =pod =head2 ccsget =for sig Signature: (indx ptr(N); indx rowids(Nnz); nzvals(Nnz); indx ix(I); missing(); [o]ixvals(I)) Lookup values in a CCS-encoded PDL by flat source index (no dataflow). Pretty much like ccsitonzi(), but returns values instead of indices. If you know that your index PDL $ix() does not refer to any missing values in the CCS-encoded matrix, then the following two calls are equivalent (modulo dataflow): $ixvals = ccsget ($ptr,$rowids,$nzvals, $ix,0); $ixvals = index($nzvals, ccsitonzi($ptr,$rowids, $ix,0)) The difference is that only the second incantation will cause subsequent changes to $ixvals to be propagated back into $nzvals. =cut *PDL::ccsget = \&ccsget; sub ccsget { my ($ptr,$rowids,$nzvals, $ix, $missing, $ixnzvals) = @_; my $nzi = ccsitonzi($ptr,$rowids, $ix,-1); my $nzi_isgood = ($nzi != -1); $ixnzvals = zeroes($nzvals->type, $ix->nelem) if (!defined($ixnzvals)); if (!all($nzi_isgood)) { my $tmp; ($tmp=$ixnzvals->where( $nzi_isgood)) .= $nzvals->index($nzi->where($nzi_isgood)); ($tmp=$ixnzvals->where(!$nzi_isgood)) .= $missing; $ixnzvals->badflag(1) if (PDL->topdl($missing)->badflag); } return $ixnzvals; } ##====================================================================== ## Vector Operations ##====================================================================== =pod =head1 Vector Operations =cut ##====================================================================== ## Vector Operations: Columns ##====================================================================== =pod =head2 ccs${OP}_cv =for sig Signature: (indx ptr(N); indx rowids(Nnz); nzvals_in(Nnz); colvec(M); [o]nzvals_out(Nnz)) Column-vector operation ${OP} on CCS-encoded PDL. Should do something like the following (without decoding the CCS matrix): ($colvec ${OP} ccsdecode(\$ptr,\$rowids,\$nzvals))->ccsencode; Missing values in the CCS-encoded PDL are not affected by this operation. ${OP} is one of the following: plus ##-- addition (alias: 'add') minus ##-- subtraction (alias: 'diff') mult ##-- multiplication (NOT matrix-multiplication) divide ##-- division (alias: 'div') modulo ##-- modulo power ##-- potentiation gt ##-- greater-than ge ##-- greater-than-or-equal lt ##-- less-than le ##-- less-than-or-equal eq ##-- equality ne ##-- inequality spaceship ##-- 3-way comparison and2 ##-- binary AND or2 ##-- binary OR xor ##-- binary XOR shiftleft ##-- left-shift shiftright ##-- right-shift =cut sub ccs_binop_compat_cv { my $ccsop = shift; return sub { $ccsop->(@_[1,2,3,4]) }; } foreach my $op (@ccs_binops) { no strict 'refs'; *{"ccs${op}_cv"} = *{"PDL::ccs${op}_cv"} = ccs_binop_compat_cv(\&{"PDL::ccs_${op}_vector_mia"}); } *ccsadd_cv = *PDL::ccsadd_cv = \&ccsplus_cv; *ccsdiff_cv = *PDL::ccsdiff_cv = \&ccsminus_cv; *ccsdiv_cv = *PDL::ccsdiv_cv = \&ccsdivide_cv; ##====================================================================== ## Vector Operations: Rows ##====================================================================== =pod =head2 ccs${OP}_rv =for sig Signature: (indx ptr(N); indx rowids(Nnz); nzvals_in(Nnz); rowvec(N); [o]nzvals_out(Nnz)) Row-vector operation ${OP} on CCS-encoded PDL. Should do something like the following (without decoding the CCS matrix): ($column->slice("*1,") ${OP} ccsdecode($ptr,$rowids,$nzvals))->ccsencode; Missing values in the CCS-encoded PDL are not effected by this operation. See ccs${OP}_cv() above for supported operations. =cut sub ccs_binop_compat_rv { my $ccsop = shift; return sub { my $ptr = shift; my ($ptri,$ptrnzi) = ccs_decode_pointer($ptr->append($_[1]->nelem)); $ccsop->($ptri, $_[1]->index($ptrnzi), @_[2,3]); }; } foreach my $op (@ccs_binops) { no strict 'refs'; *{"ccs${op}_rv"} = *{"PDL::ccs${op}_rv"} = ccs_binop_compat_rv(\&{"PDL::ccs_${op}_vector_mia"}); } *ccsadd_rv = *PDL::ccsadd_rv = \&ccsplus_rv; *ccsdiff_rv = *PDL::ccsdiff_rv = \&ccsminus_rv; *ccsdiv_rv = *PDL::ccsdiv_rv = \&ccsdivide_rv; ##------------------------------------------------------ ## Ufuncs (accumulators) ## \&ufuncsub = ccs_ufunc_compat(\&ccs_accum_sub) sub ccs_ufunc_compat { my $sub = shift; return sub { my ($ptr,$rowids,$nzvals, $M,$rowvals) = @_; my ($ixout,$valsout) = $sub->($rowids->slice("*1,"),$nzvals, 0,0); $M = $rowids->max+1 if (!defined($M)); $rowvals = zeroes($nzvals->type,$M) if (!defined($rowvals)); $rowvals->index($ixout->flat) .= $valsout; return $rowvals; }; } ## \&ufuncsub = ccs_ufunc_compat_t(\&ccs_accum_sub) sub ccs_ufunc_compat_t { my $sub = shift; return sub { my ($ptr,$rowids,$nzvals, $colvals) = @_; my ($colids,$nzix) = ccs_decode_pointer($ptr->append($nzvals->nelem)); ccs_ufunc_compat(undef,$colids,$nzvals->index($nzix), $ptr->dim(0),$colvals); }; } *ccssumover = *PDL::ccssumover = ccs_ufunc_compat (\&ccs_accum_sum); *ccssumovert = *PDL::ccssumovert = ccs_ufunc_compat_t(\&ccs_accum_sum); *ccprodover = *PDL::ccsprodover = ccs_ufunc_compat (\&ccs_accum_prod); *ccsprodovert = *PDL::ccsprodovert = ccs_ufunc_compat_t(\&ccs_accum_prod); 1; ##-- make perl happy ##====================================================================== ## Footer Administrivia ##====================================================================== ##--------------------------------------------------------------------- =pod =head1 EXAMPLES =head2 Compressed Column Format Example $a = pdl([ [10, 0, 0, 0,-2, 0], [3, 9, 0, 0, 0, 3], [0, 7, 8, 7, 0, 0], [3, 0, 8, 7, 5, 0], [0, 8, 0, 9, 9, 13], [0, 4, 0, 0, 2, -1] ]); ($ptr,$rowids,$nzvals) = ccsencode($a); print join("\n", "ptr=$ptr", "rowids=$rowids", "nzvals=$nzvals"); ... prints something like: ptr=[0 3 7 9 12 16] rowids=[ 0 1 3 1 2 4 5 2 3 2 3 4 0 3 4 5 1 4 5] nzvals=[10 3 3 9 7 8 4 8 8 7 7 9 -2 5 9 2 3 13 -1] =head2 Sparse Matrix Example ##-- create a random sparse matrix $a = random(100,100); $a *= ($a>.9); ##-- encode it ($ptr,$rowids,$nzvals) = ccsencode($a); ##-- what did we save? sub pdlsize { return PDL::howbig($_[0]->type)*$_[0]->nelem; } print "Encoding saves us ", ($saved = pdlsize($a) - pdlsize($ptr) - pdlsize($rowids) - pdlsize($nzvals)), " bytes (", (100.0*$saved/pdlsize($a)), "%)\n"; ... prints something like: Encoding saves us 71416 bytes (89.27%) =head2 Decoding Example ##-- random matrix $a = random(100,100); ##-- make an expensive copy of $a by encoding & decoding ($ptr,$rowids,$nzvals) = ccsencode($a); $a2 = ccsdecode($ptr,$rowids,$nzvals); ##-- ...and make sure it's good print all($a==$a2) ? "Decoding is good!\n" : "Nasty icky bug!\n"; =cut ##--------------------------------------------------------------------- =pod =head1 ACKNOWLEDGEMENTS Perl by Larry Wall. PDL by Karl Glazebrook, Tuomas J. Lukka, Christian Soeller, and others. Original inspiration and algorithms from the SVDLIBC C library by Douglas Rohde; which is itself based on SVDPACKC by Michael Berry, Theresa Do, Gavin O'Brien, Vijay Krishna and Sowmini Varadhan. =cut ##---------------------------------------------------------------------- =pod =head1 KNOWN BUGS Many. =cut ##--------------------------------------------------------------------- =pod =head1 AUTHOR Bryan Jurish Emoocow@cpan.orgE =head2 Copyright Policy Copyright (C) 2005-2024, Bryan Jurish. All rights reserved. This package is free software, and entirely without warranty. You may redistribute it and/or modify it under the same terms as Perl itself. =head1 SEE ALSO perl(1), PDL(3perl), PDL::SVDLIBC(3perl), PDL::CCS::Nd(3perl), SVDLIBC: http://tedlab.mit.edu/~dr/SVDLIBC/ SVDPACKC: http://www.netlib.org/svdpack/ =cut