def theta(self, S, t, vol, r):
     if t > self.T:
         return np.zeros(len(S))
     if self.op_type == 'c':
         return np.subtract(
                 np.true_divide(
                     np.multiply(-vol, 
                         np.multiply(S, norm.pdf(self.d1(S, t, vol, r)))), 
                     np.multiply(2,
                         np.power(
                             np.subtract(self.T, t), .5))),
                 np.multiply(r,
                     np.multiply(self.K,
                         np.multiply(
                             np.exp(
                                 np.multiply(-r,
                                     np.subtract(self.T, t))),
                             norm.cdf(self.d2(S, t, vol, r))))))
     else:
         return np.add(
                 np.true_divide(
                     np.multiply(-vol, 
                         np.multiply(S, norm.pdf(self.d1(S, t, vol, r)))), 
                     np.multiply(2,
                         np.power(
                             np.subtract(self.T, t), .5))),
                 np.multiply(r,
                     np.multiply(self.K,
                         np.multiply(
                             np.exp(
                                 np.multiply(-r,
                                     np.subtract(self.T, t))),
                             norm.cdf(
                                 np.multiply(-1, self.d2(S, t, vol, r)))))))

def _hist_bin_doane(x):
    """
    Doane's histogram bin estimator.

    Improved version of Sturges' formula which works better for
    non-normal data. See
    stats.stackexchange.com/questions/55134/doanes-formula-for-histogram-binning

    Parameters
    ----------
    x : array_like
        Input data that is to be histogrammed, trimmed to range. May not
        be empty.

    Returns
    -------
    h : An estimate of the optimal bin width for the given data.
    """
    if x.size > 2:
        sg1 = np.sqrt(6.0 * (x.size - 2) / ((x.size + 1.0) * (x.size + 3)))
        sigma = np.std(x)
        if sigma > 0.0:
            # These three operations add up to
            # g1 = np.mean(((x - np.mean(x)) / sigma)**3)
            # but use only one temp array instead of three
            temp = x - np.mean(x)
            np.true_divide(temp, sigma, temp)
            np.power(temp, 3, temp)
            g1 = np.mean(temp)
            return x.ptp() / (1.0 + np.log2(x.size) +
                                    np.log2(1.0 + np.absolute(g1) / sg1))
    return 0.0

def compute_IICR_n_islands(n, M, t, s=True):
    # This method evaluates the lambda function in a vector
    # of time values t.
    # If 's' is True we are in the case when two individuals where
    # sampled from the same island. If 's' is false, then the two
    # individuals where sampled from different islands.

    # Computing constants
    gamma = np.true_divide(M, n - 1)
    delta = (1 + n * gamma) ** 2 - 4 * gamma
    alpha = 0.5 * (1 + n * gamma + np.sqrt(delta))
    beta = 0.5 * (1 + n * gamma - np.sqrt(delta))

    # Now we evaluate
    x_vector = t
    if s:
        numerator = (1 - beta) * np.exp(-alpha * x_vector) + (alpha - 1) * np.exp(-beta * x_vector)
        denominator = (alpha - gamma) * np.exp(-alpha * x_vector) + (gamma - beta) * np.exp(-beta * x_vector)
    else:
        numerator = beta * np.exp(-alpha * (x_vector)) - alpha * np.exp(-beta * (x_vector))
        denominator = gamma * (np.exp(-alpha * (x_vector)) - np.exp(-beta * (x_vector)))

    lambda_t = np.true_divide(numerator, denominator)

    return lambda_t

