C++ std::log方法代碼示例

double CMT::ExponentialFunction::inverse(double data) const {
	return log(data - mEpsilon);
}

void colorize(const Image<uint32_t> &input, const Image<uint32_t> &strokes, Image<uint32_t> &output) {
    int w     = input.width();
    int h     = input.height();
    int x,y;

    int max_d          = floor(log(min(h,w))/log(2)-2);
    float scale_factor = 1.0f/( pow(2,max_d-1) );
    int padded_w       = ceil(w*scale_factor)*pow(2,max_d-1);
    int padded_h       = ceil(h*scale_factor)*pow(2,max_d-1);

    // RGB 2 YUV and padarray
    Image<float> yuv_fused(padded_w,padded_h,3);
    hl_fuse_yuv(input, strokes, yuv_fused);


    // Extract Strokes mask
    Image<float> stroke_mask(padded_w, padded_h);
    hl_nonzero(strokes,stroke_mask);

    Image<float> result(padded_w,padded_h,3);

    int n = padded_h;
    int m = padded_w;
    int k = 1;

    Tensor3d D,G,I;
    Tensor3d Dx,Dy,iDx,iDy;
    MG smk;
    G.set(n,m,k);
    D.set(n,m,k);
    I.set(n,m,k);

    int in_itr_num  = 5;
    int out_itr_num = 1;

    Dx.set(n,m,k-1);
    Dy.set(n,m,k-1);
    iDx.set(n,m,k-1);
    iDy.set(n,m,k-1);

    // Fill in label mask and luminance channels
    for ( y = 0; y<n; y++){
        for ( x = 0; x<m; x++){
            I(y,x,0) = stroke_mask(x,y);
            G(y,x,0) = yuv_fused(x,y,0);
            I(y,x,0) = !I(y,x,0);
        }
    }

    // Write output luminance
    for ( y=0; y<n; y++){
        for ( x=0; x<m; x++){
            result(x,y,0)=G(y,x,0);
        }
    }

    smk.set(n,m,k,max_d);
    smk.setI(I) ;
    smk.setG(G);
    smk.setFlow(Dx,Dy,iDx,iDy);

    // Solve chrominance
    for (int t=1; t<3; t++){
        for ( y=0; y<n; y++){
            for ( x=0; x<m; x++){
                D(y,x,0)       = yuv_fused(x,y,t);
                smk.P()(y,x,0) = yuv_fused(x,y,t);
                D(y,x,0)      *= (!I(y,x,0));
            }
        }

        smk.Div() = D ;

        Tensor3d tP2;

        for (int itr=0; itr<out_itr_num; itr++){
            smk.setDepth(max_d);
            Field_MGN(&smk, in_itr_num, 2) ;
            smk.setDepth(ceil(max_d/2));
            Field_MGN(&smk, in_itr_num, 2) ;
            smk.setDepth(2);
            Field_MGN(&smk, in_itr_num, 2) ;
            smk.setDepth(1);
            Field_MGN(&smk, in_itr_num, 4) ;
        }

        tP2 = smk.P();

        for ( y=0; y<n; y++){
            for ( x=0; x<m; x++){
                result(x,y,t) = tP2(y,x,0);
            }
        }
    }
    
    hl_yuv2rgb(result,output);
    

}

    typename return_type<T_y, T_loc, T_scale>::type
    normal_ccdf_log(const T_y& y, const T_loc& mu, const T_scale& sigma) {
      static const char* function("stan::math::normal_ccdf_log");
      typedef typename stan::partials_return_type<T_y, T_loc, T_scale>::type
        T_partials_return;

      using stan::math::check_positive;
      using stan::math::check_finite;
      using stan::math::check_not_nan;
      using stan::math::check_consistent_sizes;
      using stan::math::value_of;
      using stan::math::INV_SQRT_2;
      using std::log;
      using std::exp;

      T_partials_return ccdf_log(0.0);
      // check if any vectors are zero length
      if (!(stan::length(y)
            && stan::length(mu)
            && stan::length(sigma)))
        return ccdf_log;

      check_not_nan(function, "Random variable", y);
      check_finite(function, "Location parameter", mu);
      check_not_nan(function, "Scale parameter", sigma);
      check_positive(function, "Scale parameter", sigma);
      check_consistent_sizes(function,
                             "Random variable", y,
                             "Location parameter", mu,
                             "Scale parameter", sigma);

      OperandsAndPartials<T_y, T_loc, T_scale>
        operands_and_partials(y, mu, sigma);

      VectorView<const T_y> y_vec(y);
      VectorView<const T_loc> mu_vec(mu);
      VectorView<const T_scale> sigma_vec(sigma);
      size_t N = max_size(y, mu, sigma);
      double log_half = std::log(0.5);

      const double SQRT_TWO_OVER_PI = std::sqrt(2.0 / stan::math::pi());
      for (size_t n = 0; n < N; n++) {
        const T_partials_return y_dbl = value_of(y_vec[n]);
        const T_partials_return mu_dbl = value_of(mu_vec[n]);
        const T_partials_return sigma_dbl = value_of(sigma_vec[n]);

        const T_partials_return scaled_diff = (y_dbl - mu_dbl)
          / (sigma_dbl * SQRT_2);

        T_partials_return one_m_erf;
        if (scaled_diff < -37.5 * INV_SQRT_2)
          one_m_erf = 2.0;
        else if (scaled_diff < -5.0 * INV_SQRT_2)
          one_m_erf =  2.0 - erfc(-scaled_diff);
        else if (scaled_diff > 8.25 * INV_SQRT_2)
          one_m_erf = 0.0;
        else
          one_m_erf = 1.0 - erf(scaled_diff);

        // log ccdf
        ccdf_log += log_half + log(one_m_erf);

        // gradients
        if (contains_nonconstant_struct<T_y, T_loc, T_scale>::value) {
          const T_partials_return rep_deriv_div_sigma
            = scaled_diff > 8.25 * INV_SQRT_2
            ? std::numeric_limits<double>::infinity()
            : SQRT_TWO_OVER_PI * exp(-scaled_diff * scaled_diff)
            / one_m_erf / sigma_dbl;
          if (!is_constant_struct<T_y>::value)
            operands_and_partials.d_x1[n] -= rep_deriv_div_sigma;
          if (!is_constant_struct<T_loc>::value)
            operands_and_partials.d_x2[n] += rep_deriv_div_sigma;
          if (!is_constant_struct<T_scale>::value)
            operands_and_partials.d_x3[n] += rep_deriv_div_sigma
              * scaled_diff * stan::math::SQRT_2;
        }
      }
      return operands_and_partials.value(ccdf_log);
    }

