text:
#!/usr/bin/env python
# -*- coding:utf-8 -*-

import sys, getopt
import dpkt
import socket
import binascii
import os
import os.path
import glob
import time
import datetime
import numpy,math
import re
import platform
from orderFlowFile import *
osname=platform.system()
if osname =="Windows":
	print ("Windows tasks")
	import win_inet_pton
elif osname == "Linux":
	print ("Linux tasks")
else:
	print ("Other System tasks")

#计算Kurtosis、Skewness and Standard Error
#The first standardized moment is zero,because the first moment about the mean is always zero.
#The second standardized moment is one,because the second moment about the mean is equal to the variance σ2.
#The third standardized moment is a measure of skewness. 三阶标准矩用于定义偏度
#The fourth standardized moment refers to the kurtosis. 四阶标准矩用于定义峰度,在实际应用中，通常将峰度值做减3处理，使得正态分布的峰度0(正态分布的峰度为常数3，均匀分布的峰度（系数）为常数1.)
def calc_kurt_skew_se(dataArray):
	n=len(dataArray)
	data_mean=dataArray.mean()
	data_std=dataArray.std()
	t_kurt=0.0
	t_skew=0.0
	t_se=0.0
	if not data_std==0:
		for d in dataArray:
			t_kurt+=(d-data_mean)**4/data_std**4
			t_skew+=(d-data_mean)**3/data_std**3
		kurtosis=t_kurt/(n-1)-3
		#kurtosis=t_kurt/(n-1)
		skewness=t_skew/(n-1)
	else:
		kurtosis=-3
		#kurtosis=0
		skewness=0
	standard_error=data_std/math.sqrt(n)
	return kurtosis,skewness,standard_error

#end def calc_kurt_skew_se