def this_create_total_dict():
    input_dir='/tmp/fvs_fun/fv_deep_funneled/'
    fv_dict = dict()
    count=0;
    for name in os.listdir(input_dir):
        count=count+1
        input_file = os.path.join(input_dir,name)
        with open(input_file) as f:
            #w, h = [float(x) for x in f.readline().split()]
            #array = [[float(x) for x in line.split()] for line in f]
            s= numpy.genfromtxt(input_file, delimiter=',')
            

            word1 = "".join(re.findall("[a-zA-Z]+", name))
            word = word1[:-3]
            if word in fv_dict:
                fv_dict[word].append([numpy.true_divide(s,numpy.linalg.norm(s,ord=2))])
                #fv_dict[word]= numpy.concatenate(  (fv_dict[word],  numpy.true_divide(s,numpy.linalg.norm(s,ord=2))) ,axis=1)
            else:
                fv_dict[word]=[numpy.true_divide(s,numpy.linalg.norm(s,ord=2))]

            print count
            if count>10:
                break;
    return fv_dict

 def __call__(self, values, clip=True, out=None):
     values = _prepare(values, clip=clip, out=out)
     np.multiply(values, np.log(self.exp + 1.), out=values)
     np.exp(values, out=values)
     np.subtract(values, 1., out=values)
     np.true_divide(values, self.exp, out=values)
     return values

def thresholdColor(img):
    NTIME = time.time()

    channels = cv2.split(img)
    red = channels[2].astype(np.uint8)
    green = channels[1].astype(np.uint8)
    blue = channels[0].astype(np.uint8)

    r_g = np.true_divide(red,green)
    r_b = np.true_divide(red,blue)
    g_r = np.true_divide(green,red)
    g_b = np.true_divide(green,blue)

    
    combs = []
    for col in colors:
        (dom,minv,first,second)=color_dims[col]
        if dom == "r":
            comb = ((red > minv) & (r_g > first) & (r_b > second)).astype(np.uint8)
        elif dom == "g":
            comb = ((green > minv) & (g_r > first) & (g_b > second)).astype(np.uint8) 
        combs.append(comb)
        
    BTIME = time.time()
    rospy.loginfo("IN")
    rospy.loginfo(BTIME-NTIME)
    
    return combs

    def __call__(self, values, clip=None):

        if clip is None:
            clip = self.clip

        if isinstance(values, ma.MaskedArray):
            if clip:
                mask = False
            else:
                mask = values.mask
            values = values.filled(self.vmax)
        else:
            mask = False

        # Make sure scalars get broadcast to 1-d
        if np.isscalar(values):
            values = np.array([values], dtype=float)
        else:
            # copy because of in-place operations after
            values = np.array(values, copy=True, dtype=float)

        # Normalize based on vmin and vmax
        np.subtract(values, self.vmin, out=values)

        np.true_divide(values, self.vmax - self.vmin, out=values)

        # Clip to the 0 to 1 range
        if self.clip:
            values = np.clip(values, 0., 1., out=values)

        # Stretch values
        values = self.stretch(values, out=values, clip=False)

        # Convert to masked array for matplotlib
        return ma.array(values, mask=mask)

