src/tests/daSalome/test_aster_zzzz159a_aster_functions.py

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