def main(argv):
	try:
		opts, args = getopt.getopt(argv,'hr:i:d:m:a:s:',["rootDir=","deviceId=","dataDate=","min=","addr=","skip_pcap="])
	except getopt.GetoptError:
		print 'python buildStream.py -r <data root directory> -i <deviceId> -d <data date> -m <minnum> -a <tunaddr_prefix> -s <skip_pcap>'
		sys.exit(2)
	input_dir=""
	input_min=2
	input_tunaddr=""
	input_skip_pcap='False'
	input_deviceid=''
	input_dataDate=''
	if len(opts) < 1 :
		usage()
		sys.exit()
	for op, value in opts:
		if op in ("-r", "--rootDir"):
			input_dir = value
			print("input_dir:"+input_dir)
		elif op in("-i","--deviceId"):
			input_deviceid = value
			print("input_deviceid:%s" %(input_deviceid))
		elif op in("-d","--dataDate"):
			input_dataDate = value
			print("input_dataDate:%s" %(input_dataDate))
		elif op in("-m","--min"):
			input_min = int(value)
			print("input_min:%s" %(input_min))
		elif op in("-a","--addr"):
			input_tunaddr = value
			print("tunaddr_prefix:%s" %(input_tunaddr))
		elif op in("-s","--skip_pcap"):
			input_skip_pcap = value
			print("skip_pcap:%s" %(input_skip_pcap))
		elif op == '-h':
			usage()
			sys.exit()

	UserIP_prefix=input_tunaddr
	#最少报文数,即流的报文数小于该数则不处理
	MIN_PKTNUM=input_min
	#数据文件存放目录,包括.pcap文件和.socket文件的目录
	#WORKING_DIR='/home/liuzhen/mergedata/44cd4ccc/20160630'
	#WORKING_DIR='/home/liuzhen/mergedata/94346c37/20160816'
	#WORKING_DIR=input_dir
	#DATA_ROOT_DIR='/home/liuzhen/mergedata'
	DATA_ROOT_DIR=input_dir
	#WORKING_DIR=DATA_ROOT_DIR+'/'+input_deviceid+'/'+input_dataDate
	#WORKING_DIR=DATA_ROOT_DIR+os.path.sep+input_deviceid+os.path.sep+input_dataDate
	#构建WORKING_DIR_LIST

	list_deviceid=input_deviceid.split()
	list_dataDate=input_dataDate.split()
	print("input deviceid:%s" %(list_deviceid))
	print("input dataDate:%s" %(list_dataDate))
	WORKING_DIR_LIST=[]
	allFilesInRoot=os.listdir(DATA_ROOT_DIR)
	allDirInRoot=[]
	#print(allFilesInRoot)
	for file in allFilesInRoot:
		if file=="flows":
			continue
		abs_file=DATA_ROOT_DIR+os.path.sep+file
		if os.path.isdir(abs_file) and (not list_deviceid or file in list_deviceid):
			allDirInRoot.append(file)
	print("all deviceid:%s" %(allDirInRoot))
	for deviceidDir in allDirInRoot:
		allFilesInDeviceidDir=os.listdir(DATA_ROOT_DIR+os.path.sep+deviceidDir)
		for dataDateDir in allFilesInDeviceidDir:
			abs_dataDateDir=DATA_ROOT_DIR+os.path.sep+deviceidDir+os.path.sep+dataDateDir
			if os.path.isdir(abs_dataDateDir) and (not list_dataDate or dataDateDir in list_dataDate):
				WORKING_DIR_LIST.append(abs_dataDateDir)
	print(WORKING_DIR_LIST)
	for WORKING_DIR in WORKING_DIR_LIST:
		#组流的结果存放到数据文件根目录的flows/目录下
		#FLOW_DIR=WORKING_DIR+'/flow'
		#FLOW_DIR=DATA_ROOT_DIR+'/flows/'+input_deviceid+'_'+input_dataDate+'_flow'
		#FLOW_DIR=DATA_ROOT_DIR+os.path.sep+'flows'+os.path.sep+input_deviceid+'_'+input_dataDate+'_flow'
		print("process %s" %(WORKING_DIR))
		tmp1=os.path.split(WORKING_DIR)
		dataDate_field=tmp1[1]
		tmp2=os.path.split(tmp1[0])
		deviceId_field=tmp2[1]
		FLOW_DIR=DATA_ROOT_DIR+os.path.sep+'flows'+os.path.sep+deviceId_field+'_'+dataDate_field+'_flow'
		os.chdir(DATA_ROOT_DIR)
		if not os.path.exists(FLOW_DIR):
			os.makedirs(FLOW_DIR)

		#处理的异常数据存放到flow/目录下的exception.log文件中
		#EXCEPTION_LOG=FLOW_DIR+'/exception.log'
		EXCEPTION_LOG=FLOW_DIR+os.path.sep+'exception.log'
		exceptionLog=open(EXCEPTION_LOG,'w')
		#exceptionLog=open(EXCEPTION_LOG,'a')

		if input_skip_pcap=='False' or input_skip_pcap=='false':
			#
			print '开始处理.pcap文件进行组流'
			#pcapFileName='1_00001_20160630183417.pcap'
			#pcapFileName='1_00001_20160630070554.pcap'
			#pcapFileName='1_00001_20160816190922.pcap'
			#pcapFileNames=os.listdir(WORKING_DIR)
			#使用glob模块支持通配符,且输出的为全路径文件名
			pcapFileNames=glob.glob(WORKING_DIR+os.path.sep+'*.pcap')
			#排序,/home/liuzhen/mergedata/243aa171/20170516/1_00003_20170516003426.pcap,提取出20170516003426出来排序.
			#-------------------------------------------------------^*************^----
			#------------------------------------------------------(-19)---------(-5)--
			pcapFileNames.sort(key=lambda x:int(x[-19:-5]))
			#flowFileNames=[]
			#每个流会创建一个文件，所以注意linux同时打开文件数的限制，使用ulimit -n 10000限制可打开1万个文件
			flowFiles={}
			#需限制打开的文件数目
			MAXOPEN_NUM=300
			openFiles={}#记录文件打开使用记录，filename:count，计数规则:count=count/2+1或0，即每次更新计数时除2，文件被选中操作+1，其它+0
			notFirstOpen=[]#记录曾经打开过的文件名
			for name in pcapFileNames:
				fullName=name
				#前面改使用glob模块,无需再组合全路径已经检查是否为.pcap文件
				#fullName=os.path.join(WORKING_DIR,name)
				#if not (os.path.isfile(fullName) and fullName.endswith('.pcap')):
				#	continue
				print("process %s..." %(name))
				pcapFilePrefix=fullName.split('.')[0]
				pcapF=open(fullName,'rb') #windlows下运行出现dpkt.dpkt.NeedData ERROR,加上'rb'打开模式就没有问题
				pcap=dpkt.pcap.Reader(pcapF)

				#(ts,frame)=pcap.next()
				index=0
				skip=False
				for ts,frame in pcap:
					index=index+1
					#eth=dpkt.ethernet.Ethernet(frame) #//File encapsulation: Ethernet type
					#print `eth`
					#ip=eth.data
					try:
						ip=dpkt.ip.IP(frame) #//File encapsulation: Raw IP
						skip=False
					except Exception as ex:
						skip=True
						exceptionLog.write(pcapFilePrefix+'-'+str(index)+' PKT ERROR.'+str(ex)+'\n')
					if skip:
						continue
					srcAddr=socket.inet_ntop(socket.AF_INET,ip.src)
					dstAddr=socket.inet_ntop(socket.AF_INET,ip.dst)
					#如10.*.*.2 userIP(用户IP运行时由-a <tunaddr> 参数传入)
					#srcAddr is userIP & dstAddr is not userIP: direction=1, 流出流量
					#srcAddr is not userIP & dstAddr is userIP: direction=0, 流入流量
					direction=0
					#if srcAddr.startswith("10."):
					#	if not dstAddr.startswith("10."):
					if srcAddr.startswith(UserIP_prefix):
						if not dstAddr.startswith(UserIP_prefix):
							direction=1
							ipUser=srcAddr
							ipInternet=dstAddr
						else:
							exceptionLog.write(pcapFilePrefix+'-'+str(index)+' '+srcAddr+' '+dstAddr+'\n')
							#print(index,srcAddr,dstAddr)
							continue
					else:
						if dstAddr.startswith(UserIP_prefix):
							direction=0
							ipUser=dstAddr
							ipInternet=srcAddr
						else:
							exceptionLog.write(pcapFilePrefix+'-'+str(index)+' '+srcAddr+' '+dstAddr+'\n')
							#print(srcAddr,dstAddr)
							continue
					#ip.p //17: UDP    6: TCP	1:ICMP
					if ip.p==17:
						protocol='UDP'
					elif ip.p==6:
						protocol='TCP'
					else:
						exceptionLog.write(pcapFilePrefix+'-'+str(index)+' NOT UDP|TCP ip.p='+str(ip.p)+'\n')
						#print('pass ip.p=%s' %(ip.p))
						continue
					try:
						pkt=ip.data
						payload=pkt.data
						sport=pkt.sport
						dport=pkt.dport
						skip=False
					except Exception as ex:
						skip=True
						exceptionLog.write(pcapFilePrefix+'-'+str(index)+' PKT ERROR.'+str(ex)+'\n')
					if skip:
						continue
					if direction==1:
						portUser=pkt.sport
						portInternet=pkt.dport
					else:
						portUser=pkt.dport
						portInternet=pkt.sport
					timeStr=time.strftime("%Y-%m-%d %H:%M:%S",time.localtime(ts))
					tt=`ts`.split('.')
					#如果载荷长度为0
					#if len(payload)==0:
					#	exceptionLog.write('%s %s.%s %s %s %s %s %s %s %s\n' %(index,timeStr,tt[1],ipUser,portUser,ipInternet,portInternet,protocol,direction,binascii.hexlify(str(payload))))
					#else:
						#print('%s.%s' %(timeStr,tt[1]),ipUser,portUser,ipInternet,portInternet,protocol,direction,binascii.hexlify(str(payload)))
					flowKey=ipUser+'-'+str(portUser)+'-'+protocol+'-'+ipInternet+'-'+str(portInternet)
					if flowFiles.has_key(flowKey):
						flowFile=flowFiles.get(flowKey)
					else:
						#检查打开文件是否超过阈值,超过一个阈值则要选择计数最小的打开文件关闭，然后然后打开新的文件
						#openFiles记录了打开了哪些文件，并对每个文件打开访问过程进行计数，每打开一次计数/2再加1，每没打开一次计数/2
						if len(openFiles)>MAXOPEN_NUM:
							mincount=1
							for openFileName,openFileCount in openFiles.items():
								if openFileCount<=mincount:
									mincount=openFileCount
									min_FileName=openFileName
								#endif
							#end for openFileName
							flowFiles[min_FileName].close()		#关闭选中的文件
							del flowFiles[min_FileName]			#删除打开文件描述
							del openFiles[min_FileName]			#删除打开文件记录
						#end if len(openFiles)
						#flowFile=open(FLOW_DIR+'/'+flowKey+'.flow','w');
						if flowKey in notFirstOpen:
							flowFile=open(FLOW_DIR+os.path.sep+flowKey+'.flow','a');
						else:
							notFirstOpen.append(flowKey)
							flowFile=open(FLOW_DIR+os.path.sep+flowKey+'.flow','w');
							flowFile.write('#pkt_position(pcapFile-index)|packet_direction(1:output,0:input)|ip.len|payload.len|timestamp|timestamp-distance(2016-01-01 0:0:0)|payload\n')
						flowFiles[flowKey]=flowFile
						openFiles[flowKey]=0.0			#初始计数值为0.0，注意用小数

					#pktKey: pcapFilePrefix-index  example: 1_00001_20160630062731-1
					#pktKey: pcapFilePrefix-index
					#pktKey direction bytes timestamp payload
					d1 = datetime.datetime.strptime(timeStr, '%Y-%m-%d %H:%M:%S')
					d2 = datetime.datetime.strptime('2016-01-01 0:0:0', '%Y-%m-%d %H:%M:%S')
					delta = d1 - d2
					ts_distance=delta.total_seconds()+float('0.'+tt[1])
					flowFile.write('%s|%s|%s|%s|%s|%s|%s\n' %(pcapFilePrefix+'-'+str(index),direction,ip.len,len(payload),timeStr+'.'+tt[1],`ts_distance`,binascii.hexlify(str(payload))))
					#更新文件计数
					for name,count in openFiles.items():
						if name==flowKey:
							openFiles[name]=count/2.0+1		#打开的文件，计数/2+1
						else:
							openFiles[name]=count/2.0		#打开文件外的其他文件，计数/2
				pcapF.close()
				#end for ts,frame in pcap
			#
			#end all pcap file
			#
			#关闭所有打开的flow文件
			for key,flowFile in flowFiles.items():
				flowFile.close()
			#调用orderFowFile.py里面的函数,按报文时间大小顺序调整文件里面报文顺序
			OVERWRITE="True"
			#FLOW_DIR
			orderFlowFile(FLOW_DIR,OVERWRITE)
			exceptionLog.close()
		else:
			print '跳过.pcap文件组流'
		#
		#开始处理.socket文件,记录每个socket流同应用的对应关系
		#
		#flowAppMapList:{flowkey--->[[name,created,destroyed],[name,created,destroyed]]}
		#appList:flowAppMapList[flowkey][#index]    flowAppMapList[flowkey].append([])     del flowAppMapList[flowkey][#index]
		#name:flowAppMapList[flowkey][#index][0]
		#created:flowAppMapList[flowkey][#index][1]
		#destroyed:flowAppMapList[flowkey][#index][2]
		flowApp={}
		socketFiles=os.listdir(WORKING_DIR)
		for name in socketFiles:
			fullName=os.path.join(WORKING_DIR,name)
			if not (os.path.isfile(fullName) and fullName.endswith('.socket')):
				continue

			soctetfile=open(fullName,'r')
			#soctetfile=open(WORKING_DIR+'/94346c3720160816160830.socket','r')
			#2016-08-18 07:00:34 | 192.168.1.2, 173.252.74.22, 53917, 443, Facebook, TCP, created at 1471503634
			#2016-08-18 07:00:42 | 192.168.1.2, 98.136.223.38, 35251, 443, Yahoo Mail, TCP, destroyed at 1471503642 (created at 1471503610)
			for line in soctetfile:
				l1=line.split('|')
				#if len(l1) != 2 or '0.0.0.0' in line:
				if len(l1) != 2:
					continue
				l2=l1[1].split(',')
				ipUser=l2[0].strip()
				ipInternet=l2[1].strip()
				portUser=l2[2].strip()
				portInternet=l2[3].strip()
				appName=l2[4].strip()
				protocol=l2[5].strip()
				ts=l2[6].strip()
				if ts.startswith('destroyed'):
					tt=ts.split(' ')
					destroyed=tt[2]
					created=tt[5][0:len(tt[5])-1]
				else:
					tt=ts.split(' ')
					destroyed='-1'
					created=tt[2]
				#socket文件与pcap文件的ipUser不一致，socket文件的ipuser是手机获取到的地址，pcap文件的ipuser是tun接口地址10.*.*.2
				#flowkey=ipUser+'-'+portUser+'-'+protocol+'-'+ipInternet+'-'+portInternet
				flowkey='-'+portUser+'-'+protocol+'-'+ipInternet+'-'+portInternet
				if flowApp.has_key(flowkey):
					#oldName = flowApp[flowkey]
					#if oldName != appName:
					#	print('flowkey:%s, old appName:%s-->new appName:%s' %(flowkey,oldName,appName))
					found='false'
					for index in range(len(flowApp[flowkey])):#遍历流对应的已获取的应用列表
						#flowApp[flowkey][index][0]------name
						#flowApp[flowkey][index][1]------created
						#flowApp[flowkey][index][2]------destroyed
						if appName==flowApp[flowkey][index][0]:#存在与本条socket信息相同应用名称的列表
							if destroyed=='-1':#本条socket信息为created类型，一般满足这个就不会满足前面appName==flowApp[flowkey][index][0]这个条件
								if not flowApp[flowkey][index][2]=='-1' and eval(created)>=eval(flowApp[flowkey][index][1]) and eval(created)<=eval(flowApp[flowkey][index][2]):
									found=True#本条socket创建信息已包含在已有的信息中，可丢弃
								if flowApp[flowkey][index][2]=='-1' and eval(created)==eval(flowApp[flowkey][index][1]):
									found=True#本条socket创建信息重复，可丢弃
								#else:
									#需要添加，flowApp[flowkey].append([appName,created,destroyed])
							else:#本条socket信息为destroyed类型
								if flowApp[flowkey][index][2]=='-1':#存在与本条socket信息对应的created类型信息
									if created==flowApp[flowkey][index][1]:#存在与本条socket信息对应的created类型信息，且创建时间相同
										found=True
										flowApp[flowkey][index][2]=destroyed
									else:#存在与本条socket信息对应的created类型信息，但创建时间不相同
										if eval(created)<=eval(flowApp[flowkey][index][1]) and eval(destroyed)>=eval(flowApp[flowkey][index][1]):#已存在的created类型的socket信息的创建时间已包含在本条socket信息中，将存在的socket信息更新
											found=True
											flowApp[flowkey][index][1]=created
											flowApp[flowkey][index][2]=destroyed
										#else:#已存在的created类型的socket信息的创建时间未包含在本条socket信息中
											#需要添加，flowApp[flowkey].append([appName,created,destroyed])
								else:#不存在与本条socket信息对应的created类型信息（存在的是destroyed类型信息），一般应该已存在一条created类型信息
									if eval(created)<=eval(flowApp[flowkey][index][1]) and eval(destroyed)>=eval(flowApp[flowkey][index][2]):#存在与本条socket信息类似的destroyed类型信息，且时间范围包含在本条socket信息内，将存在的socket信息更新
										found=True
										flowApp[flowkey][index][1]=created
										flowApp[flowkey][index][2]=destroyed
									else:
										if (eval(created)<eval(flowApp[flowkey][index][1]) and eval(destroyed)>eval(flowApp[flowkey][index][1]) and eval(destroyed)<eval(flowApp[flowkey][index][2])) or (eval(created)>eval(flowApp[flowkey][index][1]) and eval(created)<eval(flowApp[flowkey][index][2]) and eval(destroyed)>eval(flowApp[flowkey][index][2])):
										#存在与本条socket信息类似的destroyed类型信息，但时间范围存在交叉，socket数据错误！！
											found=True
											print 'ERROR:'+line
										#else:存在与本条socket信息类似的destroyed类型信息，且时间范围没有交集，需要添加
											#flowApp[flowkey].append([appName,created,destroyed])
					#END for index in range(len(flowApp[flowkey]))
					if not found:
						flowApp[flowkey].append([appName,created,destroyed])
				else:
					flowApp[flowkey]=[]
					flowApp[flowkey].append([appName,created,destroyed])

			soctetfile.close()
		#.socket文件处理完毕

		#测试输出flowApp处理结果
		#{flowkey--->[[name,created,destroyed],[name,created,destroyed]]}
		'''
		print "socket文件缓冲结果\nKey\tAppName\tCreated\tDestroyed"
		for flowkey in flowApp:
			print flowkey
			for index in range(len(flowApp[flowkey])):
				print "\t"+flowApp[flowkey][index][0]+"\t"+flowApp[flowkey][index][1]+"\t"+flowApp[flowkey][index][2]
		'''
		#
		#开始处理前面处理组流生成的流文件
		#

		#ssl
		## Server Hello with certificate | Client Hello
		## This allows SSL 3.X, which includes TLS 1.0, known internally as SSL 3.1
		#^(.?.?\x16\x03.*\x16\x03|.?.?\x01\x03\x01?.*\x0b)
		#python
		regex_str='^((..)?(..)?1603(..)*1603|(..)?(..)?0103(01)?(..)*0b)'
		regex_ssl=re.compile(regex_str)

		FF_EXCEPTION_LOG=FLOW_DIR+os.path.sep+'flowfeature_exception.log'
		ff_exceptionLog=open(FF_EXCEPTION_LOG,'w')
		#print('#flowkey,总报文数,字节数,报文大小(最小,最大,平均,标准差),报文到达时间间隔(最小,最大,平均,标准差),流持续时间,appName')
		#flowFeatureFile=open(FLOW_DIR+'/flow.feature','w')
		flowFeatureFile=open(FLOW_DIR+os.path.sep+'flow.feature','w')
		flowFeatureFile.write('#flowkey(设备IP<1>,设备端口<2>,协议<3>,对端IP<4>,对端端口<5>),created<6>,destroyed<7>,appName<8>,encrypted_tag<9>,总报文数<10-1>,字节数<11-2>,报文大小(最小<12-3>,最大<13-4>,平均<14-5>,标准差<15-6>,峰度<16-7>,偏度<17-8>,标准误差<18-9>),报文到达时间间隔(最小<19-10>,最大<20-11>,平均<21-12>,标准差<22-13>),流持续时间<23-14>,IN总报文数<24-1>,IN字节数<25-2>,IN报文大小(最小<26-3>,最大<27-4>,平均<28-5>,标准差<29-6>,峰度<30-7>,偏度<31-8>,标准误差<32-9>),IN报文到达时间间隔(最小<33-10>,最大<34-11>,平均<35-12>,标准差<36-13>),IN流持续时间<37-14>,OUT总报文数<38-1>,OUT字节数<39-2>,OUT报文大小(最小<40-3>,最大<41-4>,平均<42-5>,标准差<43-6>,峰度<44-7>,偏度<45-8>,标准误差<46-9>),OUT报文到达时间间隔(最小<47-10>,最大<48-11>,平均<49-12>,标准差<50-13>),OUT流持续时间<51-14>\n')
		flowFiles=os.listdir(FLOW_DIR)
		for name in flowFiles:
			fullName=os.path.join(FLOW_DIR,name)
			if not (os.path.isfile(fullName) and fullName.endswith('.flow')):
				continue
			#print 'process flowfile '+fullName
			flowKey=name[0:len(name)-5]
			flowFile=open(fullName,'r')
			skip=flowKey.find('-')
			flowkey_skipIpUser=flowKey[skip:]

			payloadBytes_all=[]
			ipBytes_all=[]
			pktTSs_all=[]
			in_payloadBytes_all=[]
			in_ipBytes_all=[]
			in_pktTSs_all=[]
			out_payloadBytes_all=[]
			out_ipBytes_all=[]
			out_pktTSs_all=[]
			encrypted_tag_all=[]#'normal'
			if not flowApp.has_key(flowkey_skipIpUser) and 'UDP' in flowkey_skipIpUser:
				tt=flowkey_skipIpUser.split('-')
				#-35286-UDP-192.168.1.1-1900(10.77.2.2-35286-UDP-192.168.1.1-1900)
				#.socket文件中类似0.0.0.0的UDP记录只按用户端口匹配，如下面就按-35286-UDP-0.0.0.0-0匹配
				#2016-09-25 04:14:00 | 0.0.0.0, 0.0.0.0, 35286, 0, 芒果TV, UDP, created at 1474776840
				specialKey='-'+tt[1]+'-UDP-0.0.0.0-0'
				flowkey_skipIpUser=specialKey
			if flowApp.has_key(flowkey_skipIpUser):
				for index in range(len(flowApp[flowkey_skipIpUser])):
					payloadBytes_all.append([])
					ipBytes_all.append([])
					pktTSs_all.append([])
					in_payloadBytes_all.append([])
					in_ipBytes_all.append([])
					in_pktTSs_all.append([])
					out_payloadBytes_all.append([])
					out_ipBytes_all.append([])
					out_pktTSs_all.append([])
					encrypted_tag_all.append('normal')
			else:
				ff_exceptionLog.write('flow %s( %s) cannot found appName' %(flowKey,flowkey_skipIpUser))
				#print 'flow '+flowKey+'( '+flowkey_skipIpUser+') cannot found appName'
				continue
			#pkt_position(pcapFile-index)|packet_direction(1:output,0:input)|ip.len|payload.len|timestamp|timestamp-distance(2016-01-01 0:0:0)|payload
			#flow文件内容：/home/liuzhen/mergedata/44cd4ccc/20160615/20160615-17625|1|65|13|2016-06-15 23:14:56.842539|14426096.842539|0c000401083d10d45abecfe95a
			for line in flowFile:
				if line.startswith('#'):
					continue
				tt=line.split('|')
				tttime=tt[4].split('.')
				pkttimestamp=time.mktime(time.strptime(tttime[0],"%Y-%m-%d %H:%M:%S"))+float('0.'+tttime[1]) #秒数表示的时间
				#flowApp[flowkey][index][0]------name
				#flowApp[flowkey][index][1]------created
				#flowApp[flowkey][index][2]------destroyed

				selectedIndex=-1
				minDistance=-1
				for index in range(len(flowApp[flowkey_skipIpUser])):
					#比较报文时间同流时间，选择时间离created最近的为最合适的应用名
					created=flowApp[flowkey_skipIpUser][index][1]
					destroyed=flowApp[flowkey_skipIpUser][index][2]
					#print str(pkttimestamp)+' '+flowApp[flowkey_skipIpUser][index][0]+'('+created+'---'+destroyed+')'
					if minDistance<0 or minDistance>abs(pkttimestamp-eval(created)):
						minDistance=abs(pkttimestamp-eval(created))
						selectedIndex=index
						appName=flowApp[flowkey_skipIpUser][index][0]

				if selectedIndex<0:
					continue
				if tt[1]=='1':
					isOutFlow=True
				else:
					isOutFlow=False
				ipBytes_all[selectedIndex].append(int(tt[2]))
				payloadBytes_all[selectedIndex].append(int(tt[3]))
				if isOutFlow:
					out_ipBytes_all[selectedIndex].append(int(tt[2]))
					out_payloadBytes_all[selectedIndex].append(int(tt[3]))
				else:
					in_ipBytes_all[selectedIndex].append(int(tt[2]))
					in_payloadBytes_all[selectedIndex].append(int(tt[3]))

				pktTSs_all[selectedIndex].append(pkttimestamp)
				if isOutFlow:
					out_pktTSs_all[selectedIndex].append(pkttimestamp)
				else:
					in_pktTSs_all[selectedIndex].append(pkttimestamp)
				match=regex_ssl.match(tt[6])
				if match:
					encrypted_tag_all[selectedIndex]='SSL'

			for index in range(len(flowApp[flowkey_skipIpUser])):
				#flowApp[flowkey][index][0]------name
				#flowApp[flowkey][index][1]------created
				#flowApp[flowkey][index][2]------destroyed
				appName=flowApp[flowkey_skipIpUser][index][0]
				created=flowApp[flowkey_skipIpUser][index][1]
				destroyed=flowApp[flowkey_skipIpUser][index][2]
				encrypted_tag=encrypted_tag_all[index]
				#不考虑方向统计
				pktTSArray=numpy.array(pktTSs_all[index])
				pktTSArray.sort()
				pktIntervals=[]
				for i in range(1,len(pktTSArray)):
					pktIntervals.append(pktTSArray[i]-pktTSArray[i-1])
				#流入方向统计t
				in_pktTSArray=numpy.array(in_pktTSs_all[index])
				in_pktTSArray.sort()
				in_pktIntervals=[]
				for i in range(1,len(in_pktTSArray)):
					in_pktIntervals.append(in_pktTSArray[i]-in_pktTSArray[i-1])
				#流出方向统计
				out_pktTSArray=numpy.array(out_pktTSs_all[index])
				out_pktTSArray.sort()
				out_pktIntervals=[]
				for i in range(1,len(out_pktTSArray)):
					out_pktIntervals.append(out_pktTSArray[i]-out_pktTSArray[i-1])

				#只有MIN_PKTNUM报文的流忽略
				if len(pktTSs_all[index])<MIN_PKTNUM or len(in_pktTSs_all[index])<MIN_PKTNUM or len(out_pktTSs_all[index])<MIN_PKTNUM:
					#print 'flow '+flowKey+'( '+flowkey_skipIpUser+') in pkts num:'+str(len(in_pktTSs_all[index]))+'\tout pkts num:'+str(len(out_pktTSs_all[index]))
					ff_exceptionLog.write('flow %s( %s) IN pkts num:%s\tOUT pkts num:%s\n' %(flowKey,flowkey_skipIpUser,len(in_pktTSs_all[index]),len(out_pktTSs_all[index])))
					continue
				#narray=numpy.array(L),narray.min(),narray.max(),narray.mean(),narray.sum(),narray.var(),narray.std()
				pktIntervalArray=numpy.array(pktIntervals)
				ipByteArray=numpy.array(ipBytes_all[index])
				payloadByteArray=numpy.array(payloadBytes_all[index])
				in_pktIntervalArray=numpy.array(in_pktIntervals)
				in_ipByteArray=numpy.array(in_ipBytes_all[index])
				in_payloadByteArray=numpy.array(in_payloadBytes_all[index])
				out_pktIntervalArray=numpy.array(out_pktIntervals)
				out_ipByteArray=numpy.array(out_ipBytes_all[index])
				out_payloadByteArray=numpy.array(out_payloadBytes_all[index])

				#flowkey,总报文数,报文字节数,报文大小(最小,最大,平均,标准差),报文到达时间间隔(最小,最大,平均,标准差),流持续时间
				pktCount=len(payloadByteArray)
				totalPayloadBytes=payloadByteArray.sum()
				minPayload=payloadByteArray.min()
				maxPayload=payloadByteArray.max()
				meanPayload=payloadByteArray.mean()
				stdPayload=payloadByteArray.std()
				totalPktBytes=ipByteArray.sum()
				minPkt=ipByteArray.min()
				maxPkt=ipByteArray.max()
				meanPkt=ipByteArray.mean()
				stdPkt=ipByteArray.std()
				kurtosisPkt,skewnessPkt,standardErrorPkt=calc_kurt_skew_se(ipByteArray)
				minInterval=pktIntervalArray.min()
				maxInterval=pktIntervalArray.max()
				meanInterval=pktIntervalArray.mean()
				stdInterval=pktIntervalArray.std()
				duration=pktIntervalArray.sum()
				#流入方向统计
				in_pktCount=len(in_payloadByteArray)
				in_totalPayloadBytes=in_payloadByteArray.sum()
				in_minPayload=in_payloadByteArray.min()
				in_maxPayload=in_payloadByteArray.max()
				in_meanPayload=in_payloadByteArray.mean()
				in_stdPayload=in_payloadByteArray.std()
				in_totalPktBytes=in_ipByteArray.sum()
				in_minPkt=in_ipByteArray.min()
				in_maxPkt=in_ipByteArray.max()
				in_meanPkt=in_ipByteArray.mean()
				in_stdPkt=in_ipByteArray.std()
				in_kurtosisPkt,in_skewnessPkt,in_standardErrorPkt=calc_kurt_skew_se(in_ipByteArray)
				in_minInterval=in_pktIntervalArray.min()
				in_maxInterval=in_pktIntervalArray.max()
				in_meanInterval=in_pktIntervalArray.mean()
				in_stdInterval=in_pktIntervalArray.std()
				in_duration=in_pktIntervalArray.sum()
				#流出方向统计
				out_pktCount=len(out_payloadByteArray)
				out_totalPayloadBytes=out_payloadByteArray.sum()
				out_minPayload=out_payloadByteArray.min()
				out_maxPayload=out_payloadByteArray.max()
				out_meanPayload=out_payloadByteArray.mean()
				out_stdPayload=out_payloadByteArray.std()
				out_totalPktBytes=out_ipByteArray.sum()
				out_minPkt=out_ipByteArray.min()
				out_maxPkt=out_ipByteArray.max()
				out_meanPkt=out_ipByteArray.mean()
				out_stdPkt=out_ipByteArray.std()
				out_kurtosisPkt,out_skewnessPkt,out_standardErrorPkt=calc_kurt_skew_se(out_ipByteArray)
				out_minInterval=out_pktIntervalArray.min()
				out_maxInterval=out_pktIntervalArray.max()
				out_meanInterval=out_pktIntervalArray.mean()
				out_stdInterval=out_pktIntervalArray.std()
				out_duration=out_pktIntervalArray.sum()

				#print('%s|%s|%s|%s|%s|%s|%s|%s|%s|%s|%s|%s' %(flowKey,pktCount,totalBytes,minPayload,maxPayload,meanPayload,stdPayload,minInterval,maxInterval,meanInterval,stdInterval,duration))
				'''
				skip=flowKey.find('-')
				flowkey_skipIpUser=flowKey[skip:]
				if flowApp.has_key(flowkey_skipIpUser):
					appName=flowApp[flowkey_skipIpUser]
				elif 'UDP' in flowKey:
					tt=flowKey.split('-')
					#-35286-UDP-192.168.1.1-1900(10.77.2.2-35286-UDP-192.168.1.1-1900)
					#.socket文件中类似0.0.0.0的UDP记录只按用户端口匹配，如下面就按-35286-UDP-0.0.0.0-0匹配
					#2016-09-25 04:14:00 | 0.0.0.0, 0.0.0.0, 35286, 0, 芒果TV, UDP, created at 1474776840
					if flowApp.has_key('-'+tt[1]+'-UDP-0.0.0.0-0'):
						appName=flowApp['-'+tt[1]+'-UDP-0.0.0.0-0']
					elif tt[1]=='53' or tt[4]=='53':
						appName='DNS'
					else:
						appName='NOT FOUND'
				else:
					appName='NOT FOUND'
				'''

				flowFeatureFile.write('%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s,%s\n' %(flowKey.replace('-',','),created,destroyed,appName,encrypted_tag,pktCount,totalPktBytes,minPkt,maxPkt,meanPkt,stdPkt,kurtosisPkt,skewnessPkt,standardErrorPkt,minInterval,maxInterval,meanInterval,stdInterval,duration,in_pktCount,in_totalPktBytes,in_minPkt,in_maxPkt,in_meanPkt,in_stdPkt,in_kurtosisPkt,in_skewnessPkt,in_standardErrorPkt,in_minInterval,in_maxInterval,in_meanInterval,in_stdInterval,in_duration,out_pktCount,out_totalPktBytes,out_minPkt,out_maxPkt,out_meanPkt,out_stdPkt,out_kurtosisPkt,out_skewnessPkt,out_standardErrorPkt,out_minInterval,out_maxInterval,out_meanInterval,out_stdInterval,out_duration))
			flowFile.close()
			#if appName == 'NOT FOUND' or appName.startswith('NOTFOUND'):
			#	print('remove file:%s' %(fullName))
			#	os.remove(fullName)

		#print "process completed"
		flowFeatureFile.close()
		ff_exceptionLog.close()

