針對Ｃ語言的自動微分系統及其應用

(中國科學院 a.軟件研究所 并行計算實驗室；b.研究生院；c.計算機網絡信息中心， 北京 100190)



摘 要：采用切線性模式和代碼轉換策略，開發了C語言自動微分轉換系統（DTC），用于牛頓法求解非線性方程中Jacobi矩陣—向量乘積計算。介紹系統計算模型、功能、特色，并討論系統的設計與實現技術，包括編譯技術、微分代碼轉換及輸入/輸出（I/O）相關分析。最后給出了幾個具有說服力的測試與應用。

關鍵詞：自動微分；切線性模式；Jacobi矩陣

中圖分類號：TP312 文獻標志碼：A

文章編號：10013695(2009)01015504



Differentiation transforming system in C and its applications

ZHANG Chunhuia，b， CHENG Qiangc， CAO Jianwena

(a.Parallel Computing Laboratory， Institute of Software， b.Graduate School， c.Computer Network Information Center， Chinese Academy of Sciences， Beijing 100190， China)

Abstract:This paper employed the tangent linear model and the strategy of sourcetosource transformation， developed differentiation transforming system in C (DTC)to calculate the Jacobianvector product in the solution of nonlinear equations with Newton method. Then introduced the calculation model， functionality， and discussed features of the system， and the techniques of design and implementation of the system， including compiler technology， differential code transformation， and input/output (I/O) dependence analysis. Last presented some persuasive testing and applications.

Key words:automatic differentiation; tangent linear model; Jacobian



0 引言

在大規?？茖W計算的很多領域，如參數敏感性分析、變分資料同化、最優化設計和反問題等，需要計算函數不同形式的導數。自動微分(automatic differentiation，AD)[1]是一種計算函數導數的方法，它基于鏈式求導法則，不存在差分方法的截斷誤差，且應用靈活、開發代價相對較小。任何程序無論多么復雜，都可以分解為一系列順序執行的基本計算單元，如運算操作(加法、乘法等)和運算函數(三角函數等)，而形成從程序輸入變量變換到輸出變量的執行路徑。自動微分方法即是反復應用鏈式法則，逐次求出每個計算單元的局部偏導數，而最終得到所需的函數偏導數。自動微分方法計算函數導數的兩種基本模式分別為正向模式(forward mode)和反向模式(reverse mode)。其中，正向模式是按照從輸入變量到輸出變量的方向應用鏈式法則逐次計算每個計算單元的局部偏導數；而反向模式則是按相反的從輸出變量到輸入變量方向逐次計算的。兩種計算模式在不同的問題中各有優勢。例如運用伴隨模式(典型的一階反向模式)可以原計算常數倍的時間代價計算函數梯度，相應存儲代價較大；切線性模式(典型的一階正向模式)則適用于計算函數的Jacobi矩陣—向量乘積等，相應存儲代價較小。

有兩種基本的自動微分軟件實現策略，分別為算子重載(operator overloading)和代碼轉換(source transformation)。其中，算子重載策略是重載計算程序中的算術運算符及內部數學函數，其重載函數執行原運算操作及相應的局部偏導數計算，此方法適用于支持算子重載的語言，如C++和Fortran 90等；代碼轉換策略則是利用編譯技術，將原計算程序轉換為新的計算函數導數的程序代碼，它可用于任何程序語言。運用算子重載策略的AD軟件典型代表有ADOLC[2]、FADBAD/TADIFF[3]等；而運用代碼轉換策略的AD軟件典型代表有ADIFOR[4]、ADIC[5]等。

開發針對C語言的自動微分工具最大困難在于處理C語言靈活的語法。C語言不同于Fortran77等程序語言，它的數據結構更加復雜。例如，指針是C語言的重要特色，它的應用十分廣泛，如與數組配合使用、動態分配內存等，但指針的運用很容易出錯，這需要在開發針對C語言的編譯工具時準確處理指針的各種用法。再比如，C語言中定義了結構體（structure）和共用體（union）等結構，結構體和共用體類型域表中的成員也可以為結構體變量、共用體變量或指針、數組等，也就是它們可以具有復雜的嵌套結構，這需要編譯器對這些結構的處理具有通用性。另外，C語言允許用戶用typedef聲明新的類型名來代替已有類型名，即用typedef定義類型。通常，實現C語言的編譯工具本身是件很困難的事，而現有較好的C語言編譯工具又非常龐大。針對C語言的自動微分工具不僅要作為有效處理C語言的編譯工具，而且要把更多關注點放在將編譯工具適應于自動微分算法特點、簡潔高效地完成微分代碼轉換上面。所以，如何實現結構小巧卻能正確處理絕大部分C語言語法，并且更加適合自動微分代碼轉換的AD工具是開發者的一項挑戰。

