#!/usr/bin/env perl

#  Simple binary search for computing x = f^{-1}(y) such that
#  |y - f(x)| < desired_precision

#  Take as input an interval [l, u] and a target quantile, q
#  Precondition F(l) < q and F(u) > q
#  Do a binary search for the value F(m) such that |F(m) - q| < desired_precision
#  Output the final F(m) value as well as all the intermediate (v,F(v)) pairs

use strict;
use warnings;
use Getopt::Long;
use locale;
use File::Basename;
use File::Path;
use File::Find;
use FileHandle;
STDOUT->autoflush(1);

#############     Option Defaults    #########

my $maxIters  = 25;
my $lower     = "";
my $upper     = "";
my $quantile  = "";
my $precision = 0.0000001;
my $execStr   = "./gs -m 50 -n 60 -p 400";

my $defaultIntervalLen = 0.01;

my $prependQuantile = 1;

my $verbose = 0;
my $commentSymbol="%";
my $help=0;

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

#############         Usage          #########
my $baseName=basename($0);
sub usage{
   print STDERR (
   "Description: \n" .
   "Usage: $baseName [options]\n" .
   "   -help   Print this usage message.\n" .
   "\n"
   );
}

my $originalArgsString="@ARGV";
############      Read Options       #########
my $validOpts=GetOptions(
  "lower=f"     =>\$lower,
  "upper=f"     =>\$upper,
  "quantile=f"  =>\$quantile,
  "precision=f" =>\$precision,
  "execstr=s"   =>\$execStr,
  "prependq!"   =>\$prependQuantile,
  "defint=f"    =>\$defaultIntervalLen,
  "maxiters=i"  =>\$maxIters,
  "verbose!"    =>\$verbose,
  "help!"       =>\$help
);

##############################################
if($help) {
  usage(); exit(0);
}

# sub numerically { $a <=> $b; }

sub error{
    my $errString=shift @_;
    print STDERR "$errString\n";
    usage(); exit -1;
}


error("") if(!$validOpts);

if($verbose) {
    print"
lower     = $lower
upper     = $upper
quantile  = $quantile
precision = $precision
execStr   = $execStr
";
}

error("Need at least lower bound for quantile estimation") if $lower eq "";

my $l = $lower;
my $u = $upper;
$u = $l + $defaultIntervalLen if $upper eq "";

my $diff = 2 * $precision; # arbitrary starting value as long as > precision

sub F {
    my ($v) = @_;
    my $Fv = qx{$execStr -x $v};
    return $Fv;
}

my $mi = $maxIters;
my $final = 0;
while($mi > 0) {
    # calculate F(m), m midpoint of [l,u] range
    my $m  = ($l + $u)/2.0;
    my $Fm = F($m);
    chomp $Fm;
    $diff = abs($Fm - $quantile);
    print "$quantile\t" if $prependQuantile;
    print "$m\t$Fm";
    my $final = 1 if $diff < $precision;
    print "\tFINAL" if $final;
    print "\n";
    last if $final;
    if($Fm > $quantile) {
        $u = $m;
    } else {
        $l = $m
    }
    --$mi;
}