File Coverage

blib/lib/PDL/Algorithm/Center.pm

Criterion	Covered	Total	%
statement	214	229	93.4
branch	56	82	68.2
condition	24	40	60.0
subroutine	33	34	97.0
pod	2	2	100.0
total	329	387	85.0

package PDL::Algorithm::Center;

# ABSTRACT: Various methods of finding the center of a sample

use strict;
use warnings;

require v5.10;

use feature 'state';

our $VERSION = '0.15';

use Carp;

use Try::Tiny;
use Safe::Isa;
use Ref::Util qw< is_arrayref is_ref is_coderef is_hashref >;

use Hash::Wrap ( { -as => '_wrap_hash' } );

use PDL::Algorithm::Center::Failure ':all';

use PDL::Algorithm::Center::Types -all;
use Types::Standard qw[ Optional ArrayRef Undef CodeRef Num ];
use Types::Common::Numeric -types;
use Type::Params qw[ compile_named ];

use PDL::Lite ();

use Exporter 'import';

our @EXPORT_OK = qw[ sigma_clip iterate ];

# as of Types::Standard 1.003_003, Bool coerces any value into a boolean,
# which we don't want.
use constant Bool => Types::Standard::Bool->no_coercions();

sub _weighted_mean_center {
    my ( $coords, $mask, $weight, $total_weight ) = @_;

    my $wmask = $mask * $weight;

    $total_weight //= $wmask->dsum;

    iteration_empty_failure->throw( 'weighted mean center: all elements excluded or sum(weight) == 0' )
      if $total_weight == 0;

    return ( $coords * $wmask->dummy( 0 ) )->xchg( 0, 1 )->dsumover / $total_weight;
}

sub _distance {

    my ( $prev, $current ) = @_;

    return sqrt( ( ( $prev->center - $current->center )**2 )->dsum );
}

sub _sigma_clip_initialize {    ## no critic (Subroutines::ProhibitManyArgs)

    my ( $init_clip, undef, $coords, $mask, $weight, $current, $work ) = @_;

    # initialize the sigma_clip specific fields the docs promise will be there
    $current->{dist} = undef;
    $current->{clip} = $init_clip;

    my $r2 = $work->{r2} = PDL->null;
    $r2 .= ( ( $coords - $current->center )**2 )->dsumover;

    $mask *= ( $r2 <= $init_clip**2 )
      if defined $init_clip;

    my $wmask = $work->{wmask} = $mask * $weight;

    $current->total_weight( $wmask->dsum );
    $current->nelem( $mask->sum );

    if ( $current->total_weight == 0 ) {
        $current->{sigma} = undef;
        iteration_empty_failure->throw( 'sigma_clip initialize: all elements excluded or sum(weight) == 0',
        );
    }

    $current->{sigma} = sqrt( ( $wmask * $r2 )->dsum / $current->total_weight );

    return;
}

sub _sigma_clip_calc_wmask {    ## no critic (Subroutines::ProhibitManyArgs)

    my ( $nsigma, $coords, $mask, $weight, $iter, $work ) = @_;

    my $r2 = $work->{r2};
    $r2 .= ( ( $coords - $iter->center )**2 )->dsumover;

    $iter->clip( $nsigma * $iter->sigma );

    $mask *= $r2 < $iter->clip**2;

    $iter->total_weight( ( $mask * $weight )->dsum );
    $iter->nelem( $mask->sum );

    my $wmask = $work->{wmask};
    $wmask .= $mask * $weight;

    iteration_empty_failure->throw( 'sigma_clip calc_wmask: all elements excluded or sum(weight) == 0' )
      if $iter->total_weight == 0;

    $iter->sigma( sqrt( ( $wmask * $r2 )->dsum / $iter->total_weight ) );

    return;
}

sub _sigma_clip_is_converged {    ## no critic (Subroutines::ProhibitManyArgs)

    my ( undef, $dtol, undef, undef, undef, $prev, $current ) = @_;

    $current->{dist} = undef;

    # stop if standard deviations and centers haven't changed

    if ( $current->sigma == $prev->sigma
        && PDL::all( $current->center == $prev->center ) )
    {

        $current->dist( _distance( $prev, $current ) )
          if defined $dtol;

        return 1;
    }

    # or, if a tolerance was defined, stop if distance from old
    # to new centers is less than the tolerance.
    if ( defined $dtol ) {
        $current->dist( _distance( $prev, $current ) );
        return 1 if $current->dist <= $dtol;
    }

    return 0;
}

sub _sigma_clip_log_iteration {

    my $iter = shift;

    # iter n clip sigma x y
    # xxxx xxxxxxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx

    my $ncoords = $iter->center->nelem;
    if ( $iter->iter == 0 ) {

        my $fmt = '%4s %7s' . ' %10s' x ( 3 + $ncoords ) . "\n";

        printf $fmt, @$_
          for [ qw{ iter  nelem    weight      clip      sigma   }, map { 'q$_' } 1 .. $ncoords, ],
          [ qw{ ---- ------- ---------- ---------- ----------}, ( '----------' ) x $ncoords, ];
    }

    my @fmt = ( '%4d', '%7d', ( '%10.6g' ) x ( 3 + $ncoords ) );
    $fmt[3] = '%10s' if !defined $iter->clip;

    printf(
        join( q{ }, @fmt ) . "\n",
        $iter->iter,  $iter->nelem, $iter->total_weight, $iter->clip // 'undef',
        $iter->sigma, $iter->center->list,
    );
}



























































































































































































































































































sub _to_json {
    my $obj = shift;

    require Data::Visitor::Tiny;
    require Data::Clone;
    require Scalar::Util;

    # remove outer blessed hash and clone everything except for the
    # objects.
    my $hash = Data::Clone::clone( { %{$obj} } );

    # now replace the included objects
    Data::Visitor::Tiny::visit(
        $hash,
        sub {
            return unless Scalar::Util::blessed( my $object = ${ $_[1] } );
            ${ $_[1] }
              = defined( $object->can( 'TO_JSON' ) )
              ? $object->TO_JSON
              : $object->isa( 'PDL' ) ? $object->dims == 0
                  ? $object->unpdl->[0]
                  : $object->unpdl
              : "$object";
        },
        $hash,
    );
    return $hash;
}

use Hash::Wrap ( {
        -as    => '_new_iteration',
        -class => 'PDL::Algorithm::Center::Iteration',
        -clone => sub {
            my $hash = shift;

            ## no critic (BuiltinFunctions::ProhibitComplexMappings)
            return {
                map {
                    my $value = $hash->{$_};
                    $value = $value->copy if $value->$_isa( 'PDL' );
                    ( $_, $value )
                } keys %$hash,
            };
        },
        -methods => { TO_JSON => \&_to_json },
    },
    {
        -as      => '_return_iterate_results',
        -class   => 'PDL::Algorithm::Center::Iterate::Results',
        -methods => { TO_JSON => \&_to_json },
    } );

sub sigma_clip {

    state $check = compile_named(
        center      => Optional [ ArrayRef [ Num | Undef ] | Center | CodeRef ],
        clip        => Optional [PositiveNum],
        coords      => Optional [Coords],
        dtol        => PositiveNum,
        iterlim     => Optional [PositiveInt],
        log         => Optional [ Bool | CodeRef ],
        mask        => Optional [ Undef | Piddle_min1D_ne ],
        save_mask   => Optional [Bool],
        save_weight => Optional [Bool],
        nsigma      => PositiveNum,
        weight      => Optional [ Undef | Piddle_min1D_ne ],
    );

    my $opt;
    my @argv = @_;
    try {
        my %opt = %{ $check->( @argv ); };
        $opt = _wrap_hash( \%opt );
    }
    catch {
        parameter_failure->throw( $_ );
    };

    $opt->{iterlim} //= 10;

    if ( defined $opt->{log} && !is_coderef( $opt->log ) ) {
        $opt->{log} = $opt->log ? \&_sigma_clip_log_iteration : undef;
    }

    #---------------------------------------------------------------

    # now, see what kind of data we have, and ensure that all dimensions
    # are consistent

    if ( defined $opt->{coords} ) {

        for my $name ( 'mask', 'weight' ) {

            my $value = $opt->{$name};
            next unless defined $value;

            parameter_failure->throw( "<$name> must be a 1D piddle if  is specified" )
              if $value->ndims != 1;

            my $nelem_c = $value->getdim( -1 );
            my $nelem_p = $opt->coords->getdim( -1 );

            parameter_failure->throw(
                "number of elements in <$name> ($nelem_p) ) must be the same as in  ($nelem_c)", )
              if $nelem_c != $nelem_p;
        }

    }

    elsif ( defined $opt->{weight} ) {

        $opt->{coords} = $opt->weight->ndcoords( PDL::indx );

        if ( defined $opt->{mask} ) {
            parameter_failure->throw( "mask must have same shape as weight\n" )
              if $opt->mask->shape != $opt->weight->shape;

            $opt->mask( $opt->mask->flat );
        }

        $opt->weight( $opt->weight->flat );
    }

    else {

        parameter_failure->throw( 'must specify one of  or ' );
    }

    my ( $ndims ) = $opt->coords->dims;

    if ( defined $opt->{center} && is_arrayref( $opt->center ) ) {

        my $icenter = PDL->pdl( @{ $opt->center } );

        parameter_failure->throw( " must have $ndims elements" )
          unless $icenter->nelem == $ndims;

        my $defined = PDL->pdl( map { defined } @{ $opt->center } );

        if ( $defined->not->any ) {

            $icenter = $icenter->where( $defined )->sever;

            $opt->{center} = sub {

                my ( $coords, $wmask, $weight ) = @_;
                my $center = _weighted_mean_center( $coords, $wmask, $weight );
                $center->where( $defined ) .= $icenter;

                return $center;
            };
        }
    }
    else {

        $opt->{center} //= \&_weighted_mean_center;
    }

    my $nsigma = delete $opt->{nsigma};
    $opt->{calc_wmask} //= sub {
        _sigma_clip_calc_wmask( $nsigma, @_ );
        return;
    };

    $opt->{calc_center} //= sub {
        my ( $coords, $mask, $weight, $iter ) = @_;

        _weighted_mean_center( $coords, $mask, $weight, $iter->total_weight );
    };

    my ( $clip, $dtol ) = delete @{$opt}{ 'clip', 'dtol' };
    $opt->{initialize} = sub {
        _sigma_clip_initialize( $clip, $dtol, @_ );
    };

    $opt->{is_converged} = sub {
        _sigma_clip_is_converged( $clip, $dtol, @_ );
    };

    delete @{$opt}{ grep { !defined $opt->{$_} } keys %$opt };

    iterate( %$opt );
}
















































































































































































































































































































































































































































































































































