LXF138:GNU Scientific Library

Материал из Linuxformat.

GNU Scientific Library

На­уч­ная биб­лио­те­ка GNU

Ва­си­лий Оло­ни­чев на­по­ми­на­ет, что су­ще­ст­ви­тель­ное «ком­пь­ю­тер» про­ис­хо­дит от гла­го­ла «вы­чис­лять», и пред­ла­га­ет под­хо­дя­щую биб­лио­те­ку для этих са­мых вы­чис­ле­ний.

Бу­ду­чи ас­пи­ран­том, я бе­реж­но брал с пол­ки и нес к диско­во­ду раз­ме­ром с хо­ро­шую тум­боч­ку па­кет дисков ем­ко­стью 7 МБ и ве­сом ки­ло­грам­мов эдак 5. На нем на­хо­ди­лась «На­уч­ная биб­лио­те­ка Фор­тра­на» – как мне тогда ка­за­лось, са­мое глав­ное, ра­ди че­го и сде­ла­ны ком­пь­ю­те­ры. С тех пор мно­гое из­менилось: кое-кто уже и не по­доз­ре­ва­ет, что сло­ва «ком­пь­ю­тер» и «каль­ку­ля­тор» – си­нонимы. Они ду­ма­ют, что каль­ку­ля­тор – что­бы счи­тать, а ком­пь­ю­тер – что­бы бро­дить по Ин­тернету, иг­рать в иг­руш­ки и иногда на­би­рать что-нибудь в тек­сто­вом про­цес­со­ре. Хо­тя, как ни стран­но, и се­го­дня встре­ча­ют­ся чу­да­ки, пы­таю­щие­ся ис­поль­зо­вать ком­пь­ю­тер по его пер­во­на­чаль­но­му на­зна­чению. Имен­но для них и пред­на­зна­ча­ет­ся биб­лио­те­ка GSL – GNU Scientific Library (http://www.gnu.org/software/gsl), прав­да, напи­са­на она уже на C.

GSL соз­да­на, под­дер­жи­ва­ет­ся и раз­ви­ва­ет­ся в основ­ном уче­ны­ми-фи­зи­ка­ми из На­цио­наль­ной ла­бо­ра­то­рии в Лос-Ала­мо­се (США) и в пер­вую оче­редь для сво­их соб­ствен­ных нужд. Но рас­про­стра­ня­ет­ся она по ли­цен­зии GPL и по­это­му доступ­на всем. Об­нов­ления про­ис­хо­дят при­бли­зи­тель­но раз в пол­го­да. До­ку­мен­та­ция к биб­лио­те­ке про­сто пре­крас­ная, но толь­ко на анг­лий­ском язы­ке. При­ве­дем из нее па­ру аб­за­цев, даю­щих об­щее пред­став­ление о GSL:

• «Биб­лио­те­ка спро­ек­ти­ро­ва­на так, что раз­лич­ные ал­го­рит­мы мо­гут быть под­клю­че­ны или из­менены во вре­мя ис­полнения про­грам­мы, без ее пе­ре­ком­пи­ля­ции.
• Ин­тер­фейс доста­точ­но про­стой, и его мож­но ис­поль­зо­вать в язы­ках бо­лее вы­со­ко­го уров­ня, та­ких как GNU Guile или Python.
• Биб­лио­те­ка мо­жет быть ис­поль­зо­ва­на в мно­го­по­точ­ных при­ло­жениях без вся­ких ог­раничений.
• Там, где это воз­мож­но, функ­ции осно­ва­ны на на­деж­ных об­ще­доступ­ных па­ке­тах Фор­тра­на, та­ких как FFTPACK и QUADPACK.
• Биб­лио­те­ка не име­ет за­ви­си­мо­стей от дру­гих па­ке­тов и по­это­му без про­блем ком­пи­ли­ру­ет­ся и уста­нав­ли­ва­ет­ся».
• В GSL реа­ли­зо­ва­ны прак­ти­че­ски все ал­го­рит­мы вы­чис­ли­тель­ной ма­те­ма­ти­ки: ра­бо­та с по­ли­но­ма­ми, линей­ная ал­геб­ра, чис­лен­ное ин­тег­ри­ро­вание и диф­фе­рен­ци­ро­вание, по­иск корней, ре­шение систем диф­фе­рен­ци­аль­ных уравнений, ин­тер­по­ля­ция, ап­прок­си­ма­ция, ста­ти­сти­ка, оп­ти­ми­за­ция и так да­лее.

Сплай­ны

За­чем, спра­ши­ва­ет­ся, в на­ши дни са­мо­му пи­сать рас­чет­ные про­грам­мы на C, ес­ли есть Scilab, Octave, R и про­чие? Во-пер­вых, фи­зи­ки не шу­тят, когда принима­ют­ся за вы­чис­ления, и по­это­му все­гда есть по­треб­ность в эф­фек­тив­ных про­грам­мах. Во-вто­рых, есть за­да­чи, где, по­ми­мо ин­тен­сив­ных рас­че­тов, тре­бу­ет­ся доступ к систем­ным ре­сур­сам – в ча­ст­но­сти, сред­ствам меж­про­цесс­но­го взаи­мо­дей­ствия. Та­кими за­да­ча­ми мо­гут быть при­ло­жения ре­аль­но­го вре­мени, сде­лан­ные в сти­ле «Unix way» и вы­пол­няю­щие­ся на про­мыш­лен­ных ком­пь­ю­те­рах. И, на­конец, есть пре­ве­ли­кое мно­же­ство тол­ко­вых ин­женеров и эко­но­ми­стов, пред­по­чи­таю­щих все де­лать свои­ми ру­ка­ми, и для них GSL под­хо­дит иде­аль­но. Да и за­да­чи встре­ча­ют­ся са­мые раз­но­об­раз­ные: для одних хо­ро­шо под­хо­дит Scilab, а для дру­гих луч­ше напи­сать про­грам­му.

Рас­смот­рим при­мер. Пусть управ­ляю­щая систе­ма снима­ет по­ка­зания с дат­чи­ка, под­клю­чен­но­го че­рез пре­об­ра­зо­ва­тель–уси­ли­тель–АЦП. Ста­ти­че­ская ха­рак­те­ри­сти­ка дат­чи­ка нелиней­на. В этом слу­чае удоб­но ис­поль­зо­вать ин­тер­по­ля­цию сплай­на­ми. GSL со­дер­жит боль­шой на­бор средств для этих це­лей, бу­к­валь­но на все слу­чаи жизни. В про­стей­шем слу­чае, как наш, мож­но ог­раничить­ся ис­поль­зо­ванием сле­дую­щих функ­ций:

• под­го­то­ви­тель­ные

gsl_interp_accel * gsl_interp_accel_alloc (void)
 gsl_spline * gsl_spline_alloc (const gsl interp type * T, size t size)
 int gsl_spline_init (gsl spline *spline, const double xa[], const double ya[], size t size)

• ос­нов­ная ра­бо­чая

double gsl_spline_eval (const gsl spline * spline, double x, gsl interp accel * acc)

• за­вер­шаю­щие

void gsl_interp_accel_free (gsl interp accel* acc)
 void gsl_spline_free (gsl spline * spline)

Пер­вая из них соз­да­ет в па­мя­ти объ­ект-ак­се­ле­ра­тор. Это сво­его ро­да ите­ра­тор, со­хра­няю­щий пре­ды­ду­щие об­ра­ще­ния и по­зво­ляю­щий ус­ко­рить по­сле­дую­щие вы­чис­ле­ния. Сле­дую­щая функ­ция соз­да­ет в па­мя­ти сам сплайн. Сплай­ны мо­гут быть са­мых раз­но­об­раз­ных ти­пов. Кста­ти, этот тот са­мый слу­чай, ко­гда ал­го­ритм мо­жет быть из­ме­нен во вре­мя ра­бо­ты про­грам­мы.

Дан­ные функ­ции са­ми про­сят­ся, что­бы их ин­кап­су­ли­ро­ва­ли в  класс, кон­ст­рук­тор ко­то­ро­го по­лу­ча­ет в ка­че­ст­ве па­ра­мет­ра спе­ци­фи­ка­цию тек­сто­во­го фай­ла, со­дер­жа­ще­го гра­дуи­ро­воч­ную кри­вую в ви­де «вход АЦП – зна­че­ние фи­зи­че­ской ве­ли­чи­ны, из­ме­ряе­мой дат­чи­ком». Он, а так­же про­стень­кая про­грам­ма, его ис­поль­зую­щая, име­ют­ся на LXFDVD. А вот фраг­мент ко­да это­го клас­са:

class Sensor{
 protected:
     int N;
     double *x;
     double *y;
     gsl_interp_accel *acc;
     gsl_spline *spline;
 public:
     Sensor(char* f_calibr){
   // Чтение данных из текстового файла с гра дуировочной кривой
        ...
        x=new double[N*sizeof(double)];
        y=new double[N*sizeof(double)];
        ...
        acc = gsl_interp_accel_alloc();
        spline = gsl_spline_alloc(gsl_interp_cspline, N);
        gsl_spline_init(spline, x, y, N);
     }
     ~Sensor(){
        gsl_spline_free(spline);
        gsl_interp_accel_free(acc);
        delete[] x;
        delete[] y;
     }
     double inline getVal(double xi){
        return gsl_spline_eval(spline, xi, acc);
     }
 };

Здесь ис­поль­зу­ет­ся сплайн ти­па gsl_interp_cspline, со­от­вет­ст­вую­щий «клас­си­че­ско­му» ку­би­че­ско­му сплай­ну, у ко­то­ро­го вто­рая про­из­вод­ная в пер­вой и по­след­ней точ­ке ин­тер­ва­ла рав­на 0.

В про­ек­те не­об­хо­ди­мо ука­зать на биб­лио­те­ки GSL, ко­то­рые в мо­ей сис­те­ме раз­ме­ща­ют­ся в ка­та­ло­ге /usr/local/lib. Ес­ли вы ра­бо­тае­те в Code::Blocks, то на вклад­ке Linker settings сле­ду­ет пря­мо ука­зать спе­ци­фи­ка­цию биб­лио­те­ки, на­при­мер /usr/local/lib/lbgsl.so. Ес­ли вы ис­поль­зуе­те Qmake, то в файл про­ек­та сле­ду­ет до­ба­вить сле­дую­щую стро­ку: LIBS += -L/usr/local/lib -lgsl. Кро­ме этой ос­нов­ной биб­лио­те­ки, ес­ли вы пря­мо или опо­сре­до­ван­но ис­поль­зуе­те мат­ри­цы, по­тре­бу­ет­ся до­ба­вить биб­лио­те­ку libgslcblas.

Ме­тод наи­мень­ших квад­ра­тов

При­ве­дем еще один при­мер из об­лас­ти про­мав­то­ма­ти­ки, в ко­то­ром по­ка­жем при­ме­не­ние GSL для иден­ти­фи­ка­ции объ­ек­та управ­ле­ния ме­то­дом наи­мень­ших квад­ра­тов. Пол­ный текст про­грам­мы так­же име­ет­ся на дис­ке, при­ла­гае­мом к жур­на­лу, а здесь пред­став­лен ее фраг­мент.

int M; // по­ря­док объ­ек­та
   int N; // кол-во строк в ис­поль­зуе­мых век­то­рах и мат­ри­це
   // вспо­мо­га­тель­ные мас­си­вы для раз­ме­ще­ния
   // счи­тан­ных из фай­ла экс­пе­ри­мен­таль­ных дан­ных
   std::vector<double> yV,uV,tV;
   gsl_matrix *F;
   gsl_vector *T,*R,*S, *Y;
 
   // считывание данных эксперимента из текстового файла
   ...
   N=yV.size()­M;
   // выде ление памя ти под век торы и мат рицы
   F=gsl_matrix_calloc(N,2*M);
   T = gsl_vector_calloc(2*M);
   R = gsl_vector_calloc(2*M);
   S = gsl_vector_calloc(N);
   Y = gsl_vector_calloc(N);
 
   // заполнение вектора Y и матрицы F данными эксперимента
   for(int i=0; i<N; i++) gsl_vector_set(Y,i,yV[i+M]);
   for(int i=0; i<N; i++){
       for(int j=0;j<M;j++){
          gsl_matrix_set(F,i,j,(­1.0)*yV[i+M­j­1]);
       }
   }
   for(int i=0; i<N;i++){
       for(int j=0;j<M;j++){
          gsl_matrix_set(F,i,j+M,uV[i+M­j­1]);
       }
   }
 
   // определение параметров объекта
   gsl_linalg_QR_decomp(F,T);
   gsl_linalg_QR_lssolve(F,T,Y,R,S);
 
   // вывод результатов
   for(int j=0; j<M+M; j++)
       printf(“%f\n”,gsl_vector_get(R,j));
 
   // освобождение динамической памяти
   gsl_matrix_free(F);
   gsl_vector_free(T);
   gsl_vector_free(R);
   gsl_vector_free(S);
   gsl_vector_free(Y);

Здесь функ­ции gsl_matrix_calloc() и gsl_vector_calloc() вы­де­ля­ют па­мять под мат­ри­цы и век­то­ры за­дан­ных раз­ме­ров, а gsl_matrix_free() и gsl_vector_free(), на­обо­рот, ос­во­бо­ж­да­ют па­мять. Функ­ции gsl_matrix_set() и gsl_vector_set() по­зво­ля­ют за­да­вать зна­че­ния эле­мен­тов мат­риц и век­то­ров. Не­труд­но до­га­-дать­ся, что су­ще­ст­ву­ют сим­мет­рич­ные функ­ции gsl_matrix_get() и gsl_vector_get(). Оп­ре­де­ле­ние па­ра­мет­ров мо­де­ли ме­то­дом наи­мень­ших квад­ра­тов осу­ще­ст­в­ля­ет­ся при по­мо­щи функ­ций gsl_linalg_QR_decomp() и gsl_linalg_QR_lssolve(). Пер­вая вы­пол­ня­ет пред­ва­ри­тель­ную де­ком­по­зи­цию мат­ри­цы F, а вто­рая – ре­ша­ет за­да­чу. Век­то­ры T и S – вспо­мо­га­тель­ные. Век­тор Y со­дер­жит зна­че­ния на вы­хо­де объ­ек­та, а мат­ри­ца F – зна­че­ния на вхо­де объ­ек­та и зна­че­ния на вы­хо­де объ­ек­та в пре­ды­ду­щие мо­мен­ты вре­ме­ни. Ре­ше­ние, со­дер­жа­щее ко­эф­фи­ци­ен­ты чис­ли­те­ля и зна­ме­на­те­ля раз­но­ст­ной пе­ре­да­точ­ной функ­ции, по­ме­ща­ют­ся в век­тор R.

Лю­бой спе­циа­лист без тру­да под­бе­рет де­ся­ток-дру­гой за­дач из сво­ей пред­мет­ной об­лас­ти, тре­бую­щих вы­чис­ле­ний, для ко­то­рых про­грам­му про­сто так с хо­ду не на­пи­шешь. Да да­же ес­ли и на­пи­шешь, то она, ско­рее все­го, бу­дет не­эф­фек­тив­на и не­пред­ска­зуе­ма во не­штат­ных си­туа­ци­ях. Здесь, как и во мно­гих дру­гих слу­ча­ях, сво­бод­ный код по­мо­га­ет ре­шить все про­бле­мы.