#! /usr/local/bin/perl
################################################################################
#
#  Let's Do Statistics
#
#  copyright (c)
#    Fritz Lipmann Institute Jena, Genome Analysis Group, 2005-2006, 2010, 2013
#    Karol Szafranski, 2007-2008
#    Karol Szafranski at UPenn Philadelphia, Center for Bioinformatics, 2004
#    Institute of Molecular Biotechnology Jena, Dept. Genome Analysis, 1999-2004
#  author
#    Karol Szafranski, karol.szafranski@leibniz-fli.de
#
################################################################################
#
#  DESCRIPTION
#
# - See function &usage for description of command line syntax
#
# - each function comes along with a description at the beginning of the code
#   block
#
################################################################################
#
#  FUNCTIONS, DATA
#
# - MAIN
#   %GlobStore
#   $ProgFile,$ProgFstump
#   %ProgParam
#   $ProgMode,%ProgOpt,@ProgArg
#
# - usage help, command line arguments, basic I/O
#   &usage
#   &AddSwitch
#   &PrepOstump
#
# - sample metrics, approximations
#   &ProgDensDistrib
#   &ProgSampleQuant
#   &ProgSampleMeas
#   &ProgRegressLin
#
# - statistical tests
#   &ProgAnova
#   &ProgChisqr4
#   &ProgChisqrk2
#   &ProgPermutOrder
#   &ProgPermutCrowd
#
# - likelihood distributions
#   &ProgBinom
#  (&ProgGaussLimit) -> &Math::Statist::GaussLimit
#  (&ProgGaussQuant) -> &Math::Statist::GaussQuant
#
# - random values, simulations
#   &ProgRandGauss
#
################################################################################


# global constants and variables
# also used by modules as a unique global data anchor
our %GlobStore;

# include path(s), includes
use strict; #use warnings;  # OK 20080123
use Cwd qw(realpath);
BEGIN {
  unshift @INC, grep{$_} split(/:+/,$ENV{KPERLPATH}||$ENV{PERLPATH}||'');
  my ($mypath) = realpath(__FILE__)=~m/(^.*)\//;
  push @INC, $mypath;
}
use MainLib::StrRegexp qw(&TimeStr);
use MainLib::Data;
use MainLib::Path;
use MainLib::Cmdline qw(&GetoptsNArgs &QueryConfirm &GetWriteHandle);
use MainLib::File;
use MainLib::Misc;
use Math::kCalc;
use Math::kRandom;
use Math::Statist;
use Math::Plot2D;
use Math::PlotImg;  # this is depending on library GD
use database::DbPlain qw(&PlainToTable);
use database::Table;


# script ID
# - program name as specified on the command line
our $ProgFile = ( split('/',__FILE__) )[-1];
our $ProgFstump=$ProgFile; $ProgFstump=~s/\.\w{1,4}$//;

# global constants (esp. default values)
our %ProgParam;
$ProgParam{shotgun}{overlap} = 45;
$ProgParam{default}{OutImgWidth} = 640;
$ProgParam{default}{OutImgRelhigh} = 0.5;
$ProgParam{default}{PermutOrderAlpha} = 0.05;
$ProgParam{default}{PermutN} = 1000;
$ProgParam{default}{ProgMode} = 'sample';


# manage I/O #####################################

# organiz I/O handles
&Unbuffer();


# command line interface #########################
# NOTE:
# - &GetoptsNArgs (via &AddSwitch) modifies global variables:
#   $ProgMode @ProgArg %ProgOpt
# - @ProgArg may be pre-filled in &GetoptsNArgs via &AddSwitch (option -fofn)

# arguments, switches, default subprogram
our $ProgMode = undef;
our @ProgArg = ();
our %ProgOpt = ();
@ProgArg = &GetoptsNArgs();
$ProgMode ||= $ProgParam{default}{ProgMode};

# eventually open LOG file
if ($ProgOpt{-log}) {
  $ProgOpt{LogFile} = ($ProgOpt{-log} ne '1')? $ProgOpt{-log} : undef;
  &LogOpen (-file=>$ProgOpt{LogFile}, -stamp=>$ProgFstump, -prog=>"$ProgFile -$ProgMode");
}
END {
  $ProgOpt{-log} and &LogClose();
}


# work flow manifold #############################

# chain to program mode without input argument(s)
if (0) { }
elsif (!@ARGV or $ProgMode=~m/^h(elp)?$/i) { &usage() }

# ensure input argument(s)
unless (@ProgArg) {
  die "ERROR: input arguments missing\n";
}

# chain to program mode (with input argument(s))
if (0) { }
elsif ($ProgMode =~ m/^anova$/i) {
  &ProgAnova (@ProgArg);
}
elsif ($ProgMode =~ m/^binom$/i) {
  &ProgBinom (@ProgArg);
}
elsif ($ProgMode =~ m/^chisqr4$/i) {
  &ProgChisqr4 (@ProgArg);
}
elsif ($ProgMode =~ m/^chisqrk2$/i) {
  &ProgChisqrk2 (@ProgArg);
}
elsif ($ProgMode =~ m/^DensDistrib$/i) {
  &ProgDensDistrib ($ProgArg[0]);
}
elsif ($ProgMode =~ m/^DensDistrib(Weighted|Wgt)$/i) {
  $ProgOpt{-SampleWgt} = 1;
  &ProgDensDistrib ($ProgArg[0]);
}
elsif ($ProgMode =~ m/^GaussLimit$/i) {
  printf "%s\n", &GaussLimit ($ProgArg[0], -debug=>$ProgOpt{-debug});
}
elsif ($ProgMode =~ m/^GaussQuant$/i) {
  printf "%s\n", &GaussQuant ($ProgArg[0], -double=>$ProgOpt{-var}{'double-sided'}, -debug=>$ProgOpt{-debug});
}
elsif ($ProgMode =~ m/^permutorder$/i) {
  &ProgPermutOrder (@ProgArg);
}
elsif ($ProgMode =~ m/^permutcrowd$/i) {
  &ProgPermutCrowd (@ProgArg);
}
elsif ($ProgMode =~ m/^RandGauss$/i) {
  &ProgRandGauss (@ProgArg);
}
elsif ($ProgMode =~ m/^RandInt$/i) {
  my $n = $ProgArg[1] || 1;
  for (my $i=0; $i<$n; ++$i) {
    printf "%s\n", &RandInt ($ProgArg[0]||1);
  }
}
elsif ($ProgMode =~ m/^RegressLin$/i) {
  &ProgRegressLin ($ProgArg[0]);
}
elsif ($ProgMode =~ m/^Sample(Weighted|Wgt)?$/i) {
  $ProgOpt{-SampleWgt} = length($1||'') ? 1 : undef;
  &ProgSampleMeas ($ProgArg[0]);
}
elsif ($ProgMode =~ m/^SampleQuant(Wgt)?$/i) {
  $ProgOpt{-SampleWgt} = length($1||'') ? 1 : undef;
  &ProgSampleQuant ($ProgArg[0]);
}
else {
  die "ERROR: unknown program mode or switch '$ProgMode'\n";
}

# exit script successfully
# cf. END blocks!
exit 0;


################################################################################
# usage help, command line arguments, basic I/O
################################################################################


