Rev 4875 | Rev 4879 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed
| Rev 4875 | Rev 4877 | ||
|---|---|---|---|
| Line 365... | Line 365... | ||
| 365 | Vector z(nrv); |
365 | Vector z(nrv); |
| 366 | this->x_to_z(z); |
366 | this->x_to_z(z); |
| 367 | 367 | ||
| 368 | // 3) DzDmean and DzDstdv = a vector of zeros and then:
|
368 | // 3) DzDmean and DzDstdv = a vector of zeros and then:
|
| 369 | double DzDmean = 0.0; |
369 | double DzDmean = 0.0; |
| 370 | static NormalRV aStandardNormalRV(1,0.0,1.0); |
- | |
| 371 | RandomVariable *theRV = theReliabilityDomain->getRandomVariablePtr(rvTag); |
370 | RandomVariable *theRV = theReliabilityDomain->getRandomVariablePtr(rvTag); |
| 372 | if (theRV == 0) { |
371 | if (theRV == 0) { |
| 373 | opserr << "NatafProbTransf::meanSensitivityOf_x_to_u -- r.v. with tag " << rvTag |
372 | opserr << "NatafProbTransf::meanSensitivityOf_x_to_u -- r.v. with tag " << rvTag |
| 374 | << " not found in reliability domain" << endln; |
373 | << " not found in reliability domain" << endln; |
| 375 | }
|
374 | }
|
| 376 | 375 | ||
| 377 | int rvIndex = theReliabilityDomain->getRandomVariableIndex(rvTag); |
376 | int rvIndex = theReliabilityDomain->getRandomVariableIndex(rvTag); |
| 378 | - | ||
| 379 | if (strcmp(theRV->getType(),"NORMAL")==0) { |
- | |
| 380 | double sigma = theRV->getStdv(); |
- | |
| 381 | DzDmean = -1.0 / sigma; |
- | |
| 382 | }
|
- | |
| 383 | else if (strcmp(theRV->getType(),"LOGNORMAL")==0) { |
- | |
| 384 | double mean = fabs(theRV->getMean()); // more here for negative lognormal? |
- | |
| 385 | double stdv = theRV->getStdv(); |
- | |
| 386 | - | ||
| 387 | double a = mean*mean+stdv*stdv; |
- | |
| 388 | DzDmean = 0.5*(-2.0*mean*mean*log(a)+4.0*mean*mean*log(mean)-3.0*stdv*stdv*log(a)+4.0*stdv*stdv*log(mean)+2.0*stdv*stdv*log(fabs(x(rvIndex)))) |
- | |
| 389 | /(pow(log(a)-2.0*log(mean),1.5)*mean*a); |
- | |
| 390 | - | ||
| 391 | // double z = ( log ( fabs(x(rvIndex)) ) - lambda ) / zeta; // more here for negative lognormal?
|
- | |
| 392 | // double e = (stdv/mean)*(stdv/mean);
|
- | |
| 393 | // double d = 1 + e;
|
- | |
| 394 | // double f = 1.0 / (mean*d*zeta);
|
- | |
| 395 | // DzDmean = f*(-d-e+e*z/zeta);
|
- | |
| 396 | }
|
- | |
| 397 | else if (strcmp(theRV->getType(),"UNIFORM")==0) { |
- | |
| 398 | double pz = 0.39894228048*exp(-z(rvIndex)*z(rvIndex)/2.0); |
- | |
| 399 | Vector paramTemp = theRV->getParameters(); |
- | |
| 400 | double a = paramTemp(0); |
- | |
| 401 | double b = paramTemp(1); |
- | |
| 402 | DzDmean = -1.0/(pz*(b-a)); |
- | |
| 403 | }
|
- | |
| 404 | else { |
- | |
| 405 | opserr << "WARNING: Cannot compute reliability sensitivity results for " << endln |
- | |
| 406 | << " type of random variable number " << rvTag << endln; |
- | |
| 407 | //DzDmean = 0.0;
|
- | |
| 408 | //DzDmean = aStandardNormalRV.getInverseCDFvalue(theRV->getCDFMeanSensitivity(x(rvIndex)));
|
- | |
| 409 | DzDmean = aStandardNormalRV.getInverseCDFvalue(theRV->getCDFMeanSensitivity()); |
- | |
| 410 | }
|
- | |
| - | 377 | DzDmean = theRV->getCDFMeanSensitivity(); |
|
| - | 378 | opserr << "DzDmean is " << DzDmean << endln; |
|
| 411 | 379 | ||
| 412 | // 4) The hardest part: DinverseLowerCholeskyDmean
|
380 | // 4) The hardest part: DinverseLowerCholeskyDmean
|
| 413 | // DinverseLowerCholeskyDstdv
|
381 | // DinverseLowerCholeskyDstdv
|
| 414 | // This is possible, but a bit tedious
|
382 | // This is possible, but a bit tedious
|
| 415 | // The procedure would be:
|
383 | // The procedure would be:
|
| Line 427... | Line 395... | ||
| 427 | OrigLowerCholesky(i,j) = lapackA[jnrv+i]; |
395 | OrigLowerCholesky(i,j) = lapackA[jnrv+i]; |
| 428 | }
|
396 | }
|
| 429 | Matrix OrigInverseLowerCholesky(nrv,nrv); |
397 | Matrix OrigInverseLowerCholesky(nrv,nrv); |
| 430 | OrigLowerCholesky.Invert(OrigInverseLowerCholesky); |
398 | OrigLowerCholesky.Invert(OrigInverseLowerCholesky); |
| 431 | 399 | ||
| 432 | RandomVariable *aRandomVariable; |
- | |
| 433 | aRandomVariable = theReliabilityDomain->getRandomVariablePtr(rvTag); |
- | |
| 434 | double stdv = aRandomVariable->getStdv(); |
- | |
| - | 400 | double stdv = theRV->getStdv(); |
|
| 435 | double h = stdv/200.0; |
401 | double h = stdv/200.0; |
| 436 | setCorrelationMatrix(rvTag, 0, h); // rvTag or rvIndex ... does it matter? -- MHS |
- | |
| - | 402 | setCorrelationMatrix(rvTag, 0, h); |
|
| 437 | 403 | ||
| 438 | MatrixOperations someMatrixOperations(*correlationMatrix); |
404 | MatrixOperations someMatrixOperations(*correlationMatrix); |
| 439 | - | ||
| 440 | int result = someMatrixOperations.computeCholeskyAndItsInverse(); |
405 | int result = someMatrixOperations.computeCholeskyAndItsInverse(); |
| 441 | if (result < 0) { |
406 | if (result < 0) { |
| 442 | opserr << "NatafProbabilityTransformation::NatafProbabilityTransformation() - could not" << endln |
407 | opserr << "NatafProbabilityTransformation::NatafProbabilityTransformation() - could not" << endln |
| 443 | << " compute the Cholesky decomposition and its inverse " << endln |
408 | << " compute the Cholesky decomposition and its inverse " << endln |
| 444 | << " for the correlation matrix." << endln; |
409 | << " for the correlation matrix." << endln; |
| Line 449... | Line 414... | ||
| 449 | DinverseLowerCholeskyDmean.addMatrix(1.0, PerturbedInverseLowerCholesky, -1/h); |
414 | DinverseLowerCholeskyDmean.addMatrix(1.0, PerturbedInverseLowerCholesky, -1/h); |
| 450 | setCorrelationMatrix(0, 0, 0.0); |
415 | setCorrelationMatrix(0, 0, 0.0); |
| 451 | 416 | ||
| 452 | // Return the final result (the four factors)
|
417 | // Return the final result (the four factors)
|
| 453 | //return ( DinverseLowerCholeskyDmean * z + OrigInverseLowerCholesky * DzDmean );
|
418 | //return ( DinverseLowerCholeskyDmean * z + OrigInverseLowerCholesky * DzDmean );
|
| 454 | - | ||
| 455 | Vector returnVector(z); |
419 | Vector returnVector(z); |
| 456 | returnVector.addMatrixVector(0.0, DinverseLowerCholeskyDmean, z, 1.0); |
420 | returnVector.addMatrixVector(0.0, DinverseLowerCholeskyDmean, z, 1.0); |
| 457 | //returnVector.addMatrixVector(1.0, OrigInverseLowerCholesky, DzDmean, 1.0);
|
421 | //returnVector.addMatrixVector(1.0, OrigInverseLowerCholesky, DzDmean, 1.0);
|
| 458 | // Unroll for lower triangular
|
422 | // Unroll for lower triangular
|
| 459 | for (int i = rvIndex; i < nrv; i++) { |
423 | for (int i = rvIndex; i < nrv; i++) { |
| Line 479... | Line 443... | ||
| 479 | Vector z(nrv); |
443 | Vector z(nrv); |
| 480 | this->x_to_z(z); |
444 | this->x_to_z(z); |
| 481 | 445 | ||
| 482 | // 3) DzDmean and DzDstdv = a vector of zeros and then:
|
446 | // 3) DzDmean and DzDstdv = a vector of zeros and then:
|
| 483 | double DzDstdv = 0.0; |
447 | double DzDstdv = 0.0; |
| 484 | static NormalRV aStandardNormalRV(1,0.0,1.0); |
- | |
| 485 | RandomVariable *theRV = theReliabilityDomain->getRandomVariablePtr(rvTag); |
448 | RandomVariable *theRV = theReliabilityDomain->getRandomVariablePtr(rvTag); |
| 486 | if (theRV == 0) { |
449 | if (theRV == 0) { |
| 487 | opserr << "NatafProbTransf::stdvSensitivityOf_x_to_u -- r.v. with tag " << rvTag |
450 | opserr << "NatafProbTransf::stdvSensitivityOf_x_to_u -- r.v. with tag " << rvTag |
| 488 | << " not found in reliability domain" << endln; |
451 | << " not found in reliability domain" << endln; |
| 489 | }
|
452 | }
|
| 490 | 453 | ||
| 491 | int rvIndex = theReliabilityDomain->getRandomVariableIndex(rvTag); |
454 | int rvIndex = theReliabilityDomain->getRandomVariableIndex(rvTag); |
| 492 | - | ||
| 493 | if (strcmp(theRV->getType(),"NORMAL")==0) { |
- | |
| 494 | double mu = theRV->getMean(); |
- | |
| 495 | double sigma = theRV->getStdv(); |
- | |
| 496 | DzDstdv = - (x(rvIndex)-mu) / (sigma*sigma); |
- | |
| 497 | }
|
- | |
| 498 | else if (strcmp(theRV->getType(),"LOGNORMAL")==0) { |
- | |
| 499 | double mean = fabs(theRV->getMean()); // more here for negative lognormal? |
- | |
| 500 | double stdv = theRV->getStdv(); |
- | |
| 501 | - | ||
| 502 | double a = mean*mean+stdv*stdv; |
- | |
| 503 | DzDstdv = 0.5*stdv*(log(a)-2.0*log(fabs(x(rvIndex))))/(pow(log(a)-2.0*log(mean),1.5)*a); |
- | |
| 504 | - | ||
| 505 | // double z = ( log ( fabs(x(rvIndex)) ) - lambda ) / zeta; // more here for negative lognormal?
|
- | |
| 506 | // double e = (stdv/mean)*(stdv/mean);
|
- | |
| 507 | // double d = 1 + e;
|
- | |
| 508 | // double f = 1.0 / (mean*d*zeta);
|
- | |
| 509 | // DzDstdv = stdv*((1.0-z/zeta)*f/mean);
|
- | |
| 510 | }
|
- | |
| 511 | else if (strcmp(theRV->getType(),"UNIFORM")==0) { |
- | |
| 512 | double pz = 0.39894228048*exp(-z(rvIndex)*z(rvIndex)/2.0); |
- | |
| 513 | Vector paramTemp = theRV->getParameters(); |
- | |
| 514 | double a = paramTemp(0); |
- | |
| 515 | double b = paramTemp(1); |
- | |
| 516 | double DzDmean = -1.0/(pz*(b-a)); |
- | |
| 517 | double e = -DzDmean/(b-a); |
- | |
| 518 | DzDstdv = 1.732050807*(a+b-2.0*x(rvIndex))*e; |
- | |
| 519 | }
|
- | |
| 520 | else { |
- | |
| 521 | opserr << "WARNING: Cannot compute reliability sensitivity results for " << endln |
- | |
| 522 | << " type of random variable number " << rvTag << endln; |
- | |
| 523 | //DzDstdv = 0.0;
|
- | |
| 524 | //DzDstdv = aStandardNormalRV.getInverseCDFvalue(theRV->getCDFStdvSensitivity(x(rvIndex)));
|
- | |
| 525 | DzDstdv = aStandardNormalRV.getInverseCDFvalue(theRV->getCDFStdvSensitivity()); |
- | |
| 526 | }
|
- | |
| - | 455 | DzDstdv = theRV->getCDFStdvSensitivity(); |
|
| - | 456 | opserr << "DzDstdv is " << DzDstdv << endln; |
|
| 527 | 457 | ||
| 528 | // 4) The hardest part: DinverseLowerCholeskyDmean
|
458 | // 4) The hardest part: DinverseLowerCholeskyDmean
|
| 529 | // DinverseLowerCholeskyDstdv
|
459 | // DinverseLowerCholeskyDstdv
|
| 530 | // This is possible, but a bit tedious
|
460 | // This is possible, but a bit tedious
|
| 531 | // The procedure would be:
|
461 | // The procedure would be:
|
| Line 543... | Line 473... | ||
| 543 | OrigLowerCholesky(i,j) = lapackA[jnrv+i]; |
473 | OrigLowerCholesky(i,j) = lapackA[jnrv+i]; |
| 544 | }
|
474 | }
|
| 545 | Matrix OrigInverseLowerCholesky(nrv,nrv); |
475 | Matrix OrigInverseLowerCholesky(nrv,nrv); |
| 546 | OrigLowerCholesky.Invert(OrigInverseLowerCholesky); |
476 | OrigLowerCholesky.Invert(OrigInverseLowerCholesky); |
| 547 | 477 | ||
| 548 | RandomVariable *aRandomVariable; |
- | |
| 549 | aRandomVariable = theReliabilityDomain->getRandomVariablePtr(rvTag); |
- | |
| 550 | double stdv = aRandomVariable->getStdv(); |
- | |
| - | 478 | double stdv = theRV->getStdv(); |
|
| 551 | double h = stdv/200.0; |
479 | double h = stdv/200.0; |
| 552 | setCorrelationMatrix(0, rvTag, h); // rvTag or rvIndex ... does it matter? -- MHS |
- | |
| - | 480 | setCorrelationMatrix(0, rvTag, h); |
|
| 553 | 481 | ||
| 554 | MatrixOperations someMatrixOperations(*correlationMatrix); |
482 | MatrixOperations someMatrixOperations(*correlationMatrix); |
| 555 | - | ||
| 556 | int result = someMatrixOperations.computeCholeskyAndItsInverse(); |
483 | int result = someMatrixOperations.computeCholeskyAndItsInverse(); |
| 557 | if (result < 0) { |
484 | if (result < 0) { |
| 558 | opserr << "NatafProbabilityTransformation::NatafProbabilityTransformation() - could not" << endln |
485 | opserr << "NatafProbabilityTransformation::NatafProbabilityTransformation() - could not" << endln |
| 559 | << " compute the Cholesky decomposition and its inverse " << endln |
486 | << " compute the Cholesky decomposition and its inverse " << endln |
| 560 | << " for the correlation matrix." << endln; |
487 | << " for the correlation matrix." << endln; |
| Line 593... | Line 520... | ||
| 593 | // Initialize correlation matrix
|
520 | // Initialize correlation matrix
|
| 594 | correlationMatrix->Zero(); |
521 | correlationMatrix->Zero(); |
| 595 | 522 | ||
| 596 | // Put 'ones' on the diagonal
|
523 | // Put 'ones' on the diagonal
|
| 597 | for ( int i=0 ; i<nrv ; i++ ) |
524 | for ( int i=0 ; i<nrv ; i++ ) |
| 598 | {
|
- | |
| 599 | (*correlationMatrix)(i,i) = 1.0; |
525 | (*correlationMatrix)(i,i) = 1.0; |
| 600 | }
|
- | |
| 601 | 526 | ||
| 602 | // Get number of correlation coefficients
|
527 | // Get number of correlation coefficients
|
| 603 | int numberOfCorrelationCoefficients = |
528 | int numberOfCorrelationCoefficients = |
| 604 | theReliabilityDomain->getNumberOfCorrelationCoefficients(); |
529 | theReliabilityDomain->getNumberOfCorrelationCoefficients(); |
| 605 | 530 | ||