讓我們使用 C99 和 C++11 完成常見的數(shù)據(jù)科學(xué)任務(wù)。

雖然 Python 和 R 之類的語(yǔ)言在數(shù)據(jù)科學(xué)中越來(lái)越受歡迎，但是 C 和 C++ 對(duì)于高效的數(shù)據(jù)科學(xué)來(lái)說(shuō)是一個(gè)不錯(cuò)的選擇。在本文中，我們將使用 C99 和 C++11 編寫一個(gè)程序，該程序使用 Anscombe 的四重奏數(shù)據(jù)集，下面將對(duì)其進(jìn)行解釋。

我在一篇涉及 Python 和 GNU Octave 的文章中寫了我不斷學(xué)習(xí)編程語(yǔ)言的動(dòng)機(jī)，值得大家回顧。這里所有的程序都需要在命令行上運(yùn)行，而不是在圖形用戶界面(GUI)上運(yùn)行。完整的示例可在 polyglot_fit 存儲(chǔ)庫(kù)中找到。

編程任務(wù)

你將在本系列中編寫的程序：

• 從 CSV 文件中讀取數(shù)據(jù)
• 用直線插值數(shù)據(jù)(即 f(x)=m ⋅ x + q)
• 將結(jié)果繪制到圖像文件

這是許多數(shù)據(jù)科學(xué)家遇到的普遍情況。示例數(shù)據(jù)是 Anscombe 的四重奏的第一組，如下表所示。這是一組人工構(gòu)建的數(shù)據(jù)，當(dāng)擬合直線時(shí)可以提供相同的結(jié)果，但是它們的曲線非常不同。數(shù)據(jù)文件是一個(gè)文本文件，其中的制表符用作列分隔符，前幾行作為標(biāo)題。該任務(wù)將僅使用第一組(即前兩列)。

C 語(yǔ)言的方式

C 語(yǔ)言是通用編程語(yǔ)言，是當(dāng)今使用最廣泛的語(yǔ)言之一(依據(jù) TIOBE 指數(shù)、RedMonk 編程語(yǔ)言排名、編程語(yǔ)言流行度指數(shù)和 GitHub Octoverse 狀態(tài) 得來(lái))。這是一種相當(dāng)古老的語(yǔ)言(大約誕生在 1973 年)，并且用它編寫了許多成功的程序(例如 Linux 內(nèi)核和 Git 僅是其中的兩個(gè)例子)。它也是最接近計(jì)算機(jī)內(nèi)部運(yùn)行機(jī)制的語(yǔ)言之一，因?yàn)樗苯佑糜诓僮鲀?nèi)存。它是一種編譯語(yǔ)言;因此，源代碼必須由編譯器轉(zhuǎn)換為機(jī)器代碼。它的標(biāo)準(zhǔn)庫(kù)很小，功能也不多，因此人們開發(fā)了其它庫(kù)來(lái)提供缺少的功能。

我最常在數(shù)字運(yùn)算中使用該語(yǔ)言，主要是因?yàn)槠湫阅?。我覺(jué)得使用起來(lái)很繁瑣，因?yàn)樗枰芏鄻影宕a，但是它在各種環(huán)境中都得到了很好的支持。C99 標(biāo)準(zhǔn)是最新版本，增加了一些漂亮的功能，并且得到了編譯器的良好支持。

我將一路介紹 C 和 C++ 編程的必要背景，以便初學(xué)者和高級(jí)用戶都可以繼續(xù)學(xué)習(xí)。

安裝

要使用 C99 進(jìn)行開發(fā)，你需要一個(gè)編譯器。我通常使用 Clang，不過(guò) GCC 是另一個(gè)有效的開源編譯器。對(duì)于線性擬合，我選擇使用 GNU 科學(xué)庫(kù)。對(duì)于繪圖，我找不到任何明智的庫(kù)，因此該程序依賴于外部程序：Gnuplot。該示例還使用動(dòng)態(tài)數(shù)據(jù)結(jié)構(gòu)來(lái)存儲(chǔ)數(shù)據(jù)，該結(jié)構(gòu)在伯克利軟件分發(fā)版(BSD)中定義。