TEST(prob_transform, lb_f) {
  EXPECT_FLOAT_EQ(log(3.0 - 2.0), stan::prob::lb_free(3.0,2.0));
  EXPECT_FLOAT_EQ(1.7, stan::prob::lb_free(1.7, -std::numeric_limits<double>::infinity()));
}

    typename return_type<T_y, T_loc, T_scale>::type
    lognormal_log(const T_y& y, const T_loc& mu, const T_scale& sigma) {
      static const char* function("stan::math::lognormal_log");
      typedef typename stan::partials_return_type<T_y, T_loc, T_scale>::type
        T_partials_return;

      using stan::is_constant_struct;
      using stan::math::check_not_nan;
      using stan::math::check_finite;
      using stan::math::check_positive_finite;
      using stan::math::check_nonnegative;
      using stan::math::check_consistent_sizes;
      using stan::math::value_of;
      using stan::math::include_summand;


      // check if any vectors are zero length
      if (!(stan::length(y)
            && stan::length(mu)
            && stan::length(sigma)))
        return 0.0;

      // set up return value accumulator
      T_partials_return logp(0.0);

      // validate args (here done over var, which should be OK)
      check_not_nan(function, "Random variable", y);
      check_nonnegative(function, "Random variable", y);
      check_finite(function, "Location parameter", mu);
      check_positive_finite(function, "Scale parameter", sigma);
      check_consistent_sizes(function,
                             "Random variable", y,
                             "Location parameter", mu,
                             "Scale parameter", sigma);

      VectorView<const T_y> y_vec(y);
      VectorView<const T_loc> mu_vec(mu);
      VectorView<const T_scale> sigma_vec(sigma);
      size_t N = max_size(y, mu, sigma);

      for (size_t n = 0; n < length(y); n++)
        if (value_of(y_vec[n]) <= 0)
          return LOG_ZERO;

      OperandsAndPartials<T_y, T_loc, T_scale>
        operands_and_partials(y, mu, sigma);

      using stan::math::square;
      using std::log;
      using stan::math::NEG_LOG_SQRT_TWO_PI;
      using std::log;


      VectorBuilder<include_summand<propto, T_scale>::value,
                    T_partials_return, T_scale> log_sigma(length(sigma));
      if (include_summand<propto, T_scale>::value) {
        for (size_t n = 0; n < length(sigma); n++)
          log_sigma[n] = log(value_of(sigma_vec[n]));
      }

      VectorBuilder<include_summand<propto, T_y, T_loc, T_scale>::value,
                    T_partials_return, T_scale> inv_sigma(length(sigma));
      VectorBuilder<include_summand<propto, T_y, T_loc, T_scale>::value,
                    T_partials_return, T_scale> inv_sigma_sq(length(sigma));
      if (include_summand<propto, T_y, T_loc, T_scale>::value) {
        for (size_t n = 0; n < length(sigma); n++)
          inv_sigma[n] = 1 / value_of(sigma_vec[n]);
      }
      if (include_summand<propto, T_y, T_loc, T_scale>::value) {
        for (size_t n = 0; n < length(sigma); n++)
          inv_sigma_sq[n] = inv_sigma[n] * inv_sigma[n];
      }

      VectorBuilder<include_summand<propto, T_y, T_loc, T_scale>::value,
                    T_partials_return, T_y> log_y(length(y));
      if (include_summand<propto, T_y, T_loc, T_scale>::value) {
        for (size_t n = 0; n < length(y); n++)
          log_y[n] = log(value_of(y_vec[n]));
      }

      VectorBuilder<!is_constant_struct<T_y>::value,
                    T_partials_return, T_y> inv_y(length(y));
      if (!is_constant_struct<T_y>::value) {
        for (size_t n = 0; n < length(y); n++)
          inv_y[n] = 1 / value_of(y_vec[n]);
      }

      if (include_summand<propto>::value)
        logp += N * NEG_LOG_SQRT_TWO_PI;

      for (size_t n = 0; n < N; n++) {
        const T_partials_return mu_dbl = value_of(mu_vec[n]);

        T_partials_return logy_m_mu(0);
        if (include_summand<propto, T_y, T_loc, T_scale>::value)
          logy_m_mu = log_y[n] - mu_dbl;

        T_partials_return logy_m_mu_sq = logy_m_mu * logy_m_mu;
        T_partials_return logy_m_mu_div_sigma(0);
        if (contains_nonconstant_struct<T_y, T_loc, T_scale>::value)
//.........這裏部分代碼省略.........

 inline typename boost::math::tools::promote_args<T>::type
 log2(const T a) {
   using std::log;
   return log(a) / LOG_2;
 }

    typename return_type<T_y,T_loc,T_scale,T_shape>::type
    skew_normal_log(const T_y& y, const T_loc& mu, const T_scale& sigma, 
                    const T_shape& alpha) {
      static const char* function("stan::prob::skew_normal_log");
      typedef typename stan::partials_return_type<T_y,T_loc,
                                                  T_scale,T_shape>::type 
        T_partials_return;

      using std::log;
      using stan::is_constant_struct;
      using stan::math::check_positive;
      using stan::math::check_finite;
      using stan::math::check_not_nan;
      using stan::math::check_consistent_sizes;
      using stan::math::value_of;
      using stan::prob::include_summand;

      // check if any vectors are zero length
      if (!(stan::length(y) 
            && stan::length(mu) 
            && stan::length(sigma)
            && stan::length(alpha)))
        return 0.0;

      // set up return value accumulator
      T_partials_return logp(0.0);

      // validate args (here done over var, which should be OK)
      check_not_nan(function, "Random variable", y);
      check_finite(function, "Location parameter", mu);
      check_finite(function, "Shape parameter", alpha);
      check_positive(function, "Scale parameter", sigma);
      check_consistent_sizes(function,
                             "Random variable", y,
                             "Location parameter", mu,
                             "Scale parameter", sigma,
                             "Shape paramter", alpha);

      // check if no variables are involved and prop-to
      if (!include_summand<propto,T_y,T_loc,T_scale,T_shape>::value)
        return 0.0;
      
      // set up template expressions wrapping scalars into vector views
      agrad::OperandsAndPartials<T_y, T_loc, T_scale, T_shape> 
        operands_and_partials(y, mu, sigma, alpha);

      using boost::math::erfc;
      using boost::math::erf;

      VectorView<const T_y> y_vec(y);
      VectorView<const T_loc> mu_vec(mu);
      VectorView<const T_scale> sigma_vec(sigma);
      VectorView<const T_shape> alpha_vec(alpha);
      size_t N = max_size(y, mu, sigma, alpha);

      VectorBuilder<true, T_partials_return, T_scale> inv_sigma(length(sigma));
      VectorBuilder<include_summand<propto,T_scale>::value,
                    T_partials_return, T_scale> log_sigma(length(sigma));
      for (size_t i = 0; i < length(sigma); i++) {
        inv_sigma[i] = 1.0 / value_of(sigma_vec[i]);
        if (include_summand<propto,T_scale>::value)
          log_sigma[i] = log(value_of(sigma_vec[i]));
      }

      for (size_t n = 0; n < N; n++) {
        // pull out values of arguments
        const T_partials_return y_dbl = value_of(y_vec[n]);
        const T_partials_return mu_dbl = value_of(mu_vec[n]);
        const T_partials_return sigma_dbl = value_of(sigma_vec[n]);
        const T_partials_return alpha_dbl = value_of(alpha_vec[n]);

        // reusable subexpression values
        const T_partials_return y_minus_mu_over_sigma 
          = (y_dbl - mu_dbl) * inv_sigma[n];
        const double pi_dbl = stan::math::pi();

        // log probability
        if (include_summand<propto>::value)
          logp -=  0.5 * log(2.0 * pi_dbl);
        if (include_summand<propto, T_scale>::value)
          logp -= log(sigma_dbl);
        if (include_summand<propto,T_y, T_loc, T_scale>::value)
          logp -= y_minus_mu_over_sigma * y_minus_mu_over_sigma / 2.0;
        if (include_summand<propto,T_y,T_loc,T_scale,T_shape>::value)
          logp += log(erfc(-alpha_dbl * y_minus_mu_over_sigma 
                           / std::sqrt(2.0)));

        // gradients
        T_partials_return deriv_logerf 
          = 2.0 / std::sqrt(pi_dbl) 
          * exp(-alpha_dbl * y_minus_mu_over_sigma / std::sqrt(2.0) 
                * alpha_dbl * y_minus_mu_over_sigma / std::sqrt(2.0))
          / (1 + erf(alpha_dbl * y_minus_mu_over_sigma 
                     / std::sqrt(2.0)));
        if (!is_constant_struct<T_y>::value)
          operands_and_partials.d_x1[n] 
            += -y_minus_mu_over_sigma / sigma_dbl 
            + deriv_logerf * alpha_dbl / (sigma_dbl * std::sqrt(2.0)) ;
        if (!is_constant_struct<T_loc>::value)
          operands_and_partials.d_x2[n] 
//.........這裏部分代碼省略.........

template <typename Scalar> Scalar log2(Scalar v) { using std::log; return log(v)/log(Scalar(2)); }

    typename return_type<T_y,T_loc,T_scale>::type
    double_exponential_ccdf_log(const T_y& y, const T_loc& mu, 
                               const T_scale& sigma) {
      static const char* function
        = "stan::prob::double_exponential_ccdf_log(%1%)";
      
      using stan::math::check_finite;
      using stan::math::check_not_nan;
      using stan::math::check_positive;
      using stan::math::check_consistent_sizes;
      using stan::math::value_of;

      double ccdf_log(0.0);

      // check if any vectors are zero length
      if (!(stan::length(y) 
            && stan::length(mu) 
            && stan::length(sigma)))
        return ccdf_log;

      if(!check_not_nan(function, y, "Random variable", &ccdf_log))
        return ccdf_log;
      if(!check_finite(function, mu, "Location parameter", &ccdf_log))
        return ccdf_log;
      if(!check_finite(function, sigma, "Scale parameter", &ccdf_log))
        return ccdf_log;
      if(!check_positive(function, sigma, "Scale parameter", &ccdf_log))
        return ccdf_log;
      if (!(check_consistent_sizes(function, y, mu, sigma,
                                   "Random variable", "Location parameter", 
                                   "Scale Parameter", &ccdf_log)))
        return ccdf_log;
      
      using std::log;
      using std::exp;

      agrad::OperandsAndPartials<T_y, T_loc, T_scale> 
        operands_and_partials(y, mu, sigma);

      VectorView<const T_y> y_vec(y);
      VectorView<const T_loc> mu_vec(mu);
      VectorView<const T_scale> sigma_vec(sigma);
      const double log_half = std::log(0.5);
      size_t N = max_size(y, mu, sigma);

      for (size_t n = 0; n < N; n++) {
        const double y_dbl = value_of(y_vec[n]);
        const double mu_dbl = value_of(mu_vec[n]);
        const double sigma_dbl = value_of(sigma_vec[n]);
        const double scaled_diff = (y_dbl - mu_dbl) / sigma_dbl;
        const double inv_sigma = 1.0 / sigma_dbl;
        if(y_dbl < mu_dbl) {
          //log ccdf
          ccdf_log += log(1.0 - 0.5 * exp(scaled_diff));

          //gradients
          const double rep_deriv = 1.0 / (2.0 * exp(-scaled_diff) - 1.0);
        if (!is_constant_struct<T_y>::value)
          operands_and_partials.d_x1[n] -= rep_deriv * inv_sigma;
        if (!is_constant_struct<T_loc>::value)
          operands_and_partials.d_x2[n] += rep_deriv * inv_sigma;
        if (!is_constant_struct<T_scale>::value)
          operands_and_partials.d_x3[n] += rep_deriv * scaled_diff 
            * inv_sigma;
        }
        else {
          // log ccdf
          ccdf_log += log_half - scaled_diff;

          // gradients
          if (!is_constant_struct<T_y>::value)
            operands_and_partials.d_x1[n] -= inv_sigma;
          if (!is_constant_struct<T_loc>::value)
            operands_and_partials.d_x2[n] += inv_sigma;
          if (!is_constant_struct<T_scale>::value)
            operands_and_partials.d_x3[n] += scaled_diff * inv_sigma;
        }
      }
      return operands_and_partials.to_var(ccdf_log);
    }

 /**
  * Returns the base two logarithm of the argument (C99, C++11).
  *
  * The function is defined by:
  *
  * <code>log2(a) = log(a) / std::log(2.0)</code>.
  *
  * @param[in] u argument
  * @return base two logarithm of argument
  */
 inline double log2(double u) {
   using std::log;
   return log(u) / LOG_2;
 }

void BaseReverseMode<Base>::reverse_local_computation(size_t ind_num, size_t dep_num) {
    using std::sin;
    using std::cos;
    using std::acos;
    using std::sqrt;
    using std::pow;
    using std::log;
    using std::exp;

    DerivativeInfo<locint, Base> info;
    if (!trace) {
        warning_NoTraceSet();
    }
    if (ind_num != trace->get_num_ind()) {
        warning_NumberInconsistent("independent", ind_num, trace->get_num_ind());
    }
    if (dep_num != trace->get_num_dep()) {
        warning_NumberInconsistent("dependent", dep_num, trace->get_num_dep());
    }
    size_t ind_count = trace->get_num_ind();
    size_t dep_count = trace->get_num_dep();

    locint res;
    double coval;
    trace->init_reverse();
    opbyte op = trace->get_next_op_r();

    while (op != start_of_tape) {
        info.clear();
        info.opcode = op;
        switch (op) {
        case start_of_tape:
        case end_of_tape:
            break;
        case assign_ind:
            res = trace->get_next_loc_r();;
            trace->get_next_val_r();
            if (ind_count == 0) {
                // TODO(warning)
            }
            ind_count--;
            indep_index_map[res] = ind_count;
            break;
        case assign_dep:
            res = trace->get_next_loc_r();
            dep_value[res] = trace->get_next_val_r();
            init_dep_deriv(res);
            if (dep_count == 0) {
                // TODO(warning)
            }
            dep_count--;
            dep_index_map[res] = dep_count;
            break;
        case assign_param:
            info.r = trace->get_next_loc_r();
            trace->get_next_param_r();
            break;
        case assign_d:
            info.r = trace->get_next_loc_r();
            trace->get_next_coval_r();
            break;
        case assign_a:
            info.r = trace->get_next_loc_r();
            info.x = trace->get_next_loc_r();
            info.dx = 1.0;
            break;
        case comp_eq:
        case comp_lt:
            trace->get_next_loc_r();
            trace->get_next_loc_r();
            trace->get_next_coval_r();
            break;
        case eq_plus_a:
        case plus_a_a:
            info.r = trace->get_next_loc_r();
            info.y = trace->get_next_loc_r();
            info.x = trace->get_next_loc_r();
            info.dx = 1.0;
            info.dy = 1.0;
            PSEUDO_BINARY
            break;
        case eq_plus_d:
        case plus_d_a:
            info.r = trace->get_next_loc_r();
            info.x = trace->get_next_loc_r();
            trace->get_next_coval_r();
            info.dx = 1.0;
            break;
        case eq_minus_a:
        case minus_a_a:
            info.r = trace->get_next_loc_r();
            info.y = trace->get_next_loc_r();
            info.x = trace->get_next_loc_r();
            info.dx = 1.0;
            info.dy = -1.0;
            PSEUDO_BINARY
            break;
        case minus_d_a:
            info.r = trace->get_next_loc_r();
            info.x = trace->get_next_loc_r();
//.........這裏部分代碼省略.........

    typename return_type<T_y, T_scale_succ, T_scale_fail>::type
    beta_ccdf_log(const T_y& y, const T_scale_succ& alpha,
                  const T_scale_fail& beta) {
      typedef typename stan::partials_return_type<T_y, T_scale_succ,
                                                  T_scale_fail>::type
        T_partials_return;

      // Size checks
      if ( !( stan::length(y) && stan::length(alpha)
              && stan::length(beta) ) )
        return 0.0;

      // Error checks
      static const char* function("stan::math::beta_cdf");

      using stan::math::check_positive_finite;
      using stan::math::check_not_nan;
      using stan::math::check_nonnegative;
      using stan::math::check_less_or_equal;
      using boost::math::tools::promote_args;
      using stan::math::check_consistent_sizes;
      using stan::math::value_of;

      T_partials_return ccdf_log(0.0);

      check_positive_finite(function, "First shape parameter", alpha);
      check_positive_finite(function, "Second shape parameter", beta);
      check_not_nan(function, "Random variable", y);
      check_nonnegative(function, "Random variable", y);
      check_less_or_equal(function, "Random variable", y, 1);
      check_consistent_sizes(function,
                             "Random variable", y,
                             "First shape parameter", alpha,
                             "Second shape parameter", beta);

      // Wrap arguments in vectors
      VectorView<const T_y> y_vec(y);
      VectorView<const T_scale_succ> alpha_vec(alpha);
      VectorView<const T_scale_fail> beta_vec(beta);
      size_t N = max_size(y, alpha, beta);

      OperandsAndPartials<T_y, T_scale_succ, T_scale_fail>
        operands_and_partials(y, alpha, beta);

      // Compute CDF and its gradients
      using stan::math::inc_beta;
      using stan::math::digamma;
      using stan::math::lbeta;
      using std::pow;
      using std::exp;
      using std::log;
      using std::exp;

      // Cache a few expensive function calls if alpha or beta is a parameter
      VectorBuilder<contains_nonconstant_struct<T_scale_succ,
                                                T_scale_fail>::value,
                    T_partials_return, T_scale_succ, T_scale_fail>
        digamma_alpha_vec(max_size(alpha, beta));
      VectorBuilder<contains_nonconstant_struct<T_scale_succ,
                                                T_scale_fail>::value,
                    T_partials_return, T_scale_succ, T_scale_fail>
        digamma_beta_vec(max_size(alpha, beta));
      VectorBuilder<contains_nonconstant_struct<T_scale_succ,
                                                T_scale_fail>::value,
                    T_partials_return, T_scale_succ, T_scale_fail>
        digamma_sum_vec(max_size(alpha, beta));

      if (contains_nonconstant_struct<T_scale_succ, T_scale_fail>::value) {
        for (size_t i = 0; i < N; i++) {
          const T_partials_return alpha_dbl = value_of(alpha_vec[i]);
          const T_partials_return beta_dbl = value_of(beta_vec[i]);

          digamma_alpha_vec[i] = digamma(alpha_dbl);
          digamma_beta_vec[i] = digamma(beta_dbl);
          digamma_sum_vec[i] = digamma(alpha_dbl + beta_dbl);
        }
      }

      // Compute vectorized CDFLog and gradient
      for (size_t n = 0; n < N; n++) {
        // Pull out values
        const T_partials_return y_dbl = value_of(y_vec[n]);
        const T_partials_return alpha_dbl = value_of(alpha_vec[n]);
        const T_partials_return beta_dbl = value_of(beta_vec[n]);
        const T_partials_return betafunc_dbl = exp(lbeta(alpha_dbl, beta_dbl));

        // Compute
        const T_partials_return Pn = 1.0 - inc_beta(alpha_dbl, beta_dbl, y_dbl);

        ccdf_log += log(Pn);

        if (!is_constant_struct<T_y>::value)
          operands_and_partials.d_x1[n] -= pow(1-y_dbl, beta_dbl-1)
            * pow(y_dbl, alpha_dbl-1) / betafunc_dbl / Pn;

        T_partials_return g1 = 0;
        T_partials_return g2 = 0;

        if (contains_nonconstant_struct<T_scale_succ, T_scale_fail>::value) {
          stan::math::grad_reg_inc_beta(g1, g2, alpha_dbl, beta_dbl, y_dbl,
//.........這裏部分代碼省略.........

bool fabric_t::build_lid_map(bool determine_lmc)
{
  ///Always start clean
  clear_lidmap();
  
  {
    const lmc_t max_lmc_lid = lmc > 0 ? (1 << lmc) - 1 : 0;
    
    /**
    * Walk every entity and build lid map
    */
    for(
      entities_t::iterator 
        itr = entities.begin(),
        eitr = entities.end();
      itr != eitr;
      ++itr
    )
    {
      assert(itr->second.lid() > 0);
      entity_t &entity = itr->second;
      const lid_t blid = entity.lid();
      assert(blid > 0);
      
      if(entity.get_type() == port_type::HCA)
        for(lmc_t i = 0; i <= max_lmc_lid; ++i)
        {
#ifndef NDEBUG
          std::cerr << "set HCA lid " << entity.label() << "(" << itr->first << std::hex << ") = " << regex::string_cast_uint(blid + i) << std::endl;
#endif
          assert(lidmap.find(blid + i) == lidmap.end());
          lidmap[blid + i] = &entity;
        }
      else 
        if(entity.get_type() == port_type::TCA) 
        { ///Switchs do not get a second LID
#ifndef NDEBUG
          {
            entitiesmap_lid_t::iterator itr = lidmap.find(blid);
            if(itr != lidmap.end())
            {
              std::cerr << "attempting fabric lmc = " << regex::string_cast_uint(lmc) << std::endl;
              std::cerr << "found existing port " << itr->second->label() <<  " on lid " << blid << std::endl;
              std::cerr << "was going to set port " << entity.label() <<  " on lid " << blid << std::endl;
              abort();
            }
          }
          std::cerr << "set TCA lid " << entity.label() << "(" << itr->first << std::hex << ") = " << regex::string_cast_uint(blid) << std::endl;
#endif
          lidmap[blid] = &entity;
        }
      else
        abort(); ///unknown port type?
    }
  }
  
  /**
    * attempt to determine LMC value of the subnet
    * this can be done with reasonable accuracy since
    * all lmc lid values are sequential for lmc > 0
    * 
    * this is the brute force solution O(ports)*lmc
    * basically, walk every port and see if there are any other lid+lmc
    * until the smalled lid, lid+lmc*, lid sequence is found 
    * then use lmc=log2(found)
    * 
    * LIDs = BASELID to BASELID + (2^LMC - 1)
    */
  if(determine_lmc)
  {
    using std::log;
    const lmc_t current_lmc = lmc;
    assert(portmap.size() > 1);
    
    ///Start off assuming max LMC value
    lmc_t max_lmc_lid = (1 << MAX_LMC_VALUE) - 1;
    entitiesmap_lid_t::const_iterator lid_end = lidmap.end();
    
    for(
      portmap_guidport_t::const_iterator
        itr = portmap.begin(),
        eitr = portmap.end();
      itr != eitr && max_lmc_lid > 0;
      ++itr
    )
      ///Only search LIDs of HCAs
      if(itr->second->type == port_type::HCA)
      {
#ifndef NDEBUG
        std::cerr << "search port " << itr->second->label() << std::endl;
#endif

        ///walk until highest seen lmc value offset
        for(lmc_t i = 1; i <= max_lmc_lid; ++i)
        {
          assert(itr->second->lid > 0);
          assert(lidmap.find(itr->second->lid) != lid_end);
          
          ///is there lid on base lid + lmc offset
          if(lidmap.find(itr->second->lid + i) != lid_end)
//.........這裏部分代碼省略.........

bool CMT::Mixture::train(
	const MatrixXd& data,
	const Parameters& parameters,
	const Component::Parameters& componentParameters)
{
	if(data.rows() != dim())
		throw Exception("Data has wrong dimensionality.");

	if(parameters.initialize && !initialized())
		initialize(data, parameters, componentParameters);

	ArrayXXd logJoint(numComponents(), data.cols());
	Array<double, Dynamic, 1> postSum;
	Array<double, 1, Dynamic> logLik;
	ArrayXXd post;
	ArrayXXd weights;
	double avgLogLoss = numeric_limits<double>::infinity();
	double avgLogLossNew;

	for(int i = 0; i < parameters.maxIter; ++i) {
		// compute joint probability of data and assignments (E)
		#pragma omp parallel for
		for(int k = 0; k < numComponents(); ++k)
			logJoint.row(k) = mComponents[k]->logLikelihood(data) + log(mPriors[k]);

		// compute normalized posterior (E)
		logLik = logSumExp(logJoint);

		// average negative log-likelihood in bits per component
		avgLogLossNew = -logLik.mean() / log(2.) / dim();

		if(parameters.verbosity > 0)
			cout << setw(6) << i << setw(14) << setprecision(7) << avgLogLossNew << endl;

		// test for convergence
		if(avgLogLoss - avgLogLossNew < parameters.threshold)
			return true;
		avgLogLoss = avgLogLossNew;

		// compute normalized posterior (E)
		post = (logJoint.rowwise() - logLik).exp();
		postSum = post.rowwise().sum();
		weights = post.colwise() / postSum;

		// optimize prior weights (M)
		if(parameters.trainPriors) {
			mPriors = postSum / data.cols() + parameters.regularizePriors;
			mPriors /= mPriors.sum();
		}

		// optimize components (M)
		if(parameters.trainComponents) {
			#pragma omp parallel for
			for(int k = 0; k < numComponents(); ++k)
				mComponents[k]->train(data, weights.row(k), componentParameters);
		} else {
			return true;
		}
	}

	if(parameters.verbosity > 0)
		cout << setw(6) << parameters.maxIter << setw(14) << setprecision(7) << evaluate(data) << endl;

	return false;
}

    typename return_type<T_y, T_scale>::type
    rayleigh_log(const T_y& y, const T_scale& sigma) {
      static const char* function("stan::math::rayleigh_log");
      typedef typename stan::partials_return_type<T_y, T_scale>::type
        T_partials_return;

      using std::log;
      using stan::is_constant_struct;
      using stan::math::check_positive;
      using stan::math::check_not_nan;
      using stan::math::check_consistent_sizes;
      using stan::math::value_of;
      using stan::math::include_summand;
      using std::log;

      // check if any vectors are zero length
      if (!(stan::length(y) && stan::length(sigma)))
        return 0.0;

      // set up return value accumulator
      T_partials_return logp(0.0);

      // validate args (here done over var, which should be OK)
      check_not_nan(function, "Random variable", y);
      check_positive(function, "Scale parameter", sigma);
      check_positive(function, "Random variable", y);
      check_consistent_sizes(function,
                             "Random variable", y,
                             "Scale parameter", sigma);

      // check if no variables are involved and prop-to
      if (!include_summand<propto, T_y, T_scale>::value)
        return 0.0;

      // set up template expressions wrapping scalars into vector views
      OperandsAndPartials<T_y, T_scale> operands_and_partials(y, sigma);

      VectorView<const T_y> y_vec(y);
      VectorView<const T_scale> sigma_vec(sigma);
      size_t N = max_size(y, sigma);

      VectorBuilder<true, T_partials_return, T_scale> inv_sigma(length(sigma));
      VectorBuilder<include_summand<propto, T_scale>::value,
                    T_partials_return, T_scale> log_sigma(length(sigma));
      for (size_t i = 0; i < length(sigma); i++) {
        inv_sigma[i] = 1.0 / value_of(sigma_vec[i]);
        if (include_summand<propto, T_scale>::value)
          log_sigma[i] = log(value_of(sigma_vec[i]));
      }

      for (size_t n = 0; n < N; n++) {
        // pull out values of arguments
        const T_partials_return y_dbl = value_of(y_vec[n]);

        // reusable subexpression values
        const T_partials_return y_over_sigma = y_dbl * inv_sigma[n];

        static double NEGATIVE_HALF = -0.5;

        // log probability
        if (include_summand<propto, T_scale>::value)
          logp -= 2.0 * log_sigma[n];
        if (include_summand<propto, T_y>::value)
          logp += log(y_dbl);
        // if (include_summand<propto, T_y, T_scale>::value)
        logp += NEGATIVE_HALF * y_over_sigma * y_over_sigma;

        // gradients
        T_partials_return scaled_diff = inv_sigma[n] * y_over_sigma;
        if (!is_constant_struct<T_y>::value)
          operands_and_partials.d_x1[n] += 1.0 / y_dbl - scaled_diff;
        if (!is_constant_struct<T_scale>::value)
          operands_and_partials.d_x2[n]
            += y_over_sigma * scaled_diff - 2.0 * inv_sigma[n];
      }
      return operands_and_partials.to_var(logp, y, sigma);
    }