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