TMB Documentation: R_inla.hpp Source File

R_inla.hpp

1 // Copyright (C) 2015 Hans Skaug

2 // License: GPL-2

3 

16  namespace R_inla {

17 using namespace Eigen;

18 using namespace tmbutils;

19 

21 template<class Type>

22  struct spde_t{

23  SparseMatrix<Type> M0; // G0 eqn (10) in Lindgren

24  SparseMatrix<Type> M1; // G1 eqn (10) in Lindgren

25  SparseMatrix<Type> M2; // G2 eqn (10) in Lindgren

26  spde_t(SEXP x){ /* x = List passed from R */

27  M0 = asSparseMatrix<Type>(getListElement(x,"M0"));

28  M1 = asSparseMatrix<Type>(getListElement(x,"M1"));

29  M2 = asSparseMatrix<Type>(getListElement(x,"M2"));

30 }

31 };

32 

34 template<class Type>

35   SparseMatrix<Type> Q_spde(spde_t<Type> spde, Type kappa){

36  Type kappa_pow2 = kappa*kappa;

37  Type kappa_pow4 = kappa_pow2*kappa_pow2;

38 

39  return kappa_pow4*spde.M0 + Type(2.0)*kappa_pow2*spde.M1 + spde.M2; // M0=G0, M1=G1, M2=G2

40 }

41 

43 template<class Type>

44  struct spde_aniso_t{

45  int n_s;

46  int n_tri;

47  vector<Type> Tri_Area;

48  matrix<Type> E0;

49  matrix<Type> E1;

50  matrix<Type> E2;

51  matrix<int> TV;

52  SparseMatrix<Type> G0;

53  SparseMatrix<Type> G0_inv;

54  spde_aniso_t(SEXP x){ /* x = List passed from R */

55  n_s = CppAD::Integer(asVector<Type>(getListElement(x,"n_s"))[0]);

56  n_tri = CppAD::Integer(asVector<Type>(getListElement(x,"n_tri"))[0]);

57  Tri_Area = asVector<Type>(getListElement(x,"Tri_Area"));

58  E0 = asMatrix<Type>(getListElement(x,"E0"));

59  E1 = asMatrix<Type>(getListElement(x,"E1"));

60  E2 = asMatrix<Type>(getListElement(x,"E2"));

61  TV = asMatrix<int>(getListElement(x,"TV"));

62  G0 = asSparseMatrix<Type>(getListElement(x,"G0"));

63  G0_inv = asSparseMatrix<Type>(getListElement(x,"G0_inv"));

64 

65 }

66 };

67 

68 

70 template<class Type>

71   SparseMatrix<Type> Q_spde(spde_aniso_t<Type> spde, Type kappa, matrix<Type> H){

72 

73  int i;

74  Type kappa_pow2 = kappa*kappa;

75  Type kappa_pow4 = kappa_pow2*kappa_pow2;

76 

77  int n_s = spde.n_s;

78  int n_tri = spde.n_tri;

79  vector<Type> Tri_Area = spde.Tri_Area;

80  matrix<Type> E0 = spde.E0;

81  matrix<Type> E1 = spde.E1;

82  matrix<Type> E2 = spde.E2;

83  matrix<int> TV = spde.TV;

84  SparseMatrix<Type> G0 = spde.G0;

85  SparseMatrix<Type> G0_inv = spde.G0_inv;

86 

87  //Type H_trace = H(0,0)+H(1,1);

88  //Type H_det = H(0,0)*H(1,1)-H(0,1)*H(1,0);

89  SparseMatrix<Type> G1_aniso(n_s,n_s);

90  SparseMatrix<Type> G2_aniso(n_s,n_s);

91  // Calculate adjugate of H

92  matrix<Type> adj_H(2,2);

93  adj_H(0,0) = H(1,1);

94  adj_H(0,1) = -1 * H(0,1);

95  adj_H(1,0) = -1 * H(1,0);

96  adj_H(1,1) = H(0,0);

97  // Calculate new SPDE matrices

98 

99  // Calculate G1 - pt. 1

100  array<Type> Gtmp(n_tri,3,3);

101  for(i=0; i<n_tri; i++){

102  // 1st line: E0(i,) %*% adjH %*% t(E0(i,)), etc.

103  Gtmp(i,0,0) = (E0(i,0)*(E0(i,0)*adj_H(0,0)+E0(i,1)*adj_H(1,0)) + E0(i,1)*(E0(i,0)*adj_H(0,1)+E0(i,1)*adj_H(1,1))) / (4*Tri_Area(i));

104  Gtmp(i,0,1) = (E1(i,0)*(E0(i,0)*adj_H(0,0)+E0(i,1)*adj_H(1,0)) + E1(i,1)*(E0(i,0)*adj_H(0,1)+E0(i,1)*adj_H(1,1))) / (4*Tri_Area(i));

105  Gtmp(i,0,2) = (E2(i,0)*(E0(i,0)*adj_H(0,0)+E0(i,1)*adj_H(1,0)) + E2(i,1)*(E0(i,0)*adj_H(0,1)+E0(i,1)*adj_H(1,1))) / (4*Tri_Area(i));

106  Gtmp(i,1,1) = (E1(i,0)*(E1(i,0)*adj_H(0,0)+E1(i,1)*adj_H(1,0)) + E1(i,1)*(E1(i,0)*adj_H(0,1)+E1(i,1)*adj_H(1,1))) / (4*Tri_Area(i));

107  Gtmp(i,1,2) = (E2(i,0)*(E1(i,0)*adj_H(0,0)+E1(i,1)*adj_H(1,0)) + E2(i,1)*(E1(i,0)*adj_H(0,1)+E1(i,1)*adj_H(1,1))) / (4*Tri_Area(i));

108  Gtmp(i,2,2) = (E2(i,0)*(E2(i,0)*adj_H(0,0)+E2(i,1)*adj_H(1,0)) + E2(i,1)*(E2(i,0)*adj_H(0,1)+E2(i,1)*adj_H(1,1))) / (4*Tri_Area(i));

109  }

110  // Calculate G1 - pt. 2

111  for(i=0; i<n_tri; i++){

112  G1_aniso.coeffRef(TV(i,1),TV(i,0)) = G1_aniso.coeffRef(TV(i,1),TV(i,0)) + (Gtmp(i,0,1));

113  G1_aniso.coeffRef(TV(i,0),TV(i,1)) = G1_aniso.coeffRef(TV(i,0),TV(i,1)) + (Gtmp(i,0,1));

114  G1_aniso.coeffRef(TV(i,2),TV(i,1)) = G1_aniso.coeffRef(TV(i,2),TV(i,1)) + (Gtmp(i,1,2));

115  G1_aniso.coeffRef(TV(i,1),TV(i,2)) = G1_aniso.coeffRef(TV(i,1),TV(i,2)) + (Gtmp(i,1,2));

116  G1_aniso.coeffRef(TV(i,2),TV(i,0)) = G1_aniso.coeffRef(TV(i,2),TV(i,0)) + (Gtmp(i,0,2));

117  G1_aniso.coeffRef(TV(i,0),TV(i,2)) = G1_aniso.coeffRef(TV(i,0),TV(i,2)) + (Gtmp(i,0,2));

118  G1_aniso.coeffRef(TV(i,0),TV(i,0)) = G1_aniso.coeffRef(TV(i,0),TV(i,0)) + (Gtmp(i,0,0));

119  G1_aniso.coeffRef(TV(i,1),TV(i,1)) = G1_aniso.coeffRef(TV(i,1),TV(i,1)) + (Gtmp(i,1,1));

120  G1_aniso.coeffRef(TV(i,2),TV(i,2)) = G1_aniso.coeffRef(TV(i,2),TV(i,2)) + (Gtmp(i,2,2));

121  }

122  G2_aniso = G1_aniso * G0_inv * G1_aniso;

123 

124  return kappa_pow4*G0 + Type(2.0)*kappa_pow2*G1_aniso + G2_aniso;

125 }

126 

127 } // end namespace R_inla

128 

129