999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

和聲搜索的半監督聚類研究與應用

2012-07-25 11:06:36王華秋
計算機工程與設計 2012年7期

王華秋

(重慶理工大學 計算機科學與工程學院,重慶400054)

0 引 言

由于計算資源有限,對于應用程序開發者而言,計算性能是一項關注點。開發者總是致力于如何寫出高效的科學計算代碼以最優地利用超級計算機或其它的計算資源,然而,他們常常由于缺乏一些必要的信息而無法讓代碼最優化,這些信息包括:體系結構、內存機制、處理器網格等。這些復雜而隱藏的因素決定了程序的性能,比如:內存使用、I/O、編譯器、操作系統等等。

用于高性能計算機的性能分析工具就可以幫助使用者寫出高效的代碼。這些工具能夠提供程序性能評測值,并能發現一些瓶頸。瓶頸是指在程序執行時,運行時間由于無效資源的占用而浪費掉了。只要識別出瓶頸,開發者就可以修改程序以提高運行效率。正如現如今有許多并行計算硬件和并行編程語言一樣,性能分析工具也是層出不窮,其中有代表性的是如下幾種:Paradyn,TAU,Periscope,Vampir,KOJAK,SCALASCA,mpiP。

聚類機制為這些工具提供了一種基本的多變量統計分析的方法,所有處理器的Performance和Region這兩個屬性對應的Serverity值被提取出來,組成一個由處理器個數為行,Serverity為列的高維矩陣。我們用聚類算法對這些處理器進行分類,以檢測每一類處理器在不同代碼區域出現性能瓶頸的概率分布,以便對每一類處理器進行合理的資源調度。用于GENE代碼分析的處理器有1024個,這些處理器產生的性能瓶頸維度有39列之多,這樣會產生1024*39的矩陣,如何評價這些處理器的性能,因此對這樣的性能矩陣進行劃分,使其產生少數具有相似性能的處理器集合將是非常有意義的,聚類后的處理器將有利于發現隱藏在巨大數據集后的并行系統整體運行趨勢。

如今,存在著諸如劃分、層次、密度、網格等聚類算法。劃分聚類是最常用的一種方法,也是許多啟發式搜索算法首選的優化目標,這些啟發式搜索算法已廣泛地用于聚類中心的發現中,這些算法能在問題域內找到所有可能解,而且能避免陷入局部最優解。這些算法例如克隆選擇[1]、蟻群 優 化[2]和 粒 子 群 優 化[3]。和 聲 算 法 是 Geem 在2001年提出一種由音樂靈感激發的優化算法[4],受到廣泛關注,和聲搜索算法已經用于許多工程和科學領域,比如:圖像分割[5]、Web文本聚類[6]、水處理[7]。因此本文也用和聲搜索算法實現了快速準確地聚類,無論數據分布如何、形狀如何、大小如何,性能好的搜索劃分算法都可以對其聚類。

眾所周知,每個都數據包含一定數量的維度 (或者稱為屬性),聚類算法就是將這些多維數據進行劃分,在劃分之前并不知道需要劃分成為幾部分,只是將每一類內部變量盡量聚集在一起,而類與類之間的距離盡量遠離,兩點之間的差異用相似度或者距離計算。聚類的這種不確定性正是在于聚類學習是無監督的過程。

基本的劃分算法認為所有的維對于聚類同等重要,事實上,一些維存在冗余、無關、甚至易產生誤導的異常點。因此特征選擇技術出現了,這是一種用于高維數據聚類的最佳特征子集選擇技術[8]。