在 Fedora 中安裝很容易：

sudo dnf install clang gnuplot gsl gsl-devel 

代碼注釋

在 C99 中，注釋的格式是在行的開頭放置 //，行的其它部分將被解釋器丟棄。另外，/* 和 */ 之間的任何內(nèi)容也將被丟棄。

// 這是一個(gè)注釋，會(huì)被解釋器忽略 
/* 這也被忽略 */ 

必要的庫(kù)

庫(kù)由兩部分組成：

• 頭文件，其中包含函數(shù)說(shuō)明
• 包含函數(shù)定義的源文件

頭文件包含在源文件中，而庫(kù)文件的源文件則鏈接到可執(zhí)行文件。因此，此示例所需的頭文件是：

// 輸入/輸出功能 
#include <stdio.h> 
// 標(biāo)準(zhǔn)庫(kù) 
#include <stdlib.h> 
// 字符串操作功能 
#include <string.h> 
// BSD 隊(duì)列 
#include <sys/queue.h> 
// GSL 科學(xué)功能 
#include <gsl/gsl_fit.h> 
#include <gsl/gsl_statistics_double.h> 

主函數(shù)

在 C 語(yǔ)言中，程序必須位于稱為主函數(shù) main() 的特殊函數(shù)內(nèi)：

int main(void) { 
    ... 
} 

這與上一教程中介紹的 Python 不同，后者將運(yùn)行在源文件中找到的所有代碼。

定義變量

在 C 語(yǔ)言中，變量必須在使用前聲明，并且必須與類型關(guān)聯(lián)。每當(dāng)你要使用變量時(shí)，都必須決定要在其中存儲(chǔ)哪種數(shù)據(jù)。你也可以指定是否打算將變量用作常量值，這不是必需的，但是編譯器可以從此信息中受益。 以下來(lái)自存儲(chǔ)庫(kù)中的 fitting_C99.c 程序：

const char *input_file_name = "anscombe.csv"; 
const char *delimiter = "\t"; 
const unsigned int skip_header = 3; 
const unsigned int column_x = 0; 
const unsigned int column_y = 1; 
const char *output_file_name = "fit_C99.csv"; 
const unsigned int N = 100; 

C 語(yǔ)言中的數(shù)組不是動(dòng)態(tài)的，從某種意義上說(shuō)，數(shù)組的長(zhǎng)度必須事先確定(即，在編譯之前)：

int data_array[1024]; 

由于你通常不知道文件中有多少個(gè)數(shù)據(jù)點(diǎn)，因此請(qǐng)使用單鏈列表。這是一個(gè)動(dòng)態(tài)數(shù)據(jù)結(jié)構(gòu)，可以無(wú)限增長(zhǎng)。幸運(yùn)的是，BSD 提供了鏈表。這是一個(gè)示例定義：

struct data_point { 
    double x; 
    double y; 
    SLIST_ENTRY(data_point) entries; 
}; 
SLIST_HEAD(data_list, data_point) head = SLIST_HEAD_INITIALIZER(head); 
SLIST_INIT(&head); 

該示例定義了一個(gè)由結(jié)構(gòu)化值組成的 data_point 列表，該結(jié)構(gòu)化值同時(shí)包含 x 值和 y 值。語(yǔ)法相當(dāng)復(fù)雜，但是很直觀，詳細(xì)描述它就會(huì)太冗長(zhǎng)了。

打印輸出

要在終端上打印，可以使用 printf() 函數(shù)，其功能類似于 Octave 的 printf() 函數(shù)(在第一篇文章中介紹)：

printf("#### Anscombe's first set with C99 ####\n"); 

printf() 函數(shù)不會(huì)在打印字符串的末尾自動(dòng)添加換行符，因此你必須添加換行符。第一個(gè)參數(shù)是一個(gè)字符串，可以包含傳遞給函數(shù)的其他參數(shù)的格式信息，例如：

printf("Slope: %f\n", slope); 

讀取數(shù)據(jù)

