| 357 | |
| 358 | |
| 359 | void SIR_scheme::roughen_minmax (ColMatrix& P, Float K) const |
| 360 | /* Roughening |
| 361 | * Uses algorithm from Ref[1] using max-min in each state of P |
| 362 | * K is scaling factor for roughening noise |
| 363 | * unique_samples is unchanged as roughening is used to postprocess observe resamples |
| 364 | * Numerical collapse of P |
| 365 | * P with very small or zero range result in minimal roughening |
| 366 | * Exceptions: |
| 367 | * none |
| 368 | * unchanged: P |
| 369 | */ |
| 370 | { |
| 371 | using namespace std; |
| 372 | // Scale Sigma by constant and state dimensions |
| 373 | Float SigmaScale = K * pow (Float(P.size2()), -1/Float(x_size)); |
| 374 | |
| 375 | // Find min and max states in all P, precond P not empty |
| 376 | Vec xmin(x_size); noalias(xmin) = column(P,0); |
| 377 | Vec xmax(x_size); noalias(xmax) = xmin; |
| 378 | ColMatrix::iterator2 pi = P.begin2(); |
| 379 | while (pi != P.end2()) // Loop includes 0 to simplify code |
| 380 | { |
| 381 | Vec::iterator mini = xmin.begin(); |
| 382 | Vec::iterator maxi = xmax.begin(); |
| 383 | |
| 384 | #ifdef BOOST_UBLAS_NO_NESTED_CLASS_RELATION |
| 385 | for (ColMatrix::iterator1 xpi = begin(pi, ublas::iterator2_tag()); xpi != end(pi, ublas::iterator2_tag()); ++xpi) |
| 386 | #else |
| 387 | for (ColMatrix::iterator1 xpi = pi.begin(); xpi != pi.end(); ++xpi) |
| 388 | #endif |
| 389 | { |
| 390 | if (*xpi < *mini) *mini = Float(*xpi); // ISSUE mixed type proxy assignment |
| 391 | if (*xpi > *maxi) *maxi = Float(*xpi); |
| 392 | ++mini; ++maxi; |
| 393 | } |
| 394 | ++pi; |
| 395 | } |
| 396 | // Roughening st.dev max-min |
| 397 | Vec rootq(x_size); |
| 398 | rootq = xmax; |
| 399 | noalias(rootq) -= xmin; |
| 400 | rootq *= SigmaScale; |
| 401 | // Apply roughening predict based on scaled variance |
| 402 | DenseVec n(x_size); |
| 403 | for (pi = P.begin2(); pi != P.end2(); ++pi) { |
| 404 | random.normal (n); // independent zero mean normal |
| 405 | // multiply elements by std dev |
| 406 | for (DenseVec::iterator ni = n.begin(); ni != n.end(); ++ni) { |
| 407 | *ni *= rootq[ni.index()]; |
| 408 | } |
| 409 | ColMatrix::Column Pi(P,pi.index2()); |
| 410 | noalias(n) += Pi; // add to P |
| 411 | Pi = n; |
| 412 | } |
| 413 | } |
| 414 | |
| 415 | |
| 416 | /* |