特征選擇不僅提高了聚類的性能,而且對于簡化了聚類模型。特征選擇一般是在聚類之前完成的過程,但是也有一些算法將其集成于聚類算法中,即在聚類的同時進行特征選擇,文獻 [8-9]就是將特征選擇和EM聚類結合在一起,在其影響下,文獻 [10]提出了一種利用和聲搜索算法的聚類和特征選擇同步進行的方法,該方法將特征選擇作為和聲搜索聚類過程的一部分,每次搜索的目標不僅包括聚類中心,還包括被選擇出來的屬性 (維),而且為了獲得最佳的聚類中心數量,每次搜索需要不斷改變聚類的中心數量。這種集成模式的搜索采用的評價函數則是聚類誤差和屬性 (維)相對中心的偏差兩部分之和。這些改進的確提高了劃分聚類的性能,但是同時帶來了一些問題:①算法復雜度較高,隨機的聚類中心和隨機的特征屬性的組合將是非常龐大的,這對于高維度大數據而言,其算法復雜度是可想而知的;②性能評價函數不太合理,將聚類誤差和屬性偏差結合在一起不利于評價總誤差究竟是由于聚類中心數量不適當而帶來的誤差,還是由于屬性特征選擇不合理而帶來的偏差,從而造成每次迭代無法正確調整3個目標值:聚類中心數、聚類中心、屬性維度。這將對算法的準確度產生不可忽視的影響。

以上的做法是對高維數據的屬性 (特征)進行選擇,沒有被選上的屬性將被忽略或刪除,進一步,有的學者提出屬性加權的思想對特征選擇進行泛化,即是為每一維設置一個0到1之間的權值,而不是直接保留或刪除某些屬性 (維)。換句話說,如果某一屬性的特征比較明顯,那么它的權值就較高,從而對聚類結果的影響就較大;反之,特征不顯著,則權值較小,而影響較小。這樣看來,特征選擇就是為每一維設置合理的權值。近年來,許多學者致力于研究這種特征加權的方法以提高聚類質量。

Huang[11]提出了一種用于Kmeans的特征屬性自動加權的機制,他們利用迭代公式不斷修正權值,而Kmeans算法的有效性并沒有被影響,然而,該方法需要使用者主觀指定一些參數值,而且,該方法獲得的權值無法顯著區分屬性是否具有代表性和不相關性。

文獻 [12]提出了用優化的方法對屬性進行加權,首先測量不同特征屬性的分散度,然后構造了凸Kmeans算法以決定最優權值,同時保證簇內部平均分散度最小以及簇之間的平均分散度最大,在實踐中,對于所有屬性進行加權并不適合高維數據聚類。文獻 [13]提出了COSA方法用于對子屬性進行層次聚類。但是求解這個優化函數需要面臨巨大的算法復雜度。受到這樣的啟發,一些學者將這種特征加權的優化方法用于子空間聚類[14-15]。

根據以往研究的經驗,結合需要解決的問題,本文提出了一種采用了和聲搜索的半監督聚類算法,該算法克服了文獻 [10]提出的算法的那兩點問題,而且采用了一種簡單可行的特征屬性加權的方法是區別各個維在聚類中的不同作用,而不是忽略某些屬性維,采用評價函數對不同的聚類中心個數進行評價,從而保證算法收斂時能夠得到最佳聚類中心。最后本文將這種算法用于解決并行計算中的處理器聚類問題。

1 算法設計

按照引言的論述,我們要解決的問題有三點:①特征屬性的加權,即數據維加權;②采用和聲搜索進行數據聚類;③選擇最優的聚類中心數。因此我們把算法分為三步完成:

(1)在聚類之前先進行特征屬性的加權。

任何一維的數據如果過于集中,數據之間失去的區別,就不利于聚類了,這樣的維所代表的屬性就應該具有較小的權值;反之,如果這一維數據比較散亂,數據之間便于區分,那么這一維對聚類結果影響較大,應該分配較大的權值。本文就是按照這樣的思路對每一維屬性進行加權

這樣wi就代表了第i維對于聚類的權值,這里stdi是第i維數據的標準偏差,meani是第i維數據的平均值,stdi標準偏差代表了每一維的變化性,用標準偏差除以平均值則可以將每一維的數據分散性歸一化,使得所有維之間的權值具有可比性。

(2)利用和聲搜索算法計算聚類中心。

聚類評價函數值作為和聲算法的適應度,并采用一種新的調音概率公式改進和聲算法,不斷調整和聲庫直到算法收斂,直至聚類評價函數值的改變量小于要求或者達到迭代次數。

