Aster/Cata/cataSTA10/Macro/recal.py

   1 #@ MODIF recal Macro  DATE 26/05/2010   AUTEUR ASSIRE A.ASSIRE
   2 # -*- coding: iso-8859-1 -*-
   3 #            CONFIGURATION MANAGEMENT OF EDF VERSION
   4 # ======================================================================
   5 # COPYRIGHT (C) 1991 - 2002  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
  22 #___________________________________________________________________________
  23 #
  24 #           MODULE DE CALCUL DISTRIBUE POUR MACR_RECAL
  25 #
  26 #  Utilisable en mode EXTERNE, voir les flags avec "python recal.py -h"
  27 #___________________________________________________________________________
  28
  29
  30 import os
  31 import sys
  32 import shutil
  33 import tempfile
  34 import glob
  35 import math
  36 import copy
  37 import re
  38 import platform
  39 from math import log10
  40
  41 import numpy as NP
  42
  43
  44 # Importation de commandes Aster
  45 try:
  46    import aster
  47    import Macro
  48    from Accas import _F
  49 except ImportError:
  50    pass
  51
  52 include_pattern = "# -->INCLUDE<--"
  53 debug = False
  54
  55 # -------------------------------------------------------------------------------
  56 # -------------------------------------------------------------------------------
  57 def get_absolute_path(path):
  58    """Retourne le chemin absolu en suivant les liens éventuels.
  59    """
  60    if os.path.islink(path):
  61       path = os.path.realpath(path)
  62    res = os.path.normpath(os.path.abspath(path))
  63    return res
  64
  65 # -------------------------------------------------------------------------------
  66 #if os.environ.has_key('bibpytdir'): sys.path.append( os.environ['bibpytdir'] )
  67
  68 # recupere "bibpyt" à partir de "bibpyt/Macro/recal.py"
  69 sys.path.append(get_absolute_path(os.path.join(sys.argv[0], '..', '..')))
  70
  71 try:
  72    from Utilitai.Utmess import UTMESS
  73 except:
  74    def UTMESS(code='I', txt='',valk='', vali='', valr=''):
  75        print txt, valk, vali, valr
  76        if code=='F': sys.exit()
  77
  78
  79 # # -------------------------------------------------------------------------------
  80 # def find_parameter(content, param):
  81 #    """
  82 #    Return the lowest index in content where param is found and
  83 #    the index of the end of the command.
  84 #    """
  85 #    pos, endpos = -1, -1
  86 #    re_start = re.compile('^ *%s *\=' % re.escape(param), re.M)
  87 #    mat_start = re_start.search(content)
  88 #    if mat_start is not None:
  89 #       pos = mat_start.start()
  90 #       endpos = search_enclosed(content, pos)
  91 #    return pos, endpos
  92
  93
  94
  95 # -------------------------------------------------------------------------------
  96 def find_parameter(content, param):
  97     """
  98     Supprime les parametres du fichier de commande
  99     """
 100     re_start = re.compile('^ *%s *\=' % re.escape(param), re.M)
 101     l=[]
 102     for line in content.split('\n'):
 103        mat_start = re_start.search(line)
 104        if mat_start is None: l.append(line)
 105     return '\n'.join(l)
 106
 107
 108 # -------------------------------------------------------------------------------
 109 def Affiche_Param(para, val):
 110     """Affiche les parametres
 111     """
 112     t = []
 113     for p, v in zip(para, val):
 114         t.append( "     %s : %s" % ( p.ljust(9), v) )
 115     return '\n'.join(t)
 116
 117
 118 # -------------------------------------------------------------------------------
 119 def make_include_files(UNITE_INCLUDE, calcul, parametres):
 120    """  Module permettant de generer les fichiers a inclure (mode INCLUSION)
 121    """
 122
 123 #    # Importation de commandes Aster
 124 #    try:
 125 #       import aster
 126 #       import Macro
 127 #       from Accas import _F
 128 #       from Cata.cata import *
 129 #    except ImportError:
 130 #       raise "Le mode INCLUSION doit etre lance depuis Aster"
 131
 132    try:
 133        ASTER_ROOT = os.path.join(aster.repout, '..')
 134        sys.path.append(os.path.join(ASTER_ROOT, 'ASTK', 'ASTK_SERV', 'lib'))
 135        sys.path.append(os.path.join(ASTER_ROOT, 'lib', 'python%s.%s' % (sys.version_info[0], sys.version_info[1] ) , 'site-packages'))
 136    except: pass
 137    try:
 138        from asrun.utils import find_command, search_enclosed
 139    except Exception, e:
 140        print e
 141        UTMESS('F','RECAL0_99')
 142
 143
 144    # ----------------------------------------------------------------------------
 145    # Preparation des fichiers
 146    # ----------------------------------------------------------------------------
 147    liste_reponses = []
 148    for reponse in [ x[0] for x in calcul ]:
 149       if not reponse in liste_reponses: liste_reponses.append(reponse)
 150
 151    try:
 152        old = "fort.%s"     % UNITE_INCLUDE
 153        pre = "fort.%s.pre" % UNITE_INCLUDE
 154        new = "fort.%s.new" % UNITE_INCLUDE
 155
 156        # Lecture du fichier
 157        f=open(old, 'r')
 158        newtxt = f.read()
 159        f.close()
 160
 161        # On retire la commande DEBUT
 162        pos, endpos = find_command(newtxt, "DEBUT")
 163        if endpos!=-1: newtxt = newtxt[endpos+1:]
 164        if newtxt[0]==';': newtxt = newtxt[1:]  # Bug dans find_command si la commande se termine par un ";"
 165
 166        # On retire les parametres
 167        list_params = [x[0] for x in parametres]
 168        for param in list_params:
 169            newtxt = find_parameter(newtxt, param)
 170
 171        # Isole la partie a inclure si elle est specifiee
 172        n = newtxt.find(include_pattern)
 173        pretxt = None
 174        if n!=-1:
 175            pretxt = newtxt[:n]
 176            pretxt = "# -*- coding: iso-8859-1 -*-\n" + pretxt
 177            # Ecriture du nouveau fichier
 178            fw=open(pre, 'w')
 179            fw.write(pretxt)
 180            fw.close()
 181            newtxt = newtxt[n+len(include_pattern):]
 182
 183        # Retire la commande FIN
 184        pos, endpos = find_command(newtxt, "FIN")
 185        if pos!=-1: newtxt = newtxt[:pos]
 186
 187        # Ajoute un global pour ramener les courbes dans l'espace Aster
 188        newtxt = "global %s\n" % ','.join(liste_reponses) + newtxt
 189
 190        # Ajoute un encodage pour eviter les erreurs dues aux accents (ssna110a par exemple)
 191        newtxt = "# -*- coding: iso-8859-1 -*-\n" + newtxt
 192
 193        # Ecriture du nouveau fichier
 194        fw=open(new, 'w')
 195        fw.write(newtxt)
 196        fw.close()
 197    except Exception, e:
 198        raise e
 199
 200    return
 201
 202
 203 # -------------------------------------------------------------------------------
 204 def mes_concepts(list_concepts=[],base=None):
 205    """ Fonction qui liste les concepts créés """
 206    for e in base.etapes:
 207       if e.nom in ('INCLUDE','MACR_RECAL',) :
 208          list_concepts=list(mes_concepts(list_concepts=list_concepts,base=e))
 209       elif (e.sd != None) and (e.parent.nom=='INCLUDE') :
 210          nom_concept=e.sd.get_name()
 211          if not(nom_concept in list_concepts):
 212             list_concepts.append( nom_concept )
 213    return tuple(list_concepts)
 214
 215
 216 # -------------------------------------------------------------------------------
 217 def detr_concepts(self):
 218      liste_concepts=mes_concepts(base=self.parent)
 219      for e in liste_concepts:
 220         nom = string.strip(e)
 221         DETRUIRE( OBJET =self.g_context['_F'](CHAINE = nom), INFO=2)
 222         if self.jdc.g_context.has_key(nom) : del self.jdc.g_context[nom]
 223      del(liste_concepts)
 224
 225
 226 # -------------------------------------------------------------------------------
 227 def get_tables(tables_calc, tmp_repe_table, prof):
 228    """ Recupere les resultats Aster (Table Aster -> numpy)
 229    """
 230    assert (tables_calc is not None)
 231    assert (tmp_repe_table is not None)
 232
 233    # Import du module lire_table
 234    if os.environ.has_key('ASTER_ROOT'):
 235       version = prof['version'][0]
 236       bibpyt = os.path.join(os.environ['ASTER_ROOT'], version, 'bibpyt')
 237       sys.path.append(bibpyt)
 238       for mdl in glob.glob(os.path.join(bibpyt, '*')):
 239          sys.path.append(os.path.join(os.environ['ASTER_ROOT'], version, 'bibpyt', mdl))
 240    try:
 241       from lire_table_ops import lecture_table
 242    except:
 243       UTMESS('F','RECAL0_23')
 244
 245    reponses = tables_calc
 246    Lrep=[]
 247    for i in range(len(reponses)):
 248       _fic_table = tmp_repe_table + os.sep + "fort."+str(int(100+i))
 249
 250       try:
 251          f=open(_fic_table,'r')
 252          texte=f.read()
 253          f.close()
 254       except Exception, err:
 255          ier=1
 256          UTMESS('F','RECAL0_24',valk=str(err))
 257
 258       try:
 259          table_lue = lecture_table(texte, 1, ' ')
 260          list_para = table_lue.para
 261          tab_lue   = table_lue.values()
 262       except Exception, err:
 263          ier=1
 264       else:
 265          ier=0
 266
 267       if ier!=0 : UTMESS('F','RECAL0_24',valk=str(err))
 268
 269       F = table2numpy(tab_lue, list_para, reponses, i)
 270       Lrep.append(F)
 271
 272
 273    return Lrep
 274
 275
 276 # --------------------------------------------------------------------------------------------------
 277 def table2numpy(tab_lue, list_para, reponses, i):
 278    """  Extraction des resultats depuis la table Aster
 279    """
 280    try:
 281        nb_val = len(tab_lue[ list_para[0] ])
 282        F = NP.zeros((nb_val,2))
 283        for k in range(nb_val):
 284          F[k][0] = tab_lue[ str(reponses[i][1]) ][k]
 285          F[k][1] = tab_lue[ str(reponses[i][2]) ][k]
 286    except Exception, err:
 287        UTMESS('F','RECAL0_24',valk=str(err))
 288    return F
 289
 290
 291 # --------------------------------------------------------------------------------------------------
 292 def Ecriture_Fonctionnelle(output_file, type_fonctionnelle, fonctionnelle):
 293
 294    try:    os.remove(output_file)
 295    except: pass
 296
 297    f=open(output_file, 'w')
 298    if type_fonctionnelle == 'vector':
 299       try:    fonctionnelle = fonctionnelle.tolist()
 300       except: pass
 301       fonctionnelle = str(fonctionnelle).replace('[','').replace(']','').replace('\n', ' ')
 302    f.write( str(fonctionnelle) )
 303    f.close()
 304
 305
 306 # --------------------------------------------------------------------------------------------------
 307 def Ecriture_Derivees(output_file, derivees):
 308
 309    try:    os.remove(output_file)
 310    except: pass
 311
 312    # on cherche a imprimer la gradient calcule a partir de Fcalc
 313    if type(derivees) in [list, tuple]:
 314        t = []
 315        for l in derivees:
 316           l = str(l).replace('[', '').replace(']', '')
 317           t.append( l )
 318        txt = '\n'.join(t)
 319
 320    # On cherche a imprimer la matrice des sensibilite (A ou A_nodim)
 321    elif type(derivees) == NP.ndarray:
 322        t = []
 323        a = derivees
 324        for c in range(len(a[0,:])):
 325            l = a[:,c].tolist()
 326            l = str(l).replace('[', '').replace(']', '')
 327            t.append( l )
 328        txt = '\n'.join(t)
 329
 330    else: raise "Wrong type for gradient !"
 331
 332    # Ecriture
 333    f=open(output_file, 'w')
 334    f.write(txt)
 335    f.close()
 336
 337
 338
 339 # --------------------------------------------------------------------------------------------------
 340 # --------------------------------------------------------------------------------------------------
 341 class CALCULS_ASTER:
 342    """
 343       Classe gérant les calculs Aster (distribues ou include)
 344    """
 345
 346    # ---------------------------------------------------------------------------
 347    def __init__(self,
 348
 349        # MACR_RECAL inputs are optional here (if passed to self.run methods)
 350        parametres          = None,
 351        calcul              = None,
 352        experience          = None,
 353        LANCEMENT        = 'DISTRIBUTION',
 354        jdc                 = None,
 355                ):
 356
 357        self.parametres         = parametres
 358        self.calcul             = calcul
 359        self.experience         = experience
 360        #self.eval_esclave       = mode_esclave
 361        self.LANCEMENT       = LANCEMENT
 362        self.UNITE_ESCL         = None
 363        self.UNITE_INCLUDE      = None
 364        self.ASTER_ROOT         = None
 365
 366        self.jdc                = jdc
 367
 368        self.list_params        = [x[0] for x in parametres]
 369        self.list_params.sort()
 370
 371        # Valable uniquement pour le mode INCLUDE
 372        self.pre    = None
 373        self.pretxt = None
 374        self.new    = None
 375        self.newtxt = None
 376
 377        # Mode dynamique desactive par defaut
 378        self.SetDynamiqueMode(None, None)
 379
 380
 381    # ---------------------------------------------------------------------------------------------------------
 382    def SetDynamiqueMode(self, DYNAMIQUE, graph_mac):
 383        self.DYNAMIQUE = DYNAMIQUE
 384        self.graph_mac = graph_mac
 385
 386
 387    # ---------------------------------------------------------------------------------------------------------
 388    # ---------------------------------------------------------------------------------------------------------
 389    def run(self,
 390
 391                 # Current estimation
 392                 X0,
 393                 dX             = None,
 394
 395                 # Code_Aster installation
 396                 ASTER_ROOT     = None,
 397                 as_run         = None,
 398
 399                 # General
 400                 resudir        = None,
 401                 clean          = True,
 402                 info           = None,
 403                 NMAX_SIMULT   = None,
 404
 405                 # Study
 406                 export         = None,
 407
 408                 # MACR_RECAL inputs
 409                 parametres     = None,
 410                 calcul         = None,
 411                 experience     = None,
 412         ):
 413
 414         # Current estimation
 415         self.X0             = X0
 416         self.dX             = dX
 417
 418         # Code_Aster installation
 419         self.ASTER_ROOT     = ASTER_ROOT
 420         self.as_run         = as_run
 421
 422         # General
 423         self.resudir        = resudir
 424         self.clean          = clean
 425         self.info           = info
 426         if not NMAX_SIMULT: NMAX_SIMULT = 0
 427         self.NMAX_SIMULT   = NMAX_SIMULT
 428
 429         # Study
 430         self.export         = export
 431
 432         # MACR_RECAL inputs
 433         if parametres:   self.parametres     = parametres
 434         if calcul:       self.calcul         = calcul
 435         if experience:   self.experience     = experience
 436
 437         parametres  = self.parametres
 438         calcul      = self.calcul
 439         experience  = self.experience
 440
 441         list_params = self.list_params
 442
 443         if dX: CalcGradient = True
 444         else:  CalcGradient = False
 445         self.CalcGradient   = CalcGradient
 446
 447         self.list_diag      = []
 448
 449         # Pour le moment on conserve un seul fichier
 450         self.UNITE_INCLUDE  = self.UNITE_ESCL
 451
 452
 453         # ----------------------------------------------------------------------------
 454         # Liste de tous les jeux de parametres (initial + differences finies)
 455         # ----------------------------------------------------------------------------
 456         list_val = []
 457
 458         # Dictionnaire des parametres du point courant
 459         dic = dict( zip( list_params, X0 ) )
 460         list_val.append( dic )
 461
 462         # Calcul du gradient (perturbations des differences finies)
 463         if CalcGradient:
 464             UTMESS('I','RECAL0_16')
 465             # Dictionnaires des parametres des calculs esclaves
 466             for n in range(1,len(dX)+1):
 467                l = [0] * len(dX)
 468                l[n-1] = dX[n-1]
 469                X = [ X0[i] * (1+l[i]) for i in range(len(dX)) ]
 470                dic = dict( zip( list_params, X ) )
 471                list_val.append( dic )
 472
 473
 474         # ----------------------------------------------------------------------------
 475         # Aiguillage vers INCLUDE
 476         # ----------------------------------------------------------------------------
 477         if self.LANCEMENT == 'INCLUSION':
 478            UTMESS('I','RECAL0_29', valk=self.LANCEMENT)
 479            fonctionnelle, gradient = self.run_include(list_val)
 480
 481
 482         # ----------------------------------------------------------------------------
 483         # Aiguillage vers ASRUN distribue
 484         # ----------------------------------------------------------------------------
 485         elif self.LANCEMENT == 'DISTRIBUTION':
 486            UTMESS('I','RECAL0_29', valk=self.LANCEMENT)
 487            fonctionnelle, gradient = self.run_distrib(list_val)
 488
 489
 490         # ----------------------------------------------------------------------------
 491         # Erreur d'aiguillage
 492         # ----------------------------------------------------------------------------
 493         else:
 494            raise "Erreur : mode %s inconnu!" % self.LANCEMENT
 495
 496
 497         #sys.exit()
 498         # ----------------------------------------------------------------------------
 499         # Sortie
 500         # ----------------------------------------------------------------------------
 501         return fonctionnelle, gradient
 502
 503
 504
 505    # ---------------------------------------------------------------------------------------------------------
 506    # ---------------------------------------------------------------------------------------------------------
 507    def run_include(self,list_val):
 508      """  Module permettant de lancer N+1 calculs via un mecanisme d'include
 509      """
 510
 511 #      # Importation de commandes Aster
 512 #      try:
 513 #         import aster
 514 #         import Macro
 515 #         from Accas import _F
 516 #         from Cata import cata
 517 #         from Cata.cata import *
 518 #      except ImportError:
 519 #         raise "Simu_point_mat doit etre lance depuis Aster"
 520
 521      try:
 522          import aster
 523          import Macro
 524          from Cata import cata
 525          from Cata.cata import OPER, MACRO
 526          from Accas import _F
 527
 528          # Declaration de toutes les commandes Aster
 529          import cata
 530          for k,v in cata.__dict__.items() :
 531            #print k,v
 532            if isinstance(v, (OPER, MACRO)):
 533              #print k,v
 534              #self.jdc.current_context[k]= v
 535              exec("from Cata.cata import %s" % k)
 536          #self.jdc.current_context['_F']=cata.__dict__['_F']
 537      except Exception, e:
 538          raise "Le mode INCLUDE doit etre lance depuis Aster : \nErreur : " % e
 539
 540
 541      list_params = self.list_params
 542      calcul      = self.calcul
 543      reponses    = self.calcul
 544
 545 # AA : marche pas car on peut oublier des courbes, tant pis on refait des extract en trop..
 546 #      liste_reponses = []
 547 #      for reponse in [ x[0] for x in calcul ]:
 548 #          if not reponse in liste_reponses: liste_reponses.append(reponse)
 549
 550      liste_reponses = [ x[0] for x in calcul ]
 551
 552
 553      # ----------------------------------------------------------------------------
 554      # Boucle sur les N+1 calculs
 555      # ----------------------------------------------------------------------------
 556      Lcalc = []
 557      for i in range(len(list_val)):
 558          params = list_val[i]
 559
 560
 561          # ----------------------------------------------------------------------------
 562          # Affectation des valeurs des parametres
 563          # ----------------------------------------------------------------------------
 564          for nompara in list_params:
 565              valpara = params[nompara]
 566              exec( "%s=%s" % (nompara, valpara) )    #  YOUN__ = X0[0], DSDE__ = X0[1], ...
 567
 568
 569          # ----------------------------------------------------------------------------
 570          # Affichage des parametres du calcul courant
 571          # ----------------------------------------------------------------------------
 572          tpara = Affiche_Param(list_params, [ params[x] for x in list_params] )
 573          if i==0:  UTMESS('I', 'RECAL0_67', valk=tpara)
 574          else:     UTMESS('I', 'RECAL0_68', valk=(tpara, list_params[i-1]) )
 575
 576
 577          # ----------------------------------------------------------------------------
 578          # Lancement du calcul (par un include)
 579          # ----------------------------------------------------------------------------
 580          new = "fort.%s.new" % self.UNITE_INCLUDE
 581          execfile(new)
 582
 583
 584          # ----------------------------------------------------------------------------
 585          # On considere que le job est OK s'il ne s'est pas plante dans le except precedent..
 586          # ----------------------------------------------------------------------------
 587          self.list_diag.append("OK")
 588
 589
 590          # ----------------------------------------------------------------------------
 591          # Extraction des tables
 592          # ----------------------------------------------------------------------------
 593          Lrep=[]
 594          for i in range(len(liste_reponses)):
 595              reponse = liste_reponses[i]
 596              DETRUIRE(OBJET=_F(CHAINE='VEXTR___'), ALARME='NON', INFO=1)  # Faudrait proteger ce nom ici (VEXTR___ peut etre deja utilise dans l'etude)
 597              exec( "VEXTR___ = %s.EXTR_TABLE()" % reponse)
 598              list_para = VEXTR___.para
 599              tab_lue   = VEXTR___.values()
 600              F = table2numpy(tab_lue, list_para, reponses, i)
 601              Lrep.append(F)
 602
 603
 604
 605          Lcalc.append( Lrep )
 606
 607
 608          # Destruction des concepts Aster
 609          liste_concepts = self.jdc.g_context.keys()
 610          for c in liste_concepts:
 611              DETRUIRE(OBJET=_F(CHAINE=c), ALARME='NON', INFO=1);
 612
 613          #detr_concepts(self.jdc)  # marche pas !
 614          #sys.exit()
 615
 616
 617      # ----------------------------------------------------------------------------
 618      # Calcul de la fonctionnelle et du gradient
 619      # ----------------------------------------------------------------------------
 620      if debug: print "AA4/Lcalc=", Lcalc
 621      fonctionnelle, gradient = self.calc2fonc_gradient(Lcalc)
 622
 623
 624      # ----------------------------------------------------------------------------
 625      # Save all calculated responses
 626      self.Lcalc = Lcalc
 627      # ----------------------------------------------------------------------------
 628
 629
 630      return fonctionnelle, gradient
 631
 632
 633
 634
 635
 636
 637    # ---------------------------------------------------------------------------------------------------------
 638    # ---------------------------------------------------------------------------------------------------------
 639    def run_distrib(self, list_val):
 640         """ Module permettant de lancer N+1 calculs avec le module de calculs distribues d'asrun
 641         """
 642
 643         # ----------------------------------------------------------------------------
 644         # Parametres
 645         # ----------------------------------------------------------------------------
 646
 647         # Code_Aster installation
 648         ASTER_ROOT     = self.ASTER_ROOT
 649         as_run         = self.as_run
 650
 651         # General
 652         resudir        = self.resudir
 653         clean          = self.clean
 654         info           = self.info
 655
 656         # Study
 657         export         = self.export
 658
 659         # MACR_RECAL inputs
 660         parametres     = self.parametres
 661         calcul         = self.calcul
 662         experience     = self.experience
 663
 664         parametres     = self.parametres
 665         calcul         = self.calcul
 666         experience     = self.experience
 667
 668         CalcGradient   = self.CalcGradient
 669         NMAX_SIMULT   = self.NMAX_SIMULT
 670
 671
 672         # ----------------------------------------------------------------------------
 673         # Import des modules python d'ASTK
 674         # ----------------------------------------------------------------------------
 675         if not ASTER_ROOT:
 676             try:    ASTER_ROOT = os.path.join(aster.repout, '..')
 677             except: pass
 678         try:
 679             sys.path.append(os.path.join(ASTER_ROOT, 'ASTK', 'ASTK_SERV', 'lib'))
 680             sys.path.append(os.path.join(ASTER_ROOT, 'lib', 'python%s.%s' % (sys.version_info[0], sys.version_info[1] ) , 'site-packages'))
 681         except: pass
 682         try:
 683             from asrun.run          import AsRunFactory
 684             from asrun.profil       import ASTER_PROFIL
 685             from asrun.common_func  import get_hostrc
 686             from asrun.utils        import get_timeout
 687             from asrun.parametric   import is_list_of_dict
 688             from asrun.thread       import Dispatcher
 689             from asrun.distrib      import DistribParametricTask
 690         except Exception, e:
 691             print e
 692             UTMESS('F','RECAL0_99')
 693
 694
 695         assert is_list_of_dict(list_val)
 696         nbval = len(list_val)
 697
 698         # ----------------------------------------------------------------------------
 699         # Generation des etudes esclaves
 700         # ----------------------------------------------------------------------------
 701         sys.argv = ['']
 702         run = AsRunFactory()
 703         run.options['debug_stderr'] = True  # pas d'output d'executions des esclaves dans k'output maitre
 704
 705         # Master profile
 706         prof = ASTER_PROFIL(filename=export)
 707         tmp_param = tempfile.mkdtemp()
 708         try:    username = prof.param['username'][0]
 709         except: username = os.environ['LOGNAME']
 710         try:    noeud    = prof.param['noeud'][0]
 711         except: noeud    = platform.uname()[1]
 712         tmp_param = "%s@%s:%s" % ( username, noeud, tmp_param)
 713         prof.Set('R', {'type' : 'repe', 'isrep' : True, 'ul' : 0, 'compr' : False, 'path' : tmp_param })
 714         if info>=2: print prof
 715
 716         # Si batch n'est pas possible, on bascule en interactif
 717         if prof.param['mode'][0]=='batch' and run.get('batch')=='non':
 718            UTMESS('I','RECAL0_28',valk=noeud)
 719            prof.param['mode'][0] = 'interactif'
 720
 721         # result directories
 722         if resudir:
 723             if not os.path.isdir(resudir):
 724                 try:    os.mkdir(resudir)
 725                 except:
 726                     if info>=1: UTMESS('A','RECAL0_82',valk=resudir)
 727                     resudir = None
 728         if not resudir:
 729             # Par defaut, dans un sous-repertoire du repertoire d'execution
 730             pref = 'tmp_macr_recal_'
 731             # On cherche s'il y a un fichier hostfile pour placer les fichiers dans un repertoire partage
 732             l_fr = getattr(prof, 'data')
 733             l_tmp = l_fr[:]
 734             for dico in l_tmp:
 735                if dico['type']=='hostfile':
 736                   pref = os.environ['HOME'] + os.sep + 'tmp_macr_recal_'
 737                   break
 738             # Si batch alors on place les fichiers dans un repertoire partage
 739             if prof['mode'][0]=='batch': pref = os.environ['HOME'] + os.sep + 'tmp_macr_recal_'
 740
 741             resudir = tempfile.mkdtemp(prefix=pref)
 742         flashdir = os.path.join(resudir,'flash')
 743         if info>=1: UTMESS('I','RECAL0_81',valk=resudir)
 744
 745         prof.WriteExportTo( os.path.join(resudir, 'master.export') )
 746
 747         # get hostrc object
 748         hostrc = get_hostrc(run, prof)
 749
 750         # timeout before rejected a job
 751         timeout = get_timeout(prof)
 752
 753
 754         # Ajout des impressions de tables a la fin du .comm
 755         t = []
 756         reponses = calcul
 757         for i in range(len(reponses)):
 758             _ul = str(int(100+i))
 759             num_ul = '99'
 760
 761             # Pour la dynamique la table avec la matrice MAC a un traitement different
 762             if self.DYNAMIQUE:
 763                if ('MAC' in reponses[i][2]):
 764                        t.append( self.ajout_post_mac( reponses[i] ) )
 765
 766             try:    os.remove( tmp_macr_recal+os.sep+"REPE_TABLE"+os.sep+"fort."+_ul )
 767             except: pass
 768
 769             t.append("\n# Recuperation de la table : " + str(reponses[i][0]) + "\n")
 770             t.append("DEFI_FICHIER(UNITE=" + num_ul + ", FICHIER='" + os.path.join('.', 'REPE_OUT', 'fort.'+_ul) + "',);\n" )
 771             t.append("IMPR_TABLE(TABLE="+str(reponses[i][0])+", FORMAT='ASTER', UNITE="+num_ul+", INFO=1, FORMAT_R='E30.20',);\n")
 772             t.append("DEFI_FICHIER(ACTION='LIBERER', UNITE="+num_ul+",);\n")
 773
 774
 775         # number of threads to follow execution
 776         numthread = 1
 777
 778
 779         # ----------------------------------------------------------------------------
 780         # Executions des etudes esclaves
 781         # ----------------------------------------------------------------------------
 782         # ----- Execute calcutions in parallel using a Dispatcher object
 783         # elementary task...
 784         task = DistribParametricTask(run=run, prof=prof, # IN
 785                                      hostrc=hostrc,
 786                                      nbmaxitem=self.NMAX_SIMULT, timeout=timeout,
 787                                      resudir=resudir, flashdir=flashdir,
 788                                      keywords={'POST_CALCUL': '\n'.join(t)},
 789                                      info=info,
 790                                      nbnook=0, exec_result=[])            # OUT
 791         # ... and dispatch task on 'list_tests'
 792         etiq = 'calc_%%0%dd' % (int(log10(nbval)) + 1)
 793         labels = [etiq % (i+1) for i in range(nbval)]
 794         couples = zip(labels, list_val)
 795
 796         if info>=2: print couples
 797         execution = Dispatcher(couples, task, numthread=numthread)
 798
 799         # ----------------------------------------------------------------------------
 800         # Liste des diagnostics
 801         # ----------------------------------------------------------------------------
 802         d_diag = {}
 803         for result in task.exec_result:
 804             #print result
 805             label = result[0]
 806             diag  = result[2]
 807             if len(result) >= 8: output_filename = os.path.join('~', 'flasheur', str(result[7]))
 808             else:                output_filename = ''
 809             d_diag[label] = diag
 810             if not diag[0:2] in ['OK', '<A']:
 811               if not diag in ['<F>_COPY_ERROR']:
 812                   UTMESS('A', 'RECAL0_70', valk=(label, output_filename))
 813
 814                   # Affichage de l'output
 815                   try:
 816                      f=open(output_filename, 'r')
 817                      print f.read()
 818                      f.close()
 819                   except: pass
 820
 821
 822         if not d_diag:
 823                 UTMESS('F', 'RECAL0_71', valk=resudir)
 824         self.list_diag = [ d_diag[label] for label in labels ]
 825
 826         # ----------------------------------------------------------------------------
 827         # Arret si tous les jobs ne se sont pas deroules correctement
 828         # ----------------------------------------------------------------------------
 829         iret = 0
 830         if task.nbnook > 0:
 831            iret = 4
 832         if iret:
 833            UTMESS('A', 'RECAL0_71', valk=resudir)
 834            run.Sortie(iret)
 835
 836
 837
 838         # ----------------------------------------------------------------------------
 839         # Recuperation des tables calculees
 840         # ----------------------------------------------------------------------------
 841         Lcalc = []
 842         i=0
 843         for c in labels:
 844             tbl = get_tables(tables_calc=calcul, tmp_repe_table=os.path.join(resudir, c, 'REPE_OUT'), prof=prof)
 845             Lcalc.append( tbl )  # On stocke sous la forme d'une liste de numpy
 846             i+=1
 847
 848
 849         # ----------------------------------------------------------------------------
 850         # Calcul de la fonctionnelle et du gradient
 851         # ----------------------------------------------------------------------------
 852         if debug: print "AA4/Lcalc=", Lcalc
 853         fonctionnelle, gradient = self.calc2fonc_gradient(Lcalc)
 854
 855
 856         # ----------------------------------------------------------------------------
 857         # Clean result directories
 858         # ----------------------------------------------------------------------------
 859         if clean: shutil.rmtree(resudir, ignore_errors=True)
 860
 861
 862         # ----------------------------------------------------------------------------
 863         # Save all calculated responses
 864         # ----------------------------------------------------------------------------
 865         self.Lcalc = Lcalc
 866
 867         return fonctionnelle, gradient
 868
 869
 870    # ---------------------------------------------------------------------------------------------------------
 871    # ---------------------------------------------------------------------------
 872    def calc2fonc_gradient(self, Lcalc):
 873         """  Calculs de la fonctionnelle et du gradient a partir des tables calculees
 874         """
 875
 876         #print "AA1/Lcalc=", Lcalc
 877
 878         info         = self.info
 879         CalcGradient = self.CalcGradient
 880
 881         # ----------------------------------------------------------------------------
 882         # Recuperation des tables calculees
 883         # ----------------------------------------------------------------------------
 884         seq_FX   = []
 885         seq_FY   = []
 886         seq_DIMS = []
 887         lst_iter = []
 888         for i in range(len(Lcalc)):
 889             tbl = Lcalc[i]
 890             FX = []
 891             FY = []
 892             ldims = []
 893             for array in tbl:
 894                  FX.extend([ x[0] for x in array ])
 895                  FY.extend([ x[1] for x in array ])
 896                  ldims.append(len(array))
 897             # Agregation des resultats
 898             seq_FX.append(FX)
 899             seq_FY.append(FY)
 900             seq_DIMS.append(ldims)
 901             lst_iter.append(i)
 902
 903
 904         # ----------------------------------------------------------------------------
 905         # Fonctionnelle
 906         # ----------------------------------------------------------------------------
 907         # Calcul maitre (point X0)
 908         idx0 = lst_iter.index(0)   # index (les calculs arrivent-ils dans le desordre?)
 909         FY_X0 = seq_FY[idx0]
 910         fonctionnelle = FY_X0
 911
 912
 913         # ----------------------------------------------------------------------------
 914         # Procedure d'assemblage du gradient (une liste de liste)
 915         # ----------------------------------------------------------------------------
 916         gradient = []
 917         if CalcGradient:
 918             for n in range(len(lst_iter))[1:]:
 919                 idx = lst_iter.index(n)
 920                 FY   = seq_FY[idx]
 921                 col = [ (y-x) for x, y in zip(FY, FY_X0) ]
 922                 gradient.append(col)
 923                 #print 'Calcul numero: %s - Diagnostic: %s' % (n, self.list_diag[idx])
 924                 if info>=1: UTMESS('I', 'RECAL0_74', valk=(str(n), self.list_diag[idx]) )
 925
 926         # ----------------------------------------------------------------------------
 927         # Affichages
 928         # ----------------------------------------------------------------------------
 929         if info>=2:
 930             UTMESS('I', 'RECAL0_72', valk=str(fonctionnelle))
 931             import pprint
 932             if CalcGradient:
 933                 UTMESS('I', 'RECAL0_73')
 934                 pprint.pprint(gradient)
 935
 936         return fonctionnelle, gradient
 937
 938
 939    # ---------------------------------------------------------------------------------------------------------
 940    # ---------------------------------------------------------------------------
 941    def find_parameter0(self, content, param):
 942        """
 943        Return the lowest index in content where param is found and
 944        the index of the end of the command.
 945        """
 946        if not self.ASTER_ROOT:
 947            try:    ASTER_ROOT = os.path.join(aster.repout, '..')
 948            except: pass
 949        try:
 950            sys.path.append(os.path.join(ASTER_ROOT, 'ASTK', 'ASTK_SERV', 'lib'))
 951            sys.path.append(os.path.join(ASTER_ROOT, 'lib', 'python%s.%s' % (sys.version_info[0], sys.version_info[1] ) , 'site-packages'))
 952        except: pass
 953        try:
 954            from asrun.utils        import search_enclosed
 955        except Exception, e:
 956            print e
 957            UTMESS('F','RECAL0_99')
 958
 959        pos, endpos = -1, -1
 960        re_start = re.compile('^ *%s *\=' % re.escape(param), re.M)
 961        mat_start = re_start.search(content)
 962        if mat_start is not None:
 963           pos = mat_start.start()
 964           endpos = search_enclosed(content, pos)
 965        return pos, endpos
 966
 967
 968    # ---------------------------------------------------------------------------------------------------------
 969    # ---------------------------------------------------------------------------
 970    def find_parameter(self, content, param):
 971        """
 972        Supprime les parametres du fichier de commande
 973        """
 974        re_start = re.compile('^ *%s *\=' % re.escape(param), re.M)
 975        l=[]
 976        for line in content.split('\n'):
 977           mat_start = re_start.search(line)
 978           if mat_start is None: l.append(line)
 979        return '\n'.join(l)
 980
 981
 982    # ---------------------------------------------------------------------------------------------------------
 983    # ---------------------------------------------------------------------------
 984    def ajout_post_mac(self, reponse):
 985       """
 986          Ajoute un bloc a la fin de l'esclave pour l'affichage des MAC pour l'appariement manuel
 987       """
 988       txt = []
 989       txt.append( "from Macro.reca_mac import extract_mac_array, get_modes, fenetre_mac\n" )
 990       txt.append( "_mac = extract_mac_array("+str(reponse[0])+")\n" )
 991       txt.append( "l_mac=[]\n" )
 992       txt.append( "nb_freq=_mac.shape[1]\n" )
 993       if (self.DYNAMIQUE['APPARIEMENT_MANUEL']=='OUI' and self.graph_mac):
 994           txt.append( "frame =fenetre_mac(" + self.DYNAMIQUE['MODE_EXP']+"," + self.DYNAMIQUE['MODE_CALC']+",_mac)\n" )
 995           txt.append( "list_exp,list_num =frame.get_list()\n" )
 996           txt.append( "for i in range(nb_freq): l_mac.append(_mac[int(list_num[i])-1,int(list_exp[i])-1])\n" )
 997       else:
 998           txt.append( "for i in range(nb_freq): l_mac.append(_mac[i,i])\n" )
 999       txt.append( "DETRUIRE(CONCEPT=_F(NOM="+str(reponse[0])+"),)\n" )
