import sys
import math
import random
import os
import shutil
import scipy
import scipy.misc

## set of size 2^i, adding M parity constraints
## returns average number of surviving solutions
def mean(i,M):
	return pow(2,(i - M))
	
def w_star(L,i):
	j = 1
	tmp =  scipy.misc.comb(L,j)
	while (tmp<pow(2.0,i)-1):
		j = j +1
		tmp =tmp + scipy.misc.comb(L,j)
	#print  "wstar=",j, " tot=",tmp, " S=",pow(2.0,i)
	return j
		

## upper bound on the power, from ICML paper and markov chain analysis
def po_rw(L,i,M,t, verbose=True):
	tmp = 0.0
	tot_conf = 0.0
	wstar=w_star(L,i)
	if verbose:
		print "wstar=",wstar
	## fraction of 1 in the matrix
	pos_fract = t
	
	for w in range(1,wstar):
		base = 0.5 + 0.5*((1.0-2*pos_fract)**w)
		tmp = tmp + scipy.misc.comb(L,w)* (base** M)
		tot_conf += scipy.misc.comb(L,w)
		if verbose:
			print (w,(base** M))
			print scipy.misc.comb(L,w)
		#print  "power=",tmp
	
	## add rest at maximum distance
	w=wstar
	base = 0.5 + 0.5*((1.0-2*pos_fract)**w)
	if verbose:
		print (w,(base** M))
		print (pow(2.0,i)-tot_conf)
	tmp = tmp + (pow(2.0,i)-1-tot_conf)* (base** M)
	
	tmp = tmp *pow(2.0,i-M)
	tmp = tmp +pow(2.0,i-M)
	if verbose:
		print  "power=",tmp,"(",tot_conf,")"
	return tmp
		

## variance	= power - mean^2
# L binary variables
# set of size 2^i solutions
# adding M xors
# each variable added with probability t
# returns upper bound on the variance on the number of surviving solutions
def v(L,i,M,t):
	return po_rw(L,i,M,t) - mean(i,M)*mean(i,M)

#####################################
## computes minimum value of f that
##################################### 
def f_star(L,M, c=2,step=0.01, max_failure_prob=0.25):
	i = M + c
	#if M<c:
	#	return 0.5
	for f in range(1,int(math.ceil(0.5/step))+1):
		sigmas2 = v(L,i,M,step*f)
		print sigmas2
		#cantelli
		failure_prob = sigmas2 / (sigmas2 + mean(i,M)*mean(i,M))
		if failure_prob<=max_failure_prob:
			##print "variance ", sigmas2
			return step*f
	return 0.5


print f_star(100,50,max_failure_prob=0.45,step=0.01)