Rev 3402 | Go to most recent revision | Details | Last modification | View Log | RSS feed
| Rev | Author | Line No. | Line |
|---|---|---|---|
| 3352 | fmk | 1 | /* ****************************************************************** ** |
| 2 | ** OpenSees - Open System for Earthquake Engineering Simulation ** |
||
| 3 | ** Pacific Earthquake Engineering Research Center ** |
||
| 4 | ** ** |
||
| 5 | ** ** |
||
| 6 | ** (C) Copyright 2001, The Regents of the University of California ** |
||
| 7 | ** All Rights Reserved. ** |
||
| 8 | ** ** |
||
| 9 | ** Commercial use of this program without express permission of the ** |
||
| 10 | ** University of California, Berkeley, is strictly prohibited. See ** |
||
| 11 | ** file 'COPYRIGHT' in main directory for information on usage and ** |
||
| 12 | ** redistribution, and for a DISCLAIMER OF ALL WARRANTIES. ** |
||
| 13 | ** ** |
||
| 14 | ** Developed by: ** |
||
| 15 | ** Frank McKenna (fmckenna@ce.berkeley.edu) ** |
||
| 16 | ** Gregory L. Fenves (fenves@ce.berkeley.edu) ** |
||
| 17 | ** Filip C. Filippou (filippou@ce.berkeley.edu) ** |
||
| 18 | ** ** |
||
| 19 | ** Reliability module developed by: ** |
||
| 20 | ** Terje Haukaas (haukaas@ce.berkeley.edu) ** |
||
| 21 | ** Armen Der Kiureghian (adk@ce.berkeley.edu) ** |
||
| 22 | ** ** |
||
| 23 | ** ****************************************************************** */ |
||
| 24 | |||
| 25 | // $Revision: 1.1 $ |
||
| 26 | // $Date: 2008-03-13 22:22:25 $ |
||
| 27 | // $Source: /usr/local/cvs/OpenSees/SRC/reliability/analysis/analysis/ImportanceSamplingAnalysis.cpp,v $ |
||
| 28 | |||
| 29 | // |
||
| 30 | // Written by Terje Haukaas (haukaas@ce.berkeley.edu) |
||
| 31 | // |
||
| 32 | |||
| 33 | #include <ImportanceSamplingAnalysis.h> |
||
| 34 | #include <ReliabilityDomain.h> |
||
| 35 | #include <ReliabilityAnalysis.h> |
||
| 36 | #include <LimitStateFunction.h> |
||
| 37 | #include <LimitStateFunctionIter.h> |
||
| 38 | #include <ProbabilityTransformation.h> |
||
| 39 | #include <NatafProbabilityTransformation.h> |
||
| 40 | #include <GFunEvaluator.h> |
||
| 41 | #include <BasicGFunEvaluator.h> |
||
| 42 | #include <RandomNumberGenerator.h> |
||
| 43 | #include <RandomVariable.h> |
||
| 44 | #include <NormalRV.h> |
||
| 45 | #include <Vector.h> |
||
| 46 | #include <Matrix.h> |
||
| 47 | #include <MatrixOperations.h> |
||
| 48 | #include <NormalRV.h> |
||
| 49 | #include <math.h> |
||
| 50 | #include <stdlib.h> |
||
| 51 | #include <string.h> |
||
| 52 | |||
| 53 | #include <fstream> |
||
| 54 | #include <iomanip> |
||
| 55 | #include <iostream> |
||
| 56 | using std::ifstream; |
||
| 57 | using std::ios; |
||
| 58 | using std::setw; |
||
| 59 | using std::setprecision; |
||
| 60 | using std::setiosflags; |
||
| 61 | |||
| 62 | |||
| 63 | ImportanceSamplingAnalysis::ImportanceSamplingAnalysis( ReliabilityDomain *passedReliabilityDomain, |
||
| 64 | ProbabilityTransformation *passedProbabilityTransformation, |
||
| 65 | GFunEvaluator *passedGFunEvaluator, |
||
| 66 | RandomNumberGenerator *passedRandomNumberGenerator, |
||
| 67 | Tcl_Interp *passedInterp, |
||
| 68 | int passedNumberOfSimulations, |
||
| 69 | double passedTargetCOV, |
||
| 70 | double passedSamplingStdv, |
||
| 71 | int passedPrintFlag, |
||
| 72 | TCL_Char *passedFileName, |
||
| 73 | Vector *pStartPoint, |
||
| 74 | int passedAnalysisTypeTag) |
||
| 75 | :ReliabilityAnalysis() |
||
| 76 | { |
||
| 77 | theReliabilityDomain = passedReliabilityDomain; |
||
| 78 | theProbabilityTransformation = passedProbabilityTransformation; |
||
| 79 | theGFunEvaluator = passedGFunEvaluator; |
||
| 80 | theRandomNumberGenerator = passedRandomNumberGenerator; |
||
| 81 | interp = passedInterp; |
||
| 82 | numberOfSimulations = passedNumberOfSimulations; |
||
| 83 | targetCOV = passedTargetCOV; |
||
| 84 | samplingStdv = passedSamplingStdv; |
||
| 85 | printFlag = passedPrintFlag; |
||
| 86 | strcpy(fileName,passedFileName); |
||
| 87 | startPoint = pStartPoint; |
||
| 88 | analysisTypeTag = passedAnalysisTypeTag; |
||
| 89 | } |
||
| 90 | |||
| 91 | |||
| 92 | |||
| 93 | |||
| 94 | ImportanceSamplingAnalysis::~ImportanceSamplingAnalysis() |
||
| 95 | { |
||
| 96 | |||
| 97 | } |
||
| 98 | |||
| 99 | |||
| 100 | |||
| 101 | int |
||
| 102 | ImportanceSamplingAnalysis::analyze(void) |
||
| 103 | { |
||
| 104 | |||
| 105 | // Alert the user that the simulation analysis has started |
||
| 106 | opserr << "Sampling Analysis is running ... " << endln; |
||
| 107 | |||
| 108 | // Declaration of some of the data used in the algorithm |
||
| 109 | double gFunctionValue; |
||
| 110 | int result, I, k = 1, seed = 1; |
||
| 111 | double det_covariance, phi, h, q; |
||
| 112 | int numRV = theReliabilityDomain->getNumberOfRandomVariables(); |
||
| 113 | int numLsf = theReliabilityDomain->getNumberOfLimitStateFunctions(); |
||
| 114 | Vector startValues(numRV); |
||
| 115 | Vector x(numRV); |
||
| 116 | Vector z(numRV); |
||
| 117 | Vector u(numRV); |
||
| 118 | Vector randomArray(numRV); |
||
| 119 | static NormalRV aStdNormRV(1,0.0,1.0,0.0); |
||
| 120 | bool failureHasOccured = false; |
||
| 121 | |||
| 122 | det_covariance = pow(samplingStdv, numRV); |
||
| 123 | |||
| 124 | // Pre-compute some factors to minimize computations inside simulation loop |
||
| 125 | static const double twopi = 2.0*acos(-1.0); |
||
| 126 | double factor1 = 1.0 / ( pow(twopi,0.5*numRV)); |
||
| 127 | //double factor2 = 1.0 / ( pow(twopi,0.5*numRV) * sqrt(det_covariance) ); |
||
| 128 | double factor2 = factor1 / sqrt(det_covariance); |
||
| 129 | |||
| 130 | |||
| 131 | Vector sum_q(numLsf); |
||
| 132 | Vector sum_q_squared(numLsf); |
||
| 133 | double var_qbar; |
||
| 134 | double pfIn; |
||
| 135 | double CovIn; |
||
| 136 | char restartFileName[256]; |
||
| 137 | sprintf(restartFileName,"%s_%s","restart",fileName); |
||
| 138 | |||
| 139 | if (analysisTypeTag == 1) { |
||
| 140 | // Possible read data from file if this is a restart simulation |
||
| 141 | if (printFlag == 2) { |
||
| 142 | ifstream inputFile( restartFileName, ios::in ); |
||
| 143 | inputFile >> k; |
||
| 144 | inputFile >> seed; |
||
| 145 | if (k == 1 && seed == 1) { |
||
| 146 | } |
||
| 147 | else { |
||
| 148 | for (int i=1; i<=numLsf; i++) { |
||
| 149 | inputFile >> pfIn; |
||
| 150 | if (pfIn > 0.0) { |
||
| 151 | failureHasOccured = true; |
||
| 152 | } |
||
| 153 | inputFile >> CovIn; |
||
| 154 | sum_q(i-1) = pfIn*k; |
||
| 155 | var_qbar = (CovIn*pfIn)*(CovIn*pfIn); |
||
| 156 | if (k < 1.0e-6) { |
||
| 157 | opserr << "WARNING: Zero number of samples read from restart file" << endln; |
||
| 158 | } |
||
| 159 | sum_q_squared(i-1) = k*( k*var_qbar + pow(sum_q(i-1)/k,2.0) ); |
||
| 160 | } |
||
| 161 | k++; |
||
| 162 | } |
||
| 163 | inputFile.close(); |
||
| 164 | } |
||
| 165 | } |
||
| 166 | |||
| 167 | |||
| 168 | // Transform start point into standard normal space, |
||
| 169 | // unless it is the origin that is to be sampled around |
||
| 170 | Vector startPointY(numRV); |
||
| 171 | |||
| 172 | if (startPoint == 0) { |
||
| 173 | // Do nothing; keep it as zero |
||
| 174 | } |
||
| 175 | else { |
||
| 176 | /* |
||
| 177 | result = theProbabilityTransformation->set_x(*startPoint); |
||
| 178 | if (result < 0) { |
||
| 179 | opserr << "ImportanceSamplingAnalysis::analyze() - could not " << endln |
||
| 180 | << " set the x-vector for xu-transformation. " << endln; |
||
| 181 | return -1; |
||
| 182 | } |
||
| 183 | result = theProbabilityTransformation->transform_x_to_u(); |
||
| 184 | if (result < 0) { |
||
| 185 | opserr << "ImportanceSamplingAnalysis::analyze() - could not " << endln |
||
| 186 | << " transform x to u. " << endln; |
||
| 187 | return -1; |
||
| 188 | } |
||
| 189 | startPointY = theProbabilityTransformation->get_u(); |
||
| 190 | */ |
||
| 191 | result = theProbabilityTransformation->transform_x_to_u(*startPoint, startPointY); |
||
| 192 | if (result < 0) { |
||
| 193 | opserr << "ImportanceSamplingAnalysis::analyze() - could not " << endln |
||
| 194 | << " transform x to u. " << endln; |
||
| 195 | return -1; |
||
| 196 | } |
||
| 197 | } |
||
| 198 | |||
| 199 | // Initial declarations |
||
| 200 | Vector cov_of_q_bar(numLsf); |
||
| 201 | Vector q_bar(numLsf); |
||
| 202 | Vector variance_of_q_bar(numLsf); |
||
| 203 | Vector responseVariance(numLsf); |
||
| 204 | Vector responseStdv(numLsf); |
||
| 205 | Vector g_storage(numLsf); |
||
| 206 | Vector sum_of_g_minus_mean_squared(numLsf); |
||
| 207 | Matrix crossSums(numLsf,numLsf); |
||
| 208 | Matrix responseCorrelation(numLsf,numLsf); |
||
| 209 | char myString[60]; |
||
| 210 | Vector pf(numLsf); |
||
| 211 | Vector cov(numLsf); |
||
| 212 | double govCov = 999.0; |
||
| 213 | //Vector temp1; |
||
| 214 | double temp2, denumerator; |
||
| 215 | bool FEconvergence; |
||
| 216 | |||
| 217 | |||
| 218 | // Prepare output file |
||
| 219 | ofstream resultsOutputFile( fileName, ios::out ); |
||
| 220 | |||
| 221 | |||
| 222 | bool isFirstSimulation = true; |
||
| 223 | while( (k<=numberOfSimulations && govCov>targetCOV || k<=2) ) { |
||
| 224 | |||
| 225 | // Keep the user posted |
||
| 226 | if (printFlag == 1 || printFlag == 2) { |
||
| 227 | opserr << "Sample #" << k << ":" << endln; |
||
| 228 | } |
||
| 229 | |||
| 230 | |||
| 231 | // Create array of standard normal random numbers |
||
| 232 | if (isFirstSimulation) { |
||
| 233 | result = theRandomNumberGenerator->generate_nIndependentStdNormalNumbers(numRV,seed); |
||
| 234 | } |
||
| 235 | else { |
||
| 236 | result = theRandomNumberGenerator->generate_nIndependentStdNormalNumbers(numRV); |
||
| 237 | } |
||
| 238 | seed = theRandomNumberGenerator->getSeed(); |
||
| 239 | if (result < 0) { |
||
| 240 | opserr << "ImportanceSamplingAnalysis::analyze() - could not generate" << endln |
||
| 241 | << " random numbers for simulation." << endln; |
||
| 242 | return -1; |
||
| 243 | } |
||
| 244 | randomArray = theRandomNumberGenerator->getGeneratedNumbers(); |
||
| 245 | |||
| 246 | // Compute the point in standard normal space |
||
| 247 | //u = startPointY + chol_covariance * randomArray; |
||
| 248 | u = startPointY; |
||
| 249 | u.addVector(1.0, randomArray, samplingStdv); |
||
| 250 | |||
| 251 | // Transform into original space |
||
| 252 | /* |
||
| 253 | result = theProbabilityTransformation->set_u(u); |
||
| 254 | if (result < 0) { |
||
| 255 | opserr << "ImportanceSamplingAnalysis::analyze() - could not set the u-vector for xu-transformation. " << endln; |
||
| 256 | return -1; |
||
| 257 | } |
||
| 258 | |||
| 259 | |||
| 260 | result = theProbabilityTransformation->transform_u_to_x(); |
||
| 261 | if (result < 0) { |
||
| 262 | opserr << "ImportanceSamplingAnalysis::analyze() - could not transform u to x. " << endln; |
||
| 263 | return -1; |
||
| 264 | } |
||
| 265 | x = theProbabilityTransformation->get_x(); |
||
| 266 | */ |
||
| 267 | result = theProbabilityTransformation->transform_u_to_x(u, x); |
||
| 268 | if (result < 0) { |
||
| 269 | opserr << "ImportanceSamplingAnalysis::analyze() - could not transform u to x. " << endln; |
||
| 270 | return -1; |
||
| 271 | } |
||
| 272 | |||
| 273 | // Evaluate limit-state function |
||
| 274 | FEconvergence = true; |
||
| 275 | result = theGFunEvaluator->runGFunAnalysis(x); |
||
| 276 | if (result < 0) { |
||
| 277 | // In this case a failure happened during the analysis |
||
| 278 | // Hence, register this as failure |
||
| 279 | FEconvergence = false; |
||
| 280 | } |
||
| 281 | |||
| 282 | |||
| 283 | LimitStateFunctionIter &lsfIter = theReliabilityDomain->getLimitStateFunctions(); |
||
| 284 | LimitStateFunction *theLimitStateFunction; |
||
| 285 | // Loop over number of limit-state functions |
||
| 286 | //for (int lsf=0; lsf<numLsf; lsf++ ) { |
||
| 287 | while ((theLimitStateFunction = lsfIter()) != 0) { |
||
| 288 | int lsf = theLimitStateFunction->getIndex(); |
||
| 289 | int lsfTag = theLimitStateFunction->getTag(); |
||
| 290 | |||
| 291 | |||
| 292 | // Set tag of "active" limit-state function |
||
| 293 | theReliabilityDomain->setTagOfActiveLimitStateFunction(lsfTag); |
||
| 294 | |||
| 295 | // set namespace variable for tcl procedures |
||
| 296 | Tcl_SetVar2Ex(interp,"RELIABILITY_lsf",NULL,Tcl_NewIntObj(lsfTag),TCL_NAMESPACE_ONLY); |
||
| 297 | |||
| 298 | // Get value of limit-state function |
||
| 299 | result = theGFunEvaluator->evaluateG(x); |
||
| 300 | if (result < 0) { |
||
| 301 | opserr << "ImportanceSamplingAnalysis::analyze() - could not tokenize limit-state function. " << endln; |
||
| 302 | return -1; |
||
| 303 | } |
||
| 304 | gFunctionValue = theGFunEvaluator->getG(); |
||
| 305 | if (!FEconvergence) { |
||
| 306 | gFunctionValue = -1.0; |
||
| 307 | } |
||
| 308 | |||
| 309 | |||
| 310 | |||
| 311 | // ESTIMATION OF FAILURE PROBABILITY |
||
| 312 | if (analysisTypeTag == 1) { |
||
| 313 | |||
| 314 | // Collect result of sampling |
||
| 315 | if (gFunctionValue < 0.0) { |
||
| 316 | I = 1; |
||
| 317 | failureHasOccured = true; |
||
| 318 | } |
||
| 319 | else { |
||
| 320 | I = 0; |
||
| 321 | } |
||
| 322 | |||
| 323 | |||
| 324 | // Compute values of joint distributions at the u-point |
||
| 325 | phi = factor1 * exp( -0.5 * (u ^ u) ); |
||
| 326 | //temp1 = inv_covariance ^ (u-startPointY); |
||
| 327 | //temp2 = temp1 ^ (u-startPointY); |
||
| 328 | temp2 = 0.0; |
||
| 329 | for (int i = 0; i < numRV; i++) { |
||
| 330 | double uy = u(i)-startPointY(i); |
||
| 331 | temp2 += uy*uy; |
||
| 332 | } |
||
| 333 | temp2 /= samplingStdv*samplingStdv; |
||
| 334 | h = factor2 * exp( -0.5 * temp2 ); |
||
| 335 | |||
| 336 | |||
| 337 | // Update sums |
||
| 338 | q = I * phi / h; |
||
| 339 | sum_q(lsf) = sum_q(lsf) + q; |
||
| 340 | sum_q_squared(lsf) = sum_q_squared(lsf) + q*q; |
||
| 341 | |||
| 342 | |||
| 343 | |||
| 344 | if (sum_q(lsf) > 0.0) { |
||
| 345 | // Compute coefficient of variation (of pf) |
||
| 346 | q_bar(lsf) = sum_q(lsf)/k; |
||
| 347 | variance_of_q_bar(lsf) = ( sum_q_squared(lsf)/k - (sum_q(lsf)/k)*(sum_q(lsf)/k)) / k; |
||
| 348 | if (variance_of_q_bar(lsf) < 0.0) |
||
| 349 | variance_of_q_bar(lsf) = 0.0; |
||
| 350 | cov_of_q_bar(lsf) = sqrt(variance_of_q_bar(lsf)) / q_bar(lsf); |
||
| 351 | } |
||
| 352 | |||
| 353 | } |
||
| 354 | else if (analysisTypeTag == 2) { |
||
| 355 | // ESTIMATION OF RESPONSE STATISTICS |
||
| 356 | |||
| 357 | // Now q=g and q_bar=mean |
||
| 358 | q = gFunctionValue; |
||
| 359 | failureHasOccured = true; |
||
| 360 | |||
| 361 | sum_q(lsf) = sum_q(lsf) + q; |
||
| 362 | sum_q_squared(lsf) = sum_q_squared(lsf) + q*q; |
||
| 363 | |||
| 364 | g_storage(lsf) = gFunctionValue; |
||
| 365 | |||
| 366 | if (sum_q(lsf) > 0.0) { |
||
| 367 | |||
| 368 | // Compute coefficient of variation (of mean) |
||
| 369 | //q_bar(lsf) = 1.0/(double)k * sum_q(lsf); |
||
| 370 | q_bar(lsf) = sum_q(lsf)/k; |
||
| 371 | //variance_of_q_bar(lsf) = 1.0/(double)k * |
||
| 372 | // ( 1.0/(double)k * sum_q_squared(lsf) - (sum_q(lsf)/(double)k)*(sum_q(lsf)/(double)k)); |
||
| 373 | variance_of_q_bar(lsf) = ( sum_q_squared(lsf)/k - (sum_q(lsf)/k)*(sum_q(lsf)/k) ) / k; |
||
| 374 | if (variance_of_q_bar(lsf) < 0.0) { |
||
| 375 | variance_of_q_bar(lsf) = 0.0; |
||
| 376 | } |
||
| 377 | cov_of_q_bar(lsf) = sqrt(variance_of_q_bar(lsf)) / q_bar(lsf); |
||
| 378 | |||
| 379 | // Compute variance and standard deviation |
||
| 380 | if (k > 1) |
||
| 381 | //responseVariance(lsf) = 1.0/((double)k-1) * ( sum_q_squared(lsf) - 1.0/((double)k) * sum_q(lsf) * sum_q(lsf) ); |
||
| 382 | responseVariance(lsf) = ( sum_q_squared(lsf) - sum_q(lsf)/k * sum_q(lsf) ) / (k-1); |
||
| 383 | else |
||
| 384 | responseVariance(lsf) = 1.0; |
||
| 385 | |||
| 386 | if (responseVariance(lsf) <= 0.0) { |
||
| 387 | opserr << "ERROR: Response variance of limit-state function number "<< lsf |
||
| 388 | << " is zero! " << endln; |
||
| 389 | } |
||
| 390 | else { |
||
| 391 | responseStdv(lsf) = sqrt(responseVariance(lsf)); |
||
| 392 | } |
||
| 393 | } |
||
| 394 | } |
||
| 395 | else if (analysisTypeTag == 3) { |
||
| 396 | // Store g-function values to file (one in each column) |
||
| 397 | //sprintf(myString,"%12.6e",gFunctionValue); |
||
| 398 | resultsOutputFile << setiosflags(ios::scientific) << setprecision(6) << gFunctionValue << " "; |
||
| 399 | resultsOutputFile.flush(); |
||
| 400 | } |
||
| 401 | else { |
||
| 402 | opserr << "ERROR: Invalid analysis type tag found in sampling analysis." << endln; |
||
| 403 | } |
||
| 404 | |||
| 405 | // Keep the user posted |
||
| 406 | if ( (printFlag == 1 || printFlag == 2) && analysisTypeTag != 3) { |
||
| 407 | sprintf(myString," GFun #%d, estimate:%15.10f, cov:%15.10f",lsfTag,q_bar(lsf),cov_of_q_bar(lsf)); |
||
| 408 | opserr << myString << endln; |
||
| 409 | } |
||
| 410 | } |
||
| 411 | |||
| 412 | // Now all the limit-state functions have been looped over |
||
| 413 | |||
| 414 | |||
| 415 | if (analysisTypeTag == 3) { |
||
| 416 | resultsOutputFile << endln; |
||
| 417 | } |
||
| 418 | |||
| 419 | |||
| 420 | // Possibly compute correlation coefficient |
||
| 421 | if (analysisTypeTag == 2) { |
||
| 422 | |||
| 423 | for (int i=0; i<numLsf; i++) { |
||
| 424 | for (int j=i+1; j<numLsf; j++) { |
||
| 425 | |||
| 426 | //crossSums(i,j) = crossSums(i,j) + g_storage(i) * g_storage(j); |
||
| 427 | crossSums(i,j) = g_storage(i) * g_storage(j); |
||
| 428 | |||
| 429 | //denumerator = (sum_q_squared(i)-1.0/(double)k*sum_q(i)*sum_q(i)) |
||
| 430 | //*(sum_q_squared(j)-1.0/(double)k*sum_q(j)*sum_q(j)); |
||
| 431 | denumerator = (sum_q_squared(i)-sum_q(i)/k*sum_q(i))*(sum_q_squared(j)-sum_q(j)/k*sum_q(j)); |
||
| 432 | |||
| 433 | if (denumerator <= 0.0) |
||
| 434 | responseCorrelation(i,j) = 0.0; |
||
| 435 | else |
||
| 436 | responseCorrelation(i,j) = (crossSums(i,j)-sum_q(i)/k*sum_q(j)) / sqrt(denumerator); |
||
| 437 | } |
||
| 438 | } |
||
| 439 | } |
||
| 440 | |||
| 441 | |||
| 442 | // Compute governing coefficient of variation |
||
| 443 | if (!failureHasOccured) { |
||
| 444 | govCov = 999.0; |
||
| 445 | } |
||
| 446 | else { |
||
| 447 | govCov = 0.0; |
||
| 448 | for (int mmmm=0; mmmm<numLsf; mmmm++) { |
||
| 449 | if (cov_of_q_bar(mmmm) > govCov) { |
||
| 450 | govCov = cov_of_q_bar(mmmm); |
||
| 451 | } |
||
| 452 | } |
||
| 453 | } |
||
| 454 | |||
| 455 | |||
| 456 | // Make sure the cov isn't exactly zero; that could be the case if only failures |
||
| 457 | // occur in cases where the 'q' remains 1 |
||
| 458 | if (govCov == 0.0) { |
||
| 459 | govCov = 999.0; |
||
| 460 | } |
||
| 461 | |||
| 462 | |||
| 463 | // Print to the restart file, if requested. |
||
| 464 | if (printFlag == 2) { |
||
| 465 | ofstream outputFile( restartFileName, ios::out ); |
||
| 466 | outputFile << k << endln; |
||
| 467 | outputFile << seed << endln; |
||
| 468 | for (int lsf=0; lsf<numLsf; lsf++ ) { |
||
| 469 | sprintf(myString,"%15.10f %15.10f",q_bar(lsf),cov_of_q_bar(lsf)); |
||
| 470 | outputFile << myString << " " << endln; |
||
| 471 | } |
||
| 472 | outputFile.close(); |
||
| 473 | } |
||
| 474 | |||
| 475 | // Increment k (the simulation number counter) |
||
| 476 | k++; |
||
| 477 | isFirstSimulation = false; |
||
| 478 | |||
| 479 | } |
||
| 480 | |||
| 481 | // Step 'k' back a step now that we went out |
||
| 482 | k--; |
||
| 483 | opserr << endln; |
||
| 484 | |||
| 485 | |||
| 486 | if (analysisTypeTag != 3) { |
||
| 487 | |||
| 488 | if (!failureHasOccured) { |
||
| 489 | opserr << "WARNING: Failure did not occur for any of the limit-state functions. " << endln; |
||
| 490 | } |
||
| 491 | |||
| 492 | LimitStateFunctionIter &lsfIter = theReliabilityDomain->getLimitStateFunctions(); |
||
| 493 | LimitStateFunction *theLimitStateFunction; |
||
| 494 | //for (int lsf=1; lsf<=numLsf; lsf++ ) { |
||
| 495 | while ((theLimitStateFunction = lsfIter()) != 0) { |
||
| 496 | int lsfIndex = theLimitStateFunction->getIndex(); |
||
| 497 | int lsf = theLimitStateFunction->getTag(); |
||
| 498 | |||
| 499 | if ( q_bar(lsfIndex) == 0.0 ) { |
||
| 500 | |||
| 501 | resultsOutputFile << "#######################################################################" << endln; |
||
| 502 | resultsOutputFile << "# SAMPLING ANALYSIS RESULTS, LIMIT-STATEFUNCDTION NUMBER " |
||
| 503 | <<setiosflags(ios::left)<<setprecision(1)<<setw(4)<<lsf <<" #" << endln; |
||
| 504 | resultsOutputFile << "# #" << endln; |
||
| 505 | resultsOutputFile << "# Failure did not occur, or zero response! #" << endln; |
||
| 506 | resultsOutputFile << "# #" << endln; |
||
| 507 | resultsOutputFile << "#######################################################################" << endln << endln << endln; |
||
| 508 | } |
||
| 509 | else { |
||
| 510 | |||
| 511 | // Some declarations |
||
| 512 | double beta_sim, pf_sim, cov_sim; |
||
| 513 | int num_sim; |
||
| 514 | |||
| 515 | |||
| 516 | // Set tag of "active" limit-state function |
||
| 517 | theReliabilityDomain->setTagOfActiveLimitStateFunction(lsf); |
||
| 518 | |||
| 519 | |||
| 520 | // Store results |
||
| 521 | if (analysisTypeTag == 1) { |
||
| 522 | beta_sim = -aStdNormRV.getInverseCDFvalue(q_bar(lsfIndex)); |
||
| 523 | pf_sim = q_bar(lsfIndex); |
||
| 524 | cov_sim = cov_of_q_bar(lsfIndex); |
||
| 525 | num_sim = k; |
||
| 526 | theLimitStateFunction->SimulationReliabilityIndexBeta = beta_sim; |
||
| 527 | theLimitStateFunction->SimulationProbabilityOfFailure_pfsim = pf_sim; |
||
| 528 | theLimitStateFunction->CoefficientOfVariationOfPfFromSimulation = cov_sim; |
||
| 529 | theLimitStateFunction->NumberOfSimulations = num_sim; |
||
| 530 | } |
||
| 531 | |||
| 532 | |||
| 533 | // Print results to the output file |
||
| 534 | if (analysisTypeTag == 1) { |
||
| 535 | resultsOutputFile << "#######################################################################" << endln; |
||
| 536 | resultsOutputFile << "# SAMPLING ANALYSIS RESULTS, LIMIT-STATE FUNCTION NUMBER " |
||
| 537 | <<setiosflags(ios::left)<<setprecision(1)<<setw(4)<<lsf <<" #" << endln; |
||
| 538 | resultsOutputFile << "# #" << endln; |
||
| 539 | resultsOutputFile << "# Reliability index beta: ............................ " |
||
| 540 | <<setiosflags(ios::left)<<setprecision(5)<<setw(12)<<beta_sim |
||
| 541 | << " #" << endln; |
||
| 542 | resultsOutputFile << "# Estimated probability of failure pf_sim: ........... " |
||
| 543 | <<setiosflags(ios::left)<<setprecision(5)<<setw(12)<<pf_sim |
||
| 544 | << " #" << endln; |
||
| 545 | resultsOutputFile << "# Number of simulations: ............................. " |
||
| 546 | <<setiosflags(ios::left)<<setprecision(5)<<setw(12)<<num_sim |
||
| 547 | << " #" << endln; |
||
| 548 | resultsOutputFile << "# Coefficient of variation (of pf): .................. " |
||
| 549 | <<setiosflags(ios::left)<<setprecision(5)<<setw(12)<<cov_sim |
||
| 550 | << " #" << endln; |
||
| 551 | resultsOutputFile << "# #" << endln; |
||
| 552 | resultsOutputFile << "#######################################################################" << endln << endln << endln; |
||
| 553 | } |
||
| 554 | else { |
||
| 555 | resultsOutputFile << "#######################################################################" << endln; |
||
| 556 | resultsOutputFile << "# SAMPLING ANALYSIS RESULTS, LIMIT-STATE FUNCTION NUMBER " |
||
| 557 | <<setiosflags(ios::left)<<setprecision(1)<<setw(4)<<lsf <<" #" << endln; |
||
| 558 | resultsOutputFile << "# #" << endln; |
||
| 559 | resultsOutputFile << "# Estimated mean: .................................... " |
||
| 560 | <<setiosflags(ios::left)<<setprecision(5)<<setw(12)<<q_bar(lsfIndex) |
||
| 561 | << " #" << endln; |
||
| 562 | resultsOutputFile << "# Estimated standard deviation: ...................... " |
||
| 563 | <<setiosflags(ios::left)<<setprecision(5)<<setw(12)<<responseStdv(lsfIndex) |
||
| 564 | << " #" << endln; |
||
| 565 | resultsOutputFile << "# #" << endln; |
||
| 566 | resultsOutputFile << "#######################################################################" << endln << endln << endln; |
||
| 567 | } |
||
| 568 | } |
||
| 569 | } |
||
| 570 | |||
| 571 | if (analysisTypeTag == 2) { |
||
| 572 | resultsOutputFile << "#######################################################################" << endln; |
||
| 573 | resultsOutputFile << "# RESPONSE CORRELATION COEFFICIENTS #" << endln; |
||
| 574 | resultsOutputFile << "# #" << endln; |
||
| 575 | if (numLsf <=1) { |
||
| 576 | resultsOutputFile << "# Only one limit-state function! #" << endln; |
||
| 577 | } |
||
| 578 | else { |
||
| 579 | resultsOutputFile << "# gFun gFun Correlation #" << endln; |
||
| 580 | resultsOutputFile.setf(ios::fixed, ios::floatfield); |
||
| 581 | LimitStateFunctionIter lsfIterI = theReliabilityDomain->getLimitStateFunctions(); |
||
| 582 | LimitStateFunctionIter lsfIterJ = theReliabilityDomain->getLimitStateFunctions(); |
||
| 583 | LimitStateFunction *lsfI, *lsfJ; |
||
| 584 | //for (i=0; i<numLsf; i++) { |
||
| 585 | while ((lsfI = lsfIterI()) != 0) { |
||
| 586 | int iTag = lsfI->getTag(); |
||
| 587 | int i = lsfI->getIndex(); |
||
| 588 | lsfIterJ.reset(); |
||
| 589 | //for (int j=i+1; j<numLsf; j++) { |
||
| 590 | while ((lsfJ = lsfIterJ()) != 0) { |
||
| 591 | int jTag = lsfJ->getTag(); |
||
| 592 | int j = lsfJ->getIndex(); |
||
| 593 | resultsOutputFile << "# " <<setw(3)<< iTag <<" "<<setw(3)<< jTag <<" "; |
||
| 594 | if (responseCorrelation(i,j)<0.0) { resultsOutputFile << "-"; } |
||
| 595 | else { resultsOutputFile << " "; } |
||
| 596 | resultsOutputFile <<setprecision(7)<<setw(11)<<fabs(responseCorrelation(i,j)); |
||
| 597 | resultsOutputFile << " #" << endln; |
||
| 598 | } |
||
| 599 | } |
||
| 600 | } |
||
| 601 | resultsOutputFile << "# #" << endln; |
||
| 602 | resultsOutputFile << "#######################################################################" << endln << endln << endln; |
||
| 603 | } |
||
| 604 | |||
| 605 | } |
||
| 606 | |||
| 607 | |||
| 608 | |||
| 609 | |||
| 610 | |||
| 611 | // Print summary of results to screen |
||
| 612 | opserr << "Simulation Analysis completed." << endln; |
||
| 613 | |||
| 614 | // Clean up |
||
| 615 | resultsOutputFile.close(); |
||
| 616 | |||
| 617 | return 0; |
||
| 618 | } |
||
| 619 |