1000       txt.append( str(reponse[0]) + "=CREA_TABLE(LISTE=(_F(PARA='NUME_ORDRE',LISTE_I=range(1,nb_freq+1),),_F(PARA='MAC',LISTE_R=l_mac,),),)\n" )
1001       return '\n'.join(txt)
1002
1003
1004 # --------------------------------------------------------------------------------------------------
1005 # --------------------------------------------------------------------------------------------------
1006 class CALC_ERROR:
1007    """
1008       Classe gérant l'erreur par rapport aux donnees experimentales, la matrice des sensibilites
1009    """
1010    # ---------------------------------------------------------------------------
1011    def __init__(self, experience, X0, calcul, poids=None, objective_type='vector', info=0, unite_resu=None):
1012
1013        if poids is None:
1014             poids = NP.ones(len(experience))
1015        self.experience     = experience
1016        self.X0             = X0
1017        self.calcul         = calcul
1018        self.poids          = poids
1019        self.objective_type = objective_type
1020        self.INFO           = info
1021        self.unite_resu     = unite_resu
1022
1023        from Macro import reca_interp, reca_algo
1024        self.test_convergence   = reca_algo.test_convergence
1025        self.calcul_gradient    = reca_algo.calcul_gradient
1026        self.Simul              = reca_interp.Sim_exp(self.experience, self.poids)
1027        try:    self.Dim   = reca_algo.Dimension(copy.copy(self.X0))
1028        except: self.Dim   = reca_algo.Dimension(copy.copy(self.X0), None)  # gere l'ancienne version de MACR_RECAL
1029        #self.Dim   = reca_algo.Dimension(copy.copy(self.X0))
1030
1031        self.F               = None
1032        self.L_J_init        = None
1033        self.L_J             = None
1034        self.J_init          = None
1035        self.J               = None
1036        self.L_init          = None
1037        self.erreur          = None
1038        self.norme           = None
1039        self.A               = None
1040        self.A_nodim         = None
1041        self.norme_A_        = None
1042        self.norme_A_nodim   = None
1043
1044        if info>=3: self.debug = True
1045        else:       self.debug = False
1046        #if debug: self.debug = True
1047        self.debug = True
1048
1049
1050    # ---------------------------------------------------------------------------
1051    def CalcError(self, Lcalc):
1052
1053        self.F = Lcalc[0]
1054        if self.L_init is None:    self.L_init   = copy.copy(self.F)
1055
1056        self.L_J, self.erreur = self.Simul.multi_interpole(self.F, self.calcul)
1057        if self.L_J_init is None:  self.L_J_init = copy.copy(self.L_J)
1058
1059        self.J = self.Simul.norme_J( copy.copy(self.L_J_init), copy.copy(self.L_J) )
1060        if self.J_init is None:      self.J_init   = copy.copy(self.J)
1061
1062        # norme de l'erreur
1063        self.norme = NP.sum( [x**2 for x in self.erreur] )
1064
1065        if self.debug:
1066            print "AA1/F=", self.F
1067            print "AA1/calcul=", self.calcul
1068            print "AA1/L_J=", self.L_J
1069            print "AA1/erreur=", self.erreur
1070            print "AA1/L_J_init=", self.L_J_init
1071            print "AA1/J=", self.J
1072            print "AA1/norme de l'erreur=", self.norme
1073            print "AA1/norme de J (fonctionnelle)=", str(self.J)
1074
1075        if self.INFO>=1:
1076            UTMESS('I', 'RECAL0_30')
1077
1078        if self.objective_type=='vector':
1079            if self.INFO>=1: UTMESS('I', 'RECAL0_35', valr=self.norme)
1080            return self.erreur
1081        else:
1082            if self.INFO>=1: UTMESS('I', 'RECAL0_36', valr=self.norme)
1083            return self.norme
1084
1085
1086    # ---------------------------------------------------------------------------
1087    def CalcSensibilityMatrix(self, Lcalc, val, dX=None, pas=None):
1088
1089       """
1090          Calcul de F(X0) et de tous les F(X0+h)
1091          Formation de la matrice des sensibilites A
1092          N+1 calculs distribues
1093       """
1094
1095       if not dX and not pas: raise "Need 'dX' or 'pas' parameter."
1096       if     dX and     pas: raise "Need 'dX' or 'pas' parameter, not both."
1097       if pas: dX = len(val)*[pas]
1098       if len(dX) != len(val): raise "Error : 'dX' and 'val' parameters aren't compatible (lenght are not equal).\ndX = %s\nval = %s" % (str(dx), str(val))
1099
1100       reponses  = self.calcul
1101       resu_exp  = self.experience
1102       len_para  = len(val)  # initialement len(self.para)
1103
1104
1105       # Erreur de l'interpolation de F_interp : valeur de F interpolée sur les valeurs experimentales
1106       F = Lcalc[0]
1107       F_interp = self.Simul.multi_interpole_sensib(F, reponses)  #F_interp est une liste contenant des tab num des reponses interpolés
1108
1109
1110       # Creation de la liste des matrices de sensibilités
1111       L_A=[]
1112       for i in range(len(reponses)):
1113           L_A.append(NP.zeros((len(resu_exp[i]),len(val))) )
1114
1115       for k in range(len(val)):   # pour une colone de A (dim = nb parametres)
1116
1117           F_perturbe = Lcalc[k+1]
1118
1119           # Erreur de l'interpolation de F_perturb : valeur de F (perturbée) interpolée sur les valeurs experimentales
1120           F_perturbe_interp = self.Simul.multi_interpole_sensib(F_perturbe, reponses)
1121
1122           # Calcul de L_A (matrice sensibilité des erreurs sur F interpolée)
1123           h = val[k]*dX[k]
1124           for j in range(len(reponses)):
1125              for i in range(len(resu_exp[j])):
1126                 if NP.all(h != 0.):
1127                    L_A[j][i,k] = -1*(F_interp[j][i] - F_perturbe_interp[j][i])/h
1128                 else:
1129                    if self.unite_resu:
1130                        fic=open(os.getcwd()+'/fort.'+str(unite_resu),'a')
1131                        fic.write('\n Probleme de division par zéro dans le calcul de la matrice de sensiblité')
1132                        fic.write('\n Le parametre '+para[k]+'est nul ou plus petit que la précision machine')
1133                        fic.close()
1134                    UTMESS('F','RECAL0_45',valk=para[k])
1135                    return
1136
1137       # On construit la matrice de sensiblité sous forme d'un tab num
1138       dim =[]
1139       for i in range(len(L_A)):
1140          dim.append(len(L_A[i]))
1141       dim_totale = NP.sum(dim)
1142       a=0
1143       self.A_nodim = NP.zeros((dim_totale,len(val)))
1144       for n in range(len(L_A)):
1145          for k in range(len(val)):
1146             for i in range(dim[n]):
1147                self.A_nodim[i+a][k] = L_A[n][i,k]
1148          a=dim[n]
1149
1150       del(L_A)
1151
1152
1153       self.A = self.Dim.adim_sensi( copy.copy(self.A_nodim) )
1154
1155       # Si on n'est pas encore passe par CalcError...
1156       if self.erreur is None:
1157           self.erreur = self.CalcError(Lcalc)
1158       self.gradient_init = self.calcul_gradient(self.A, self.erreur)  #utile pour le test de convergence, on prend les valeurs dimensionnées
1159       self.residu = self.test_convergence(self.gradient_init, self.erreur, self.A, NP.zeros(len(self.gradient_init)))
1160
1161       if self.debug:
1162           print "AA1/erreur=", self.erreur
1163           print "AA1/residu=", self.residu
1164           print "AA1/A_nodim=", self.A_nodim
1165           print "AA1/A=", self.A
1166
1167
1168       if self.objective_type=='vector':
1169           return self.erreur, self.residu, self.A_nodim, self.A
1170       else:
1171           # norme de l'erreur
1172           self.norme = NP.dot(self.erreur, self.erreur)**0.5
1173           self.norme_A_nodim = NP.zeros( (1,len_para))
1174           self.norme_A       = NP.zeros( (1,len_para))
1175           for c in range(len(self.A[0,:])):
1176               norme_A_nodim = 0
1177               norme_A       = 0
1178               for l in range(len(self.A[:,0])):
1179                    norme_A_nodim += self.A_nodim[l,c] * self.A_nodim[l,c]
1180                    norme_A       += self.A[l,c] * self.A[l,c]
1181               self.norme_A_nodim[0,c] = math.sqrt( norme_A_nodim )
1182               self.norme_A[0,c] = math.sqrt( norme_A )
1183           if self.debug:
1184               print "AA1/norme_A_nodim=", self.norme_A_nodim
1185               print "AA1/norme_A=", self.norme_A
1186           return self.erreur, self.residu, self.norme_A_nodim, self.norme_A
1187
1188
1189
1190
1191
1192
1193 # ----------------------------------------------------------------------------------------------------------------------
1194 # ----------------------------------------------------------------------------------------------------------------------
1195 if __name__ == '__main__':
1196
1197     # Execution via YACS ou en externe
1198     isFromYacs = globals().get('ASTER_ROOT', None)
1199
1200
1201     # ------------------------------------------------------------------------------------------------------------------
1202     #                               Execution depuis YACS
1203     # ------------------------------------------------------------------------------------------------------------------
1204     if isFromYacs:
1205         # Execution depuis YACS : les parametres sont deja charges en memoire
1206
1207         # ----------------------------------------------------------------------------
1208         # Parametres courant
1209         X0 = globals().get('X0', [ 80000.,  1000., 30. ])
1210         dX = globals().get('dX', [ 0.001, 0.001, 0.0001])
1211         # ----------------------------------------------------------------------------
1212
1213         # ----------------------------------------------------------------------------
1214         # Parametres
1215         os.environ['ASTER_ROOT'] = ASTER_ROOT
1216         if debug:
1217             clean = False
1218             info  = 1
1219         else:
1220             clean = True
1221             info  = 0
1222         # ----------------------------------------------------------------------------
1223
1224
1225     # ------------------------------------------------------------------------------------------------------------------
1226     #                               Execution en mode EXTERNE
1227     # ------------------------------------------------------------------------------------------------------------------
1228     else:
1229         # Execution en mode EXTERNE : on doit depouiller les parametres de la ligne de commande
1230
1231
1232         from optparse import OptionParser, OptionGroup
1233
1234         p = OptionParser(usage='usage: %s fichier_export [options]' % sys.argv[0])
1235
1236         # Current estimation
1237         p.add_option('--input',             action='store',       dest='input',             type='string',                                       help='Chaine de texte contenant les parametres')
1238         p.add_option('--input_step',        action='store',       dest='input_step',        type='string',                                       help='Chaine de texte contenant les pas de discretisation des differences finies')
1239         p.add_option('--input_file',        action='store',       dest='input_file',        type='string',   default='input.txt',                help='Fichier contenant les parametres')
1240         p.add_option('--input_step_file',   action='store',       dest='input_step_file',   type='string',                                       help='Fichier contenant les pas de discretisation des differences finies')
1241
1242         # Outputs
1243         p.add_option('--output',            action='store',       dest='output',            type='string',   default='output.txt',               help='fichier contenant la fonctionnelle')
1244         p.add_option('--output_grad',       action='store',       dest='output_grad',       type='string',   default='grad.txt',                 help='fichier contenant le gradient')
1245
1246         # Code_Aster installation
1247         p.add_option('--aster_root',        action='store',       dest='aster_root',        type='string',                                       help="Chemin d'installation d'Aster")
1248         p.add_option('--as_run',            action='store',       dest='as_run',            type='string',                                       help="Chemin vers as_run")
1249
1250         # General
1251         p.add_option('--resudir',           action='store',       dest='resudir',           type='string',                                       help="Chemin par defaut des executions temporaires d'Aster")
1252         p.add_option("--noclean",           action="store_false", dest="clean",                              default=True,                       help="Erase temporary Code_Aster execution directory")
1253         p.add_option('--info',              action='store',       dest='info',              type='int',      default=1,                          help="niveau de message (0, [1], 2)")
1254         p.add_option('--sources_root',      action='store',       dest='SOURCES_ROOT',      type='string',                                       help="Chemin par defaut des surcharges Python")
1255         #p.add_option('--slave_computation', action='store',       dest='slave_computation', type='string',   default='distrib',                  help="Evaluation de l'esclave ([distrib], include)")
1256
1257         # MACR_RECAL parameters
1258         p.add_option('--objective',         action='store',       dest='objective',         type='string',   default='fcalc',                    help="Fonctionnelle ([fcalc]/[error])")
1259         p.add_option('--objective_type',    action='store',       dest='objective_type',    type='string',   default='vector',                   help="type de la fonctionnelle (float/[vector])")
1260         p.add_option('--gradient_type',     action='store',       dest='gradient_type' ,    type='string',   default='no',                       help="calcul du gradient par Aster ([no]/normal/adim)")
1261
1262         # MACR_RECAL inputs
1263         p.add_option('--mr_parameters',     action='store',       dest='mr_parameters',     type='string',   default='N_MR_Parameters.py',       help="Fichier de parametres de MACR_RECAL : parametres, calcul, experience")
1264         p.add_option('--study_parameters',  action='store',       dest='study_parameters',  type='string',                                       help="Fichier de parametre de l'etude : export")
1265         p.add_option('--parameters',        action='store',       dest='parameters',        type='string',                                       help="Fichier de parametres")
1266
1267         options, args = p.parse_args()
1268
1269
1270         # Study : .export file
1271         if args: export =  args[0]
1272         else:
1273            liste = glob.glob('*.export')
1274            export = liste[0]
1275         if not os.path.isfile(export): raise "Export file : is missing!"
1276
1277
1278         # Code_Aster installation
1279         ASTER_ROOT = None
1280         if options.aster_root:                  ASTER_ROOT = options.aster_root
1281         elif os.environ.has_key('ASTER_ROOT'):  ASTER_ROOT = os.environ['ASTER_ROOT']
1282         if not ASTER_ROOT: raise "ASTER_ROOT is missing! Set it by --aster_root flag or environment variable ASTER_ROOT"
1283         if not os.path.isdir(ASTER_ROOT): raise "Wrong directory for ASTER_ROOT : %s" % ASTER_ROOT
1284         os.environ['ASTER_ROOT'] = ASTER_ROOT
1285 #         sys.path.append(get_absolute_path(os.path.join(ASTER_ROOT, 'STA10.1', 'bibpyt' )))
1286 #         from Utilitai.Utmess import UTMESS
1287
1288         if options.as_run:          as_run = options.as_run
1289         else:                       as_run = os.path.join(ASTER_ROOT, 'bin', 'as_run')
1290
1291
1292         # General
1293         if options.resudir: resudir = options.resudir
1294         clean = options.clean
1295
1296 #         if   options.info == 0: info = False
1297 #         elif options.info == 1: info = False
1298 #         elif options.info == 2: info = True
1299         info = options.info
1300
1301         # Import des modules supplementaires
1302         if options.SOURCES_ROOT:
1303              if not os.path.isdir(options.SOURCES_ROOT): raise "Wrong directory for sources_root : %s" % options.SOURCES_ROOT
1304              else:
1305                  sys.path.insert(0, options.SOURCES_ROOT)
1306                  sys.path.insert(0, os.path.join(options.SOURCES_ROOT, 'sources'))
1307
1308
1309         # MACR_RECAL inputs
1310         if options.mr_parameters:
1311             try:
1312                 if info>=1: print "Read MR parameters file : %s" % options.mr_parameters
1313                 execfile(options.mr_parameters)
1314             except: raise "Wrong file for MR Parameters: %s" % options.mr_parameters
1315         else: raise "MR Parameters file needed ! Use --mr_parameters flag"
1316         parametres = globals().get('parametres',  None)
1317         calcul     = globals().get('calcul',      None)
1318         experience = globals().get('experience',  None)
1319         poids      = globals().get('poids',       None)
1320
1321         if not parametres:  raise "MR Parameters file need to define 'parametres' variable"
1322         if not calcul:      raise "MR Parameters file need to define 'calcul' variable"
1323         if type(parametres)  != list: raise "Wrong type for 'parametres' variable in MR parameters file : %s"  % options.mr_parameters
1324         if type(calcul)      != list: raise "Wrong type for 'calcul' variable in MR parameters file : %s"      % options.mr_parameters
1325
1326         if options.objective == 'error':
1327              if type(experience) != list: raise "For error objective output, the 'experience' variable must be a list of arrays"
1328              if type(poids) not in [list, tuple, NP.ndarray]: raise "The 'poids' variable must be a list or an array"
1329              if len(poids) != len(experience): raise "'experience' and 'poids' lists must have the same lenght"
1330
1331
1332         # MACR_RECAL parameters
1333         objective      = options.objective
1334         objective_type = options.objective_type
1335         gradient_type  = options.gradient_type
1336
1337
1338         # X0 : read from commandline flag or from file
1339         if not os.path.isfile(options.input_file): options.input_file = None
1340         if not (options.input or  options.input_file): raise "Missing input parameters"
1341         if     (options.input and options.input_file): raise "Error : please use only one choice for input parameters definition"
1342
1343         if options.input_file:
1344             try:
1345                 f = open(options.input_file, 'r')
1346                 options.input = f.read()
1347                 f.close()
1348             except:
1349                 raise "Can't read input parameters file : %s" % options.input_file
1350
1351         # Extract X0 from text
1352         try:
1353             txt = options.input.strip()
1354             txt = txt.replace(',', ' ')
1355             txt = txt.replace(';', ' ')
1356             X0 = [ float(x) for x in txt.split() ]
1357             if type(X0) != list: raise "Wrong string for input parameters : %s" % options.input
1358         except:
1359             raise "Can't decode input parameters string : %s.\n It should be a comma separated list." % options.input
1360
1361
1362         # dX : read from commandline flag or from file
1363         dX = None
1364         if options.gradient_type == 'no':
1365             if (options.input_step or  options.input_step_file): raise "You must set 'gradient_type' to another choice than 'no' or remove input step parameters from commandline"
1366         else:
1367             if not (options.input_step or  options.input_step_file): raise "Missing input step parameters"
1368             if     (options.input_step and options.input_step_file): raise "Error : please use only one choice for input step parameters definition"
1369
1370             if options.input_step_file:
1371                 try:
1372                     f = open(options.input_step_file, 'r')
1373                     options.input_step = f.read()
1374                     f.close()
1375                 except:
1376                     raise "Can't read file for discretisation step : %s" % options.input_step_file
1377
1378             # Extract dX from text
1379             try:
1380                 txt = options.input_step.strip()
1381                 txt = txt.replace(',', ' ')
1382                 txt = txt.replace(';', ' ')
1383                 dX = [ float(x) for x in txt.split() ]
1384                 if type(dX) != list: raise "Wrong string for discretisation step : %s" % options.input_step
1385             except:
1386                 raise "Can't decode input parameters string : %s.\n It should be a comma separated list." % options.input_step
1387
1388
1389
1390
1391     # ------------------------------------------------------------------------------------------------------------------
1392     #                               Execution des calculs (N+1 calculs distribues si dX est fourni)
1393     # ------------------------------------------------------------------------------------------------------------------
1394
1395     # Repertoire contenant les resultats des calculs Aster (None = un rep temp est cree)
1396     resudir = globals().get('resudir', None)
1397
1398     # Affichage des parametres
1399     lpara = [x[0] for x in parametres]
1400     lpara.sort()
1401     if info >=1:
1402        lpara = [x[0] for x in parametres]
1403        lpara.sort()
1404        print "Calcul avec les parametres : \n%s" % Affiche_Param(lpara, X0)
1405
1406     C = CALCULS_ASTER(
1407                 # MACR_RECAL inputs
1408                 parametres          = parametres,
1409                 calcul              = calcul,
1410                 experience          = experience,
1411                      )
1412
1413     fonctionnelle, gradient = C.run(
1414                 # Current estimation
1415                 X0                  = X0,
1416                 dX                  = dX,
1417
1418                 # Code_Aster installation
1419                 ASTER_ROOT          = ASTER_ROOT,
1420                 as_run              = as_run,
1421
1422                 # General
1423                 resudir             = resudir,
1424                 clean               = clean,
1425                 info                = info,
1426
1427                 # Study
1428                 export              = export,
1429
1430 #                 # MACR_RECAL inputs
1431 #                 parametres          = parametres,
1432 #                 calcul              = calcul,
1433 #                 experience          = experience,
1434     )
1435
1436     # ------------------------------------------------------------------------------------------------------------------
1437     #                               Calcul de l'erreur par rapport aux donnees experimentale
1438     # ------------------------------------------------------------------------------------------------------------------
1439     if not isFromYacs:        # Execution en mode EXTERNE uniquement
1440
1441         # Calcul de l'erreur par rapport aux donnees experimentale
1442         if objective == 'error':
1443             E = CALC_ERROR(
1444                 experience          = experience,
1445                 X0                  = X0,
1446                 calcul              = calcul,
1447                 poids               = poids,
1448                 objective_type      = objective_type,
1449                 info=info,
1450             )
1451
1452             erreur                      = E.CalcError(C.Lcalc)
1453             erreur, residu, A_nodim, A  = E.CalcSensibilityMatrix(C.Lcalc, X0, dX=dX, pas=None)
1454
1455             fonctionnelle = erreur
1456             if   gradient_type == 'normal': gradient = A
1457             elif gradient_type == 'adim':   gradient = A_nodim
1458             else: raise "??"
1459
1460
1461
1462     # ------------------------------------------------------------------------------------------------------------------
1463     #                               Ecriture des resultats
1464     # ------------------------------------------------------------------------------------------------------------------
1465     if not isFromYacs:        # Execution en mode EXTERNE uniquement
1466
1467         # Fonctionnelle
1468         if options.objective_type == 'float':
1469            fonctionnelle = math.sqrt( NP.sum( [x**2 for x in fonctionnelle] ) )
1470         Ecriture_Fonctionnelle(output_file=options.output, type_fonctionnelle=options.objective_type, fonctionnelle=fonctionnelle)
1471
1472         # Gradient
1473         if gradient: Ecriture_Derivees(output_file=options.output_grad, derivees=gradient)
1474
1475
1476
1477     # ------------------------------------------------------------------------------------------------------------------
1478     #                               Affichages
1479     # ------------------------------------------------------------------------------------------------------------------
1480     if info>=2:
1481         print "\nFonctionnelle au point X0: \n%s" % str(fonctionnelle)
1482         import pprint
1483         if dX:
1484            print "\nGradient au point X0:"
1485            pprint.pprint(gradient)