Python scipy 模块，arange() 实例源码

def _plotFMeasures(fstepsize=.1,  stepsize=0.0005, start = 0.0, end = 1.0):
    """Plots 10 fmeasure Curves into the current canvas."""
    p = sc.arange(start, end, stepsize)[1:]
    for f in sc.arange(0., 1., fstepsize)[1:]:
        points = [(x, _fmeasureCurve(f, x)) for x in p
                  if 0 < _fmeasureCurve(f, x) <= 1.5]
        try:
            xs, ys = zip(*points)
            curve, = pl.plot(xs, ys, "--", color="gray", linewidth=0.8)  # , label=r"$f=%.1f$"%f) # exclude labels, for legend
            # bad hack:
            # gets the 10th last datapoint, from that goes a bit to the left, and a bit down
            datapoint_x_loc = int(len(xs)/2)
            datapoint_y_loc = int(len(ys)/2)
            # x_left = 0.05
            # y_left = 0.035
            x_left = 0.035
            y_left = -0.02
            pl.annotate(r"$f=%.1f$" % f, xy=(xs[datapoint_x_loc], ys[datapoint_y_loc]), xytext=(xs[datapoint_x_loc] - x_left, ys[datapoint_y_loc] - y_left), size="small", color="gray")
        except Exception as e:
            print e 

#colors = "gcmbbbrrryk"
#colors = "yyybbbrrrckgm"  # 7 is a prime, so we'll loop over all combinations of colors and markers, when zipping their cycles

def bowtie_py(modis_img,p,cs):
    print 'removie bowtie effect from image... '
    stripwidth=10000/cs
    #Loop over every x coordinate of the image
    for x in sp.arange(modis_img.shape[1]):
        #Loop over every sanning strip
        overlap=sp.polyval(p,x).round() #get the overlap from the polynom
        if overlap > 0:
            for y in sp.arange(stripwidth,modis_img.shape[0],stripwidth):
                #cut out the current part:
                strippart=modis_img[y:y+stripwidth,x]
                #delete the upper and lower few pixels of the strippart:
                strippart=strippart[int(overlap/2.):-int(round(overlap/2.))]
                #Interpolat to stipwidth length
                f=interp1d(sp.arange(0,len(strippart)),strippart)
                strippart=f(sp.linspace(0,len(strippart)-1,stripwidth))
                #replace the current sick part in the image by the new healthy one
                modis_img[y:y+stripwidth,x]=strippart
    print 'done'
    return modis_img

def KMV_f(E,D,T,r,sigmaE):
    n=10000
    m=2000
    diffOld=1e6     # a very big number
    for i in sp.arange(1,10):
        for j in sp.arange(1,m):
            A=E+D/2+i*D/n
            sigmaA=0.05+j*(1.0-0.001)/m
            d1 = (log(A/D)+(r+sigmaA*sigmaA/2.)*T)/(sigmaA*sqrt(T))
            d2 = d1-sigmaA*sqrt(T)
            diff4E= (A*N(d1)-D*exp(-r*T)*N(d2)-E)/A  # scale by assets
            diff4A= A/E*N(d1)*sigmaA-sigmaE          # a small number already
            diffNew=abs(diff4E)+abs(diff4A)
            if diffNew<diffOld:
               diffOld=diffNew
               output=(round(A,2),round(sigmaA,4),round(diffNew,5))
    return output
#

def up_and_out_call(s0,x,T,r,sigma,n_simulation,barrier):
    n_steps=100. 
    dt=T/n_steps 
    total=0 
    for j in sp.arange(0, n_simulation): 
        sT=s0 
        out=False
        for i in range(0,int(n_steps)): 
            e=sp.random.normal() 
            sT*=sp.exp((r-0.5*sigma*sigma)*dt+sigma*e*sp.sqrt(dt)) 
            if sT>barrier: 
               out=True 
        if out==False: 
            total+=bsCall(s0,x,T,r,sigma) 
    return total/n_simulation 
#

def up_and_out_call(s0,x,T,r,sigma,n_simulation,barrier):
    n_steps=100. 
    dt=T/n_steps 
    total=0 
    for j in sp.arange(0, n_simulation): 
        sT=s0 
        out=False
        for i in range(0,int(n_steps)): 
            e=sp.random.normal() 
            sT*=sp.exp((r-0.5*sigma*sigma)*dt+sigma*e*sp.sqrt(dt)) 
            if sT>barrier: 
               out=True 
        if out==False: 
            total+=bsCall(s0,x,T,r,sigma) 
    return total/n_simulation 
#