(3)選擇最優的聚類中心數。

不斷改變聚類中心數,利用第二步計算其聚類中心和評價函數值,選擇評價函數值最小對應的聚類中心作為最后的輸出。

由此可以看出,第二步事實上是本算法的關鍵,下面我們詳細介紹這一步計算過程。

2 和聲搜索聚類算法研究

和Kmeans、Kmediods、FCM等算法一樣,和聲搜索聚類算法同樣離不開聚類中心反復 “迭代”,只不過這時的“迭代”演變成了 “演奏”。

2.1 初始化和設置參數

現已知樣本

式中:n——樣本數量,d——樣本的維度,每一維數據的上界 Uxj∈d和下界 Lxj∈d。

(1)數據歸一化:由于各維數據代表的意義不同,造成量綱不統一,這樣不利于計算彼此的距離。由于在屬性加權時我們

已經對每一維的數據計算了標準差和平均值了,這里就采用Z-Score歸一化方法,由下式計算

(2)聚類需要確定的參數是:我們需要定義聚類質量評價函數,這個評價應該體現聚類的原則:簇內樣本盡量靠近該簇中心,簇和另一個簇盡量遠離。因此我們把盡量小的簇內距離除以盡量大的簇間距離作為評價函數 (也稱為適應度)

式中:C——聚類中心的集合,k——聚類中心的個數,Ci和Cj分是第i個和第j個聚類中心,nj是屬于第j類的樣本數量。這樣我們聚類的目標就是通過和聲搜索到使得f(C,k)達到最小的C和k,C用向量的方式存儲,C=(C1,C2,…,Ck)= (c11,…,c1d,…,ck1,…,ckd),例如:第一個聚類中心就是c11,c12,…,c1d,第k個聚類中心就是ck1,ck2,…,ckd。

(3)和聲搜索算法也有一些參數需要確定:和聲存儲庫的行數 (HMS=20),和聲存儲庫的列數 (COL=k*d),和聲存儲庫的選擇概率 (HMCR=0.9),調音概率(PAR=0.25),隨機帶寬 (bw=0.2),音演奏次數 (即迭代次數ITERATION=100)或者停止條件 (適應度變化量小于一個小數DELTA<0.01)。

2.2 初始化和聲存儲庫

在這一步中,和聲存儲庫HM矩陣最初由有上界Uxj∈d和下界 Lxj∈d限制的隨機數組成

調用聚類質量評價函數計算各行和音的適應度,得到適應度列向量F= [f1,f2,…,fHMS]T,記錄其中最佳和最差的適應度及其對應的和聲變量。

和聲存儲庫的作用就是在每次迭代時,由此庫產生出新的和聲變量,并根據這新變量不斷更新和聲存儲庫,使之收斂到全局最優解。

2.3 演奏新的和音

一個新的和音向量

這里,1≤i≤k,1≤j≤d,這個新和音向量有3種途徑產生:①和音存儲庫中選擇;②從和音庫中選擇后調音;③在允許范圍內隨機產生。產生新的和音的過程稱之為 “演奏”。

在第一種方法中,第一個類中心的第一維變量c′11將從和音存儲庫的第一列向量中隨機選擇,以此類推可以得到其它中心的各維變量cij’。

HMCR的變化范圍是0到1,一般為0.6~0.9之間,本文為0.9,這是決定是否從和音存儲庫HM中選擇性的變量的概率,換句話說,(1-HMCR)就是決定是否從上下限范圍內隨機產生新的變量的概率。選擇操作如下

對于每一個從和聲存儲庫中得到的變量還要決定是否調音,PAR就是初始調音的概率,一般為0.1~0.4之間,本文為0.25,概率 (1-PAR)就是不需要調音的概率。調音操作如下

所謂調音,就是在原來的變量的基礎上增加或減少一個隨機帶寬量bw

上述3種方法產生的新的和音將被返回并用于更新和音存儲庫。