現(xiàn)在來(lái)到了困難的部分……有一些用 C 語(yǔ)言解析 CSV 文件的庫(kù)，但是似乎沒(méi)有一個(gè)庫(kù)足夠穩(wěn)定或流行到可以放入到 Fedora 軟件包存儲(chǔ)庫(kù)中。我沒(méi)有為本教程添加依賴項(xiàng)，而是決定自己編寫此部分。同樣，討論這些細(xì)節(jié)太啰嗦了，所以我只會(huì)解釋大致的思路。為了簡(jiǎn)潔起見，將忽略源代碼中的某些行，但是你可以在存儲(chǔ)庫(kù)中找到完整的示例代碼。

首先，打開輸入文件：

FILE* input_file = fopen(input_file_name, "r"); 

然后逐行讀取文件，直到出現(xiàn)錯(cuò)誤或文件結(jié)束：

while (!ferror(input_file) && !feof(input_file)) { 
    size_t buffer_size = 0; 
    char *buffer = NULL; 
    
    getline(&buffer, &buffer_size, input_file); 
    ... 
} 

getline() 函數(shù)是 POSIX.1-2008 標(biāo)準(zhǔn)新增的一個(gè)不錯(cuò)的函數(shù)。它可以讀取文件中的整行，并負(fù)責(zé)分配必要的內(nèi)存。然后使用 strtok() 函數(shù)將每一行分成字元token。遍歷字元，選擇所需的列：

char *token = strtok(buffer, delimiter); 
while (token != NULL) 
{ 
    double value; 
    sscanf(token, "%lf", &value); 
    if (column == column_x) { 
        x = value; 
    } else if (column == column_y) { 
        y = value; 
    } 
    column += 1; 
    token = strtok(NULL, delimiter); 
} 

最后，當(dāng)選擇了 x 和 y 值時(shí)，將新數(shù)據(jù)點(diǎn)插入鏈表中：

struct data_point *datum = malloc(sizeof(struct data_point)); 
datum->x = x; 
datum->y = y; 
SLIST_INSERT_HEAD(&head, datum, entries); 

malloc() 函數(shù)為新數(shù)據(jù)點(diǎn)動(dòng)態(tài)分配(保留)一些持久性內(nèi)存。

擬合數(shù)據(jù)

GSL 線性擬合函數(shù) gslfitlinear() 期望其輸入為簡(jiǎn)單數(shù)組。因此，由于你將不知道要?jiǎng)?chuàng)建的數(shù)組的大小，因此必須手動(dòng)分配它們的內(nèi)存：

const size_t entries_number = row - skip_header - 1; 
double *x = malloc(sizeof(double) * entries_number); 
double *y = malloc(sizeof(double) * entries_number); 

然后，遍歷鏈表以將相關(guān)數(shù)據(jù)保存到數(shù)組：

SLIST_FOREACH(datum, &head, entries) { 
    const double current_x = datum->x; 
    const double current_y = datum->y; 
    x[i] = current_x; 
    y[i] = current_y; 
    i += 1; 
} 

現(xiàn)在你已經(jīng)處理完了鏈表，請(qǐng)清理它。要總是釋放已手動(dòng)分配的內(nèi)存，以防止內(nèi)存泄漏。內(nèi)存泄漏是糟糕的、糟糕的、糟糕的(重要的話說(shuō)三遍)。每次內(nèi)存沒(méi)有釋放時(shí)，花園侏儒都會(huì)找不到自己的頭：

while (!SLIST_EMPTY(&head)) { 
    struct data_point *datum = SLIST_FIRST(&head); 
    SLIST_REMOVE_HEAD(&head, entries); 
    free(datum); 
} 

終于，終于!你可以擬合你的數(shù)據(jù)了：

gsl_fit_linear(x, 1, y, 1, entries_number, 
               &intercept, &slope, 
               &cov00, &cov01, &cov11, &chi_squared); 
const double r_value = gsl_stats_correlation(x, 1, y, 1, entries_number); 
printf("Slope: %f\n", slope); 
printf("Intercept: %f\n", intercept); 
printf("Correlation coefficient: %f\n", r_value); 

繪圖

你必須使用外部程序進(jìn)行繪圖。因此，將擬合數(shù)據(jù)保存到外部文件：

