Vsebina

Ta modul uporablja:

https://github.com/leskovec/pyC_part.II

Struktura je zastavljena takole:

napišemo majhno C++ »znanstveno jedro«, ki uporablja zunanje knjižnice (recimo Eigen + GSL)
jedro prevedemo kot običajno C++ knjižnico (da ga lahko kliče tudi C++ program)
isto knjižnico ovijemo s pybind11, da jo lahko kliče tudi Python
z majhno Python skripto uvozimo modul in pogledamo vrnjene podatke

Ključna učna poanta ni samo »Eigen in GSL delujeta«, ampak:

Ko vežete knjižnično kodo, navadno želite vezati svoj majhen API, ne cele knjižnice. Tukaj je naš API ena funkcija get_my_eigvals(n), ki vrne majhen result strukt.

1) C++ API: glava, ki določa kaj izpostavimo

Glava knjižnice lib/lib.h definira:

majhen strukturiran tip result, ki hrani izhode v STL vsebnikih
funkcijo, ki jo želimo klicati iz C++ in Pythona: get_my_eigvals(n)

#pragma once

#include <iostream>

// eigen: matrix and eigenvalues
#include <Eigen/Dense>
#include <Eigen/Eigenvalues>

// gsl: random number generator
#include <GSL/gsl_rng.h>

// struct to store the result
struct result {
    std::vector<double> eigval;
    std::vector<std::vector<double>> matrix;
};

// wrapper function that does all the work
result get_my_eigvals(const int n);

Zakaj `result` uporablja `std::vector`

Python ne ve vedno kako delati s posamičnimi vsebniki; pozna Eigen ampak ne GSL recimo. Zna pa STL vsebnike avtomatsko pretvoriti v Python vsebnike.

Zato hranimo:

eigval kot std::vector<double> -> v Pythonu postane list[float]
matrix kot std::vector<std::vector<double>> -> v Pythonu postane »seznam seznamov«

To je zelo pogosta strategija, prijazna začetnikom:

numeriko izvajate v Eigen/GSL, rezultat pa na meji serializirate v STL vsebnike.

2) Numerično jedro: Eigen + GSL v `lib/lib.cxx`

Spodaj je implementacija. Tu se dejansko zgodi »znanstveni« del.

#include "lib.h"


result get_my_eigvals(const int n){
    Eigen::MatrixXd A = Eigen::MatrixXd::Zero(n,n);

    // use GSL to generate random numbers
    gsl_rng * r = gsl_rng_alloc (gsl_rng_taus);
    gsl_rng_set(r, 1);

    // populate A with random numbers from GSL
    for (int i=0; i<n; i++){
        for (int j=0; j<n; j++){
            A(i,j) = gsl_rng_uniform (r);
        }
    }

    // make a self-adjoint matrix
    Eigen::MatrixXd B = A.transpose()*A;

    // solve the eigensystem
    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigsys(B);

    // store the eigenvalues and the matrix to result and return
    // note, python won't know what to do with Eigen Vectors and Matrices,
    // but it can handle stl containers like std::vector

    result res;
    res.eigval.resize(n);
    for(int i=0; i<n; ++i){
      res.eigval[i] = eigsys.eigenvalues()(i);
    }
    res.matrix.resize(n);
    for(int i=0; i<n; ++i){
      res.matrix[i].resize(n);
      for(int j=0; j<n; ++j){
        res.matrix[i][j] = B(i,j);
      }
    }

    return res;
}

Kaj koda počne po vrstnem redu (informativno):

Korak A - sestavimo matriko `A` velikosti `n×n`

Eigen::MatrixXd A = Eigen::MatrixXd::Zero(n,n);

Začnemo z ničlami, nato A napolnimo z naključnimi vrednostmi.

Korak B - generiranje naključnih števil z GSL

gsl_rng * r = gsl_rng_alloc (gsl_rng_taus);
gsl_rng_set(r, 1);
A(i,j) = gsl_rng_uniform (r);

gsl_rng_* ustvari generator naključnih števil
gsl_rng_uniform(r) vrne vrednosti v intervalu [0, 1)
seme je fiksno (1), zato je primer ponovljiv (didaktično pomembno)

Korak C - zgradimo simetrično pozitivno semidefinitno matriko `B = Aᵀ A`

Eigen::MatrixXd B = A.transpose()*A;

To je trik s prosojnic: Aᵀ A je vedno simetrična in pozitivno semidefinitna matrika, zato je primerna za self-adjoint eigensolver.

Korak D - rešimo problem lastnih vrednosti

Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigsys(B);

Ker je B simetrična, uporabimo SelfAdjointEigenSolver.

Korak E - rezultate prekopiramo v STL vsebnike za Python

Koda nato prekopira:

lastne vrednosti v res.eigval
matriko B v res.matrix

To je ključen »most«. Da, to je kopiranje, vendar je zelo eksplicitno, lahko razumljivo in zanesljivo.

3) Gradnja knjižnice: `lib/CMakeLists.txt`

Kodo Eigen/GSL najprej obravnavamo kot knjižnico. To dosežemo z add_library(lib ...) + add_subdirectory(lib). Datoteka CMake za knjižnico:

target_sources(lib PUBLIC
    "${CMAKE_CURRENT_SOURCE_DIR}/lib.h"
    "${CMAKE_CURRENT_SOURCE_DIR}/lib.cxx"
)

find_package(Eigen3 3.4 REQUIRED NO_MODULE)
target_link_libraries(lib PUBLIC Eigen3::Eigen)

find_package(GSL REQUIRED)
target_link_libraries(lib PUBLIC GSL::gsl GSL::gslcblas)

