void gmm_objective_b()

in src/cpp/modules/tapenade/gmm/gmm_b.c [375:474]

• 99 lines of code
• 12 McCabe index (conditional complexity)

void gmm_objective_b(int d, int k, int n, const double *alphas, double *
        alphasb, const double *means, double *meansb, const double *icf, 
        double *icfb, const double *x, Wishart wishart, double *err, double *
        errb) {
    int ix, ik;
    /* TFIX */
    const double CONSTANT = -n*d*0.5*log(2*PI);
    int icf_sz = d*(d+1)/2;
    double *Qdiags;
    double *Qdiagsb;
    double result1;
    double result1b;
    int ii1;
    Qdiagsb = (double *)malloc(d*k*sizeof(double));
    for (ii1 = 0; ii1 < d*k; ++ii1)
        Qdiagsb[ii1] = 0.0;
    Qdiags = (double *)malloc(d*k*sizeof(double));
    double *sum_qs;
    double *sum_qsb;
    sum_qsb = (double *)malloc(k*sizeof(double));
    for (ii1 = 0; ii1 < k; ++ii1)
        sum_qsb[ii1] = 0.0;
    sum_qs = (double *)malloc(k*sizeof(double));
    double *xcentered;
    double *xcenteredb;
    xcenteredb = (double *)malloc(d*sizeof(double));
    for (ii1 = 0; ii1 < d; ++ii1)
        xcenteredb[ii1] = 0.0;
    xcentered = (double *)malloc(d*sizeof(double));
    double *Qxcentered;
    double *Qxcenteredb;
    Qxcenteredb = (double *)malloc(d*sizeof(double));
    for (ii1 = 0; ii1 < d; ++ii1)
        Qxcenteredb[ii1] = 0.0;
    Qxcentered = (double *)malloc(d*sizeof(double));
    double *main_term;
    double *main_termb;
    main_termb = (double *)malloc(k*sizeof(double));
    for (ii1 = 0; ii1 < k; ++ii1)
        main_termb[ii1] = 0.0;
    main_term = (double *)malloc(k*sizeof(double));
    preprocess_qs_nodiff(d, k, icf, &(sum_qs[0]), &(Qdiags[0]));
    double slse = 0.;
    double slseb = 0.0;
    for (ix = 0; ix < n; ++ix)
        for (ik = 0; ik < k; ++ik) {
            pushReal8Array(xcentered, d); /* TFIX */
            subtract_nodiff(d, &(x[ix*d]), &(means[ik*d]), &(xcentered[0]));
            pushReal8Array(Qxcentered, d); /* TFIX */
            Qtimesx_nodiff(d, &(Qdiags[ik*d]), &(icf[ik*icf_sz + d]), &(
                           xcentered[0]), &(Qxcentered[0]));
            result1 = sqnorm_nodiff(d, &(Qxcentered[0]));
            pushReal8(main_term[ik]);
            main_term[ik] = alphas[ik] + sum_qs[ik] - 0.5*result1;
        }
    double lse_alphas;
    double lse_alphasb;
    slseb = *errb;
    lse_alphasb = -(n*(*errb));
    result1b = *errb;
    *errb = 0.0;
    log_wishart_prior_b(d, k, wishart, &(sum_qs[0]), &(sum_qsb[0]), &(Qdiags[0
                        ]), &(Qdiagsb[0]), icf, icfb, result1b);
    for (ii1 = 0; ii1 < k; ii1++) /* TFIX */
        alphasb[ii1] = 0.0;
    log_sum_exp_b(k, alphas, alphasb, lse_alphasb);
    for (ii1 = 0; ii1 < d * k; ii1++) /* TFIX */
        meansb[ii1] = 0.0;
    for (ix = n-1; ix > -1; --ix) {
        result1b = slseb;
        log_sum_exp_b(k, &(main_term[0]), &(main_termb[0]), result1b);
        for (ik = k-1; ik > -1; --ik) {
            popReal8(&(main_term[ik]));
            alphasb[ik] = alphasb[ik] + main_termb[ik];
            sum_qsb[ik] = sum_qsb[ik] + main_termb[ik];
            result1b = -(0.5*main_termb[ik]);
            main_termb[ik] = 0.0;
            sqnorm_b(d, &(Qxcentered[0]), &(Qxcenteredb[0]), result1b);
            popReal8Array(Qxcentered, d); /* TFIX */
            Qtimesx_b(d, &(Qdiags[ik*d]), &(Qdiagsb[ik*d]), &(icf[ik*icf_sz + 
                      d]), &(icfb[ik*icf_sz + d]), &(xcentered[0]), &(
                      xcenteredb[0]), &(Qxcentered[0]), &(Qxcenteredb[0]));
            popReal8Array(xcentered, d); /* TFIX */
            subtract_b(d, &(x[ix*d]), &(means[ik*d]), &(meansb[ik*d]), &(
                       xcentered[0]), &(xcenteredb[0]));
        }
    }
    preprocess_qs_b(d, k, icf, icfb, &(sum_qs[0]), &(sum_qsb[0]), &(Qdiags[0])
                    , &(Qdiagsb[0]));
    free(main_term);
    free(main_termb);
    free(Qxcentered);
    free(Qxcenteredb);
    free(xcentered);
    free(xcenteredb);
    free(sum_qs);
    free(sum_qsb);
    free(Qdiags);
    free(Qdiagsb);
}