調音是避免和聲算法出現 “早熟”的重要措施,在最初的和聲算法中,并沒有對調音概率PAR和隨機帶寬bw做調整,因此需要很多次迭代才能調出一個不同于其它的和聲。當庫中的和聲趨同時,我們以一個遠大于通常的調音概率PAR執行一次操作,同時增加隨機帶寬量bw,這樣就能夠隨機產生新的音符,從而使整個和音庫擺脫 “早熟”。

當和聲庫中的最佳適應度fbest與平均適應度favg滿足如下關系

我們就認為和聲庫具有趨同的趨勢,就適當提高調音概率PAR和隨機帶寬bw

這樣的改進可以增加算法的隨機搜索的幾率,提高算法的全局收斂性,但是考慮到收斂時間會延長,因此PAR和bw的增加系數均不太大。

2.4 更新和音存儲庫

使用新和音向量c′=cij′= (c11′,…,c1d′,…,ck1′,…,ckd′)作為聚類中心,調用聚類質量評價函數計算c′的適應度f′,如果新的和音向量c′在適應度上優于和音存儲庫中的最差的向量,那么和音存儲庫就需要更新了,新的和音向量將取代那個最差的向量。如果這個新的向量同時也是和音存儲庫中最優的,那么保存這個新向量的適應度為和音存儲庫的最佳適應度,這時的最佳適應度對應的變量就是最終聚集成為k個中心時的聚類結果。

2.5 檢查終止條件

如果滿足終止條件,結束聚類中心的搜索,否則,轉到3.3節。終止條件就是:①迭代次數達到最大值;②評價函數值 (適應度)的變化量小于限定值。

3 算法復雜度分析

為了分析整個算法的復雜性,我們介紹一下整個算法的流程,如圖1所示。

圖1 算法流程

我們把整個算法分為了若干段,每一段完成一定的功能,也消耗一定的時間,外部循環讓這些程序段反復迭代直至收斂。從內部開始分析,和聲搜索聚類中心需要的時間為:R(初始化)+k(中心數)×d(維度)×S(式 (4)的演奏+式 (5)和式 (6)的調音)+F(式 (3)的適應度計算n×k)+U (更新和聲庫)+P(式 (7)和式 (8)的大概率調音),由于有了尋找最佳聚類中心數的半監督過程,該算法重復了次和聲搜索聚類,另外還有式 (1)的維度加權和式(2)數據歸一化,因此對于整個半監督聚類而言,整個算法復雜度為。除了適應度公式之外,和聲搜索的這些公式都不太復雜,計算速度可以得到保證。

PSO聚類的空間復雜度為O (k*C* (N+1)),時間復雜度為O (k*m*C*n),總體復雜度為O (k*C*(N+1)+k*m*C*n)。其中k為粒子個數,C為需要類別數,N為數據維數,n為樣本數,m為各類樣本的均值向量數。由于FCM沒有采用多個粒子搜索聚類中心,因此其算法復雜度為O (C* (N+1)+m*C*n)。文獻[6]采用了和聲搜索和Kmeans的結合,文獻 [10]在和聲搜索聚類時考慮了特征選擇,因此二者的算法復雜度均為:O(W+* (R+k*d*S+k* (N+1)+m*k*n+n*k+U))。從理論上看,除FCM在復雜度上比本文提出的算法有優勢外,其余的算法復雜度都比本文的較高。

4 算法性能比較

我們把本文的HSW和FCM、PSO聚類、文獻 [6]和文獻 [10]算法進行比較,測試數據來自4種Benchmark數據,分別是:分為3類有150個4維數據的Iris,分為6類有871個3維數據的Vowel,分為2類有768個8維數據的Diabetes,分為3類有178個13維數據的Wine。由于上述幾種算法的性能均與初始值有關,因此我們每種算法分別測試10次,取其平均性能參數進行比較。

對于每種數據集測試時性能參數取:①測試中的最佳聚類中心的聚類評價函數值,②測試中的平均聚類評價函數值,③測試中的最差聚類評價函數值,④獲得最佳聚類中心的平均收斂時間。

