/*******************************************************************************
 *
 *  subroutines for LDPC codes study, version 3.2.3
 *  Copyright (C) 2003-2011 Misha Stepanov
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 ******************************************************************************/

#ifdef LDPCC_RANDOM
  #include <gsl/gsl_rng.h>
  #include <gsl/gsl_randist.h>
#endif

#ifdef LDPCC_LP
  #include <glpk.h>
#endif

#define LDPCC_BIG_VALUE 1.e+10

#define ERROR(M) \
  { fprintf(stderr, M); exit(1); }

#define ALLOCATE(A, N, T, P, S) \
  { if (((A) = (P)malloc((N) * sizeof(T))) == NULL) \
     { fprintf(stderr, "malloc ERROR: %d %s\n", (int)N, S); exit(1); } }

typedef double real;

/*===type_structures==========================================================*/

typedef struct {
  int bits; /* length of the codeword (in bits) */
  int checks; /* number of parity checks */

  int **Hi, **Ha; /* parity check matrix arrays */
  int *punctured; /* 0 / 1 --- bit (is not) / (is) punctured */
  int *message; /* encoded message, the values are +/- 1 */
  real *h; /* log-likelihoods (magnetic fields) h_i */
  real *m; /* decoding output (magnetization) m_i */

  real **eta, **mu, **lambda, **old_eta, *pre_eta; /* messages */

  int HiHa_grown, ID_grown;
    /* 0 / 1 --- the corresponding arrays (are not) / (are) generated */
#ifdef LDPCC_LP
  LPX *lp; /* linear programming structure */
  int LP_grown; /* 0 / 1 --- is (is not) generated */
#endif
 } code;

/*----------------------------------------------------------------------------*/

#define Gaussian_channel 0
#define AWGN_channel 0
#define Laplacian_channel 1
#define BSC_channel 8
#define BEC_channel 9
#define isotropic_erasure_channel 19

typedef struct {
  real SNR; /* the ratio of *amplitudes* of signal and noise,
               for AWGN differs from standard one by SNR -> SNR^2 */
  real bit_damage_probability;

#ifdef LDPCC_RANDOM
  gsl_rng *RNG;
  int RNG_grown;
#endif

  int type; /* variants are Gaussian_channel, ..., see definitions above */
 } channel;

/*----------------------------------------------------------------------------*/

#define sum_product 0
#define min_sum 1
#define bit_flipping 2

/* M.G. Stepanov, M. Chertkov, Improving convergence of Belief
   Propagation decoding */
#define lambda_relaxation 1
#define naive_relaxation 2

typedef struct {
  int n_iter; /* number of iterations in iterative decoding */
/*  real (*node)(); */
  int WCC; /* != 0 --- with codeword check on each iteration */
  int relaxed; real tau; /* relaxation, slowing of iterations */
 } iterative_decoder;

/*===channel_noise============================================================*/

#ifdef LDPCC_RANDOM

void kill_RNG(channel *pC)
 {
  if (pC->RNG_grown)
   {
    gsl_rng_free(pC->RNG);
    pC->RNG_grown = 0;
   }
 }

void grow_RNG(channel *pC)
 {
  if (pC->RNG_grown) kill_RNG(pC);

  gsl_rng_env_setup();
  pC->RNG = gsl_rng_alloc(gsl_rng_ranlxd2);
  gsl_rng_set(pC->RNG, (unsigned long int)time(NULL));

  pC->RNG_grown = 1;
 }

void toss_noise(channel *pC, real *xi, const int bits)
 {
  int i, sign;

  for (i = 0; i < bits; i++)
    switch (pC->type)
     {
      case Gaussian_channel:
        xi[i] = gsl_ran_gaussian(pC->RNG, 1. / (pC->SNR));
        break;
      case Laplacian_channel:
        sign = 2 * gsl_rng_uniform_int(pC->RNG, 1) - 1;
        xi[i] = sign * gsl_ran_exponential(pC->RNG, 1. / (pC->SNR));
        break;
      case BSC_channel:
        xi[i] = 0.;
        if (gsl_rng_uniform(pC->RNG) < pC->bit_damage_probability)
          xi[i] = 2.;
        break;
      case BEC_channel:
        xi[i] = 0.;
        if (gsl_rng_uniform(pC->RNG) < pC->bit_damage_probability)
          xi[i] = 1.;
        break;
      case isotropic_erasure_channel:
        xi[i] = 1. + gsl_ran_gaussian(pC->RNG, 1. / (pC->SNR));
        break;
      default: /* Gaussian channel */
        xi[i] = gsl_ran_gaussian(pC->RNG, 1. / (pC->SNR));
        break;
     }
 }

#endif

/*===general_decoding_subroutines=============================================*/

void send_msg_w2b_calc_h(code *pH, channel *pC, real *xi, const int verbose)
/* initializing h_i --- sending messages from world to bits */
/*   h[i] = (1/2) log(p(x_i=+|x_out) / p(x_i=-|x_out)) */
 {
  int i;
  real SNR, SNR2, *h, amp_h;

  h = pH->h;
  SNR = pC->SNR; SNR2 = SNR * SNR;
  if (pC->type == 8)
   {
    amp_h = pC->bit_damage_probability;
    amp_h = (0.5 * log((1. - amp_h) / amp_h));
   }

  for (i = 0; i < pH->bits; i++)
    if ((pH->punctured)[i] == 0)
      switch (pC->type)
       {
        case Gaussian_channel:
          h[i] = ((real)((pH->message)[i]) - xi[i]) * SNR2; break;
        case Laplacian_channel:
          h[i] = ((real)((pH->message)[i]) - xi[i]) * SNR; break;
        case BSC_channel:
          if (xi[i] < 1.) h[i] = amp_h; else h[i] = -amp_h;
          h[i] *= ((real)((pH->message)[i]));
          break;
        case BEC_channel:
          if (xi[i] > 0.5) h[i] = 0.;
           else h[i] = (real)((pH->message)[i]) * LDPCC_BIG_VALUE;
          break;
        case isotropic_erasure_channel:
          h[i] = ((real)((pH->message)[i]) - xi[i]) * SNR; break;
          break;
        default: /* Gaussian channel */
          h[i] = ((real)((pH->message)[i]) - xi[i]) * SNR2; break;
       }
     else h[i] = 0.;

  for (i = 0; i < pH->bits; i++) (pH->m)[i] = h[i];
  if (verbose)
   {
    printf("# MSG_W->B\n0 ");
    for (i = 0; i < pH->bits; i++) printf("%e ", h[i]);
    printf("\n"); fflush(stdout);
   }
 }

int check_codeword(code *pH)
  /* checks if the decoding output m_i gives a codeword */
 {
  int a, I, flag, itmp;

  for (flag = 1, a = 0; a < (pH->checks); a++)
   {
    for (itmp = 1, I = 1; I <= (pH->Hi)[a][0]; I++)
     {
      if ((pH->m)[(pH->Hi)[a][I]] < 0.) itmp = -itmp;
      if ((pH->m)[(pH->Hi)[a][I]] == 0.) itmp = 0;
     }
    if (itmp <= 0) { flag = 0; break; }
   }

  return flag;
 }

int check_error(code *pH)
 {
  int i, flag;

  for (flag = 0, i = 0; i < (pH->bits); i++)
    if ((((real)(pH->message)[i]) * (pH->m)[i]) <= 0.)
      { flag = 1; break; }

  return flag;
 }

/*===iterative_decoding=======================================================*/

void kill_iterative_decoder(code *pH)
 {
  int i;

  if (pH->ID_grown)
   {
    if (pH->HiHa_grown == 0)
      ERROR("Iterative decoder is grown but not the code!\n")

    free(pH->pre_eta);
    for (i = 0; i < (pH->bits); i++)
      { free((pH->eta)[i]); free((pH->mu)[i]); free((pH->lambda)[i]); }
    free(pH->eta); free(pH->mu); free(pH->lambda);

    pH->ID_grown = 0;
   }
 }

void grow_iterative_decoder(code *pH)
 {
  int i;

  if (pH->HiHa_grown == 0)
   {
    if (pH->ID_grown)
      ERROR("Iterative decoder is grown but not the code!\n")
    else return;
   }
  kill_iterative_decoder(pH);

  ALLOCATE(pH->eta, pH->bits, real *, real **, "eta")
  for (i = 0; i < (pH->bits); i++)
    ALLOCATE((pH->eta)[i], (pH->Ha)[i][0] + 1, real, real *, "eta[]")

  ALLOCATE(pH->mu, pH->bits, real *, real **, "mu")
  for (i = 0; i < (pH->bits); i++)
    ALLOCATE((pH->mu)[i], (pH->Ha)[i][0] + 1, real, real *, "mu[]")

  ALLOCATE(pH->lambda, pH->bits, real *, real **, "lambda")
  for (i = 0; i < (pH->bits); i++)
    ALLOCATE((pH->lambda)[i], (pH->Ha)[i][0] + 1, real, real *, "lambda[]")

  ALLOCATE(pH->old_eta, pH->bits, real *, real **, "old_eta")
  for (i = 0; i < (pH->bits); i++)
    ALLOCATE((pH->old_eta)[i], (pH->Ha)[i][0] + 1, real, real *, "old_eta[]")

  ALLOCATE(pH->pre_eta, pH->bits, real, real *, "pre_eta")

  pH->ID_grown = 1;
 }