sub usage {
  print "\n";
  print <<END_USAGE;
DESCRIPTION
 $ProgFile is desired for sequencing result folder and file handling.

COMMAND LINE SYNTAX
 $ProgFile  -<ModeSwitch> [-<OptionalSwitch> ...] <Arg1> [<Arg2> ...]

Arguments
---------
 See function descriptions.

path arguments:
 Path "-" resolves to STDIN or STDOUT, depending on the context.

Input Table File Format
-----------------------
 Table input files are expected to have TAB-delimited format by standard.
 Any other delimiter concept can be specified through option -DelimCol.

ModeSwitch (case-insensitive)
-----------------------------
<none>            default ModeSwitch -$ProgParam{default}{ProgMode} if program arguments are given.
                  Otherwise like ModeSwitch -help.
-anova            perform one-way ANOVA test
                  Arg1        data table, one sample per line, values separated
                              by whitespaces
-binom            calculate binomial table
                  Arg1        range / step size for n value, like '20..40:2'.
                  Arg2        range / step size for p value, like '0.02..0.80:0.02'
                  Arg3?       range / step size for x value, like '20..40',
                              default: complete range of x in 0..n
                              This will result in the specified range of x
                              displayed only. The tables, of course will be 
                              calculated based on the complete possible range
                              of x.
                  -accur=n
-chisqr4          chi square for four-field table
                  Arg1-4      fields of table
-chisqrk2         chi square for k*2-field table
                  Arg1        path of table file
-DensDistrib      case density distribution (table) for a sample
                  Arg1        path of file containing array of values
                  --xmin=F    minimal x value for lower interval beyond data
-DensDistribWgt   same as -DensDistrib, but 2nd data column of input table
                  specifies weight for each data entry
-GaussLimit       determine quantile limit as abs(delta-my)/sigma in Gaussian
                  distribution
                  Arg1        quantile
-GaussQuant       Gaussian quantile of a given multiple of sigma
                  Arg1        multiple of sigma
                  --double-sided=B
                              quantile for double-sided testing, default one-
                              sided
-h(elp)           output command line syntax description and exit
-PermutCrowd      significance of item crowding, assessed using permutation
                  test
                  Arg1++      path of table file
                  --permut=N  number of permutations, default $ProgParam{default}{PermutN}
-PermutOrder      significance of item order in repeated sampling, assessed
                  using permutation test
                  Arg1++      path of table file
                  --alpha=F   alpha value, default $ProgParam{default}{PermutOrderAlpha}
                  --permut=N  number of permutations, default $ProgParam{default}{PermutN}
-RandGauss        random values following normal-distribution
                  Arg1        my
                  Arg2        sigma
                  Arg3?       size of sample, default: 1
-RandInt          random integer number
                  Arg1        maximum value, default: 1
                  Arg2?       size of sample, default: 1
-RegressLin       linear regression on a sample of (x,y) data pixels
                  Arg1        path of table file containing data pixels
                              col0=x, col1=y, col2=weight (optional)
                  -OutImg     prepare image of data and regression lines
                  --origin=B  use origin in displayed regression: 0 (=NO),
                              1 (=YES), default: both. In the generated plot
                              graph (switch -OutImg), the regression variants
                              are color-coded as blue (no origin constraint) and
                              green (no origin constraint).
                  --varaxis=S data axis which is associated with the variation,
                              possible: x, y (default)
-sample           standard calculations on a sample
                  Arg1        path of file containing array of values
                  -OutImg     prepare image of the case density distribution
                  --climax=N  direction of cumulation in conjunction with
                              -OutImg, default: +1 = with increasing x;
                              possible: -1 with decreasing x
                  --ImgGauss=B
                              add Gauss approximation function to distribution
                              plot, default: do not
                  --smoothws=F
                              smoothen sample plot in created image by calc-
                              ulating the mean value in a sliding window with
                              the specified size
                  --xmin=F    minimal x value for lower interval beyond data of
                              case density distribution plot
-SampleWgt        same as -sample, but 2nd data column specifies weight for
                  each data entry
-SampleQuant      quantile plot (table) for a sample, i.e. a plot of the
                  cumulative case density distribution (cmp. mode switch
                  -DensDistrib)
                  Arg1        path of file containing array of values
                  --abs=B     plot cumulation of absolute number of cases rather
                              than the default normalized case frequency
                              (summing to 1.0)
                  --climax=N  direction of cumulation, default: +1 = with
                              increasing x; possible: -1 with decreasing x
-SampleQuantWgt   same as -SampleQuant, but 2nd data column specifies weight
                  for each data entry

OptionalSwitch (case-insensitive)
---------------------------------
switch argument types: B:=boolean, F:=floating point/scientific, N:=integer,
S:=string, X:=varying type.

-accur=n          numerical accuracy (number of digits) in probability value
                  calculations
-debug(=N)        print debug protocol to STDOUT/STDERR. A debug depth value
                  may be specified via switch parameter N in syntax '-debug=N'.
-DelimCol=S       regexp for data field splitting, default: '\\t'. Note that
                  the regexp shall not contain bracketed expressions.
-format=S         other than the default format scheme
-log(=S)          redirect STDOUT and STDERR to LOG file
                  S           log file path, default path worked out
                              automatically from built-in directives.
-OutDir=S         directory for multi-file output. This switch overrides any
                  directory statement provided with switch -OutStump.
-OutImg(=S)       force program to produce image output and (optionally)
                  specify output path
-OutImgRelhigh=F  relative height of an output image in relation to its width,
                  default $ProgParam{default}{OutImgRelhigh}
-OutImgTransp     turn an image's background transparent
-OutImgWidth=N    define pixel width of an output image (plot data field),
                  default: $ProgParam{default}{OutImgWidth}
-OutStump=S       path stump for multi-file output. A default is derived from
                  input file names in most cases.
-redo=N           redo batch size, for simulations etc.
-v(erbose)        print extended action protocol to STDOUT
--*               program mode-specific switches. See the descriptions there.
                  Case-sensitive!

Environment Variables
---------------------
 \$PERLPATH        primary search path for Perl package look-up
END_USAGE
  print "\n";
  exit 0;
}