這次比較實驗在Intel Pentium Dual-Core CPU T2310 1.46GHz 1.46GHz,1G內存的計算機上運行,上述算法的代碼使用Eclipse 3.5.2的Java語言編寫。

Iris數據集測試與性能比較結果見表1。

由于Iris數據量最少,維度也低,被分為3類,因此各種算法聚類效果都比較理想。

Vowel數據集測試與性能比較結果見表2。

表1 Iris數據集測試與性能比較結果

表2 Vowel數據集測試與性能比較結果

由于Vowel數據量最多,但是維度不高,被分為6類,類之間的差別比較大,因此本文算法聚類效果也比較好。Diabetes數據集測試與性能比較結果見表3。

表3 Diabetes數據集測試與性能比較結果

由于Diabetes數據量比較大,維度比較高,卻僅分為2類,因此類之間的差別比較小,因此幾種算法聚類效果都不太好。

Wine數據集測試與性能比較結果見表4。

表4 Wine數據集測試與性能比較結果

Wine和Iris類似,數據量少,也分為3類,但維度最高,為計算增加了難度,各種算法聚類效果仍然比較可以接受。

從上面的比較可以看出,本文這種算法聚類質量比較高,而且耗時較少。雖然FCM算法速度更快,但是誤差較大。因此我們采用本文的算法進行并行處理器性能聚類。

5 并行處理器性能聚類分析

為了分析和識別大規模科學計算中各處理器性能瓶頸。我們這里用回旋電磁數值實驗代碼 (gyrokinetic electromagnetic numerical experiment,GENE)進行并行計算,GENE代碼通過迭代求解具有5維相空間的非線性回旋方程,用于識別電磁等離子聚變產生的渦流特征。GENE可以運行在擁有大量處理器的超級計算機上,代碼用MPI并行化的Fortran95實現。數據采集實驗在ALTIX超級計算機上進行,GENE代碼其中的1024核中運行。通過測試可以收集GENE的不同代碼區域的各種MPI屬性值和可伸縮性,比如:節點性能,負載能力,信息組織方式,等等。示例數據如下所示:

<property cluster="false"ID="7-511">

<name>Excessive MPI communication time</name>

<context FileID="7"FileName="velo.f"RFL="33"Config="256x1"Region="CALL_REGION"RegionId="7-556">

<execObj process="247"thread="0"/>

</context>

<severity>54.36</severity>

<confidence>1.00</confidence>

<addInfo>

<CallTime>60866</CallTime>

<PhaseTime>111959</PhaseTime>

</addInfo>

</property>

<property cluster="false"ID="1-17">

<name>IA64Pipeline Stall Cycles</name>

<context FileID="5"FileName="main.f"RFL="63"Config="256x1"Region="USER_REGION"RegionId="5-584">

<execObj process="29"thread="0"/>

</context>

<severity>40.40</severity>

<confidence>1.00</confidence>

<addInfo>

</addInfo>

</property>

我們首先將這些數據進行分析,對于處理器有其編號process唯一確定,由于ID唯一表示了該CPU出現性能瓶頸的名稱,而RegionID則唯一標示了GENE代碼中出現性能瓶頸的區域位置,ID和RegionID的結合就可以確定GENE的哪一部分代碼運行時出現了什么樣的性能問題,而嚴重程度就是serverity。因此上述XML文件中,需要提取出來進行用于聚類的屬性為:process,ID,RegionID,serverity。我們對這些屬性進行構造,就是由給定的屬性構造和添加新的屬性,以有利于聚類。比如,我們根據屬性ID和RegionID可以構造Performance屬性。通過這種組合屬性,屬性構造可以發現關于數據屬性間聯系的丟失信息,這對知識發現有用的。于是我們構建了如下1024行39列的數據矩陣。

性能分析數據見表5。

表5 性能分析數據

利用這些數據,我們運用HSC算法對這1024個處理器進行聚類分析,將瓶頸性能近似的處理器歸為同一類,聚類的結果就是這1024個處理器被自動分為8類。

并行性能聚類結果見表6。

表6 并行性能聚類結果