int iterative_decoding(code *pH, iterative_decoder *pD, const int verbose)
 {
  int it, bits, checks, i, j, a, J, A, B, sign;
  real *h, **eta, **mu, **lambda, **old_eta, *pre_eta, *m;
  real abs, min_abs;

  if (pH->ID_grown == 0)
    ERROR("Iterative decoder is not grown!\n")

/* making variables more local */
  bits = pH->bits; checks = pH->checks;
  h = pH->h; m = pH->m; eta = pH->eta; mu = pH->mu;
  lambda = pH->lambda; old_eta = pH->old_eta; pre_eta = pH->pre_eta;
  
/* setting eta_ia = mu_ia = lambda_ia = 0 */
  for (i = 0; i < bits; i++)
    for (A = 1; A <= (pH->Ha)[i][0]; A++) mu[i][A] = 0.;
  if (pD->relaxed)
    for (i = 0; i < bits; i++)
      for (A = 1; A <= (pH->Ha)[i][0]; A++)
        { eta[i][A] = 0.; lambda[i][A] = 0.; }

/* iterations, main cycle of the decoder */
  for (it = 0; it < (pD->n_iter); it++)
   {
    if (pD->WCC) if (check_codeword(pH)) break;

/* making backup of eta_ia */
    if (pD->relaxed)
      for (i = 0; i < bits; i++) for (A = 1; A <= (pH->Ha)[i][0]; A++)
        old_eta[i][A] = eta[i][A];

/* sending messages from bits to checks */
    for (i = 0; i < bits; i++) for (A = 1; A <= (pH->Ha)[i][0]; A++)
      for (eta[i][A] = h[i], B = 1; B <= (pH->Ha)[i][0]; B++)
        if (B != A) eta[i][A] += mu[i][B];

    if (pD->relaxed == lambda_relaxation)
     {
/* calculating lambda_ia from eta_ia, old_eta_ia */
      for (i = 0; i < bits; i++) for (A = 1; A <= (pH->Ha)[i][0]; A++)
        lambda[i][A] += ((pD->tau) * (eta[i][A] - old_eta[i][A]));
/* calculating eta_ia from lambda_ia */
      for (i = 0; i < bits; i++)
       {
        for (pre_eta[i] = 0., A = 1; A <= (pH->Ha)[i][0]; A++)
          pre_eta[i] += lambda[i][A];
        pre_eta[i] /= (real)((pH->Ha)[i][0] - 1.);
       }
      for (i = 0; i < bits; i++) for (A = 1; A <= (pH->Ha)[i][0]; A++)
        eta[i][A] = pre_eta[i] - lambda[i][A];
     }
    if (pD->relaxed == naive_relaxation)
      for (i = 0; i < bits; i++) for (A = 1; A <= (pH->Ha)[i][0]; A++)
        eta[i][A] = (pD->tau) * eta[i][A] + (1. - (pD->tau)) * old_eta[i][A];

    if (verbose > 1)
     {
      printf("### MSG_B->C\n");
      for (i = 0; i < bits; i++)
       {
        if (i != 0) printf("| ");
        for (A = 1; A <= (pH->Ha)[i][0]; A++) printf("%e ", eta[i][A]);
       }
      printf("\n"); fflush(stdout);
     }

/* sending messages from checks to bits */
    for (i = 0; i < bits; i++) for (A = 1; A <= (pH->Ha)[i][0]; A++)
      /* a = Ha[i][A] */
     {
      a = (pH->Ha)[i][A];
      sign = 1; min_abs = -1.;
      for (J = 1; J <= (pH->Hi)[a][0]; J++) if ((pH->Hi)[a][J] != i)
        /* j = Hi[a][J] */
       {
        j = (pH->Hi)[a][J];
        for (B = 1; B <= (pH->Ha)[j][0]; B++) if ((pH->Ha)[j][B] == a) break;
        abs = eta[j][B];
        if (abs < 0.) { sign = -sign; abs = -abs; }
        if ((abs < min_abs) || (min_abs < 0.)) min_abs = abs;
       }
      mu[i][A] = sign * min_abs;
     }

    if (verbose > 1)
     {
      printf("### MSG_C->B\n");
      for (i = 0; i < bits; i++)
       {
        if (i != 0) printf("| ");
        for (A = 1; A <= (pH->Ha)[i][0]; A++) printf("%e ", mu[i][A]);
       }
      printf("\n"); fflush(stdout);
     }

/* calculating m_i */
    for (i = 0; i < bits; i++)
      for (m[i] = h[i], A = 1; A <= (pH->Ha)[i][0]; A++) m[i] += mu[i][A];
    if (verbose)
     {
      printf("# MSG_B->W\n%d ", it + 1);
      for (i = 0; i < bits; i++) printf("%e ", m[i]);
      printf("\n"); fflush(stdout);
     }
   }

  return it;
 }

