def pdf_x(L, l, delta, dv):
    if -L / 2 + l / 2 + dv < delta < L / 2 - l / 2 - dv:
        return 2 * dv / l * ln(L / (L - l))
    if  L / 2 - l / 2 - dv < delta < L / 2 - l / 2 + dv:
        return 1 / 2. + 1 / l * ((delta - dv) * ln(L - l) + (delta + dv) * (1 - ln(2 * (delta + dv))) + 2 * dv * ln(L) - L / 2.)
    if (delta > L / 2. - l / 2. + dv) and (delta > L / 2. - l / 2. - dv):
        return 1 / l * (2 * dv * (1 + ln(L / 2.)) - (delta + dv) * ln(delta + dv) + (delta - dv) * ln(delta - dv))

def _get_I_total(R0, d, t):
    A = np.ln(R0)
    B = np.ln(1+ d) 
    mu = 0.5 * A / B
    I = 0.5 * np.exp( 0.25 * A**2 / B * np.sqrt( np.pi / B))
    I = I * np.sqrt(B) * (np.erf(t - mu)- np.erf(-mu))
    return I

def GKB_1(u_zref, zref, h, LAI, Wfol, Ta, pa):
   """Same as FKB_1, but then for spatial in- and output"""
   
   # Constants
   C_d = 0.2   # foliage drag coefficient
   C_t = 0.05  # heat transfer coefficient
   k = 0.41     # Von Karman constant
   Pr = 0.7    # Prandtl number
   hs = 0.009  # height of soil roughness obstacles (0.009-0.024)
   
   # Calculations
   Wsoil = 1.0 - Wfol
   h = ifthenelse(Wfol == 0.0, hs, h)
   
   z0 = 0.136 * h   # Brutsaert (1982)
   u_h0 = u_zref * ln(2.446) / ln ((zref - 0.667 * h)/z0) # wind speed at canopy height
   ust2u_h = 0.32 - 0.264 / exp(15.1 * C_d * LAI)
   ustarh = ust2u_h * u_h0
   nu0 = 1.327E-5 * (101325.0/pa) * (Ta / 273.15 + 1.0) ** 1.81 # kinematic viscosity
   n_h = C_d * LAI / (2.0 * ust2u_h ** 2.0)
   # First term
   F1st = ifthenelse(pcrne(n_h, 0.0), k * C_d / (4.0 * C_t * ust2u_h * (1.0 - exp(pcrumin(n_h)/2.0))) * Wfol ** 2.0, 0.0)
   # Second term
   S2nd = k * ust2u_h * 0.136 * Pr ** (2.0/3.0) * sqrt(ustarh * h / nu0) * Wfol ** 2.0 * Wsoil ** 2.0
   # Third term
   T3rd = (2.46 * (u_zref * k / ln(zref/hs) * hs / nu0) ** 0.25 - ln(7.4)) * Wsoil ** 2.0
   
   return F1st + S2nd + T3rd

def pdf_x_0(L, l, delta, dv):
    if -L / 2 + l / 2 + dv < delta < L / 2 - l / 2 - dv:
        return 1 / l * ln(L / (L - l))
    if  L / 2 - l / 2 - dv < delta < L / 2 - l / 2 + dv:
        return 0
    if (delta > L / 2. - l / 2. + dv) and (delta > L / 2. - l / 2. - dv):
        return 1 / l * ln(L / (2 * delta))

def probability_cut_nooverlaps(lc, lf, delta):
    """
        Probability that we cut fiber of nooverlapped fibers
    """
    # type I and III
    if (lf <= lc / 2.0 - delta) and (lf <= lc / 2.0 + delta):
        # print "Varianta I or III ",
        return lc / lf * (ln(lc / (lc - lf))) - 1
    # type II
    if (lf < lc / 2.0 + delta) and (lf > lc / 2.0 - delta):
        # print "Varianta II ",
        return (
            1
            / 2.0
            / lf
            * (
                (lc + 2.0 * delta) * ln(2.0 / (lc + 2.0 * delta))
                - (ln(lc - lf) + 1) * (lc - 2.0 * delta)
                + 2.0 * lc * ln(lc)
            )
        )
    # type IV
    if (lf >= lc / 2.0 + delta) and (lf >= lc / 2.0 - delta):
        # print "Varianta IV ",
        return 1 + 1 / lf * (
            -lc + (lc / 2.0 - delta) * ln((lc + 2 * delta) / (lc - 2 * delta)) + lc * ln(2 * lc / (lc + 2 * delta))
        )

