#! perl -w

=pod

=head1 NAME

SimpleHistogram - create one-dimensional histograms from neuron data

=head1 SYNOPSIS

SimpleHistogram allows the creation of one-dimensional histograms that plot
a given value (or range of values) against the number of occurrences.

=head1 USAGE

Run SimpleHistogram without parameters to see the available command line
options.

=head1 AUTHOR

Gunther Schmidl <gschmidl@gmx.at>

=cut

use strict;
use warnings;
use POSIX;
no warnings 'recursion';

my $VERSION = "1.0";

use Neuron::SWC::Read;

use Neuron::Morphometrics::Bifurcation;
use Neuron::Morphometrics::Segment;
use Neuron::Morphometrics::Tools;

my (@arr, %hdr, %opts);






#==============================================================================
#
# Subroutine for running all the statistics necessary to set up segments with
# the required data
#
#==============================================================================

sub setup()
{
    my $temp;

    $temp = minFragmentation(\@arr)          if $opts{'m'} eq "fragmentation";
    $temp = contraction(\@arr)               if $opts{'m'} eq "contraction";
    $temp = minSegmentEuclidianLength(\@arr) if $opts{'m'} eq "elength";
    $temp = minSegmentPathLength(\@arr)      if $opts{'m'} eq "length";
    $temp = avgDaughterRatio(\@arr)          if $opts{'m'} eq "kidratio";
    $temp = bifurcationTorqueAngle(\@arr)    if $opts{'m'} eq "btangle";
    $temp = compartmentAmplitudeAngle(\@arr) if $opts{'m'} eq "caangle";
    $temp = diameterRatio(\@arr)             if $opts{'m'} eq "dratio";
    $temp = powerCorrection(\@arr)           if $opts{'m'} eq "pc";
    $temp = rall(\@arr)                      if $opts{'m'} eq "rall";
    $temp = segmentAmplitudeAngle(\@arr)     if $opts{'m'} eq "saangle";
    $temp = tilt(\@arr)                      if $opts{'m'} eq "tilt";
}






#==============================================================================
#
# Subroutine for running all the statistics necessary to set up segments with
# the required data
#
#==============================================================================

sub varlist
{
    print<<"EOV";

usage: SimpleHistogram [-l<n> -h<n> -w<var>] [-s<n> | -z<n>] [-t1]
                       [-m<var> | -m?] <inputfile>

    possible values for -m:

    SEGMENTS
    contraction     segment contraction
    distance        path distance from soma
    edistance       euclidian distance from soma
    elength         euclidian length
    endd            final diameter
    fragmentation   number of compartments
    initiald        initial diameter
    length          path length
    order           branch order
    taper           segment taper

    BIFURCATIONS
    btangle         bifurcation torque angle
    caangle         compartment amplitude angle
    dratio          parent/child ratio
    kidratio        daughter ratio
    pc              power correction
    rall            Rall's value
    saangle         segment amplitude angle
    tilt            tilt
EOV
    exit;
}






#==============================================================================
#
# Options
#
#==============================================================================

sub usage
{
    print<<"EOU";

usage: SimpleHistogram [-l<n> -h<n> -w<var>] [-s<n> | -z<n>] [-t1]
                       [-m<var> | -m?] <inputfile>

    required parameters:
    -m <var>        measurement (use -m? for a full list)
    inputfile       data in SWC format

    optional parameters:
    -w <var>        variable for limit (order, distance, or edistance) *
    -l <n>          lower limit for variable *
    -h <n>          upper limit for variable *
    -s <n>          number of bins **
    -z <n>          size of bins **
    -t1             generate tab-delimited instead of comma-delimited output
     * these three parameters must be used together
    ** these two parameters can not be used together
EOU
    exit;
}