# add program switches to global table (hash)
#
# INTERFACE
# - argument 1:  switch argument without leading '-'
#
# - global data:
#   $ProgMode
#   %ProgOpt  switch data which gets processed here
#
# DESCRIPTION
# - this function gets called by &MainLib::Misc::GetoptsNArgs
# - switch arguments are tested for validity. Arguments are parsed with highest
#   possible tolerance. This way, syntax errors can reported in accordance to
#   the actual switch, rather than reporting ANY syntax error.
#
sub AddSwitch {
  my $switch = shift;

  # optional switches
  if ($switch =~ m/^accur=(.+)$/i) {
    $ProgOpt{-accur} = $1;
    return;
  }
  if ($switch =~ m/^debug(=(\d+))?$/i) {
    $ProgOpt{-debug} = defined($2) ? int($2) : 1;
    return;
  }
  if ($switch =~ m/^DelimCol=(.+)$/i) {
    $ProgOpt{-DelimCol} = $1;
    return;
  }
  if ($switch =~ m/^format=(.+)$/i) {
    $ProgOpt{-format} = $1;
    return;
  }
  if ($switch =~ m/^log(=(.*))?$/i) {
    $ProgOpt{-log} = $2 ? &PathExpand($2) : 1;
    return;
  }
  if ($switch =~ m/^OutDir=(.+)$/i) {
    $ProgOpt{-OutDir} = &PathExpand ($1);
    unless (-d $ProgOpt{-OutDir}) {
      if (-e $ProgOpt{-OutDir}) {
        die sprintf "ERROR: output destination exists, but is not a directory: %s\n", $ProgOpt{-OutDir}||"''";
      }
      
      # this dialogue will also appear if there's no output do be done to any
      # file/directory
      else {
        if (int(grep{ $_ eq '-' }@ARGV)) {
          printf STDERR "creating non-existent output directory (skip dialogue in STDIN input mode)\n";
          mkdir ($ProgOpt{-OutDir});
        } else {
          printf STDERR "output directory does not exist, create?";
          if (&QueryConfirm()) {
            mkdir ($ProgOpt{-OutDir});
          }
          else { exit 1 }
        }
      }
    }
    return;
  }
  if ($switch =~ m/^OutImg(=(.+))?$/i) {
    $ProgOpt{-OutImg} = $2 ? &PathExpand($2) : '';
    return;
  }
  if ($switch =~ m/^OutImgRelhigh=([0-9\.eE+-]+)?$/i) {
    $ProgOpt{-OutImg} ||= '';
    $ProgOpt{-OutImgRelhigh} = $1;
    return;
  }
  if ($switch =~ m/^OutImgTransp(ar)?$/i) {
    $ProgOpt{-OutImg} ||= '';
    $ProgOpt{-OutImgTransp} = 1;
    return;
  }
  if ($switch =~ m/^OutImgWidth=(\d+)?$/i) {
    $ProgOpt{-OutImg} ||= '';
    $ProgOpt{-OutImgWidth} = $1;
    return;
  }
  if ($switch =~ m/^OutSt[au]mp=(.+)$/i) {
    $ProgOpt{-OutStump} = $1;
    return;
  }
  if ($switch =~ m/^redo(=(\d+))$/i) {
    $ProgOpt{-redo} = $2;
    return;
  }
  if ($switch =~ m/^v(erbose)?$/i) {
    $ProgOpt{-verbose} = 1;
    return;
  }
  if ($switch =~ m/^(?:-|var=)(\w+)[,=](.+)$/i) {
    $ProgOpt{-var}{$1} = $2;
    return;
  }

  # program mode switches
  if (defined $ProgMode) {
    die sprintf "ERROR: multiple specification of program mode or unknown switch, %s and %s\n",
      '-'.($ProgMode||"''"), '-'.($switch||"''");
  }
  else {
    $ProgMode = $switch;
    return;
  }
}


# work out output path base from source file and/or switches
#
# INTERFACE
# - global options:
#   -debug       [STD]
#   -OutDir      [STD]
#   -OutStump    [STD]
# - return val:  output path base
#
sub PrepOstump {
  my ($debug);
  my ($PathStamp, $PathRef);

  # function parameters
  $debug = $ProgOpt{-debug};

  # prepare output path base
  # 1st: $ProgOpt{-OutStump}
  # 2nd: $ProgArg[0]
  unless ($PathStamp = $ProgOpt{-OutStump}) {
    $PathRef = $ProgArg[0];
    $PathStamp = &PathChgSuffix (&PathExpand($PathRef,-pipe=>'<'), '', -last=>1);
  }

  # determine effective directory
  $ProgOpt{-OutDir} and $PathStamp =
    $ProgOpt{-OutDir} .'/'. &PathSplit($PathStamp)->{name};

  # final refinement: expand to nice, rooted path
  $PathStamp = &PathExpand ($PathStamp);

  # return path
  return $PathStamp;
}


################################################################################
# sample metrics, approximations
################################################################################


# case density distribution for a sample
#
# INTERFACE
# - argument 1: path of file that contains array of values
#
# - global options:
#   -debug      print debug protocol to STDERR
#   -SampleWgt  sample data is weighted, means: 2nd data column specifies
#               weight for each data entry
#
# DESCRIPTION
# - this is nothing more than a wrapper to library function &DistribEmpir
#
sub ProgDensDistrib {
  my ($PathTab) = @_;
  my $debug = $ProgOpt{-debug};
  my $dbg2  = $debug ? $debug-1 : undef;
  my $bWgt = $ProgOpt{-SampleWgt};

  # main data
  my $DataTabType = $bWgt ? 'AA' : 'A1';
  my $pData = &PlainToTable ($PathTab, -TabType=>$DataTabType,
    -DelimCol=>$ProgOpt{-DelimCol}, -comments=>1, -debug=>$dbg2);
  my $ValN;
  unless ($ValN = int(@$pData)) {
    die sprintf "ERROR: don't see any data in file %s\n", $PathTab||"''";
  }
  # is data really weighted?
  if ($bWgt and int (grep{ $$_[1]>0 }@$pData) != $ValN) {
    printf STDERR "WARNING: data does not seem to be weighted correctly\n";
    if (grep{ $$_[1]<0 }@$pData) {
      print  STDERR "ERROR: weight < 0 is not allowed\n"
        . "  analysing data unweighted\n";
      $pData = &TableConvert ('AA', 'A1', $pData, -debug=>$dbg2);
      undef $bWgt;
    }
  }

  # empirical distribution function ("atomic resolution histogram")
  my $pDistrib = &DistribEmpir ($pData,
    exists($ProgOpt{-var}{xmin})?(-RangeMin=>$ProgOpt{-var}{xmin}):());
  if (!$pDistrib or !int(@$pDistrib)) {
    print  STDERR "ERROR in \&DistribEmpir\n";
    &DistribEmpir ($pData, -debug=>1,
      exists($ProgOpt{-var}{xmin})?(-RangeMin=>$ProgOpt{-var}{xmin}):());
    exit 1;
  }
  foreach (@$pDistrib) {
    print  "$_->[0]\t$_->[1]\n";
  }
}


# quantile plot (cumulative case density distribution) for a sample
#
# INTERFACE
# - argument 1: path of file containing array of values
#
# - global options:
#   -debug      print debug protocol to STDERR
#   -SampleWgt  sample data is weighted, means: 2nd data column specifies
#               weight for each data entry
#   --abs       cumulation of absolute case number
#   --climax    direction of cumulation
#
sub ProgSampleQuant {
  my ($PathTab) = @_;
  my $debug = $ProgOpt{-debug};
  my $dbg2  = $debug ? $debug-1 : undef;
  my $bWgt = $ProgOpt{-SampleWgt};
  my $climax = (($ProgOpt{-var}{climax}||1)>0) ? +1 : -1;

  # main data
  my $DataTabType = $bWgt ? 'AA' : 'A1';
  my $pData = &PlainToTable ($PathTab, -TabType=>$DataTabType,
    -DelimCol=>$ProgOpt{-DelimCol}, -comments=>1, -debug=>$dbg2);
  my ($ValN, $ValWgtsum);
  unless ($ValN = int(@$pData)) {
    die sprintf "ERROR: don't see any data in file %s\n", $PathTab||"''";
  }
  # is data really weighted?
  if ($bWgt and int (grep{ $$_[1]>0 }@$pData) != $ValN) {
    printf STDERR "WARNING: data does not seem to be weighted correctly\n";
    if (grep{ $$_[1]<0 }@$pData) {
      print  STDERR "ERROR: weight < 0 is not allowed\n"
        . "  analysing data unweighted\n";
      $pData = &TableConvert ('AA', 'A1', $pData, -debug=>$dbg2);
      undef $bWgt;
    }
  }
  $ValWgtsum = $bWgt ? &Sum(map{ $$_[1] }@$pData) : $ValN;
  # force data to weighted format (= "AA")
  if (! $bWgt) {
    $pData = [map { [$_,1] } @$pData];
  }

  # empirical case density distribution ("atomic resolution histogram")
  my (%DistribDens, @DistribDist, $DistFlank);
  my $bFirst = 1;
  my $ValLast;
  foreach (sort { $a->[0]<=>$b->[0] } @$pData) {
    if (!$bFirst) {
      if ((my $dist=$_->[0]-$ValLast) > 0) {
        push @DistribDist, $dist;
      }
    }
    $DistribDens{$_->[0]} ||= 0;
    $DistribDens{$_->[0]} += $_->[1];
    $bFirst = 0;
    $ValLast = $_->[0];
  }
  $DistFlank = &Min(@DistribDist) / 2;
  my $fnormal = $ProgOpt{-var}{abs} ? 1 : $ValWgtsum;
  my $CumulWgt = ($climax>0) ? 0.0 : $ValWgtsum;
  foreach (sort { $a<=>$b } keys %DistribDens) {
    printf "%s\t%s\n", $_-$DistFlank, $CumulWgt/$fnormal;
    printf "%s\t%s\n", $_+$DistFlank, ($CumulWgt+=$DistribDens{$_}*$climax)/$fnormal;
  }
}