開發采取切線性模式的自動微分工具的困難在于確保所生成的切線性模式微分代碼的可靠性。在大規模計算中，切線性模式微分代碼的可靠性不能保證，所以需要檢測切線性微分代碼的計算精度。切線性模式計算精度的檢測需要知道程序中定義的函數，其各個參數是輸入變量、輸出變量或是既為輸入又為輸出的變量（輸入/輸出變量）。這個函數參數的I/O信息對于準確生成測試程序十分關鍵。雖然數學上所有變量（無論是自變量還是因變量）都可以被等同地視為輸入/輸出變量，但在程序中，將函數參數都作為輸入/輸出變量將帶來很大的內存開銷，特別是在大型計算中不可行。所以，判斷函數參數的輸入/輸出屬性，即I/O相關分析非常必要?，F有的AD工具能夠自動生成微分模式可靠性檢測程序，或者具有I/O相關分析功能的很少。這也是開發者在實現AD工具中的一項薄弱環節。

目前，針對C語言，采取代碼轉換策略實現的AD工具僅有美國Argonne國家實驗室開發的ADIC。ADIC采用一種形式的正向模式計算模型，適用于計算Jacobi矩陣或者Jacobi矩陣的一部分。ADIC的開發借助了語法分析器Sage++，其AIF和SparsLinC部件富有特色。

1 DTC系統介紹

1.1 計算模型

DTC計算函數一階導數采取切線性模式的計算模型。將程序分解為一系列基本運算的復合，形成從輸入變量到輸出變量的計算路徑。依據一階微分形式不變性，切線性模式沿著這條計算路徑順序地計算各個中間變量的微分，直到最終計算得到輸出變量的微分。若程序計算的函數為多元向量函數，則DTC的這一切線性模式執行后，可得到函數對應Jacobi矩陣與向量的乘積。

設程序計算的向量函數為F。其中，函數F：Rn|→Rm，X=（x1，x2，…，xn）T對應程序的輸入變量，Y=(y1，y2，…，ym)T對應程序的輸出變量。F對應的Jacobi矩陣F′存在。

切線性模式計算模型可表示為

dY=F′dX(1)

其中：dX=(dx1，dx2，…，dxn)T；dY=(dy1，dy2，…，dym)T。

可以看到，給定適當的dx1，dx2，…，dxn，即可直接求得程序計算的函數F對應的Jacobi矩陣—向量積。特別地，對于m=1的情況，即可計算出函數沿（dx1，dx2，…，dxn）方向的方向導數。

1.2 系統功能

1.2.1 生成切線性模式

DTC處理的基本程序對象是函數，其默認接口為TYPE Func(TYPE X_IN， TYPE X_OUT， TYPE X_IN_OUT， TYPE N)。其中：“Func”表示定義的函數名；“TYPE X_IN、TYPE X_OUT、TYPE X_IN_OUT、TYPE NDIM”表示函數參數列表；X_IN、X_OUT和X_IN_OUT分別表示一系列輸入變量、輸出變量和輸入/輸出變量；N表示一系列整型或字符型變量；TYPE表示對應的類型，如double型、int型、指針型、數組型、結構體類型等。相應地，DTC生成的切線性模式代碼函數接口為TYPE diff_Func(TYPE diff_X_IN， TYPE X_IN， TYPE diff_X_OUT， TYPE X_OUT， TYPE diff_X_IN_OUT， TYPE X_IN_OUT， TYPE N)。其中，每個浮點型參數緊跟在它對應的微分擾動之后，微分擾動用相應參數加上默認前綴“diff_”表示。

1.2.2 生成測試代碼