def alt_results(self, samples, kplanets):
        titles = sp.array(["Amplitude","Period","Longitude", "Phase","Eccentricity", 'Acceleration', 'Jitter', 'Offset', 'MACoefficient', 'MATimescale', 'Stellar Activity'])
        namen = sp.array([])
        ndim = kplanets * 5 + self.nins*2*(self.MOAV+1) + self.totcornum + 1

        RESU = sp.zeros((ndim, 5))
        for k in range(kplanets):
            namen = sp.append(namen, [titles[i] + '_'+str(k) for i in range(5)])
        namen = sp.append(namen, titles[5])  # for acc
        for i in range(self.nins):
            namen = sp.append(namen, [titles[ii] + '_'+str(i+1) for ii in sp.arange(2)+6])
            for c in range(self.MOAV):
                namen = sp.append(namen, [titles[ii] + '_'+str(i+1) + '_'+str(c+1) for ii in sp.arange(2)+8])
        for h in range(self.totcornum):
            namen = sp.append(namen, titles[-1]+'_'+str(h+1))

        alt_res = map(lambda v: (v[2], v[3]-v[2], v[2]-v[1], v[4]-v[2], v[2]-v[0]),
                      zip(*np.percentile(samples, [2, 16, 50, 84, 98], axis=0)))
        logdat = '\nAlternative results with uncertainties based on the 2nd, 16th, 50th, 84th and 98th percentiles of the samples in the marginalized distributions'
        logdat = '\nFormat is like median +- 1-sigma, +- 2-sigma'
        for res in range(ndim):
            logdat += '\n'+namen[res]+'     : '+str(alt_res[res][0])+' +- '+str(alt_res[res][1:3]) +'    2%   +- '+str(alt_res[res][3:5])
            RESU[res] = sp.percentile(samples, [2, 16, 50, 84, 98], axis=0)[:, res]
        print(logdat)
        return RESU

def test_fetch_one_column():
    _urlopen_ref = datasets.mldata.urlopen
    try:
        dataname = 'onecol'
        # create fake data set in cache
        x = sp.arange(6).reshape(2, 3)
        datasets.mldata.urlopen = mock_mldata_urlopen({dataname: {'x': x}})

        dset = fetch_mldata(dataname, data_home=tmpdir)
        for n in ["COL_NAMES", "DESCR", "data"]:
            assert_in(n, dset)
        assert_not_in("target", dset)

        assert_equal(dset.data.shape, (2, 3))
        assert_array_equal(dset.data, x)

        # transposing the data array
        dset = fetch_mldata(dataname, transpose_data=False, data_home=tmpdir)
        assert_equal(dset.data.shape, (3, 2))
    finally:
        datasets.mldata.urlopen = _urlopen_ref

def kde(data, N=None, MIN=None, MAX=None):

    # Parameters to set up the mesh on which to calculate
    N = 2**14 if N is None else int(2**sci.ceil(sci.log2(N)))
    if MIN is None or MAX is None:
        minimum = min(data)
        maximum = max(data)
        Range = maximum - minimum
        MIN = minimum - Range/10 if MIN is None else MIN
        MAX = maximum + Range/10 if MAX is None else MAX

    # Range of the data
    R = MAX-MIN

    # Histogram the data to get a crude first approximation of the density
    M = len(data)
    DataHist, bins = sci.histogram(data, bins=N, range=(MIN,MAX))
    DataHist = DataHist/M
    DCTData = scipy.fftpack.dct(DataHist, norm=None)

    I = [iN*iN for iN in xrange(1, N)]
    SqDCTData = (DCTData[1:]/2)**2

    # The fixed point calculation finds the bandwidth = t_star
    guess = 0.1
    try:
        t_star = scipy.optimize.brentq(fixed_point, 0, guess, 
                                       args=(M, I, SqDCTData))
    except ValueError:
        print 'Oops!'
        return None

    # Smooth the DCTransformed data using t_star
    SmDCTData = DCTData*sci.exp(-sci.arange(N)**2*sci.pi**2*t_star/2)
    # Inverse DCT to get density
    density = scipy.fftpack.idct(SmDCTData, norm=None)*N/R
    mesh = [(bins[i]+bins[i+1])/2 for i in xrange(N)]
    bandwidth = sci.sqrt(t_star)*R

    density = density/sci.trapz(density, mesh)
    return bandwidth, mesh, density

def iterate_minibatches(self,X, y, batchsize, shuffle=False):
        assert len(X) == len(y)
        if shuffle:
            indices = sp.arange(len(X))
            sp.random.shuffle(indices)
        for start_idx in range(0, len(X) - batchsize + 1, batchsize):
            if shuffle:
                excerpt = indices[start_idx:start_idx + batchsize]
            else:
                excerpt = slice(start_idx, start_idx + batchsize)
            yield X[excerpt], y[excerpt]