const double step_x = ((max_x + 1) - (min_x - 1)) / N; 
for (unsigned int i = 0; i < N; i += 1) { 
    const double current_x = (min_x - 1) + step_x * i; 
    const double current_y = intercept + slope * current_x; 
    fprintf(output_file, "%f\t%f\n", current_x, current_y); 
} 

用于繪制兩個(gè)文件的 Gnuplot 命令是：

plot 'fit_C99.csv' using 1:2 with lines title 'Fit', 'anscombe.csv' using 1:2 with points pointtype 7 title 'Data' 

結(jié)果

在運(yùn)行程序之前，你必須編譯它：

clang -std=c99 -I/usr/include/ fitting_C99.c -L/usr/lib/ -L/usr/lib64/ -lgsl -lgslcblas -o fitting_C99 

這個(gè)命令告訴編譯器使用 C99 標(biāo)準(zhǔn)、讀取 fitting_C99.c 文件、加載 gsl 和 gslcblas 庫(kù)、并將結(jié)果保存到 fitting_C99。命令行上的結(jié)果輸出為：

#### Anscombe's first set with C99 #### 
Slope: 0.500091 
Intercept: 3.000091 
Correlation coefficient: 0.816421 

這是用 Gnuplot 生成的結(jié)果圖像：

C++11 方式

C++ 語(yǔ)言是一種通用編程語(yǔ)言，也是當(dāng)今使用的最受歡迎的語(yǔ)言之一。它是作為 C 的繼承人創(chuàng)建的(誕生于 1983 年)，重點(diǎn)是面向?qū)ο蟪绦蛟O(shè)計(jì)(OOP)。C++ 通常被視為 C 的超集，因此 C 程序應(yīng)該能夠使用 C++ 編譯器進(jìn)行編譯。這并非完全正確，因?yàn)樵谀承O端情況下它們的行為有所不同。 根據(jù)我的經(jīng)驗(yàn)，C++ 與 C 相比需要更少的樣板代碼，但是如果要進(jìn)行面向?qū)ο箝_發(fā)，語(yǔ)法會(huì)更困難。C++11 標(biāo)準(zhǔn)是最新版本，增加了一些漂亮的功能，并且基本上得到了編譯器的支持。

由于 C++ 在很大程度上與 C 兼容，因此我將僅強(qiáng)調(diào)兩者之間的區(qū)別。我在本部分中沒(méi)有涵蓋的任何部分，則意味著它與 C 中的相同。

安裝

這個(gè) C++ 示例的依賴項(xiàng)與 C 示例相同。 在 Fedora 上，運(yùn)行：

sudo dnf install clang gnuplot gsl gsl-devel 

必要的庫(kù)

庫(kù)的工作方式與 C 語(yǔ)言相同，但是 include 指令略有不同：

#include <cstdlib> 
#include <cstring> 
#include <iostream> 
#include <fstream> 
#include <string> 
#include <vector> 
#include <algorithm> 
extern "C" { 
#include <gsl/gsl_fit.h> 
#include <gsl/gsl_statistics_double.h> 
} 

由于 GSL 庫(kù)是用 C 編寫的，因此你必須將這個(gè)特殊情況告知編譯器。

定義變量

與 C 語(yǔ)言相比，C++ 支持更多的數(shù)據(jù)類型(類)，例如，與其 C 語(yǔ)言版本相比，string 類型具有更多的功能。相應(yīng)地更新變量的定義：

const std::string input_file_name("anscombe.csv"); 

對(duì)于字符串之類的結(jié)構(gòu)化對(duì)象，你可以定義變量而無(wú)需使用 = 符號(hào)。

打印輸出

你可以使用 printf() 函數(shù)，但是 cout 對(duì)象更慣用。使用運(yùn)算符 << 來(lái)指示要使用 cout 打印的字符串(或?qū)ο?：

std::cout << "#### Anscombe's first set with C++11 ####" << std::endl; 
... 
std::cout << "Slope: " << slope << std::endl; 
std::cout << "Intercept: " << intercept << std::endl; 
std::cout << "Correlation coefficient: " << r_value << std::endl; 

