package Statistics::QMethod::QuasiNormalDist;

use warnings;
use strict;
use POSIX;
use base 'Exporter';
use Carp;
our @EXPORT = qw/get_q_dist/;
our $VERSION = '0.01';

# dispatch table for dealing with rounding errors theoretically this
# could all be memoised

my %difs = (
    0 => sub { return },
    -2 => sub {
        my ($res, $N) = @_;
        $res->[floor($#$res/2)] = $res->[floor($#$res/2)] + 2;
    },
    '-1' => sub {
        my ($res, $N) = @_;
        $res->[floor($#$res/2)] = $res->[floor($#$res/2)] + 1;
    },
    1 => sub {
        my ($res, $N) = @_;
        $res->[floor($#$res/2)] = $res->[floor($#$res/2)]   + 1 ;
        $res->[floor($#$res/2)-1] = $res->[floor($#$res/2)] - 2 ;
        $res->[floor($#$res/2)+1] = $res->[floor($#$res/2)] - 2 ;

    },
    2 => sub {
        my ($res, $N) = @_;
        $res->[floor($#$res/2)-1] = $res->[floor($#$res/2)] - 1 ;
        $res->[floor($#$res/2)+1] = $res->[floor($#$res/2)] - 1 ;
    },
    3 => sub {
        my ($res, $N) = @_;
        $res->[floor($#$res/2)-1] = $res->[floor($#$res/2)] - 2 ;
        $res->[floor($#$res/2)+1] = $res->[floor($#$res/2)] - 2 ;
        $res->[floor($#$res/2)] = $res->[floor($#$res/2)]   + 1 ;
    },
    4 => sub {
        my ($res, $N) = @_;
        $res->[floor($#$res/2)-1] = $res->[floor($#$res/2)] - 2 ;
        $res->[floor($#$res/2)+1] = $res->[floor($#$res/2)] - 2 ;
    },
    5 => sub {
        my ($res, $N) = @_;
        $res->[floor($#$res/2)-1] = $res->[floor($#$res/2)] - 2 ;
        $res->[floor($#$res/2)+1] = $res->[floor($#$res/2)] - 2 ;
        $res->[floor($#$res/2)] = $res->[floor($#$res/2)]   - 1 ;
    },
    6 => sub {
        my ($res, $N) = @_;
        $res->[floor($#$res/2)-1] = $res->[floor($#$res/2)] - 3 ;
        $res->[floor($#$res/2)] = $res->[floor($#$res/2)]   - 3 ;
    },
    7 => sub {
        my ($res, $N) = @_;
        $res->[floor($#$res/2)-1] = $res->[floor($#$res/2)] - 3 ;
        $res->[floor($#$res/2)+1] = $res->[floor($#$res/2)] - 3 ;
        $res->[floor($#$res/2)] = $res->[floor($#$res/2)]   - 1 ;
    },
    8 => sub {
        my ($res, $N) = @_;
        $res->[floor($#$res/2)-1] = $res->[floor($#$res/2)] - 3 ;
        $res->[floor($#$res/2)+1] = $res->[floor($#$res/2)] - 3 ;
        $res->[floor($#$res/2)] = $res->[floor($#$res/2)]   - 2 ;
    },
    9 => sub {
        my ($res, $N) = @_;
        $res->[floor($#$res/2)-1] = $res->[floor($#$res/2)] - 4 ;
        $res->[floor($#$res/2)+1] = $res->[floor($#$res/2)] - 4 ;
        $res->[floor($#$res/2)] = $res->[floor($#$res/2)]   - 1 ;
    },
    10 => sub {
        my ($res, $N) = @_;
        $res->[floor($#$res/2)-1] = $res->[floor($#$res/2)] - 4 ;
        $res->[floor($#$res/2)+1] = $res->[floor($#$res/2)] - 4 ;
        $res->[floor($#$res/2)] = $res->[floor($#$res/2)]   - 2 ;
    },
    11 => sub {
        my ($res, $N) = @_;
        $res->[floor($#$res/2)-1] = $res->[floor($#$res/2)] - 4 ;
        $res->[floor($#$res/2)+1] = $res->[floor($#$res/2)] - 4 ;
        $res->[floor($#$res/2)] = $res->[floor($#$res/2)]   - 3 ;
    },
    12 => sub {
        my ($res, $N) = @_;
        $res->[floor($#$res/2)-1] = $res->[floor($#$res/2)] - 5 ;
        $res->[floor($#$res/2)+1] = $res->[floor($#$res/2)] - 5 ;
        $res->[floor($#$res/2)] = $res->[floor($#$res/2)]   - 2 ;
    },
    13 => sub {
        my ($res, $N) = @_;
        $res->[floor($#$res/2)-1] = $res->[floor($#$res/2)] - 5 ;
        $res->[floor($#$res/2)+1] = $res->[floor($#$res/2)] - 5 ;
        $res->[floor($#$res/2)] = $res->[floor($#$res/2)]   - 3 ;
    },
    14 => sub {
        my ($res, $N) = @_;
        $res->[floor($#$res/2)-1] = $res->[floor($#$res/2)] - 5 ;
        $res->[floor($#$res/2)+1] = $res->[floor($#$res/2)] - 5 ;
        $res->[floor($#$res/2)] = $res->[floor($#$res/2)]   - 4 ;
    },
    15 => sub {
        my ($res, $N) = @_;
        $res->[floor($#$res/2)-1] = $res->[floor($#$res/2)] - 6 ;
        $res->[floor($#$res/2)+1] = $res->[floor($#$res/2)] - 6 ;
        $res->[floor($#$res/2)] = $res->[floor($#$res/2)]   - 3 ;
    },
    16 => sub {
        my ($res, $N) = @_;
        $res->[floor($#$res/2)-1] = $res->[floor($#$res/2)] - 6 ;
        $res->[floor($#$res/2)+1] = $res->[floor($#$res/2)] - 6 ;
        $res->[floor($#$res/2)] = $res->[floor($#$res/2)]   - 4 ;
    },
    17 => sub {
        my ($res, $N) = @_;
        $res->[floor($#$res/2)-1] = $res->[floor($#$res/2)] - 6 ;
        $res->[floor($#$res/2)+1] = $res->[floor($#$res/2)] - 6 ;
        $res->[floor($#$res/2)] = $res->[floor($#$res/2)]   - 5 ;
    },

);

sub get_q_dist {
    my $N = shift;
    die "do not support q samples of $N" if $N < 5 || $N > 500;
    my @weights = qw/0.0013 0.023 0.16  0.22 0.22 0.22 0.16 0.023 0.0013/;
    my @res = map { sprintf "%.0f", $_ * $N} @weights;
    @res = grep { $_ } @res; # strip the trailing and leading zeroes
    my $sum = _sum(\@res);
    my $dif = $sum-$N;
    eval  {
        $difs{$dif}->(\@res, $N);
    };
    if ($@) {
        warn"oops, dif was $dif\n"
    };
    _massage_tails(\@res,$N);
    _massage_centre(\@res, $N);
    _check_dist(\@res,$N);
    return \@res;
}

sub _sum {
    my $ary = shift;
    carp "needs an array referenc" if ref $ary ne 'ARRAY';
    my $sum = 0;
    map { $sum = $sum+$_} @$ary;
    return $sum;
}

sub _massage_centre {
    my ($res, $N) = @_;
    my $mid = floor($#$res/2);
    if ($res->[$mid-1] == $res->[$mid] && $res->[$mid] == $res->[$mid+1]) {
        $res->[$mid-1]--;
        $res->[$mid+1]--;
        $res->[$mid] +=2;
    if ($N == 5) { # special case
        shift @$res;
        pop @$res;
        $res->[0]=1;
        $res->[$#$res]=1;
    }
    }
    elsif ($res->[$mid] < $res->[$mid + 1 ]) {
        $res->[$mid-1]--;
        $res->[$mid+1]--;
        $res->[$mid] +=2;
    }
}

sub _massage_tails {
    my ($res, $N) = @_;
    if ($res->[0] !=1) {
        if ($res->[0] == 2) {
            $res->[0]--;
            $res->[1]++;
            $res->[$#$res]--;
            $res->[$#$res-1]++;

        }
        elsif ($res->[0] >= 3) {
            $res->[0]--;
            $res->[$#$res]--;
            push @$res,1;
            unshift @$res,1;
        }
    }
}

sub _check_dist {
    my ($res, $N) = @_;
    my $mid = floor($#$res/2);
    my @top = @$res[0 .. $mid];
    my @bottom = reverse @$res;
    @bottom = @$res[0 .. $mid];
    # check dist is symmetrical
    die "$N distrib is not symmetrical" if join ('', @top) != join ('',@bottom);
    my $n = $bottom[0];
    foreach (@bottom[1 .. $#bottom]) {
        # check that each half of the distribution is quasi-normal
        die "the distribution is not quasi normal for $N" if $_ < $n;
    }
}

=head1 Statistics::QMethod::QuasiNormalDist

=head1 SUMMARY

    use Statistics::QMethod::QuasiNormalDist;
    my @dist = get_q_dist(50); # quasi normalo distribution for a 50
                               # item Q sample

Q methodology requires the generation of a quasi normal distribution
in order to provide a ranking method for subjects in a Q sort.  This
module generates these for Q samples between 5 and 500 items long.

=head1 EXPORTS

    get_q_dist($N)

where N is the number of items in a q sample.


=head1 REFERENCES

Stephenson, W. (1953). The study of behaviour: q.technique and its
methodology. Chicago: University of Chicago Press.

=head1 Author

Kieren Diment <zarquon@cpan.org>

=head1 LICENSE

This software can be redistributed under the same conditions as perl itself.

=cut

1;