# main calculations for a sample
#
# INTERFACE
# - argument 1: path of file containing array of values
#
# - global options:
#   -debug      print debug protocol to STDERR
#   -OutImg     [STD]
#   -OutImgWidth  [STD]
#   -SampleWgt  sample data is weighted, means: 2nd data column specifies
#               weight for each data entry
#   --climax    direction of cumulation
#   --smoothws  window size for smoothening
#
sub ProgSampleMeas {

  # function parameters
  my ($PathTab) = @_;
  my $debug = $ProgOpt{-debug};
  my $dbg2  = $debug ? $debug-1 : undef;
  my $bWgt = $ProgOpt{-SampleWgt};
  my %path; $path{stamp} = &PrepOstump();
  my %img;
  $img{flag} = int (defined($ProgOpt{-OutImg}) or $ProgOpt{-OutImgWidth}
    or $ProgOpt{-OutImgType});
  $img{width} = $ProgOpt{-OutImgWidth} || $ProgParam{default}{OutImgWidth};
  $img{StepNum} = $img{width} * 0.5;
  $img{gauss} = $ProgOpt{-var}{ImgGauss} ? 1 : 0;

  #############################################################################
  # calculate statistics, create/write report
  my %result;
  my ($pPlot,$xLast,$xCurr);

  # main data
  $result{time} = &TimeStr();
  $result{SrcFile} = &PathExpand ($PathTab, -pipe=>'<');
  my $DataTabType = $bWgt ? 'AA' : 'A1';
  $result{data} = &PlainToTable ($PathTab, -TabType=>$DataTabType,
    -DelimCol=>$ProgOpt{-DelimCol}, -comments=>1, -debug=>$dbg2);
  my $ValN = int @{$result{data}};
  unless ($ValN) {
    die sprintf "ERROR: don't see any data in file %s\n", $PathTab||"''";
  }
  # is data really weighted?
  if ($bWgt and int (grep{ $$_[1] }@{$result{data}}) != int(@{$result{data}})) {
    printf STDERR "WARNING: some data lines contain zero weights, turning to unweighted mode\n";
    $result{data} = &TableConvert ('AA', 'A1', $result{data}, -debug=>$dbg2);
    undef $bWgt;
  }

  # sample metrics
  my @DataSorted;
  if ($bWgt) {
    @DataSorted = sort { $$a[0] <=> $$b[0]; } @{$result{data}};
    @{$result{DataRange}}{'off','end'} = ($DataSorted[0][0], $DataSorted[-1][0]);
    $result{metrics} = &SampleMetricsWgt ($result{data}, -median=>1, -debug=>$dbg2);
  } else {
    @DataSorted = sort { $a <=> $b; } @{$result{data}};
    @{$result{DataRange}}{'off','end'} = ($DataSorted[0], $DataSorted[-1]);
    $result{metrics} = &SampleMetrics ($result{data}, -median=>1, -debug=>$dbg2);
  }
  $result{metrics}{min} = $result{DataRange}{off};
  $result{metrics}{max} = $result{DataRange}{end};

  # empirical case density distribution ("atomic resolution histogram")
  $result{distribution} = &DistribEmpir (\@DataSorted,
    exists($ProgOpt{-var}{xmin})?(-RangeMin=>$ProgOpt{-var}{xmin}):(),
    #exists($ProgOpt{-var}{xmin})?(-RangeMin=>$ProgOpt{-var}{xmax}):(),
    );
  if (!$result{distribution} or !int(@{$result{distribution}})) {
    print  STDERR "ERROR in \&DistribEmpir\n";
    &DistribEmpir (\@DataSorted, -debug=>1,
      exists($ProgOpt{-var}{xmin})?(-RangeMin=>$ProgOpt{-var}{xmin}):());
    exit 1;
  }
  # eventually smoothen empirical distribution function
  # - $result{distribution}[-1][0], $result{distribution}[0][0] are not equal to
  #   $result{DataRange}{off}, $result{DataRange}{end}
  $img{StepSize} = ($result{distribution}[-1][0] - $result{distribution}[0][0]) / $img{StepNum};
  if ($ProgOpt{-var}{smoothws}) {
    $pPlot = Math::Plot2D->new($result{distribution});
    $result{distribution} = $pPlot->SmoothPlot ($img{StepSize},
      -window=>$ProgOpt{-var}{smoothws}, -debug=>$dbg2);
    if ($ProgOpt{-var}{xmin}) {
      $pPlot = Math::Plot2D->new($result{distribution});
      $result{distribution} = $pPlot->FlipWing ($ProgOpt{-var}{xmin}, -1,
        -debug=>$dbg2);
    }
#    if ($ProgOpt{-var}{xmax}) {
#      $pPlot = Math::Plot2D->new($result{distribution});
#      $result{distribution} = $pPlot->FlipWing ($ProgOpt{-var}{xmax}, +1,
#        -debug=>$dbg2);
#    }
  }

  # Gaussian approximation
  if ($img{flag} and $img{gauss}) {
    $xCurr = $result{distribution}[0][0];
    while ($xCurr <= $result{distribution}[-1][0]) {
      push @{$result{DistribGauss}},
        [ $xCurr, &GaussVal($result{metrics}{mean},$result{metrics}{s},$xCurr) ];
      $xCurr += $img{StepSize};
    }
  }

  # cumulative frequency function
  $pPlot = Math::Plot2D->new($result{distribution});
  $xCurr = $xLast = $pPlot->Xmin();
  my $integr_dir = (($ProgOpt{-var}{climax}||1)>0) ? +1 : -1;
  my $integral = (($ProgOpt{-var}{climax}||1)>0) ? 0 : 1;
  while ($xCurr <= $pPlot->Xmax()) {
    $integral += $integr_dir * $pPlot->Integral($xLast, $xCurr);
    push @{$result{DistribCumul}}, [ $xCurr, $integral ];
    $xLast = $xCurr;
    $xCurr += $img{StepSize};
  }

  # output result data, verbose &GetWriteHandle
  $path{data} = $path{stamp} .'_sample.dat';
  if (my $hOutDat = &GetWriteHandle($path{data})) {
    &DataPrint (\%result, -handle=>$hOutDat, -NoAddr=>1);
  }

  #############################################################################
  # graphical presentation

  # prepare graphs
  if ($img{flag}) {
    my $hOutDat;

    # graph of empirical distribution function
    my %graph = (
      BgTranspar => $ProgOpt{-OutImgTransp},
      plot => [
        { DimPixel     => { x=>$img{width} },
          HeightRel    => $ProgOpt{-OutImgRelhigh} || $ProgParam{default}{OutImgRelhigh},
          DataType     => 'AA',
          data         => $result{distribution},
          DataRange    => {
              exists($ProgOpt{-var}{xmin}) ? (x=>[$ProgOpt{-var}{xmin}]) : (),
              y => [ 0 ]
            },
          ReprType     => 'line',
        },
        $img{gauss} ?
        { DataType     => 'AA',
          data         => $result{DistribGauss},
          DataRangeRef => 0,
          ReprType     => 'line',
          ReprColor    => 'blue',
        } : (),
        ],
      scale => [
        { location => 'x',
          PlotNum  => 0,
        },
        { location => 'left',
          PlotNum  => 0,
        },
        { location => 'right',
          PlotNum  => 0,
        },
        ],
      );

    # save image and plot graph data, verbose &GetWriteHandle
    $path{data} = $path{stamp} .'_distrib_img.dat';
    if ($hOutDat = &GetWriteHandle($path{data})) {
      &DataPrint (\%graph, -handle=>$hOutDat, -NoAddr=>1);
    }
    $path{img} = $path{stamp} .'_distrib.png';
    if (&Graph (\%graph, -save=>$path{img}, -debug=>$dbg2)) {
      printf "writing file %s\n", $path{img}||"''";
    } else {
      printf STDERR "%s. ERROR: failure in plotting distribution graph to %s\n", &MySub, $path{img}||"''";
    }

    # graph of cumulative frequency function
    %graph = (
      BgTranspar => $ProgOpt{-OutImgTransp},
      plot => [
        { DimPixel  => { x=>$img{width} },
          HeightRel => 0.7,
          DataType  => 'AA',
          data      => $result{DistribCumul},
          DataRange => {
              exists($ProgOpt{-var}{xmin}) ? (x=>[$ProgOpt{-var}{xmin}]) : (),
              y => [ 0, 1 ],
            },
          ReprType  => 'line',
        },
        ],
      scale => [
        { location => 'x',
          PlotNum  => 0,
        },
        { location => 'left',
          PlotNum  => 0,
        },
        { location => 'right',
          PlotNum  => 0,
        },
        ],
      );

    # save image and plot graph data
    $path{data} = $path{stamp} .'_cumfreq_img.dat';
    $hOutDat = FileHandle->new($path{data},'w');
    &DataPrint (\%graph, -handle=>$hOutDat, -NoAddr=>1);
    printf "writing file %s\n", $path{data}||"''";
    $path{img} = $path{stamp} .'_cumfreq.png';
    if (&Graph (\%graph, -save=>$path{img}, -debug=>$dbg2)) {
      printf "writing file %s\n", $path{img}||"''";
    } else {
      printf STDERR "%s. ERROR: failure in plotting cumulative frequency graph to %s\n", &MySub, $path{img}||"''";
    }
  }
}


