Aster/Cata/cataSTA9/Macro/post_dyna_alea_ops.py

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