// transform/fmllr-raw-test.cc

// Copyright  2009-2011 Microsoft Corporation
//            2013  Johns Hopkins University (author: Daniel Povey)

// See ../../COPYING for clarification regarding multiple authors
//
// 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
//
// THIS CODE IS PROVIDED *AS IS* BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
// KIND, EITHER EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION ANY IMPLIED
// WARRANTIES OR CONDITIONS OF TITLE, FITNESS FOR A PARTICULAR PURPOSE,
// MERCHANTABLITY OR NON-INFRINGEMENT.
// See the Apache 2 License for the specific language governing permissions and
// limitations under the License.

#include "util/common-utils.h"
#include "gmm/diag-gmm.h"
#include "transform/fmllr-diag-gmm.h"
#include "transform/fmllr-raw.h"

namespace kaldi {


void InitRandomGmm (DiagGmm *gmm_in) {
  int32 num_gauss = 5 + rand () % 4;
  int32 dim = 6 + Rand() % 5;
  DiagGmm &gmm(*gmm_in);
  gmm.Resize(num_gauss, dim);
  Matrix<BaseFloat> inv_vars(num_gauss, dim),
      means(num_gauss, dim);
  Vector<BaseFloat> weights(num_gauss);
  for (int32 i = 0; i < num_gauss; i++) {
    for (int32 j = 0; j < dim; j++) {
      inv_vars(i, j) = Exp(RandGauss() * (1.0 / (1 + j)));
      means(i, j) = RandGauss() * (1.0 / (1 + j));
    }
    weights(i) = Exp(RandGauss());
  }
  weights.Scale(1.0 / weights.Sum());
  gmm.SetWeights(weights);
  gmm.SetInvVarsAndMeans(inv_vars, means);
  gmm.ComputeGconsts();
}

void UnitTestFmllrRaw(bool use_offset) {
  using namespace kaldi;
  DiagGmm gmm;
  InitRandomGmm(&gmm);
  int32 model_dim =  gmm.Dim();

  int32 raw_dim = 5 + Rand() % 3;
  int32 num_splice = 1 + Rand() % 5;
  while (num_splice * raw_dim < model_dim) {
    num_splice++;
  }

  int32 full_dim = num_splice * raw_dim;
  int32 npoints = raw_dim*(raw_dim+1)*10;

  Matrix<BaseFloat> rand_points(npoints, full_dim);
  rand_points.SetRandn();

  Matrix<BaseFloat> lda_mllt(full_dim, full_dim + (use_offset ? 1 : 0)); // This is the full LDA+MLLT
  // matrix.  TODO: test with offset.
  lda_mllt.SetRandn();

  FmllrRawAccs accs(raw_dim, model_dim, lda_mllt);

  BaseFloat prev_objf_impr;
  for (int32 iter = 0; iter < 4; iter++) {

    for (int32 i = 0; i < npoints; i++) {
      SubVector<BaseFloat> sample(rand_points, i);
      accs.AccumulateForGmm(gmm, sample, 1.0);
    }

    Matrix<BaseFloat> fmllr_mat(raw_dim, raw_dim + 1);
    fmllr_mat.SetUnit(); // sets diagonal elements to one.

    FmllrRawOptions opts;
    BaseFloat objf_impr, count;
    accs.Update(opts, &fmllr_mat, &objf_impr, &count);

    KALDI_ASSERT(objf_impr > 0.0);

    if (iter != 0) {
      // This is not something provable, but is always true
      // in practice.
      KALDI_ASSERT(objf_impr < prev_objf_impr);
    }
    prev_objf_impr = objf_impr;


    // Now transform the raw features.
    for (int32 splice = 0; splice < num_splice; splice++) {
      SubMatrix<BaseFloat> raw_feats(rand_points,
                                     0, npoints,
                                     splice * raw_dim, raw_dim);
      for (int32 t = 0; t < npoints; t++) {
        SubVector<BaseFloat> this_feat(raw_feats, t);
        ApplyAffineTransform(fmllr_mat, &this_feat);
      }
    }
    accs.SetZero();
  }
}


}  // namespace kaldi ends here

int main() {
  kaldi::g_kaldi_verbose_level = 5;

  for (int i = 0; i < 2; i++) {  // did more iterations when first testing...
    kaldi::UnitTestFmllrRaw(i % 2 == 0);
  }
  std::cout << "Test OK.\n";
}