# linear regression
#
# INTERFACE
# - argument 1: path of file containing array of value pairs
#
# - global options:
#   -debug      print debug protocol to STDERR
#   -OutImg     [STD]
#   -OutImgWidth  [STD]
#
sub ProgRegressLin {
  my ($PathTab) = @_;

  # function parameters
  my $debug = $ProgOpt{-debug};
  my $dbg2  = $debug ? $debug-1 : undef;
  my %path; $path{stamp} = &PrepOstump();
  my %img; $img{flag} = int (defined($ProgOpt{-OutImg}) or $ProgOpt{-OutImgWidth}
    or $ProgOpt{-OutImgType});

  #############################################################################
  # mathematics part
  my %result;

  # load data
  $result{time} = &TimeStr();
  $result{SrcFile} = &PathExpand ($PathTab, -pipe=>'<');
  $result{data} = &PlainToTable ($PathTab, -TabType=>'AA',
    -DelimCol=>$ProgOpt{-DelimCol}, -comments=>1, -debug=>$dbg2);
  unless (int @{$result{data}}) {
    die sprintf "ERROR: don't see any data in file %s\n", $PathTab||"''";
  }

  # do regressions
  $result{lin0} = &RegressLin0 ($result{data});
  $result{lin} = &RegressLin ($result{data});

  # residual deviation squares
  foreach my $model (qw(lin0 lin)) {
    my $devsqr=0;
    foreach (@{$result{data}}) {
      my $e = $result{$model}{ay} * $_->[0] + $result{$model}{by};
      $devsqr += ($_->[1]-$e) **2;
    }
    $result{$model}{resid_vary} = $devsqr / ($result{lin0}{n}-1);
    $devsqr=0;
    foreach (@{$result{data}}) {
      my $e = ($_->[1]-$result{$model}{bx}) / $result{$model}{ax};
      $devsqr += ($_->[0]-$e) **2;
    }
    $result{$model}{resid_varx} = $devsqr / ($result{lin0}{n}-1);
  }

  # output result data, verbose &GetWriteHandle
  $path{data} = $path{stamp} .'_RegressLin.dat';
  if (my $hOutDat = &GetWriteHandle($path{data})) {
    &DataPrint ({ constraint_lin0=>$result{lin0}, constraint_lin=>$result{lin} },
      -handle=>$hOutDat, -NoAddr=>1);
  }

  #############################################################################
  # graphical presentation

  # prepare graphs
  if ($img{flag}) {
    $img{varaxis}=$ProgOpt{-var}{varaxis}||'y';

    # graph of empirical distribution function
    my %graph = (
      BgTranspar => $ProgOpt{-OutImgTransp},
      plot => [
          { DimPixel  => { x=>$img{width} },
            HeightRel => 0.7,
            DataType  => 'AA',
            data      => $result{data},
            #DataRange => { y=>[0] },
            ReprType  => 'pixel',
            ReprSizePixel => 3,
          },
  #       we need:
  #       - regression lines for both regression subtypes (0, normal)
  #       - coordinates of the cross points of the displayed data range
          (!exists($ProgOpt{-var}{origin}) or $ProgOpt{-var}{origin}==0) ?
          { DataType => 'AA',
            data     => [
                 [&Min(map{$_->[0]}@{$result{data}}),&Min(map{$_->[0]}@{$result{data}})*$result{lin}{'a'.$img{varaxis}}+$result{lin}{'b'.$img{varaxis}}],
                 [&Max(map{$_->[0]}@{$result{data}}),&Max(map{$_->[0]}@{$result{data}})*$result{lin}{'a'.$img{varaxis}}+$result{lin}{'b'.$img{varaxis}}],
               ],
            DataRangeRef => 0,
            ReprType  => 'line',
            ReprColor => 'blue',
          } : (),
          (!exists($ProgOpt{-var}{origin}) or $ProgOpt{-var}{origin}==1) ?
          { DataType => 'AA',
            data     => [
                 [&Min(map{$_->[0]}@{$result{data}}),&Min(map{$_->[0]}@{$result{data}})*$result{lin0}{'a'.$img{varaxis}}+$result{lin0}{'b'.$img{varaxis}}],
                 [&Max(map{$_->[0]}@{$result{data}}),&Max(map{$_->[0]}@{$result{data}})*$result{lin0}{'a'.$img{varaxis}}+$result{lin0}{'b'.$img{varaxis}}],
               ],
            DataRangeRef => 0,
            ReprType  => 'line',
            ReprColor => 'green',
          } : (),
        ],
      scale => [
          { location => 'x',
            PlotNum  => 0,
          },
          { location => 'left',
            PlotNum  => 0,
          },
          { location => 'right',
            PlotNum  => 0,
          },
        ],
      );

    # save image and plot graph data
    $path{data} = $path{stamp} .'_distrib_img.dat';
    if (my $hOutDat = &GetWriteHandle($path{data})) {
      printf "writing file %s\n", $path{data}||"''";
      &DataPrint (\%graph, -handle=>$hOutDat, -NoAddr=>1);
    }
    $path{img} = $path{stamp} .'_distrib.png';
    if (&Graph (\%graph, -save=>$path{img}, -debug=>$dbg2)) {
      printf "writing file %s\n", $path{img}||"''";
    } else {
      printf STDERR "%s. ERROR: failure in plotting distribution graph to %s\n", &MySub, $path{img}||"''";
    }
  }
}


################################################################################
# statistical tests
################################################################################


# one-way ANOVA test
#
# INTERFACE
# - argument 1: table file
#
# - global options:
#   -debug      print debug protocol to STDERR
#
# DESCRIPTION
# - variables for empirical measures used as in
#   http://helios.bto.ed.ac.uk/bto/statistics/tress6.html
#
sub ProgAnova {
  my ($PathTab) = @_;
  my $debug = $ProgOpt{-debug};
  my $dbg2  = $debug ? $debug-1 : undef;
  require Statistics::Distributions;

  # load data
  my $pVal = &PlainToTable ($PathTab, -TabType=>'AA',
    -DelimCol=>$ProgOpt{-DelimCol}, -comments=>1, -debug=>$dbg2);
  if (int(@$pVal) < 2) {
    die sprintf "ERROR: less than two data lines in file %s\n", $PathTab||"''";
  }

  # derive line sums (@lsum), sum of column A ($asum), and total sum ($n)
  my (@v,@eval,$bErrComma);
  foreach my $pLn (@$pVal) {
    if (int(@$pLn) < 2) {
      die sprintf "ERROR: data line %d contains less than two values\n", int(@eval)+1;
    }
    if (!$bErrComma and grep{ m/,/ } @$pLn) {
      warn sprintf "WARNING: data line %d contains comma(s)\n", int(@eval)+1;
      $bErrComma=1;
    }
    push @v, @$pLn;
    my $met = &SampleMetrics($pLn);
    $met->{A} = &Sum(map{ $_**2 } @$pLn);
    $met->{B} = &Sum(@$pLn)**2 / int(@$pLn);
    push @eval, $met;
  }

  # report values
  my $ni_max = &Max(map{ int(@$_) } @$pVal);
  for (my $i=0; $i<$ni_max; ++$i) {
    printf "replic%d", $i+1;
    foreach my $pLn (@$pVal) {
      printf "\t%s", $pLn->[$i];
    }
    print  "\n";
  }
  # report calculi
  foreach my $calc (qw(n A B DevSq var)) {
    printf "$calc";
    for (my $i=0; $i<int(@$pVal); ++$i) {
      printf "\t%s", $eval[$i]->{$calc};
    }
    printf "\t%s\n", &Sum(map{ $_->{$calc} } @eval);
  }
  print  "\n";
  # sum of squares, square of total by n
  my $A = &Sum(map{ $_->{A} } @eval);
  my $B = &Sum(map{ $_->{B} } @eval);
  # grand total squared, divided by n
  my $D = &Sum(@v)**2 / int(@v);
  # degrees of freedom
  my $df_inter = int(@$pVal)-1;
  my $df_intra = &Sum(map{ $_->{n}-1 } @eval);
  print  "D\t$D\n";
  printf "df_inter\t%s\n", $df_inter;
  printf "df_intra\t%s\n", $df_intra;

  # degrees of freedom
  print  "\n";
  my $Var_inter = ($B-$D)/$df_inter;
  my $Var_intra = ($A-$B)/$df_intra;
  printf "Var_inter=(B-D)/df_inter\t%s\n", $Var_inter;
  printf "Var_intra=(A-B)/df_intra\t%s\n", $Var_intra;
  my $fprob = &Statistics::Distributions::fprob($df_inter,$df_intra,$Var_inter/$Var_intra);
  printf "F,P\t%s\t%.2e%s\n", $Var_inter/$Var_intra,
    $fprob, ($fprob>0.01)? sprintf(' = %.3f',$fprob):'';
}


# chi square for four-field table
#
# INTERFACE
# - argument 1-4: table values
#
# - global options:
#   -debug      print debug protocol to STDERR
#
sub ProgChisqr4 {
  my (@val) = @_;
  if (int(@val)<4) {
    warn sprintf "WARNING: few arguments (%d), 4 expected\n", int(@val);
  }
  require Statistics::Distributions;

  # derive line sums (@lsum[1..2]) and total sum ($n)
  my @lsum = ($val[0]+$val[1], $val[2]+$val[3]);   ### well, not further used
  my $n = &Sum (@lsum);

  # calculate chi square
  my $chisqr = $n * ($val[0]*$val[3] - $val[1]*$val[2])**2
    / ($val[0]+$val[1]) / ($val[2]+$val[3])
    / ($val[0]+$val[2]) / ($val[1]+$val[3])
    ;

  # p-value, report
  my $chisprob = &Statistics::Distributions::chisqrprob(1,$chisqr);
  printf "chi^2,P\t%s\t%.2e%s\n", $chisqr,
    $chisprob, ($chisprob>0.01)? sprintf(' = %.3f',$chisprob):'';
}


# chi square for k*2-field table
#
# INTERFACE
# - argument 1: table file, two values per line
#
# - global options:
#   -debug      print debug protocol to STDERR
#
sub ProgChisqrk2 {
  my ($PathTab) = @_;
  my $debug = $ProgOpt{-debug};
  my $dbg2  = $debug ? $debug-1 : undef;
  require Statistics::Distributions;

  # load data
  my $pVal = &PlainToTable ($PathTab, -TabType=>'AA',
    -DelimCol=>$ProgOpt{-DelimCol}, -comments=>1, -debug=>$dbg2);
  if (int(@$pVal) < 2) {
    die sprintf "ERROR: less than two data lines in file %s\n", $PathTab||"''";
  }

  # derive line sums (@lsum), sum of column A ($asum), and total sum ($n)
  my (@lsum,$asum);
  foreach my $pLn (@$pVal) {
    if (int(@$pLn) < 2) {
      die sprintf "ERROR: data line %d contains less than two values\n", int(@lsum)+1;
    }
    push @lsum, $pLn->[0]+$pLn->[1];
    $asum += $pLn->[0];
  }
  my $n = &Sum (@lsum);

  # calculate chi square
  my $chisqr = $n**2 / $asum / ($n-$asum);
  my $m=0;
  for (my $i=0; $i<int(@lsum); ++$i) {
    $m += $pVal->[$i][0]**2 / $lsum[$i];
  }
  $m -= $asum**2 / $n;
  $chisqr *= $m;

  # p-value, report
  my $chisprob = &Statistics::Distributions::chisqrprob(int(@$pVal)-1,$chisqr);
  printf "chi^2,f,P\t%s\t%d\t%.2e%s\n", $chisqr, int(@$pVal)-1,
    $chisprob, ($chisprob>0.01)? sprintf(' = %.3f',$chisprob):'';
}


# significance of item order in repeated sampling (by permutation test)
#
# INTERFACE
# - argument 1:  table file with list of items
# - argument 2+: table file with list of items
#
sub ProgPermutOrder {
  my (@PathTab) = @_;
  my $alpha = $ProgOpt{-var}{alpha} || $ProgParam{default}{PermutOrderAlpha};

  # load data
  my @val;
  foreach my $f (@PathTab) {
    push @val, my $pTab = [];
    my $h = FileHandle->new($f) or die "ERROR: unable to read file $f\n";
    while (<$h>) {
      if (m/^./) {
        chomp;
        push @$pTab, $_;
      }
    }
    printf "# file %s with %d items\n", $f, int(@$pTab);
  }
  if (!int(@val)) { return }

  # empirical position distributions
  my $pInfoIs = &_PermutOrderStat(\@val);

  # position distributions from permutations
  my $nperm = $ProgOpt{-var}{permut} || $ProgParam{default}{PermutN};
  printf "# number of permutations %d\n", $nperm;
  my %InfoCumul;
  for (my $j=0; $j<=$nperm; ++$j) {
    my $pPermut = &PermutOrderPermut(\@val);
  #&DataPrint ($pPermut); exit;
    my $pInfo = &_PermutOrderStat($pPermut);
  #&DataPrint ($pInfo); exit;
    foreach my $i (keys %$pInfoIs) {
      foreach my $v (qw(mean s DevSq)) {
        push @{$InfoCumul{$i}{$v}||=[]}, $pInfo->{$i}{$v};
      }
    }
  }

  # report distributions
  my $sAlpha = sprintf('%.2f',$alpha*100);
  $sAlpha =~ s/\./p/;
  printf "#\n#item\tn_pos\tmean_pos\tsd_pos\tdevsqr_pos\tcalc_mean\tcalc_mean_sd\tcalc_mean_trusted%sperc\tcalc_sd\tcalc_sd_sd\tcalc_sd_trusted%sperc\tcalc_devsqr\tcalc_devsqr_sd\n",
    $sAlpha, $sAlpha;
  foreach my $i (sort{ $pInfoIs->{$a}{mean}<=>$pInfoIs->{$b}{mean} } keys %$pInfoIs) {
    my $pMeanMetrics = &SampleMetrics($InfoCumul{$i}{mean});
    $pMeanMetrics->{alpha} = do{
      my @v=sort{ $a<=>$b } @{$InfoCumul{$i}{mean}};
      $v[$alpha*($nperm-1)];
      };
    $pMeanMetrics->{'95'} = do{
      my @v=sort{ $a<=>$b } @{$InfoCumul{$i}{mean}};
      $v[0.95*($nperm-1)];
      };
    my $pSdMetrics = &SampleMetrics($InfoCumul{$i}{s});
    $pSdMetrics->{alpha} = do{
      my @v=sort{ $a<=>$b } @{$InfoCumul{$i}{s}};
      $v[$alpha*($nperm-1)];
      };
    $pSdMetrics->{'95'} = do{
      my @v=sort{ $a<=>$b } @{$InfoCumul{$i}{s}};
      $v[0.95*($nperm-1)];
      };
    my $pDevsqrMetrics = &SampleMetrics($InfoCumul{$i}{DevSq});
    my $accurorder = (int(@{$val[0]})<100)? 3 : 2+log(int(@{$val[0]}))/log(10);
    printf "%s\t%d\t%.*f\t%.*f\t%.2f\t%.3f\t%.3f\t%.3f-%.3f\t%.3f\t%.3f\t%.3f-%.3f\t%.2f\n", $i,
      $pInfoIs->{$i}{n}, $accurorder, $pInfoIs->{$i}{mean}, $accurorder, $pInfoIs->{$i}{s}, $pInfoIs->{$i}{DevSq},
      $pMeanMetrics->{mean}, $pMeanMetrics->{s}, $pMeanMetrics->{alpha}, $pMeanMetrics->{'95'},
      $pSdMetrics->{mean}, $pSdMetrics->{s}, $pSdMetrics->{alpha}, $pSdMetrics->{'95'},
      $pDevsqrMetrics->{mean};
  }
}

# evaluation of item order in samples
#
# INTERFACE
# - argument 1: reference to 2D array of table->item
# - return val: reference to statistics info (hash of info hashes)
#
sub _PermutOrderStat {
  my ($pData) = @_;

  # sample relative positions
  my %aPos;
  foreach my $pTab (@$pData) {
    my $sz=$#$pTab;
    for (my $i=0; $i<=$sz; ++$i) {
      push @{$aPos{$pTab->[$i]}||=[]}, $i/$sz;
    }
  }

  # evaluate distributions of positions
  my %dic;
  foreach my $i (keys %aPos) {
    $dic{$i} = &SampleMetrics($aPos{$i});
  }

  return \%dic;
}

# randomize item order in samples
#
# INTERFACE
# - argument 1: reference to 2D array of table->item
# - return val: reference to randomized 2D array of table->item
#
sub PermutOrderPermut {
  my ($pData) = @_;
  my @val;
  foreach my $pTab (@$pData) {
    push @val, &RandArrayOrder($pTab);
  }
  return \@val;
}


# significance of item order in repeated sampling (by permutation test)
#
# INTERFACE
# - argument 1+: table file with list of items
#
sub ProgPermutCrowd {
  my (@PathTab) = @_;

  # load data
  my @val;
  foreach my $f (@PathTab) {
    my $n=0;
    my $h = FileHandle->new($f) or die "ERROR: unable to read file $f\n";
    while (<$h>) {
      if (m/^./) {
        chomp;
        push @val,$_;
        ++$n;
      }
    }
    printf "# file %s with %d items\n", $f, $n;
  }

  # empirical position distributions
  my $pInfoIs = &_PermutCrowdStat(\@val);
  my $fIsScore = sqrt(
    &Sum(map{ $_->{devsqr} } values %$pInfoIs) / (int(@val)-int(keys %$pInfoIs))
    ) * 1.00001;

  # position distributions from permutations
  my $nperm = $ProgOpt{-var}{permut} || $ProgParam{default}{PermutN};
  printf "# real skew score %f\n", $fIsScore;
  printf "# number of permutations %d\n", $nperm;
  my $nout=0;
  for (my $j=0; $j<=$nperm; ++$j) {
    my $pPermut = &RandArrayOrder(\@val);
    my $pInfo = &_PermutCrowdStat($pPermut);
    my $dssum = &Sum(map{ $_->{devsqr} } values %$pInfo);
    my $score = sqrt($dssum/(int(@val)-int(keys %$pInfoIs)));
    if ($score<$fIsScore) { ++$nout }
  }
  printf "# quantile probability %f\n", $nout/$nperm;
}

# evaluation of item order in samples
#
# INTERFACE
# - argument 1: reference to 2D array of table->item
# - return val: reference to statistics info (hash of info hashes, per item)
#
sub _PermutCrowdStat {
  my ($pData) = @_;

  # sample relative positions
  my %aPos;
  my $sz=$#$pData;
  for (my $i=0; $i<=$sz; ++$i) {
    push @{$aPos{$pData->[$i]}||=[]}, $i/$sz;
  }

  # evaluate distributions of positions
  my %dic;
  foreach my $k (keys %aPos) {
    $dic{$k} = {
      n => int(@{$aPos{$k}}),
      mean => my $m=&Mean(@{$aPos{$k}}),
      };
    $dic{$k}{devsqr} = &Sum(map{ ($_-$m)**2 } @{$aPos{$k}}),
  }

  return \%dic;
}


################################################################################
# likelihood distributions
################################################################################


# calculate binomial table
#
# INTERFACE
# - argument 1:  range of value n
# - argument 2:  range of value p
# - argument 3?: range for x values, e.g. '..150' meaning:
#                range from default to 150.
#                This takes effect inside this code.
#
# - global options:
#   -debug       [STD]
#
sub ProgBinom {
  my ($ArgN,$ArgP,$ArgX) = @_;
  my $debug = $ProgOpt{-debug};
  my $dbg2  = $debug ? $debug-1 : undef;
  my $accur = $ProgOpt{-accur} || 5;

  # preliminarily resolve value range of x
  my %val;
  if ($ArgX =~ m/^(\d*)\.\.(\d*)$/) {
    $val{x} = { min=>$1, max=>$2 };
  } else {
    $val{x} = { min=>"0", max=>'' };
  }
  # resolve value range of n
  if ($ArgN =~ m/^\d+$/) { @{$val{n}} = ($ArgN); }
  elsif ($ArgN=~m/^(\d+)\.\.(\d+)(:(\d+))?$/ and $3||=1 and $3>0) {
    for (my $i=$1+0; $i<=$2; $i+=$3) {
      push @{$val{n}}, $i;
    }
  }
  unless (@{$val{n}}) {
    printf STDERR "%s. ERROR: unable to resolve range argument for n: %s\n", &MySub, $ArgN;
    return undef;
  }
  # resolve value range of p
  if ($ArgP =~ m/^[\d.]+([eE][+-]?[\d]+)?$/) { @{$val{p}} = ($ArgP); }
  elsif ($ArgP =~ m/^([\d.]+([eE][+-]?[\d]+)?)\.\.([\d.]+([eE][+-]?[\d]+)?):([\d.]+([eE][+-]?[\d]+)?)$/
  and $5 > 0) {
    for (my $i=$1+0; $i<=($3+0.1*$5); $i+=$5) {
      push @{$val{p}}, $i;
    }
  }
  unless (@{$val{p}}) {
    printf STDERR "%s. ERROR: unable to resolve range argument p: %s\n", &MySub, $ArgP;
    return undef;
  }
  if (grep { $_<=0 or $_>=1 } @{$val{p}}) {
    printf STDERR "%s. ERROR: bad range argument p: %s%s\n", &MySub, $ArgP,
      $debug ? (' => '.join(' ',@{$val{p}}) ) : '';
    return undef;
  }
  # debug arguments: p, n, x
  if ($debug) {
    printf STDERR "%s. p range: %s\n", &MySub, join(',',@{$val{p}});
    printf STDERR "%s. n range: %s\n", &MySub, join(',',@{$val{n}});
    printf STDERR "%s. x range: %s..%s\n", &MySub, $val{x}{min}||"0", $val{x}{max}||"n";
  }

  # output table header
  printf "# %s %s: binomial likelihood distributions\n", $ProgFile, $ProgMode;
  printf "# date/time: %s\n", &TimeStr();
  print  "# range of n: $ArgN\n";
  print  "# range of p: $ArgP\n";
  print  "#\n";

  # loop over values of n
  my $pTable;
  foreach my $ItN (@{$val{n}}) {

    # determine effective range of X
    $val{x}{min_eff} = length($val{x}{min}) ? int($val{x}{min}) : 0;
    $val{x}{max_eff} = length($val{x}{min}) ? int($val{x}{max}) : $ItN;

    # calculate binomial distribution
    foreach my $ItP (@{$val{p}}) {
      my ($pTableProb,$pTablePsi) = &BinomRow ($ItN, $ItP,
        length($val{x}{min}) ? (-ForceMinX=>$val{x}{min}) : (),
        length($val{x}{max}) ? (-ForceMaxX=>$val{x}{max}) : (),
        -accur=>$ProgOpt{-accur}, -debug=>$dbg2);
      $$pTable{probab}{$ItP} = $pTableProb;
      $$pTable{psi}{$ItP} = $pTablePsi;
    }

    # output table of case probabilities
    print  "# table of P values\n";
    printf "# n: %d\n", $ItN;
    print  "# range of p: $ArgP\n";
    print  "#\n";
    print  "# column labels:\n";
    printf "# X\t%s\n", join ("\t", map { "p=$_" } @{$val{p}});
    for (my $ItX=$val{x}{min_eff}; $ItX<=$val{x}{max_eff}; ++$ItX) {
      printf "%d", $ItX;
      foreach my $ItP (@{$val{p}}) {
        printf "\t%.*f", $accur, $$pTable{probab}{$ItP}[$ItX];
      }
      print "\n";
    }

    # output table of cumulative case probabilities
    print  "# table of Psi values\n";
    printf "# n: %d\n", $ItN;
    print  "# range of p: $ArgP\n";
    print  "#\n";
    print  "# column labels:\n";
    printf "# X\t%s\n", join ("\t", map { "p=$_" } @{$val{p}});
    for (my $ItX=$val{x}{min_eff}; $ItX<=$val{x}{max_eff}; ++$ItX) {
      printf "%d", $ItX;
      foreach my $ItP (@{$val{p}}) {
        printf "\t%.*f", $accur, $$pTable{psi}{$ItP}[$ItX];
      }
      print "\n";
    }
  }  # end loop over values of n
}


################################################################################
# random values, simulations
################################################################################


# random values following normal-distribution
#
# INTERFACE
# - argument 1:  my value
# - argument 2:  sigma value
# - argument 3?: sample size, default: 1
#
# - global options:
#   -debug       [STD]
#
sub ProgRandGauss {
  my ($ArgMy,$ArgSigma,$ArgN) = @_;
  my $debug = $ProgOpt{-debug};
  my $dbg2  = $debug ? $debug-1 : undef;
  $ArgN ||= 1;

  # loop over values of n
  for (my $i=0; $i<$ArgN; ++$i) {

    # need two components:
    # - sign
    # - homogeneously random value in range [0,1[
    my $ValRand = (rand() - 0.5) * 2;

    # evaluate to position in Gaussian distribution
    my $ValTimesS = &GaussLimit (abs($ValRand), -debug=>$dbg2);

    # output
    printf "%s\n", $ArgMy + &Sign($ValRand) * $ValTimesS * $ArgSigma;
  }
}
# $Id: Statist.pl,v 1.58 2018/06/05 18:02:56 szafrans Exp $
