exit.py－20170818

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Aug 14 11:01:48 2017

@author: vicky
"""

from numpy import *
import numpy as np
import pandas as pd
import math 
import matplotlib.pyplot as plt

[n,p]=np.shape(y); #输入y=原始数据，t=时间点

#波峰识别算法
#------------------------平滑处理-------------------------
#f＝y的四个区间累计
f=np.zeros((n,p), dtype=np.int)
for i in range(n-3): 
    for j in range (0,p):
        f[i+3,j]=y[i,j]+y[i+1,j]+y[i+2,j]+y[i+3,j]

#-----------------------取f的二阶差分----------------------
#f的一,二，三阶差分df1,df2
df1=np.zeros((n,p), dtype=np.int)
df2=np.zeros((n,p), dtype=np.int)
#df3=np.zeros((n,p), dtype=np.int)
for j in range (p):
    for i in range (n-1):
        df1[i+1,j]=f[i+1,j]-f[i,j]
    for i in range (n-2):
        df2[i+2,j]=f[i+2,j]-2*f[i+1,j]+f[i,j]

#-----------------------所有波峰判断------------------------
#所有波峰fmax=f的极大值且df2的局部极小值
fmax_all=np.zeros((n,p), dtype=np.int)
for j in range(p):
    for i in range(2,n-2):
        if f[i,j]>=max(f[i-1,j],f[i+1,j]): #and 
#            (df2[i,j]<=min(df2[i-1,j],df2[i+1,j]) or 
#             df2[i+1,j]<=min(df2[i,j],df2[i+2,j]) or 
#             df2[i-1,j]<=min(df2[i-2,j],df2[i,j]))):
            fmax_all[i,j]=f[i,j]

#------------------------时间段处理-------------------------
#ffmax=fmax去0，且8点以后,且>100;ttmax=ffmax对应的时间点
ffmax_all=np.zeros((n,p), dtype=np.int)
ttmax_all=np.zeros((n,p), dtype=np.int)
for j in range(0,p):  
    k=0
    for i in range(97,n):
        if (fmax_all[i,j]!=0 and fmax_all[i,j]>100):
            ffmax_all[k,j]=fmax_all[i,j]
            ttmax_all[k,j]=i
            k=k+1

ffmax_all=ffmax_all[~np.all(ffmax_all == 0, axis=1)]
ttmax_all=ttmax_all[~np.all(ttmax_all == 0, axis=1)]

#--------------------------去除假波峰-----------------------
[n1,p1]=np.shape(ffmax_all)
rffmax_all=np.zeros(ffmax_all.shape) #rffmax＝波峰的增长率
for j in range(0,p1):
    for i in range(0,n1):
        if ffmax_all[i,j]!=0 and ffmax_all[i-1,j]!=0 :
            rffmax_all[i,j]=(ffmax_all[i,j]-ffmax_all[i-1,j])/ffmax_all[i-1,j]

#rffmax[1,:]=0 #从第2个波峰起开始判断
#fffmax=去除rffmax大于2.223的假波峰，起点波峰之后的波峰,tttmax=对应时间点
fffmax_all=np.zeros((n,p), dtype=np.int)
tttmax_all=np.zeros((n,p), dtype=np.int)
for j in range(0,p):
    tmp=np.array(np.where(rffmax_all[:,j]>2.223))
    k=0
    if size(tmp)>0:
        for i in range(tmp[0,0],n1):
            fffmax_all[k,j]=ffmax_all[i,j]
            tttmax_all[k,j]=ttmax_all[i,j]
            k=k+1
    else:
        fffmax_all[:,j]=0
        tttmax_all[:,j]=0
fffmax_all=fffmax_all[~np.all(fffmax_all==0,axis=1)]
tttmax_all=tttmax_all[~np.all(tttmax_all==0,axis=1)]        
[n2,p2]=np.shape(fffmax_all)

#---------------------------图例----------------------------
#td=range(0,288)
#for j in range(p):
#    fig = plt.plot(td,f[:,j])
##   plt.savefig('/Users/vicky/Documents/地图/演唱会/pythonfig/'+np.str(j)+'.jpg')
#    plt.show()
  
#plt.plot(df1[:,0])         
#plt.plot(df2[:,0]) 

#-------------------------散场斜率---------------------------
#te_beg=zeros((1,p),dtype=int) #te_beg=末端波峰位置，即看作散场开始时间
#kt=zeros((1,p),dtype=int) 
#for j in range(p):
#    for i in range(n1):
#        if rffmax_all[i,j]>2.223:
#            kt[0,j]=i
#    for i in range(kt[0,j],n1):
#        if rffmax_all[i,j]<-0.9:
#            te_beg[0,j]=ttmax_all[i-1,j]

##单位根检验
#import statsmodels.tsa.stattools as ts
#for j in range(5):
#    result = ts.adfuller(df1[:,j])
#    print(j)
#    if result[0]<result[4]['1%']:
#       print ('没有单位根')    
#    else:
#       print('有单位根')

te_end=zeros((1,p),dtype=int) #te_end=散场结束时间，前后一阶查分平稳
for j in range(p):
    if tttmax_all[0,j]!=0:
      for i in range(tttmax_all[0,j],n):
        if max(abs(df1[i-5:i+5,j]))<=50 and f[i,j]<500:
            te_end[0,j]=i
            break
#    else: 
#        te_end[0,j]=287
        
te_beg=zeros((1,p),dtype=int) #te_beg=末端波峰位置，即看作散场开始时间
for j in range(p):
    for i in range(n1):
        if ttmax_all[i,j]!=0:
            if ttmax_all[i,j]+5<te_end[0,j]:
                te_beg[0,j]=ttmax_all[i,j]
       
td=range(0,n)
for j in range(p):
    fig=plt.plot(td,f[:,j],te_beg[0,j],f[te_beg[0,j],j],'x',te_end[0,j],f[te_end[0,j],j],'o')
    plt.title("fig"+np.str(j))
    plt.show()           

j=26
fig=plt.plot(td,f[:,j],te_beg[0,j],f[te_beg[0,j],j],'x',te_end[0,j],f[te_end[0,j],j],'o')
    
rdfmean=zeros((1,p))
rdf=df1/f
s=0
k=0
for j in range(p):            
    if te_beg[0,j]>0 and te_end[0,j]>0:
        rdfmean[0,j]=np.mean(rdf[te_beg[0,j]+1:te_end[0,j],j],axis=0)
    if dfmean[0,j]<0:
       s=s+rdfmean[0,j]
       k=k+1
m=s/k #平均下降率，每20分钟下降m个人
print(m,k)

#[n4,p4]=np.shape(y4)    
#f4=np.zeros((n4,p4), dtype=np.int)
#for i in range(n4-3): 
#    for j in range (0,p4):
#        f4[i+3,j]=y4[i,j]+y4[i+1,j]+y4[i+2,j]+y4[i+3,j]
#
#td2=range(0,n4)
#for j in range(p4):
#    fig=plt.plot(td2,f4[:,j])
#    plt.title("fig"+np.str(j))
#    plt.show()        
#
#t1=np.array([259,259,240])
#t2=np.array([336,330,302])
#t3=np.array([287,287,287])
##tt=np.hstack((t,t))
#
#for j in range(p4):
#    fig=plt.plot(td2,f4[:,j],t1[j],f4[t1[j],j],'x',t2[j],f4[t2[j],j],'o',t3[j],f4[t3[j],j],'v',td2,ff4[:,j])
#    plt.show()
#
#f4mean=zeros((1,p4))
#for j in range(p4):
#    f4mean[0,j]=np.mean(f4[36:97,j])
#
#ff4=f4
#for j in range(p4):
#    for i in range(t1[j],t3[j]):
#        if ff4[i,j]>f4mean[0,j]:
#            ff4[i,j]=ff4[i-1,j]*(1+m)