text: // Copyright 2013 by Chris Dyer
  //
  // Licensed under the Apache License, Version 2.0 (the "License");
  // you may not use this file except in compliance with the License.
  // You may obtain a copy of the License at
  //
  //     http://www.apache.org/licenses/LICENSE-2.0
  //
  // Unless required by applicable law or agreed to in writing, software
  // distributed under the License is distributed on an "AS IS" BASIS,
  // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  // See the License for the specific language governing permissions and
  // limitations under the License.
  //
  
  #include <iostream>
  #include <cstdlib>
  #include <cmath>
  #include <utility>
  #include <fstream>
  #include <getopt.h>
  #include <sstream>
  
  #include "src/corpus.h"
  #include "src/ttables.h"
  #include "src/da.h"
  
  using namespace std;
  
  struct PairHash {
    size_t operator()(const pair<short, short>& x) const {
      return (unsigned short) x.first << 16 | (unsigned) x.second;
    }
  };
  
  Dict d; // integerization map
  
  void ParseLine(const string& line,
                 vector<unsigned>* src,
                 vector<unsigned>* trg) {
    static const unsigned kDIV = d.Convert("|||");
    vector<unsigned> tmp;
    src->clear();
    trg->clear();
    d.ConvertWhitespaceDelimitedLine(line, kDIV, &tmp);
    unsigned i = 0;
    while (i < tmp.size() && tmp[i] != kDIV) {
      src->push_back(tmp[i]);
      ++i;
    }
    if (i < tmp.size() && tmp[i] == kDIV) {
      ++i;
      for (; i < tmp.size(); ++i)
        trg->push_back(tmp[i]);
    }
  }
  
  string input;
  string conditional_probability_filename = "";
  string input_model_file = "";
  double mean_srclen_multiplier = 1.0;
  int is_reverse = 0;
  int ITERATIONS = 5;
  int favor_diagonal = 0;
  double beam_threshold = -4.0;
  double prob_align_null = 0.08;
  double diagonal_tension = 4.0;
  int optimize_tension = 0;
  int variational_bayes = 0;
  double alpha = 0.01;
  int no_null_word = 0;
  size_t thread_buffer_size = 10000;
  bool force_align = false;
  int print_scores = 0;
  struct option options[] = {
      {"input",             required_argument, 0,                  'i'},
      {"reverse",           no_argument,       &is_reverse,        1  },
      {"iterations",        required_argument, 0,                  'I'},
      {"favor_diagonal",    no_argument,       &favor_diagonal,    1  },
      {"force_align",       required_argument, 0,                  'f'},
      {"mean_srclen_multiplier", required_argument, 0,             'm'},
      {"beam_threshold",    required_argument, 0,                  't'},
      {"p0",                required_argument, 0,                  'q'},
      {"diagonal_tension",  required_argument, 0,                  'T'},
      {"optimize_tension",  no_argument,       &optimize_tension,  1  },
      {"variational_bayes", no_argument,       &variational_bayes, 1  },
      {"alpha",             required_argument, 0,                  'a'},
      {"no_null_word",      no_argument,       &no_null_word,      1  },
      {"conditional_probabilities", required_argument, 0,          'p'},
      {"thread_buffer_size", required_argument, 0,                 'b'},
      {0,0,0,0}
  };
  
  bool InitCommandLine(int argc, char** argv) {
    while (1) {
      int oi;
      int c = getopt_long(argc,
                          argv,
                          "i:rI:df:m:t:q:T:ova:Np:b:s",
                          options,
                          &oi);
      if (c == -1) break;
      cerr << "ARG=" << (char)c << endl;
      switch(c) {
        case 'i': input = optarg; break;
        case 'r': is_reverse = 1; break;
        case 'I': ITERATIONS = atoi(optarg); break;
        case 'd': favor_diagonal = 1; break;
        case 'f': force_align = 1; conditional_probability_filename = optarg; break;
        case 'm': mean_srclen_multiplier = atof(optarg); break;
        case 't': beam_threshold = atof(optarg); break;
        case 'q': prob_align_null = atof(optarg); break;
        case 'T': favor_diagonal = 1; diagonal_tension = atof(optarg); break;
        case 'o': optimize_tension = 1; break;
        case 'v': variational_bayes = 1; break;
        case 'a': alpha = atof(optarg); break;
        case 'N': no_null_word = 1; break;
        case 'p': conditional_probability_filename = optarg; break;
        case 'b': thread_buffer_size = atoi(optarg); break;
        case 's': print_scores = 1; break;
        default: return false;
      }
    }
    if (input.size() == 0) return false;
    return true;
  }
  
  void UpdateFromPairs(const vector<string>& lines, const int lc, const int iter,
      const bool final_iteration, const bool use_null, const unsigned kNULL,
      const double prob_align_not_null, double* c0, double* emp_feat,
      double* likelihood, TTable* s2t, vector<string>* outputs) {
    if (final_iteration) {
      outputs->clear();
      outputs->resize(lines.size());
    }
    double emp_feat_ = 0.0;
    double c0_ = 0.0;
    double likelihood_ = 0.0;
  #pragma omp parallel for schedule(dynamic) reduction(+:emp_feat_,c0_,likelihood_)
    for (int line_idx = 0; line_idx < static_cast<int>(lines.size());
        ++line_idx) {
      vector<unsigned> src, trg;
      ParseLine(lines[line_idx], &src, &trg);
      if (is_reverse)
        swap(src, trg);
      if (src.size() == 0 || trg.size() == 0) {
        cerr << "Error in line " << lc << "\n" << lines[line_idx] << endl;
        //return 1;
      }
      ostringstream oss; // collect output in last iteration
      vector<double> probs(src.size() + 1);
      bool first_al = true;  // used when printing alignments
      double local_likelihood = 0.0;
      for (unsigned j = 0; j < trg.size(); ++j) {
        const unsigned& f_j = trg[j];
        double sum = 0;
        double prob_a_i = 1.0 / (src.size() + use_null);  // uniform (model 1)
        if (use_null) {
          if (favor_diagonal)
            prob_a_i = prob_align_null;
          probs[0] = s2t->prob(kNULL, f_j) * prob_a_i;
          sum += probs[0];
        }
        double az = 0;
        if (favor_diagonal)
          az = DiagonalAlignment::ComputeZ(j + 1, trg.size(), src.size(),
              diagonal_tension) / prob_align_not_null;
        for (unsigned i = 1; i <= src.size(); ++i) {
          if (favor_diagonal)
            prob_a_i = DiagonalAlignment::UnnormalizedProb(j + 1, i, trg.size(),
                src.size(), diagonal_tension) / az;
          probs[i] = s2t->prob(src[i - 1], f_j) * prob_a_i;
          sum += probs[i];
        }
        if (final_iteration) {
          double max_p = -1;
          int max_index = -1;
          if (use_null) {
            max_index = 0;
            max_p = probs[0];
          }
          for (unsigned i = 1; i <= src.size(); ++i) {
            if (probs[i] > max_p) {
              max_index = i;
              max_p = probs[i];
            }
          }
          if (max_index > 0) {
            if (first_al)
              first_al = false;
            else
              oss << ' ';
            if (is_reverse)
              oss << j << '-' << (max_index - 1);
            else
              oss << (max_index - 1) << '-' << j;
          }
        } else {
          if (use_null) {
            double count = probs[0] / sum;
            c0_ += count;
            s2t->Increment(kNULL, f_j, count);
          }
          for (unsigned i = 1; i <= src.size(); ++i) {
            const double p = probs[i] / sum;
            s2t->Increment(src[i - 1], f_j, p);
            emp_feat_ += DiagonalAlignment::Feature(j, i, trg.size(), src.size()) * p;
          }
        }
        local_likelihood += log(sum);
      }
      likelihood_ += local_likelihood;
      if (final_iteration) {
        if (print_scores) {
          double log_prob = Md::log_poisson(trg.size(), 0.05 + src.size() * mean_srclen_multiplier);
          log_prob += local_likelihood;
          oss << " ||| " << log_prob;
        }
        oss << endl;
        (*outputs)[line_idx] = oss.str();
      }
    }
    *emp_feat += emp_feat_;
    *c0 += c0_;
    *likelihood += likelihood_;
  }
  
  inline void AddTranslationOptions(vector<vector<unsigned> >& insert_buffer,
      TTable* s2t) {
    s2t->SetMaxE(insert_buffer.size()-1);
  #pragma omp parallel for schedule(dynamic)
    for (unsigned e = 0; e < insert_buffer.size(); ++e) {
      for (unsigned f : insert_buffer[e]) {
        s2t->Insert(e, f);
      }
      insert_buffer[e].clear();
    }
  }
  
  void InitialPass(const unsigned kNULL, const bool use_null, TTable* s2t,
      double* n_target_tokens, double* tot_len_ratio,
      vector<pair<pair<short, short>, unsigned>>* size_counts) {
    ifstream in(input.c_str());
    if (!in) {
      cerr << "Can't read " << input << endl;
    }
    unordered_map<pair<short, short>, unsigned, PairHash> size_counts_;
    vector<vector<unsigned>> insert_buffer;
    size_t insert_buffer_items = 0;
    vector<unsigned> src, trg;
    string line;
    bool flag = false;
    int lc = 0;
    cerr << "INITIAL PASS " << endl;
    while (true) {
      getline(in, line);
      if (!in)
        break;
      lc++;
      if (lc % 1000 == 0) { cerr << '.'; flag = true; }
      if (lc %50000 == 0) { cerr << " [" << lc << "]\n" << flush; flag = false; }
      ParseLine(line, &src, &trg);
      if (is_reverse)
        swap(src, trg);
      if (src.size() == 0 || trg.size() == 0) {
        cerr << "Error in line " << lc << "\n" << line << endl;
      }
      *tot_len_ratio += static_cast<double>(trg.size()) / static_cast<double>(src.size());
      *n_target_tokens += trg.size();
      if (use_null) {
        for (const unsigned f : trg) {
          s2t->Insert(kNULL, f);
        }
      }
      for (const unsigned e : src) {
        if (e >= insert_buffer.size()) {
          insert_buffer.resize(e+1);
        }
        for (const unsigned f : trg) {
          insert_buffer[e].push_back(f);
        }
        insert_buffer_items += trg.size();
      }
      if (insert_buffer_items > thread_buffer_size * 100) {
        insert_buffer_items = 0;
        AddTranslationOptions(insert_buffer, s2t);
      }
      ++size_counts_[make_pair<short, short>(trg.size(), src.size())];
    }
    for (const auto& p : size_counts_) {
      size_counts->push_back(p);
    }
    AddTranslationOptions(insert_buffer, s2t);
  
    mean_srclen_multiplier = (*tot_len_ratio) / lc;
    if (flag) {
      cerr << endl;
    }
    cerr << "expected target length = source length * " << mean_srclen_multiplier << endl;
  }
  
  int main(int argc, char** argv) {
    if (!InitCommandLine(argc, argv)) {
      cerr << "Usage: " << argv[0] << " -i file.fr-en\n"
           << " Standard options ([USE] = strongly recommended):\n"
           << "  -i: [REQ] Input parallel corpus\n"
           << "  -v: [USE] Use Dirichlet prior on lexical translation distributions\n"
           << "  -d: [USE] Favor alignment points close to the monotonic diagonoal\n"
           << "  -o: [USE] Optimize how close to the diagonal alignment points should be\n"
           << "  -r: Run alignment in reverse (condition on target and predict source)\n"
           << "  -c: Output conditional probability table\n"
           << " Advanced options:\n"
           << "  -I: number of iterations in EM training (default = 5)\n"
           << "  -q: p_null parameter (default = 0.08)\n"
           << "  -N: No null word\n"
           << "  -a: alpha parameter for optional Dirichlet prior (default = 0.01)\n"
           << "  -T: starting lambda for diagonal distance parameter (default = 4)\n"
           << "  -s: print alignment scores (alignment ||| score, disabled by default)\n";
      return 1;
    }
    const bool use_null = !no_null_word;
    if (variational_bayes && alpha <= 0.0) {
      cerr << "--alpha must be > 0\n";
      return 1;
    }
    const double prob_align_not_null = 1.0 - prob_align_null;
    const unsigned kNULL = d.Convert("<eps>");
    TTable s2t, t2s;
    vector<pair<pair<short, short>, unsigned>> size_counts;
    double tot_len_ratio = 0;
    double n_target_tokens = 0;
  
    if (force_align) {
      ifstream in(conditional_probability_filename.c_str());
      s2t.DeserializeLogProbsFromText(&in, d);
      ITERATIONS = 0; // don't do any learning
    } else {
      InitialPass(kNULL, use_null, &s2t, &n_target_tokens, &tot_len_ratio, &size_counts);
      s2t.Freeze();
    }
  
    for (int iter = 0; iter < ITERATIONS; ++iter) {
      const bool final_iteration = (iter == (ITERATIONS - 1));
      cerr << "ITERATION " << (iter + 1) << (final_iteration ? " (FINAL)" : "") << endl;
      ifstream in(input.c_str());
      if (!in) {
        cerr << "Can't read " << input << endl;
        return 1;
      }
      double likelihood = 0;
      const double denom = n_target_tokens;
      int lc = 0;
      bool flag = false;
      string line;
      double c0 = 0;
      double emp_feat = 0;
      vector<string> buffer;
      vector<string> outputs;
      while(true) {
        getline(in, line);
        if (!in) break;
        ++lc;
        if (lc % 1000 == 0) { cerr << '.'; flag = true; }
        if (lc %50000 == 0) { cerr << " [" << lc << "]\n" << flush; flag = false; }
        buffer.push_back(line);
  
        if (buffer.size() >= thread_buffer_size) {
          UpdateFromPairs(buffer, lc, iter, final_iteration, use_null, kNULL,
              prob_align_not_null, &c0, &emp_feat, &likelihood, &s2t, &outputs);
          if (final_iteration) {
            for (const string& output : outputs) {
              cout << output;
            }
          }
          buffer.clear();
        }
      } // end data loop
      if (buffer.size() > 0) {
        UpdateFromPairs(buffer, lc, iter, final_iteration, use_null, kNULL,
            prob_align_not_null, &c0, &emp_feat, &likelihood, &s2t, &outputs);
        if (final_iteration) {
          for (const string& output : outputs) {
            cout << output;
          }
        }
        buffer.clear();
      }
  
      // log(e) = 1.0
      double base2_likelihood = likelihood / log(2);
  
      if (flag) {
        cerr << endl;
      }
      emp_feat /= n_target_tokens;
      cerr << "  log_e likelihood: " << likelihood << endl;
      cerr << "  log_2 likelihood: " << base2_likelihood << endl;
      cerr << "     cross entropy: " << (-base2_likelihood / denom) << endl;
      cerr << "        perplexity: " << pow(2.0, -base2_likelihood / denom) << endl;
      cerr << "      posterior p0: " << c0 / n_target_tokens << endl;
      cerr << " posterior al-feat: " << emp_feat << endl;
      //cerr << "     model tension: " << mod_feat / toks << endl;
      cerr << "       size counts: " << size_counts.size() << endl;
      if (!final_iteration) {
        if (favor_diagonal && optimize_tension && iter > 0) {
          for (int ii = 0; ii < 8; ++ii) {
            double mod_feat = 0;
  #pragma omp parallel for reduction(+:mod_feat)
            for(size_t i = 0; i < size_counts.size(); ++i) {
              const pair<short,short>& p = size_counts[i].first;
              for (short j = 1; j <= p.first; ++j)
                mod_feat += size_counts[i].second * DiagonalAlignment::ComputeDLogZ(j, p.first, p.second, diagonal_tension);
            }
            mod_feat /= n_target_tokens;
            cerr << "  " << ii + 1 << "  model al-feat: " << mod_feat << " (tension=" << diagonal_tension << ")\n";
            diagonal_tension += (emp_feat - mod_feat) * 20.0;
            if (diagonal_tension <= 0.1) diagonal_tension = 0.1;
            if (diagonal_tension > 14) diagonal_tension = 14;
          }
          cerr << "     final tension: " << diagonal_tension << endl;
        }
        if (variational_bayes)
          s2t.NormalizeVB(alpha);
        else
          s2t.Normalize();
      }
    }
    if (!force_align && !conditional_probability_filename.empty()) {
      cerr << "conditional probabilities: " << conditional_probability_filename << endl;
      s2t.ExportToFile(conditional_probability_filename.c_str(), d, beam_threshold);
    }
    if (force_align) {
      istream* pin = &cin;
      if (input != "-" && !input.empty())
        pin = new ifstream(input.c_str());
      istream& in = *pin;
      string line;
      vector<unsigned> src, trg;
      int lc = 0;
      double tlp = 0;
      while(getline(in, line)) {
        ++lc;
        ParseLine(line, &src, &trg);
        for (auto s : src) cout << d.Convert(s) << ' ';
        cout << "|||";
        for (auto t : trg) cout << ' ' << d.Convert(t);
        cout << " |||";
        if (is_reverse)
          swap(src, trg);
        if (src.size() == 0 || trg.size() == 0) {
          cerr << "Error in line " << lc << endl;
          return 1;
        }
        double log_prob = Md::log_poisson(trg.size(), 0.05 + src.size() * mean_srclen_multiplier);
  
        // compute likelihood
        for (unsigned j = 0; j < trg.size(); ++j) {
          unsigned f_j = trg[j];
          double sum = 0;
          int a_j = 0;
          double max_pat = 0;
          double prob_a_i = 1.0 / (src.size() + use_null);  // uniform (model 1)
          if (use_null) {
            if (favor_diagonal) prob_a_i = prob_align_null;
            max_pat = s2t.safe_prob(kNULL, f_j) * prob_a_i;
            sum += max_pat;
          }
          double az = 0;
          if (favor_diagonal)
            az = DiagonalAlignment::ComputeZ(j+1, trg.size(), src.size(), diagonal_tension) / prob_align_not_null;
          for (unsigned i = 1; i <= src.size(); ++i) {
            if (favor_diagonal)
              prob_a_i = DiagonalAlignment::UnnormalizedProb(j + 1, i, trg.size(), src.size(), diagonal_tension) / az;
            double pat = s2t.safe_prob(src[i-1], f_j) * prob_a_i;
            if (pat > max_pat) { max_pat = pat; a_j = i; }
            sum += pat;
          }
          log_prob += log(sum);
          if (true) {
            if (a_j > 0) {
              cout << ' ';
              if (is_reverse)
                cout << j << '-' << (a_j - 1);
              else
                cout << (a_j - 1) << '-' << j;
            }
          }
        }
        tlp += log_prob;
        cout << " ||| " << log_prob << endl << flush;
      } // loop over test set sentences
      cerr << "TOTAL LOG PROB " << tlp << endl;
    }
    return 0;
  }
  