sub getopt
{
    # get options
    my $argumentative = $_[0];
    my ($first,$rest);

    while (@ARGV && ($_ = $ARGV[0]) =~ /^-(.)(.*)/)
    {
        ($first,$rest) = ($1,$2);

        if (index($argumentative,$first) >= 0)
        {
            if ($rest ne '') { shift(@ARGV); }
            else { shift(@ARGV); $rest = shift(@ARGV); }
            $opts{$first} = $rest;
        }
        else
        {
            if ($rest ne '') { $ARGV[0] = "-$rest"; }
            else { shift(@ARGV); }
        }
    }

    # help requested?
    usage if !defined($opts{'m'});
    varlist if $opts{'m'} eq '?';

    print "Initialising...";

    # check option syntax
    my $test = 0;

    # -w, -l, -h must be used together, or not at all
    $test++ if defined($opts{'w'}); $test++ if defined($opts{'l'});
    $test++ if defined($opts{'h'});
    usage if ($test != 0 && $test != 3);

    # -s, -z must not be used together
    usage if (defined($opts{'s'}) && defined($opts{'z'}));

    # one and only one filename
    usage if !defined($ARGV[0]);
    usage if defined($ARGV[1]);

    -e $ARGV[0] || die "Error: $ARGV[0] does not exist.";

    %hdr = readHeader($ARGV[0]);
    @arr = readData($ARGV[0]);

    # set standard values if undefined
    $opts{'w'} = 'order' unless defined $opts{'w'};
    $opts{'l'} = 0 unless defined $opts{'l'};
    $opts{'h'} = maxBranchOrder(\@arr) unless defined $opts{'h'};

    $opts{'s'} = 1 unless defined($opts{'s'});
}






#==============================================================================
#
# Main program starts here
#
#==============================================================================

# print pretty header
print "SimpleHistogram $VERSION (Neuron $NEURON_VERSION)\n";

# get options and load input file
getopt('lhwsmtz');
setup();
print " done.\n";

# declare variables
my ($min, $max) = (65000, -1);
my (@worklist, @bins);
my ($numbins, $binsize);

# get maximum and minimum

if ($opts{'m'} eq "contraction" || $opts{'m'} eq "elength" || $opts{'m'} eq 
    "endd" || $opts{'m'} eq "edistance" || $opts{'m'} eq "distance" || 
    $opts{'m'} eq "taper" || $opts{'m'} eq "length" || $opts{'m'} eq "order" || 
    $opts{'m'} eq "initiald" || $opts{'m'} eq "fragmentation")
{
    foreach my $seg (@SEGLIST)
    {
        # check boundaries
        next if (@{$seg}[0]->{"$opts{'w'}"} < $opts{'l'}) ||
                (@{$seg}[0]->{"$opts{'w'}"} > $opts{'h'});
    
        # pop values on a work list, figure out minimum and maximum
        my $temp = @{$seg}[0]->{"$opts{'m'}"};
        push @worklist, $temp;
        $min = $temp if $temp < $min;
        $max = $temp if $temp > $max;
    }
}
else
{
    foreach my $fork (@FORKLIST)
    {
        # check boundaries
        next if ($fork->{"$opts{'w'}"} < $opts{'l'}) ||
                ($fork->{"$opts{'w'}"} > $opts{'h'});
    
        # pop values on a work list, figure out minimum and maximum
        my $temp = $fork->{"$opts{'m'}"};
        push @worklist, $temp;
        $min = $temp if $temp < $min;
        $max = $temp if $temp > $max;
    }    
}    

$max++ if $min == 0;

# set number of bins and bin size
if (defined($opts{'z'}))
{
    $binsize = $opts{'z'};
    $numbins = ($max - $min) / $binsize;
}
else
{
    $numbins = $opts{'s'};
    $binsize = ($max - $min) / $numbins;
}

# fill in unused bins
for (my $i = 0; $i <= $numbins; $i++)
{
    $bins[$i] = 0;
}

# get samples
print "Sampling between $opts{'w'} of $opts{'l'} and $opts{'h'} ";
print "(", ceil($numbins), " bins of size $binsize).\n";

# sort values into bins
my ($bin, $val);

while (@worklist)
{
    $val = pop(@worklist);
    $bin = ceil($numbins) - 1;
    while ($bin >= 0)
    { 
        $bins[$bin]++ if $val >= (($bin * $binsize) + $min);
        last if $val >= (($bin * $binsize) + $min);
        $bin--;
    }
}

# print to file
my $output_file = $ARGV[0];
my $DELIMITER;

if (defined($opts{'t'}))
{
    $output_file =~ s/(.*)\..*/$1.txt/;
    $DELIMITER = "\t";
}
else
{
    $output_file =~ s/(.*)\..*/$1.csv/;
    $DELIMITER = ",";
}

print "Output file is $output_file.\n";

open FILEOUT, ">$output_file" or die "Could not open output file: $!";

for (my $i = $min; $i < $numbins; $i++)
{
    if ($i == $numbins - 1)
    {
        print FILEOUT "[", $i * $binsize, ";$max]${DELIMITER}$bins[$i]\n";
    }
    else
    {
        print FILEOUT "[", $i * $binsize, ";", ($i + 1) * $binsize;
        print FILEOUT "[${DELIMITER}$bins[$i]\n";
    }
}

close FILEOUT;

print "All done.";