1 /*  _______________________________________________________________________
2 
3     PECOS: Parallel Environment for Creation Of Stochastics
4     Copyright (c) 2011, Sandia National Laboratories.
5     This software is distributed under the GNU Lesser General Public License.
6     For more information, see the README file in the top Pecos directory.
7     _______________________________________________________________________ */
8 
9 #include "pecos_stat_util.hpp"
10 #include "RandomVariable.hpp"
11 
12 static const char rcsId[]="@(#) $Id: pecos_stat_util.cpp 4768 2007-12-17 17:49:32Z mseldre $";
13 
14 namespace Pecos {
15 
16 
int_range_to_xy_pdf(int l_bnd,int u_bnd,RealArray & x_val,RealArray & y_val)17 void int_range_to_xy_pdf(int l_bnd,        int u_bnd,
18 			 RealArray& x_val, RealArray& y_val)
19 {
20   int i, num_params = u_bnd - l_bnd + 1;
21   x_val.resize(num_params);  y_val.assign(num_params, 1.);
22   for (i=0; i<num_params; ++i)
23     x_val[i] = (Real)(l_bnd + i);
24 }
25 
26 
bins_to_xy_cdf(const RealRealMap & h_bin_prs,RealArray & x_val,RealArray & y_val)27 void bins_to_xy_cdf(const RealRealMap& h_bin_prs,
28 		    RealArray& x_val, RealArray& y_val)
29 {
30   size_t i, num_params = h_bin_prs.size(), last_index = num_params - 1;
31   RRMCIter cit = h_bin_prs.begin();
32 
33   // Note: LHS continuous linear accumulates CDF with first y=0 and last y=1
34   x_val.resize(num_params);  y_val.resize(num_params);
35   for (i=0; i<num_params; ++i, ++cit)
36     x_val[i] = cit->first;
37   y_val[0] = 0.;  cit = h_bin_prs.begin();
38   for (i=0; i<last_index; ++i, ++cit)
39     y_val[i+1] = y_val[i] + cit->second * (x_val[i+1] - x_val[i]);
40   // normalize if necessary (h_bin_pairs should have normalized counts)
41   Real& cdf_last = y_val[last_index];
42   if (cdf_last != 1.) {
43     for (i=1; i<last_index; ++i)
44       y_val[i] /= cdf_last;
45     cdf_last = 1.;
46   }
47 #ifdef DEBUG
48   for (i=0; i<num_params; ++i)
49     PCout << "hbuv cdf: x_val[" << i << "] is " << x_val[i]
50 	  << " y_val[" << i << "] is " << y_val[i] << '\n';
51 #endif // DEBUG
52 }
53 
54 
intervals_to_xy_cdf(const RealRealPairRealMap & ci_bpa,RealArray & x_val,RealArray & y_val)55 void intervals_to_xy_cdf(const RealRealPairRealMap& ci_bpa,
56 			 RealArray& x_val, RealArray& y_val)
57 {
58   RealArray prob_dens;
59   intervals_to_xy_pdf(ci_bpa, x_val, prob_dens);
60   size_t i, num_params = x_val.size(), last_index = num_params - 1;
61 
62   // LHS continuous linear accumulates CDF with first y=0 and last y=1:
63   // put the densities in a cumulative format necessary for LHS histograms.
64   y_val.resize(num_params);
65   y_val[0] = 0.;
66   Real pdf;
67   for (i=1; i<num_params; ++i) {
68     pdf = prob_dens[i-1];
69     y_val[i] = (pdf > 0.) ? y_val[i-1] + pdf * (x_val[i] - x_val[i-1]) :
70                             y_val[i-1] + 1.e-4; // handle case of a gap
71   }
72   // normalize if necessary
73   Real& cdf_last = y_val[last_index];
74   if (cdf_last != 1.) {
75     for (i=1; i<last_index; ++i)
76       y_val[i] /= cdf_last;
77     cdf_last = 1.;
78   }
79 #ifdef DEBUG
80   for (i=0; i<num_params; ++i)
81     PCout << "ciuv cdf: x_val[" << i << "] is " << x_val[i]
82 	  << " y_val[" << i << "] is " << y_val[i] << '\n';
83 #endif // DEBUG
84 }
85 
86 
87 /*
88 Note: these are not active since RandomVariable type is insufficient to
89       define the source variable type in the presence of probability
90       transformations (e.g., from design/epistemic/state to std uniform).
91 
92 void design_state_subset(const std::vector<RandomVariable>& random_vars,
93 			 BitArray& subset, size_t start_set, size_t num_set)
94 {
95   size_t i, num_rv = random_vars.size(), end_set = start_set + num_set;
96   subset.resize(num_rv, false); // init bits to false
97   for (i=start_set; i<end_set; ++i)
98     // activate design + state vars
99     switch (random_vars[i].type()) {
100     case CONTINUOUS_RANGE:    case DISCRETE_RANGE: case DISCRETE_SET_INT:
101     case DISCRETE_SET_STRING: case DISCRETE_SET_REAL:
102       subset.set(i); break;
103     }
104 }
105 
106 
107 void uncertain_subset(const std::vector<RandomVariable>& random_vars,
108 		      BitArray& subset)
109 {
110   size_t i, num_rv = random_vars.size();
111   subset.resize(num_rv, true); // init bits to true
112   for (i=0; i<num_rv; ++i)
113     // deactivate complement of uncertain vars
114     switch (random_vars[i].type()) {
115     case CONTINUOUS_RANGE:    case DISCRETE_RANGE: case DISCRETE_SET_INT:
116     case DISCRETE_SET_STRING: case DISCRETE_SET_REAL:
117       subset.reset(i); break;
118     }
119 }
120 
121 
122 void aleatory_uncertain_subset(const std::vector<RandomVariable>& random_vars,
123 			       BitArray& subset)
124 {
125   size_t i, num_rv = random_vars.size();
126   subset.resize(num_rv, true); // init bits to true
127   for (i=0; i<num_rv; ++i)
128     // deactivate complement of aleatory uncertain vars
129     switch (random_vars[i].type()) {
130     case CONTINUOUS_RANGE:    case DISCRETE_RANGE: case DISCRETE_SET_INT:
131     case DISCRETE_SET_STRING: case DISCRETE_SET_REAL:
132     case CONTINUOUS_INTERVAL_UNCERTAIN: case DISCRETE_INTERVAL_UNCERTAIN:
133     case DISCRETE_UNCERTAIN_SET_INT:    case DISCRETE_UNCERTAIN_SET_STRING:
134     case DISCRETE_UNCERTAIN_SET_REAL:
135       subset.reset(i); break;
136     }
137 }
138 
139 
140 void epistemic_uncertain_subset(const std::vector<RandomVariable>& random_vars,
141 				BitArray& subset)
142 {
143   size_t i, num_rv = random_vars.size();
144   subset.resize(num_rv, false); // init bits to false
145   for (i=0; i<num_rv; ++i)
146     // activate epistemic uncertain vars
147     switch (random_vars[i].type()) {
148     case CONTINUOUS_INTERVAL_UNCERTAIN: case DISCRETE_INTERVAL_UNCERTAIN:
149     case DISCRETE_UNCERTAIN_SET_INT:    case DISCRETE_UNCERTAIN_SET_STRING:
150     case DISCRETE_UNCERTAIN_SET_REAL:
151       subset.set(i); break;
152     }
153 }
154 */
155 
156 } // namespace Pecos
157