圖2中的橫坐標是處理器的編號,縱坐標是類別。根據這樣的聚類,可以重新分配GENG的代碼區域到處理器上,取得更好的計算性能。

例如:p1,p20,p140,p230,p663,p915,…,被歸為一類,這一類的中心按照Serverity的大小排序為:54.36,52.22,43.36,43.17,40.71,20.34,18.46,7.71,…,按照Serverity的大小排序,這一類主要的瓶頸問題就是:

(1)Serverity值54.36:CALL_REGION 程序段 (7-556)出現了Excessive MPI communication time性能問題(7-511);

(2)Serverity值52.22:CALL_REGION 程序段 (7-556)出現了Excessive MPI time in receive due to late sender性能問題 (7-518);

圖2 處理器聚類結果

以此類推。有了這些性能特征,程序管理員或開發者就可以針對這任何一類處理器的任務重新做出分配或調整。

6 結束語

將隨機啟發式搜索算法之一的和聲搜索算法運用于聚類中,僅僅保留了聚類的評價函數,改變了聚類的方式,從而提高了聚類的質量,本文嘗試對其做了一些改進,和其它的算法相比,本文的HSW的半監督聚類有如下優點:

(1)通過半監督學習,算法能夠自動找到正確或合理的聚類中心數;

(2)通過特征屬性加權,改變了每一維對聚類結果的影響,提高了聚類準確性;

(3)在算法運行后期采用大概率調音保證了算法不會收斂到局部極小。

進一步的工作:這種算法畢竟是一種隨機啟發式搜索最優聚類中心算法,對于大數據集或高維數據聚類而言,收斂速度始終是一個問題,是否可以將代表不同聚類中心數的和聲庫分配給不同的進程或線程分布式計算其聚類中心?這種啟發式隨機搜索算法具有天然并行性,因此本文認為分布式和聲搜索聚類算法是可行的,必然提高收斂速度。

[1]LUO Yin-sheng,LI Ren-hou,ZHANG Wei-xi.Clustering algorithm based on clone selection theory [J].Control and Decision,2006,20 (11):1261-1264 (in Chinese). [羅印升,李人厚,張維璽.一種基于克隆選擇的聚類算法 [J].控制與決策,2005,20 (11):1261-1264.]

[2]XU Xiaohua,CHEN Ling.An adaptive ant clustering algorithm [J].Journal of Software,2006,17 (9):1884-1889(in Chinese). [徐曉華,陳崚.一種自適應的螞蟻聚類算法[J].軟件學報,2006,17 (9):1884-1889.]

[3]LI Shuai,WANG Xin-jun,GAO Dan-dan.PSO clustering algorithm based on internal spatial characteristic [J].Computer Engineering,2009,35 (5):197-199 (in Chinese). [李帥,王新軍,高丹丹.基于內部空間特性的PSO聚類算法 [J].計算機工程,2009,35 (5):197-199.]

[4]Santos Coelho L d,de Andrade Bernert D L.An improved harmony search algorithm for synchronization of discrete-time chaotic systems[J].Chaos,Solitons and Fractals,2009,41(15):2526-2532.

[5]Osama Alia,Rajeswari Mandava,Mohd Aziz.A hybrid harmony search algorithm for MRI brain segmentation [J].Evolutionary Intelligence,2011,4 (1):31-49.

[6]Forsati R,Meybodi M R,Mahdavi M,et al.Hybridization of K-means and harmony search methods for web page clustering[C].International Conference on Web Intelligence and Intelligent Agent Technology,Sydney,Australia:IEEE Computer Society,2008:329-335.

[7]Ayvaz M T.Simultaneous determination of aquifer parameters and zone structures with fuzzy c-means clustering and meta-heuristic harmony search algorithm [J].Advances in Water Resources,2007,30 (11):2326-2338.

[8]Weiguo S,Xiaohui L,Fairhurst M.A niching memetic algorithm for simultaneous clustering and feature selection [J].IEEE Trans on Knowledge and Data Engineering,2008,20(7):868-879.