DTC可以自動生成檢測切線性模式微分代碼正確性的程序。對于切線性模式微分代碼正確性的檢測，準確的方法是檢測切線性模式計算的函數每個輸出變量在任意點梯度的每一項是否能很好地近似當增量的范數足夠小時差分方法計算的相應值。然而，由于當輸入變量與輸出變量數目很大時，時間與空間開銷很大，這樣的檢測方法在實際中不可行。DTC中采用的檢測方法是在隨機點檢測切線性模式微分計算結果與差分方法計算結果是否符合。

設程序計算的向量函數為F：Rn|→Rm，若它在給定點X0處可微，則Jacobi矩陣F′(X0)存在且

‖F(X0+ΔX)-F(X0)‖/‖F′(X0)ΔX‖→1，‖ΔX‖→0（2）

DTC用統計的方法檢測切線性模式正確性，其生成的相應測試代碼接口為TYPE check_diff_Func(TYPE X_IN， TYPE X_OUT， TYPE X_IN_OUT， TYPE N)。其中，函數參數與原函數參數是一致的。為了確定程序中計算F的函數，其參數分別對應于F的自變量還是因變量，需要確定函數參數的I/O屬性，即進行I/O相關分析。在微分代碼檢測程序中，調用原函數和切線性模式函數前后需要將輸入/輸出變量進行保存和恢復。

1.3 系統特色

1.3.1 靜態相關分析

DTC的靜態相關分析主要體現在I/O相關分析方面。DTC可以記錄當前函數所有參數的I/O屬性，通過深層遞歸I/O相關分析，計算得到函數內參數最終的I/O屬性。

DTC亦可以對程序中全局變量、函數參數及局部變量進行簡潔準確的數據類型分析。這是通過DTC的變量表結構及相關操作完成的。

1.3.2 結構化的微分代碼轉換

DTC采用清晰的程序結構實現微分代碼轉換，所生成的切線性模式代碼中，微分擾動變量和原變量、微分計算語句和原計算語句總是成對出現。切線性模式代碼保持了原模式代碼的結構，如變量定義、循環結構、復合結構等，即采取了語句到語句、結構到結構的微分代碼轉換。

下面是DTC代碼的例子。它包含了一個函數定義，且沒有其他函數調用。斜壓一層半線性化準地轉位勢渦度方程模式eng.c源程序代碼如下所示：

/* This model is the linear baroclinic vortex equation model and \"z\" is the stream function. */

/* Author: Yinglai Jia. Revised to C version by Chunhui Zhang */

#include 〈math.h〉

/* Input parameters:z，dx，dy Output parameters:sek*/

int ix=513，jy=76;

void eng(double*sek， double z[ix][jy]， double dx[jy]，double dy)

{

double u，v，ek;

int i， j;

*sek = 0.0;

for(j=2;j<=jy-1;j++)

for(i=2;i<=ix-1;i++)

{

u=-1.0*(z[i][j+1]-z[i][j-1])/2.0/dx[j];

v=(z[i+1][j]-z[i-1][j])/2.0/dy;

ek=0.5*(u*u+v*v);

*sek=*sek+ek;

}

*sek=*sek/(double)(jy-2)/(double)(ix-2);

return;

}

DTC系統生成的對應切線性程序diff_eng.c的代碼如下所示：

#include \"dtcmath.h\"

#include 〈math.h〉

int ix = 513， jy = 76;

void diff_eng(double *diff_sek， double *sek， double diff_z[ix][jy]， double z[ix][jy]， double diff_dx[jy]， double dx[jy]， doublediff_dy， double dy)

{

double diff_u， u， diff_v， v， diff_ek， ek;

int ix = 513， jy = 75;

int i， j;

{ *diff_sek = 0.0;

*sek = 0.0; }

for(j=2;j<=jy-1;j++)

for(i=2;i<=ix-1;i++)

{

{diff_u=(-1.0*(diff_z[i][j+1]-diff_z[i][j-1])/

2.0*dx[j]-diff_dx[j]*-1.0*(z[i][j+1]-z[i][j-1])/2.0) /(dx[j]*dx[j]);

u=-1.0*(z[i][j+1]-z[i][j-1])/2.0/dx[j]; }

{diff_v=((diff_z[i+1][j]-diff_z[i-1][j])/2.0*dy-diff_dy*(z[i+1][j]-z[i-1][j])/2.0)/(dy*dy);

v=(z[i+1][j]-z[i-1][j])/2.0/dy;}

{diff_ek=0.5*(2.0*u*diff_u+2.0*v*diff_v);

ek=0.5*(u*u+v*v);}

{ *diff_sek=*diff_sek+diff_ek;

*sek=*sek + ek;}

}

{*diff_sek=*diff_sek/(double)(jy-2)/(double)(ix-2);

*sek=*sek/(double)(jy-2)/(double)(ix-2);}

return;

}