def usage():
	print 'python buildStream.py -r <root data directory> -i <deviceId> -d <dataDate>  -m <minnum> -a <tunaddr_prefix>'
	print 'buildStream.py usage:'
	print '-r, --rootDir: .pcap & .socket file root directory'
	print '-i, --deviceId: device id'
	print '-d, --dataDate: data date'
	print '-m, --min: min pkt num'
	print '-a, --addr: tunaddr_prefix'
	print '-s, --skip_pcap: skip_pcap'
	print 'example:'
	print '\tpython buildStream.py -r /home/liuzhen/mergedata -m 2 -a 10. -s True'
	print '\tpython buildStream.py -r /home/liuzhen/mergedata -i 44cd4ccc -m 2 -a 10. -s True'
	print '\tpython buildStream.py -r /home/liuzhen/mergedata -i "44cd4ccc 66aed7a1" -m 2 -a 10. -s True'
	print '\tpython buildStream.py -r /home/liuzhen/mergedata -d 20160821 -m 2 -a 10. -s True'
	print '\tpython buildStream.py -r /home/liuzhen/mergedata -d "20160821 20160921" -m 2 -a 10. -s True'
	print '\tpython buildStream.py -r /home/liuzhen/mergedata -i 44cd4ccc -d 20160821 -m 2 -a 10. -s True'
	print '\tpython buildStream.py -r /home/liuzhen/mergedata -i "44cd4ccc 66aed7a1" -d "20160821 20160921" -m 2 -a 10. -s True'
	print '\tpython buildStream.py -r /home/liuzhen/mergedata -i 44cd4ccc -d 20160821 -m 2 -a 10. -s False'
	print 'result:\tplease see directory /home/liuzhen/mergedata/flows/44cd4ccc_20160821_flow/'

if __name__=='__main__':
	main(sys.argv[1:])


text:
"""MCMC Test Rig for COVID-19 UK model"""
import argparse
import os
from time import perf_counter

import tqdm
import yaml
import h5py
import numpy as np
import tensorflow as tf
import tensorflow_probability as tfp

from covid import config
from covid.impl.util import compute_state
from covid.impl.mcmc import UncalibratedLogRandomWalk, random_walk_mvnorm_fn
from covid.impl.event_time_mh import UncalibratedEventTimesUpdate
from covid.impl.occult_events_mh import UncalibratedOccultUpdate, TransitionTopology
from covid.impl.gibbs import DeterministicScanKernel, GibbsStep, flatten_results
from covid.impl.multi_scan_kernel import MultiScanKernel

import model_spec

# Uncomment to block GPU use
# os.environ["CUDA_VISIBLE_DEVICES"] = "-1"
if tf.test.gpu_device_name():
    print("Using GPU")
else:
    print("Using CPU")

tfd = tfp.distributions
tfb = tfp.bijectors
DTYPE = config.floatX