/*===LP_decoding==============================================================*/

#ifdef LDPCC_LP

void kill_lp_structure(code *pH)
 {
  if (pH->LP_grown)
   {
    if (pH->HiHa_grown == 0)
      ERROR("Linear programming structure is grown but not the code!\n")

    glp_delete_prob(pH->lp);

    pH->LP_grown = 0;
   }
 }

void grow_lp_decoder(code *pH)
 {
  int i, a, flag, rows, row, *ind;
  real rflag, *val;

  if (pH->HiHa_grown == 0)
   {
    if (pH->LP_grown)
      ERROR("Linear programming decoder is grown but not the code!\n")
    else return;
   }
  kill_lp_structure(pH);

  pH->lp = glp_create_prob();
  glp_set_obj_dir(pH->lp, GLP_MAX);

  glp_add_cols(pH->lp, pH->bits);
  for (i = 0; i < (pH->bits); i++)
    glp_set_col_bnds(pH->lp, i + 1, GLP_DB, -1., 1.);

  ALLOCATE(ind, (pH->bits) + 1, int, int *, "ind")
  ALLOCATE(val, (pH->bits) + 1, real, real *, "val")

  for (rows = 0, a = 0; a < (pH->checks); a++) if ((pH->Hi)[a][0] > 0) 
    rows += (1 << ((pH->Hi)[a][0] - 1));
  glp_add_rows(pH->lp, rows);

  for (row = 1, a = 0; a < (pH->checks); a++)
    switch ((pH->Hi)[a][0])
     {
      case 0: break;
      case 1: ind[1] = (pH->Hi)[a][1] + 1; val[1] = 1.;
              glp_set_row_bnds(pH->lp, row, GLP_FX, 1., 1.);
              glp_set_mat_row(pH->lp, row, 1, ind, val);
              row++; break;
      case 2: ind[1] = (pH->Hi)[a][1] + 1; val[1] = 1.;
              ind[2] = (pH->Hi)[a][2] + 1; val[2] = -1.;
              glp_set_row_bnds(pH->lp, row, GLP_FX, 0., 0.);
              glp_set_mat_row(pH->lp, row, 2, ind, val);
              row++; break;
      default: for (i = 1; i <= (pH->Hi)[a][0]; i++)
                 { ind[i] = (pH->Hi)[a][i] + 1; val[i] = 1.; }
               for (flag = 0; flag <= (pH->Hi)[a][0]; )
                {
                 for (rflag = 1., i = 1; i <= (pH->Hi)[a][0]; i++)
                   rflag *= val[i];

                 if (rflag < 0.) /* not a local codeword */
                  {
                   glp_set_row_bnds(pH->lp, row, GLP_UP,
                                    0., ((real)(pH->Hi)[a][0]) - 2.);
                   glp_set_mat_row(pH->lp, row, (pH->Hi)[a][0], ind, val);
                   row++;
                  }

                 for (flag = 1; flag <= (pH->Hi)[a][0]; flag++)
                  {
                   val[flag] *= -1.;
                   if (val[flag] < 0.) break;
                  }
                }
               break;
     }

  free(ind); free(val);

  pH->LP_grown = 1;
 }

void grow_lp_cone_cross_section(code *pH)
 {
  int i, a, rows, row, *ind;
  real *val;

  if (pH->HiHa_grown == 0)
   {
    if (pH->LP_grown)
      ERROR("[Small polytope] cone is grown but not the code!\n")
    else return;
   }
  kill_lp_structure(pH);

  pH->lp = glp_create_prob();
  glp_set_obj_dir(pH->lp, GLP_MAX);

  glp_add_cols(pH->lp, pH->bits);
  for (i = 0; i < (pH->bits); i++)
    glp_set_col_bnds(pH->lp, i + 1, GLP_DB, 0., 1.);

  ALLOCATE(ind, (pH->bits) + 1, int, int *, "ind")
  ALLOCATE(val, (pH->bits) + 1, real, real *, "val")

  for (rows = 1, a = 0; a < (pH->checks); a++) if ((pH->Hi)[a][0] > 0) 
    { if ((pH->Hi)[a][0] != 2) rows += (pH->Hi)[a][0]; else rows++; }
  glp_add_rows(pH->lp, rows);

  for (i = 1; i <= pH->bits; i++) { ind[i] = i; val[i] = 1.; }
  glp_set_row_bnds(pH->lp, 1, GLP_FX, 1., 1.);
  glp_set_mat_row(pH->lp, 1, pH->bits, ind, val);
  for (row = 2, a = 0; a < (pH->checks); a++)
    switch ((pH->Hi)[a][0])
     {
      case 0: break;
      case 1: ind[1] = (pH->Hi)[a][1] + 1; val[1] = 1.;
              glp_set_row_bnds(pH->lp, row, GLP_FX, 0., 0.);
              glp_set_mat_row(pH->lp, row, 1, ind, val);
              row++; break;
      case 2: ind[1] = (pH->Hi)[a][1] + 1; val[1] = 1.;
              ind[2] = (pH->Hi)[a][2] + 1; val[2] = -1.;
              glp_set_row_bnds(pH->lp, row, GLP_FX, 0., 0.);
              glp_set_mat_row(pH->lp, row, 2, ind, val);
              row++; break;
      default: for (i = 1; i <= (pH->Hi)[a][0]; i++)
                 { ind[i] = (pH->Hi)[a][i] + 1; val[i] = 1.; }
               for (i = 1; i <= (pH->Hi)[a][0]; i++)
                {
                 val[i] = -1.;
                 glp_set_row_bnds(pH->lp, row, GLP_LO, 0., 0.);
                 glp_set_mat_row(pH->lp, row, (pH->Hi)[a][0], ind, val);
                 row++; val[i] = 1.;
                }
               break;
     }

  free(ind); free(val);

  pH->LP_grown = 1;
 }

