Расчет Якобиана с использованием Stan Math дает 0 - PullRequest
Новая
       27

Расчет Якобиана с использованием Stan Math дает 0

0 голосов
/ 29 июня 2019

Я пытаюсь вычислить якобиан векторной функции. Модель, с которой я работаю, является версией Lotka-Volterra-Model (Hunter-Prey-Modell).

Я реализовал функцию оценки как функтор. 

Matrix<var, Dynamic, 1> LotkaVolterra::operator() (Matrix<var, Dynamic, 1> v) {
    Matrix<var, Dynamic, 1> result(v.size()); //prepare result vector
    double sum;
    for (int i = 0; i < v.size(); i++) { //calculations
        sum = bi.at(i);
        for (int j = 0; j < numOfEq; j++) {
            sum += v(j).val() * aij.at(i*numOfEq + j);
        }
        result(i) = sum*v(i).val();
    }
    return result;
}

В общем, я держался довольно близко к примеру, приведенному здесь (стр. 13). Так что мой фактический расчет якобиана выглядит примерно так же: 

void LotkaVolterra::jacobi(LotkaVolterra lv, int numEq, std::vector<double> v) {
        if (v.size()!= numEq) {std::cout << "ERROR";}

        //convert input to variable type
        cout << "x_var: " << endl;

        Matrix<var, Dynamic, 1> x_var(v.size()); //inputs as variables
        for (int i = 0; i < numEq; i++) {
            x_var(i) = v.at(i);
        }
        cout << x_var << endl;

        cout << "f_x_var:" << endl;

        //compute f(x) from variable type as variable  type
        Matrix<var, Dynamic, 1> f_x_var = lv(x_var);

        cout << f_x_var << endl;

        //convert f(x) from variable type back to double
        cout << "f_x: "<< endl;

        Matrix<double, Dynamic, 1> f_x(f_x_var.size());
        for (int i = 0; i< f_x_var.size(); i++) {
            f_x(i) = f_x_var(i).val();
        }
        cout << f_x << endl;
        Matrix<double, Dynamic, Dynamic>J(f_x_var.size(), x_var.size());
        for (int i = 0; i < f_x_var.size(); ++i) {
            if (i > 0) stan::math::set_zero_all_adjoints();
            f_x_var(i).grad();
            for (int j = 0; j < x_var.size(); ++j) {
                J(i,j) = x_var(j).adj();            
            }
        }
        std::cout << "J" << endl << J << endl;
}

К сожалению, якобиан J состоит только из нулей. Вывод на консоль выглядит следующим образом: 

x_var: 
0.6
0.4
0.5
0.5
0.3
f_x_var:
 0.5919
 0.3948
0.49261
0.49367
0.295758
f_x: 
 0.5919
 0.3948
0.49261
0.49367
0.295758
J
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0

Может кто-нибудь указать, что я делаю не так?

...