if __name__ == "__main__":

    # Read in settings
    parser = argparse.ArgumentParser()
    parser.add_argument(
        "-c", "--config", default="example_config.yaml", help="configuration file",
    )
    args = parser.parse_args()
    print("Loading config file:", args.config)

    with open(args.config, "r") as f:
        config = yaml.load(f, Loader=yaml.FullLoader)

    covar_data = model_spec.read_covariates(config["data"])

    # We load in cases and impute missing infections first, since this sets the
    # time epoch which we are analysing.
    cases = model_spec.read_cases(config["data"]["reported_cases"])

    # Impute censored events, return cases
    events = model_spec.impute_censored_events(cases)

    # Initial conditions are calculated by calculating the state
    # at the beginning of the inference period
    #
    # Imputed censored events that pre-date the first I-R events
    # in the cases dataset are discarded.  They are only used to
    # to set up a sensible initial state.
    state = compute_state(
        initial_state=tf.concat(
            [covar_data["N"], tf.zeros_like(events[:, 0, :])], axis=-1
        ),
        events=events,
        stoichiometry=model_spec.STOICHIOMETRY,
    )
    start_time = state.shape[1] - cases.shape[1]
    initial_state = state[:, start_time, :]
    events = events[:, start_time:, :]
    num_metapop = covar_data["N"].shape[0]

    ########################################################
    # Build the model, and then construct the MCMC kernels #
    ########################################################
    model = model_spec.CovidUK(
        covariates=covar_data,
        initial_state=initial_state,
        initial_step=0,
        num_steps=events.shape[1],
    )

    # Full joint log posterior distribution
    # $\pi(\theta, \xi, y^{se}, y^{ei} | y^{ir})$
    def logp(theta, xi, events):
        return model.log_prob(
            dict(beta1=theta[0], beta2=theta[1], gamma=theta[2], xi=xi, seir=events,)
        )

    # Build Metropolis within Gibbs sampler
    #
    # Kernels are:
    #     Q(\theta, \theta^\prime)
    #     Q(\xi, \xi^\prime)
    #     Q(Z^{se}, Z^{se\prime}) (partially-censored)
    #     Q(Z^{ei}, Z^{ei\prime}) (partially-censored)
    #     Q(Z^{se}, Z^{se\prime}) (occult)
    #     Q(Z^{ei}, Z^{ei\prime}) (occult)
    def make_theta_kernel(scale, bounded_convergence, name):
        return GibbsStep(
            0,
            tfp.mcmc.MetropolisHastings(
                inner_kernel=UncalibratedLogRandomWalk(
                    target_log_prob_fn=logp,
                    new_state_fn=random_walk_mvnorm_fn(scale, p_u=bounded_convergence),
                )
            ),
            name=name,
        )

    def make_xi_kernel(scale, bounded_convergence, name):
        return GibbsStep(
            1,
            tfp.mcmc.RandomWalkMetropolis(
                target_log_prob_fn=logp,
                new_state_fn=random_walk_mvnorm_fn(scale, p_u=bounded_convergence),
            ),
            name=name,
        )

    def make_partially_observed_step(
        target_event_id, prev_event_id=None, next_event_id=None, name=None
    ):
        return GibbsStep(
            2,
            tfp.mcmc.MetropolisHastings(
                inner_kernel=UncalibratedEventTimesUpdate(
                    target_log_prob_fn=logp,
                    target_event_id=target_event_id,
                    prev_event_id=prev_event_id,
                    next_event_id=next_event_id,
                    initial_state=initial_state,
                    dmax=config["mcmc"]["dmax"],
                    mmax=config["mcmc"]["m"],
                    nmax=config["mcmc"]["nmax"],
                )
            ),
            name=name,
        )

    def make_occults_step(prev_event_id, target_event_id, next_event_id, name):
        return GibbsStep(
            2,
            tfp.mcmc.MetropolisHastings(
                inner_kernel=UncalibratedOccultUpdate(
                    target_log_prob_fn=logp,
                    topology=TransitionTopology(
                        prev_event_id, target_event_id, next_event_id
                    ),
                    cumulative_event_offset=initial_state,
                    nmax=config["mcmc"]["occult_nmax"],
                    t_range=(events.shape[1] - 21, events.shape[1]),
                    name=name,
                )
            ),
            name=name,
        )

    # MCMC tracing functions
    def trace_results_fn(_, results):
        """Returns log_prob, accepted, q_ratio"""

        def is_accepted(result):
            if hasattr(result, "is_accepted"):
                return tf.cast(result.is_accepted, DTYPE)
            return is_accepted(result.inner_results)

        def f(result):
            log_prob = result.proposed_results.target_log_prob
            accepted = is_accepted(result)
            q_ratio = result.proposed_results.log_acceptance_correction
            if hasattr(result.proposed_results, "extra"):
                proposed = tf.cast(result.proposed_results.extra, log_prob.dtype)
                return tf.concat([[log_prob], [accepted], [q_ratio], proposed], axis=0)
            return tf.concat([[log_prob], [accepted], [q_ratio]], axis=0)

        def recurse(f, list_or_atom):
            if isinstance(list_or_atom, list):
                return [recurse(f, x) for x in list_or_atom]
            return f(list_or_atom)

        return recurse(f, results)

    # Build MCMC algorithm here.  This will be run in bursts for memory economy
    @tf.function(autograph=False, experimental_compile=True)
    def sample(n_samples, init_state, sigma_theta, sigma_xi):
        with tf.name_scope("main_mcmc_sample_loop"):

            init_state = init_state.copy()

            kernel = DeterministicScanKernel(
                [
                    make_theta_kernel(sigma_theta, 1.0, "theta_kernel"),
                    make_xi_kernel(sigma_xi, 1.0, "xi_kernel"),
                    MultiScanKernel(
                        config["mcmc"]["num_event_time_updates"],
                        DeterministicScanKernel(
                            [
                                make_partially_observed_step(0, None, 1, "se_events"),
                                make_partially_observed_step(1, 0, 2, "ei_events"),
                                make_occults_step(None, 0, 1, "se_occults"),
                                make_occults_step(0, 1, 2, "ei_occults"),
                            ]
                        ),
                    ),
                ],
                name="gibbs_kernel",
            )

            samples, results = tfp.mcmc.sample_chain(
                n_samples, init_state, kernel=kernel, trace_fn=trace_results_fn
            )

            return samples, results

    ####################################
    # Construct bursted MCMC loop here #
    ####################################

    # MCMC Control
    NUM_BURSTS = config["mcmc"]["num_bursts"]
    NUM_BURST_SAMPLES = config["mcmc"]["num_burst_samples"]
    NUM_EVENT_TIME_UPDATES = config["mcmc"]["num_event_time_updates"]
    THIN_BURST_SAMPLES = NUM_BURST_SAMPLES // config["mcmc"]["thin"]
    NUM_SAVED_SAMPLES = THIN_BURST_SAMPLES * NUM_BURSTS

    # RNG stuff
    tf.random.set_seed(2)

    current_state = [
        np.array([0.85, 0.3, 0.25], dtype=DTYPE),
        np.zeros(model.model["xi"]().event_shape[-1], dtype=DTYPE),
        events,
    ]

    # Output Files
    posterior = h5py.File(
        os.path.expandvars(config["output"]["posterior"]),
        "w",
        rdcc_nbytes=1024 ** 2 * 400,
        rdcc_nslots=100000,
        libver="latest",
    )
    event_size = [NUM_SAVED_SAMPLES] + list(current_state[2].shape)

    posterior.create_dataset("initial_state", data=initial_state)

    # Ideally we insert the inference period into the posterior file
    # as this allows us to post-attribute it to the data.  Maybe better
    # to simply save the data into it as well.
    posterior.create_dataset("config", data=yaml.dump(config))
    theta_samples = posterior.create_dataset(
        "samples/theta",
        [NUM_SAVED_SAMPLES, current_state[0].shape[0]],
        dtype=np.float64,
    )
    xi_samples = posterior.create_dataset(
        "samples/xi", [NUM_SAVED_SAMPLES, current_state[1].shape[0]], dtype=np.float64,
    )
    event_samples = posterior.create_dataset(
        "samples/events",
        event_size,
        dtype=DTYPE,
        chunks=(32, 32, 32, 1),
        compression="szip",
        compression_opts=("nn", 16),
    )

    output_results = [
        posterior.create_dataset("results/theta", (NUM_SAVED_SAMPLES, 3), dtype=DTYPE,),
        posterior.create_dataset("results/xi", (NUM_SAVED_SAMPLES, 3), dtype=DTYPE,),
        posterior.create_dataset(
            "results/move/S->E", (NUM_SAVED_SAMPLES, 3 + num_metapop), dtype=DTYPE,
        ),
        posterior.create_dataset(
            "results/move/E->I", (NUM_SAVED_SAMPLES, 3 + num_metapop), dtype=DTYPE,
        ),
        posterior.create_dataset(
            "results/occult/S->E", (NUM_SAVED_SAMPLES, 6), dtype=DTYPE
        ),
        posterior.create_dataset(
            "results/occult/E->I", (NUM_SAVED_SAMPLES, 6), dtype=DTYPE
        ),
    ]
    posterior.swmr_mode = True

    print("Initial logpi:", logp(*current_state))

    theta_scale = tf.constant(
        [
            [1.12e-3, 1.67e-4, 1.61e-4],
            [1.67e-4, 7.41e-4, 4.68e-5],
            [1.61e-4, 4.68e-5, 1.28e-4],
        ],
        dtype=DTYPE,
    )
    theta_scale = theta_scale * 0.2 / theta_scale.shape[0]

    xi_scale = tf.eye(current_state[1].shape[0], dtype=DTYPE)
    xi_scale = xi_scale * 0.0001 / xi_scale.shape[0]

    # We loop over successive calls to sample because we have to dump results
    #   to disc, or else end OOM (even on a 32GB system).
    # with tf.profiler.experimental.Profile("/tmp/tf_logdir"):
    for i in tqdm.tqdm(range(NUM_BURSTS), unit_scale=NUM_BURST_SAMPLES):
        samples, results = sample(
            NUM_BURST_SAMPLES,
            init_state=current_state,
            sigma_theta=theta_scale,
            sigma_xi=xi_scale,
        )
        current_state = [s[-1] for s in samples]
        s = slice(i * THIN_BURST_SAMPLES, i * THIN_BURST_SAMPLES + THIN_BURST_SAMPLES)
        idx = tf.constant(range(0, NUM_BURST_SAMPLES, config["mcmc"]["thin"]))
        theta_samples[s, ...] = tf.gather(samples[0], idx)
        xi_samples[s, ...] = tf.gather(samples[1], idx)
        # cov = np.cov(
        #     np.log(theta_samples[: (i * NUM_BURST_SAMPLES + NUM_BURST_SAMPLES), ...]),
        #     rowvar=False,
        # )
        print(current_state[0].numpy(), flush=True)
        # print(cov, flush=True)
        # if (i * NUM_BURST_SAMPLES) > 1000 and np.all(np.isfinite(cov)):
        #     theta_scale = 2.38 ** 2 * cov / 2.0

        start = perf_counter()
        event_samples[s, ...] = tf.gather(samples[2], idx)
        end = perf_counter()

        flat_results = flatten_results(results)
        for i, ro in enumerate(output_results):
            ro[s, ...] = tf.gather(flat_results[i], idx)

        posterior.flush()
        print("Storage time:", end - start, "seconds")
        print(
            "Acceptance theta:",
            tf.reduce_mean(tf.cast(flat_results[0][:, 1], tf.float32)),
        )
        print(
            "Acceptance xi:", tf.reduce_mean(tf.cast(flat_results[1][:, 1], tf.float32))
        )
        print(
            "Acceptance move S->E:",
            tf.reduce_mean(tf.cast(flat_results[2][:, 1], tf.float32)),
        )
        print(
            "Acceptance move E->I:",
            tf.reduce_mean(tf.cast(flat_results[3][:, 1], tf.float32)),
        )
        print(
            "Acceptance occult S->E:",
            tf.reduce_mean(tf.cast(flat_results[4][:, 1], tf.float32)),
        )
        print(
            "Acceptance occult E->I:",
            tf.reduce_mean(tf.cast(flat_results[5][:, 1], tf.float32)),
        )

    print(f"Acceptance theta: {output_results[0][:, 1].mean()}")
    print(f"Acceptance xi: {output_results[1][:, 1].mean()}")
    print(f"Acceptance move S->E: {output_results[2][:, 1].mean()}")
    print(f"Acceptance move E->I: {output_results[3][:, 1].mean()}")
    print(f"Acceptance occult S->E: {output_results[4][:, 1].mean()}")
    print(f"Acceptance occult E->I: {output_results[5][:, 1].mean()}")

    posterior.close()

ext: py

reason: import numpy


text:
# -*- coding: utf-8 -*-
# Copyright 2021, SERTIT-ICube - France, https://sertit.unistra.fr/
# This file is part of sertit-utils project
#     https://github.com/sertit/sertit-utils
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
""" Script testing raster function (with rasterio) """
import os
import shutil
import tempfile

import numpy as np
import pytest
import rasterio
import shapely

from CI.SCRIPTS.script_utils import rasters_path, s3_env
from sertit import ci, rasters_rio, vectors