def this_create_total_dict(fromPercent, toPercent):
    input_dir="/tmp/fvs_fun/fv_deep_funneled/"
    total = len(os.listdir(input_dir))
    start_index = int(fromPercent*total)
    end_index =  int(toPercent*total)
    fv_dict = dict()
    count=0;
    
    for name in os.listdir(input_dir)[start_index:end_index]:
        count=count+1
        input_file = os.path.join(input_dir,name)
        
        with open(input_file) as f:
            #w, h = [float(x) for x in f.readline().split()]
            #array = [[float(x) for x in line.split()] for line in f]
            s= numpy.genfromtxt(input_file, delimiter=',')
            word1 = "".join(re.findall("[a-zA-Z]+", name))
            word = word1[:-3]
            if word in fv_dict:
                fv_dict[word].append([numpy.true_divide(s,numpy.linalg.norm(s,ord=2))])
                #fv_dict[word]= numpy.concatenate(  (fv_dict[word],  numpy.true_divide(s,numpy.linalg.norm(s,ord=2))) ,axis=1)
            else:
                fv_dict[word]=[numpy.true_divide(s,numpy.linalg.norm(s,ord=2))]
            print count
        #if len(fv_dict.keys())>50:
        #    break;
    return fv_dict

    def general_work(self, input_items, output_items):
        in0 = input_items[0]
	flags = input_items[1]
        out = output_items[0]
	if len(in0) !=len(flags):
	  print "TX:input length buffer DOESN'T EQUAL"
	len_in = min(len(in0), len(flags))
	if self.state == 0:
	  for i in range(len_in):
	    if flags[i] == 1:
	      print "TX: FLAG FOUND!!!", i
              self.state =1
	      self.consume(0,i+1)
	      self.consume(1,i+1)
              return 0
        if self.state == 1:
	  Y = in0[-self.pilot_length:0]
	  corr = np.true_divide(np.dot(Y.transpose(),self.pilot_seq.conj()),self.pilot_length)
	  A = np.dot(corr,corr.transpose().conj())
	  B = np.true_divide(np.dot(Y.transpose(),Y.conj()),self.pilot_length)
	  Omega_hat = B-A
	  Omega = Omega_hat - self.NHAT + np.identity(self.nt)
	  Sigma = np.dot(np.linalg.pinv(Omega)-np.linalg.pinv(B),self.weight) #nt x nt
	  #=====================debugging msg========================
	  print "Pilot detected! TX sync Sigma ="
	  print Sigma
          out[:]=Sigma.reshape(self.nt*self.nt)
	  self.consume(0,len_in)
	  self.consume(1,len_in)
          return 1

def compare_cdf_g(model, n_obs=10000, case='T2s',
                  t_vector=np.arange(0, 100, 0.1), path2ms='./'):
    # Do a graphical comparison by plotting the theoretical cdf and the 
    # empirical pdf
    T_list_ms = np.true_divide(model.T_list,2)
    M_list = model.M_list
    if case == 'T2s':
        cmd = create_ms_command_T2(n_obs, model.n, T_list_ms, M_list, 'same')
        theor_cdf = model.cdf_T2s
    elif case == 'T2d':
        cmd = create_ms_command_T2(n_obs, model.n, T_list_ms, M_list, 'disctint')
        theor_cdf = model.cdf_T2d
    else:
        return 1
    
    obs = simulate_T2_ms(cmd, path2ms)
    obs = np.array(obs) * 2
    
    delta = t_vector[-1] - t_vector[-2]
    bins = t_vector
    
    f_obs = np.histogram(obs, bins=bins)[0]
    cum_f_obs = [0] + list(f_obs.cumsum())
    F_obs = np.true_divide(np.array(cum_f_obs), n_obs)
    F_theory = [theor_cdf(t) for t in bins]    
    
    # Now we plot
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.plot(bins, F_obs, '-r', label = 'Empirical cdf')
    ax.plot(bins, F_theory, '-b', label = 'Theoretical cdf')
    
    plt.legend()
    plt.show()

 def __init__(self, n=10, M=0.1):
     T2_StSI.__init__(self, n, M)
     [A, B, a, alpha, beta, E, AplusB, AminusB] = \
     T2_StSI.compute_constants(self, n, M)
     [self.a, self.alpha, self.beta] = [a, alpha, beta]
     self.gamma = np.true_divide(M, n-1)
     self.c = np.true_divide(self.gamma, beta-alpha)