從這個例子可以看到，DTC的切線性模式代碼生成是一個語句到語句、結構到結構的微分代碼轉換過程。具體地，函數參數表和函數局部變量聲明中，每個浮點型變量，即一些文獻中提到的活躍變量（active variable）都緊隨在其對應微分擾動變量定義或聲明之后，而每個賦值語句緊隨在其對應微分語句之后。

2 DTC的設計與實現

2.1 系統設計

DTC借助詞法分析工具Flex和語法分析工具Bison開發，可以運行在Linux/UNIX平臺上。DTC的邏輯結構分為五個主要部分，分別為詞法分析、語法分析、靜態相關分析、AD函數庫和微分代碼轉換，如圖1所示。為了幫助用戶更有效地使用DTC系統生成微分代碼，下面介紹DTC的運行機制。

1)詞法分析 DTC使用詞法分析器生成器Flex[6]進行詞法分析，Flex是詞法分析程序Lex的實現。DTC接受用戶給定的原模式程序(將要計算微分的程序)，首先進行詞法分析，識別出C程序的單詞符號串。

2)語法分析 DTC使用語法分析器生成器Bison[7]進行語法分析，Bison是GNU的語法分析工具Yacc的升級。DTC的詞法分析得到的單詞符號串接下來進行語法分析，按照相應的語法規則進行規定的語義動作。

3)數據相關分析 DTC在詞法分析和語法分析的基礎上進行I/O相關分析和數據類型分析，作為微分代碼轉換、生成切線性模式程序和測試程序的必要步驟。為了處理C程序中各種變量的類型，特別是復雜類型如指針、數組、結構體、聯合體、typedef定義類型等，DTC使用富有特色的變量表簡潔有效地進行全局變量、函數參數及局部變量的記錄和查找，對于變量的復雜嵌套類型亦可準確處理。

4)AD函數庫 DTC的AD函數庫包含了奇異函數微分的等價數學定義及實現、一些有用接口及通用測試代碼等。對于計算程序中出現的奇異函數，如fabs、sqrt等，DTC的處理是調用其數學庫函數（在dtcmath.h中定義），如sign、diff_sqrt等，進行等價的微分代碼轉換。對奇異函數的處理在其他文獻中也稱做開關問題。

5）微分代碼轉換 在語法分析的語義動作中給出切線性模式微分代碼轉換方法，將原模式程序通過DTC的詞法分析、語法分析和相應微分轉換后，即可得到切線性模式程序；同時，在數據相關分析的基礎上，按照1.2.2節的方法即可生成切線性模式代碼可靠性的自動檢測程序。

2.2 詞法分析與語法分析

DTC使用詞法分析器生成器Flex和語法分析器生成器Bison進行詞法分析與語法分析，其工作流程如圖2所示。DTC系統給出Flex源程序和Bison源程序。其中，Flex源程序中給出描述C詞法結構的正則表達式和相應詞法動作，通過Flex工具自動轉換為詞法分析程序yylex()；Bison源程序給出C語法對應的語法規則和相應語義動作，通過Bison工具自動轉換為語法分析程序yyparse()。這里，語義動作實際上給出了恰當的微分代碼轉換方法。原模式輸入串通過詞法分析程序yylex()產生能被語法分析程序yyparse()識別的單詞符號串(終結符)，最后生成切線性模式輸出串。

2.3 變量表

DTC使用變量表詳細記錄程序中定義的全局變量、函數、函數參數、函數內的局部變量、用戶定義結構等。DTC的變量表定義了函數表、變量散列表、結構表、類型表等，如圖3所示。其中，函數表記錄了程序中定義的函數屬性。變量散列表分為全局變量散列表和函數局部變量散列表，分別記錄全局變量和函數局部變量，其變量插入散列表的散列函數選擇為除余函數。結構表記錄程序中定義的結構體（struct）、聯合體（union）和枚舉類型（enum），表中每個成員指向類型項，它記錄了C的復雜類型，即指針型、數組型、struct/union型、枚舉型和typedef型。其中，類型項的union域表示變量類型可以為這些復雜類型的其中一種。這里，類型項的每種復雜類型均被合理定義，如struct/union型中更詳細地記錄了結構體/聯合體的成員個數以及成員表，其他復雜類型類似。這樣，復雜類型可以嵌套定義，如結構體成員表中的成員也可以為結構體等復雜類型。最后，類型表記錄結構表定義以外的類型，包括基本類型以及指針、數組和用戶定義類型。

