Aster/Cata/cataSTA9/Macro/post_dyna_alea_ops.py

   1 #@ MODIF post_dyna_alea_ops Macro  DATE 25/04/2008   AUTEUR ZENTNER I.ZENTNER
   2 # -*- coding: iso-8859-1 -*-
   3 #            CONFIGURATION MANAGEMENT OF EDF VERSION
   4 # ======================================================================
   5 # COPYRIGHT (C) 1991 - 2006  EDF R&D                  WWW.CODE-ASTER.ORG
   6 # THIS PROGRAM IS FREE SOFTWARE; YOU CAN REDISTRIBUTE IT AND/OR MODIFY
   7 # IT UNDER THE TERMS OF THE GNU GENERAL PUBLIC LICENSE AS PUBLISHED BY
   8 # THE FREE SOFTWARE FOUNDATION; EITHER VERSION 2 OF THE LICENSE, OR
   9 # (AT YOUR OPTION) ANY LATER VERSION.
  10 #
  11 # THIS PROGRAM IS DISTRIBUTED IN THE HOPE THAT IT WILL BE USEFUL, BUT
  12 # WITHOUT ANY WARRANTY; WITHOUT EVEN THE IMPLIED WARRANTY OF
  13 # MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. SEE THE GNU
  14 # GENERAL PUBLIC LICENSE FOR MORE DETAILS.
  15 #
  16 # YOU SHOULD HAVE RECEIVED A COPY OF THE GNU GENERAL PUBLIC LICENSE
  17 # ALONG WITH THIS PROGRAM; IF NOT, WRITE TO EDF R&D CODE_ASTER,
  18 #    1 AVENUE DU GENERAL DE GAULLE, 92141 CLAMART CEDEX, FRANCE.
  19 # ======================================================================
  20
  21 import random
  22 import Numeric
  23 from sets  import Set
  24 from types import ListType, TupleType
  25 from math  import pi,sqrt,log,exp
  26
  27 EnumTypes = (ListType, TupleType)
  28
  29
  30 def post_dyna_alea_ops(self,INTE_SPEC, FRAGILITE,TITRE,INFO,**args):
  31    import aster
  32    from Accas               import _F
  33    from Utilitai.Utmess     import UTMESS
  34    from Utilitai.t_fonction import t_fonction
  35    from Utilitai.Table      import Table
  36
  37    commande='POST_DYNA_ALEA'
  38
  39    ier = 0
  40    # La macro compte pour 1 dans la numerotation des commandes
  41    self.set_icmd(1)
  42
  43    # Le concept sortant (de type table_sdaster ou dérivé) est tab
  44    self.DeclareOut('tabout', self.sd)
  45
  46    # On importe les definitions des commandes a utiliser dans la macro
  47    # Le nom de la variable doit etre obligatoirement le nom de la commande
  48    CREA_TABLE    = self.get_cmd('CREA_TABLE')
  49    CALC_TABLE    = self.get_cmd('CALC_TABLE')
  50    IMPR_TABLE    = self.get_cmd('IMPR_TABLE')
  51    RECU_FONCTION = self.get_cmd('RECU_FONCTION')
  52    IMPR_FONCTION = self.get_cmd('IMPR_FONCTION')
  53    DEFI_LIST_REEL = self.get_cmd('DEFI_LIST_REEL')
  54    DEFI_FONCTION  = self.get_cmd('DEFI_FONCTION')
  55    CALC_FONCTION  = self.get_cmd('CALC_FONCTION')
  56
  57
  58 #  ------------------------------------------------------------------
  59 #---------algorithme d'optimisation pour le  maximum de vraisemblance
  60    def vrais(x):
  61       am=x[0]
  62       assert am >0.000, 'optimize.py: am negatif'
  63       beta=x[1]
  64       assert am >0.000, 'optimize.py: beta negatif'
  65       res=1.0
  66       for k in range(Nbval):
  67          ai=liste_indic[k]
  68          xi=float(liste_def[k])
  69          val=log(ai/am)
  70          pfa=normcdf(val/beta)
  71          f0=pfa**xi*(1.-pfa)**(1-xi)
  72          res=res*f0
  73       return -res
  74
  75    def boot_vrais(x):
  76       am=x[0]
  77       beta=x[1]
  78       res=1.0
  79       for k in range(Nbval):
  80          ai=liste_indic[list_rand[k]]
  81          xi=float(liste_def[list_rand[k]])
  82          val=log(ai/am)
  83          pfa=normcdf(val/beta)
  84          f0=pfa**xi*(1.-pfa)**(1-xi)
  85          res=res*f0
  86       return -res
  87
  88 #  ------------------------------------------------------------------
  89 #  OPTION FRAGILITE
  90 # ------------------------------------------------------------------
  91    if FRAGILITE !=None :
  92       from Utilitai.optimize   import fmin
  93       from Utilitai.stats   import normcdf
  94
  95       if FRAGILITE['LIST_PARA'] != None :
  96          liste_a = FRAGILITE['LIST_PARA'].VALE.get()
  97       elif FRAGILITE['VALE'] != None :
  98          liste_a =FRAGILITE['VALE']
  99
 100
 101       Nba=len(liste_a)
 102       lpfa=[]
 103       tab2 = FRAGILITE['TABL_RESU'].EXTR_TABLE()
 104       dicta = tab2.values()
 105
 106       if dicta.has_key('DEFA') :
 107          liste_def = dicta['DEFA']
 108       else:
 109          UTMESS('F','TABLE0_1',valk=('DEFA'))
 110       if dicta.has_key('PARA_NOCI') :
 111         liste_indic = dicta['PARA_NOCI']
 112       else:
 113         UTMESS('F','TABLE0_1',valk=('PARA_NOCI'))
 114
 115       Nbval=len(liste_indic)
 116
 117       test1 = Numeric.equal(None,liste_indic)
 118       test2 = Numeric.equal(None,liste_def)
 119       if test1 >=1 or test2 >=1:
 120          UTMESS('F','TABLE0_15')
 121
 122       # estimation paramètres
 123       x0 = [FRAGILITE['AM_INI'],FRAGILITE['BETA_INI']]
 124       xopt = fmin(vrais,x0)
 125
 126       texte='PARAMETRES Am, beta ESTIMES : '+str(xopt)+'\n'
 127       aster.affiche('MESSAGE',texte)      #print 'parametres Am, beta estimes: ', xopt
 128
 129       #courbe de fragilité
 130       vec_a=Numeric.array(liste_a)
 131       vecval=(Numeric.log(vec_a/xopt[0]))/xopt[1]
 132       for m in range(Nba):
 133          lpfa.append(normcdf(vecval[m]))
 134
 135       # table sortie
 136
 137       mcfact=[]
 138       if  TITRE !=None :
 139            mcfact.append(_F(PARA= 'TITRE' , LISTE_K= TITRE  ))
 140
 141       mcfact.append(_F(PARA= 'AM' ,LISTE_R=xopt[0] ))
 142       mcfact.append(_F(PARA= 'BETA' ,LISTE_R=xopt[1] ))
 143       mcfact.append(_F(PARA= 'PARA_NOCI' ,LISTE_R =liste_a  ))
 144       mcfact.append(_F(PARA= 'PFA' ,LISTE_R = lpfa ))
 145
 146    #print 'fractiles a calculer par bootstrap : ', FRAGILITE['FRACTILE']
 147
 148
 149       # si calcul de fractiles (intervalles de confiance) par bootstrap
 150
 151       if FRAGILITE['FRACTILE']!= None :
 152          if INFO==2 :
 153             texte='FRACTILES A CALCULER PAR BOOTSTRAP '+ str(FRAGILITE['FRACTILE']) +'\n'
 154             aster.affiche('MESSAGE',texte)
 155          if FRAGILITE['NB_TIRAGE']!= None :
 156             Nboot = FRAGILITE['NB_TIRAGE']
 157             if Nboot > Nbval :
 158                UTMESS('F','PROBA0_11')     #assert Nboot <= Nbval , 'ERREUR: nombre de tirages demandes trop grand'
 159          else:
 160             Nboot = Nbval
 161
 162          list_fonc = []
 163          lfract =FRAGILITE['FRACTILE']
 164          __F1=[None]*Nbval
 165          __ABS=[None]*Nbval
 166          __ORDO=[None]*Nbval
 167
 168          for kb in range(Nboot) : #in range(Nbval)
 169
 170             lpfa = []
 171             list_rand = []
 172
 173             for kb2 in range(Nbval) :
 174                list_rand.append(random.randint(0,Nbval-1))
 175
 176             xopt = fmin(boot_vrais,x0)
 177             if INFO==2 :
 178                texte1='BOOTSTRAP TIRAGE '+ str(kb+1)
 179                texte2='  PARAMETRES Am, beta ESTIMES : '+str(xopt)+'\n'
 180                aster.affiche('MESSAGE',texte1) #print 'bootstrap tirage', kb+1, ', -  parametres Am, beta estimes: ', xopt
 181                aster.affiche('MESSAGE',texte2)
 182             vecval=(Numeric.log(vec_a/xopt[0]))/xopt[1]
 183             for m in range(Nba):
 184                lpfa.append(normcdf(vecval[m]))
 185
 186             __ABS[kb]=DEFI_LIST_REEL( VALE = liste_a  );
 187             __ORDO[kb]=DEFI_LIST_REEL( VALE = lpfa );
 188
 189             __F1[kb]=DEFI_FONCTION(  NOM_PARA='PGAZ',
 190                                      NOM_RESU = 'PFA',
 191                                      VALE_PARA = __ABS[kb],
 192                                      VALE_FONC = __ORDO[kb],);
 193             list_fonc.append(__F1[kb],)
 194
 195
 196          #__FRACTILE = [None]*len(lfract)
 197          liste = [None]*len(lfract)
 198          for kb in range(len(lfract)):
 199             __FRACTILE=CALC_FONCTION(FRACTILE=_F(FONCTION=(list_fonc),
 200                              FRACT=lfract[kb]), );
 201             liste[kb]= __FRACTILE.Ordo()
 202             mcfact.append(_F(PARA= str(lfract[kb]) ,LISTE_R =liste[kb]  ))
 203
 204
 205       #   fin FRAGILITE
 206       tabout = CREA_TABLE(LISTE=mcfact,TITRE = 'POST_DYNA_ALEA concept : '+self.sd.nom)
 207
 208 #  ------------------------------------------------------------------
 209
 210
 211 #  ------------------------------------------------------------------
 212 #  OPTION INTESPEC
 213 # ------------------------------------------------------------------
 214    if INTE_SPEC !=None :
 215
 216       TOUT_ORDRE = args['TOUT_ORDRE']
 217       NUME_VITE_FLUI=args['NUME_VITE_FLUI']
 218
 219       NUME_ORDRE_I=args['NUME_ORDRE_I']
 220       NOEUD_I=args['NOEUD_I']
 221       OPTION=args['OPTION']
 222
 223       MOMENT=args['MOMENT']
 224
 225       intespec=INTE_SPEC.EXTR_TABLE()
 226       # pour la clarté !
 227       NUME_VITE_FLUI_present = 'NUME_VITE_FLUI' in intespec.para
 228       NUME_ORDRE_I_present   = 'NUME_ORDRE_I'   in intespec.para
 229       NOEUD_I_present        = 'NOEUD_I'        in intespec.para
 230
 231       # table résultat
 232       tabres = Table(titr='POST_DYNA_ALEA concept : %s' % self.sd.nom)
 233
 234 #     ------------------------------------------------------------------
 235 #     Liste des moments spectraux
 236 #     repérer le type de l'interspectre et son nom
 237 #                   1- concept interspectre
 238 #                   2- table de table d interspectre
 239
 240       if NUME_VITE_FLUI_present :
 241          if TOUT_ORDRE != None :
 242             jvite = list(Set(intespec.NUME_VITE_FLUI.not_none_values()))
 243             jvite.sort()
 244          else :
 245             jvite=[NUME_VITE_FLUI,]
 246       else :
 247          jvite  =[None]
 248
 249 #     ------------------------------------------------------------------
 250 #     Repérer les couples d'indices selectionnés
 251 #     vérification de l'égalité du nombre d indices en i et j
 252
 253       if NUME_ORDRE_I!=None :
 254         l_ind_i=NUME_ORDRE_I
 255         l_ind_j=args['NUME_ORDRE_J']
 256         if type(l_ind_i) not in EnumTypes : l_ind_i=[l_ind_i]
 257         if type(l_ind_j) not in EnumTypes : l_ind_j=[l_ind_j]
 258         if len(l_ind_i)!=len(l_ind_j) :
 259            UTMESS('F','PROBA0_8')
 260         # paramètres fixes de la table
 261         tabres.add_para(['NUME_ORDRE_I','NUME_ORDRE_J'], 'I')
 262       elif NOEUD_I!=None :
 263         l_ind_i=NOEUD_I
 264         l_ind_j=args['NOEUD_J']
 265         l_cmp_i=args['NOM_CMP_I']
 266         l_cmp_j=args['NOM_CMP_J']
 267         if type(l_cmp_i) not in EnumTypes : l_cmp_i=[l_cmp_i]
 268         if type(l_cmp_j) not in EnumTypes : l_cmp_j=[l_cmp_j]
 269         if type(l_ind_i) not in EnumTypes : l_ind_i=[l_ind_i]
 270         if type(l_ind_j) not in EnumTypes : l_ind_j=[l_ind_j]
 271         if len(l_ind_i)!=len(l_ind_j) :
 272            UTMESS('F','PROBA0_8')
 273         if len(l_cmp_i)!=len(l_cmp_j) :
 274            UTMESS('F','PROBA0_9')
 275         if len(l_ind_i)!=len(l_cmp_i) :
 276            UTMESS('F','PROBA0_10')
 277         # paramètres fixes de la table
 278         tabres.add_para(['NOEUD_I','NOEUD_J','NOM_CMP_I','NOM_CMP_J'], 'K8')
 279
 280 #     ------------------------------------------------------------------
 281 #     Cas de tous les indices centraux
 282
 283       elif OPTION!=None :
 284          if NUME_VITE_FLUI_present :
 285                intespec = intespec.NUME_VITE_FLUI == jvite[0]
 286
 287          if NUME_ORDRE_I_present :
 288             imode = list(Set(intespec.NUME_ORDRE_I.not_none_values()))
 289             l_ind_i=imode
 290             l_ind_j=imode
 291             # paramètres fixes de la table
 292             tabres.add_para(['NUME_ORDRE_I','NUME_ORDRE_J'], 'I')
 293          else :
 294             l_ind_i = intespec.NOEUD_I.values()
 295             l_ind_j = intespec.NOEUD_J.values()
 296             if len(l_ind_i) != len(l_ind_j) :
 297                UTMESS('F','PROBA0_8')
 298             l_cmp_i = intespec.NOM_CMP_I.values()
 299             l_cmp_j = intespec.NOM_CMP_J.values()
 300             if (len(l_ind_i) != len(l_cmp_i) or len(l_ind_j) != len(l_cmp_j)) :
 301                UTMESS('F','PROBA0_10')
 302             l_l=zip(zip(l_ind_i,l_cmp_i),zip(l_ind_j,l_cmp_j))
 303             l_ind_i=[]
 304             l_ind_j=[]
 305             l_cmp_i=[]
 306             l_cmp_j=[]
 307             for ai,aj in l_l :
 308                 if ai==aj :
 309                    l_ind_i.append(ai[0])
 310                    l_ind_j.append(aj[0])
 311                    l_cmp_i.append(ai[1])
 312                    l_cmp_j.append(aj[1])
 313             # paramètres fixes de la table
 314             tabres.add_para(['NOEUD_I','NOEUD_J','NOM_CMP_I','NOM_CMP_J'], 'K8')
 315
 316       if jvite[0]!=None :
 317          tabres.add_para('NUME_VITE_FLUI', 'I')
 318
 319
 320 #     ------------------------------------------------------------------
 321 #     Liste des moments spectraux
 322
 323       l_moments=[0,1,2,3,4]
 324       if MOMENT!=None :
 325          l_moments.extend(list(MOMENT))
 326          l_moments=list(Set(l_moments))
 327
 328 #     ------------------------------------------------------------------
 329 #     Boucle sur les fonctions
 330
 331       if NOEUD_I_present :
 332          l_ind=zip(l_ind_i,l_ind_j, l_cmp_i,l_cmp_j)
 333       else :
 334          l_ind=zip(l_ind_i, l_ind_j )
 335
 336
 337       # pour la présentation de la table finale, on stocke le nbre de paramètres "initiaux"
 338       nbpara0 = len(tabres.para)
 339
 340       for vite in jvite :
 341         if INFO==2 :
 342            texte='POUR LA MATRICE INTERSPECTRALE '+INTE_SPEC.nom+'\n'
 343            aster.affiche('MESSAGE',texte)
 344         for ind in l_ind :
 345            dlign = {}
 346            mcfact=[]
 347            if vite!=None :
 348              dlign['NUME_VITE_FLUI'] = vite
 349              mcfact.append(_F(NOM_PARA='NUME_VITE_FLUI',VALE_I=vite))
 350            if NOEUD_I_present :
 351              i_foncstat = ind[0] == ind[1] and  ind[2] == ind[3]
 352              dlign['NOEUD_I'], dlign['NOEUD_J'], dlign['NOM_CMP_I'], dlign['NOM_CMP_J'] = \
 353                   ind[0], ind[1], ind[2], ind[3]
 354              mcfact.append(_F(NOM_PARA='NOEUD_I',VALE_K=ind[0]))
 355              mcfact.append(_F(NOM_PARA='NOEUD_J',VALE_K=ind[1]))
 356              mcfact.append(_F(NOM_PARA='NOM_CMP_I',VALE_K=ind[2]))
 357              mcfact.append(_F(NOM_PARA='NOM_CMP_J',VALE_K=ind[3]))
 358              if INFO==2 :
 359                 aster.affiche('MESSAGE','INDICES :'+ind[0]+' - '+ind[1])
 360                 aster.affiche('MESSAGE','INDICES :'+ind[2]+' - '+ind[3]+'\n')
 361            else :
 362              i_foncstat = ind[0] == ind[1]
 363              dlign['NUME_ORDRE_I'], dlign['NUME_ORDRE_J'] = ind[0], ind[1]
 364              mcfact.append(_F(NOM_PARA='NUME_ORDRE_I',VALE_I=ind[0]))
 365              mcfact.append(_F(NOM_PARA='NUME_ORDRE_J',VALE_I=ind[1]))
 366              if INFO==2 :
 367                 aster.affiche('MESSAGE','INDICES :'+str(ind[0])+' - '\
 368                                                    +str(ind[1])+'\n')
 369
 370            __fon1=RECU_FONCTION(TABLE        = INTE_SPEC,
 371                                 NOM_PARA_TABL= 'FONCTION_C',
 372                                 FILTRE       = mcfact, )
 373
 374            val  = __fon1.Valeurs()
 375            fvalx= Numeric.array(val[0])
 376            fvaly= Numeric.array(val[1])
 377            frez = fvalx[0]
 378
 379            # -- moments spectraux
 380
 381            val_mom={}
 382            for i_mom in l_moments :
 383                trapz     = Numeric.zeros(len(fvaly),Numeric.Float)
 384                trapz[0]  = 0.
 385                valy      = fvaly*(2*pi*fvalx)**i_mom
 386                trapz[1:] = (valy[1:]+valy[:-1])/2*(fvalx[1:]-fvalx[:-1])
 387                prim_y    = Numeric.cumsum(trapz)
 388                val_mom[i_mom] = prim_y[-1]
 389            for i_mom in l_moments :
 390              chmo='LAMBDA_'+str(i_mom).zfill(2)
 391              dlign[chmo] = val_mom[i_mom]
 392
 393         #--- si auto-spectre:
 394            if i_foncstat:
 395               # test si le spectre est bien à valeurs positives
 396               if min(fvaly) < 0.0 :
 397                  aster.affiche('MESSAGE', str(ind)+'\n')
 398                  UTMESS('F','MODELISA9_95')
 399               # -- fonctions statistiques
 400               if NUME_VITE_FLUI or frez >= 0. :
 401                   # -- cas NUME_VITE_FLUI, seule la partie positive du spectre est utilisée
 402                   # -- Il faut donc doubler lambda  pour calculer le bon écart type
 403                   dlign['ECART'] = sqrt(val_mom[0]*2.)
 404               else :
 405                   dlign['ECART'] = sqrt(val_mom[0])
 406               if abs(val_mom[2])>=1e-20 :
 407                     dlign['NB_EXTREMA_P_S'] = 1./pi*sqrt(val_mom[4]/val_mom[2])
 408               if abs(val_mom[0])>=1e-20 :
 409                  dlign['NB_PASS_ZERO_P_S'] = 1./pi*sqrt(val_mom[2]/val_mom[0])
 410                  dlign['FREQ_APPAR'] = 0.5*dlign['NB_PASS_ZERO_P_S']
 411                  if abs(val_mom[4])>=1e-20 :
 412                     dlign['FACT_IRRE'] = sqrt( val_mom[2]*val_mom[2]/val_mom[0]/val_mom[4])
 413
 414            # ajoute la ligne à la Table
 415            tabres.append(dlign)
 416
 417 #--- construction de la table produite
 418
 419       # tri des paramètres
 420       ord_para = tabres.para[nbpara0:]
 421       ord_para.sort()
 422       ord_para = tabres.para[:nbpara0] + ord_para
 423       dprod = tabres[ord_para].dict_CREA_TABLE()
 424
 425       tabout = CREA_TABLE(**dprod)
 426
 427    return ier
 428