def montecarlo(self,N):
     ppf = MFnLineArray(xdata=self.hui_cdf, ydata=self.x)
     fragment_lengths = []
     for i in range(N):
         fragment_length = ppf.get_value(np.random.rand(1)) + ppf.get_value(np.random.rand(1))
         fragment_lengths.append(fragment_length / 2.)
     return np.array(fragment_lengths)

 def ppf( self, x ):
     self.check()
     if len( self.cdf_values ) != 0:
         xx, cdf = self.x_values, self.cdf_values
     else:
         xx, cdf = self.x_values, self.cdf( self.x_values )
     ppf_mfn_line = MFnLineArray( xdata = cdf, ydata = xx )
     return ppf_mfn_line.get_values( x )

 def strains( self ):
     mfn = MFnLineArray()
     mfn.xdata, mfn.ydata = self.values
     strains_fn = frompyfunc( mfn.get_diff, 1, 1 )
     strains = strains_fn( mfn.xdata )
     strains[0] = strains[1]
     strains[-2] = strains[-1]
     return strains

 def _get_cached_pdf( self ):
     if len( self.pdf_values ) == 0:
         cdf_line = MFnLineArray( xdata = self.x_values, ydata = self.cdf_values )
         pdf = []
         for x in self.x_values:
             pdf.append( cdf_line.get_diff( x ) )
         return MFnLineArray( xdata = self.x_values, ydata = pdf )
     else:
         return MFnLineArray( xdata = self.x_values, ydata = self.pdf_values )

 def get_epsf_x(self, w):
     '''
     evaluates the mean fiber strain profile at given load
     '''
     self.CB_model.w = w
     if self.CB_model.Ll > self.CB_model.Lr:
         epsf_interp = MFnLineArray(xdata=-self.CB_model._x_arr[::-1], ydata=self.CB_model._epsf_arr[::-1])
     else:
         epsf_interp = MFnLineArray(xdata=self.CB_model._x_arr, ydata=self.CB_model._epsf_arr)
     return epsf_interp.get_values(self.x)

 def get_epsm_x(self, w):
     """
     evaluates the matrix strain profile at given load
     """
     self.CB_model.w = w
     if self.CB_model.Ll > self.CB_model.Lr:
         epsm_interp = MFnLineArray(xdata=-self.CB_model._x_arr[::-1], ydata=self.CB_model._epsm_arr[::-1])
     else:
         epsm_interp = MFnLineArray(xdata=self.CB_model._x_arr, ydata=self.CB_model._epsm_arr)
     return epsm_interp.get_values(self.x)

 def w_x_results( self, w_arr, x ):
     epsm = np.zeros( ( len( w_arr ), len( x ) ) )
     mu_epsf = np.zeros( ( len( w_arr ), len( x ) ) )
     sigma_c = []
     for i, w in enumerate( w_arr ):
         self.model.w = w
         epsm_line = MFnLineArray( xdata = self.x_arr, ydata = self.epsm_arr )
         mu_epsf_line = MFnLineArray( xdata = self.x_arr, ydata = self.mu_epsf_arr )
         epsm[i, :] = epsm_line.get_values( x )
         mu_epsf[i, :] = mu_epsf_line.get_values( x )
         sigma_c.append( self.sigma_c )
     return epsm, mu_epsf, np.array( sigma_c )

    def _get_pdf_array( self ):

        # get the normed histogram implemented in the base class YMBHist
        h, b = self.normed_hist
        b_cen = ( b[1:] + b[:-1] ) / 2.
        mask = ( h != 0. )
        h = h[mask]
        b_cen = b_cen[mask]
        b = hstack( [min( b ), b_cen, max( b )] )
        h = hstack( [0, h, 0] )
        h = h / trapz( h, b )
        p = MFnLineArray( xdata=b, ydata=h )
        return p.get_values( self.x_array )

 def _get__epsm_arr(self):
     if len(self.sorted_reinf_lst[0]) != 0 and len(self.sorted_reinf_lst[1]) != 0:
         epsm_cont_interp = MFnLineArray(xdata=self.cont_fibers.x_arr, ydata=self.cont_fibers.epsm_arr)
         epsm_short_interp = MFnLineArray(xdata=self.short_fibers.x_arr, ydata=self.short_fibers.epsm_arr)
         added_epsm_cont = self.cont_fibers.epsm_arr + epsm_short_interp.get_values(self.cont_fibers.x_arr) 
         added_epsm_short = self.short_fibers.epsm_arr + epsm_cont_interp.get_values(self.short_fibers.x_arr) 
         sorted_unique_idx = np.unique(np.hstack((self.cont_fibers.x_arr, self.short_fibers.x_arr)), return_index=True)[1]
         return np.hstack((added_epsm_cont, added_epsm_short))[sorted_unique_idx]
     elif len(self.sorted_reinf_lst[0]) != 0:
         return self.cont_fibers.epsm_arr
     elif len(self.sorted_reinf_lst[1]) != 0:
         self.short_fibers.w = self.w
         return self.short_fibers.epsm_arr

    def get_df_average( self, n_points ):
        '''derive the average phi-function based on all entries 
        in damage_function_list
        '''

        def get_y_average( self, x_average ): 
            '''get the y-values from the mfn-functions in df_list for 
            'x_average' and return the average.
            Note that the shape of 'mfn.xdata' does not necessarily needs to be equal in all
            'DamageFunctionEntries' as the number of steps used for calibration or the adaptive 
            refinement in 'tloop' might have been different for each case.
            '''  
            y_list = [ self.damage_function_list[i].damage_function.get_value( x_average ) \
                       for i in range(len( self.damage_function_list )) ]
            return sum(y_list) / len(y_list)
    
        get_y_average_vectorized = frompyfunc( get_y_average, 2, 1 )
        
        mfn = MFnLineArray()
        
        # take the smallest value of the strains for the average function. Beyond this value 
        # the average does not make sense anymore because it depends on the arbitrary number
        # of entries in the df_list  
        #
        xdata_min = min( self.damage_function_list[i].damage_function.xdata[-1] \
                         for i in range( len( self.damage_function_list ) ) )
        
        # number of sampling point used for the average phi function
        #
        mfn.xdata = linspace( 0., xdata_min, num = n_points )

        # get the corresponding average ydata values
        #
        mfn.ydata = self.get_y_average_vectorized( mfn.xdata )

        return mfn

 def get_sigma_m_x_input( self, sigma ):
     self.apply_load( sigma )
     line = MFnLineArray( xdata = self.x_arr,
                         ydata = self.epsm_arr )
     return line.get_values( self.x_input )

 def get_L( self, Ll, Lr ):
     self.result_values
     BC_line = MFnLineArray( xdata = self.BC_range, ydata = self.BC_range, extrapolate = 'constant' )
     return BC_line.get_values( [Ll, Lr] )

 def _mfn_line_array_target_default(self):
     mfn = MFnLineArray( xdata = self.xdata_target, ydata = self.ydata_target )
     mfn.data_changed = True
     return mfn

 def eta(self, psi):
     psi_line = MFnLineArray(xdata=self.psi_line.ydata,
                             ydata=self.psi_line.xdata)
     return psi_line.get_values(psi, k=1)

                       3.75724285,  3.84216067,  3.92147416,  3.997584  ,  4.07213255,  4.1462654,
                       4.22079085,  4.29628189,  4.37314373,  4.4516601 ,  4.53202557,  4.61436863,
                       4.69876847,  4.78526735,  4.87387975,  4.96459931,  5.0574041 ,  5.15226056,
                       5.24912666,  5.3479542 ,  5.44869074,  5.55128092,  5.65566766,  5.76179299,
                       5.86959871,  5.97902693,  6.09002044,  0.74043954,  0.74661105,  0.75304179,
                       0.75972494,  0.7666537 ,  0.77382139,  0.7812214 ,  0.78884722,  0.79669249,
                       0.80475096,  0.81301651,  0.82148317,  0.83014509,  0.83899658,  0.84803208,
                       0.85724617,  0.86663356,  0.87618913,  0.88590785,  0.89578486,  0.90581541,
                       0.91599489,  0.9263188 ,  0.93678277,  0.94738257,  0.95811405,  0.9689732,
                       0.9689732 ])



   
    # interpolation function for fitting data
    mfn_line_array_fit = MFnLineArray()
    mfn_line_array_fit.set( xdata = xdata_fit, ydata = ydata_fit )
    mfn_line_array_fit.data_changed = True
    
    
    
    # get the current strain:
    eps_app_tn1 = copy( fitter.tloop.U_k )
    print 'eps_app_tn1', eps_app_tn1    
    # get the current sctx:
    sctx = fitter.tloop.tstepper.sctx
    print 'sctx', sctx    

    # for fitting use default method 'get_value' of 'MFnLineArray' as 'phi_fn':
    fitter.mats2D_eval.polar_fn.phi_fn.mfn.get_value = MFnLineArray().get_value

 def pdf(self, e, r, L):
     cdf_line = MFnLineArray(xdata=e, ydata=self.cdf(e, r, L))
     return cdf_line.get_diffs(e)