void lp_decoding(code *pH)
 {
  int i;

  if (pH->LP_grown == 0)
    { printf("Linear programming decoder is not grown!\n"); return; }

  for (i = 0; i < pH->bits; i++)
    glp_set_obj_coef(pH->lp, i + 1, (pH->h)[i]);

  glp_simplex(pH->lp, NULL);

  for (i = 0; i < pH->bits; i++)
    (pH->m)[i] = glp_get_col_prim(pH->lp, i + 1);
 }

#ifdef LDPCC_RANDOM

/* Finding low weight pseudocodewords using LP decoding and the
 * procedure described in
 *   M. Chertkov, M.G. Stepanov, An efficient pseudocodeword search
 *   algorithm for linear programming decoding of LDPC codes, IEEE
 *   Transactions on Information Theory 54 (4) 1514-1520 (2008).
 *     [arxiv: cs.IT/0601113] */

real lp_search(code *pH, channel *pC, real *xi)
 {
  int i, flag;
  real sxi, sxi2, old_sxi2;

  for (flag = 0, old_sxi2 = 1.e+10, sxi2 = 1.e+9; sxi2 != old_sxi2; )
   {
    if (flag == 0) toss_noise(pC, xi, pH->bits);
    send_msg_w2b_calc_h(pH, pC, xi, 0);
    lp_decoding(pH);

    for (flag = 0, i = 0; i < (pH->bits); i++)
      if ((((real)(pH->message)[i]) * (pH->m)[i]) <= 0.)
        { flag = 1; break; }

    if (flag)
     {
      old_sxi2 = sxi2;
      for (sxi = 0., sxi2 = 0., i = 0; i < (pH->bits); i++)
       {
        xi[i] = 1. - (pH->m)[i];
        sxi += fabs(xi[i]); sxi2 += (xi[i] * xi[i]);
       }

      for (i = 0; i < (pH->bits); i++) xi[i] *= (sxi / sxi2);
      for (sxi2 = 0., i = 0; i < (pH->bits); i++) sxi2 += (xi[i] * xi[i]);
     }
   }

  return sxi2;
 }

real lp_search2(code *pH, channel *pC, real *xi)
 {
  int i;
  real sxi, sxi2, old_sxi2;

  toss_noise(pC, xi, pH->bits);
  for (old_sxi2 = 1.e+10, sxi2 = 1.e+9; sxi2 != old_sxi2; )
   {
    for (i = 0; i < pH->bits; i++) (pH->h)[i] = xi[i];
    lp_decoding(pH);

    old_sxi2 = sxi2;
    for (sxi = 0., sxi2 = 0., i = 0; i < pH->bits; i++)
     {
      xi[i] = (pH->m)[i];
      sxi += fabs(xi[i]); sxi2 += (xi[i] * xi[i]);
     }

    for (i = 0; i < pH->bits; i++) xi[i] *= (sxi / sxi2);
    for (sxi2 = 0., i = 0; i < pH->bits; i++) sxi2 += (xi[i] * xi[i]);

/*    printf(" sxi = %22.16e\n", sxi);
    for (i = 0; i < pH->bits; i++) if (xi[i] < -1.e-10) printf("#");

    printf("sxi2 = %22.16e\n", sxi2); fflush(stdout); */
   }

  return sxi2;
 }

#endif

#endif

/*===parity_check_matrix======================================================*/

