op_chi2rnd_meat.hpp 3.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199
  1. // Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au)
  2. // Copyright 2008-2016 National ICT Australia (NICTA)
  3. //
  4. // Licensed under the Apache License, Version 2.0 (the "License");
  5. // you may not use this file except in compliance with the License.
  6. // You may obtain a copy of the License at
  7. // http://www.apache.org/licenses/LICENSE-2.0
  8. //
  9. // Unless required by applicable law or agreed to in writing, software
  10. // distributed under the License is distributed on an "AS IS" BASIS,
  11. // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  12. // See the License for the specific language governing permissions and
  13. // limitations under the License.
  14. // ------------------------------------------------------------------------
  15. //! \addtogroup op_chi2rnd
  16. //! @{
  17. template<typename T1>
  18. inline
  19. void
  20. op_chi2rnd::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_chi2rnd>& in)
  21. {
  22. arma_extra_debug_sigprint();
  23. typedef typename T1::elem_type eT;
  24. const Proxy<T1> P(in.m);
  25. if(P.is_alias(out) == false)
  26. {
  27. op_chi2rnd::apply_noalias(out, P);
  28. }
  29. else
  30. {
  31. Mat<eT> tmp;
  32. op_chi2rnd::apply_noalias(tmp, P);
  33. out.steal_mem(tmp);
  34. }
  35. }
  36. template<typename T1>
  37. inline
  38. void
  39. op_chi2rnd::apply_noalias(Mat<typename T1::elem_type>& out, const Proxy<T1>& P)
  40. {
  41. arma_extra_debug_sigprint();
  42. #if defined(ARMA_USE_CXX11)
  43. {
  44. typedef typename T1::elem_type eT;
  45. op_chi2rnd_varying_df<eT> generator;
  46. const uword n_rows = P.get_n_rows();
  47. const uword n_cols = P.get_n_cols();
  48. out.set_size(n_rows, n_cols);
  49. eT* out_mem = out.memptr();
  50. if(Proxy<T1>::use_at == false)
  51. {
  52. const uword N = P.get_n_elem();
  53. typename Proxy<T1>::ea_type Pea = P.get_ea();
  54. for(uword i=0; i<N; ++i)
  55. {
  56. out_mem[i] = generator( Pea[i] );
  57. }
  58. }
  59. else
  60. {
  61. for(uword col=0; col < n_cols; ++col)
  62. for(uword row=0; row < n_rows; ++row)
  63. {
  64. (*out_mem) = generator( P.at(row,col) ); ++out_mem;
  65. }
  66. }
  67. }
  68. #else
  69. {
  70. out.reset();
  71. arma_stop_logic_error("chi2rnd(): C++11 compiler required");
  72. }
  73. #endif
  74. }
  75. template<typename eT>
  76. inline
  77. void
  78. op_chi2rnd::fill_constant_df(Mat<eT>& out, const eT df)
  79. {
  80. arma_extra_debug_sigprint();
  81. #if defined(ARMA_USE_CXX11)
  82. {
  83. if(df > eT(0))
  84. {
  85. typedef std::mt19937_64 motor_type;
  86. typedef std::mt19937_64::result_type seed_type;
  87. typedef std::chi_squared_distribution<eT> distr_type;
  88. motor_type motor; motor.seed( seed_type(arma_rng::randi<int>()) );
  89. distr_type distr(df);
  90. const uword N = out.n_elem;
  91. eT* out_mem = out.memptr();
  92. for(uword i=0; i<N; ++i)
  93. {
  94. out_mem[i] = eT( distr(motor) );
  95. }
  96. }
  97. else
  98. {
  99. out.fill( Datum<eT>::nan );
  100. }
  101. }
  102. #else
  103. {
  104. out.reset();
  105. arma_ignore(df);
  106. arma_stop_logic_error("chi2rnd(): C++11 compiler required");
  107. }
  108. #endif
  109. }
  110. //
  111. #if defined(ARMA_USE_CXX11)
  112. template<typename eT>
  113. inline
  114. op_chi2rnd_varying_df<eT>::~op_chi2rnd_varying_df()
  115. {
  116. arma_extra_debug_sigprint();
  117. }
  118. template<typename eT>
  119. inline
  120. op_chi2rnd_varying_df<eT>::op_chi2rnd_varying_df()
  121. {
  122. arma_extra_debug_sigprint();
  123. typedef std::mt19937_64::result_type seed_type;
  124. motor.seed( seed_type(arma_rng::randi<int>()) );
  125. }
  126. template<typename eT>
  127. inline
  128. eT
  129. op_chi2rnd_varying_df<eT>::operator()(const eT df)
  130. {
  131. arma_extra_debug_sigprint();
  132. // as C++11 doesn't seem to provide a way to explicitly set the parameter
  133. // of an existing chi_squared_distribution object,
  134. // we need to create a new object each time
  135. if(df > eT(0))
  136. {
  137. std::chi_squared_distribution<eT> distr(df);
  138. return eT( distr(motor) );
  139. }
  140. else
  141. {
  142. return Datum<eT>::nan;
  143. }
  144. }
  145. #endif
  146. //! @}