// TestGenericMatrixOde.cpp
//
// Computing exp(At) as:
//
//	dB/dt = A*B for t > 0
//  B(0) = C (identity matrix), where C is a given matrix
//
// (C) Datasim Education BV 2018
//

#ifdef _MSC_VER
#define _CRT_SECURE_NO_WARNINGS
#endif

#include <iostream>
#include <map>
#include <vector>
#include <string>
#include <fstream>

#include <boost/numeric/odeint.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>

namespace ublas = boost::numeric::ublas;

using value_type = double;
using state_type = boost::numeric::ublas::matrix<value_type>;

class MatrixOde
{
private:
	// dB/dt = A*B, B(0) = C;
	ublas::matrix<value_type> A_;
	ublas::matrix<value_type> C_;
public:

	MatrixOde(const ublas::matrix<value_type>& A, const ublas::matrix<value_type>& IC)
		: A_(A), C_(IC) {}

    void operator()( const state_type &x , state_type &dxdt , double t ) const
    {

        for( std::size_t i=0 ; i < x.size1();++i )
        {
            for( std::size_t j=0 ; j < x.size2(); ++j )
            {
				dxdt(i, j) = 0.0;
				for (std::size_t k = 0; k < x.size2(); ++k)
				{
					dxdt(i, j) += A_(i,k)*x(k,j);
				}
			}
		}
 	}

};


void write_out(const state_type& U, const double t )
{

	{
		// put your wishes here
	}
}

int main()
{

	namespace Bode = boost::numeric::odeint;

 	std::size_t N = 2;
	ublas::identity_matrix<value_type> C(N);
	ublas::matrix<value_type> A(N,N);
	//A(0, 0) = A(1, 1) = 3.0; A(0, 1) = 1.0; A(1, 0) = 0.0;
	A(0, 0) = 6.0;  A(0, 1) = -1.0; A(1,0) = 9.0; A(1, 1) = 0.0;
/*	A(0, 0) = 0.0;  A(0, 1) = 1.0; A(0, 2) = 2.0; 
	A(1, 0) = 0.0;  A(0, 1) = 0.0; A(0, 2) = -1.0;
	A(2, 0) = 0.0;  A(2, 1) = 0.0; A(2, 2) = 0.0;*/

	ublas::matrix<value_type> B(N,N);

	// Initialise the solution matrix B
	B = C;

	MatrixOde ode(A, C);
	double L = 0.0;
	double T = 3.0;
	double dt = 0.001;

	std::size_t steps = Bode::integrate(ode, B, L, T, dt, write_out); 

	std::cout << "Value at time: " << T << '\n';
	std::cout << B << '\n';

	// Computing the exact matrix
	ublas::matrix<value_type> check(N, N);
	double e = std::exp(3.0*T);
	check(0, 0) = (1.0 + 3.0*T)*e;
	check(0, 1) = -T*e;
	check(1, 0) = 9.0*T*e;
	check(1, 1) = (1.0 - 3.0*T)*e;
	std::cout << check << '\n';

	// Sanity check; 1X1 matrix case dB/t = B on (0,5), B(0) = I
	{ // dB/dt = A*B, B(0) = C

		std::size_t N = 1;
		ublas::identity_matrix<value_type> C(N);
		ublas::matrix<value_type> A(N, N);
		A(0, 0) = 1.0;
		ublas::matrix<value_type> B(N, N);

		// Initialise the solution matrix B
		B = C;

		MatrixOde ode(A, C);
		double L = 0.0;
		double T = 5.0;
		double dt = 0.001;

		std::size_t steps = Bode::integrate(ode, B, L, T, dt, write_out);

		std::cout << "Value at time: " << T << '\n';
		std::cout << B << ", " << std::exp(5.0) <<  '\n'; // 148.413
	}

    return 0;
}

# https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html

import numpy as np
from scipy.integrate import odeint


def deriv(A, t, Ab):
    return np.dot(Ab, A)


Ab = np.array([[-0.25,    0,    0],
               [ 0.25, -0.2,    0],
               [    0,  0.2, -0.1]])

time = np.linspace(0, 25, 101)
A0 = np.array([10, 20, 30])

MA = odeint(deriv, A0, time, args=(Ab,))