讀取數(shù)據(jù)

該方案與以前相同。將打開文件并逐行讀取文件，但語(yǔ)法不同：

std::ifstream input_file(input_file_name); 
while (input_file.good()) { 
    std::string line; 
    getline(input_file, line); 
    ... 
} 

使用與 C99 示例相同的功能提取行字元。代替使用標(biāo)準(zhǔn)的 C 數(shù)組，而是使用兩個(gè)向量。向量是 C++ 標(biāo)準(zhǔn)庫(kù)中對(duì) C 數(shù)組的擴(kuò)展，它允許動(dòng)態(tài)管理內(nèi)存而無(wú)需顯式調(diào)用 malloc()：

std::vector<double> x; 
std::vector<double> y; 
// Adding an element to x and y: 
x.emplace_back(value); 
y.emplace_back(value); 

擬合數(shù)據(jù)

要在 C++ 中擬合，你不必遍歷列表，因?yàn)橄蛄靠梢员ＷC具有連續(xù)的內(nèi)存。你可以將向量緩沖區(qū)的指針直接傳遞給擬合函數(shù)：

gsl_fit_linear(x.data(), 1, y.data(), 1, entries_number, 
               &intercept, &slope, 
               &cov00, &cov01, &cov11, &chi_squared); 
const double r_value = gsl_stats_correlation(x.data(), 1, y.data(), 1, entries_number); 
std::cout << "Slope: " << slope << std::endl; 
std::cout << "Intercept: " << intercept << std::endl; 
std::cout << "Correlation coefficient: " << r_value << std::endl; 

繪圖

使用與以前相同的方法進(jìn)行繪圖。 寫入文件：

const double step_x = ((max_x + 1) - (min_x - 1)) / N; 
for (unsigned int i = 0; i < N; i += 1) { 
    const double current_x = (min_x - 1) + step_x * i; 
    const double current_y = intercept + slope * current_x; 
    output_file << current_x << "\t" << current_y << std::endl; 
} 
output_file.close(); 

然后使用 Gnuplot 進(jìn)行繪圖。

結(jié)果

在運(yùn)行程序之前，必須使用類似的命令對(duì)其進(jìn)行編譯：

clang++ -std=c++11 -I/usr/include/ fitting_Cpp11.cpp -L/usr/lib/ -L/usr/lib64/ -lgsl -lgslcblas -o fitting_Cpp11 

命令行上的結(jié)果輸出為：

#### Anscombe's first set with C++11 #### 
Slope: 0.500091 
Intercept: 3.00009 
Correlation coefficient: 0.816421 

這就是用 Gnuplot 生成的結(jié)果圖像：

結(jié)論

本文提供了用 C99 和 C++11 編寫的數(shù)據(jù)擬合和繪圖任務(wù)的示例。由于 C++ 在很大程度上與 C 兼容，因此本文利用了它們的相似性來(lái)編寫了第二個(gè)示例。在某些方面，C++ 更易于使用，因?yàn)樗糠譁p輕了顯式管理內(nèi)存的負(fù)擔(dān)。但是其語(yǔ)法更加復(fù)雜，因?yàn)樗肓藶? OOP 編寫類的可能性。但是，仍然可以用 C 使用 OOP 方法編寫軟件。由于 OOP 是一種編程風(fēng)格，因此可以在任何語(yǔ)言中使用。在 C 中有一些很好的 OOP 示例，例如 GObject 和 Jansson庫(kù)。

對(duì)于數(shù)字運(yùn)算，我更喜歡在 C99 中進(jìn)行，因?yàn)樗恼Z(yǔ)法更簡(jiǎn)單并且得到了廣泛的支持。直到最近，C++11 還沒(méi)有得到廣泛的支持，我傾向于避免使用先前版本中的粗糙不足之處。對(duì)于更復(fù)雜的軟件，C++ 可能是一個(gè)不錯(cuò)的選擇。

你是否也將 C 或 C++ 用于數(shù)據(jù)科學(xué)?在評(píng)論中分享你的經(jīng)驗(yàn)。