#!/usr/bin/perl


## Script written by Elisenda Feliu, January 2010 ##

## This script computes the P-value of finding at least l near-native decoys when selecting p decoys from a pool of d decoys containing n near-native ones. ##

use strict;
use warnings;

use Getopt::Long;

select(STDERR);

my ($help,$ndec,$n0,$pick,$atleast);

&GetOptions("d=s"   => \$ndec,
             "n=s" => \$n0,
             "p=s" => \$pick,
	     "l=s" => \$atleast,
             "h"     => \$help       );

unless (defined $ndec && defined $n0 && defined $pick && not defined $help){
   printf "** This script computes the P-value of finding at least l near-native decoys when selecting p decoys from a pool of d decoys containing n near-native ones.**\nUsage:\n\t./compute_prob.pl -d [Total number of decoys] -n [Number of near-native decoys in the set] -p [Number of decoys to pick] -l [Minimum number of hits]\n\nOption specification is required for all options but -l which defaults to 1.
          \n";
   exit;
}

if(not defined $atleast){ $atleast=1;}
my ($k,$prob,$min);


# Compute the probability of at least one hit at p predictions, for this type of family (depends on n0)
if($n0!=0){ 
	$min=&min($n0,$pick);
	$prob=0;
	foreach $k ($atleast..$min){   # prob of l guesses among p predictions. l must be lower than n0 and p
		$prob += &probability($n0,$ndec-$n0,$k,$pick);
	}
}

print  "The probability of getting at least $atleast hit taking $pick decoys from a set of $ndec from which there are $n0 near-native decoys is:\n$prob \n";

sub probability {
	my ($n0,$n1,$k,$m)=@_;
	my ($prob,$j);
	if($k>$m) {$prob=0; return($prob); }
	if($k>$n0){ $prob=0; return($prob);}
	if($k<$m-$n1){$prob=0; return($prob);}
	$prob=1;
	for($j=$n0;$j>=$n0-$k+1;$j=$j-1){
		$prob=$prob*$j;	
	}
	for($j=$n0+$n1;$j>=$n1+$k+1;$j=$j-1){
		$prob=$prob*($j-$m)/$j;	
	}
	for($j=$n1+$k;$j>=$n1+1;$j=$j-1){
		$prob=$prob/$j;
	}

 	for($j=$m;$j>=$m-$k+1;$j=$j-1){
 		$prob=$prob*$j;	
 	}
	$prob=$prob/&factorial($k);
	return($prob);	
}

sub factorial {
	my ($n)=@_;
	if($n<0){ print "Negative factorial \n"; exit}
	my $i;
	my $prod;
	if($n==0 || $n==1){
		$prod=1;
	}
	else{
		$prod=1;
		for($i=2;$i<=$n;$i++){ 	$prod=$prod*$i;}
	}
	return($prod);
}

sub min {
	my @n=@_;
	my $min;
	if($n[0]>$n[1]){ $min=$n[1];}
	else{ $min=$n[0];}
	return($min);
}