from mathkit.mfn.mfn_line.mfn_line import MFnLineArray
from numpy import linspace, frompyfunc, array
from math import sin, cos

sigma_max = 0.4

xdata = linspace( 0., 1., 10000 )
fnydata = lambda x:  1.0+0.0001*sin(19.0*x)*cos(200*x) * sigma_max * 10000

fnydata_vec = frompyfunc( fnydata, 1, 1 )
ydata = array( fnydata_vec( xdata ), dtype = float )

#print ydata

mfn = MFnLineArray( xdata = xdata, ydata = ydata )
mfn.configure_traits()

 def _get_matrix_strength(self):
     # evaluates a random field
     # realization and creates a spline representation
     rf = self.random_field.random_field
     rf_spline = MFnLineArray(xdata=self.random_field.xgrid, ydata=rf)
     return rf_spline.get_values(self.x_arr)

    def Pw(self, w):
        mask = self.cached_pdfs[3] < w
        dx2 = (self.x[-1]-self.x[-2])**2
        return np.sum(self.cached_pdfs[2] * mask * dx2)
    
    def CDF_w(self, w_arr):
        cdfs = []
        for w in w_arr:
            cdfs.append(self.Pw(w))
        return cdfs

if __name__ == '__main__':
    ac = AnalyticalCracks(l0=1., d=0.007, tau=0.1, sigma0=2200., s=2.0,
                          rho=5.0, c=1.0, x=np.linspace(0.0, 5.0, 500))
    for s in [1., 3., 5.0]:
        ac.s = s
        pdf_x = ac.p_x(s, ac.x)
        pdf_x = pdf_x / np.trapz(pdf_x, ac.x)
        cdf_x = np.hstack((0., cumtrapz(pdf_x, ac.x)))
        sf = MFnLineArray(xdata=cdf_x + 1e-12 * ac.x, ydata=ac.x)
        rand_vals1 = sf.get_values(np.random.rand(10000))
        rand_vals2 = sf.get_values(np.random.rand(10000))
        cracks = ac.w_func(rand_vals1, rand_vals2)
        plt.hist(cracks, bins=40, normed=True, label=str(s), cumulative=True)
        w_arr = np.linspace(0.0,7.0,100)
        cdf_w = ac.CDF_w(w_arr)
        plt.plot(w_arr, cdf_w, color='black', lw=2)
    plt.ylim(0,1.2)
    plt.legend()
    plt.show()