Skip to content
Snippets Groups Projects
Commit a0820765 authored by patavirt's avatar patavirt
Browse files

Cleanup + edit

parent ecbca958
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/action.hpp
+ 32
− 8
View file @ a0820765
......@@ -14,6 +14,10 @@
#include "common.hpp"
#include "array.hpp"
#define BLOCK_00(mat) ((mat).template block<2,2>(0,0))
#define BLOCK_01(mat) ((mat).template block<2,2>(0,2))
#define BLOCK_10(mat) ((mat).template block<2,2>(2,0))
#define BLOCK_11(mat) ((mat).template block<2,2>(2,2))
template <typename Scalar, typename ScalarW = typename std::remove_const<Scalar>::type>
inline Matrix<ScalarW,4,4>
......@@ -21,7 +25,7 @@ get_Q_matrix(const array::ArrayView<Scalar, array::Shape<2,2,2> >& Q)
{
auto g = Q.part(0u).matrix();
auto gt = Q.part(1u).matrix();
Matrix<double,2,2> I;
Matrix<Complex,2,2> I;
Matrix<ScalarW,2,2> N, Nt;
Matrix<ScalarW,4,4> Qm;
......@@ -31,10 +35,10 @@ get_Q_matrix(const array::ArrayView<Scalar, array::Shape<2,2,2> >& Q)
N = (I + g*gt).inverse();
Nt = (I + gt*g).inverse();
Qm.template block<2,2>(0,0) = N * (I - g * gt);
Qm.template block<2,2>(0,2) = 2 * N * g;
Qm.template block<2,2>(2,0) = 2 * Nt * gt;
Qm.template block<2,2>(2,2) = -Nt * (I - gt * g);
BLOCK_00(Qm) = N * (I - g * gt);
BLOCK_01(Qm) = 2 * N * g;
BLOCK_10(Qm) = 2 * Nt * gt;
BLOCK_11(Qm) = -Nt * (I - gt * g);
return Qm;
}
......@@ -42,6 +46,10 @@ get_Q_matrix(const array::ArrayView<Scalar, array::Shape<2,2,2> >& Q)
template <typename Scalar, typename Shape, typename Shape2, typename ScalarW = typename std::remove_const<Scalar>::type>
inline ScalarW S_2(const array::ArrayView<Scalar, Shape>& Q, const array::ArrayView<const Complex, Shape2>& U, double h)
{
using array::Array;
using array::Shape;
using array::Dynamic;
const int N_i[2] = {0, 1};
const int N_j[2] = {1, 0};
const size_t invk[4] = {3, 2, 1, 0};
......@@ -50,9 +58,18 @@ inline ScalarW S_2(const array::ArrayView<Scalar, Shape>& Q, const array::ArrayV
size_t nx = Q.dim(0);
size_t ny = Q.dim(1);
Array<ScalarW, Shape<Dynamic, Dynamic, 4, 4> > Qmat({nx, ny, 4, 4});
Array<char, Shape<Dynamic, Dynamic> > initQ({nx, ny});
for (size_t i = 0; i < nx; ++i) {
for (size_t j = 0; j < ny; ++j) {
Qmat.part(i,j).matrix() = get_Q_matrix(Q.part(i,j));
}
}
for (size_t i = 0; i < nx; ++i) {
for (size_t j = 0; j < ny; ++j) {
auto Q1 = get_Q_matrix(Q.part(i,j));
auto Q1 = Qmat.part(i,j).matrix();
/* Sum over neighbors */
for (size_t k = 0; k < 2; ++k) {
......@@ -71,11 +88,18 @@ inline ScalarW S_2(const array::ArrayView<Scalar, Shape>& Q, const array::ArrayV
size_t i2 = i + di;
size_t j2 = j + dj;
auto Q2 = get_Q_matrix(Q.part(i2,j2));
auto Q2 = Qmat.part(i2,j2).matrix();
auto U12 = U.part(i,j,k).matrix();
auto U21 = U.part(i2,j2,invk[k]).matrix();
auto M = (Q1 - U12*Q2*U21);
Eigen::Matrix<ScalarW,4,4> M;
/* M = Q1 - U12 * Q * U21; enforce Nambu block diagonal U */
BLOCK_00(M) = BLOCK_00(U12) * BLOCK_00(Q2) * BLOCK_00(U21) - BLOCK_00(Q1);
BLOCK_01(M) = BLOCK_00(U12) * BLOCK_01(Q2) * BLOCK_11(U21) - BLOCK_01(Q1);
BLOCK_10(M) = BLOCK_11(U12) * BLOCK_10(Q2) * BLOCK_00(U21) - BLOCK_10(Q1);
BLOCK_11(M) = BLOCK_11(U12) * BLOCK_11(Q2) * BLOCK_11(U21) - BLOCK_11(Q1);
S += 2/(4*h) * (M * M).trace();
}
}
......
src/common.hpp
+ 17
− 114
View file @ a0820765
......@@ -21,55 +21,16 @@
#include "vector.hpp"
typedef std::complex<double> Complex;
typedef CppAD::AD<double> ADDouble;
typedef std::complex<ADDouble> ADComplex;
class AdeptDouble
{
typedef adept::adouble T;
public:
AdeptDouble() : v_() {}
AdeptDouble(const double& v) : v_(v) {}
AdeptDouble(const T& v) : v_(v) {}
AdeptDouble(const T&& v) : v_(std::move(v)) {}
T& value() { return v_; }
const T& value() const { return v_; }
AdeptDouble& operator+=(const AdeptDouble& other)
{ v_ += other.v_; return *this; }
AdeptDouble& operator-=(const AdeptDouble& other)
{ v_ -= other.v_; return *this; }
AdeptDouble& operator*=(const AdeptDouble& other)
{ v_ *= other.v_; return *this; }
AdeptDouble& operator/=(const AdeptDouble& other)
{ v_ /= other.v_; return *this; }
AdeptDouble operator-() const { return std::move(AdeptDouble(-v_)); }
private:
adept::adouble v_;
};
AdeptDouble operator+(const AdeptDouble &a, const AdeptDouble &b) { return std::move(adept::adouble(a.value() + b.value())); }
AdeptDouble operator-(const AdeptDouble &a, const AdeptDouble &b) { return std::move(adept::adouble(a.value() - b.value())); }
AdeptDouble operator*(const AdeptDouble &a, const AdeptDouble &b) { return std::move(adept::adouble(a.value() * b.value())); }
AdeptDouble operator/(const AdeptDouble &a, const AdeptDouble &b) { return std::move(adept::adouble(a.value() / b.value())); }
typedef std::complex<AdeptDouble> AdeptComplex;
typedef CppAD::AD<Complex> ADComplex;
template <typename Scalar_, int Rows_, int Cols_>
using Matrix = Eigen::Matrix<Scalar_, Rows_, Cols_, Eigen::RowMajor>;
namespace Eigen
{
/* Eigen binop maps: (double, ADDouble) -> ADDouble */
template<typename BinaryOp> struct ScalarBinaryOpTraits<ADDouble,double,BinaryOp> { typedef ADDouble ReturnType; };
template<typename BinaryOp> struct ScalarBinaryOpTraits<double,ADDouble,BinaryOp> { typedef ADDouble ReturnType; };
/* Eigen binop maps: (double/complex, ADComplex) -> ADComplex */
template<typename BinaryOp> struct ScalarBinaryOpTraits<ADComplex,int,BinaryOp> { typedef ADComplex ReturnType; };
template<typename BinaryOp> struct ScalarBinaryOpTraits<int,ADComplex,BinaryOp> { typedef ADComplex ReturnType; };
template<typename BinaryOp> struct ScalarBinaryOpTraits<ADComplex,double,BinaryOp> { typedef ADComplex ReturnType; };
template<typename BinaryOp> struct ScalarBinaryOpTraits<double,ADComplex,BinaryOp> { typedef ADComplex ReturnType; };
template<typename BinaryOp> struct ScalarBinaryOpTraits<ADComplex,Complex,BinaryOp> { typedef ADComplex ReturnType; };
......@@ -77,78 +38,20 @@ namespace Eigen
}
namespace CppAD {
ADComplex adcomplex(const ADComplex&& cpx)
{
return {cpx};
}
ADComplex adcomplex(const std::complex<double>&& cpx)
{
return {cpx};
}
ADComplex operator*(const ADComplex& a, const double& b) { return a * adcomplex(b); }
ADComplex operator*(const double& a, const ADComplex& b) { return adcomplex(a) * b; }
ADComplex operator*(const ADComplex& a, const Complex& b) { return a * adcomplex(b); }
ADComplex operator*(const Complex& a, const ADComplex& b) { return adcomplex(a) * b; }
ADComplex operator+(const ADComplex& a, const double& b) { return a + adcomplex(b); }
ADComplex operator+(const double& a, const ADComplex& b) { return adcomplex(a) + b; }
ADComplex operator+(const ADComplex& a, const Complex& b) { return a + adcomplex(b); }
ADComplex operator+(const Complex& a, const ADComplex& b) { return adcomplex(a) + b; }
ADComplex operator-(const ADComplex& a, const double& b) { return a - adcomplex(b); }
ADComplex operator-(const double& a, const ADComplex& b) { return adcomplex(a) - b; }
ADComplex operator-(const ADComplex& a, const Complex& b) { return a - adcomplex(b); }
ADComplex operator-(const Complex& a, const ADComplex& b) { return adcomplex(a) - b; }
ADComplex operator/(const ADComplex& a, const double& b) { return a / adcomplex(b); }
ADComplex operator/(const double& a, const ADComplex& b) { return adcomplex(a) / b; }
ADComplex operator/(const ADComplex& a, const Complex& b) { return a / adcomplex(b); }
ADComplex operator/(const Complex& a, const ADComplex& b) { return adcomplex(a) / b; }
ADComplex adcomplex(const ADComplex&& cpx) { return {cpx}; }
ADComplex adcomplex(const Complex&& cpx) { return {cpx}; }
#define COMMON_HPP_DECLARE_BINOP(op) \
ADComplex operator op(const ADComplex& a, const double& b) { return a op adcomplex(b); } \
ADComplex operator op(const double& a, const ADComplex& b) { return adcomplex(a) op b; } \
ADComplex operator op(const ADComplex& a, const Complex& b) { return a op adcomplex(b); } \
ADComplex operator op(const Complex& a, const ADComplex& b) { return adcomplex(a) op b; }
COMMON_HPP_DECLARE_BINOP(+);
COMMON_HPP_DECLARE_BINOP(-);
COMMON_HPP_DECLARE_BINOP(*);
COMMON_HPP_DECLARE_BINOP(/);
#undef COMMON_HPP_DECLARE_BINOP
}
namespace Eigen
{
/* Eigen binop maps: (double, ADDouble) -> ADDouble */
template<typename BinaryOp> struct ScalarBinaryOpTraits<AdeptDouble,double,BinaryOp> { typedef AdeptDouble ReturnType; };
template<typename BinaryOp> struct ScalarBinaryOpTraits<double,AdeptDouble,BinaryOp> { typedef AdeptDouble ReturnType; };
/* Eigen binop maps: (double/complex, AdeptComplex) -> AdeptComplex */
template<typename BinaryOp> struct ScalarBinaryOpTraits<AdeptComplex,double,BinaryOp> { typedef AdeptComplex ReturnType; };
template<typename BinaryOp> struct ScalarBinaryOpTraits<double,AdeptComplex,BinaryOp> { typedef AdeptComplex ReturnType; };
template<typename BinaryOp> struct ScalarBinaryOpTraits<AdeptComplex,Complex,BinaryOp> { typedef AdeptComplex ReturnType; };
template<typename BinaryOp> struct ScalarBinaryOpTraits<Complex,AdeptComplex,BinaryOp> { typedef AdeptComplex ReturnType; };
}
AdeptComplex adcomplex(const AdeptComplex&& cpx)
{
return {cpx};
}
AdeptComplex adcomplex(const std::complex<double>&& cpx)
{
return {cpx};
}
AdeptComplex operator*(const AdeptComplex& a, const double& b) { return a * adcomplex(b); }
AdeptComplex operator*(const double& a, const AdeptComplex& b) { return adcomplex(a) * b; }
AdeptComplex operator*(const AdeptComplex& a, const Complex& b) { return a * adcomplex(b); }
AdeptComplex operator*(const Complex& a, const AdeptComplex& b) { return adcomplex(a) * b; }
AdeptComplex operator+(const AdeptComplex& a, const double& b) { return a + adcomplex(b); }
AdeptComplex operator+(const double& a, const AdeptComplex& b) { return adcomplex(a) + b; }
AdeptComplex operator+(const AdeptComplex& a, const Complex& b) { return a + adcomplex(b); }
AdeptComplex operator+(const Complex& a, const AdeptComplex& b) { return adcomplex(a) + b; }
AdeptComplex operator-(const AdeptComplex& a, const double& b) { return a - adcomplex(b); }
AdeptComplex operator-(const double& a, const AdeptComplex& b) { return adcomplex(a) - b; }
AdeptComplex operator-(const AdeptComplex& a, const Complex& b) { return a - adcomplex(b); }
AdeptComplex operator-(const Complex& a, const AdeptComplex& b) { return adcomplex(a) - b; }
AdeptComplex operator/(const AdeptComplex& a, const double& b) { return a / adcomplex(b); }
AdeptComplex operator/(const double& a, const AdeptComplex& b) { return adcomplex(a) / b; }
AdeptComplex operator/(const AdeptComplex& a, const Complex& b) { return a / adcomplex(b); }
AdeptComplex operator/(const Complex& a, const AdeptComplex& b) { return adcomplex(a) / b; }
#endif
src/core.cpp
+ 19
− 45
View file @ a0820765
......@@ -23,23 +23,23 @@ private:
double h_;
typedef Vector<size_t> SizeVector;
typedef Vector<double> DoubleVector;
typedef Vector<Complex> ComplexVector;
typedef std::vector<bool> BoolVector;
typedef CppAD::sparse_rc<SizeVector> Sparsity;
CppAD::ADFun<double> f_;
CppAD::ADFun<Complex> f_;
CppAD::sparse_hes_work hes_work_;
Sparsity hes_pattern_;
CppAD::sparse_rcv<SizeVector, DoubleVector> hes_subset_;
CppAD::sparse_rcv<SizeVector, ComplexVector> hes_subset_;
bool computed_ = false;
pybind11::object py_sparse_matrix_type_;
static Vector<double> ravel(const QType& Q)
static Vector<Complex> ravel(const QType& Q)
{
double *data = reinterpret_cast<double *>(const_cast<Complex *>(Q.data()));
size_t size = 2 * Q.size();
return std::move(Vector<double>(data, size));
Complex *data = const_cast<Complex *>(Q.data());
size_t size = Q.size();
return std::move(Vector<Complex>(data, size));
}
public:
......@@ -57,12 +57,12 @@ public:
return S_2(Q, U_, h_);
}
CppAD::ADFun<double> eval_ad(const QType& Qval)
CppAD::ADFun<Complex> eval_ad(const QType& Qval)
{
if (Qval.dim(0) != nx_ || Qval.dim(1) != ny_)
throw std::out_of_range("Q array wrong size");
Array<ADComplex, Shape<Dynamic,Dynamic,2,2,2>, std::vector<ADDouble> >
Array<ADComplex, Shape<Dynamic,Dynamic,2,2,2> >
Q({nx_, ny_, 2, 2, 2});
auto a = Qval.reshape<Dynamic>({Qval.size()});
......@@ -84,8 +84,8 @@ public:
* more elements than there are variables.
*/
Array<ADDouble, Shape<1> > res;
res(0u) = S_2(Q, U_, h_).real();
Array<ADComplex, Shape<1> > res;
res(0u) = S_2(Q, U_, h_);
return {Q.storage(), res.storage()};
}
......@@ -93,6 +93,7 @@ public:
void compute(const QType& Q)
{
f_ = eval_ad(Q);
f_.optimize();
/*
* Hessian sparsity pattern
......@@ -102,7 +103,7 @@ public:
hes_work_.clear();
/* Right-multiplier is identity matrix */
size_t n = 2*Q.size();
size_t n = Q.size();
size_t nr = n;
size_t nc = n;
size_t nnz_in = n;
......@@ -118,6 +119,7 @@ public:
BoolVector s = {true};
f_.for_jac_sparsity(pattern_in, false, false, false, hes_pattern_);
f_.rev_hes_sparsity(s, false, false, hes_pattern_);
f_.size_forward_set(0);
hes_subset_ = hes_pattern_;
......@@ -126,24 +128,24 @@ public:
auto grad(const QType& Q)
{
Vector<double> x = ravel(Q);
Vector<Complex> x = ravel(Q);
if (!computed_)
compute(Q);
auto jac = f_.Jacobian(x);
return Array<double, Shape<Dynamic>, Vector<double> >(std::move(jac), {jac.size()});
return Array<Complex, Shape<Dynamic>, Vector<Complex> >(std::move(jac), {jac.size()});
}
auto hess(const QType& Q)
{
Vector<double> x = ravel(Q);
Vector<Complex> x = ravel(Q);
if (!computed_)
compute(Q);
/* Hessian */
DoubleVector w = {1.0};
ComplexVector w = {1.0};
std::string coloring = "cppad.symmetric";
f_.sparse_hes(x, w, hes_subset_, hes_pattern_, coloring, hes_work_);
......@@ -151,39 +153,12 @@ public:
Array<size_t, Shape<Dynamic>, SizeVector> row(std::move(SizeVector(hes_subset_.row())), {nnz});
Array<size_t, Shape<Dynamic>, SizeVector> col(std::move(SizeVector(hes_subset_.col())), {nnz});
Array<double, Shape<Dynamic>, DoubleVector> val(std::move(DoubleVector(hes_subset_.val())), {nnz});
Array<Complex, Shape<Dynamic>, ComplexVector> val(std::move(ComplexVector(hes_subset_.val())), {nnz});
return py_sparse_matrix_type_(
pybind11::make_tuple(std::move(val), pybind11::make_tuple(std::move(row), std::move(col))),
pybind11::make_tuple(x.size(), x.size()));
}
auto grad_adept(const QType& Q)
{
if (Q.dim(0) != nx_ || Q.dim(1) != ny_)
throw std::out_of_range("Q array wrong size");
Vector<double> x = ravel(Q);
size_t n = x.size();
adept::Stack stack;
Array<AdeptComplex, Shape<Dynamic,Dynamic,2,2,2>, std::vector<AdeptDouble> >
Qad({nx_, ny_, 2, 2, 2});
Array<AdeptDouble, Shape<1> > res;
adept::set_values(reinterpret_cast<adept::adouble *>(Qad.storage().data()), n, x.data());
stack.new_recording();
res(0u) = S_2(Qad, U_, h_).real();
stack.independent(reinterpret_cast<adept::adouble *>(Qad.storage().data()), n);
stack.dependent(reinterpret_cast<adept::adouble *>(res.data()), 1);
Array<double, Shape<Dynamic> > jac({n});
stack.jacobian(jac.data());
return jac;
}
};
PYBIND11_MODULE(_core, m)
......@@ -200,6 +175,5 @@ PYBIND11_MODULE(_core, m)
.def("grad", &Action::grad)
.def("hess", &Action::hess)
.def("compute", &Action::compute)
.def("grad_adept", &Action::grad_adept)
;
}
src/xtest.f95
+ 17
− 4
View file @ a0820765
......@@ -2,15 +2,19 @@ program main
use adjac
implicit none
type(adjac_complexan), dimension(50,50,2,2,2) :: Q
double complex, dimension(50,50,4,4,4) :: U
type(adjac_complexan), dimension(:,:,:,:,:), allocatable :: Q
double complex, dimension(:,:,:,:,:), allocatable :: U
type(adjac_complexan), dimension(1) :: res
double complex, dimension(1, 50*50*2*2*2) ::jac
double complex, dimension(:, :), allocatable ::jac
double complex :: v
integer i1, i2, i3, i4, i5, p, p2
allocate(Q(50,50,2,2,2))
allocate(U(50,50,4,4,4))
allocate(jac(1, 50*50*2*2*2))
call adjac_reset()
p = 1
......@@ -122,6 +126,7 @@ contains
logical, dimension(:,:), allocatable :: initQ
type(adjac_complexan), dimension(4,4) :: M, Q1, Q2
double complex, dimension(4,4) :: U12, U21
integer :: i, j, k, i2, j2, nx, ny
......@@ -164,7 +169,15 @@ contains
Q1 = Qmat(i,j,:,:)
Q2 = Qmat(i2,j2,:,:)
M = Q1 - matmul(U(i,j,k,:,:), matmul(Q2, U(i2,j2,invk(k),:,:)))
U12 = U(i,j,k,:,:)
U21 = U(i2,j2,invk(k),:,:)
U12(1:2,3:4) = 0
U12(3:4,1:2) = 0
U21(1:2,3:4) = 0
U21(3:4,1:2) = 0
M = Q1 - matmul(U12, matmul(Q2, U21))
res = res + 2/(4*h) * trace(matmul(M, M))
end do
end do
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment