#!perl -w

=pod

=head1 NAME

Neuron::Statistics::StdDev - Calculate standard deviation

=head1 SYNOPSIS

Neuron::Statistics::StdDev provides a function to calculate the standard
deviation of a number of morphometrics.

=head1 FUNCTIONS

In all the following functions, the optional $parameters I<min> and I<max> define
limits in the $parameter specified by I<which>: branch order ("order"), path
distance ("distance") or euclidian distance ("edistance"); for example,
numSegments(\@myneuron, 2, 4, order) would only consider segments of branch
orders 2 through 4, inclusively. If $parameters are given at all, both I<min>
and I<max> must be given.

=over 0

=item *
stdDev(I<neuron>, I<$parameter> [, I<min>, I<max>, I<which>]);

return standard deviation of I<$parameter>.

=back

=head1 AUTHOR

Gunther Schmidl <gschmidl@gmx.at>

=cut

package Neuron::Statistics::StdDev;
use strict;
use warnings;

use Neuron::Morphometrics::Bifurcation;
use Neuron::Morphometrics::Segment;
use Neuron::Morphometrics::Tools;

# prepare variable export

use vars qw(@ISA @EXPORT @EXPORT_OK $VERSION);
require Exporter;
@ISA = qw(Exporter);
@EXPORT = qw(stdDev);
@EXPORT_OK = qw();
$VERSION = 1.00;






#==============================================================================
#
# stdDev(neuron, parameter [, min, max, which]) - return standard deviation of
# a parameter from the following list:
#
# btangle         bifurcation torque angle     (bifurcation)
# caangle         compartment amplitude angle  (bifurcation)
# contraction     segment contraction          (segment)
# distance        path distance from soma      (segment)
# dratio          parent/child ratio           (bifurcation)
# edistance       euclidian distance from soma (segment)
# elength         euclidian length             (segment)
# endd            final diameter               (segment)
# fragmentation   number of compartments       (segment)
# initiald        initial diameter             (segment)
# kidratio        daughter ratio               (bifurcation)
# length          path length                  (segment)
# order           branch order                 (segment)
# pc              power correction             (bifurcation)
# rall            Rall's value                 (bifurcation)
# saangle         segment amplitude angle      (bifurcation)
# taper           segment taper                (segment)
# tilt            tilt                         (bifurcation)
#
#==============================================================================

sub stdDev
{
    # get $parameters
    my @arr = @{$_[0]};
    my $param = $_[1];
    my $min = (defined($_[2]) ? $_[2] : 0);
    my $max = (defined($_[3]) ? $_[3] : maxBranchOrder(\@arr));
    my $wch = (defined($_[4]) ? $_[4] : "order");

    my ($temp, @worklist);
    my ($mean, $stddev, $count) = (0, 0, 0);

    # run a morphometric function if necessary
    for ($param)
    {
     /btangle/       and do
        { $temp = bifurcationTorqueAngle(\@arr, $min, $max, $wch); last; };
     /caangle/       and do
        { $temp = compartmentAmplitudeAngle(\@arr, $min, $max, $wch); last; };
     /contraction/   and do
        { $temp = contraction(\@arr, $min, $max, $wch); last; };
     /dratio/        and do
        { $temp = diameterRatio(\@arr, $min, $max, $wch); last; };
     /elength/       and do
        { $temp = avgSegmentEuclidianLength(\@arr, $min, $max, $wch); last; };
     /fragmentation/ and do
        { $temp = avgFragmentation(\@arr, $min, $max, $wch); last; };
     /length/        and do
        { $temp = avgSegmentPathLength(\@arr, $min, $max, $wch); last; };
     /kidratio/      and do
        { $temp = avgDaughterRatio(\@arr, $min, $max, $wch); last; };
     /pc/            and do
        { $temp = powerCorrection(\@arr, $min, $max, $wch); last; };
     /rall/          and do
        { $temp = rall(\@arr, $min, $max, $wch); last; };
     /saangle/       and do
        { $temp = segmentAmplitudeAngle(\@arr, $min, $max, $wch); last; };
     /tilt/          and do
        { $temp = tilt(\@arr, $min, $max, $wch); last; };
    }

    # put required elements on work list
    if ($param eq "contraction" || $param eq "elength" || $param eq "endd" || $param
    eq "edistance" || $param eq "distance" || $param eq "taper" || $param eq
    "length" || $param eq "order" || $param eq "initiald" || $param eq
    "fragmentation")
    {
        foreach my $seg (@SEGLIST)
        {
            # check boundaries
            next if (@{$seg}[0]->{"$wch"} < $min) ||
                    (@{$seg}[0]->{"$wch"} > $max);
            my $tmp = @{$seg}[0]->{"$param"};
            push @worklist, $tmp unless ($tmp == -1);
        }
    }
    else
    {
        foreach my $fork (@FORKLIST)
        {
            # check boundaries
            next if ($fork->{"$wch"} < $min) ||
                    ($fork->{"$wch"} > $max);
            my $tmp = $fork->{"$param"};
            push @worklist, $tmp unless ($tmp == -1);
        }
    }
    
    # calculate mean
    foreach my $val (@worklist)
    {
        $count++;
        $mean += $val;
    }
    return 0 if ($count == 0);
    $mean /= $count;
    
    # calculate standard deviation
    foreach my $val (@worklist)
    {
        $stddev += (($val - $mean) ** 2);
    }
    $stddev = sqrt($stddev / $count);
    
    return $stddev;
}