def loadData(path="../data/",k=5,log='add',pca_n=0,SEED=34):
	from pandas import DataFrame, read_csv
	from numpy import log as ln
	from sklearn.cross_validation import KFold
	from sklearn.preprocessing import LabelEncoder
	from sklearn.preprocessing import StandardScaler
	train = read_csv(path+"train.csv")
	test = read_csv(path+"test.csv")
	id = test.id
	target = train.target
	encoder = LabelEncoder()
	target_nnet = encoder.fit_transform(target).astype('int32')
	feat_names = [x for x in train.columns if x.startswith('feat')]
	train = train[feat_names].astype(float)
	test = test[feat_names]
	if log == 'add':
		for v in train.columns:
			train[v+'_log'] = ln(train[v]+1)
			test[v+'_log'] = ln(test[v]+1)
	elif log == 'replace':
		for v in train.columns:
			train[v] = ln(train[v]+1)
			test[v] = ln(test[v]+1)      
	if pca_n > 0:
		from sklearn.decomposition import PCA
		pca = PCA(pca_n)
		train = pca.fit_transform(train)
		test = pca.transform(test)
	scaler = StandardScaler()
	scaler.fit(train)
	train = DataFrame(scaler.transform(train),columns=['feat_'+str(x) for x in range(train.shape[1])])
	test = DataFrame(scaler.transform(test),columns=['feat_'+str(x) for x in range(train.shape[1])])
	cv = KFold(len(train), n_folds=k, shuffle=True, random_state=SEED)
	return train, test, target, target_nnet, id, cv, encoder

def Ar(z, lf):
    #z is the dist between the 
    lf = np.float(lf)
    #######################
            #int#
    #######################
    
    if type(z) == int:
        if z == 0:
            return 1
        if z == lf / 2:
            return 0
        return lf * (lf - 2. * z + 2. * z * (.6931471806 + ln((1. / lf) * z))) / (lf - 2. * z) ** 2
    
    #######################
            #ARRAY#
    #######################
    
    res = 2 * z / lf * (ln(2 * z / lf) - 1) + 1
    if any(z == lf / 2):
        z = list(z)
        ind_h = z.index(lf / 2)
        res[ind_h] = 0
        z = np.array(z)
    if any(z == 0):
        z = list(z)
        ind_z = z.index(0) 
        res[ind_z] = 1 
        z = np.array(z)
    return res 

def PSIma(f, g):
   a = 0.33
   b = 0.41
   pi = 3.141592654
   tangens = scalar(atan((2.0 * g - 1.0) / sqrt(3.0))) * pi /180
   tangens = ifthenelse(tangens > pi/2.0, tangens - 2.0 * pi, tangens) 
   PSIma = ln(a + f) - 3.0 * b * f ** (1.0 / 3.0) + b * a ** (1.0 / 3.0) / 2.0 * ln((1 + g) ** 2.0 / (1.0 - g + sqrt(g))) + sqrt(3.0) * b * a ** (1.0 / 3.0) * tangens
   return PSIma

def Cw(hi, L, z0, z0h):
   alfa = 0.12
   beta = 125.0
   C0 = (alfa / beta) * hi
   C1 = pcrumin(z0h) / L
   C11 = -alfa * hi / L
   C21 = hi / (beta * z0)
   C22 = -beta * z0 / L
   C = ifthenelse(z0 < C0, pcrumin(ln(alfa)) + PSIh_y(C11) - PSIh_y(C1), ln(C21) + PSIh_y(C22) - PSIh_y(C1))
   Cw = ifthenelse(C < 0.0, 0.0, C) # This results from unfortunate parameter combination!
   return Cw

