schain Commit - r1683:256acf43d60a · Jicamarca Repository

23

from numpy import NaN

23

from numpy import NaN

24

from scipy.optimize.optimize import OptimizeWarning

24

from scipy.optimize.optimize import OptimizeWarning

25

warnings.filterwarnings('ignore')

25

warnings.filterwarnings('ignore')

26

import pdb

26

27

28

SPEED_OF_LIGHT = 299792458

29

SPEED_OF_LIGHT = 299792458

2732

dataOut.nIncohInt= int(dataOut.nIncohInt)

2733

dataOut.nIncohInt= int(dataOut.nIncohInt)

2733

dataOut.nProfiles = int(dataOut.nProfiles)

2734

dataOut.nProfiles = int(dataOut.nProfiles)

2734

dataOut.nCohInt = int(dataOut.nCohInt)

2735

dataOut.nCohInt = int(dataOut.nCohInt)

2735

print('sale',dataOut.nProfiles,dataOut.nHeights)

2736

#dataOut.nFFTPoints=nprofs

2736

#dataOut.nFFTPoints=nprofs

2737

#dataOut.normFactor = nprofs

2737

#dataOut.normFactor = nprofs

2738

dataOut.channelList = channelList

2738

dataOut.channelList = channelList

2741

vmax = (300000000/49920000.0/2) / (dataOut.ippSeconds)

2741

vmax = (300000000/49920000.0/2) / (dataOut.ippSeconds)

2742

#dataOut.ippSeconds=ipp

2742

#dataOut.ippSeconds=ipp

2743

absc = vmax*( numpy.arange(nProf,dtype='float')-nProf/2.)/nProf

2743

absc = vmax*( numpy.arange(nProf,dtype='float')-nProf/2.)/nProf

2744

print('sale 2',dataOut.ippSeconds,M,N)

2745

print('Empieza procesamiento offline')

2746

if path != None:

2744

if path != None:

2747

sys.path.append(path)

2745

sys.path.append(path)

2748

self.library = importlib.import_module(file)

2746

self.library = importlib.import_module(file)

2786

jvelr = numpy.zeros(nHeights, dtype = 'float')

2784

jvelr = numpy.zeros(nHeights, dtype = 'float')

2787

hvalid = [0]

2785

hvalid = [0]

2788

coh2 = abs(dataOut.data_cspc[i,1:nProf,:])**2/(dataOut.data_spc[0+i*2,1:nProf-0,:]*dataOut.data_spc[1+i*2,1:nProf-0,:])

2786

coh2 = abs(dataOut.data_cspc[i,1:nProf,:])**2/(dataOut.data_spc[0+i*2,1:nProf-0,:]*dataOut.data_spc[1+i*2,1:nProf-0,:])

2787

2789

for h in range(nHeights):

2788

for h in range(nHeights):

2790

smooth = clean_num_aver[i+1,h]

2789

smooth = clean_num_aver[i+1,h]

2791

signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth

2790

signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth

2821

snr0 = numpy.sum(signal0/n0)/(nProf-1)

2820

snr0 = numpy.sum(signal0/n0)/(nProf-1)

2822

snr1 = numpy.sum(signal1/n1)/(nProf-1)

2821

snr1 = numpy.sum(signal1/n1)/(nProf-1)

2823

if snr0 > snrth and snr1 > snrth and clean_num_aver[i+1,h] > 0 :

2822

if snr0 > snrth and snr1 > snrth and clean_num_aver[i+1,h] > 0 :

2824

#Covariance Matrix

2823

#Covariance Matrix

2825

D = numpy.diag(d**2)

2824

D = numpy.diag(d**2)

2826

ind = 0

2825

ind = 0

2827

for pairs in listComb:

2826

for pairs in listComb:

2846

LT=L.T

2845

LT=L.T

2847

2846

2848

dp = numpy.dot(LT,d)

2847

dp = numpy.dot(LT,d)

2849

#Initial values

2848

#Initial values

2850

data_spc = dataOut.data_spc[coord,:,h]

2849

data_spc = dataOut.data_spc[coord,:,h]

2851

w = data_spc/data_spc

2850

w = data_spc/data_spc

2852

if filec != None:

2851

if filec != None:

2853

w = self.weightf.weightfit(w,tini.tm_year,tini.tm_yday,index,h,i)

2852

w = self.weightf.weightfit(w,tini.tm_year,tini.tm_yday,index,h,i)

2854

if (h>6)and(error1[3]<25):

2853

if (h>6)and(error1[3]<25):

2855

p0 = dataOut.data_param[i,:,h-1]

2854

p0 = dataOut.data_param[i,:,h-1].copy()

2856

else:

2855

else:

2857

p0 = numpy.array(self.library.initialValuesFunction(data_spc*w, constants))# sin el i(data_spc, constants, i)

2856

p0 = numpy.array(self.library.initialValuesFunction(data_spc*w, constants))# sin el i(data_spc, constants, i)

2858

#print("AQUI REEMPLAZAS CON EL fd0",fd0)

2859

p0[3] = fd0

2857

p0[3] = fd0

2858

2860

if filec != None:

2859

if filec != None:

2861

p0 = self.weightf.Vrfit(p0,tini.tm_year,tini.tm_yday,index,h,i)

2860

p0 = self.weightf.Vrfit(p0,tini.tm_year,tini.tm_yday,index,h,i)

2862

2861

2863

#print("minP0-ANTES DE OPTIMIZE")

2864

#print(p0)

2865

try:

2862

try:

2866

#Least Squares

2863

#Least Squares

2867

minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)

2864

minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)

2868

print("VERIFICAR SI ESTA VARIANDO minp--------------------------------", minp[3])

2869

print("COMPARACION ---- p0[3],fd0",p0[3],fd0)

2870

#minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))

2865

#minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))

2871

#Chi square error

2866

#Chi square error

2872

error0 = numpy.sum(infodict['fvec']**2)/(2*N)

2867

error0 = numpy.sum(infodict['fvec']**2)/(2*N)

2882

p0 = numpy.array(self.library.initialValuesFunction(data_spc, constants))

2877

p0 = numpy.array(self.library.initialValuesFunction(data_spc, constants))

2883

minp = p0*numpy.nan

2878

minp = p0*numpy.nan

2884

error0 = numpy.nan

2879

error0 = numpy.nan

2885

error1 = p0*numpy.nan

2880

error1 = p0*numpy.nan

2886

if dataOut.data_param is None:

2881

if dataOut.data_param is None:

2887

dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan

2882

dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan

2888

dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan

2883

dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan

2889

print("CLASE SF minp- ?????",minp)

2890

dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))

2884

dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))

2891

dataOut.data_param[i,:,h] = minp

2885

dataOut.data_param[i,:,h] = minp

2886

2892

for ht in range(nHeights-1) :

2887

for ht in range(nHeights-1) :

2893

smooth = coh_num_aver[i+1,ht] #datc[0,ht,0,beam]

2888

smooth = coh_num_aver[i+1,ht] #datc[0,ht,0,beam]

2894

dataOut.data_paramC[4*i,ht,1] = smooth

2889

dataOut.data_paramC[4*i,ht,1] = smooth

2927

mom1 = self.moments(doppler,signalpn1_n1,nProf)

2922

mom1 = self.moments(doppler,signalpn1_n1,nProf)

2928

dataOut.data_paramC[2+4*i,ht,0] = (mom0[0]+mom1[0])/2.

2923

dataOut.data_paramC[2+4*i,ht,0] = (mom0[0]+mom1[0])/2.

2929

dataOut.data_paramC[3+4*i,ht,0] = (mom0[1]+mom1[1])/2.

2924

dataOut.data_paramC[3+4*i,ht,0] = (mom0[1]+mom1[1])/2.

2930

2931

dataOut.data_spc = jspectra

2925

dataOut.data_spc = jspectra

2932

dataOut.spc_noise = my_noises*nProf*M

2926

dataOut.spc_noise = my_noises*nProf*M

2933

if numpy.any(proc): dataOut.spc_noise = my_noises*nProf*M

2927

if numpy.any(proc): dataOut.spc_noise = my_noises*nProf*M

3654

else :

3648

else :

3655

moments=numpy.vstack(([velRadialm[0,:]],[velRadialm[0,:]]))

3649

moments=numpy.vstack(([velRadialm[0,:]],[velRadialm[0,:]]))

3656

dataOut.moments=moments

3650

dataOut.moments=moments

3657

print("CLASE EWD Estimation",moments)

3658

# Coherent

3651

# Coherent

3659

smooth_wC = ebufc[0,:]

3652

smooth_wC = ebufc[0,:]

3660

p_w0C = rbufc[0,:]

3653

p_w0C = rbufc[0,:]

	Site-wide shortcuts
/	Use quick search box
g h	Goto home page
g g	Goto my private gists page
g G	Goto my public gists page
g 0-9	Goto bookmarked items from 0-9
n r	New repository page
n g	New gist page

	Repositories
g s	Goto summary page
g c	Goto changelog page
g f	Goto files page
g F	Goto files page with file search activated
g p	Goto pull requests page
g o	Goto repository settings
g O	Goto repository access permissions settings
t s	Toggle sidebar on some pages

             from numpy import NaN
             from scipy.optimize.optimize import OptimizeWarning
             warnings.filterwarnings('ignore')