[9]Zeng H,Cheung Y M.A new feature selection method for Gaussian mixture clustering [J].Pattern Recognition,2009,42 (2):243-250.

[10]Sarvari H,Khairdoost N,Fetanat A.Harmony search algorithm for simultaneous clustering and feature selection [C].International Conference of Soft Computing and Pattern Recognition.Paris, France:IEEE Computer Society,2010:202-207.

[11]Huang J Z,Ng M K,Rong H,et al.Automated variable weighting in k-means type clustering [J].IEEE Trans Pattern Anal Mach Intell,2005,27 (5):657-668.

[12]Yu Z,Wong H S.Genetic based k-means algorithm for selection of feature variables [C].International Conference on Pattern Recognition.Hong Kong,China:IEEE Computer Society,2006:744-747.

[13]Anil Kumar Tiwari,Lokesh Kumar Sharma,Rama Krishna G.Entropy weighting genetic k-means algorithm for subspace clustering [J].International Journal of Computer Applications,2010,7 (7):27-30.

[14]Jing L,Ng M K,Huang J Z.An entropy weighting k-means algorithm for subspace clustering of high dimensional sparse data [J].IEEE Trans on Knowledge and Data Engineering,2007,19 (8):1026-1041.

[15]Domeniconi C,Gunopulos D,Ma S,et al.Locally adaptive metrics for clustering high dimensional data [J].Data Mining and Knowledge Discovery,2007,14 (1):63-97.

主站蜘蛛池模板: 国产情侣一区| 青青操国产| 激情综合图区| 毛片免费视频| 亚洲欧美日韩成人高清在线一区| 91小视频在线| a天堂视频| 狠狠色丁香婷婷| 久久精品人人做人人综合试看| 亚洲乱码在线视频| 99久久国产精品无码| 国产精品自在线拍国产电影| 国产欧美视频综合二区| 怡春院欧美一区二区三区免费| 伊人成人在线视频| 日韩精品资源| 欧美成人a∨视频免费观看| 高清无码不卡视频| 97精品国产高清久久久久蜜芽| 91精品国产综合久久香蕉922| 特级欧美视频aaaaaa| 免费高清毛片| 99ri精品视频在线观看播放| 精品国产香蕉伊思人在线| 欧美日本不卡| 鲁鲁鲁爽爽爽在线视频观看| 波多野结衣无码AV在线| 999精品色在线观看| 国产1区2区在线观看| 国内毛片视频| 无码精品国产dvd在线观看9久| 日本欧美一二三区色视频| 成人福利视频网| 国产资源站| 一本一道波多野结衣一区二区| 欧美激情网址| 日本免费a视频| 国产极品嫩模在线观看91| 欧美日本激情| 国产精品久久精品| 国内精品一区二区在线观看| 操国产美女| 国产亚洲欧美另类一区二区| 日韩欧美中文字幕在线精品| 国产精品久久久久久久伊一| 99精品伊人久久久大香线蕉| 九九线精品视频在线观看| 国产亚洲高清在线精品99| 国产成人盗摄精品| 在线观看亚洲精品福利片| 国产99精品视频| 无码高清专区| 亚洲一区二区黄色| 亚洲中文无码av永久伊人| 亚洲AV无码久久天堂| 久久成人国产精品免费软件| 久久久国产精品免费视频| 久久久久久久蜜桃| 欧美一区二区自偷自拍视频| 国产探花在线视频| 亚洲无卡视频| 香蕉99国内自产自拍视频| 夜夜操国产| 亚洲热线99精品视频| 一级毛片在线播放免费观看| 国产成人av一区二区三区| 天堂网国产| 成年女人a毛片免费视频| 久久久黄色片| 亚洲综合经典在线一区二区| 97在线碰| 中文字幕不卡免费高清视频| 亚洲男人在线天堂| 国内熟女少妇一线天| 在线观看的黄网| 老司机午夜精品网站在线观看| 91久久夜色精品| 欧美成人免费一区在线播放| 国产亚洲精| 97亚洲色综久久精品| 乱人伦视频中文字幕在线| 国产毛片高清一级国语|