變量表的相應操作包括建表與初始化表、在表中插入元素、在表中查找元素、釋放表空間等。

DTC的變量表及相應表操作可以靈活而有效地處理各種C語言的函數定義與變量類型，包括處理C特有的復雜類型，如指針、數組、結構體/共用體類型、枚舉類型、typedef定義類型，以及嵌套定義的復雜類型。

2.4 微分代碼轉換

為了高效地進行切線性模式微分代碼轉換，DTC采取了一些簡潔的方法。例如，由于計算函數中恒為整數的變量其局部偏導數為0，DTC將程序中出現的變量分為整型(可微類型)和浮點型(不可微類型)。考慮C語法規則中的一個賦值z=f(x，y)(其中f(x，y)為一具體賦值函數)的切線性模式，可以有幾種情況：如果z為整型，則不對此賦值進行切線性模式變換(而不影響計算結果)；如果z為浮點型，則其切線性模式的變換可分別討論x與y的可微屬性來確定，如f(x，y)為基本運算(+、-、*、/)時，DTC所作的微分代碼轉換如表1所示。表中，前面有“diff_”前綴的變量表示微分擾動變量。f(x， y)為其余函數時亦是如此。DTC的這種微分代碼轉換方法可以避免所生成的微分模式程序中不必要的0項，從而帶來計算效率的提高。

表1 DTC中對基本運算的微分代碼轉換

原模式

切線性模式

x(整型)y(整型)x(浮點型)y(整型)x(整型)y(浮點型)x(浮點型)y(浮點型)

x+y0.0diff_xdiff_ydiff_x+diff_y

x-y0.0diff_xdiff_ydiff_x-diff_y

x*y0.0diff_x*yx*diff_yx*diff_y+y*diff_x

x/y(y!=0)0.0diff_x/y-x*diff_y/(y*y)(y*diff_x-x*diff_y)/(y*y)

2.5 I/O相關分析

下面給出I/O屬性的定義和I/O相關分析算法。

定義 I/O屬性。指定的程序片斷中，一個變量有三種I/O屬性，分別為輸入屬性(IN)、輸出屬性(OUT)、輸入/輸出屬性(IN_OUT)。其中，變量為輸入(IN)變量，指變量被引用一次或多次并且未改變其值；變量為輸出(OUT)變量，指變量被引用之前改變了其值；變量為輸入/輸出(IN_OUT)變量，指變量被引用一次或多次之后改變了其值。

一個變量可以在不同的程序片斷中具有不同的I/O屬性。特別地，DTC中考察函數參數在函數中的I/O屬性。I/O相關分析是在程序片斷內通過一個變量已知的上文I/O屬性以及當前賦值表達式(assignment expression)此變量的I/O屬性，判斷程序運行到此變量的I/O屬性，直到最終得到此變量在程序片斷內的IO屬性。為了算法實現上的方便，增加未知I/O屬性(UNKNOWN)。

算法 (I/O相關分析)。已知變量的上文I/O屬性及其在當前賦值表達式中的I/O屬性，通過表2的基本規則可以得到程序運行到此變量的I/O屬性。注意到如果變量當前I/O屬性為OUT或IN_OUT，則它將一直保持這個I/O屬性。當判斷程序中定義的函數其參數在函數中的I/O屬性時，可通過表2的規則按照每個賦值表達式依次進行計算，而得到參數在函數中的最終I/O屬性。特別地，當賦值表達式含有函數引用時，可遞歸地進行參數I/O屬性的計算。

表2 I/O屬性計算的基本規則

I/O屬性

當前賦值表達式中變量I/O屬性

UNKNOWNINOUTIN_OUT

變量上文I/O屬性UNKNOWNUNKNOWNINOUT

IN_OUT

ININININ_OUTIN_OUT