def normalize(ndarray_list):
	""" Normalize each ndarray and scale all
	> All ndarray must have same dimensions!

	Parameters :
	- ndarray_list : ndarray list of images


	Returns :
	- ndarray : ndarray in wich each pixel match to pixels normalize

	"""
	liste = []
	lenght = len(ndarray_list)
	print "Taille de la liste d'image à traiter : "+ str(lenght)
	# 1) Normalize each flat field, i.e. divide it by its mean entry value
	for i in range(lenght):
		print "Normalize 1/2 : " + str(i+1) + "/" + str(lenght)
		mean = np.mean(ndarray_list[i]) # Find the mean of the ndarray
		print "Moyenne de la " + str(i+1) + " image : " + str(mean)
		liste.append(mean)
		ndarray_list[i] = np.true_divide(ndarray_list[i], mean) # Divide ndarray by its mean, normalizing it
	# 2) Scale all of the fields’ means so that their individual averages are equal to one another
	meanofmean = sum(liste) / len(liste)  # Find the mean of the total set of ndarray_list
	print "Moyenne global de toutes les images : " + str(meanofmean)
	for i in range(lenght):
		print "Normalize 2/2 with : "+ str(i+1) + "/" + str(lenght)
		ndarray_list[i] = np.multiply(ndarray_list[i], np.true_divide(meanofmean, np.mean(ndarray_list[i]))) # Divide
	return ndarray_list

def sigmoid_inverse(z):
    # z = 0.998*z + 0.001
    assert(np.max(z) <= 1.0 and np.min(z) >= 0.0)
    z = 0.998*z + 0.001
    return np.log(np.true_divide(
        1.0, (np.true_divide(1.0, z) - 1)
    ))

def ArtistsRate(JCR, period, istest=False):    
    global constants
    cols = ["id", "reviews_avgnote", "reviews_all", "playlisters_logged", "starers_logged", "listeners_logged", "downloaders_logged"]
    if period == 'total': cols += ['days_since_publication']
    COLUMNS = JCR.getColumns(cols)

    
    #REVIEWS RATE
    reviews_reduceunder = constants[period]['reviews_reduceunder']    
    #for album in JCR.iterAlbumJoinArtist() 
    #bestreviewed
    COLUMNS['reviews_rate'] = computeBayesAvg(COLUMNS["reviews_avgnote"], COLUMNS['reviews_all'], constants[period]['reviews_bayes_const'], reviews_reduceunder) 

    #LIKESRATE
    likes = COLUMNS["playlisters_logged"] + 1.5*COLUMNS["starers_logged"] #*1.5 to bring the 2 values-set to about the same level avg level 
                                            
    ratios =  np.true_divide(likes, COLUMNS["listeners_logged"]+1)
    likes_reduceunder = constants[period]['likes_reduceunder']
    COLUMNS['likes_rate'] = computeBayesAvg(ratios, COLUMNS["listeners_logged"]+1, constants[period]['likes_bayes_const'], likes_reduceunder)

    #FINAL RATE COMBINING BY RANK RATE
    if period=='total': 
        days = 1 + nomramlizeTo0_1(np.log2( COLUMNS['days_since_publication'] + 2 ))
        downloaders = np.true_divide(COLUMNS['downloaders_logged'], days) 
    else: downloaders = COLUMNS['downloaders_logged']
    
    COLUMNS['rate'] = getRanks(COLUMNS['likes_rate']) + 1.5*getRanks(COLUMNS['reviews_rate']) + getRanks(COLUMNS['downloaders_logged']) 
    
    
    if istest: 
        return COLUMNS
    else: 
        return {'id': COLUMNS['id'], 'rate': COLUMNS['rate']}

def persp(winsize, nf, dim=3):
    """
    Assumptions: 2D: Device coordinates (eye coordinates) have
    (0,0) at the top left corner. z is just normalized.
    This is because in 2D I like to pretend the screen is a big bitmap.

    3D: I am not going to bother explaining what is going on here.
        It is relatively simple but this was incredibly painful to implement
        on top of poorly written code.
        It is dull mathematics. One thing to note: If you are using your own
        geometry shader you MUST do some sort of clipping plane check
        on the output, because if you don't you will get weird artifacts around the near
        plane because of the asymptotic behavior.

    :param winsize: width, height
    :return: device coordinates normalized to the x,y axis
    """

    w = np.true_divide(2, winsize[0])
    h = np.true_divide(2, winsize[1])
    n, f = nf
    if dim == 3:
        return np.array([[w*n, 0, 0, 0],
                [0, h*n, 0, 0],
                [0, 0, -(n+f)/float(f-n), -2*n*f/float(f-n)],
                [0, 0, -1, 0]], dtype="float32")

    else:
        return [[2*w, 0, -1],
                [0, -2*h, 1],
                [0, 0, 1]]