def le(L, l):
    '''
        Solve embedded length l_e of fibes (including null values) (integral)
    '''
    if L < l:
        #print 'very short specimen',
        return L / 4.
    if L < 2. * l:
        #print 'short specimen',
        return -L / 4. - L / 4. * ln(L / 2.) + L / 4. * ln(L) + 1 / 2. * l * (1 - L / (2. * l)) + 1 / 4. * L * ln(L / (2. * l)) + l / 4.
    if L >= 2. * l:
        #print 'long specimen',
        return -1 / 4. * l - 1 / 4. * L * ln(L - l) + 1 / 4. * L * ln(L)

def prob( lc, lf, delta ):
    #type I and III
    if ( lf <= lc / 2. - delta ) and ( lf <= lc / 2. + delta ):
        print "Varianta I or III ",
        return lc / lf * ( ln( lc / ( lc - lf ) ) ) - 1
    #type II
    if ( lf < lc / 2. + delta ) and ( lf > lc / 2. - delta ):
        print "Varianta II ",
        return 1 / 2. / lf * ( ( lc + 2. * delta ) * ln( 2. / ( lc + 2. * delta ) ) - ( ln( lc - lf ) + 1 ) * ( lc - 2. * delta ) + 2. * lc * ln( lc ) )
    #type IV
    if ( lf >= lc / 2. + delta ) and ( lf >= lc / 2. - delta ):
        print "Varianta IV ",
        return  1 + 1 / lf * ( -lc + ( lc / 2. - delta ) * ln( ( lc + 2 * delta ) / ( lc - 2 * delta ) ) + lc * ln( 2 * lc / ( lc + 2 * delta ) ) )

def u_pbl(NDVI):
   """Calculates Planetary Boundary Layer wind speed [m s-1] from NDVI
   
   NDVI Input PCRaster NDVI map (scalar, ratio between 0 and 1)"""
   
   z0m = 0.005 + 0.5 * (nd_mid/nd_max) ** 2.5
   assert z0m >= 0.0
   fc = ((nd_mid - nd_min) / nd_df) ** 2.0    # fractional vegetation cover == Wfol (-)
   assert fc >= 0.0
   h = z0m / 0.136                            # total height of vegetation (m)
   d = 2.0/3.0 * h			      # zero plane displacement (m)
   u_c = ln((z_pbl - d) / z0m) / ln((z_ms - d) / z0m)
   u_pbl = u_s * u_c
   return u_pbl, z0m, d, fc, h

 def _get_results(self):
     """dividing the interval <a,b>,
     returns aprox x, error estimation, No. of steps ..."""
     int = [self.a, self.b]
     if self.f(self.a) * self.f(self.b) > 0:
         print "None or more than 1 roots in selected interval"
     else:
         while abs(int[0] - int[1])/2 > self.error:
             if self.f(int[0]) * self.f((int[0] + int[1])*0.5) < 0:
                 int.insert(1,(int[0]+int[1])*0.5), int.pop()
             else:
                 int.insert(1,(int[0]+int[1])*0.5), int.pop(0) 
         return [(int[0] + int[1])/2,
                 abs(int[0] - int[1])/2,
                 (ln(2)-ln( (abs( int[0] - int[1] ) /2)/(self.b-self.a)))/ln(2)]

def calc_bethe_entropy(B):
    """Calculate the Bethe entropy given beliefs B"""
    Sbethe = 0
    for roti,rotj in rotedges:
        try:
            if B[(roti,rotj)]>0:
                Sbethe -= B[(roti,rotj)]* ln(B[(roti,rotj)])
        except RuntimeError:
            pass  # can't find edge. . .
    sumqi = sum([len(res2partners[resid])-1 for resid in resids])
    blogb = 0
    for roti in eg.GetVertexIDs():
        if B[roti]>0:
            blogb += B[roti]* ln(B[roti])
    Sbethe += sumqi*blogb
    return Sbethe