@s3_env
def test_rasters_rio():
    """Test raster functions"""
    # Rasters
    raster_path = rasters_path().joinpath("raster.tif")
    raster_masked_path = rasters_path().joinpath("raster_masked.tif")
    raster_cropped_path = rasters_path().joinpath("raster_cropped.tif")
    raster_sieved_path = rasters_path().joinpath("raster_sieved.tif")
    raster_to_merge_path = rasters_path().joinpath("raster_to_merge.tif")
    raster_merged_gtiff_path = rasters_path().joinpath("raster_merged.tif")

    # Vectors
    mask_path = rasters_path().joinpath("raster_mask.geojson")
    extent_path = rasters_path().joinpath("extent.geojson")
    footprint_path = rasters_path().joinpath("footprint.geojson")
    if shapely.__version__ >= "1.8a1":
        vect_truth_path = rasters_path().joinpath("vector.geojson")
        diss_truth_path = rasters_path().joinpath("dissolved.geojson")
    else:
        print("USING OLD VECTORS")
        vect_truth_path = rasters_path().joinpath("vector_old.geojson")
        diss_truth_path = rasters_path().joinpath("dissolved_old.geojson")
    nodata_truth_path = rasters_path().joinpath("nodata.geojson")
    valid_truth_path = rasters_path().joinpath("valid.geojson")

    # Create tmp file
    # VRT needs to be build on te same disk
    with tempfile.TemporaryDirectory() as tmp_dir:
        # tmp_dir = rasters_path().joinpath("OUTPUT_RASTERIO")
        # os.makedirs(tmp_dir, exist_ok=True)

        # Get Extent
        extent = rasters_rio.get_extent(raster_path)
        truth_extent = vectors.read(extent_path)
        ci.assert_geom_equal(extent, truth_extent)

        # Get Footprint
        footprint = rasters_rio.get_footprint(raster_path)
        truth_footprint = vectors.read(footprint_path)
        ci.assert_geom_equal(footprint, truth_footprint)

        # Read
        raster, meta = rasters_rio.read(raster_path)
        with rasterio.open(str(raster_path)) as dst:
            raster_1, meta1 = rasters_rio.read(dst, resolution=20)
            raster_2, _ = rasters_rio.read(dst, resolution=[20, 20])
            raster_3, _ = rasters_rio.read(dst, size=(meta1["width"], meta1["height"]))
            raster_4, _ = rasters_rio.read(raster_path, indexes=[1])
            with pytest.raises(ValueError):
                rasters_rio.read(dst, resolution=[20, 20, 20])

            assert raster.shape == (dst.count, dst.height, dst.width)
            assert meta["crs"] == dst.crs
            assert meta["transform"] == dst.transform
            np.testing.assert_array_equal(raster_1, raster_2)
            np.testing.assert_array_equal(raster_1, raster_3)
            np.testing.assert_array_equal(raster, raster_4)  # 2D array

            # Write
            raster_out = os.path.join(tmp_dir, "test.tif")
            rasters_rio.write(raster, meta, raster_out)
            assert os.path.isfile(raster_out)

            # Mask
            raster_masked_out = os.path.join(tmp_dir, "test_mask.tif")
            mask = vectors.read(mask_path)
            mask_arr, mask_meta = rasters_rio.mask(dst, mask.geometry)
            rasters_rio.write(mask_arr, mask_meta, raster_masked_out)

            # Crop
            raster_cropped_out = os.path.join(tmp_dir, "test_crop.tif")
            crop_arr, crop_meta = rasters_rio.crop(dst, mask)
            rasters_rio.write(crop_arr, crop_meta, raster_cropped_out)

            # Sieve
            sieve_out = os.path.join(tmp_dir, "test_sieved.tif")
            sieve_arr, sieve_meta = rasters_rio.sieve(
                raster, meta, sieve_thresh=20, connectivity=4
            )
            rasters_rio.write(sieve_arr, sieve_meta, sieve_out, nodata=255)

            # Collocate
            coll_arr, coll_meta = rasters_rio.collocate(
                meta, raster, meta
            )  # Just hope that it doesnt crash
            assert coll_meta == meta

            # Merge GTiff
            raster_merged_gtiff_out = os.path.join(tmp_dir, "test_merged.tif")
            rasters_rio.merge_gtiff(
                [raster_path, raster_to_merge_path],
                raster_merged_gtiff_out,
                method="max",
            )

            # Vectorize
            val = 2
            vect = rasters_rio.vectorize(raster_path)
            vect_val = rasters_rio.vectorize(raster_path, values=val)
            vect_val_diss = rasters_rio.vectorize(
                raster_path, values=val, dissolve=True
            )
            vect_val_disc = rasters_rio.vectorize(
                raster_path, values=[1, 255], keep_values=False
            )
            vect_truth = vectors.read(vect_truth_path)
            diss_truth = vectors.read(diss_truth_path)
            ci.assert_geom_equal(vect, vect_truth)
            ci.assert_geom_equal(vect_val, vect_truth.loc[vect_truth.raster_val == val])
            ci.assert_geom_equal(vect_val_diss, diss_truth)
            ci.assert_geom_equal(
                vect_val_disc, vect_truth.loc[vect_truth.raster_val == val]
            )

            # Get valid vec
            valid_vec = rasters_rio.get_valid_vector(raster_path)
            valid_truth = vectors.read(valid_truth_path)
            ci.assert_geom_equal(valid_vec, valid_truth)

            # Get nodata vec
            nodata_vec = rasters_rio.get_nodata_vector(raster_path)
            nodata_truth = vectors.read(nodata_truth_path)
            ci.assert_geom_equal(nodata_vec, nodata_truth)

        # Tests
        ci.assert_raster_equal(raster_path, raster_out)
        ci.assert_raster_equal(raster_masked_out, raster_masked_path)
        ci.assert_raster_equal(raster_cropped_out, raster_cropped_path)
        ci.assert_raster_equal(sieve_out, raster_sieved_path)
        ci.assert_raster_equal(raster_merged_gtiff_out, raster_merged_gtiff_path)

@s3_env
@pytest.mark.skipif(
    shutil.which("gdalbuildvrt") is None,
    reason="Only works if gdalbuildvrt can be found.",
)
def test_vrt():
    raster_merged_vrt_path = rasters_path().joinpath("raster_merged.vrt")
    raster_to_merge_path = rasters_path().joinpath("raster_to_merge.tif")
    raster_path = rasters_path().joinpath("raster.tif")

    with tempfile.TemporaryDirectory() as tmp_dir:
        # Merge VRT
        raster_merged_vrt_out = os.path.join(tmp_dir, "test_merged.vrt")
        rasters_rio.merge_vrt(
            [raster_path, raster_to_merge_path], raster_merged_vrt_out
        )
        ci.assert_raster_equal(raster_merged_vrt_out, raster_merged_vrt_path)

def test_dim():
    """Test on BEAM-DIMAP function"""
    dim_path = rasters_path().joinpath("DIM.dim")
    assert rasters_rio.get_dim_img_path(dim_path) == rasters_path().joinpath(
        "DIM.data", "dim.img"
    )

def test_bit():
    """Test bit arrays"""
    np_ones = np.ones((1, 2, 2), dtype=np.uint16)
    ones = rasters_rio.read_bit_array(np_ones, bit_id=0)
    zeros = rasters_rio.read_bit_array(np_ones, bit_id=list(np.arange(1, 15)))
    assert (np_ones == ones).all()
    for arr in zeros:
        assert (np_ones == 1 + arr).all()

@s3_env
def test_dem_fct():
    """ Test DEM fct, ie. slope and hillshade"""
    # Paths IN
    dem_path = rasters_path().joinpath("dem.tif")
    hlsd_path = rasters_path().joinpath("hillshade.tif")
    slope_path = rasters_path().joinpath("slope.tif")
    slope_r_path = rasters_path().joinpath("slope_r.tif")
    slope_p_path = rasters_path().joinpath("slope_p.tif")

    # Create tmp file
    with tempfile.TemporaryDirectory() as tmp_dir:
        # Path OUT
        hlsd_path_out = os.path.join(tmp_dir, "hillshade_out.tif")
        slope_path_out = os.path.join(tmp_dir, "slope.tif")
        slope_r_path_out = os.path.join(tmp_dir, "slope_r.tif")
        slope_p_path_out = os.path.join(tmp_dir, "slope_p.tif")

        # Compute
        hlsd, meta = rasters_rio.hillshade(dem_path, 34.0, 45.2)
        rasters_rio.write(hlsd, meta, hlsd_path_out)

        slp, meta = rasters_rio.slope(dem_path)
        rasters_rio.write(slp, meta, slope_path_out)

        slp_r, meta = rasters_rio.slope(dem_path, in_pct=False, in_rad=True)
        rasters_rio.write(slp_r, meta, slope_r_path_out)

        slp_p, meta = rasters_rio.slope(dem_path, in_pct=True)
        rasters_rio.write(slp_p, meta, slope_p_path_out)

        # Test
        ci.assert_raster_almost_equal(hlsd_path, hlsd_path_out, decimal=4)
        ci.assert_raster_almost_equal(slope_path, slope_path_out, decimal=4)
        ci.assert_raster_almost_equal(slope_r_path, slope_r_path_out, decimal=4)
        ci.assert_raster_almost_equal(slope_p_path, slope_p_path_out, decimal=4)

max_stars_repo_name: sertit/sertit-utils

max_issues_repo_name: sertit/sertit-utils

max_forks_repo_name: sertit/sertit-utils


text:
import os
import utils
import config
import traceback
import logging.config
from luna import LunaExcepion

logging.config.fileConfig("logging.conf")
logger = logging.getLogger()

if __name__ == '__main__':
    try:
        logger.info("------ start ------")
        utils.lock()

        from keras.models import Sequential
        from keras.layers import Conv2D, MaxPooling2D
        from keras.layers import Activation, Dropout, Flatten, Dense
        from keras.preprocessing.image import ImageDataGenerator
        from keras.utils.vis_utils import plot_model
        import numpy as np

        if not os.path.exists(config.result_dir):
            os.mkdir(config.result_dir)

        # モデルを構築
        model = Sequential()
        model.add(Conv2D(32, (3, 3), input_shape=(150, 150, 3)))
        model.add(Activation('relu'))
        model.add(MaxPooling2D(pool_size=(2, 2)))

        model.add(Conv2D(32, (3, 3)))
        model.add(Activation('relu'))
        model.add(MaxPooling2D(pool_size=(2, 2)))

        model.add(Conv2D(64, (3, 3)))
        model.add(Activation('relu'))
        model.add(MaxPooling2D(pool_size=(2, 2)))

        model.add(Flatten())
        model.add(Dense(64))
        model.add(Activation('relu'))
        model.add(Dropout(0.5))
        model.add(Dense(1))
        model.add(Activation('sigmoid'))

        model.compile(
            loss='binary_crossentropy',
            optimizer='adam',
            metrics=['accuracy']
        )

        #model.summary()
        plot_model(model, to_file='model.png')

        # 訓練データとバリデーションデータを生成するジェネレータを作成
        train_datagen = ImageDataGenerator(
            rescale=1.0 / 255,
            shear_range=0.2,
            zoom_range=0.2,
            horizontal_flip=True
        )
        test_datagen = ImageDataGenerator(rescale=1.0 / 255)
        train_generator = train_datagen.flow_from_directory(
            config.train_dir,
            target_size=(150, 150),
            batch_size=32,
            class_mode='binary'
        )
        validation_generator = test_datagen.flow_from_directory(
            config.validation_dir,
            target_size=(150, 150),
            batch_size=32,
            class_mode='binary'
        )
        print(train_generator.class_indices)

        # 訓練
        history = model.fit_generator(
            generator=train_generator,
            steps_per_epoch=int(np.floor(2000/32)),
            epochs=50,
            validation_data=validation_generator,
            validation_steps=int(np.floor(800/32))
        )
        utils.plot_history(history)

        # 結果を保存
        model.save(os.path.join(config.result_dir, 'scratch_model.h5'))
        model.save_weights(os.path.join(config.result_dir, 'scratch_weights.h5'))
        utils.save_history(history, os.path.join(config.result_dir, 'scratch_history.txt'))

    except (KeyboardInterrupt, SystemExit):
        utils.unlock()
        utils.error(config.syserr)
    except LunaExcepion as e:
        utils.error(e.value)
        if (e.value == config.locked):
            exit()
            logger.info("------ end ------")
    except Exception as e:
        logger.error(e)
        logger.error(traceback.format_exc())
        utils.error(config.syserr)
    utils.unlock()
    logger.info("------ end ------")

max_stars_repo_name: mrm-xiefan/lunania-ai

max_issues_repo_name: kevinelly/lunania-ai

max_forks_repo_name: kevinelly/lunania-ai


text:
import numpy as np
import torch

# 相交矩形的面积
def intersect(box_a, box_b):
    """计算两组矩形两两之间相交区域的面积
    Args:
        box_a: (tensor) bounding boxes, Shape: [A, 4].
        box_b: (tensor) bounding boxes, Shape: [B, 4].
    Return:
      (tensor) intersection area, Shape: [A, B].
    """
    A = box_a.size(0)
    B = box_b.size(0)
    max_xy = torch.min(box_a[:, 2:].unsqueeze(1).expand(A, B, 2),
                       box_b[:, 2:].unsqueeze(0).expand(A, B, 2))
    min_xy = torch.max(box_a[:, :2].unsqueeze(1).expand(A, B, 2),
                       box_b[:, :2].unsqueeze(0).expand(A, B, 2))
    inter = torch.clamp((max_xy - min_xy), min=0)
    return inter[:, :, 0] * inter[:, :, 1]