def fit(self, X):
        n_samples, n_features = X.shape
        n_classes = self.n_classes
        max_iter = self.max_iter
        tol = self.tol

        rand_center_idx = sprand.permutation(n_samples)[0:n_classes]
        center = X[rand_center_idx].T
        responsilibity = sp.zeros((n_samples, n_classes))

        for iter in range(max_iter):
            # E step
            dist = sp.expand_dims(X, axis=2) - sp.expand_dims(center, axis=0)
            dist = spla.norm(dist, axis=1)**2
            min_idx = sp.argmin(dist, axis=1)
            responsilibity.fill(0)
            responsilibity[sp.arange(n_samples), min_idx] = 1

            # M step
            center_new = sp.dot(X.T, responsilibity) / sp.sum(responsilibity, axis=0)
            diff = center_new - center
            print('K-Means: {0:5d} {1:4e}'.format(iter, spla.norm(diff) / spla.norm(center)))
            if (spla.norm(diff) < tol * spla.norm(center)):
                break

            center = center_new

        self.center = center.T
        self.responsibility = responsilibity

        return self

def bowtie_polynom(modis_img,cs,folder):
    print 'Determine overlap pattern... '
    sw=10000/cs #stripwidth
    overlaplist=[]#define list to store number of overlapped lines
    #devide in parts with a width of 40 pixel
    for i in sp.arange(0,modis_img.shape[1]-40,40):
        part=modis_img[:,i:i+39]
        #search in every scanning strip
        samples=[]
        for j in sp.arange(sw-2,part.shape[0]-sw,sw):
            target=part[j-1:j+1,:] #cut out a target, which overlapped counter-part shall be found
            searchwindow=part[j+2:j+sw+2] #,: cut out the window, where the overlapped counter part might be located
            #start the search
            c=[] #calculate correlation coefficients of every given offset from 3 to 11
            for offset in sp.arange(3,sw/2+1):
                imgpart=searchwindow[offset-3:offset-1] #,: cut out image, which has to be compared with the target
                c.append(sp.corrcoef(imgpart.flatten(),target.flatten())[0,1])#calculate correlatoin coefficient
            c=sp.array(c)
            overl=sp.ndimage.measurements.maximum_position(c)[0]+3 #find the overlap with the highes correlation coefficient
            samples.append([overl,c.max()]) #attach overlap and correlation coefficient to the sample list
        samples=sp.array(samples)
        #print i, samples[:,1].mean()
        if samples[:,1].mean() > 0.9: #chek the mean correlation coefficient:
            #print('Bowtie Correlation high - removing effect')
            overlaplist.append([i+20,samples[:,0].mean()]) #save result, if correlation coefficient is high
            #print(overlaplist)
            o=sp.array(overlaplist)
            X=o[:,0]
            overlap=o[:,1]
            #Calculate a second order Polynom to describe the overlap
            p=sp.polyfit(X,overlap,2)
            #print 'done, Overlap polynom: '+str(p)
        else:
            #print('low Bowtie correlation')
            p = [1.,  1.,  1.]
            #overlaplist.append([i+20,1])
            #os.system('rm -r '+folder)
            #print('scene deleted') 
    return p

def LPSD_f(returns, Rf):
    y=returns[returns-Rf<0]  
    m=len(y)
    total=0.0
    for i in sp.arange(m):
        total+=(y[i]-Rf)**2
    return total/(m-1)

def LPSD_f(returns, Rf):
    y=returns[returns-Rf<0]  
    m=len(y)
    total=0.0
    for i in sp.arange(m):
        total+=(y[i]-Rf)**2
    return total/(m-1)

# function 3: estimate Sortino

def BIS_f(R,s,n):
    R=R0
    for i in sp.arange(0,n):
        deltaR=z[i]*s/sp.sqrt(2.)
        logR=sp.log(R)
        R=sp.exp(logR+deltaR)
        output.append(round(R,5))
    return output 
#

def binomialCall(s,x,T,r,sigma,n=100): 
    deltaT = T /n
    u = exp(sigma * sqrt(deltaT)) 
    d = 1.0 / u
    a = exp(r * deltaT)
    p = (a - d) / (u - d)
    v = [[0.0 for j in sp.arange(i + 1)]  for i in sp.arange(n + 1)] 
    for j in sp.arange(n+1):
        v[n][j] = max(s * u**j * d**(n - j) - x, 0.0) 
    for i in sp.arange(n-1, -1, -1):
        for j in sp.arange(i + 1):
            v[i][j]=exp(-r*deltaT)*(p*v[i+1][j+1]+(1.0-p)*v[i+1][j]) 
    return v[0][0]