def calc_meanfield_entropy(B):
    """Calculate the mean-field entropy given beliefs B"""
    Smeanfield = 0
    for roti in eg.GetVertexIDs():
        if B[roti]>0:
            Smeanfield -= B[roti] * ln(B[roti])
    return Smeanfield

def significance(signal, background, min_bkg=0, highstat=True):

    if isinstance(signal, (list, tuple)):
        signal = sum(signal)
    if isinstance(background, (list, tuple)):
        background = sum(background)
    sig_counts = np.array(list(signal.y()))
    bkg_counts = np.array(list(background.y()))
    # reverse cumsum
    S = sig_counts[::-1].cumsum()[::-1]
    B = bkg_counts[::-1].cumsum()[::-1]
    exclude = B < min_bkg
    with np.errstate(divide='ignore', invalid='ignore'):
        if highstat:
            # S / sqrt(S + B)
            sig = np.ma.fix_invalid(np.divide(S, np.sqrt(S + B)),
                fill_value=0.)
        else:
            # sqrt(2 * (S + B) * ln(1 + S / B) - S)
            sig = np.sqrt(2 * (S + B) * np.ln(1 + S / B) - S)
    bins = list(background.xedges())[:-1]
    max_bin = np.argmax(np.ma.masked_array(sig, mask=exclude))
    max_sig = sig[max_bin]
    max_cut = bins[max_bin]
    return sig, max_sig, max_cut

 def get_eps_x_reinf (self, crack_x):
     self.cbs_allocation(abs(crack_x[0]))
     self.weakest_cb_check()
     if  crack_x[0]  not in self.unsorted_crack_list:
         self.unsorted_crack_list.append(crack_x[0])
         le_array, phi_array, f_array , tau_array = self.fiber_data_tuple[2][self.crack_list.index(abs(crack_x[0]))] , self.fiber_data_tuple[1][self.crack_list.index(abs(crack_x[0]))], self.fiber_data_tuple[3][self.crack_list.index(abs(crack_x[0]))], self.fiber_data_tuple[4][self.crack_list.index(abs(crack_x[0]))]
         diff = self.sum_of_fiberforces(crack_x, le_array, phi_array , f_array, tau_array)
         self.diff_list.append(diff)
     index_lists = self.unsorted_crack_list.index(crack_x[0])
     F = self.P - self.diff_list[index_lists]
     #Ar_fibers = self.active_fibers_list[index_lists] * Pi * self.r **2
     #########TESTPLOT#######################
     Af = len(le_array) * Pi * self.r ** 2 * 2
     plt.plot(crack_x, self.ar_list[index_lists] / Af)
     lf = self.fiber_length
     z = crack_x
     r = self.r
     testAr = 2 * Pi * r ** 2 * (-2 * z + z * ln(2 * z / lf) + lf) / lf
     #testAcki = (10 * [10 - 2. * z + 2. * z * (.6931471806 + ln((1 / 10) * z))]) / (10 - 2. * z) ** 2
     #plt.plot(crack_x, testAcki / (2 * Pi * r ** 2))
     plt.plot(crack_x, testAr / (2 * Pi * r ** 2))
     plt.plot(crack_x, Ar_alc(z, self.fiber_length))
     #maxAr = np.max(self.ar_list[index_lists])
     #minAr = np.min(self.ar_list[index_lists])
     #x_gerade = [0, 4.5]
     #y_gerade = [maxAr, minAr]
     #plt.plot(x_gerade, y_gerade)
     #plt.plot()
     plt.show()
     ################################
     Ar_fibers = self.ar_list[index_lists]
     eps = F / self.Er / Ar_fibers
     return eps * H(eps)