def jaccard(box_a, box_b):
    """计算两组矩形两两之间的iou
    Args:
        box_a: (tensor) bounding boxes, Shape: [A, 4].
        box_b: (tensor) bounding boxes, Shape: [B, 4].
    Return:
        ious: (tensor) Shape: [A, B]
    """
    inter = intersect(box_a, box_b)
    area_a = ((box_a[:, 2]-box_a[:, 0]) *
              (box_a[:, 3]-box_a[:, 1])).unsqueeze(1).expand_as(inter)  # [A, B]
    area_b = ((box_b[:, 2]-box_b[:, 0]) *
              (box_b[:, 3]-box_b[:, 1])).unsqueeze(0).expand_as(inter)  # [A, B]
    union = area_a + area_b - inter
    return inter / union  # [A, B]

def _matrix_nms(bboxes, cate_labels, cate_scores, kernel='gaussian', sigma=2.0):
    """Matrix NMS for multi-class bboxes.
    Args:
        bboxes (Tensor): shape (n, 4)
        cate_labels (Tensor): shape (n), mask labels in descending order
        cate_scores (Tensor): shape (n), mask scores in descending order
        kernel (str):  'linear' or 'gaussian'
        sigma (float): std in gaussian method
    Returns:
        Tensor: cate_scores_update, tensors of shape (n)
    """
    n_samples = len(cate_labels)
    if n_samples == 0:
        return []

    # 计算一个n×n的IOU矩阵，两组矩形两两之间的IOU
    iou_matrix = jaccard(bboxes, bboxes)   # shape: [n_samples, n_samples]
    iou_matrix = iou_matrix.triu(diagonal=1)   # 只取上三角部分

    # label_specific matrix.
    cate_labels_x = cate_labels.expand(n_samples, n_samples)   # shape: [n_samples, n_samples]
    # 第i行第j列表示的是第i个预测框和第j个预测框的类别id是否相同。我们抑制的是同类的预测框。
    label_matrix = (cate_labels_x == cate_labels_x.transpose(1, 0)).float().triu(diagonal=1)   # shape: [n_samples, n_samples]

    # IoU compensation
    # 非同类的iou置为0，同类的iou保留。逐列取最大iou
    compensate_iou, _ = (iou_matrix * label_matrix).max(0)   # shape: [n_samples, ]
    compensate_iou = compensate_iou.expand(n_samples, n_samples).transpose(1, 0)   # shape: [n_samples, n_samples]

    # IoU decay
    # 非同类的iou置为0，同类的iou保留。
    decay_iou = iou_matrix * label_matrix   # shape: [n_samples, n_samples]

    # matrix nms
    if kernel == 'gaussian':
        decay_matrix = torch.exp(-1 * sigma * (decay_iou ** 2))
        compensate_matrix = torch.exp(-1 * sigma * (compensate_iou ** 2))
        decay_coefficient, _ = (decay_matrix / compensate_matrix).min(0)
    elif kernel == 'linear':
        decay_matrix = (1-decay_iou)/(1-compensate_iou)
        decay_coefficient, _ = decay_matrix.min(0)
    else:
        raise NotImplementedError

    # 更新分数
    cate_scores_update = cate_scores * decay_coefficient
    return cate_scores_update

def matrix_nms(bboxes,
               scores,
               score_threshold,
               post_threshold,
               nms_top_k,
               keep_top_k,
               use_gaussian=False,
               gaussian_sigma=2.):
    inds = (scores > score_threshold)
    cate_scores = scores[inds]
    if len(cate_scores) == 0:
        return torch.zeros((1, 6), device=bboxes.device) - 1.0

    inds = inds.nonzero()
    cate_labels = inds[:, 1]
    bboxes = bboxes[inds[:, 0]]

    # sort and keep top nms_top_k
    sort_inds = torch.argsort(cate_scores, descending=True)
    if nms_top_k > 0 and len(sort_inds) > nms_top_k:
        sort_inds = sort_inds[:nms_top_k]
    bboxes = bboxes[sort_inds, :]
    cate_scores = cate_scores[sort_inds]
    cate_labels = cate_labels[sort_inds]

    # Matrix NMS
    kernel = 'gaussian' if use_gaussian else 'linear'
    cate_scores = _matrix_nms(bboxes, cate_labels, cate_scores, kernel=kernel, sigma=gaussian_sigma)

    # filter.
    keep = cate_scores >= post_threshold
    if keep.sum() == 0:
        return torch.zeros((1, 6), device=bboxes.device) - 1.0
    bboxes = bboxes[keep, :]
    cate_scores = cate_scores[keep]
    cate_labels = cate_labels[keep]

    # sort and keep keep_top_k
    sort_inds = torch.argsort(cate_scores, descending=True)
    if len(sort_inds) > keep_top_k:
        sort_inds = sort_inds[:keep_top_k]
    bboxes = bboxes[sort_inds, :]
    cate_scores = cate_scores[sort_inds]
    cate_labels = cate_labels[sort_inds]

    cate_scores = cate_scores.unsqueeze(1)
    cate_labels = cate_labels.unsqueeze(1).float()

    # pred = torch.cat([cate_labels, cate_scores, bboxes], 1)
    # return pred
    return bboxes, cate_scores, cate_labels.type(torch.int32)

# yolo_box
# https://github.com/miemie2013/Pytorch-PPYOLO/blob/master/model/head.py
def yolo_box(conv_output, anchors, stride, num_classes, scale_x_y, im_size, clip_bbox, conf_thresh, use_gpu=False):
    conv_output = conv_output.permute(0, 2, 3, 1)
    conv_shape       = conv_output.shape
    batch_size       = conv_shape[0]
    output_size      = conv_shape[1]
    anchor_per_scale = len(anchors)
    conv_output = conv_output.reshape((batch_size, output_size, output_size, anchor_per_scale, 5 + num_classes))
    conv_raw_dxdy = conv_output[:, :, :, :, 0:2]
    conv_raw_dwdh = conv_output[:, :, :, :, 2:4]
    conv_raw_conf = conv_output[:, :, :, :, 4:5]
    conv_raw_prob = conv_output[:, :, :, :, 5: ]

    rows = torch.arange(0, output_size, dtype=torch.float32, device=conv_raw_dxdy.device)
    cols = torch.arange(0, output_size, dtype=torch.float32, device=conv_raw_dxdy.device)
    rows = rows[np.newaxis, np.newaxis, :, np.newaxis, np.newaxis].repeat((1, output_size, 1, 1, 1))
    cols = cols[np.newaxis, :, np.newaxis, np.newaxis, np.newaxis].repeat((1, 1, output_size, 1, 1))
    offset = torch.cat([rows, cols], dim=-1)
    offset = offset.repeat((batch_size, 1, 1, anchor_per_scale, 1))
    # Grid Sensitive
    pred_xy = (scale_x_y * torch.sigmoid(conv_raw_dxdy) + offset - (scale_x_y - 1.0) * 0.5 ) * stride

    # _anchors = torch.Tensor(anchors, device=conv_raw_dxdy.device)   # RuntimeError: legacy constructor for device type: cpu was passed device type: cuda, but device type must be: cpu
    _anchors = torch.Tensor(anchors)
    if use_gpu and torch.cuda.is_available():
        _anchors = _anchors.cuda()
    pred_wh = (torch.exp(conv_raw_dwdh) * _anchors)

    pred_xyxy = torch.cat([pred_xy - pred_wh / 2, pred_xy + pred_wh / 2], dim=-1)   # 左上角xy + 右下角xy
    pred_conf = torch.sigmoid(conv_raw_conf)
    # mask = (pred_conf > conf_thresh).float()
    pred_prob = torch.sigmoid(conv_raw_prob)
    pred_scores = pred_conf * pred_prob
    # pred_scores = pred_scores * mask
    # pred_xyxy = pred_xyxy * mask

    # paddle中实际的顺序
    # pred_xyxy = pred_xyxy.permute(0, 3, 1, 2, 4)
    # pred_scores = pred_scores.permute(0, 3, 1, 2, 4)

    pred_xyxy = pred_xyxy.reshape((batch_size, output_size*output_size*anchor_per_scale, 4))
    pred_scores = pred_scores.reshape((batch_size, pred_xyxy.shape[1], num_classes))
    im_size = torch.from_numpy(im_size)
    _im_size_h = im_size[:, 0:1]
    _im_size_w = im_size[:, 1:2]
    _im_size = torch.cat([_im_size_w, _im_size_h], 1)
    _im_size = _im_size.unsqueeze(1)
    _im_size = _im_size.repeat((1, pred_xyxy.shape[1], 1))

    pred_x0y0 = pred_xyxy[:, :, 0:2] / output_size / stride * _im_size
    pred_x1y1 = pred_xyxy[:, :, 2:4] / output_size / stride * _im_size
    if clip_bbox:
        x0 = pred_x0y0[:, :, 0:1]
        y0 = pred_x0y0[:, :, 1:2]
        x1 = pred_x1y1[:, :, 0:1]
        y1 = pred_x1y1[:, :, 1:2]
        x0 = torch.where(x0 < 0, x0 * 0, x0)
        y0 = torch.where(y0 < 0, y0 * 0, y0)
        x1 = torch.where(x1 > _im_size[:, :, 0:1], _im_size[:, :, 0:1], x1)
        y1 = torch.where(y1 > _im_size[:, :, 1:2], _im_size[:, :, 1:2], y1)
        pred_xyxy = torch.cat([x0, y0, x1, y1], -1)
    else:
        pred_xyxy = torch.cat([pred_x0y0, pred_x1y1], -1)
    return pred_xyxy, pred_scores


text:
from capacityreq import cap
import numpy as np
train_data=np.load('embedded_train_data.npy')
train_labels=np.load('train_labels.npy')
batch_size=int(10000)
indy=np.random.randint(0, train_labels.shape[0], size=batch_size)
train_labels=train_labels[indy]
train_data=train_data[indy]
bestmemeqcap=9999999999
best_feat_nums=None
num_feats=10
for i in range(1000):

	feat_nums=np.random.randint(0, train_data.shape[1], size=num_feats)
	featurized=train_data[:,feat_nums]
	print(featurized.shape, train_labels.shape)
	train_csv=np.concatenate((featurized,train_labels.reshape(train_labels.shape[0],1)),axis=1)
	print(feat_nums)
	memeqcap=cap(train_csv)
	if memeqcap<bestmemeqcap:
		bestmemeqcap=memeqcap
		best_feat_nums=feat_nums
print(best_feat_nums)
print(bestmemeqcap)


text:
import os
import cv2
import dlib
import numpy as np
import argparse
from wide_resnet import WideResNet

# Get Pose Predictor from dlib
# This allows us to detect landmark points in faces and understand the pose/angle of the face
shape_predictor = dlib.shape_predictor('face_recognition/shape_predictor_68_face_landmarks.dat')
# Get the face recognition model
# This is what gives us the face encodings (numbers that identify the face of a particular person)
face_recognition_model = dlib.face_recognition_model_v1('face_recognition/dlib_face_recognition_resnet_model_v1.dat')
# This is the tolerance for face comparisons
# The lower the number - the stricter the comparison
# To avoid false matches, use lower value
# To avoid false negatives (i.e. faces of the same person doesn't match), use higher value
# 0.5-0.6 works well
TOLERANCE = 0.6

known_faces = []

# This function takes a list of known faces
def compare_face_encodings(faces):
    # Finds the difference between each known face and the given face (that we are comparing)
    # Calculate norm for the differences with each known face
    # Return an array with True/Face values based on whether or not a known face matched with the given face
    # A match occurs when the (norm) difference between a known face and the given face is less than or equal to the TOLERANCE value
    if len(faces) > 0:
        face = faces[0]
        for known_face in known_faces:
            val = (np.linalg.norm(known_face - face, axis=1) <= TOLERANCE)
            if val.any():
                return True

    return False

def get_args():
    parser = argparse.ArgumentParser(description="This script detects faces from web cam input, "
                                                 "and estimates age and gender for the detected faces.",
                                     formatter_class=argparse.ArgumentDefaultsHelpFormatter)
    parser.add_argument("--weight_file", type=str, default=None,
                        help="path to weight file (e.g. weights.18-4.06.hdf5)")
    parser.add_argument("--depth", type=int, default=16,
                        help="depth of network")
    parser.add_argument("--width", type=int, default=8,
                        help="width of network")
    args = parser.parse_args()
    return args

def draw_label(image, point, label, font=cv2.FONT_HERSHEY_SIMPLEX,
               font_scale=1, thickness=2):
    size = cv2.getTextSize(label, font, font_scale, thickness)[0]
    x, y = point
    cv2.rectangle(image, (x, y - size[1]), (x + size[0], y), (255, 0, 0), cv2.FILLED)
    cv2.putText(image, label, point, font, font_scale, (255, 255, 255), thickness)

# This function will take an image and return its face encodings using the neural network
def get_face_encodings(image, detected_faces):
    # Load image using scipy
    #image = scipy.misc.imread(path_to_image)
    # Detect faces using the face detector
    #detected_faces = face_detector(image, 1)
    # Get pose/landmarks of those faces
    # Will be used as an input to the function that computes face encodings
    # This allows the neural network to be able to produce similar numbers for faces of the same people, regardless of camera angle and/or face positioning in the image
    shapes_faces = [shape_predictor(image, face) for face in detected_faces]
    # For every face detected, compute the face encodings
    return [np.array(face_recognition_model.compute_face_descriptor(image, face_pose, 1)) for face_pose in shapes_faces]