void kill_H_matrix(code *pH)
 {
  int i, a;

  if (pH->HiHa_grown)
   {
    if (pH->ID_grown) kill_iterative_decoder(pH);
#ifdef LDPCC_LP
    if (pH->LP_grown) kill_lp_structure(pH);
#endif

    free(pH->h); free(pH->m); free(pH->punctured); free(pH->message);
    for (i = 0; i < pH->bits; i++) free((pH->Ha)[i]); free(pH->Ha);
    for (a = 0; a < pH->checks; a++) free((pH->Hi)[a]); free(pH->Hi);

    pH->HiHa_grown = 0;
   }
 }

void read_H_matrix(const char *filename, code *pH)
 /*
  * The format of the file is the following:
  * --- begin ----------- | 3 2   |
  * #bits #checks         |       | <- example of the H-matrix file
  * matrix from 0s and 1s | 1 1 0 |    for the (3, 1) repetition code
  * --- end ------------- | 0 1 1 |
  * matrix form: #checks rows, #bits columns
  */
 {
  int i, a, c, *ni, *na;
  FILE *in;

  if (pH->ID_grown) kill_iterative_decoder(pH);
#ifdef LDPCC_LP
  if (pH->LP_grown) kill_lp_structure(pH);
#endif
  if (pH->HiHa_grown) kill_H_matrix(pH);

  in = fopen(filename, "r");
  fscanf(in, "%d %d", &i, &a); pH->bits = i; pH->checks = a;

  ALLOCATE(na, pH->bits, int, int *, "na")
  ALLOCATE(ni, pH->checks, int, int *, "ni")
  for (a = 0; a < (pH->checks); a++)
    for (i = 0;i < (pH->bits); i++)
      { fscanf(in, "%d", &c); if (c) { na[i]++; ni[a]++; } }
  fclose(in);

  ALLOCATE(pH->Hi, pH->checks, int *, int **, "Hi")
  for (a = 0; a < pH->checks; a++)
   {
    ALLOCATE((pH->Hi)[a], ni[a] + 1, int, int *, "Hi[]")
    (pH->Hi)[a][0] = ni[a];
   }
  ALLOCATE(pH->Ha, pH->bits, int *, int **, "Ha")
  for (i = 0; i < pH->bits; i++)
   {
    ALLOCATE((pH->Ha)[i], na[i] + 1, int, int *, "Ha[]")
    (pH->Ha)[i][0] = na[i];
   }

  for (i = 0; i< (pH->bits); i++) na[i] = 0;
  for (a = 0; a< (pH->checks); a++) ni[a] = 0;
  in = fopen(filename, "r");
  fscanf(in, "%d %d", &i, &a);
  for (a = 0; a < (pH->checks); a++) for (i = 0; i < (pH->bits); i++)
   {
    fscanf(in, "%d", &c);
    if (c)
      { na[i]++; (pH->Ha)[i][na[i]] = a; ni[a]++; (pH->Hi)[a][ni[a]] = i; }
   }
  fclose(in); free(na); free(ni);

  ALLOCATE(pH->punctured, pH->bits, int, int *, "punctured")
  for (i = 0; i < (pH->bits); i++) (pH->punctured)[i] = 0;

  ALLOCATE(pH->message, pH->bits, int, int *, "message")
  for (i = 0; i < (pH->bits); i++) (pH->message)[i] = 1;

  ALLOCATE(pH->h, pH->bits, real, real *, "h")
  ALLOCATE(pH->m, pH->bits, real, real *, "m")

  pH->HiHa_grown = 1;
 }

