Aster/Cata/cataSTA10/Macro/post_dyna_alea_ops.py

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