def main():
    args = get_args()
    depth = args.depth
    k = args.width
    weight_file = args.weight_file

    if not weight_file:
        weight_file = os.path.join("pretrained_models", "weights.18-4.06.hdf5")

    # for face detection
    detector = dlib.get_frontal_face_detector()

    # load model and weights
    img_size = 64
    model = WideResNet(img_size, depth=depth, k=k)()
    model.load_weights(weight_file)

    # capture video
    cap = cv2.VideoCapture(0)
    cap.set(cv2.CAP_PROP_FRAME_WIDTH, 640)
    cap.set(cv2.CAP_PROP_FRAME_HEIGHT, 480)

    while True:
        # get video frame
        ret, img = cap.read()

        if not ret:
            print("error: failed to capture image")
            return -1

        input_img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
        img_h, img_w, _ = np.shape(input_img)

        # detect faces using dlib detector
        detected = detector(input_img, 1)
        faces = np.empty((len(detected), img_size, img_size, 3))

        for i, d in enumerate(detected):
            x1, y1, x2, y2, w, h = d.left(), d.top(), d.right() + 1, d.bottom() + 1, d.width(), d.height()
            xw1 = max(int(x1 - 0.4 * w), 0)
            yw1 = max(int(y1 - 0.4 * h), 0)
            xw2 = min(int(x2 + 0.4 * w), img_w - 1)
            yw2 = min(int(y2 + 0.4 * h), img_h - 1)
            cv2.rectangle(img, (x1, y1), (x2, y2), (255, 0, 0), 2)
            # cv2.rectangle(img, (xw1, yw1), (xw2, yw2), (255, 0, 0), 2)
            faces[i,:,:,:] = cv2.resize(img[yw1:yw2 + 1, xw1:xw2 + 1, :], (img_size, img_size))

        if len(detected) > 0:
            # predict ages and genders of the detected faces
            results = model.predict(faces)
            predicted_genders = results[0]
            ages = np.arange(0, 101).reshape(101, 1)
            predicted_ages = results[1].dot(ages).flatten()

        # draw results
        for i, d in enumerate(detected):
            label = "{}, {}".format(int(predicted_ages[i]),
                                    "F" if predicted_genders[i][0] > 0.5 else "M")
            draw_label(img, (d.left(), d.top()), label)
            #print(label)
            #print(get_face_encodings(input_img, detected))
            #faces = get_face_encodings(input_img, detected)
            #if not compare_face_encodings(faces):
            #    print("new face")
            #    known_faces.append(faces)

        cv2.imshow("result", img)
        key = cv2.waitKey(30)

        if key == 27:
            break

if __name__ == '__main__':
    main()

max_stars_repo_name: chicocx/age-gender-estimation

max_issues_repo_name: chicocx/age-gender-estimation

max_forks_repo_name: chicocx/age-gender-estimation


text:
import sympy as S
import mpmath as poly
import numpy as N

class NoRealRoots(Exception):
    pass

class ODE:
    """This class represents the ODEs and associated computation of time.

    """

    MAX_QUANTA = 10**-3
    NUM_TERMS = 5               # This should be adjustable

    TRIG_FUNCS = [S.sin, S.cos, S.tan, S.cot, S.sec, S.csc]
    INV_TRIG_FUNCS = [S.asin, S.acos, S.atan, S.acot, S.asec, S.acsc, S.atan2]
    HYPERBOLIC_FUNCS = [S.sinh, S.cosh, S.tanh, S.coth, S.sech, S.csch]
    INV_HYPERBOLIC_FUNCS = [S.asinh, S.acosh, S.atanh, S.acoth, S.asech,
                            S.acsch]
    EXP_LOG = [S.exp, S.ln]
    TRANSCEDENTAL_FUNCS = (TRIG_FUNCS + INV_TRIG_FUNCS + HYPERBOLIC_FUNCS +
                           INV_HYPERBOLIC_FUNCS + EXP_LOG)

    def __init__(self, env, lvalue, rvalue, qorder=1, torder=1,
                 iterations=20, vtol=0, ttol=10**-2, taylor_expand=5,
                 trans_funcs=[], simplify_poly=False, half=False):
        """The quantized state order and taylor series order by default is 1.
        The maximum number of back-stepping iterations is 20 be default.
        The tolerance by default is 10^-2. taylor_expand gives the
        number to terms that we expand transcendental function too,
        default 5. Simplify the polynomial before finding roots, can
        take a very long time. Usually simplification is needed if
        polynomial has both a numerator and a denominator.

        """
        self.env = env
        self.lvalue = lvalue
        self.rvalue = rvalue
        self.qorder = qorder
        self.torder = torder
        self.iterations = iterations
        self.vtol = vtol
        self.ttol = ttol
        ODE.NUM_TERMS = taylor_expand
        ODE.TRANSCEDENTAL_FUNCS += trans_funcs
        ODE.simplify_poly = simplify_poly
        self.half = half

    @staticmethod
    # XXX:
    def replace(expr, s, v):
        if expr == s:
            return (S.sympify(v))
        elif expr.args == ():
            return expr
        else:
            return expr.func(*[ODE.replace(a, s, v)
                               for a in expr.args])

    # XXX: Post order traversal
    @staticmethod
    def taylor_expand(expr, around=0):
        if expr.args is ():
            return expr
        args = [ODE.taylor_expand(a) for a in expr.args]
        if expr.func in ODE.TRANSCEDENTAL_FUNCS:
            if len(args) != 1:
                raise RuntimeError('Cannot create a taylor series '
                                   'approximation of: ', expr)
            else:
                # XXX: Build the polynomial for arg
                coeffs = poly.taylor(expr.func, around, ODE.NUM_TERMS)
                # print(coeffs)
                coeffs = [(S.Mul(float(a), S.Mul(*[args[0]
                                                   for i in range(c)])))
                          for c, a in enumerate(coeffs)][::-1]
                # print(coeffs)
                return S.Add(*coeffs)
        else:
            return expr.func(*args)

    def compute(self, init, time):
        # Now init is a dictionary of all the required initial values.
        slope = self.rvalue
        for k in init:
            slope = ODE.replace(slope, k, init[k])
        slope = slope.subs('t', time)
        return init[self.lvalue.args[0]] + float(slope)*time

    def _delta1(self, init):
        return init[self.lvalue.args[0]]

    def _delta2(self, init):
        # slope = ODE.replace(self.rvalue, self.lvalue.args[0], init)
        slope = self.rvalue
        for k in init:
            slope = ODE.replace(slope, k, init[k])
        t = S.Symbol('t')
        return (S.Add(init[self.lvalue.args[0]],
                      (S.Mul(slope, (t - self.env.now)))))

    def _taylor1(self, init, q, q2, quanta, count):
        def is_den(x):
            return (type(x) == S.Pow and x.args[1] == -1)

        def compute_delta(part_poly, d, dl, quanta):
            # XXX: Positive quantum, so polynomial - quanta = 0
            polynomial1 = S.Add(part_poly, -quanta)
            # XXX: Assumption that the time line is called "t"
            # print(polynomial1)
            if not ODE.simplify_poly:
                polynomial1 = (polynomial1.expand().subs('t', d))
            else:
                polynomial1 = S.simplify(polynomial1.expand().subs('t', d))
            ppoly = polynomial1
            # XXX: Taking care of numerator and denominators after
            # expansion.
            if type(polynomial1) == S.Mul:
                if not is_den(polynomial1.args[0]):
                    polynomial1 = polynomial1.args[0]  # Get just the numerator
                else:
                    polynomial1 = polynomial1.args[1]
            # print('polynomial:', polynomial1)
            # If "δ" vanishes after expansion then just return None
            if (type(polynomial1) is S.Float):
                return None
            polynomial1 = S.Poly(polynomial1)
            try:
                nsoln = N.roots(polynomial1.all_coeffs())
                nsoln = nsoln[N.isreal(nsoln)]
                nsoln = nsoln[N.where(nsoln >= 0)]
                # soln = poly.polyroots([poly.mpf(a) for
                #                        a in polynomial1.all_coeffs()])
                # print('1:', nsoln, soln)
            except S.PolynomialError as e:
                print('When trying to solve: ', ppoly)
                raise e
            # dl += [float(a) for a in soln
            #        if type(a) is poly.mpf and float(a) >= 0]
            dl += list(nsoln)
            # The second polynomial
            # XXX: Negative quantum, so polynomial + quanta = 0
            polynomial2 = S.Add(part_poly, quanta)
            # XXX: Assumption that the time line is called "t"
            if not ODE.simplify_poly:
                polynomial2 = (polynomial2.expand().subs('t', d))
            else:
                polynomial2 = S.simplify(polynomial2.expand().subs('t', d))
            ppoly = polynomial2
            # print(ppoly.args[0], ppoly.args[1])
            if type(polynomial2) == S.Mul:
                if not is_den(polynomial2.args[0]):
                    polynomial2 = polynomial2.args[0]  # Get just the numerator
                else:
                    polynomial2 = polynomial2.args[1]
            polynomial2 = S.poly(polynomial2)
            try:
                nsoln = N.roots(polynomial2.all_coeffs())
                nsoln = nsoln[N.isreal(nsoln)]
                nsoln = nsoln[N.where(nsoln >= 0)]
                # soln = poly.polyroots([poly.mpf(a) for
                #                        a in polynomial2.all_coeffs()])
                # print('2:', nsoln, soln)
            except S.PolynomialError as e:
                print('When trying to solve: ', ppoly)
                raise e
            # dl += [float(a) for a in soln
            #        if type(a) is poly.mpf and float(a) >= 0]
            dl += list(nsoln)
            return dl

        def get_d(q):
            d = S.Symbol('d', positive=True, real=True)
            # XXX: My rvalue can depend upon a whole vector os q's
            # TODO: Convert it into a taylor series
            # print(self.rvalue, q)

            # XXX: Making a taylor polynomial if it is transcendental
            # function
            slope = ODE.taylor_expand(self.rvalue)
            # print('slope: ', slope)
            # print(q)
            for k in q:
                slope = ODE.replace(slope, k, q[k]).evalf()
            # print(slope)
            # XXX: IMP CHANGE! Here I am chaning QSS to compare with a
            # constant level not qith "Q". Note that q is the slope
            # itself.
            part_poly = S.Mul(d, slope)
            # print('ppoly: ', part_poly.subs('t', 'd').expand().evalf())
            # XXX: compute_delta saolves for the roots of the polynomial
            dl = compute_delta(part_poly, d, [], quanta)
            if dl is None:
                return None     # The constant slope case
            elif dl == []:
                raise NoRealRoots('No real positive root for: ',
                                  S.Eq(part_poly.subs('t', d).expand(),
                                       quanta), '{:.2e}'.format(quanta))
            d = min(dl)
            return d

        # print('getting δ1')
        d1 = get_d(q)           # Get the future time event from QSS-1
        # print('getting δ2')
        d2 = get_d(q2)          # Get the future time event from modified QSS-2
        if d1 is None:
            return S.oo, quanta  # This is returning infinity, wrong HA
        if d2 is None:
            # d1s = '{:.2e}'.format(d1)
            # quanta = '{:.2e}'.format(quanta)
            # print('chosen Δq: %s δ: %s' % (quanta, d1s))
            return d1, quanta
        elif abs(d1 - d2) <= self.ttol:
            # d1s = '{:.2e}'.format(d1)
            # d2s = '{:.2e}'.format(d2)
            # pquanta = '{:.2e}'.format(quanta)
            # print('chosen Δq: %s δ1: %s δ2: %s' % (pquanta, d1s, d2s))
            # In this case we have satisfied εt so returning first δ
            return d1, quanta
        elif count < self.iterations:
            # If the delta step results in output that is within the
            # user defined error bounds then great. Else, half the
            # quanta and keep on doing this until number of iterations
            # is met. This is reducing the quanta in a geometric
            # progression.

            # XXX: Adaptive Stepsize Runge-Kutta Integration William H.
            # Press, and Saul A. Teukolsky
            newquanta = d1 * pow(abs(self.ttol / (d1 - d2)), 1.0/2)
            quanta = newquanta if newquanta <= quanta else 0.5*quanta
            return self._taylor1(init, q, q2, quanta, (count+1))
        else:
            raise RuntimeError('Could not find delta that satisfies '
                               'the user specified error bound: '
                               'ε: %s δ1: %s δ2: %s Q1: %s Q2: %s '
                               'Δq: %s. Increase interation count'
                               % (self.ttol, d1, d2, q, q2, quanta))

    def _taylor(self, init, q, q2, quanta):
        if self.torder == 1:    # First order taylor only supported
            # The delta step
            return self._taylor1(init, q, q2, quanta, 0)
        elif self.torder > 1:
            raise RuntimeError('Currently only first order taylor supported')

    def get_q(self, init, order):
        # First calculate the q(t) given the qorder
        if order == 1:
            q = self._delta1(init)
        elif order == 2:
            q = self._delta2(init)
        elif order > 2:
            raise RuntimeError('Curretly only upto QSS2 is supported')
        return q

    # XXX: This is the main function, which returns the future time
    # event per level crossing per variable.
    def delta(self, init, other_odes=None, quanta=MAX_QUANTA):
        """This is the main function that returns back the delta step-size.
        Arguments: The initial value of the ode. Returns: The delta
        step-size that is within the user specified error.

        """
        # These two are me XXX: Here we are building the quantized
        # states, i.e., hysterisis for qorder=1 and integration for
        # qorder-2.
        qs = {self.lvalue.args[0]: self.get_q(init, self.qorder)}
        q2s = {self.lvalue.args[0]: self.get_q(init, self.qorder+1)}
        # XXX: Building the quantized states for other odes that we
        # might depend upon, because we can have coupled ODEs.
        if other_odes is not None:
            for ode in other_odes:
                qs[ode.lvalue.args[0]] = ode.get_q(init, ode.qorder)
                q2s[ode.lvalue.args[0]] = ode.get_q(init, ode.qorder+1)
        # XXX: delta is the returned value
        delta, nquanta = self._taylor(init, qs, q2s, quanta)
        # XXX: Handling sudden jumps
        if self.half and abs(quanta) > ODE.MAX_QUANTA:
            delta, _ = self._taylor(init, qs, q2s, quanta*0.5)
        elif ((not (type(self.rvalue) is S.Float
                    or type(self.rvalue) is S.Integer) and
               float(nquanta) == quanta and abs(quanta) > ODE.MAX_QUANTA)):
            delta, _ = self._taylor(init, qs, q2s, quanta*0.5)
        return delta

    def __str__(self):
        ode = str(self.lvalue) + ' = ' + str(self.rvalue)
        ret = ' '.join([ode])
        return ret+'\n'