OUTOUTOUTOUTOUT

IN_OUTIN_OUTIN_OUTIN_OUTIN_OUT

3 測試與應用

DTC系統在Linux/UNIX平臺上進行了一些測試與應用，如表3所示。數值計算測試庫與斜壓一層半線性化準地轉位勢渦度方程模式是在曙光天闊服務器S4800A1上測試的，有關淺水波方程的應用是在IBM X40上測試的。其中，數值計算測試庫選取了實矩陣相乘、復矩陣相乘、實矩陣求逆、復矩陣求逆、Trench方法求Toeplitz矩陣逆和對稱正定矩陣的Cholesky分解六個測試（表3中前面6項）；準地轉位勢渦度方程模式選取了三個函數測試，計算網格為513×75（表3中第7~9項）；DTC系統的一項應用是求解二維全球正壓大氣淺水波方程，它運用科學計算工具箱PETSc[8]，采取“matrixfree”方式的NewtonKrylov方法求解，而DTC用于Jacobi矩陣—向量乘積的計算，取計算網格為640×320（表3中最后6項）。測試中，切線性模式的可靠性指標均達到了六位有效數字以上的精度，代碼復雜性均按照非注釋語句行數統計。

表3 DTC系統的數值測試結果

測試實例運行時間(2~10 s)

原模式切線性模式

代碼復雜性(行)

原模式切線性模式

從表中可以看到，DTC系統生成的切線性模式運行時間代價大致是原模式運行時間代價的兩倍左右。淺水波方程應用中ITA的計算由于切線性模式代碼經過了部分手工優化，其運行時間代價較小；而LF的切線性模式計算由于緩存效應，其運行時間代價較大。代碼復雜性方面切線性模式亦是原模式的兩倍左右。說明DTC系統在實際測試與應用中是高效和實用的。

4 結束語

自動微分是應用鏈式法則計算函數導數的一系列方法。對比傳統的手寫代碼、符號微分和有限差分方法，自動微分方法具有不存在截斷誤差、應用靈活且開發代價較小等特點。DTC是針對C語言、采取代碼轉換策略實現的自動微分工具。與另一個此類軟件ADIC不同，DTC采取切線性模式的計算模型，它在牛頓法求解非線性方程的Jacobi矩陣—向量乘積計算中有著重要應用。DTC運用編譯技術有效地處理C語言特有的指針與結構等，并且具有結構化微分代碼轉換、自動生成微分代碼檢測程序和靜態相關分析等特色。DTC系統在數值計算測試庫、斜壓一層半線性化準地轉位勢渦度方程模式和用PETSc實現的淺水波方程應用中的測試結果表明，其生成的切線性模式代碼可靠、高效而且實用。在下一步研究工作中，將進一步完善I/O相關分析及其他相關分析，將DTC應用到更多實際數值計算程序中，并對DTC系統作C++語言的部分擴展。

參考文獻：

［1］

GRIEWANK A.Evaluating derivatives:principles and techniques of algorithmic differentiation [M]. Philadelphia: SIAM， 2000.

［2］GRIEWANK A，JUEDES D，UTKE J.ADOLC:a package for the automatic differentiation of algorithms written in C/C++[J].ACM Trans on Mathematical Software，1996，22(2):131167.

［3］BENDTSEN C，STAUNING O.FADBAD:a flexible C++ package for automatic differentiation[R]. Lyngby: Department of Mathematical Modelling， Technical University of Denmark， 1996.

［4］BISCHOF C H， CARLE A， CORLISS G F，et al.ADIFOR: generating derivative codes from Fortran programs[J].Scientific Programming，1992，1(1): 1129.

［5］BISCHOF C H， ROH L， MAUER A.ADIC:an extensible automatic differentiation tool for ANSIC[J].Software Practice and Experience，1997，27(12):14271456.

［6］PAXSON V， ESTES W， MILLAWAY J.Flex: the fast lexical analyzer[EB/OL].(20080226)[20080328].http://www.gnu.org/software/flex/.

［7］DONNELLY C， STALLMAN R.Bison: GNU parser generator [EB/OL].(20070413)[20080328]. http://www.gnu.org/software/bison/.

［8］BALAY S， BUSCHELMAN K， EIJKHOUT V，et al.PETSc users manual[R].[S.l.]:Argonne National Laboratory， 2004.