def get_clipped_resized_masks(boxes, input_masks):
  target_size = cfg.MASK_SIZE
  num_boxes = boxes.shape[0]
  num_channels = input_masks.shape[1]
  #masks = cymask.clip_resize_mask(sp, reg2sp, boxes, target_size)
  masks = np.zeros((num_boxes, num_channels, target_size, target_size))
  totaltime=0.
  yv = np.arange(target_size).reshape(-1,1)
  xv = np.arange(target_size).reshape(1,-1)
  yv = np.true_divide(yv, float(target_size)-1.)
  xv = np.true_divide(xv, float(target_size)-1.)
  
  for k in range(num_boxes):
    box = boxes[k,:]
    xmin = box[0]
    ymin = box[1]
    w = box[2]-box[0]
    h = box[3]-box[1]
    xv1 = np.round(xmin + xv*w).astype(int)
    yv1 = np.round(ymin + yv*h).astype(int)
    for l in range(num_channels):
      M = input_masks[k,l,yv1, xv1]
      masks[k,l,:,:]=M
    #S = sp[box[1]:box[3]+1,box[0]:box[2]+1]
    #M = reg2sp[:,k]
    #M = M[S-1]
    #start = time.clock()
    #masks[k,:,:] = scipymisc.imresize(M, (target_size, target_size), 'nearest')
    #stop = time.clock()
    #totaltime+=stop-start 
  return masks

def hand_to_features(hand, convolution_function, play="river"):
    # hand: list representation of n hands e.g., ['Kc','Ac','6c','7h']
    # returns the output of the filters

    diction = {"river": 7, "turn": 6, "flop": 5}
    flush_filter = numpy.ones((1, 13))
    vertical_filter = numpy.ones((4, 1))
    horizontal_filter = numpy.ones((1, 5))

    if play not in diction.keys():
        print "play must be river, turn or flop"
        quit()
    hand_matrix = string_to_vector(hand)

    column_sums = numpy.sum(hand_matrix, axis=0)
    column_sums[column_sums > 1] = 1
    column_sums = column_sums.reshape((1, 14))
    stripped_hand_matrix = numpy.delete(hand_matrix, 0, 1)

    straight_features = convolution_function(column_sums, horizontal_filter)
    straight_features_normalized = numpy.true_divide(straight_features, 5)
    two_three_four_features = convolution_function(stripped_hand_matrix, vertical_filter)
    two_three_four_features_normalized = numpy.true_divide(two_three_four_features, 4)
    flush_features = convolution_function(stripped_hand_matrix, flush_filter).T
    flush_features_normalized = numpy.true_divide(flush_features, diction[play])

    return numpy.concatenate(
        (straight_features_normalized, flush_features_normalized, two_three_four_features_normalized), axis=1
    )

def make_bdt(data,label,weight):

    # create and save tree model

    # configure weights
    wp = len(np.where(label == 1)[0])
    wd = len(np.where(label == 0)[0])
    print 'Scale pos. weight: {}'.format(3.*np.true_divide(wd,wp))

    # setup parameters for xgboost
    param = {}
    # use logistic regression loss, use raw prediction before logistic transformation
    # since we only need the rank
    param['objective'] = 'binary:logistic'
    # scale weight of positive examples
    param['scale_pos_weight'] = 3.*np.true_divide(wd,wp)
    param['eta']               = 0.05
    param['eval_metric']       = 'error'
    param['silent']            = 1
    param['nthread']           = 6
    param['min_child_weight']  = 2
    param['max_depth']         = 10
    param['gamma']             = 0.0
    param['colsample_bytree']  = 0.8
    param['subsample']         = 0.9
    param['reg_alpha']         = 1e-5

    # boost 25 tres
    num_round = 300
    
    # make dmatrices from xgboost
    dtrain = xgb.DMatrix( data, label=label )
    bst    = xgb.train(plst, dtrain, num_round)
        
    return bst