def predict(self,xt,tau=None,confidenceMap=None):
        '''
        Function that predict the label for sample xt using the learned model
        Inputs:
            xt: the samples to be classified
        Outputs:
            y: the class
            K: the decision value for each class
        '''

        MAX = sp.finfo(sp.float64).max
        E_MAX = sp.log(MAX) # Maximum value that is possible to compute with sp.exp

        ## Get information from the data
        nt = xt.shape[0]        # Number of testing samples
        C = self.ni.shape[0]    # Number of classes

        ## Initialization
        K = sp.empty((nt,C))

        if tau is None:
            TAU=self.tau
        else:
            TAU=tau

        for c in range(C):
            invCov,logdet = self.compute_inverse_logdet(c,TAU)
            cst = logdet - 2*sp.log(self.prop[c]) # Pre compute the constant

            xtc = xt-self.mean[c,:]
            temp = sp.dot(invCov,xtc.T).T
            K[:,c] = sp.sum(xtc*temp,axis=1)+cst
            del temp,xtc        

        yp = sp.argmin(K,1)

        if confidenceMap is None:

            ## Assign the label save in classnum to the minimum value of K 
            yp = self.classnum[yp]

            return yp

        else:

            K *= -0.5
            K[K>E_MAX],K[K<-E_MAX] = E_MAX,-E_MAX
            sp.exp(K,out=K)
            K /= K.sum(axis=1).reshape(nt,1)
            K = K[sp.arange(len(K)),yp]
            #K = sp.diag(K[:,yp])

            yp = self.classnum[yp]

            return yp,K

def pt_pos(self, kplanets, boundaries, inslims, acc_lims):
        ndim = 1 + 5 * kplanets + self.nins*2*(self.MOAV+1) + self.totcornum
        pos = sp.array([sp.zeros(ndim) for i in range(self.nwalkers)])
        k = -2
        l = -2
        ll = -2  ##
        for j in range(ndim):
            if j < 5 * kplanets:
                k += 2
                if j%5==1:
                    fact = sp.absolute(boundaries[k] - boundaries[k+1]) / self.nwalkers
                else:
                    #fact = sp.absolute(boundaries[k]) / (self.nwalkers)
                    fact = (sp.absolute(boundaries[k] - boundaries[k+1]) * 2) / (5 * self.nwalkers)
                dif = sp.arange(self.nwalkers) * fact * np.random.uniform(0.9, 0.999)
                for i in range(self.nwalkers):
                    if j%5==1:
                        pos[i][j] = boundaries[k] + (dif[i] + fact/2.0)
                    else:
                        #pos[i][j] = boundaries[k] * 0.5 + (dif[i] + fact/2.0)
                        pos[i][j] = (boundaries[k+1]+3*boundaries[k])/4 + (dif[i] + fact/2.0)
            if j == 5 * kplanets:  # acc
                fact = sp.absolute(acc_lims[0] - acc_lims[1]) / self.nwalkers
                dif = sp.arange(self.nwalkers) * fact * np.random.uniform(0.9, 0.999)
                for i in range(self.nwalkers):
                    pos[i][j] = acc_lims[0] + (dif[i] + fact/2.0)

            if 5 * kplanets < j < 5*kplanets + self.nins*2*(self.MOAV+1) + 1:
                l += 2
                fact = sp.absolute(inslims[l] - inslims[l+1]) / self.nwalkers
                dif = sp.arange(self.nwalkers) * fact * np.random.uniform(0.9, 0.999)

                if (j-5*kplanets-1) % self.nins*2*(self.MOAV+1) == 0:  # ojo aqui
                    jitt_ini = sp.sort(sp.fabs(sp.random.normal(0, 1, self.nwalkers))) * 0.1
                    dif = jitt_ini * np.random.uniform(0.9, 0.999)
                for i in range(self.nwalkers):
                    pos[i][j] = inslims[l] + (dif[i] + fact/2.0)
            if self.totcornum:
                if j > 5*kplanets + self.nins*2*(self.MOAV+1):
                    fact = sp.absolute(acc_lims[0] - acc_lims[1]) / self.nwalkers

                    dif = sp.arange(self.nwalkers) * fact * np.random.uniform(0.8, 0.999)
                    for i in range(self.nwalkers):
                        pos[i][j] = acc_lims[0] + (dif[i] + fact/2.0)
                        #print(pos[i][j])

        pos = sp.array([pos for h in range(self.ntemps)])
        return pos