开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 图形图像 > 图形图像处理 > 查看源码
pdfbthmt.cc
资源名称：contourletHMT.rar　[点击查看]
上传用户：l56789
上传日期：2022-02-25
资源大小：2422k
文件大小：41k
源码类别：

图形图像处理

开发平台：

Matlab

pdfbthmt.cc：源码内容
  1. //-----------------------------------------------------------------------------
  2. // pdfbthmt.cc  -  Tying Hidden Markov Tree Models for PDFB
  3. //-----------------------------------------------------------------------------
  6. #include <stdlib.h>
  7. #include <stdio.h>
  8. #include <math.h>
  9. #include <time.h>
  10. #include <iostream>
  11. #include "pdfbthmt.hh"
  12. #include "utils.hh"
  13. #include "mex.h"
  15. using std::cerr;
  16. using std::endl;
  18. //-----------------------------------------------------------------------------
  19. THMT::THMT(int ns, int nl, int* levndir, bool zm)
  20.     :M(ns), nCh(4), nLev(nl), zeromean(zm),
  21.      model_trans(nLev, vector<matrix<double> >() ),
  22.      model_mean(nLev, vector< vector<double> >() ), 
  23.      model_stdv(nLev, vector< vector<double> >() )
  24. {
  25.     int i, J, k;
  26.     
  27.     for(i=0; i<nLev; i++){
  28.       model_trans[i] = vector< matrix<double> >((int)pow(2.0,levndir[i]
  29.  -levndir[0]), 
  30. matrix<double>(M,M));
  31.       model_mean[i] = vector< vector<double> >((int)pow(2.0,levndir[i]
  32.  -levndir[0]), 
  33.        vector<double>(M));
  34.       model_stdv[i] = vector< vector<double> >((int)pow(2.0,levndir[i]
  35.  -levndir[0]), 
  36.        vector<double>(M));
  37.     }
  39.     rnd_init_model();
  41.     // subbandtree is a tree structure used to keep track of which
  42.     // directional subband does each coefficient in a tree belong to
  43.     subbandtree = tree<int> (1, nCh, nLev, 0);
  44.     // setting up the subbandtree tree
  45.     for(J = 0; J<nLev-1; J++)
  46.       for(i=0; i<subbandtree[J].size(); i++) 
  47. if(levndir[J+1]==levndir[J])
  48.   for(k=0; k<nCh; k++)
  49.     subbandtree[J+1][i*nCh+k] = subbandtree[J][i];
  50. else
  51.   for(k=0; k<nCh; k++)
  52.     if (k<nCh/2)
  53.       subbandtree[J+1][i*nCh+k] = subbandtree[J][i]*2;
  54.     else
  55.       subbandtree[J+1][i*nCh+k] = subbandtree[J][i]*2+1;
  56. }
  59. //-----------------------------------------------------------------------------
  60. THMT::THMT(const mxArray* initmodel, int* levndir)
  61. {
  62.     register int J, MM, m, mm, n, i, k;   
  63.     double* doubleptr;
  64.     mxArray* pointer, *pointer2, *pointer3;
  65.     int numfields, numels, zmlen;
  66.     char* zm;
  68.     numfields = mxGetNumberOfFields(initmodel);
  69.     numels = mxGetNumberOfElements(initmodel);
  70.     if ((numfields != 6) && (numfields != 7))
  71.       mexErrMsgTxt("ERROR: number of fields in struct model is incorrect");
  73.     if (numels != 1)
  74.       mexErrMsgTxt("ERROR: Too many elements");
  76.     if (strcmp(mxGetFieldNameByNumber(initmodel, 0), "nstates") != 0)
  77.       mexErrMsgTxt("Field 0 has wrong name");
  79.     pointer = mxGetFieldByNumber(initmodel, 0, 0);
  80.     doubleptr = mxGetPr(pointer);
  81.     M = (int)(*doubleptr);
  83.     nCh = 4;
  85.     if (strcmp(mxGetFieldNameByNumber(initmodel, 1), "nlevels") != 0)
  86.       mexErrMsgTxt("Field 1 has wrong name");    
  88.     pointer = mxGetFieldByNumber(initmodel, 0, 1);
  89.     doubleptr = mxGetPr(pointer);
  90.     nLev = (int)*doubleptr;
  92.     if (strcmp(mxGetFieldNameByNumber(initmodel, 2), "zeromean") != 0)
  93.       mexErrMsgTxt("Field 2 has wrong name");
  94.     pointer = mxGetFieldByNumber(initmodel, 0, 2);
  95.     zmlen = mxGetM(pointer) * mxGetN(pointer) + 1; 
  96.     zm = (char*)mxCalloc(zmlen, sizeof(char));
  97.     if (mxGetString(pointer, zm, zmlen) != 0)
  98.       mexErrMsgTxt("Error: cannot get zeromean variable");
  99.     if (strcmp(zm, "yes") == 0)
  100.       zeromean = true;
  101.     else if (strcmp(zm, "no") == 0)
  102.       zeromean = false;
  103.     else
  104.       mexErrMsgTxt("Error: Zeromean can only be yes or no");
  106.     model_trans = vector<vector<matrix<double> > >(nLev, 
  107.   vector<matrix<double> >() );
  108.     model_mean = vector<vector<vector<double> > >(nLev, 
  109.  vector< vector<double> >() );
  110.     model_stdv = vector<vector<vector<double> > >(nLev, 
  111.  vector< vector<double> >() );
  112.     for(i=0; i<nLev; i++){
  113.       model_trans[i] = vector< matrix<double> >((int)pow(2.0,levndir[i]
  114.  -levndir[0]), 
  115. matrix<double>(M,M));
  116.       model_mean[i] = vector< vector<double> >((int)pow(2.0,levndir[i]
  117.  -levndir[0]), 
  118.        vector<double>(M));
  119.       model_stdv[i] = vector< vector<double> >((int)pow(2.0,levndir[i]
  120.  -levndir[0]), 
  121.        vector<double>(M));
  122.     }
  124.     if (strcmp(mxGetFieldNameByNumber(initmodel, 3), "rootprob") != 0)
  125.       mexErrMsgTxt("Field 3 has wrong name");
  126.     pointer = mxGetFieldByNumber(initmodel, 0, 3);
  127.     if (pointer == NULL) {
  128.       mexPrintf("%s%dn",
  129.      "FIELD:", 3);
  130.       mexErrMsgTxt("Above field is empty!"); 
  131.     }
  132.     doubleptr = mxGetPr(pointer);
  133.     for (m = 0; m < M; m++)
  134.          model_trans[0][0][m][0]= doubleptr[m];
  135.  
  136.     if (strcmp(mxGetFieldNameByNumber(initmodel, 4), "transprob") != 0)
  137.       mexErrMsgTxt("Field 4 has wrong name");
  138.     pointer = mxGetFieldByNumber(initmodel, 0, 4);
  139.     if (pointer == NULL) {
  140.       mexPrintf("%s%dn",
  141.      "FIELD:", 4);
  142.       mexErrMsgTxt("Above field is empty!"); 
  143.     }
  145.     // Trans probs
  146.     for (J =  1; J < nLev; J++)
  147.     { 
  148.       pointer2 = mxGetCell(pointer,J-1);
  149.       for (m = 0; m < M; m++){
  150. for (mm = 0; mm < model_trans[J].size(); mm++) {
  151.     pointer3 = mxGetCell(pointer2, mm);
  152.     doubleptr = mxGetPr(pointer3);
  153.     for (n = 0; n < M; n++)
  154.        model_trans[J][mm][m][n] = doubleptr[n*M+m];
  155. }
  156.       }
  157.     }
  159.     // Mean
  160.     if (!zeromean)
  161.     {
  162.         if (strcmp(mxGetFieldNameByNumber(initmodel, 5), "mean") != 0)
  163.   mexErrMsgTxt("Field 5 has wrong name");
  164.         pointer = mxGetFieldByNumber(initmodel, 0, 5);
  165. if (pointer == NULL) {
  166.   mexPrintf("%s%dn",
  167.     "FIELD:", 5);
  168.   mexErrMsgTxt("Above field is empty!");
  169. }
  170. for (J = 0; J < nLev; J++)
  171. {
  172.   pointer2 = mxGetCell(pointer, J);
  173.   for (mm=0; mm < model_mean[J].size(); mm++){
  174.     pointer3 = mxGetCell(pointer2, mm);
  175.     doubleptr = mxGetPr(pointer3);
  176.     for (m=0; m<M; m++)
  177.        model_mean[J][mm][m] = doubleptr[m];
  178.   }
  179. }
  180.     }
  182.     if(!zeromean){
  183.       if (strcmp(mxGetFieldNameByNumber(initmodel, 6), "stdv") != 0)
  184. mexErrMsgTxt("Field 6 has wrong name");
  185.       pointer = mxGetFieldByNumber(initmodel, 0, 6);
  186.       if (pointer == NULL) {
  187. mexPrintf("%s%dn",
  188.   "FIELD:", 6);
  189. mexErrMsgTxt("Above field is empty!");
  190.       }
  191.     }
  192.     else{
  193.       if (strcmp(mxGetFieldNameByNumber(initmodel, 5), "stdv") != 0)
  194. mexErrMsgTxt("Field 5 has wrong name");
  195.       pointer = mxGetFieldByNumber(initmodel, 0, 5);
  196.       if (pointer == NULL) {
  197. mexPrintf("%s%dn",
  198.   "FIELD:", 5);
  199. mexErrMsgTxt("Above field is empty!");
  200.       } 
  201.     }
  203.     // Standard deviation
  204.     for (J =  0; J < nLev; J++)
  205.     {
  206.       pointer2 = mxGetCell(pointer, J);
  207.       for (mm=0; mm < model_stdv[J].size(); mm++){
  208. pointer3 = mxGetCell(pointer2, mm);
  209. doubleptr = mxGetPr(pointer3);
  210. for (m = 0; m < M; m++) 
  211.     model_stdv[J][mm][m] = doubleptr[m];
  212.       }
  213.     }
  215.     subbandtree = tree<int> (1, nCh, nLev, 0);
  216.     // setting up the subbandtree tree
  217.     for(J = 0; J<nLev-1; J++)
  218.       for(i=0; i<subbandtree[J].size(); i++) 
  219. if(levndir[J+1]==levndir[J])
  220.   for(k=0; k<nCh; k++)
  221.     subbandtree[J+1][i*nCh+k] = subbandtree[J][i];
  222. else
  223.   for(k=0; k<nCh; k++)
  224.     if (k<nCh/2)
  225.       subbandtree[J+1][i*nCh+k] = subbandtree[J][i]*2;
  226.     else
  227.       subbandtree[J+1][i*nCh+k] = subbandtree[J][i]*2+1;
  228. }
  231. //-----------------------------------------------------------------------------
  232. THMT::THMT(char *filename, int* levndir, long& startpos)
  233. {
  234.     FILE *fp;
  235.     register int J, MM, m, mm, n, i, k;
  236.     char szm[10];
  237.    
  238.     fp = fopen(filename,"a+");
  239.     if (!fp) 
  240.     {
  241. cerr << "ERROR: can not open for reading " << filename << endl;
  242. exit(-1);
  243.     }
  245.     fseek(fp, startpos, SEEK_SET);
  247.     if (fscanf(fp, "nStates: %dn", &M) != 1) 
  248.     {
  249. cerr << "ERROR: problem reading nStates field of " 
  250.      << filename << endl;
  251. exit(-1);
  252.     }
  254.     nCh = 4;
  256.     if (fscanf(fp, "nLevels: %dn", &nLev) != 1) 
  257.     {
  258. cerr << "ERROR: problem reading nLevels field of " 
  259.      << filename << endl;
  260. exit(-1);
  261.     }
  263.     if (fscanf(fp, "zeroMean: %sn", &szm) != 1) 
  264.     {
  265. cerr << "ERROR: problem reading zeroMean field of " 
  266.      << filename << endl;
  267. exit(-1);
  268.     }
  270.     if (strcmp(szm, "yes") == 0)
  271. zeromean = true;
  272.     else
  273. zeromean = false;
  275.     // Allocate space for model parameters
  276.     model_trans = vector<vector<matrix<double> > >(nLev, 
  277.   vector<matrix<double> >() );
  278.     model_mean = vector<vector<vector<double> > >(nLev, 
  279.  vector< vector<double> >() );
  280.     model_stdv = vector<vector<vector<double> > >(nLev, 
  281.  vector< vector<double> >() );
  282.     for(i=0; i<nLev; i++){
  283.       model_trans[i] = vector< matrix<double> >((int)pow(2.0,levndir[i]
  284.  -levndir[0]), 
  285. matrix<double>(M,M));
  286.       model_mean[i] = vector< vector<double> >((int)pow(2.0,levndir[i]
  287.  -levndir[0]), 
  288.        vector<double>(M));
  289.       model_stdv[i] = vector< vector<double> >((int)pow(2.0,levndir[i]
  290.  -levndir[0]), 
  291.        vector<double>(M));
  292.     }
  294.     // Initial probs
  295.     fscanf(fp, "n");
  296.     for (m = 0; m < M; m++)
  297. fscanf(fp, "%lf ", &model_trans[0][0][m][0]);  
  298.     fscanf(fp, "nn");
  299.   
  300.     // Trans prob
  301.     for (J = 1; J < nLev; J++) 
  302.     {
  303.       for (m = 0; m < M; m++){
  304. for (mm = 0; mm < model_trans[J].size(); mm++)
  305.   {
  306.     for (n = 0; n < M; n++)
  307.       fscanf(fp, "%lf ", &model_trans[J][mm][m][n]);
  308.   }
  310. fscanf(fp, "n");
  311.       }
  312.     }
  313.     fscanf(fp, "n");
  315.     // Mean
  316.     if (zeromean)
  317.     {
  318.         for (J = 0; J < nLev; J++){
  319.   for (mm=0; mm<model_mean[J].size(); mm++)
  320.     for (m = 0; m < M; m++) 
  321. model_mean[J][mm][m] = 0.0;
  322. }
  323.     }
  324.     else
  325.     {
  326. for (J = 0; J < nLev; J++)
  327. {
  328.   for (mm=0; mm<model_mean[J].size(); mm++)
  329.     for (m = 0; m < M; m++) 
  330. fscanf(fp, "%lf ", &model_mean[J][mm][m]);
  331.     
  332.   fscanf(fp, "n");
  333. }
  334. fscanf(fp, "n");
  335.     }
  337.     // Standard deviation
  338.     for (J=0; J < nLev; J++)
  339.     {
  340.       for (mm=0; mm<model_stdv[J].size(); mm++)
  341. for (m = 0; m < M; m++) 
  342.     fscanf(fp, "%lf ", &model_stdv[J][mm][m]);
  343.   
  344.       fscanf(fp, "n");
  345.     }
  346.     fscanf(fp, "n");
  348.     startpos = ftell(fp);
  349.    
  350.     fclose(fp);
  352.     subbandtree = tree<int> (1, nCh, nLev, 0);
  353.     // setting up the subbandtree tree
  354.     for(J = 1; J<nLev-1; J++)
  355.       for(i=0; i<subbandtree[J].size(); i++) 
  356. if(J%2==0)
  357.   for(k=0; k<nCh; k++)
  358.     subbandtree[J+1][i*nCh+k] = subbandtree[J][i];
  359. else
  360.   for(k=0; k<nCh; k++)
  361.     if (k<nCh/2)
  362.       subbandtree[J+1][i*nCh+k] = subbandtree[J][i]*2;
  363.     else
  364.       subbandtree[J+1][i*nCh+k] = subbandtree[J][i]*2+1;
  365. }
  368. //-----------------------------------------------------------------------------
  369. void THMT::allocate_training()
  370. {
  371.     int i, J, k;
  373.     if (alpha.nlev() == 0) // avoid re-allocate
  374.       {
  375. alpha = 
  376.     tree< vector<double> > (1, nCh, nLev, vector<double>(M));
  377. beta = beta_par = 
  378.     tree< vector<double> > (1, nCh, nLev, vector<double>(M));
  379. state_prob = tree< vector<double> > (nObs, nCh, nLev, 
  380.      vector<double>(M));
  381. scaling = tree<double> (1, nCh, nLev);
  382. sum_prob = sum_mean = sum_stdv = 
  383.     vector< vector< vector<double> > >(nLev, vector< vector<double> >());
  385. sum_trans = 
  386.   vector< vector< matrix<double> > > 
  387.   (nLev, vector<matrix<double> >());
  389. for(i=0; i<nLev; i++) {
  390.     sum_prob[i] = vector< vector<double> >(model_stdv[i].size(), 
  391.    vector<double>(M));
  392.     sum_mean[i] = vector< vector<double> >(model_stdv[i].size(), 
  393.    vector<double>(M));
  394.     sum_stdv[i] = vector< vector<double> >(model_stdv[i].size(), 
  395.    vector<double>(M));
  396.     sum_trans[i] = vector< matrix<double> >(model_stdv[i].size(), 
  397.     matrix<double>(M, M)); 
  398. }
  399.     }    
  400. }  
  403. //-----------------------------------------------------------------------------
  404. void THMT::allocate_testing()
  405. {
  406.     int i, J, k;
  408.     if (beta.nlev() == 0) // avoid re-allocate
  409.     {
  410. beta = beta_par = 
  411.     tree< vector<double> > (1, nCh, nLev, vector<double>(M));
  412. scaling = tree<double> (1, nCh, nLev);
  413.     }
  414. }
  417. //-----------------------------------------------------------------------------
  418. void THMT::compute_beta(int ob_ind)
  419. {
  420.     register int J, nNode, i, m, mm, n, c;
  421.     register double o;
  422.     register double sum, prod;
  423.   
  425.     // Initialization of beta down-tree
  426.     J = nLev - 1; // finest scale
  428.     for (i = 0, nNode = beta[J].size(); i < nNode; i++)
  429.     {
  430. // o_{J,i}^{ob_ind}
  431. o = (*obs)[J][ob_ind * ipow(nCh,J) + i];
  433. mm = subbandtree[J][i];
  435. for (m = 0; m < M; m++) 
  436.   beta[J][i][m] = compute_g(o, model_mean[J][mm][m], 
  437.       model_stdv[J][mm][m]);
  438.     } // end initialization
  440.     // Rescale beta
  441.     rescale(J);
  443.     while (J > 0) 
  444.     {
  445. // Compute $beta_{i,p(i)}$
  446. for (i = 0, nNode = beta_par[J].size(); i < nNode; i++) 
  447. {   
  448.   mm = subbandtree[J][i];
  449.   for (m = 0; m < M; m++) 
  450.     {
  451.       sum = 0.0;
  453.       for (n = 0; n < M; n++)
  454. sum += model_trans[J][mm][n][m] * beta[J][i][n];
  456.       beta_par[J][i][m] = sum;
  457.     }
  458. }
  460. // Compute $beta_{p(i)}(m)$
  461. J--;
  462. for (i = 0, nNode = beta[J].size(); i < nNode; i++) 
  463. {
  464.     // o_{J,i}^{ob_ind}
  465.     o = (*obs)[J][ob_ind * ipow(nCh,J) + i];
  466.     mm = subbandtree[J][i];
  468.     for (m = 0; m < M; m++) 
  469.     {
  470. prod = compute_g(o, model_mean[J][mm][m], 
  471.  model_stdv[J][mm][m]);
  473. for (c = 0; c < nCh; c++)  // look for child of $p(i)$
  474.     prod *= beta_par[J+1][i*nCh+c][m];
  476. beta[J][i][m] = prod;
  477.     }
  478. }
  480. // Rescale $beta_{p(i)}(m)$
  481. rescale(J);
  482.     }
  483. }
  486. //-----------------------------------------------------------------------------
  487. void THMT::rescale(int J)
  488. {
  489.     register int i, m, nNode = beta[J].size();
  490.     register double sum;
  492.     // Rescale for each node in this level (J)
  493.     for (i = 0; i < nNode; i++) 
  494.     {
  495. sum = 0.0;
  497. for (m = 0; m < M; m++) 
  498.     sum += beta[J][i][m];
  499.     
  500. scaling[J][i] = sum;
  502. for (m = 0; m < M; m++) 
  503.     beta[J][i][m] /= sum;
  504.     }
  505. }
  507. //-----------------------------------------------------------------------------
  508. void THMT::compute_alpha()
  509. {
  510.     register int J, i, m, mm, n, nNode;
  511.     register double sum;
  513.     // Initialize the coarsest level
  514.     for (m = 0; m < M; m++)
  515. alpha[0][0][m] = model_trans[0][0][m][0];
  517.     for (J = 1; J < nLev; J++) 
  518.     {
  519. for (i = 0, nNode = alpha[J].size(); i < nNode; i++) 
  520. {
  521.   mm = subbandtree[J][i];
  522.   for (m = 0; m < M; m++) 
  523.     {
  524.       sum = 0.0;
  526.       for (n = 0; n < M; n++) 
  527. sum += model_trans[J][mm][m][n] *
  528.   alpha[J-1][i/nCh][n] *
  529.   beta[J-1][i/nCh][n] /
  530.   beta_par[J][i][n];
  532. alpha[J][i][m] = sum;
  533.     } // end for m
  534. } // end for i
  535.     } // end for J
  536. }
  539. //-----------------------------------------------------------------------------
  540. double THMT::compute_likelihood()
  541. {
  542.     register int J, i, nNode, m;
  543.     register double f;
  544.     register double log_scale;
  546.     for (m = 0, f = 0.0; m < M; m++)
  547. f += beta[0][0][m] * model_trans[0][0][m][0];
  549.     // Re-scale back using saved scaling factors
  550.     for (J = 0, log_scale = 0.0; J < nLev; J++)
  551. for (i = 0, nNode = scaling[J].size(); i < nNode; i++) 
  552.     log_scale += log(scaling[J][i]);
  554.     // Final result
  555.     return (log_scale + log(f));
  556. }
  559. //-----------------------------------------------------------------------------
  560. void THMT::update_model(double& delta, double& avf)
  561. {
  562.     register int J, MM, i, m, mm, n, nNode;
  563.     register int ob_ind;           // observation index
  564.     register double o, denom, prob1, prob2, newval;
  566.     // Reinitialize
  567.     delta = avf = 0.0;
  569.     for (J = 0; J < nLev; J++)
  570.       for (mm=0; mm<sum_prob[J].size(); mm++)
  571. for (m = 0; m < M; m++) 
  572. {
  573.     sum_prob[J][mm][m] = 0.0;
  574.     sum_mean[J][mm][m] = 0.0;
  575.     sum_stdv[J][mm][m] = 0.0;
  577.     for (n = 0; n < M; n++)
  578.       sum_trans[J][mm][m][n] = 0.0;
  579. }
  580.   
  581.     // For each training tree
  582.     for (ob_ind = 0; ob_ind < nObs; ob_ind++) 
  583.     {
  584. // Compute $alpha, beta$
  585. compute_beta(ob_ind);
  586. compute_alpha();
  588. // Denominator for restimated probs
  589. denom = 0.0;
  590. for (m = 0; m < M; m++)
  591.     denom += beta[0][0][m] * alpha[0][0][m];
  593. // Compute state probabilities for denoising
  594. for (J = 0; J < nLev; J++)
  595.     for (i = 0, nNode = alpha[J].size(); i < nNode; i++)
  596.       for (m = 0; m < M; m++)
  597. {
  598.   state_prob[J][ob_ind*ipow(nCh,J)+i][m] = 
  599.     alpha[J][i][m]*beta[J][i][m]/denom;
  600. }
  602. // Update total log-likelihood
  603. avf += compute_likelihood();
  605. for (J = 0; J < nLev; J++)
  606.             for (i = 0, nNode = alpha[J].size(); i < nNode; i++) 
  607.            {
  608.         // o_{J,i}^{ob_ind}
  609.         o = (*obs)[J][ob_ind * ipow(nCh,J) + i];
  611.         for (m = 0; m < M; m++) 
  612.      {
  613.     // Compute $prob(S_i | O)$
  614.             prob1 = alpha[J][i][m] * beta[J][i][m]
  615.         / denom;
  617.     mm = subbandtree[J][i];
  618.     // Summing for all trees
  619.             sum_prob[J][mm][m] += prob1;
  620.             sum_mean[J][mm][m] += o * prob1;
  621.             sum_stdv[J][mm][m] += (o - model_mean[J][mm][m]) *
  622.         (o - model_mean[J][mm][m]) *
  623.         prob1;
  624.         }
  625.            }
  627. for (J = 1; J < nLev; J++)
  628.     for (i = 0, nNode = alpha[J].size(); i < nNode; i++)
  629. for (m = 0; m < M; m++)
  630.     for (n = 0; n < M; n++) 
  631.     {
  632.         mm = subbandtree[J][i];
  633. // Compute $prob(S_i, S_{p(i) | O)$
  634. prob2 = beta[J][i][m] *
  635.   model_trans[J][mm][m][n] *
  636.   alpha[J-1][i/nCh][n] *
  637.   beta[J-1][i/nCh][n] /
  638.   beta_par[J][i][n]
  639.   / denom;
  640. sum_trans[J][mm][m][n] += prob2;
  641.     }
  643.     } // end ob_ind
  645.     // Average log-likelihood
  646.     avf = avf / double(nObs);
  648.     // Normalize and update model parameters
  649.     for (J = 0; J < nLev; J++)
  650.       for (mm = 0; mm < model_stdv[J].size(); mm++)
  651. for (m = 0; m < M; m++) 
  652. {
  653.     if (!zeromean) // Only update means for non-zeromean model
  654.     {
  655. newval = sum_mean[J][mm][m] / sum_prob[J][mm][m];
  656. delta += fabs(newval - model_mean[J][mm][m]);
  657. model_mean[J][mm][m] = newval;
  658.     }
  660.     newval = sqrt(sum_stdv[J][mm][m] / sum_prob[J][mm][m]);
  661.     delta += fabs(newval - model_stdv[J][mm][m]);
  662.     model_stdv[J][mm][m] = newval;
  663. }
  665.     //state probs 
  666.     for (m = 0; m < M; m++) 
  667.     {
  668. newval = sum_prob[0][0][m] / nObs;
  669. delta += fabs(newval - model_trans[0][0][m][0]);
  670. model_trans[0][0][m][0] = newval;
  671.     }
  673.     // And transition probs 
  674.     for (J = 1; J < nLev; J++){
  675.       for (mm=0; mm < model_trans[J].size(); mm++)
  676. for (m = 0; m < M; m++)
  677.     for (n = 0; n < M; n++) 
  678.     {
  679.       if (sum_prob[J].size()>sum_prob[J-1].size())
  680. newval = sum_trans[J][mm][m][n] / sum_prob[J-1][mm/2][n]
  681.   /(nCh/2);
  682.       else
  683. newval =sum_trans[J][mm][m][n] / sum_prob[J-1][mm][n]/nCh;
  684.       delta += fabs(newval - model_trans[J][mm][m][n]);
  685.       model_trans[J][mm][m][n] = newval;
  686.     }
  687.     }
  688. }
  690. //-----------------------------------------------------------------------------
  691. void THMT::reorder_model()
  692. {
  693.   int level, dir, state, state2, largeststate, k;
  694.   double largeststdv, tempdouble;
  696.   // for each node
  697.   for(level = 0; level < nLev; level++){
  698.     for(dir = 0; dir < model_stdv[level].size(); dir++) {
  699.       for(state=0; state<M-1; state++) {
  700. // initialize
  701. largeststdv = -1;
  702. largeststate = M;
  703. for(state2=state; state2<M; state2++){
  704.   // search for the state with the largest standard deviation
  705.   if (model_stdv[level][dir][state2]>largeststdv){
  706.     largeststdv = model_stdv[level][dir][state2];
  707.     largeststate = state2;
  708.   }
  709. }
  711. // if the current state is not the state with the largest stdv, then
  712. // need to swap the order
  713. if(largeststate != state){
  715.   // swap the order of the stdv
  716.   model_stdv[level][dir][largeststate] = model_stdv[level][dir][state];
  717.   model_stdv[level][dir][state] = largeststdv;
  719.   // swap the order of the mean
  720.   if (!zeromean) {
  721.     tempdouble = model_mean[level][dir][largeststate];
  722.     model_mean[level][dir][largeststate] = model_mean[level][dir]
  723.       [state];
  724.     model_mean[level][dir][state] = tempdouble;
  725.   }
  727.   // swap the order of the transition matrix with parent
  728.   for(k=0; k<M; k++){
  729.     tempdouble = model_trans[level][dir][largeststate][k];
  730.     model_trans[level][dir][largeststate][k] = model_trans[level][dir]
  731.       [state][k];
  732.     model_trans[level][dir][state][k] = tempdouble;
  733.   }
  735.   // swap the order of the transition matrix with children if this
  736.   // is not a leaf node
  737.   if(level != nLev-1){
  738.     if(model_trans[level].size()==model_trans[level+1].size())
  739.       for(k=0; k<M; k++){
  740. tempdouble = model_trans[level+1][dir][k][largeststate];
  741. model_trans[level+1][dir][k][largeststate] = 
  742.   model_trans[level+1][dir][k][state];
  743. model_trans[level+1][dir][k][state] = tempdouble;
  744.       }
  745.     else if (model_trans[level].size()*2==model_trans[level+1].size())
  746.       for (k=0; k<M; k++){
  747. tempdouble =model_trans[level+1][2*dir][k][largeststate];
  748. model_trans[level+1][2*dir][k][largeststate] = 
  749.   model_trans[level+1][2*dir][k][state];
  750. model_trans[level+1][2*dir][k][state] = tempdouble;
  752. tempdouble =model_trans[level+1][2*dir+1][k][largeststate];
  753. model_trans[level+1][2*dir+1][k][largeststate] = 
  754.   model_trans[level+1][2*dir+1][k][state];
  755. model_trans[level+1][2*dir+1][k][state] = tempdouble;
  756.       } 
  757.   }
  758. }
  759.       }
  760.     }   
  761.   }
  763. }
  765. //-----------------------------------------------------------------------------
  766. void THMT::batch_train(tree<double> *trainTree, double min_delta)
  767. {
  768.   
  769.     // Assign data pointer
  770.     obs = trainTree;
  771.     // Train HMT model
  772.     train_all(min_delta);
  773. }
  776. //-----------------------------------------------------------------------------
  777. void THMT::train_all(double min_delta)
  778. {
  779.     register int count = 0, J, i, m, nNode;
  780.     register double delta = min_delta;
  781.     register double avf;
  782.     register double last_avf = -10e6;
  783.  
  784.     if (obs->nlev() == 0)
  785.       mexErrMsgTxt("ERROR in THMT::train_all(): empty training data");
  786.     if ((obs->nlev() != nLev) || (obs->nch() != nCh))
  787.       mexErrMsgTxt("ERROR in THMT::train_all(): incompatible training data");
  789.     nObs = obs->nrt();
  791.     // Allocate space for training
  792.     allocate_training();
  794.     while ((delta >= min_delta) && (count++ <= MAX_ITR)) 
  795.     {
  796. update_model(delta, avf);
  798. //if (avf < last_avf){
  799. //  mexWarnMsgTxt("WARNING: Log-likelihood decreases in training!");
  800. //  break;
  801. //}
  803. last_avf = avf;
  805. //#ifdef DEBUG
  806. //mexPrintf("count = %dndelta = %fnavf = %fn", count, delta, avf);
  807. //#endif
  808.     }
  810.     // change the model so that the state 1 always has the largest variance
  811.     // state 2 the second, and so on.....
  812.     //reorder_model();
  814.     mexPrintf("done batch-train:ncount = %dndelta = %fnavf = %fn", 
  815.       count, delta, avf);
  816. }
  819. //-----------------------------------------------------------------------------
  820. double THMT::batch_test(tree<double> *testTree)
  821. {
  822.     // Assign data pointer
  824.   obs = testTree;
  825.  
  826.   // Compute average log-likelihood
  827.   return test_all();
  828. }
  831. //-----------------------------------------------------------------------------
  832. double THMT::batch_test(char *filename)
  833. {
  834.     // Read data from file
  835.     obs = new tree<double> (filename);
  837.     // Compute average log-likelihood
  838.     double avf = test_all();
  840.     // Delete data
  841.     delete obs;
  843.     return avf;
  844. }
  847. //-----------------------------------------------------------------------------
  848. double THMT::test_all()
  849. {
  850.     register int ob_ind;
  851.     register double f, sumf = 0.0;
  852.     
  853.     
  854.     if (obs->nlev() == 0)
  855.     {
  856.       //cerr << "ERROR in THMT::test_all(): empty data"
  857.       //     << endl;
  858. return 0.0;
  859.     }
  860.     
  861.    
  862.     if ((obs->nlev() != nLev) || (obs->nch() != nCh))
  863.     {
  864.       //cerr << "ERROR in THMT::test_all(): incompatible data"
  865.       //      << endl;
  866. return 0.0;
  867.     }
  868.   
  869.     nObs = obs->nrt();
  870.   
  871.     // Allocate space for training
  872.     allocate_testing();
  874.     for (ob_ind = 0; ob_ind < nObs; ob_ind++) 
  875.     {
  876. compute_beta(ob_ind);
  877. f = compute_likelihood();
  879. //#ifdef DEBUG
  880. //cout << "ob_ind = " << ob_ind << "tf = " << f << endl;
  881. //#endif
  883. sumf += f;
  884.     }
  886.     return (sumf / double(nObs));
  887. }
  890. //-----------------------------------------------------------------------------
  891. void THMT::rnd_init_model()
  892. {
  893.     const double MAX_MEAN = 100.0;
  894.     const double MAX_STDV = 100.0;
  895.     register int J, MM, m, mm, n;
  896.     double temp;
  897.     vector<double> vprob(M);
  898.     int idum = -time(NULL);    // random seed
  900.     for (J = 0; J < nLev; J++){
  901.       for (mm=0; mm<model_mean[J].size(); mm++)
  902. for (m = 0; m < M; m++) 
  903.   {
  904.     if (zeromean)
  905.       model_mean[J][mm][m] = 0.0;
  906.     else
  907.       model_mean[J][mm][m] = (2.0*ran1(idum) - 1.0) * MAX_MEAN;
  908.     model_stdv[J][mm][m]  = ran1(idum) * MAX_STDV;
  910.     ranprobs(vprob, idum);
  911.     for (n = 0; n < M; n++)
  912.       model_trans[J][mm][n][m] = vprob[n];
  913.   }
  914.     }
  915.     
  916. }
  919. //-----------------------------------------------------------------------------
  920. void THMT::generate_data(tree<double>&obs, int n)
  921. {
  922.   tree<double> aTree(1, nCh, nLev);
  923.   register int ob_ind, J;
  924.   int idum = -time(NULL);     // random seed
  926.   // Resize output tree if necessary
  927.   if ((obs.nlev() != nLev) || (obs.nch() != nCh) || (obs.nrt() != n))
  928.     obs = tree<double>(n, nCh, nLev);
  929.   
  930.   for (ob_ind = 0; ob_ind < n; ob_ind++) {
  931.     generate_one(aTree, idum);
  932.     
  933.     for (J = 0; J < nLev; J++)
  934.       copy(aTree[J].begin(), aTree[J].end(),
  935.    obs[J].begin() + ob_ind * ipow(nCh,J));
  936.   }
  937. }
  940. //-----------------------------------------------------------------------------
  941. void THMT::generate_data(char *filename, int n)
  942. {
  943.     tree<double> genTree(n, nCh, nLev);
  944.     tree<double> aTree(1, nCh, nLev);
  945.     register int ob_ind, J;
  946.     int idum = -time(NULL);     // random seed
  948.     for (ob_ind = 0; ob_ind < n; ob_ind++) {
  949. generate_one(aTree, idum);
  951. for (J = 0; J < nLev; J++)
  952.     copy(aTree[J].begin(), aTree[J].end(),
  953.  genTree[J].begin() + ob_ind * ipow(nCh,J));
  954.     }
  956.     genTree.dump(filename);
  957. }
  960. //-----------------------------------------------------------------------------
  961. void THMT::generate_one(tree<double> &aTree, int& idum)
  962. {
  963.     tree<int> states(1, nCh, nLev);
  964.     vector<double> vprob(M);
  965.     register int J, i, nNode, m, mm;
  966.     register double mean, stdv;
  968.     // Initial state
  969.     for (m = 0; m < M; m++)
  970. vprob[m] = model_trans[0][0][m][0];
  971.     
  972.     states[0][0] = ranind(vprob, idum);
  974.     mean = model_mean[0][0][states[0][0]];
  975.     stdv = model_stdv[0][0][states[0][0]];
  976.     aTree[0][0] = mean + stdv * rangas(idum);
  977.   
  978.     // All others
  979.     for (J = 1; J < nLev; J++) 
  980.     {
  981. for (i = 0, nNode = aTree[J].size(); i < nNode; i++) 
  982. {
  983.   mm = subbandtree[J][i];
  985.   // build vector prob
  986.   for (m = 0; m < M; m++)
  987.     vprob[m] = model_trans[J][mm][m][states[J-1][i/nCh]];
  989.   states[J][i] = ranind(vprob, idum);
  990.   mean = model_mean[J][mm][states[J][i]];
  991.   stdv = model_stdv[J][mm][states[J][i]];
  992.   aTree[J][i] = mean + stdv * rangas(idum);
  993. }
  994.     }
  995. }
  998. //-----------------------------------------------------------------------------
  999. void THMT::generate_one(tree<double> &aTree, int& idum, double initval)
  1000. {
  1001.     tree<int> states(1, nCh, nLev);
  1002.     vector<double> vprob(M);
  1003.     register int J, i, nNode, m,mm;
  1004.     register double mean, stdv;
  1006.     aTree[0][0] = initval;
  1008.     /********* WRONG!!!
  1009.     // Find initial states from given initial value
  1010.     for (m = 0; m < M; m++)
  1011. vprob[m] = model_trans[0][m][0] *
  1012.     compute_g(initval, model_mean[0][m], model_stdv[0][m]);
  1014.     // It is the state with highest probability
  1015.     states[0][0] = 0;
  1016.     double maxprob = vprob[0];
  1018.     for (m = 1; m < M; m++)
  1019.     {
  1020. if (vprob[m] > maxprob)
  1021. {
  1022.     states[0][0] = m;
  1023.     maxprob = vprob[m];
  1024. }
  1025.     }
  1026.     **********/
  1028.     /***** HACK (before the end of the Millennium!!!) *****/
  1029.     if (M != 2)
  1030.       mexErrMsgTxt("Only works for 2 states");
  1032.     int smallState = (model_stdv[0][0][0] < model_stdv[0][0][1]) ? 0 : 1;
  1034.     // Cumunative probability
  1035.     double cumprob = 0.0;
  1036.     for (m = 0; m < M; m++)
  1037. cumprob += model_trans[0][0][m][0] *
  1038.     Psi(fabs(initval - model_mean[0][0][m]) / model_stdv[0][0][m]);
  1040.     if (cumprob > (0.5 + 0.5 * model_trans[0][0][smallState][0]))
  1041. states[0][0] = smallState;
  1042.     else
  1043. states[0][0] = 1 - smallState; // largeState
  1045.     /***** END HACK *****/
  1047.   
  1048.     // All others
  1049.     for (J = 1; J < nLev; J++) 
  1050.     {
  1051. for (i = 0, nNode = aTree[J].size(); i < nNode; i++) 
  1052. {
  1053.   mm = subbandtree[J][i];
  1055.   // build vector prob
  1056.   for (m = 0; m < M; m++)
  1057.     vprob[m] = model_trans[J][mm][m][states[J-1][i/nCh]];
  1059.   states[J][i] = ranind(vprob, idum);
  1060.   mean = model_mean[J][mm][states[J][i]];
  1061.   stdv = model_stdv[J][mm][states[J][i]];
  1062.   aTree[J][i] = mean + stdv * rangas(idum);
  1063. }
  1064.     }
  1066. }
  1069. //-----------------------------------------------------------------------------
  1070. tree<double>* THMT::denoise(double nvar, tree<double>* source, 
  1071.     const mxArray* stateprob)
  1072. {
  1073.   double temp = 0, *doubleptr;
  1074.   int J, i, m, mm, nNode;
  1075.   mxArray* stateprobcell;
  1077.   // Read data from file
  1078.   obs = new tree<double> (*source);
  1079.   nObs = obs->nrt();
  1081.   state_prob = tree< vector<double> >(nObs, nCh, nLev, 
  1082.       vector<double>(M));
  1083.   for (J=0; J<nLev; J++){
  1084.     stateprobcell = mxGetCell(stateprob, J);
  1085.     doubleptr = mxGetPr(stateprobcell);
  1086.     for (i = 0, nNode = state_prob[J].size(); i < nNode; i++) {
  1087.       for (m = 0; m < M; m++){
  1088. state_prob[J][i][m] = doubleptr[m*nNode+i];
  1089.       }
  1090.     }
  1091.   }
  1092.  
  1093.   for (J = 0; J < nLev; J++) {
  1094.     for (i = 0, nNode = state_prob[J].size(); i < nNode; i++) {
  1095.       mm = subbandtree[J][i%ipow(nCh,J)];
  1096.       for (m = 0; m < M; m++){
  1097. temp += state_prob[J][i][m]*model_stdv[J][mm][m]*model_stdv[J][mm][m]
  1098.   /(nvar+model_stdv[J][mm][m]*model_stdv[J][mm][m]);
  1099.       }
  1100.       (*obs)[J][i] = (*obs)[J][i]*temp;
  1101.       temp = 0;
  1102.     }
  1103.   }
  1105.   return obs;
  1106. }
  1108. //-----------------------------------------------------------------------------
  1109. void THMT::dump_model(char *filename)
  1110. {
  1111.     FILE *fp;
  1112.     register int J, MM, m, mm, n;
  1114.     fp = fopen(filename,"a");
  1115.     if (!fp) 
  1116.       mexErrMsgTxt("ERROR: can not open for writing");
  1118.     fprintf(fp, "nStates: %dn", M);
  1120.     fprintf(fp, "nLevels: %dn", nLev);
  1122.     if (zeromean)
  1123. fprintf(fp, "zeroMean: yesn");
  1124.     else 
  1125. fprintf(fp, "zeroMean: non");
  1127.     // Initial probs
  1128.     fprintf(fp, "n");
  1129.     for (m = 0; m < M; m++)
  1130. fprintf(fp, "%f ", model_trans[0][0][m][0]);
  1131.     fprintf(fp, "nn");
  1133.     // Trans probs
  1134.     for (J =  1; J < nLev; J++)
  1135.     { 
  1136.       for (m = 0; m < M; m++){
  1137. for (mm = 0; mm < model_trans[J].size(); mm++) 
  1138.     for (n = 0; n < M; n++)
  1139.       fprintf(fp, "%f ", model_trans[J][mm][m][n]);
  1140. fprintf(fp, "n");
  1141.       }
  1142.       fprintf(fp, "n");
  1143.     }
  1144.     fprintf(fp, "n");
  1146.     // Mean
  1147.     if (!zeromean){
  1148. for (J = 0; J < nLev; J++)
  1149. {
  1150.   for (mm=0; mm < model_mean[J].size(); mm++)
  1151.     for (m=0; m<M; m++)
  1152.       fprintf(fp, "%f ", model_mean[J][mm][m]);
  1154.   fprintf(fp, "n");
  1155. }
  1156. fprintf(fp, "n");
  1157.     }
  1159.     // Standard deviation
  1160.     for (J =  0; J < nLev; J++)
  1161.     {
  1162.       for (mm=0; mm < model_stdv[J].size(); mm++)
  1163. for (m = 0; m < M; m++) 
  1164.     fprintf(fp, "%f ", model_stdv[J][mm][m]);
  1165.   
  1166.       fprintf(fp, "n");
  1167.     }
  1168.     fprintf(fp, "n");
  1170.     fclose(fp);
  1171. }
  1173. //-----------------------------------------------------------------------------
  1174. void THMT::dump_model_struct(const mxArray* model)
  1175. {
  1176.     double* doubleptr;
  1177.     mxArray* pointer, *pointer2, *pointer3, *assignptr;
  1178.     int numfields, numels;
  1179.     register int J, MM, m, mm, n;
  1181.     numfields = mxGetNumberOfFields(model);
  1182.     numels = mxGetNumberOfElements(model);
  1183.     if (((numfields != 6) && (zeromean)) || 
  1184. ((numfields != 7) && (!zeromean)))
  1185.       mexErrMsgTxt("ERROR: number of fields in struct model is incorrect");
  1187.     if (numels != 1)
  1188.       mexErrMsgTxt("ERROR: Too many elements");
  1190.     if (strcmp(mxGetFieldNameByNumber(model, 0), "nstates") != 0)
  1191.       mexErrMsgTxt("Field 0 has wrong name");
  1193.     pointer = mxGetFieldByNumber(model, 0, 0);
  1194.     doubleptr = mxGetPr(pointer);
  1195.     *doubleptr = (double)M;
  1197.     if (strcmp(mxGetFieldNameByNumber(model, 1), "nlevels") != 0)
  1198.       mexErrMsgTxt("Field 1 has wrong name");    
  1200.     pointer = mxGetFieldByNumber(model, 0, 1);
  1201.     doubleptr = mxGetPr(pointer);
  1202.     *doubleptr = (double)nLev;
  1204.     if (strcmp(mxGetFieldNameByNumber(model, 2), "zeromean") != 0)
  1205.       mexErrMsgTxt("Field 2 has wrong name");
  1206.     if (zeromean){
  1207.         assignptr = mxCreateString("yes");
  1208.         mxSetFieldByNumber((mxArray*)model, 0, 2, assignptr);
  1209.     }
  1210.     else {
  1211.         assignptr = mxCreateString("no");
  1212.         mxSetFieldByNumber((mxArray*)model, 0, 2, assignptr);
  1213.     }
  1215.     if (strcmp(mxGetFieldNameByNumber(model, 3), "rootprob") != 0)
  1216.       mexErrMsgTxt("Field 3 has wrong name");
  1217.     pointer = mxGetFieldByNumber(model, 0, 3);
  1218.     if (pointer == NULL) {
  1219.       mexPrintf("%s%dn",
  1220.      "FIELD:", 3);
  1221.       mexErrMsgTxt("Above field is empty!"); 
  1222.     }
  1223.     doubleptr = mxGetPr(pointer);
  1224.     for (m = 0; m < M; m++)
  1225.         doubleptr[m] = model_trans[0][0][m][0];
  1226.  
  1227.     if (strcmp(mxGetFieldNameByNumber(model, 4), "transprob") != 0)
  1228.       mexErrMsgTxt("Field 4 has wrong name");
  1229.     pointer = mxGetFieldByNumber(model, 0, 4);
  1230.     if (pointer == NULL) {
  1231.       mexPrintf("%s%dn",
  1232.      "FIELD:", 4);
  1233.       mexErrMsgTxt("Above field is empty!"); 
  1234.     }
  1236.     // Trans probs
  1237.     for (J =  1; J < nLev; J++)
  1238.     { 
  1239.       pointer2 = mxGetCell(pointer,J-1);
  1240.       for (m = 0; m < M; m++){
  1241. for (mm = 0; mm < model_trans[J].size(); mm++) {
  1242.     pointer3 = mxGetCell(pointer2, mm);
  1243.     doubleptr = mxGetPr(pointer3);
  1244.     for (n = 0; n < M; n++)
  1245.       doubleptr[n*M+m] = model_trans[J][mm][m][n];
  1246. }
  1247.       }
  1248.     }
  1250.     // Mean
  1251.     if (!zeromean)
  1252.     {
  1253.         if (strcmp(mxGetFieldNameByNumber(model, 5), "mean") != 0)
  1254.   mexErrMsgTxt("Field 5 has wrong name");
  1255.         pointer = mxGetFieldByNumber(model, 0, 5);
  1256. if (pointer == NULL) {
  1257.   mexPrintf("%s%dn",
  1258.     "FIELD:", 5);
  1259.   mexErrMsgTxt("Above field is empty!");
  1260. }
  1261. for (J = 0; J < nLev; J++)
  1262. {
  1263.   pointer2 = mxGetCell(pointer, J);
  1264.   for (mm=0; mm < model_mean[J].size(); mm++){
  1265.     pointer3 = mxGetCell(pointer2, mm);
  1266.     doubleptr = mxGetPr(pointer3);
  1267.     for (m=0; m<M; m++)
  1268.       doubleptr[m] = model_mean[J][mm][m];
  1269.   }
  1270. }
  1271.     }
  1273.     if(!zeromean){
  1274.       if (strcmp(mxGetFieldNameByNumber(model, 6), "stdv") != 0)
  1275. mexErrMsgTxt("Field 6 has wrong name");
  1276.       pointer = mxGetFieldByNumber(model, 0, 6);
  1277.       if (pointer == NULL) {
  1278. mexPrintf("%s%dn",
  1279.   "FIELD:", 6);
  1280. mexErrMsgTxt("Above field is empty!");
  1281.       }
  1282.     }
  1283.     else{
  1284.       if (strcmp(mxGetFieldNameByNumber(model, 5), "stdv") != 0)
  1285. mexErrMsgTxt("Field 5 has wrong name");
  1286.       pointer = mxGetFieldByNumber(model, 0, 5);
  1287.       if (pointer == NULL) {
  1288. mexPrintf("%s%dn",
  1289.   "FIELD:", 5);
  1290. mexErrMsgTxt("Above field is empty!");
  1291.       } 
  1292.     }
  1294.     // Standard deviation
  1295.     for (J =  0; J < nLev; J++)
  1296.     {
  1297.       pointer2 = mxGetCell(pointer, J);
  1298.       for (mm=0; mm < model_stdv[J].size(); mm++){
  1299. pointer3 = mxGetCell(pointer2, mm);
  1300. doubleptr = mxGetPr(pointer3);
  1301. for (m = 0; m < M; m++) 
  1302.     doubleptr[m] = model_stdv[J][mm][m];
  1303.       }
  1304.     }
  1305. }
  1307. //-----------------------------------------------------------------------------
  1308. void THMT::dump_state_prob(const mxArray* stateprob)
  1309. {
  1310.     int i, J, m, nNode;
  1311.     mxArray* tempptr;
  1312.     double* dataptr;
  1314.     for (J=0; J<nLev; J++){
  1315.       tempptr = mxGetCell(stateprob, J);
  1316.       dataptr = mxGetPr(tempptr);
  1317.       for (i = 0, nNode = state_prob[J].size(); i < nNode; i++) 
  1318. for (m = 0; m < M; m++)
  1319.   dataptr[m*state_prob[J].size()+i] = state_prob[J][i][m];
  1320.     }
  1321. }  
  1322. //-----------------------------------------------------------------------------
  1323. double KLD_est(THMT model1, THMT model2, int nObservations)
  1324. {
  1325.     tree<double> genTree;
  1326.     double logden1, logden2;
  1327.     
  1328.     model1.generate_data(genTree, nObservations);
  1330.     logden1 = model1.batch_test(&genTree);
  1331.     logden2 = model2.batch_test(&genTree);
  1332.     
  1333.     return (logden1 - logden2);
  1334. }
  1337. //-----------------------------------------------------------------------------
  1338. // Compute the Kullback-Leibler distance between two discrete 
  1339. // probality mass functions
  1340. double KLD_disc(const vector<double>& prob1, const vector<double>& prob2)
  1341. {
  1342.     if (prob1.size() != prob2.size())
  1343. cerr << "KLD_disc: Two probability vectors have different length"
  1344.      << endl;
  1346.     double d = 0.0;
  1348.     for (int i = 0; i < prob1.size(); i++)
  1349. if ((prob1[i] != 0.0) && (prob2[0] != 0.0))
  1350.     d += prob1[i] * log(prob1[i] / prob2[i]);
  1352.     return d;
  1353. }
  1356. //-----------------------------------------------------------------------------
  1357. // Compute the Kullback-Leibler distance between two continuous 
  1358. // Gaussian probability desity functions
  1359. double KLD_gauss(double mean1, double stdv1, double mean2, double stdv2)
  1360. {
  1361.     double r1 = stdv1 / stdv2;
  1362.     double r2 = (mean1 - mean2) / stdv2;
  1364.     return 0.5 * (-2*log(r1) - 1 + r1*r1 + r2*r2); 
  1365. }
  1368. //-----------------------------------------------------------------------------
  1369. double KLD_upb(THMT model1, THMT model2)
  1370. {
  1371.     int J, m, n, dir, maxdir;
  1372.     int M, nCh, nLev;
  1373.     double test_sum;
  1374.    
  1375.     M = model1.M;
  1376.     if (model2.M != M)
  1377. mexErrMsgTxt("KLD_upb: Incompatible models.");
  1379.     nCh = model1.nCh;
  1380.     if (model2.nCh != nCh)
  1381. mexErrMsgTxt("KLD_upb: Incompatible models.");
  1383.     nLev = model1.nLev;
  1384.     if (model2.nLev != nLev)
  1385. mexErrMsgTxt("KLD_upb: Incompatible models.");
  1387.     // assume the lowest level has the largest number of directions
  1388.     maxdir = model1.model_stdv[nLev-1].size();
  1391.     vector<vector<double> > D(maxdir, vector<double>(M));
  1392.     vector<vector<double> > d(maxdir, vector<double>(M));
  1393.     vector<double> trans1(M);
  1394.     vector<double> trans2(M);
  1395.  
  1396.     // Initial: lowest level
  1397.     J = nLev - 1;
  1398.     for (dir = 0; dir < model1.model_stdv[J].size(); dir++){
  1399.       for (m = 0; m < M; m++)
  1400. D[dir][m] = KLD_gauss(model1.model_mean[J][dir][m], 
  1401.  model1.model_stdv[J][dir][m],
  1402.  model2.model_mean[J][dir][m], 
  1403.  model2.model_stdv[J][dir][m]);
  1404.     } 
  1405.       
  1406.     // DEBUG
  1407.     //mexPrintf("%sn", "Lowest level: ");
  1408.     //for (dir = 0; dir < model1.model_stdv[J].size(); dir++)
  1409.     //  for (m = 0; m < M; m++)
  1410.     // mexPrintf("%d %d %f n", dir, m, D[dir][m]);
  1412.     // Induction:
  1413.     for (J = nLev-1; J > 0; J--) {
  1414.       for (dir = 0; dir < model1.model_stdv[J-1].size(); dir++){
  1415. for (m = 0; m < M; m++) {
  1417.   d[dir][m] = KLD_gauss(model1.model_mean[J-1][dir][m], 
  1418. model1.model_stdv[J-1][dir][m],
  1419. model2.model_mean[J-1][dir][m], 
  1420. model2.model_stdv[J-1][dir][m]);
  1421.  
  1422.   if ( model1.model_stdv[J].size() == model1.model_stdv[J-1].size()){
  1423.     for (n=0; n<M; n++){
  1424.       trans1[n] = model1.model_trans[J][dir][n][m];
  1425.       trans2[n] = model2.model_trans[J][dir][n][m];
  1426.     }
  1427.     d[dir][m] += nCh * KLD_disc(trans1,trans2);
  1428.   }
  1429.   else if (model1.model_stdv[J].size() == 
  1430.    2*model1.model_stdv[J-1].size()){
  1431.     for (n=0; n<M; n++){
  1432.       trans1[n] = model1.model_trans[J][dir*2][n][m];
  1433.       trans2[n] = model2.model_trans[J][dir*2][n][m];
  1434.     }
  1436.     d[dir][m] += nCh/2 * KLD_disc(trans1,trans2);
  1437.     for (n=0; n<M; n++){
  1438.       trans1[n] = model1.model_trans[J][dir*2+1][n][m];
  1439.       trans2[n] = model2.model_trans[J][dir*2+1][n][m];
  1440.     }
  1441.     d[dir][m] += nCh/2 * KLD_disc(trans1,trans2);
  1442.   }
  1443.   else
  1444.     mexErrMsgTxt("Error: Multiple parents for one child");
  1447.   for (n = 0; n < M; n++)
  1448.     if ( model1.model_stdv[J].size() == model1.model_stdv[J-1].size())
  1449.       d[dir][m] += nCh * model1.model_trans[J][dir][n][m] 
  1450. * D[dir][n]; 
  1451.     else {
  1452.       d[dir][m] += nCh/2 * model1.model_trans[J][dir*2][n][m] 
  1453. * D[dir*2][n];
  1454.       d[dir][m] += nCh/2 * model1.model_trans[J][dir*2+1][n][m] 
  1455. * D[dir*2+1][n];
  1456.     }
  1457. }
  1458.       }
  1459.       // DEBUG:
  1461.       //for (dir = 0; dir < model1.model_stdv[J-1].size(); dir++)
  1462.       // for (m = 0; m < M; m++)
  1463.       //   mexPrintf("%d %d %f %fn", dir, m, D[dir][m], d[dir][m]);
  1465.       // updating the temporary distance vector D
  1466.       D = d;
  1467.     }
  1468.  
  1469.     // Final:
  1470.     double dist;
  1471.     for (n=0; n<M; n++){
  1472.       trans1[n] = model1.model_trans[0][0][n][0];
  1473.       trans2[n] = model2.model_trans[0][0][n][0];
  1474.     }
  1476.     dist = KLD_disc(trans1, trans2);
  1478.     // DEBUG
  1479.     //mexPrintf("%s %f", "Final: KLD_disc = ", dist);
  1481.     for (m = 0; m < M; m++)
  1482. dist += model1.model_trans[0][0][m][0] * D[0][m];
  1484.     return dist;
  1485. }
  1487. /*********************************************************************
  1488.     TEMPLATES INSTANCIATION !!!
  1489. *************************************************************/
  1490. /*template class matrix<float>;
  1492. template class vector<double>;
  1493. template class tree<double>;
  1494. template class tree<int>;
  1495. template class tree<vector<double> >; */