def yes(features, knowledge, smooth):
    
    #calculate probability of yes
    p_y = math.log(np.true_divide(float(knowledge["nmb_y"]), (knowledge["nmb_y"] + knowledge["nmb_n"])))
    p_f_y = p_y
    #calculate probability of no
    p_n = math.log(np.true_divide(float(knowledge["nmb_n"]), (knowledge["nmb_y"] + knowledge["nmb_n"])))
    p_f_n = p_n
    for e in features:
        #calculate p(yes|features)
        if (not knowledge["frqy"].has_key(e)) and (not knowledge["frqn"].has_key(e)):
            continue
        tmp = (0 if (not knowledge["frqy"].has_key(e)) else knowledge["frqy"][e]) + smooth
        if tmp == 0:
            p_f_y =float('-inf')
        else:
            p_f_y = p_f_y + math.log(float(tmp)/(len(knowledge["frqy"]) * smooth + knowledge["total_f_y"]))
        
        #calculate p(no|features) 
        tmp = (0 if (not knowledge["frqn"].has_key(e)) else knowledge["frqn"][e]) + smooth
        if tmp == 0:
            p_f_n =float('-inf')
        else:
            p_f_n = p_f_n + math.log(float(tmp)/(len(knowledge["frqn"]) * smooth + knowledge["total_f_n"]))

    if (p_f_y == float('-inf') and p_f_n == float('-inf')):
        rlt = 'yes' if p_y > p_n else 'n0'     
    else:
        rlt = 'yes' if p_f_y > p_f_n else 'no' 
    return rlt

def run_cv(data,label,weight):

    # run cross validation on training set to get stats

    # configure weights
    wp = len(np.where(label == 1)[0])
    wd = len(np.where(label == 0)[0])
    print 'Scale pos. weight: {}'.format(3.*np.true_divide(wd,wp))

    # setup parameters for xgboost
    param = {}
    # use logistic regression loss, use raw prediction before logistic transformation
    # since we only need the rank
    param['objective'] = 'binary:logistic'
    # scale weight of positive examples
    param['scale_pos_weight'] = 3.*np.true_divide(wd,wp)
    param['eta']               = 0.05
    param['eval_metric']       = 'error'
    param['silent']            = 1
    param['nthread']           = 6
    param['min_child_weight']  = 2
    param['max_depth']         = 10
    param['gamma']             = 0.0
    param['colsample_bytree']  = 0.8
    param['subsample']         = 0.9
    param['reg_alpha']         = 1e-5

    # boost 25 tres
    num_round = 300

    test_error    = []
    test_falsepos = []
    test_falseneg = []
    scores        = np.zeros((2,len(label)))

    # get folds
    skf = StratifiedKFold(label, 10, shuffle=True)
    for i, (train, test) in enumerate(skf):
        #print train, test
        Xtrain = data[train]
        ytrain = label[train]
        wtrain = label[train]
        Xtest  = data[test]
        ytest  = label[test]
        wtest = label[test]
        # make dmatrices from xgboost
        dtrain = xgb.DMatrix( Xtrain, label=ytrain )
        dtest  = xgb.DMatrix( Xtest )
        #watchlist = [ (dtrain,'train') ]

        bst   = xgb.train(param, dtrain, num_round)
        ypred = bst.predict(dtest)
        fold_error,fold_falsepos,fold_falseneg = compute_stats(ytest,ypred)
        test_error.append(fold_error)
        test_falsepos.append(fold_falsepos)
        test_falseneg.append(fold_falseneg)
        scores[0,test] = ytest
        scores[1,test] = ypred

    return test_error,test_falsepos,test_falseneg,scores