reason: import numpy,import sympy,import mpmath


text:
"""
model class is the augmented model class that expands on the pysd model class with additional information

- allows attachment of a full_df with all the runs
- allows run with tracking for run time error tracking

Version 0.2
Update 30.07.18/sk
"""

import pandas as pd
from timeit import default_timer as timer
import numpy as np
import traceback
from pysd.py_backend.functions import Model, Time
from configparser import ConfigParser
import os

class AugModel(Model):
    """
    Augmodel is built on top of the pysd model class to add functionality
    """

    def __init__(self, model_path):
        """
        Careful with this init... model.reload() is going to initialize the model again, so stuff that needs to be kept
        beyond the reload cannot be in the init function

        :param model_path:
        """
        super(Model, self).__init__(model_path, None, None)
        self.time = Time()
        self.time.stage = 'Load'
        self.initialize()

    def set_tracking(self):
        """
        This needs to be outside of init because model reload starts resets the model but we would like to keep
        the error list and the full df dict

        essentially all the init that shouldn't be reloaded

        :return:
        """
        self.cf = ConfigParser()
        # if this is run with MP, then the current working directory is the run pipe
        # if not, then it is ok
        self.cf.read(os.path.join(os.getcwd().replace('\\run_pipe',''), '_config', 'settings.ini'))
        # precision for rounding and finding negative flows and stocks
        # unlimited precision is sometimes generating negative values where there shouldn't be any
        self.precision = self.cf['model'].getint('round_precision', fallback=8)
        self.max_ts = self.cf['model'].getint('max_ts', fallback=10)
        # list for tracking errors
        self.err_lst = []
        # full df dict initialized
        self.full_df_dict = {}

    def start_model(self):
        """
        Starting the timer for the model

        Deprecated 27.07.18/sk
        Not used 27.07.18/sk

        :return:
        """
        self.start = timer()

    def end_model(self):
        """
        Ending the timer for the model

        Deprecated 27.07.18/sk
        Not used 27.07.18/sk

        :return:
        """
        self.end = timer()
        self.elapsed = (self.end - self.start) / 60

    def create_full_df(self, base, group):
        """
        Creating the full_df by group with the base run as initial point

        :param base: base run from the model
        :param group: identifier for the full_df dict
        :return:
        """
        self.full_df_dict[group] = pd.DataFrame()
        self.add_run(base, 'base', group)

    def add_run(self, new_df, name, group):
        """
        adds run to full.df, which is the Dataframe that contains all the endogenous variables with all the runs
        full df is also used then to create the heatmaps

        :param group: identifier for the full_df dict
        :param new_df: df with the new run (only endogenous variables)
        :param name: name of the new run to be added, is added as level 1 in the multi index of the full df
        :return: full_df with added new run
        """
        new_df = pd.concat([new_df], axis=1, keys=[name])
        new_df = new_df.astype('float64', copy=False)
        self.full_df_dict[group] = pd.concat([self.full_df_dict[group], new_df], axis=1)

    @staticmethod
    def clean_trace(trace_str):
        """
        cleans the traceback string to the relevant line and cleans it, returns location of runtime error

        input is:
        Traceback (most recent call last):
          File "C:\winprog\Anaconda3\lib\site-packages\pysd\py_backend\functions.py", line 71, in cached
            assert cached.cache_t == func.__globals__['time']()
        AttributeError: 'function' object has no attribute 'cache_t'
        During handling of the above exception, another exception occurred:
        Traceback (most recent call last):
          File "C:\code\testingbattery\tb\tb_backend\model.py", line 224, in run_with_tracking
            run.reload)
          File "C:\winprog\Anaconda3\lib\site-packages\pysd\py_backend\functions.py", line 648, in run
            res = self._integrate(t_series, capture_elements, return_timestamps)
          File "C:\winprog\Anaconda3\lib\site-packages\pysd\py_backend\functions.py", line 747, in _integrate
            outputs.append({key: getattr(self.components, key)() for key in capture_elements})
          File "C:\winprog\Anaconda3\lib\site-packages\pysd\py_backend\functions.py", line 747, in <dictcomp>
            outputs.append({key: getattr(self.components, key)() for key in capture_elements})
          File "C:\winprog\Anaconda3\lib\site-packages\pysd\py_backend\functions.py", line 75, in cached
            cached.cache_val = func(*args)
          File "C:\code\testingbattery\test180711\workforce_treated.py", line 293, in pressure_to_hire
            return (task_backlog() - target_backlog()) / target_backlog()
        ZeroDivisionError: float division by zero

        Output is:
        pressure_to_hire

        :param trace_str: error trace
        :return:
        """
        lines = trace_str.splitlines()
        i = -2
        while True:
            # checks until it finds the last file, and saves i
            if lines[i].replace('\t', '').replace(' ', '').startswith('File'):
                break
            else:
                i -= 1
        # returns the last element of the line before the return
        return lines[i].rsplit(' ', 1)[-1].strip()

    @staticmethod
    def clean_type(type_str):
        """
        cleans the type of error if it's based on a warning at runtime and returns the error type

        input is:
        <class 'ZeroDivisionError'>

        output is:
        ZeroDivisionError

        :param type_str:
        :return:
        """
        return str(type_str).split('\'')[1]

    @staticmethod
    def clean_desc(desc_str):
        """
        cleans the error description and returns the description when the error is a warning at runtime

        input:
        'invalid value encountered in double_scalars'

        output:
        'Invalid Value as Result'

        :param desc_str: cleans up the description for the error file
        :return:
        """

        desc_str = str(desc_str)
        if desc_str == 'invalid value encountered in double_scalars':
            desc_str = 'Invalid Value as Result'
        elif desc_str == 'divide by zero encountered in double_scalars':
            desc_str = 'Division by Zero'
        return desc_str

    def neg_flow_check(self, run, param_dict, flow_names, test_name, change=''):
        """
        checks negative flow values unless the flow name contains the word 'net'

        :param run: pd.Dataframe(), result of the model execution
        :param param_dict: dict, parameter settings for the run
        :param flow_names: list, flows in the model
        :param test_name: str, name of the test
        :param change: str, if only one variable is changed, the variable is reported here
        :return:
        """

        errors = []
        # as type added to avoid problems when a df is interpeted with an object type column
        flow_df = run[flow_names].astype(np.float64)
        # iterating over each flow in the df
        for column in flow_df:
            # ignores the flows that have net in the name
            if 'net' not in column.lower().split(' '):
                # reports time stamps at which this error occurs (index is timestamps for run)
                # we don't have to account for the np.nan here (unlike in the run comparison)
                # because np.nan inequalities are always FALSE
                neg_lst = flow_df.index[flow_df[column].round(self.precision) < 0].tolist()
                # only the first 10 time steps are reported, to keep the error file clean
                neg_lst = neg_lst[:self.max_ts]
                if neg_lst:
                    errors.append((test_name, 'Negative Flow', neg_lst, column, param_dict, change))
        return errors

    def neg_stock_check(self, run, param_dict, stock_names, test_name, change=''):
        """
        checks negative stock values, there is no keyword for stocks that could also go negative, should be declared

        :param run: pd.Dataframe(), run results
        :param param_dict: dict, parameter settings for the run
        :param stock_names: list, name of stocks in the model
        :param test_name: str, test the run originates from
        :param change: str, if only one variable is changed, it is reported here
        :return:
        """
        errors = []
        # as type added to avoid problems when a df is interpeted with an object type column
        stock_df = run[stock_names].astype(np.float64)
        # iterating over each stock
        for column in stock_df:
            # reports time stamps at which this error occurs (index is timestamps for run)
            # we don't have to account for the np.nan here (unlike in the run comparison)
            # because np.nan inequalities are always FALSE
            neg_lst = stock_df.index[stock_df[column].round(self.precision) < 0].tolist()
            # only the first 10 time steps are reported, to keep the error file clean
            neg_lst = neg_lst[:self.max_ts]
            if neg_lst:
                errors.append((test_name, 'Negative Stock', neg_lst, column, param_dict, change))
        return errors

    def run_with_tracking(self, args):
        """
        args = run,flow_names,stock_names,test_name

        returns a tuple with three elements
        name, res, errors

        tracks errors at runtime and reports those errors

        errors are logged in the error file with the information structure:
        source, error_type, description, location, parameter, change to base
        source is the name of the test and comes from the test.init functions
        error_type is created by this routine, currently includes
            division by 0, floating point, negative stock, negative flow
        description is explanatory text (division by 0, floating point) or
            the timesteps at which the error occurs (negative stock, negative flow)
        location is the variable at which the error occured
        parameter: is the dict with all parameter settings
        change to base: is used when only one variable is changed compared to base (for easier viewing)

        :param args: tuple with run info, flow names, stock names and test name
        :return: run results and updated error list
        """
        run, flow_names, stock_names, test_name = args
        # creating the empty data frame to store the run in
        res = pd.DataFrame()
        errors = []
        # the model is run with np.errorstate in order to change the warnings to exceptions
        # and being able to capture them
        # this will break the run, thus the timestamp for exception errors cannot be retrieved
        with np.errstate(divide='raise', invalid='raise'):
            try:
                res = self.run(run.input_dict, run.return_columns, run.return_timestamps, run.initial_condition,
                               run.reload)
            except Exception as e:
                errors.append(
                    (test_name, self.clean_type(type(e)), self.clean_desc(e), self.clean_trace(traceback.format_exc()),
                     run.input_dict, run.change))
        # if an exception occurs above, the run remains empty as execution is halted, thus the run has to be executed
        # again, but outside of np.errstate to just have warnings, that way the output is always a completed run
        # if a invalid value is hard coded (e.g. in the knockout analysis), then it still raises an error, so the try
        # block takes care of that
        if res.empty:
            try:
                res = self.run(run.input_dict, run.return_columns, run.return_timestamps, run.initial_condition,
                               run.reload)
                # then the negative flows and stock checks are run on the completed run
            except Exception:
                pass
        # the testing needs to be checked at the end
        if not res.empty:
            nf_errors = self.neg_flow_check(res, run.input_dict, flow_names, test_name, run.change)
            ns_errors = self.neg_stock_check(res, run.input_dict, stock_names, test_name, run.change)
            errors.extend(nf_errors)
            errors.extend(ns_errors)
        # this has to be fixed, the returning element should be only the run object with the information attached
        # also then the run check could be done here instead of in collect results (where there's no MP) 250718/sk
        return run.name, res, errors

    def add_errors(self, errors):
        """
        Adds errors to the error list

        Deprecated 27.07.18/sk
        Not used 27.07.18/sk

        :param errors:
        :return:
        """
        self.err_lst.extend(errors)

    def check_err_list(self, source):
        """
        Error list check to make sure that at least one entry happens to avoid problems when converting to DF

        :param source: test where the check originates from
        :return:
        """
        try:
            if not self.err_lst:
                self.err_lst.append((source, 'No Error', '', '', '', ''))
        except AttributeError:
            print('run time error')


text:
import numpy as np

def project_data(X, U, K):
    # You need to return the following variables correctly.
    Z = np.zeros((X.shape[0], K))

    # ===================== Your Code Here =====================
    # Instructions: Compute the projection of the data using only the top K
    #               eigenvectors in U (first K columns).
    #               For the i-th example X[i], the projection on to the k-th
    #               eigenvector is given as follows:
    #                   x = X(i, :)';
    #                   projection_k = x' * U(:, k);
    #                   (above is octave code)
    #
    Ureduce = U[:, np.arange(K)]
    Z = np.dot(X, Ureduce)

    # ==========================================================

    return Z

max_stars_repo_name: RunningGump/coursera-ml-py

max_issues_repo_name: RunningGump/coursera-ml-py

max_forks_repo_name: RunningGump/coursera-ml-py