sub iterate {

    state $check = compile_named(
        center       => Center | CodeRef,
        initialize   => CodeRef,
        calc_center  => CodeRef,
        calc_wmask   => CodeRef,
        is_converged => CodeRef,
        coords       => Coords,
        iterlim      => PositiveInt,
        log          => Optional [CodeRef],
        mask         => Optional [Piddle1D_ne],
        weight       => Optional [Piddle1D_ne],
        save_mask    => Optional [Bool],
        save_weight  => Optional [Bool],
    );

    my $opt = _wrap_hash( $check->( @_ ) );

    $opt->{log}         //= undef;
    $opt->{save_mask}   //= 0;
    $opt->{save_weight} //= 0;
    $opt->{mask}        //= undef;
    $opt->{weight}      //= undef;

    my ( $ndims, $nelem ) = $opt->coords->dims;

    parameter_failure->throw( "<$_> must have $nelem elements" )
      for grep { defined $opt->{$_} && $opt->{$_}->nelem != $nelem } qw[ mask weight ];

    $opt->weight(
        defined $opt->weight
        ? PDL::convert( $opt->weight, PDL::double )
        : PDL->ones( PDL::double, $nelem ),
    );

    my $total_weight;

    if ( defined $opt->mask ) {
        $nelem        = $opt->mask->sum;
        $total_weight = ( $opt->mask * $opt->weight )->dsum;
    }
    else {
        $opt->mask( PDL->ones( PDL::long, $nelem ) );
        $total_weight = $opt->weight->dsum;
    }

    my $mask   = $opt->mask->copy;
    my $weight = $opt->weight->copy;

    $opt->center( $opt->center->( $opt->coords, $mask, $weight, $total_weight ) )
      if is_coderef( $opt->center );

    parameter_failure->throw( " must be a 1D piddle with $ndims elements" )
      unless is_Piddle1D( $opt->center ) && $opt->center->nelem == $ndims;

    #############################################################################

    # Iterate until convergence.

    my @iteration;

    my $work = {};

    # Set up initial state

    push @iteration,
      _new_iteration( {
          center       => $opt->center,
          total_weight => $total_weight,
          nelem        => $nelem,
          iter         => 0,
      } );

    $mask   .= $opt->mask;
    $weight .= $opt->weight;

    my $iteration = 0;
    my $converged;
    my $error;

    eval {

        $opt->initialize->( $opt->coords, $mask, $weight, $iteration[-1], $work, );

        $opt->log && $opt->log->( _new_iteration( $iteration[-1] ) );

        while ( !$converged && ++$iteration <= $opt->iterlim ) {

            my $prev = $iteration[-1];

            my $current = _new_iteration( $prev );
            push @iteration, $current;

            ++$current->{iter};

            $current->total_weight( $total_weight );
            $current->nelem( $nelem );

            $mask   .= $opt->mask;
            $weight .= $opt->weight;

            $opt->calc_wmask->( $opt->coords, $mask, $weight, $current, $work );

            $current->total_weight( ( $mask * $weight )->dsum )
              unless defined $current->total_weight;
            $current->nelem( $mask->sum )
              unless defined $current->nelem;

            iteration_empty_failure->throw( 'no elements left after clip' )
              if $current->nelem == 0;

            $current->center( $opt->calc_center->( $opt->coords, $mask, $weight, $current, $work, ) );

            $converged = $opt->is_converged->( $opt->coords, $mask, $weight, $prev, $current, $work, );

            $opt->log && $opt->log->( _new_iteration( $current ) );
        }
        1;
    } or $error = $@;

    $error
      = iteration_limit_reached_failure->new( msg => "iteration limit (@{[ $opt->iterlim ]}) reached" )
      if $iteration > $opt->iterlim;

    _return_iterate_results( {
        %{ $iteration[-1] },
        ( $opt->save_mask   ? ( mask   => $mask )   : () ),
        ( $opt->save_weight ? ( weight => $weight ) : () ),
        iterations => \@iteration,
        success    => !$error,
        error      => $error,
    } );
}

1;

#
# This file is part of PDL-Algorithm-Center
#
# This software is Copyright (c) 2017 by Smithsonian Astrophysical Observatory.
#
# This is free software, licensed under:
#
#   The GNU General Public License, Version 3, June 2007
#

__END__