def Bw(hi, L, z0):
   # constants (Brutsaert, 1999)
   alfa = 0.12
   beta = 125.0
   
   # calculations
   B0 = (alfa / beta) * hi
   B1 = -1.0 *z0 / L
   B11 = -alfa * hi / L
   B21 = hi / (beta * z0)
   B22 = -beta * z0 / L
   tempB11 = PSIm_y(B11)
   tempB1 = PSIm_y(B1)
   B = ifthenelse(z0 < B0, -1.0 * ln(alfa) + PSIm_y(B11) - PSIm_y(B1), ln(B21) + PSIm_y(B22) - PSIm_y(B1))
   Bw = ifthenelse(B < 0.0, 0.0, B) # This results from unfortunate parameter combination!
   return Bw

def FKB_1(u_zref, zref, h, LAI, Wfol, Ta, pa):
   """Initial determination of roughness length for heat transfer (non-spatial)
   KB-1 function according to Massman, 1999
   Convention of variable names:
   f_z = f(z)
   d2h = d/h
   
   u_zref Input wind speed at reference height [m s-1]
   zref Input reference height [m]
   h Input canopy height [m]
   LAI Input canopy total Leaf Area Index [-]
   Wfol Input Fractional canopy cover [-]
   Ta Input ambient temperature [degrees Celsius]
   pa Input ambient air pressure [Pa]"""

   # Constants
   C_d = 0.2   # foliage drag coefficient
   C_t = 0.01  # heat transfer coefficient
   k = 0.41     # Von Karman constant
   Pr = 0.7    # Prandtl number
   hs = 0.009  # height of soil roughness obstacles (0.009-0.024)
   
   # Calculations
   Wsoil = 1.0 - Wfol
   if Wfol == 0.0: # for bare soil take soil roughness
      h = hs
   assert Wfol >= 0.0 and Wfol <= 1.0 and Wsoil >= 0.0 and Wsoil <= 1.0
   z0 = 0.136 * h   # Brutsaert (1982)
   u_h0 = u_zref * ln(2.446) / ln ((zref - 0.667 * h) / z0) # wind speed at canopy height
   u_h0 = cellvalue(u_h0, 0, 0)
   u_h0 = u_h0[0]
   assert u_h0 >= 0.0
   ust2u_h = 0.32 - 0.264/exp(15.1 * C_d * LAI)
   ustarh = ust2u_h * u_h0
   nu0 = 1.327E-5 * (101325.0 / pa) * (Ta / 273.15 + 1.0) ** 1.81 # kinematic viscosity
   n_h = C_d * LAI / (2.0 * ust2u_h ** 2.0)
   # First term
   if n_h <> 0.0:
      F1st = k * C_d / (4.0 * C_t * ust2u_h * (1.0 - exp(pcrumin(n_h)/2.0))) * Wfol ** 2.0
   else:
      F1st = 0.0
   # Second term
   S2nd = k * ust2u_h * 0.136 * Pr ** (2.0/3.0) * sqrt(ustarh * h / nu0) * Wfol ** 2.0 * Wsoil ** 2.0
   # Third term
   T3rd = (2.46 * (u_zref * k / ln(zref/hs) * hs / nu0) ** 0.25 - ln(7.4)) * Wsoil ** 2.0
   
   return F1st + S2nd + T3rd

 def mean_approx(self,l):
     yarn = self._get_gbundle_props() 
     # bundle length
     l_b = yarn[0]
     # mean and stdev of Gaussian bundle
     mu_b = yarn[1]
     gamma_b = yarn[2]
     # No. of bundles in series 
     nb = l/l_b
     if nb == 1:
         return mu_b, mu_b, mu_b
     w = ln(nb)
     # approximation of the mean for k = (1;300) (Mirek)
     mu = mu_b + gamma_b * (-0.007 * w**3 + 0.1025 * w**2 - 0.8684 * w)
     ### approximation for the mean from extremes of Gaussian (tends to Gumbel as mb grows large)
     a = gamma_b / sqrt(2*w)
     b = mu_b + gamma_b * ( (ln(w) + ln(4*pi))/ sqrt(8 * w) - sqrt(2 * w))
     med = b + a * ln(ln(2))
     mean = b - 0.5772156649015328606 * a
     return mu, mean, med