##-*- Mode: CPerl -*-
##======================================================================
## Header Administrivia
##======================================================================
use PDL::VectorValued::Dev;
my $VERSION = '1.24.1'; ##-- update with perl-reversion from Perl::Version module
pp_setversion($VERSION);
##------------------------------------------------------
## pm headers
pp_addpm({At=>'Top'},<<'EOPM');
#use PDL::CCS::Config;
use strict;
=pod
=head1 NAME
PDL::CCS::Utils - Low-level utilities for compressed storage sparse PDLs
=head1 SYNOPSIS
use PDL;
use PDL::CCS::Utils;
##---------------------------------------------------------------------
## ... stuff happens
=cut
EOPM
## /pm additions
##------------------------------------------------------
##------------------------------------------------------
## Exports: None
#pp_export_nothing();
##------------------------------------------------------
## Includes / defines
pp_addhdr(<<'EOH');
#include "ccsutils.h"
EOH
##------------------------------------------------------
## index datatype
require "../Config.pm";
our $INDX = $PDL::CCS::Config::ccsConfig{INDX_SIG};
pp_addpm( $PDL::CCS::Config::ccsConfig{INDX_FUNCDEF} );
pp_addhdr( $PDL::CCS::Config::ccsConfig{INDX_TYPEDEF} );
##======================================================================
## C Utilities
##======================================================================
# (none)
##======================================================================
## PDL::PP Wrappers
##======================================================================
##======================================================================
## Non-missing Value Counts
##======================================================================
pp_addpm(<<'EOPM');
=pod
=head1 Non-missing Value Counts
=cut
EOPM
##------------------------------------------------------
## nnz() : get number of nonzero values
pp_def(
'nnz',
Pars => "a(N); $INDX\ [o]nnz()",
HandleBad => 1,
Code => q{
broadcastloop %{
$nnz()=0;
loop (N) %{
if ($a()!=0) ++$nnz();
%}
%}
$PDLSTATESETGOOD(nnz); /*-- output is always good --*/
},
Doc => <<'EOD',
Get number of non-zero values in a PDL $a();
For 1d PDLs, should be equivalent to:
$nnz = nelem(which($a!=0));
For k>1 dimensional PDLs, projects via number of nonzero elements
to N-1 dimensions by computing the number of nonzero elements
along the the 1st dimension.
EOD
BadDoc=> 'The output PDL $nnz() never contains BAD values.',
);
##------------------------------------------------------
## nnza() : get number of non-approximate zero values
use PDL;
my %absfunc = (
map {
my $typ = PDL->can($_);
($typ
? ($typ->()->ppsym => ($typ->()->ctype eq 'long' ? "labs" : "llabs"))
: qw())
} qw (longlong indx)
);
pp_def(
'nnza',
Pars => "a(N); eps(); $INDX\ [o]nnz();",
HandleBad => 1,
Code => '
broadcastloop %{
$nnz()=0;
loop (N) %{
types(BSUL) %{ if ( abs($a()) > $eps()) ++$nnz(); %}
'.join("\n ", map {"types($_) %{ if ($absfunc{$_}(\$a()) > \$eps()) ++\$nnz(); %}"} sort keys(%absfunc)).'
types(F) %{ if (fabsf($a()) > $eps()) ++$nnz(); %}
types(D) %{ if (fabs ($a()) > $eps()) ++$nnz(); %}
%}
%}
$PDLSTATESETGOOD(nnz); /*-- output is always good --*/
',
Doc => <<'EOD',
Like nnz() using tolerance constant $eps().
For 1d PDLs, should be equivalent to:
$nnz = nelem(which(!$a->approx(0,$eps)));
EOD
BadDoc=> 'The output PDL $nnz() never contains BAD values.',
);
##======================================================================
## Encoding
##======================================================================
pp_addpm(<<'EOPM');
=pod
=head1 Encoding Utilities
=cut
EOPM
##------------------------------------------------------
## ccs_encode_pointers() : get encoded pointer & index translation PDL
pp_def(
'ccs_encode_pointers',
Pars => "$INDX\ ix(Nnz); $INDX\ N(); $INDX\ [o]ptr(Nplus1); $INDX\ [o]ixix(Nnz);",
PMCode=> q{
sub PDL::ccs_encode_pointers {
my ($ix,$N,$ptr,$ixix) = @_;
barf("Usage: ccs_encode_pointers(ix(Nnz), N(), [o]ptr(N+1), [o]ixix(Nnz)") if (!defined($ix));
&PDL::_ccs_encode_pointers_int($ix, ($N // $ix->max+1), ($ptr //= null), ($ixix //= null));
return ($ptr,$ixix);
}
},
RedoDimsCode => q{
if ($SIZE(Nplus1) < 0) {
$SIZE(Nplus1) = $N() + 1;
}
else if ($SIZE(Nplus1) <= $N()) {
$CROAK("dimension Nplus1 (=%td) must be greater than N (=%td)", $SIZE(Nplus1), $N());
}
},
Code => q{
/*-- Local variables --*/
CCS_Indx ixval, ixval_next, ixval_prev, nzi, nzj, sizeN=$SIZE(Nplus1)-1, sizeNnz=$SIZE(Nnz);
//
/*-- Count number of NZs in each column; store in ptr[N=>ixval] --*/
loop (Nplus1) %{ $ptr()=0; %}
loop (Nnz) %{ ixval=$ix(); ++$ptr(Nplus1=>ixval); %}
//
/*-- tweak ptr(): fill each cell with the starting point of the previous row --*/
ixval_prev = sizeN-1;
$ptr(Nplus1=>sizeN) = sizeNnz - $ptr(Nplus1=>ixval_prev);
for (ixval_next=sizeN, ixval=ixval_prev; ixval > 0; ixval_next=ixval--) {
ixval_prev = ixval-1;
$ptr(Nplus1=>ixval) = $ptr(Nplus1=>ixval_next) - $ptr(Nplus1=>ixval_prev);
}
$ptr(Nplus1=>0) = 0;
//
/*-- Assign columns and values --*/
for (nzi=0; nzi < sizeNnz; nzi++) {
ixval = $ix(Nnz=>nzi);
ixval_next = ixval+1;
nzj = $ptr(Nplus1=>ixval_next)++;
$ixix(Nnz=>nzj) = nzi;
}
},
Doc => <<'EOD'
General CCS encoding utility.
Get a compressed storage "pointer" vector $ptr
for a dimension of size $N with non-missing values at indices $ix. Also returns a vector
$ixix() which may be used as an index for $ix() to align its elements with $ptr()
along the compressed dimension.
The induced vector $ix-E<gt>index($ixix) is
guaranteed to be stably sorted along dimension $N():
\forall $i,$j with 1 <= $i < $j <= $Nnz :
$ix->index($ixix)->at($i) < $ix->index($ixix)->at($j) ##-- primary sort on $ix()
or
$ixix->at($i) < $ixix->at($j) ##-- ... stable
EOD
);
##======================================================================
## Decoding
##======================================================================
pp_addpm(<<'EOPM');
=pod
=head1 Decoding Utilities
=cut
EOPM
##------------------------------------------------------
## ccs_decode_pointer() : decode a CCS-encoded pointer
pp_def(
'ccs_decode_pointer',
Pars => "$INDX ptr(Nplus1); $INDX proj(Nproj); $INDX\ [o]projix(NnzProj); $INDX\ [o]nzix(NnzProj)",
OtherPars => 'PDL_Indx nnzProj;',
PMCode=> q{
sub PDL::ccs_decode_pointer {
my ($ptr,$proj,$projix,$nzix,$nnzproj) = @_;
barf("Usage: ccs_decode_pointer(ptr(N+1), proj(Nproj), [o]projix(NnzProj), [o]nzix(NnzProj), NnzProj?")
if (!defined($ptr));
if (!defined($proj)) {
$proj = PDL->sequence(ccs_indx(), $ptr->dim(0)-1);
$nnzproj //= $ptr->at(-1);
}
$projix //= null;
$nzix //= null;
$nnzproj //= ($projix->isnull && $nzix->isnull
? ($ptr->index($proj+1)-$ptr->index($proj))->sum
: -1);
return (null,null) if (!$nnzproj);
&PDL::_ccs_decode_pointer_int($ptr,$proj, $projix,$nzix, $nnzproj);
return ($projix,$nzix);
}
},
RedoDimsCode => q{
if ($SIZE(NnzProj) < 0) {
$SIZE(NnzProj) = $COMP(nnzProj);
}
},
Code => q{
/*-- Local variables --*/
CCS_Indx ni,ni_next, nzi,nzi_next, ixi=0, sizeNproj=$SIZE(Nproj), sizeNnzProj=$SIZE(NnzProj);
loop (Nproj) %{
ni = $proj();
ni_next = ni+1;
nzi = $ptr(Nplus1=>ni);
nzi_next = $ptr(Nplus1=>ni_next);
for ( ; nzi < nzi_next && ixi < sizeNnzProj; nzi++, ixi++) {
$projix(NnzProj=>ixi) = Nproj;
$nzix(NnzProj=>ixi) = nzi;
}
%}
},
Doc=><<'EOD'
General CCS decoding utility.
Project indices $proj() from a compressed storage "pointer" vector $ptr().
If unspecified, $proj() defaults to:
sequence($ptr->dim(0) - 1)
EOD
);
##======================================================================
## Indexing
##======================================================================
pp_addpm(<<'EOPM');
=pod
=head1 Indexing Utilities
=cut
EOPM
##------------------------------------------------------
## ccs_xindex1d()
## + optimized dice_axis on 0th dimension, no pointer required
vvpp_def(
'ccs_xindex1d',
Pars => "$INDX which(Ndims,Nnz); $INDX a(Na); $INDX\ [o]nzia(NnzA); $INDX\ [o]nnza()",
OtherPars => 'PDL_Indx sizeNnzA;',
PMCode=> q{
sub PDL::ccs_xindex1d {
my ($which,$a,$nzia,$nnza) = @_;
barf("Usage: ccs_xindex2d(which(Ndims,Nnz), a(Na), [o]nzia(NnzA), [o]nnza()")
if ((grep {!defined($_)} @_[0..1]) || $which->ndims < 2 || $which->dim(0) < 1);
$nzia //= null;
$nnza //= null;
&PDL::_ccs_xindex1d_int($which,$a,$nzia,$nnza, ($nnza ? $nnza->sclr : -1));
return wantarray ? ($nzia,$nnza) : $nzia->reshape($nnza->sclr);
}
},
RedoDimsCode => q(
if ($SIZE(NnzA) < 0) {
$SIZE(NnzA) = $COMP(sizeNnzA) >= 0 ? $COMP(sizeNnzA) : $SIZE(Nnz);
}
),
Code => q{
CCS_Indx a_min=0, a_max=$SIZE(Nnz);
CCS_Indx a_lb, a_ub, a_ubmax=a_max;
CCS_Indx nnzai = 0;
#if 0
/*-- DEBUG --*/
CCS_Indx size_nnz = $SIZE(Nnz);
CCS_Indx size_na = $SIZE(Na);
CCS_Indx size_nnza = $SIZE(NnzA);
printf("Nnz=%td, Na=%td [%td:%td], NnzA=%td\n", size_nnz, size_na,$a(Na=>0),$a(Na=>size_na-1), size_nnza);
#endif
loop (Na) %{
a_ubmax = a_max;
$LB('$a()', '$which(Ndims=>0,Nnz=>$_)', 'a_min','a_max', 'a_lb',ubmaxvar=>'a_ubmax');
if ($which(Ndims=>0,Nnz=>a_lb) != $a()) { a_min=a_lb; continue; }
$LB('$a()+1', '$which(Ndims=>0,Nnz=>$_)', 'a_lb' ,'a_ubmax', 'a_ub');
if ($which(Ndims=>0,Nnz=>a_ub) <= $a()) ++a_ub;
for ( ; a_lb < a_ub && nnzai < $SIZE(NnzA); ++a_lb, ++nnzai ) {
$nzia(NnzA=>nnzai) = a_lb;
}
if (nnzai >= $SIZE(NnzA)) break;
if (a_ub < a_max) a_min = a_ub;
%}
$nnza() = nnzai;
for ( ; nnzai < $SIZE(NnzA); ++nnzai) {
$nzia(NnzA=>nnzai) = -1;
}
},
Doc=><<'EOD'
Compute indices $nzai() along dimension C<Nnz> of $which() whose initial values $which(0,$nzai)
match some element of $a(). Appropriate for indexing a sparse encoded PDL
with non-missing entries at $which()
along the 0th dimension, a la L<dice_axis(0,$a)|PDL::Slices/dice_axis>.
$which((0),) and $a() must be both sorted in ascending order.
In list context, returns a list ($nzai,$nnza), where $nnza() is the number of indices found,
and $nzai are those C<Nnz> indices. In scalar context, trims the output vector $nzai() to $nnza()
elements.
EOD
);
##------------------------------------------------------
## ccs_xindex2d()
## + Cartesian-product index
vvpp_def(
'ccs_xindex2d',
Pars => "$INDX which(Ndims,Nnz); $INDX a(Na); $INDX b(Nb); $INDX\ [o]ab(Nab); $INDX\ [o]nab()",
PMCode=> q{
sub PDL::ccs_xindex2d {
my ($which,$a,$b,$ab,$nab) = @_;
barf("Usage: ccs_xindex2d(which(2,Nnz), a(Na), b(Nb), [o]nab(), [o]ab(Nab)")
if ((grep {!defined($_)} @_[0..2]) || $which->ndims != 2 || $which->dim(0) < 2);
&PDL::_ccs_xindex2d_int($which,$a,$b, ($ab//=null), ($nab//=null));
return wantarray ? ($ab,$nab) : $ab->reshape($nab->sclr);
}
},
RedoDimsCode => q{
if ($SIZE(Nab) < 0) {
if ($PDL(nab)->nvals > 0) {
$SIZE(Nab) = $nab();
} else {
$SIZE(Nab) = ($SIZE(Na)*$SIZE(Nb)) < $SIZE(Nnz) ? ($SIZE(Na)*$SIZE(Nb)) : $SIZE(Nnz);
}
}
},
Code => q{
CCS_Indx a_min=0, a_max=$SIZE(Nnz);
CCS_Indx a_lb, a_ub, a_ubmax=a_max;
CCS_Indx b_min, b_max, b_lb;
CCS_Indx abi = 0;
#if 0
/*-- DEBUG --*/
CCS_Indx size_nnz = $SIZE(Nnz);
CCS_Indx size_na = $SIZE(Na);
CCS_Indx size_nb = $SIZE(Nb);
CCS_Indx size_nab = $SIZE(Nab);
printf("Nnz=%td, Na=%td [%td:%td], Nb=%td [%td:%td], Nab=%td\n", size_nnz, size_na,$a(Na=>0),$a(Na=>size_na-1), size_nb,$b(Nb=>0),$b(Nb=>size_nb-1), size_nab);
#endif
loop (Na) %{
a_ubmax = a_max;
$LB('$a()', '$which(Ndims=>0,Nnz=>$_)', 'a_min','a_max', 'a_lb',ubmaxvar=>'a_ubmax');
//printf("a:LB(a=%td,min=%td,max=%td)=%td --> %td (ubmax=%td)\n", $a(),a_min,a_max,a_lb, $which(Ndims=>0,Nnz=>a_lb), a_ubmax); fflush(stdout);
if ($which(Ndims=>0,Nnz=>a_lb) != $a()) { a_min=a_lb; continue; }
//
$LB('$a()+1', '$which(Ndims=>0,Nnz=>$_)', 'a_lb' ,'a_ubmax', 'a_ub');
if ($which(Ndims=>0,Nnz=>a_ub) <= $a()) ++a_ub;
//printf("a:UB(a=%td,min=%td,max=%td)=%td --> %td\n", $a(),a_lb,a_ubmax,a_ub, $which(Ndims=>0,Nnz=>a_ub)); fflush(stdout);
//
b_min = a_lb;
b_max = a_ub;
loop (Nb) %{
if (b_min >= b_max) break;
//printf("+ b:LB(a=%td,b=%td,min=%td,max=%td)=", $a(),$b(),b_min,b_max); fflush(stdout);
$LB('$b()', '$which(Ndims=>1,Nnz=>$_)', 'b_min','b_max', 'b_lb');
//printf("%td --> %td", b_lb, $which(Ndims=>1,Nnz=>b_lb));
if ($which(Ndims=>1,Nnz=>b_lb) == $b()) {
//printf(" *[%td]", abi); fflush(stdout);
$ab(Nab=>abi) = b_lb;
++abi;
++b_lb;
if (abi >= $SIZE(Nab)) break;
}
b_min = b_lb;
//printf("\n"); fflush(stdout);
%}
if (abi >= $SIZE(Nab)) break;
if (a_ub < a_max) a_min = a_ub;
%}
$nab() = abi;
for ( ; abi < $SIZE(Nab); ++abi) {
$ab(Nab=>abi) = -1;
}
},
Doc=><<'EOD'
Compute indices along dimension C<NNz> of $which() corresponding to any combination
of values in the Cartesian product of $a() and $b(). Appropriate for indexing a
2d sparse encoded PDL with non-missing entries at $which() via the ND-index piddle
$a-E<gt>slice("*1,")-E<gt>cat($b)-E<gt>clump(2)-E<gt>xchg(0,1), i.e. all pairs $ai,$bi with $ai in $a()
and $bi in $b(). $a() and $b() values must be be sorted in ascending order
In list context, returns a list ($ab,$nab), where $nab() is the number of indices found,
and $ab are those C<Nnz> indices. In scalar context, trims the output vector $ab() to $nab()
elements.
EOD
);
##======================================================================
## Debugging
##======================================================================
pp_addpm(<<'EOPM');
=pod
=head1 Debugging Utilities
=cut
EOPM
##------------------------------------------------------
## ccs_dump_which()
## + prints a text dump of an index
pp_def(
'ccs_dump_which',
Pars => "$INDX which(Ndims,Nnz)",
OtherPars => 'SV *HANDLE; char *fmt; char *fsep; char *rsep',
PMCode=> q{
sub PDL::ccs_dump_which {
my ($which,$fh,$fmt,$fsep,$rsep) = @_;
$fmt = '%td' if (!defined($fmt) || $fmt eq '');
$fsep = " " if (!defined($fsep) || $fsep eq '');
$rsep = "$/" if (!defined($rsep) || $rsep eq '');
$fh = \*STDOUT if (!defined($fh));
&PDL::_ccs_dump_which_int($which,$fh,$fmt,$fsep,$rsep);
}
},
Code => q{
CCS_Indx dimi, sizeNdims=$SIZE(Ndims);
char *fmt_str = $COMP(fmt);
char *fsep_str = $COMP(fsep);
char *rsep_str = $COMP(rsep);
PerlIO *pio;
IO *io;
/*-- get PerlIO from SV (lifted from _rasc() n PDL_SRC_ROOT/IO/Misc/misc.pd) --*/
io = sv_2io($COMP(HANDLE));
if (!io || !(pio = IoIFP(io))) {
croak("can\'t get PerlIO pointer from HANDLE");
}
loop (Nnz) %{
PerlIO_printf(pio, fmt_str, $which(Ndims=>0));
for (dimi=1; dimi<sizeNdims; dimi++) {
PerlIO_puts(pio,fsep_str);
PerlIO_printf(pio, fmt_str, $which(Ndims=>dimi));
}
PerlIO_puts(pio,rsep_str);
%}
},
Doc=><<'EOD'
Print a text dump of an index PDL to the filehandle C<HANDLE>, which default to C<STDUT>.
C<$fmt> is a printf() format to use for output, which defaults to "%td".
C<$fsep> and C<$rsep> are field-and record separators, which default to
a single space and C<$/>, respectively.
EOD
);
##======================================================================
## Footer Administrivia
##======================================================================
##------------------------------------------------------
## pm additions: footer
pp_addpm(<<'EOPM');
##---------------------------------------------------------------------
=pod
=head1 ACKNOWLEDGEMENTS
Perl by Larry Wall.
PDL by Karl Glazebrook, Tuomas J. Lukka, Christian Soeller, and others.
=cut
##----------------------------------------------------------------------
=pod
=head1 KNOWN BUGS
Probably many.
=cut
##---------------------------------------------------------------------
=pod
=head1 AUTHOR
Bryan Jurish E<lt>moocow@cpan.orgE<gt>
=head2 Copyright Policy
Copyright (C) 2007-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)
=cut
EOPM
# Always make sure that you finish your PP declarations with
# pp_done
pp_done();
##----------------------------------------------------------------------