src/tests/daSalome/test_aster_zzzz159a_aster_functions.py

   1 #!/usr/bin/python
   2 # -*- coding: iso-8859-1 -*-
   3
   4 import os, sys, shutil, tempfile, glob
   5 from math import log10
   6
   7 # Variables globales
   8 debug=False
   9 ASTER_ROOT   = ''
  10 SOURCES_ROOT = ''
  11 export = None
  12 calcul = None
  13 parametres = None
  14 python_version = ''
  15
  16 #===============================================================================
  17 def UTMESS(code='I', txt=''):
  18    print txt
  19    if code=='F': sys.exit()
  20
  21 #===============================================================================
  22 def get_tables(tables_calc,tmp_repe_table,prof):
  23    """ Recupere les resultats Aster (Table Aster -> Numeric Python)
  24    """
  25    import Numeric
  26    assert (tables_calc is not None)
  27    assert (tmp_repe_table is not None)
  28
  29    # Import du module lire_table
  30    if os.environ.has_key('ASTER_ROOT'):
  31       version = prof['version'][0]
  32       bibpyt = os.path.join(os.environ['ASTER_ROOT'], version, 'bibpyt')
  33       sys.path.append(bibpyt)
  34       for mdl in glob.glob(os.path.join(bibpyt, '*')):
  35          sys.path.append(os.path.join(os.environ['ASTER_ROOT'], version, 'bibpyt', mdl))
  36    try:
  37       from lire_table_ops import lecture_table
  38    except:
  39       UTMESS('F', "Impossible d'importer le module lire_table!")
  40
  41    reponses = tables_calc
  42    Lrep=[]
  43    _TB = [None]*len(reponses)
  44    for i in range(len(reponses)):
  45       _fic_table = tmp_repe_table + os.sep + "fort."+str(int(100+i))
  46
  47       try:
  48          file=open(_fic_table,'r')
  49          texte=file.read()
  50          file.close()
  51       except Exception, err:
  52          ier=1
  53          message = "Erreur 1!\n" + str(err)
  54          UTMESS('F', message)
  55
  56       try:
  57          table_lue = lecture_table(texte, 1, ' ')
  58          list_para = table_lue.para
  59          tab_lue   = table_lue.values()
  60       except Exception, err:
  61          ier=1
  62          message = "Erreur 2!\n" + str(err)
  63       else:
  64          ier=0
  65
  66       if ier!=0 : UTMESS('F', message)
  67
  68       try:
  69           nb_val = len(tab_lue[ list_para[0] ])
  70           F = Numeric.zeros((nb_val,2), Numeric.Float)
  71           for k in range(nb_val):
  72             F[k][0] = tab_lue[ str(reponses[i][1]) ][k]
  73             F[k][1] = tab_lue[ str(reponses[i][2]) ][k]
  74           Lrep.append(F)
  75       except Exception, err:
  76           message = "Erreur 3!\n" + str(err)
  77           UTMESS('F', message)
  78    resu_calc = Lrep
  79    from N_Parameters import debug
  80    if debug: print 'resu_calc:', resu_calc
  81
  82    return resu_calc
  83 #===============================================================================
  84
  85
  86 #===============================================================================
  87 def Calcul_Aster_Ponctuel( X0 = None ):
  88     #
  89     import numpy
  90     if type(X0) is type(numpy.matrix([])):
  91         X0 = X0.A1.tolist()
  92     else:
  93         X0 = list(X0)
  94     # ----------------------------------------------------------------------------
  95     # Parametres
  96     isFromYacs = globals().get('ASTER_ROOT', None)  # execution via YACS ou en externe
  97     if not isFromYacs:
  98         from N_Parameters import ASTER_ROOT, debug, SOURCES_ROOT, DISPLAY
  99         from N_Study_Parameters import export
 100         from N_MR_Parameters import calcul, parametres
 101     os.environ['ASTER_ROOT'] = ASTER_ROOT
 102
 103     # ----------------------------------------------------------------------------
 104     # Repertoire contenant les resultats des calculs Aster (None = un rep temp est cree)
 105     resudir = globals().get('resudir', None)
 106
 107     # ----------------------------------------------------------------------------
 108     # Parametres remis en forme
 109     list_params = [x[0] for x in parametres]
 110     list_calc = calcul
 111
 112     # ----------------------------------------------------------------------------
 113     # Procedure de calculs distribues
 114     #
 115     # Import des modules python d'ASTK
 116     astk_serv_root = os.path.join(ASTER_ROOT, 'ASTK', 'ASTK_SERV')
 117     sys.path.append(os.path.join(astk_serv_root, 'lib'))
 118     sys.path.append(os.path.join(ASTER_ROOT, 'lib', python_version, 'site-packages'))
 119     if debug:
 120         print sys.path
 121     try:
 122       from asrun.run          import AsRunFactory
 123       from asrun.profil       import ASTER_PROFIL
 124       from asrun.common_func  import get_hostrc
 125       from asrun.utils        import get_timeout
 126       from asrun.parametric   import is_list_of_dict
 127       from asrun.thread       import Dispatcher
 128       from asrun.distrib      import DistribParametricTask
 129     except Exception, e:
 130       print e
 131       UTMESS('F', "Impossible de determiner l'emplacement d'Aster ! Fixer le chemin avec la variable d'environnement ASTER_ROOT.")
 132
 133     # Import des modules supplementaires
 134     sys.path.insert(0, SOURCES_ROOT)
 135     sys.path.insert(0, os.path.join(SOURCES_ROOT, 'sources'))
 136
 137
 138     # result directories
 139     if resudir:
 140        if not os.path.isdir(resudir):
 141           try:    os.mkdir(resudir)
 142           except:
 143              UTMESS('A', "Impossible de creer le repertoire : %s. On utilise un repertoire temporaire" % resudir)
 144              resudir = None
 145     if not resudir: resudir = tempfile.mkdtemp(prefix='tmp_macr_recal_')
 146     flashdir = os.path.join(resudir,'flash')
 147     UTMESS('I', "\n       ASTER Exécution simple\n       Répertoire temporaire de résultats : %s" % resudir)
 148
 149     sys.argv = ['']
 150
 151     run = AsRunFactory()
 152
 153     prof = ASTER_PROFIL(filename=export)
 154     #prof = init_profil_from(run, prof, keep_surch=True)
 155     prof.Set('R', {
 156        'type' : 'repe', 'isrep' : True, 'ul' : 0, 'compr' : False,
 157        'path' : '/tmp/test_param' })
 158
 159     if debug: print prof
 160     prof.WriteExportTo( os.path.join(resudir, 'master.export') )
 161
 162     # get hostrc object
 163     hostrc = get_hostrc(run, prof)
 164
 165     # timeout before rejected a job
 166     timeout = get_timeout(prof)
 167
 168
 169     # list of parameters
 170     list_val = []
 171
 172     # Dictionnaire des parametres du point courant
 173     dic = dict( zip( list_params, X0 ) )
 174     list_val.append( dic )
 175
 176     assert is_list_of_dict(list_val)
 177     nbval = len(list_val)
 178
 179
 180     # Ajout des impressions de tables a la fin du .comm
 181     t = []
 182     reponses = list_calc
 183     for i in range(len(reponses)):
 184         _ul = str(int(100+i))
 185         num_ul = '99'
 186
 187         try:    os.remove( tmp_macr_recal+os.sep+"REPE_TABLE"+os.sep+"fort."+_ul )
 188         except: pass
 189
 190         t.append("\n# Recuperation de la table : " + str(reponses[i][0]) + "\n")
 191         t.append("DEFI_FICHIER(UNITE=" + num_ul + ", FICHIER='" + os.path.join('.', 'REPE_OUT', 'fort.'+_ul) + "',);\n" )
 192         t.append("IMPR_TABLE(TABLE="+str(reponses[i][0])+", FORMAT='ASTER', UNITE="+num_ul+", INFO=1, FORMAT_R='E30.20',);\n")
 193         t.append("DEFI_FICHIER(ACTION='LIBERER', UNITE="+num_ul+",);\n")
 194
 195
 196     # number of threads to follow execution
 197     numthread = 1
 198
 199     # ----- Execute calcutions in parallel using a Dispatcher object
 200     # elementary task...
 201     task = DistribParametricTask(run=run, prof=prof, # IN
 202                        hostrc=hostrc,
 203                        nbmaxitem=0, timeout=timeout,
 204                        resudir=resudir, flashdir=flashdir,
 205                        keywords={'POST_CALCUL': '\n'.join(t)},
 206                        info=1,
 207                        nbnook=0, exec_result=[])            # OUT
 208     # ... and dispatch task on 'list_tests'
 209     etiq = 'calc_%%0%dd' % (int(log10(nbval)) + 1)
 210     labels = [etiq % (i+1) for i in range(nbval)]
 211     couples = zip(labels, list_val)
 212     execution = Dispatcher(couples, task, numthread=numthread)
 213
 214     iret = 0
 215     if task.nbnook > 0:
 216         iret = 4
 217     #run.Sortie(iret)
 218
 219     # Recuperation des tables calculees
 220     seq_FX   = []
 221     seq_FY   = []
 222     seq_DIMS = []
 223     lst_DIAG = []
 224     lst_iter = []
 225     i=0
 226     for c in labels:
 227         tbl = get_tables(tables_calc=list_calc, tmp_repe_table=os.path.join(resudir, c, 'REPE_OUT'), prof=prof)
 228         FX = []
 229         FY = []
 230         ldims = []
 231         for array in tbl:
 232 #           print 'AA1:', array
 233 #           print array[0]
 234             FX.extend([ x[0] for x in array ])
 235             FY.extend([ x[1] for x in array ])
 236             ldims.append(len(array))
 237
 238         # Agregation des resultats
 239         seq_FX.append(FX)
 240         seq_FY.append(FY)
 241         seq_DIMS.append(ldims)
 242         lst_iter.append(i)
 243         i+=1
 244
 245     # Liste des diagnostics
 246     d_diag = {}
 247     for result in task.exec_result:
 248         label = result[0]
 249         diag  = result[2]
 250         d_diag[label] = diag
 251     lst_DIAG = [ d_diag[label] for label in labels]
 252
 253     if debug:
 254         print
 255         print "list_calc =",list_calc
 256         print "seq_FX    =",seq_FX
 257         print "seq_FY    =",seq_FY
 258         print "seq_DIMS  =",seq_DIMS
 259         print "lst_DIAG  =",lst_DIAG
 260         print "lst_iter  =",lst_iter
 261         print
 262
 263     # ----------------------------------------------------------------------------
 264     # Procedure d'assemblage du gradient
 265
 266     # Calcul maitre (point X0)
 267     idx0 = lst_iter.index(0)   # index (les calculs arrivent-ils dans le desordre?)
 268     FY_X0 = seq_FY[idx0]
 269     H_de_X = FY_X0
 270
 271     # Arret si tous les jobs ne se sont pas deroules correctement
 272     for diag in lst_DIAG:
 273         if not diag[0:2] in ['OK', '<A']:
 274             raise ValueError("Au moins un calcul ne s'est pas deroule normalement")
 275
 276     if debug:
 277         print "\nH_de_X (calcul ponstuel) au point X0: \n%s" % str(H_de_X)
 278
 279     return H_de_X
 280
 281 #===============================================================================
 282 def Calcul_Aster_Jacobienne( X0 = None ):
 283     #
 284     import numpy
 285     if type(X0) is type(numpy.matrix([])):
 286         X0 = X0.A1.tolist()
 287     else:
 288         X0 = list(X0)
 289     # ----------------------------------------------------------------------------
 290     FacteurdX = 1.e-4
 291     dX = globals().get('dX', [ x*FacteurdX for x in X0 ])  # execution via YACS ou en externe
 292     # dX = globals().get('dX', [ 0.1, 0.1, 0.001])  # execution via YACS ou en externe
 293     # ----------------------------------------------------------------------------
 294     # Parametres
 295     isFromYacs = globals().get('ASTER_ROOT', None)  # execution via YACS ou en externe
 296     if not isFromYacs:
 297         from N_Parameters import ASTER_ROOT, debug, SOURCES_ROOT, DISPLAY
 298         from N_Study_Parameters import export
 299         from N_MR_Parameters import calcul, parametres
 300     os.environ['ASTER_ROOT'] = ASTER_ROOT
 301
 302     # ----------------------------------------------------------------------------
 303     # Repertoire contenant les resultats des calculs Aster (None = un rep temp est cree)
 304     resudir = globals().get('resudir', None)
 305
 306     # ----------------------------------------------------------------------------
 307     # Parametres remis en forme
 308     list_params = [x[0] for x in parametres]
 309     list_calc = calcul
 310
 311     # ----------------------------------------------------------------------------
 312     # Procedure de calculs distribues
 313     #
 314     # Import des modules python d'ASTK
 315     astk_serv_root = os.path.join(ASTER_ROOT, 'ASTK', 'ASTK_SERV')
 316     sys.path.append(os.path.join(astk_serv_root, 'lib'))
 317     sys.path.append(os.path.join(ASTER_ROOT, 'lib', python_version, 'site-packages'))
 318     if debug:
 319         print sys.path
 320     try:
 321       from asrun.run          import AsRunFactory
 322       from asrun.profil       import ASTER_PROFIL
 323       from asrun.common_func  import get_hostrc
 324       from asrun.utils        import get_timeout
 325       from asrun.parametric   import is_list_of_dict
 326       from asrun.thread       import Dispatcher
 327       from asrun.distrib      import DistribParametricTask
 328     except Exception, e:
 329       print e
 330       UTMESS('F', "Impossible de determiner l'emplacement d'Aster ! Fixer le chemin avec la variable d'environnement ASTER_ROOT.")
 331
 332     # Import des modules supplementaires
 333     sys.path.insert(0, SOURCES_ROOT)
 334     sys.path.insert(0, os.path.join(SOURCES_ROOT, 'sources'))
 335
 336
 337     # result directories
 338     if resudir:
 339        if not os.path.isdir(resudir):
 340           try:    os.mkdir(resudir)
 341           except:
 342              UTMESS('A', "Impossible de creer le repertoire : %s. On utilise un repertoire temporaire" % resudir)
 343              resudir = None
 344     if not resudir: resudir = tempfile.mkdtemp(prefix='tmp_macr_recal_')
 345     flashdir = os.path.join(resudir,'flash')
 346     UTMESS('I', "\n       ASTER Exécutions multiples\n       Répertoire temporaire de résultats : %s" % resudir)
 347
 348     sys.argv = ['']
 349
 350     run = AsRunFactory()
 351
 352     prof = ASTER_PROFIL(filename=export)
 353     #prof = init_profil_from(run, prof, keep_surch=True)
 354     prof.Set('R', {
 355        'type' : 'repe', 'isrep' : True, 'ul' : 0, 'compr' : False,
 356        'path' : '/tmp/test_param' })
 357
 358     if debug: print prof
 359     prof.WriteExportTo( os.path.join(resudir, 'master.export') )
 360
 361     # get hostrc object
 362     hostrc = get_hostrc(run, prof)
 363
 364     # timeout before rejected a job
 365     timeout = get_timeout(prof)
 366
 367
 368     # list of parameters
 369     list_val = []
 370
 371     # Dictionnaire des parametres du point courant
 372     dic = dict( zip( list_params, X0 ) )
 373     list_val.append( dic )
 374
 375     # Dictionnaires des parametres des calculs esclaves (perturbations des differences finies)
 376     for n in range(1,len(dX)+1):
 377         l = [0] * len(dX)
 378         l[n-1] = dX[n-1]
 379         X = [ X0[i] + l[i] for i in range(len(dX)) ]
 380         dic = dict( zip( list_params, X ) )
 381         list_val.append( dic )
 382
 383     assert is_list_of_dict(list_val)
 384     nbval = len(list_val)
 385
 386
 387     # Ajout des impressions de tables a la fin du .comm
 388     t = []
 389     reponses = list_calc
 390     for i in range(len(reponses)):
 391         _ul = str(int(100+i))
 392         num_ul = '99'
 393
 394         try:    os.remove( tmp_macr_recal+os.sep+"REPE_TABLE"+os.sep+"fort."+_ul )
 395         except: pass
 396
 397         t.append("\n# Recuperation de la table : " + str(reponses[i][0]) + "\n")
 398         t.append("DEFI_FICHIER(UNITE=" + num_ul + ", FICHIER='" + os.path.join('.', 'REPE_OUT', 'fort.'+_ul) + "',);\n" )
 399         t.append("IMPR_TABLE(TABLE="+str(reponses[i][0])+", FORMAT='ASTER', UNITE="+num_ul+", INFO=1, FORMAT_R='E30.20',);\n")
 400         t.append("DEFI_FICHIER(ACTION='LIBERER', UNITE="+num_ul+",);\n")
 401
 402
 403     # number of threads to follow execution
 404     numthread = 1
 405
 406     # ----- Execute calcutions in parallel using a Dispatcher object
 407     # elementary task...
 408     task = DistribParametricTask(run=run, prof=prof, # IN
 409                        hostrc=hostrc,
 410                        nbmaxitem=0, timeout=timeout,
 411                        resudir=resudir, flashdir=flashdir,
 412                        keywords={'POST_CALCUL': '\n'.join(t)},
 413                        info=1,
 414                        nbnook=0, exec_result=[])            # OUT
 415     # ... and dispatch task on 'list_tests'
 416     etiq = 'calc_%%0%dd' % (int(log10(nbval)) + 1)
 417     labels = [etiq % (i+1) for i in range(nbval)]
 418     couples = zip(labels, list_val)
 419     execution = Dispatcher(couples, task, numthread=numthread)
 420
 421     iret = 0
 422     if task.nbnook > 0:
 423         iret = 4
 424     #run.Sortie(iret)
 425
 426     # Recuperation des tables calculees
 427     seq_FX   = []
 428     seq_FY   = []
 429     seq_DIMS = []
 430     lst_DIAG = []
 431     lst_iter = []
 432     i=0
 433     for c in labels:
 434         tbl = get_tables(tables_calc=list_calc, tmp_repe_table=os.path.join(resudir, c, 'REPE_OUT'), prof=prof)
 435         FX = []
 436         FY = []
 437         ldims = []
 438         for array in tbl:
 439 #           print 'AA1:', array
 440 #           print array[0]
 441             FX.extend([ x[0] for x in array ])
 442             FY.extend([ x[1] for x in array ])
 443             ldims.append(len(array))
 444
 445         # Agregation des resultats
 446         seq_FX.append(FX)
 447         seq_FY.append(FY)
 448         seq_DIMS.append(ldims)
 449         lst_iter.append(i)
 450         i+=1
 451
 452     # Liste des diagnostics
 453     d_diag = {}
 454     for result in task.exec_result:
 455         label = result[0]
 456         diag  = result[2]
 457         d_diag[label] = diag
 458     lst_DIAG = [ d_diag[label] for label in labels]
 459
 460     if debug:
 461         print
 462         print "list_calc =",list_calc
 463         print "seq_FX    =",seq_FX
 464         print "seq_FY    =",seq_FY
 465         print "seq_DIMS  =",seq_DIMS
 466         print "lst_DIAG  =",lst_DIAG
 467         print "lst_iter  =",lst_iter
 468         print "dX        =",dX
 469         print
 470
 471     # ----------------------------------------------------------------------------
 472     # Procedure d'assemblage du gradient
 473
 474     # Calcul maitre (point X0)
 475     idx0 = lst_iter.index(0)   # index (les calculs arrivent-ils dans le desordre?)
 476     FY_X0 = seq_FY[idx0]
 477     H_de_X = FY_X0
 478
 479     # Arret si tous les jobs ne se sont pas deroules correctement
 480     for diag in lst_DIAG:
 481         if not diag[0:2] in ['OK', '<A']:
 482             raise ValueError("Au moins un calcul ne s'est pas deroule normalement")
 483
 484     # Calcul du gradient (une liste de liste)
 485     Gradient_de_H_en_X = []
 486
 487     for n in range(len(lst_iter))[1:]:
 488         idx = lst_iter.index(n)
 489         FY   = seq_FY[idx]
 490         col = [ -(y-x)/dX[idx-1] for x, y in zip(FY, FY_X0) ]
 491         Gradient_de_H_en_X.append(col)
 492         if debug: print 'Calcul numero: %s - Diagnostic: %s' % (n, lst_DIAG[idx])
 493
 494     if debug:
 495         print "\nCalcul H au point X0: \n%s" % str(H_de_X)
 496         import pprint
 497         print "\nGradient au point X0:"
 498         pprint.pprint(Gradient_de_H_en_X)
 499
 500     return Gradient_de_H_en_X
 501
 502 #===============================================================================
 503 def Calcul_Aster_Adjoint( (X0, dY) ):
 504     #
 505     if 0:
 506         print
 507         print "CALCUL ADJOINT"
 508         print "X0 =",X0
 509         print "dY =",dY
 510     #
 511     import numpy
 512     #
 513     Y0 = numpy.asmatrix(dY).flatten().T
 514     #
 515     Delta_HX = Calcul_Aster_Jacobienne( X0 )
 516     Delta_HX = numpy.matrix( Delta_HX )
 517     #
 518     HtY = numpy.dot(Delta_HX, Y0)
 519     #
 520     if 0:
 521         print "dHX =",Delta_HX
 522         print "HtY =",HtY
 523         print
 524     #
 525     return HtY.A1