00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #ifndef FORWARD_ALGORITHM_H
00021 #define FORWARD_ALGORITHM_H
00022 #include <tr1/tuple>
00023 #include <tr1/array>
00024 #include <TooN/TooN.h>
00025 #include <vector>
00026 #include <cmath>
00027
00028
00029
00030
00031
00032
00033
00034 inline double ln(double x)
00035 {
00036 if(x == 0)
00037 return -1e100;
00038 else
00039 return std::log(x);
00040 }
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153 template<int States, class Btype, class Otype> std::tr1::tuple<double, TooN::Vector<Btype::NumParameters>, TooN::Matrix<Btype::NumParameters> > forward_algorithm_hessian(TooN::Matrix<States> A, TooN::Vector<States> pi, const Btype& B, const std::vector<Otype>& O, bool compute_deriv=1, bool compute_hessian=1)
00154 {
00155 using namespace TooN;
00156 using namespace std;
00157 using namespace std::tr1;
00158
00159 if(compute_hessian == 1)
00160 compute_deriv=1;
00161
00162 static const int M = Btype::NumParameters;
00163 int states = pi.size();
00164
00165
00166 vector<array<double, States> > delta(O.size());
00167
00168
00169 vector<array<Vector<M>,States > > diff_delta(O.size());
00170
00171
00172
00173 vector<array<Matrix<M>,States > > hess_delta(O.size());
00174
00175
00176
00177 for(int i=0; i < states; i++)
00178 {
00179 delta[0][i] = ln(pi[i]) + B.log(i, O[0]);
00180
00181 if(compute_deriv)
00182 diff_delta[0][i] = B.diff_log(i, O[0]);
00183
00184 if(compute_hessian)
00185 hess_delta[0][i] = B.hess_log(i, O[0]);
00186 }
00187
00188
00189
00190 for(unsigned int t=1; t < O.size(); t++)
00191 {
00192 for(int j=0; j < states; j++)
00193 {
00194 double Ztj = -HUGE_VAL;
00195 for(int i=0; i < states; i++)
00196 Ztj = max(Ztj, delta[t-1][i] + ln(A[i][j]));
00197
00198 double sum=0;
00199 for(int i=0; i < states; i++)
00200 sum += exp(delta[t-1][i] + ln(A[i][j]) - Ztj);
00201
00202 delta[t][j] = B.log(j, O[t]) + Ztj + ln(sum);
00203
00204 if(compute_deriv)
00205 {
00206 Vector<M> sum_top = Zeros;
00207 for(int i=0; i < states; i++)
00208 sum_top += diff_delta[t-1][i] * exp(delta[t-1][i] + ln(A[i][j]) - Ztj);
00209
00210 diff_delta[t][j] = B.diff_log(j, O[t]) + (sum_top) / sum;
00211
00212 if(compute_hessian)
00213 {
00214 Matrix<M> sum_top2 = Zeros;
00215 for(int i=0; i < states; i++)
00216 sum_top2 += exp(delta[t-1][i] + ln(A[i][j]) - Ztj) * ( hess_delta[t-1][i] + diff_delta[t-1][i].as_col() * diff_delta[t-1][i].as_row());
00217
00218 hess_delta[t][j] = B.hess_log(j, O[t]) + sum_top2 / sum - sum_top.as_col() * sum_top.as_row() / (sum*sum);
00219 }
00220 }
00221 }
00222 }
00223
00224
00225 double Z = -HUGE_VAL;
00226 for(int i=0; i < states; i++)
00227 Z = max(Z, delta.back()[i]);
00228
00229 double sum =0;
00230 for(int i=0; i < states; i++)
00231 sum += exp(delta.back()[i] - Z);
00232
00233 double log_prob = Z + ln(sum);
00234
00235
00236 Vector<M> diff_log = Zeros;
00237
00238
00239 for(int i=0; compute_deriv && i < states; i++)
00240 diff_log += exp(delta.back()[i] - log_prob)*diff_delta.back()[i];
00241
00242 Matrix<M> hess_log = Zeros;
00243
00244
00245 for(int i=0; compute_hessian && i < states; i++)
00246 hess_log += exp(delta.back()[i] - log_prob) * (hess_delta.back()[i] + diff_delta.back()[i].as_col() * diff_delta.back()[i].as_row());
00247
00248 hess_log -= diff_log.as_col() * diff_log.as_row();
00249
00250
00251 return make_tuple(log_prob, diff_log, hess_log);
00252 }
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263 template<int States, class Btype, class Otype> double forward_algorithm(TooN::Matrix<States> A, TooN::Vector<States> pi, const Btype& B, const std::vector<Otype>& O)
00264 {
00265 using namespace TooN;
00266 using namespace std;
00267 using namespace std::tr1;
00268
00269 int states = pi.size();
00270
00271
00272 vector<array<double, States> > delta(O.size());
00273
00274
00275
00276 for(int i=0; i < states; i++)
00277 delta[0][i] = ln(pi[i]) + B.log(i, O[0]);
00278
00279
00280
00281 for(unsigned int t=1; t < O.size(); t++)
00282 {
00283 for(int j=0; j < states; j++)
00284 {
00285 double Ztj = -HUGE_VAL;
00286 for(int i=0; i < states; i++)
00287 Ztj = max(Ztj, delta[t-1][i] + ln(A[i][j]));
00288
00289 double sum=0;
00290 for(int i=0; i < states; i++)
00291 sum += exp(delta[t-1][i] + ln(A[i][j]) - Ztj);
00292
00293 delta[t][j] = B.log(j, O[t]) + Ztj + ln(sum);
00294 }
00295 }
00296
00297
00298 double Z = -HUGE_VAL;
00299 for(int i=0; i < states; i++)
00300 Z = max(Z, delta.back()[i]);
00301
00302 double sum =0;
00303 for(int i=0; i < states; i++)
00304 sum += exp(delta.back()[i] - Z);
00305
00306 double log_prob = Z + ln(sum);
00307
00308 return log_prob;
00309 }
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321 template<int States, class Btype, class Otype> std::pair<double, TooN::Vector<Btype::NumParameters> > forward_algorithm_deriv(TooN::Matrix<States> A, TooN::Vector<States> pi, const Btype& B, const std::vector<Otype>& O)
00322 {
00323 using namespace std::tr1;
00324 double p;
00325 TooN::Vector<Btype::NumParameters> v;
00326 tie(p,v, ignore) = forward_algorithm_hessian(A, pi, B, O, 1, 0);
00327 return make_pair(p,v);
00328 }
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340 template<int States, class Btype, class Otype>
00341 std::vector<std::tr1::array<double, States> >
00342 forward_algorithm_delta(TooN::Matrix<States> A, TooN::Vector<States> pi, const Btype& B, const std::vector<Otype>& O)
00343 {
00344 using namespace TooN;
00345 using namespace std;
00346 using namespace std::tr1;
00347
00348 int states = pi.size();
00349
00350
00351 vector<array<double, States> > delta(O.size());
00352
00353
00354
00355 for(int i=0; i < states; i++)
00356 delta[0][i] = ln(pi[i]) + B.log(i, O[0]);
00357
00358
00359
00360
00361 for(unsigned int t=1; t < O.size(); t++)
00362 {
00363 for(int j=0; j < states; j++)
00364 {
00365 double Ztj = -HUGE_VAL;
00366 for(int i=0; i < states; i++)
00367 Ztj = max(Ztj, delta[t-1][i] + ln(A[i][j]));
00368
00369 double sum=0;
00370 for(int i=0; i < states; i++)
00371 sum += exp(delta[t-1][i] + ln(A[i][j]) - Ztj);
00372
00373 delta[t][j] = B.log(j, O[t]) + Ztj + ln(sum);
00374 }
00375 }
00376
00377 return delta;
00378 }
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389 template<int States, class Btype, class Otype>
00390 void forward_algorithm_delta2(TooN::Matrix<States> A, TooN::Vector<States> pi, const Btype& B, const std::vector<Otype>& O, std::vector<std::tr1::array<double, States> >& delta)
00391 {
00392 using namespace TooN;
00393 using namespace std;
00394 using namespace std::tr1;
00395
00396 int states = pi.size();
00397
00398
00399 delta.resize(O.size());
00400
00401
00402
00403 for(int i=0; i < states; i++)
00404 delta[0][i] = ln(pi[i]) + B.log(i, O[0]);
00405
00406 Matrix<States> lA;
00407 for(int r=0; r < States; r++)
00408 for(int c=0; c < States; c++)
00409 lA[r][c] = ln(A[r][c]);
00410
00411
00412
00413
00414 for(unsigned int t=1; t < O.size(); t++)
00415 {
00416 for(int j=0; j < states; j++)
00417 {
00418 double Ztj = -HUGE_VAL;
00419 for(int i=0; i < states; i++)
00420 Ztj = max(Ztj, delta[t-1][i] + lA[i][j]);
00421
00422 double sum=0;
00423 for(int i=0; i < states; i++)
00424 sum += exp(delta[t-1][i] + lA[i][j] - Ztj);
00425
00426 delta[t][j] = B.log(j, O[t]) + Ztj + ln(sum);
00427 }
00428 }
00429 }
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441 template<int States, class Btype, class Otype>
00442 std::pair<std::vector<std::tr1::array<double, States> >, std::vector<std::tr1::array<double, States> > >
00443 forward_backward_algorithm(TooN::Matrix<States> A, TooN::Vector<States> pi, const Btype& B, const std::vector<Otype>& O)
00444 {
00445 using namespace TooN;
00446 using namespace std;
00447 using namespace std::tr1;
00448
00449 int states = pi.size();
00450
00451
00452 vector<array<double, States> > delta = forward_algorithm_delta(A, pi, B, O);
00453
00454
00455
00456 vector<array<double, States> > epsilon(O.size());
00457
00458
00459 for(int i=0; i < states; i++)
00460 epsilon[O.size()-1][i] = 0;
00461
00462
00463 for(int t=O.size()-2; t >= 0; t--)
00464 {
00465 for(int i=0; i < states; i++)
00466 {
00467
00468 double Z = -HUGE_VAL;
00469
00470 for(int j=0; j < states; j++)
00471 Z = max(Z, ln(A[i][j]) + B.log(j, O[t+1]) + epsilon[t+1][j]);
00472
00473 double sum=0;
00474 for(int j= 0; j < states; j++)
00475 sum += exp(ln(A[i][j]) + B.log(j, O[t+1]) + epsilon[t+1][j] - Z);
00476
00477 epsilon[t][i] = ln(sum) + Z;
00478 }
00479 }
00480
00481 return make_pair(delta, epsilon);
00482 }
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498 template<class A, class Rng> int select_random_element(const A& v, const double scale, Rng& rng)
00499 {
00500 double total=0, choice = rng()*scale;
00501
00502 for(int i=0; i < (int)v.size(); i++)
00503 {
00504 total += v[i];
00505 if(choice <= total)
00506 return i;
00507 }
00508 return v.size()-1;
00509 }
00510
00511
00512
00513
00514
00515
00516 template<int N, class Rng> int sample_unscaled_log(std::tr1::array<double, N> a, Rng& rng)
00517 {
00518 double hi = *max_element(a.begin(), a.end());
00519 double sum=0;
00520
00521 for(unsigned int i=0; i < a.size(); i++)
00522 {
00523 a[i] = exp(a[i] - hi);
00524 sum += a[i];
00525 }
00526
00527 return select_random_element(a, sum, rng);
00528 }
00529
00530
00531
00532
00533
00534
00535
00536
00537 template<int States, class StateType, class Rng>
00538 std::vector<StateType> backward_sampling(TooN::Matrix<States> A, const std::vector<std::tr1::array<double, States> >& delta, Rng& rng)
00539 {
00540
00541 for(int r=0; r < A.num_rows(); r++)
00542 for(int c=0; c < A.num_cols(); c++)
00543 A[r][c] = ln(A[r][c]);
00544
00545 std::vector<StateType> samples(delta.size());
00546
00547 samples.back() = sample_unscaled_log<States, Rng>(delta.back(), rng);
00548
00549
00550
00551 for(int i=delta.size()-2; i >= 0; i--)
00552 {
00553 std::tr1::array<double, States> reverse_probabilities = delta[i];
00554
00555 for(int j=0; j < States; j++)
00556 reverse_probabilities[j] += A[j][samples[i+1]];
00557
00558 samples[i] = sample_unscaled_log<States, Rng>(reverse_probabilities, rng);
00559 }
00560
00561 return samples;
00562 }
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572 template<int States>
00573 std::vector<int> backward_sampling(const TooN::Matrix<States>& A, const std::vector<std::tr1::array<double, States> >& delta)
00574 {
00575 RngDrand48 d;
00576 return backward_sampling<States, int, RngDrand48>(A, delta, d);
00577 }
00578 */
00579
00580 #endif