void read_short_H_matrix(const char *filename, code *pH)
 /*
  * The format of the file is the following:
  * --- begin ----------- | 3 2   |
  * #bits #checks         |       | <- example of the H-matrix file
  * #checks rows with the | 2 0 1 |    for the (3, 1) repetition code
  *   #bits connected and | 2 1 2 |
  *   their list          |
  * --- end ------------- |
  */
 {
  int i, a, I, *ni, *na;
  FILE *in;

  if (pH->ID_grown) kill_iterative_decoder(pH);
#ifdef LDPCC_LP
  if (pH->LP_grown) kill_lp_structure(pH);
#endif
  if (pH->HiHa_grown) kill_H_matrix(pH);

  in = fopen(filename, "r");
  fscanf(in, "%d %d", &i, &a); pH->bits = i; pH->checks = a;

  ALLOCATE(na, pH->bits, int, int *, "na")
  ALLOCATE(ni, pH->checks, int, int *, "ni")

  ALLOCATE(pH->Hi, pH->checks, int *, int **, "Hi")
  for (i = 0; i < (pH->bits); i++) na[i] = 0;
  for (a = 0; a < (pH->checks); a++)
   {
    fscanf(in, "%d", &(ni[a]));
    if (ni[a] > pH->bits)
      { fprintf(stderr, "ni[a] = %d\n", ni[a]); fflush(stderr); }
    ALLOCATE((pH->Hi)[a], ni[a] + 1, int, int *, "Hi[]")
    (pH->Hi)[a][0] = ni[a];

    for (I = 1; I <= (pH->Hi)[a][0]; I++)
      { fscanf(in, "%d", &i); (pH->Hi)[a][I] = i; na[i]++; }
   }
  fclose(in); free(ni);

  ALLOCATE(pH->Ha, pH->bits, int *, int **, "Ha")
  for (i = 0; i < pH->bits; i++)
   {
    ALLOCATE((pH->Ha)[i], na[i] + 1, int, int *, "Ha[]")
    (pH->Ha)[i][0] = na[i]; na[i] = 1;

    for (a = 0; a < (pH->checks); a++)
      for (I = 1; I <= (pH->Hi)[a][0]; I++)
        if ((pH->Hi)[a][I] == i) { (pH->Ha)[i][na[i]] = a; na[i]++; break; }
   }
  free(na);

  ALLOCATE(pH->punctured, pH->bits, int, int *, "punctured")
  for (i = 0; i < (pH->bits); i++) (pH->punctured)[i] = 0;

  ALLOCATE(pH->message, pH->bits, int, int *, "message")
  for (i = 0; i < (pH->bits); i++) (pH->message)[i] = 1;

  ALLOCATE(pH->h, pH->bits, real, real *, "h")
  ALLOCATE(pH->m, pH->bits, real, real *, "m")

  pH->HiHa_grown = 1;
 }

void write_H_matrix(const char *filename, code *pH)
 {
  int i, a, A, flag;
  FILE *out;

  out = fopen(filename, "w");
  fprintf(out, "%d %d\n\n", pH->bits, pH->checks);
  for (a = 0; a < (pH->checks); a++)
   {
    for (i = 0; i < (pH->bits); i++)
     {
      for (flag = 0, A = 1; A <= (pH->Ha)[i][0]; A++)
        if ((pH->Ha)[i][A] == a) { flag = 1; break; }
      fprintf(out, "%d ", flag);
     }
    fprintf(out, "\n");
   }
  fclose(out);
 }

void write_short_H_matrix(const char *filename, code *pH)
 {
  int a, I;
  FILE *out;

  out = fopen(filename, "w");
  fprintf(out, "%d %d\n\n", pH->bits, pH->checks);
  for (a = 0; a < (pH->checks); a++)
   {
    fprintf(out, "%d ", (pH->Hi)[a][0]);
    for (I = 1; I <= (pH->Hi)[a][0]; I++)
      fprintf(out, "%d ", (pH->Hi)[a][I]);
    fprintf(out, "\n");
   }
  fclose(out);
 }

/*===diluting_parity_check_matrix=============================================*/

/* "Dendro trick" for lowering connectivity of parity checks described in
 *   M. Chertkov, M. Stepanov, Pseudo-codeword Landscape,
 *   IEEE International Symposium on Information Theory 2007
 *   (24-29 June 2007, Nice, France), pp. 1546-1550.
 *     DOI: 10.1109/ISIT.2007.4557442
 *     [arxiv: cs.IT/0701084] */

int dilute_addition(int degree)
 {
  if (degree < 4) return 0;
   else
   {
    if (degree == 4) return 1;
     else
      return ((degree >> 1) + dilute_addition((degree + 1) >>  1));
   }
 }