target_include_directories(lib PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}" 
                                       "${Eigen_INCLUDE_DIRS}"  
                                       "${GSL_INCLUDE_DIRS}")
target_compile_definitions(lib PUBLIC BASE_VERSION=1.0)

Kaj to naredi:

lib.h + lib.cxx prevede v statični cilj knjižnice lib
poišče Eigen in GSL
ju poveže v lib
izvozi include poti, da lahko vsak cilj, ki poveže lib, vključuje "lib.h"

To je struktura »kot v pravem projektu«: jedro ni ujeto samo znotraj Python vezave.

4) Čisti C++ odjemalec: `main.cxx`

Preden karkoli vežemo, repozitorij pokaže še, kako isto funkcijo uporabimo iz C++ programa.

#include "lib.h"
#include <iostream>

int main(int argc, char* argv[]){

    if (argc < 2 ){
      std::cout << "Usage: " << argv[0] << " n" << std::endl;
      std::cout << "n: number of eigenvalues to compute" << std::endl;
      return 1;
    }

    int n = atoi(argv[1]);
    std::cout<< "computing eigenvalues of a " 
             << n << "x" << n << " matrix" 
             << std::endl;
    result r = get_my_eigvals(n);
    std::cout << "eigenvalues: ";
    for(int i=0; i<n; ++i){
      std::cout<< r.eigval[i] << " ";
    } 
    std::cout << std::endl;
    return 0;
}

To je dobra praksa v znanstvenem/inženirskem računanju:

jedro preverite brez Pythona
obdržite referenčni izvršljiv program, za katerega veste, da deluje
šele nato dodajte vezavo

5) pybind11 modul: `eigvals.cxx`

Poglejmo zdaj: »eigvals.cxx« modul.

#include <pybind11/pybind11.h>
#include <pybind11/stl.h>

#include <sstream>

#include "lib.h"

namespace py = pybind11;

PYBIND11_MODULE(eigvals, m) {
  m.doc() = "get eigenvalues from a positive definite random matrix";
  // add the results struct
  py::class_<result>(m, "result")
    .def(py::init<>())
    .def_readwrite("eigval", &result::eigval)
    .def_readwrite("matrix", &result::matrix);

  // add the get my evals function
  m.def("get_my_eigvals", &get_my_eigvals, "get eigenvalues from a positive definite random matrix", py::arg("n"),py::return_value_policy::reference);
}

Razčlenimo vezavno plast.

Definicija modula

PYBIND11_MODULE(eigvals, m) { ... }

To proizvede Python modul, ki ga uvozite kot:

import eigvals

`result` strukt kot Python razred

py::class_<result>(m, "result")
  .def(py::init<>())
  .def_readwrite("eigval", &result::eigval)
  .def_readwrite("matrix", &result::matrix);

To je razredna različica m.def(...):

ustvari Pythonu viden razred result
izpostavi dva javna člana (eigval, matrix) kot Python atributa

Ker smo vključili pybind11/stl.h, se STL vsebniki avtomatsko pretvarjajo.

Izpostavitev funkcije

m.def("get_my_eigvals", &get_my_eigvals, ...);

To v Pythonu omogoči klic eigvals.get_my_eigvals(n). Ker get_my_eigvals vrne objekt result po vrednosti, pybind11 iz vrnjenega C++ objekta ustvari Pythonu viden primerek result. Vezava pri tem tudi eksplicitno določi:

py::return_value_policy::reference

Pedagoško je tu pomemben namen:

vrnemo strukturiran tip (ne samo skalarja)
Python dobi instanco vezanega razreda result
člani strukture postanejo Python vsebniki

(V produkcijski vezavi bi življenjski cikel obravnavali še strožje; za delavnico pa je to ok.)

6) Trik »brez ročnega kopiranja .so«: `source_py/test.py`

V prejšnjih modulih ste preveden .so kopirali v Python mapo. V tem modulu pa pokažemo boljšo prakso: modul namestite in ga dodajte v sys.path.

Repozitorijski test.py:

import sys

sys.path.insert(0, '/path/to/pyC_part.II/install/pymodule')

import eigvals as ev


n=4
a = ev.get_my_eigvals(n)

print("matrix:")
for i in range(n):
    print(a.matrix[i][:])

print("eigenvalues:")
print(a.eigval)

print(type(a))
print(type(a.matrix))
print(type(a.eigval))

Ključna vrstica je:

sys.path.insert(0, ".../install/pymodule")

Ta ustreza install pravilu v vrhnji CMakeLists.txt:

install(TARGETS eigvals DESTINATION ${CMAKE_INSTALL_PREFIX}/pymodule)

Potek dela je torej:

prevedi + namesti v install/ prefiks
Python uvaža iz install/pymodule/

V praksi je to precej bolj praktično, kot ročno kopiranje .so datotek.

Kaj morate odnesti iz tega modula

Ta modul je malce aljši, ker pokaže realen potek:

jedrno numerično kodo v knjižnici (lib)
dve fronti uporabe: C++ izvršljivo datoteko in Python modul
majhen, ekspliciten API (get_my_eigvals(n), ki vrne result)
STL vsebnike na meji (ker se čisto mapirajo v Python)

Ko znate projekt strukturirati tako, lahko ovijete bistveno zahtevnejše znanstvene knjižnice, ne da bi vezavna koda postala nepregledna.

Vsebina

1) C++ API: glava, ki določa kaj izpostavimo

Zakaj result uporablja std::vector

2) Numerično jedro: Eigen + GSL v lib/lib.cxx

Korak A - sestavimo matriko A velikosti n×n