+            import pdb
             SPEED_OF_LIGHT = 299792458
                         dataOut.nIncohInt= int(dataOut.nIncohInt)
                         dataOut.nProfiles = int(dataOut.nProfiles)
                         dataOut.nCohInt = int(dataOut.nCohInt)
-                        print('sale',dataOut.nProfiles,dataOut.nHeights)
                         #dataOut.nFFTPoints=nprofs
                         #dataOut.normFactor = nprofs
                         dataOut.channelList = channelList
                         vmax = (300000000/49920000.0/2) / (dataOut.ippSeconds)
                         #dataOut.ippSeconds=ipp
                         absc = vmax*( numpy.arange(nProf,dtype='float')-nProf/2.)/nProf
-                        print('sale 2',dataOut.ippSeconds,M,N)
-                        print('Empieza procesamiento offline')
                         if path != None:
                             sys.path.append(path)
                         self.library = importlib.import_module(file)
                         jvelr = numpy.zeros(nHeights, dtype = 'float')
                         hvalid = [0]
                         coh2 = abs(dataOut.data_cspc[i,1:nProf,:])**2/(dataOut.data_spc[0+i*2,1:nProf-0,:]*dataOut.data_spc[1+i*2,1:nProf-0,:])
                         for h in range(nHeights):
                             smooth = clean_num_aver[i+1,h]
                             signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth
                             snr0 = numpy.sum(signal0/n0)/(nProf-1)
                             snr1 = numpy.sum(signal1/n1)/(nProf-1)
                             if snr0 > snrth and snr1 > snrth and clean_num_aver[i+1,h] > 0 :
-                            #Covariance Matrix
+                                #Covariance Matrix
                                 D = numpy.diag(d**2)
                                 ind = 0
                                 for pairs in listComb:
                                 LT=L.T
                                 dp = numpy.dot(LT,d)
-                            #Initial values
+                                #Initial values
                                 data_spc = dataOut.data_spc[coord,:,h]
                                 w = data_spc/data_spc
                                 if filec != None:
                                    w = self.weightf.weightfit(w,tini.tm_year,tini.tm_yday,index,h,i)
                                 if (h>6)and(error1[3]<25):
-                                    p0 = dataOut.data_param[i,:,h-1]
+                                    p0 = dataOut.data_param[i,:,h-1].copy()
                                 else:
                                     p0 = numpy.array(self.library.initialValuesFunction(data_spc*w, constants))# sin el i(data_spc, constants, i)
-                                #print("AQUI REEMPLAZAS CON EL fd0",fd0)
                                 p0[3] = fd0
                                 if filec != None:
                                    p0 = self.weightf.Vrfit(p0,tini.tm_year,tini.tm_yday,index,h,i)
-                                #print("minP0-ANTES DE OPTIMIZE")
-                                #print(p0)
                                 try:
                                     #Least Squares
                                     minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
-                                    print("VERIFICAR SI ESTA VARIANDO minp--------------------------------", minp[3])
-                                    print("COMPARACION ---- p0[3],fd0",p0[3],fd0)
                                     #minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
                                     #Chi square error
                                     error0 = numpy.sum(infodict['fvec']**2)/(2*N)
                                 p0 = numpy.array(self.library.initialValuesFunction(data_spc, constants))
                                 minp = p0*numpy.nan
                                 error0 = numpy.nan
                                 error1 = p0*numpy.nan
                             if dataOut.data_param is None:
                                 dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
                                 dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
-                            print("CLASE  SF minp- ?????",minp)
                             dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
                             dataOut.data_param[i,:,h] = minp
                         for ht in range(nHeights-1) :
                             smooth = coh_num_aver[i+1,ht] #datc[0,ht,0,beam]
                             dataOut.data_paramC[4*i,ht,1] = smooth
                                 mom1 = self.moments(doppler,signalpn1_n1,nProf)
                                 dataOut.data_paramC[2+4*i,ht,0] = (mom0[0]+mom1[0])/2.
                                 dataOut.data_paramC[3+4*i,ht,0] = (mom0[1]+mom1[1])/2.
                     dataOut.data_spc = jspectra
                     dataOut.spc_noise = my_noises*nProf*M
                     if numpy.any(proc): dataOut.spc_noise = my_noises*nProf*M
                     else :
                        moments=numpy.vstack(([velRadialm[0,:]],[velRadialm[0,:]]))
                     dataOut.moments=moments
-                    print("CLASE EWD Estimation",moments)
                     # Coherent
                     smooth_wC = ebufc[0,:]
                     p_w0C = rbufc[0,:]