void dilute_H_matrix(code *pH1, code *pH2, int number, int ind[])
 {
  int i, j, I, J, a, b, A, B, *ni, *na, *map, flag;

  kill_H_matrix(pH2);

  pH2->bits = pH1->bits; pH2->checks = pH1->checks;
  for (A = 0; A < number; A++)
   {
    j = dilute_addition((pH1->Hi)[ind[A]][0]);
    pH2->bits += j; pH2->checks += j;
   }

  ALLOCATE(na, pH2->bits, int, int *, "na")
  ALLOCATE(ni, pH2->checks, int, int *, "ni")
  ALLOCATE(pH2->Hi, pH2->checks, int *, int **, "Hi")
  ALLOCATE(pH2->Ha, pH2->bits, int *, int **, "Ha")

  for (i = 0; i < (pH1->bits); i++) na[i] = (pH1->Ha)[i][0];
  for (; i < (pH2->bits); i++) na[i] = 2;

/* moving all to be diluted checks to the back, the a'th (not being 
 * diluted) check goes to map[a]'th check in a new H matrix */
  ALLOCATE(map, pH1->bits, int, int *, "map")
  for (a = 0, b = 0; a < (pH1->checks); a++)
   {
    map[a] = -1;
    for (flag = 1, A = 0; A < number; A++)
      if (a == ind[A]) { flag = 0; break; }
    if (flag)
      { map[a] = b; ni[b] = (pH1->Hi)[a][0]; b++; }
   }
  for (; b < (pH2->checks); b++) ni[b] = 3;

  for (a = 0; a < pH2->checks; a++)
   {
    ALLOCATE((pH2->Hi)[a], ni[a] + 1, int, int *, "Hi[]")
    (pH2->Hi)[a][0] = ni[a];
   }
  for (i = 0; i < pH2->bits; i++)
   {
    ALLOCATE((pH2->Ha)[i], na[i] + 1, int, int *, "Ha[]")
    (pH2->Ha)[i][0] = na[i];
   }

/* untouched part of parity check matrix */
  for (i = 0; i < (pH1->bits); i++)
    for (A = 1; A <= na[i]; A++) (pH2->Ha)[i][A] = map[(pH1->Ha)[i][A]];
  for (; i < (pH2->bits); i++) { (pH2->Ha)[i][1] = (pH2->Ha)[i][2] = -1; }

  for (a = 0; a < (pH1->checks); a++) if (map[a] != -1)
    for (I = 1; I <= ni[map[a]]; I++) (pH2->Hi)[map[a]][I] = (pH1->Hi)[a][I];

/* diluted checks */
/* I / na is now number / list of bits we are to process */
/* J / B --- first available bit / check */
  for (J = (pH1->bits), B = (pH1->checks) - number, A = 0; A < number; A++)
   {
    I = (pH1->Hi)[ind[A]][0];
    for (i = 0; i < I; i++) na[i] = (pH1->Hi)[ind[A]][i + 1];
    while (I > 3)
     {
      if (I == 4)
       {
        for (j = 0; j < 4; j++)
         {
          for (b = 1; b <= (pH2->Ha)[na[j]][0]; b++)
            if ((pH2->Ha)[na[j]][b] == -1) break;
          (pH2->Ha)[na[j]][b] = B + (j >> 1);

          (pH2->Hi)[B + (j >> 1)][(j % 2) + 1] = na[j];
         }
        (pH2->Hi)[B][3] = J; (pH2->Hi)[B + 1][3] = J;
        (pH2->Ha)[J][1] = B; (pH2->Ha)[J][2] = B + 1;
        J++; B += 2; I = 0;
       }
       else
       {
        for (j = 0; j < (I >> 1); j++, J++)
         {
          for (b = 1; b <= (pH2->Ha)[na[2 * j]][0]; b++)
            if ((pH2->Ha)[na[2 * j]][b] == -1) break;
          (pH2->Ha)[na[2 * j]][b] = B + j;
          for (b = 1; b <= (pH2->Ha)[na[2 * j + 1]][0]; b++)
            if ((pH2->Ha)[na[2 * j + 1]][b] == -1) break;
          (pH2->Ha)[na[2 * j + 1]][b] = B + j;

          (pH2->Hi)[B + j][1] = na[2 * j];
          (pH2->Hi)[B + j][2] = na[2 * j + 1];

          (pH2->Hi)[B + j][3] = J; (pH2->Ha)[J][2] = B + j;

          na[j] = J;
         }
        if ((I % 2) == 1) na[j] = na[I - 1];
        B += (I >> 1); I = ((I + 1) >> 1);
        if ((I % 2) == 1)
         {
/* flipping the list of bits to make the diluting more symmetric */
          for (j = 0; j < (I >> 1); j++)
            { i = na[j]; na[j] = na[I - j - 1]; na[I - j - 1] = i; }
         }
       }
     }
    if (I == 3)
     {
      for (j = 0; j < 3; j++)
       {
        for (b = 1; b <= (pH2->Ha)[na[j]][0]; b++)
          if ((pH2->Ha)[na[j]][b] == -1) break;
        (pH2->Ha)[na[j]][b] = B;
        (pH2->Hi)[B][j + 1] = na[j];
       }
      B++;
     }
   }

  free(na); free(ni);

  ALLOCATE(pH2->punctured, pH2->bits, int, int *, "punctured")
  for (i = 0; i < (pH1->bits); i++) (pH2->punctured)[i] = (pH1->punctured)[i];
  for (; i < (pH2->bits); i++) (pH2->punctured)[i] = 1;

  ALLOCATE(pH2->message, pH2->bits, int, int *, "message")
  for (i = 0; i < (pH2->bits); i++) (pH2->message)[i] = 1;

  ALLOCATE(pH2->h, pH2->bits, real, real *, "h")
  ALLOCATE(pH2->m, pH2->bits, real, real *, "m")

  pH2->HiHa_grown = 1;
 }

/******************************************************************************/