##-*- 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::Version;
use strict;
=pod
=head1 NAME
PDL::CCS::Ops - Low-level binary operations 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');
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
##======================================================================
##======================================================================
## Operations: Binary: ALIGN: missing-is-annihilator
##======================================================================
vvpp_def
('ccs_binop_align_block_mia',
Pars => ("\n "
.join("\n ",
"$INDX\ ixa(Ndims,NnzA); $INDX\ ixb(Ndims,NnzB); $INDX\ istate(State);",
"$INDX\ [o]nzai(NnzC); $INDX\ [o]nzbi(NnzC); $INDX\ [o]ostate(State);",
'')),
Code =>
(q(
CCS_Indx sizeNnzA=$SIZE(NnzA), sizeNnzB=$SIZE(NnzB), sizeNnzC=$SIZE(NnzC);
CCS_Indx nnzai=0, nnzbi=0,nnzbi0, nnzci=0, nnzai_nxt=0,nnzbi_nxt=0,nnzci_nxt=0;
CCS_Indx cmpme1,cmpme2;
int cmpval=0;
//
//-- initialize: parse istate() [ nnzai,nnzai_nxt, nnzbi,nnzbi_nxt, nnzci,nnzci_nxt, cmpval ]
if ($SIZE(State) >= 7) {
nnzai = $istate(State=>0);
nnzai_nxt = $istate(State=>1);
nnzbi = $istate(State=>2);
nnzbi_nxt = $istate(State=>3);
nnzci = $istate(State=>4);
nnzci_nxt = $istate(State=>5);
cmpval = $istate(State=>6);
}
//
//-- main loop: start at current nnzai,nnzbi,nnzci
for ( ; nnzai<sizeNnzA && nnzbi<sizeNnzB && nnzci_nxt<=sizeNnzC; ) {
//
//-- dispatch on cached cmpval for current index-runs (nnzai..nnzai_nxt),(nnzbi..nnzbi_nxt)
if (cmpval < 0) {
//-- CASE ixa(,ai) < ixb(,bi) : INSERT ( ixa(,ai) => (ixa(,ai) . -1) ); INCR(ai);
//-- increment ai: detect next run-length
for (nnzai=nnzai_nxt, nnzai_nxt=nnzai+1; nnzai_nxt<sizeNnzA; nnzai_nxt++) {
$CMPVEC('$ixa(NnzA=>nnzai)','$ixa(NnzA=>nnzai_nxt)','Ndims','cmpval',var1=>'cmpme1',var2=>'cmpme2');
if (cmpval != 0) break;
}
}
else if (cmpval > 0) {
//-- CASE ixa(,ai) > ixb(,bi) : INSERT ( ixb(,bi) => (-1 . ixb(,bi)) ); INCR(bi);
//-- increment bi: detect next run-length
for (nnzbi=nnzbi_nxt, nnzbi_nxt=nnzbi+1; nnzbi_nxt<sizeNnzB; nnzbi_nxt++) {
$CMPVEC('$ixb(NnzB=>nnzbi)','$ixb(NnzB=>nnzbi_nxt)','Ndims','cmpval',var1=>'cmpme1',var2=>'cmpme2');
if (cmpval != 0) break;
}
}
else {
//-- CASE ixa(,ai) == ixb(,bi) : INSERT ( ixa(,ai) => (ixa(,ai) . ixb(,bi)) ); INCR(ai); INCR(bi);
for (nnzbi0=nnzbi; nnzai<nnzai_nxt; nnzai++) {
for (nnzbi=nnzbi0; nnzbi<nnzbi_nxt; nnzbi++, nnzci++) {
$nzai(NnzC=>nnzci) = nnzai;
$nzbi(NnzC=>nnzci) = nnzbi;
}
}
//-- increment ai,bi: detect next run-lengths
for (nnzai_nxt=nnzai+1; nnzai_nxt<sizeNnzA; nnzai_nxt++) {
$CMPVEC('$ixa(NnzA=>nnzai)','$ixa(NnzA=>nnzai_nxt)','Ndims','cmpval',var1=>'cmpme1',var2=>'cmpme2');
if (cmpval != 0) break;
}
for (nnzbi_nxt=nnzbi+1; nnzbi_nxt<sizeNnzB; nnzbi_nxt++) {
$CMPVEC('$ixb(NnzB=>nnzbi)','$ixb(NnzB=>nnzbi_nxt)','Ndims','cmpval',var1=>'cmpme1',var2=>'cmpme2');
if (cmpval != 0) break;
}
}
//
//-- compare current index-run values
$CMPVEC('$ixa(NnzA=>nnzai)','$ixb(NnzB=>nnzbi)','Ndims','cmpval',var1=>'cmpme1',var2=>'cmpme2');
if (cmpval < 0) { nnzci_nxt = nnzci + (nnzai_nxt-nnzai); }
else if (cmpval > 0) { nnzci_nxt = nnzci + (nnzbi_nxt-nnzbi); }
else { nnzci_nxt = nnzci + (nnzai_nxt-nnzai)*(nnzbi_nxt-nnzbi); }
} //-- end main loop
//
//-- gobble leftovers
if (nnzci_nxt < sizeNnzC) {
nnzai = nnzai_nxt = sizeNnzA;
nnzbi = nnzbi_nxt = sizeNnzB;
nnzci_nxt = nnzci;
}
//
//-- save state
if ($SIZE(State) >= 7) {
$ostate(State=>0) = nnzai;
$ostate(State=>1) = nnzai_nxt;
$ostate(State=>2) = nnzbi;
$ostate(State=>3) = nnzbi_nxt;
$ostate(State=>4) = nnzci;
$ostate(State=>5) = nnzci_nxt;
$ostate(State=>6) = cmpval;
}
)),
Doc =>
(q{
Partially aligns a pair of lexicographically sorted index-vector lists C<$ixa()> and C<$ixb()>,
e.g. for block-wise incremental computation of binary operations over sparse index-encoded PDLs,
assuming missing indices correspond to annihilators.
On return, the vectors C<$nzai> and C<$nzbi> hold indices into C<NnzA> and C<NnzB>
respectively, and are constructed such that:
($ixa(,$nzai->slice("0:$nzci_max")) == $ixb(,$nzbi->slice("0:$nzci_max"))
At most C<NnzC> alignments are performed, and alignment ceases
as soon as any of the PDLs C<$ixa()>, C<$ixb()>, C<$nzai()>, or C<$nzbi()>
has been exhausted.
The parameters C<$istate()> and C<$ostate()> hold the state of the algorithm,
for incremental block-wise computation at the perl level. Each state PDL
is a 7-element PDL containing the following values:
INDEX LABEL DESCRIPTION
-----------------------------------------------------------------------
0 nnzai minimum offset in NnzA of current $ixa() value
1 nnzai_nxt minimum offset in NnzA of next $ixa() value
2 nnzbi minimum offset in NnzB of current $ixb() value
3 nnzbi_nxt minimum offset in NnzB of next $ixb() value
4 nnzci minimum offset in NnzC of current ($ixa(),$ixb()) value pair
5 nnzci_nxt minimum offset in NnzC of next ($ixa(),$ixb()) value pair
6 cmpval 3-way comparison value for current ($ixa(),$ixb()) value pair
For computation of the first block, $istate() can be safely set to C<zeroes(long,7)>.
Repetitions may occur in input index PDLs C<$ixa()> and C<$ixb()>.
If an index-match occurs on such a "run", I<all pairs> of matching values are
added to the output PDLs.
All alignments have been performed if:
$ostate(0)==$NnzA && $ostate(1)==$NnzB
B<WARNING:> this alignment method ignores index-vectors which are not present
in I<both> C<$ixa()> and C<$ixb()>, which is a Good Thing if your are feeding
the aligned values into an operation for which missing values are annihilators:
$missinga * $bval == ($missinga * $missingb) for each $bval \in $b, and
$aval * $missingb == ($missinga * $missingb) for each $aval \in $a
This ought to be the case for all operations if missing values are C<BAD> (see L<PDL::Bad>),
but might cause unexpected results if e.g. missing values are zero and the operation
in question is addition.
}),
); ##--/ccs_binop_align_block_mia
##======================================================================
## 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
No support for (pseudo)-threading.
=cut
##---------------------------------------------------------------------
=pod
=head1 AUTHOR
Bryan Jurish E<lt>moocow@cpan.orgE<gt>
=head2 Copyright Policy
All other parts 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();
##----------------------------------------------------------------------