From e7621e95dcd542b4ff15c26afa26abdac5b59da8 Mon Sep 17 00:00:00 2001 From: Pascale Noyret Date: Wed, 29 Nov 2006 11:17:32 +0000 Subject: [PATCH] *** empty log message *** --- Aster/Cata/Utilitai/Sensibilite.py | 78 +++ Aster/Cata/Utilitai/System.py | 217 ++++++ Aster/Cata/Utilitai/liss_enveloppe.py | 931 ++++++++++++++++++++++++++ Aster/Cata/Utilitai/optimize.py | 593 ++++++++++++++++ 4 files changed, 1819 insertions(+) create mode 100644 Aster/Cata/Utilitai/Sensibilite.py create mode 100644 Aster/Cata/Utilitai/System.py create mode 100644 Aster/Cata/Utilitai/liss_enveloppe.py create mode 100644 Aster/Cata/Utilitai/optimize.py diff --git a/Aster/Cata/Utilitai/Sensibilite.py b/Aster/Cata/Utilitai/Sensibilite.py new file mode 100644 index 00000000..9d0ad2b0 --- /dev/null +++ b/Aster/Cata/Utilitai/Sensibilite.py @@ -0,0 +1,78 @@ +#@ MODIF Sensibilite Utilitai DATE 07/03/2006 AUTEUR MCOURTOI M.COURTOIS +# -*- coding: iso-8859-1 -*- +# CONFIGURATION MANAGEMENT OF EDF VERSION +# ====================================================================== +# COPYRIGHT (C) 1991 - 2006 EDF R&D WWW.CODE-ASTER.ORG +# THIS PROGRAM IS FREE SOFTWARE; YOU CAN REDISTRIBUTE IT AND/OR MODIFY +# IT UNDER THE TERMS OF THE GNU GENERAL PUBLIC LICENSE AS PUBLISHED BY +# THE FREE SOFTWARE FOUNDATION; EITHER VERSION 2 OF THE LICENSE, OR +# (AT YOUR OPTION) ANY LATER VERSION. +# +# THIS PROGRAM IS DISTRIBUTED IN THE HOPE THAT IT WILL BE USEFUL, BUT +# WITHOUT ANY WARRANTY; WITHOUT EVEN THE IMPLIED WARRANTY OF +# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. SEE THE GNU +# GENERAL PUBLIC LICENSE FOR MORE DETAILS. +# +# YOU SHOULD HAVE RECEIVED A COPY OF THE GNU GENERAL PUBLIC LICENSE +# ALONG WITH THIS PROGRAM; IF NOT, WRITE TO EDF R&D CODE_ASTER, +# 1 AVENUE DU GENERAL DE GAULLE, 92141 CLAMART CEDEX, FRANCE. +# ====================================================================== + +""" +""" + +from types import StringTypes +import aster +from Utilitai.Utmess import UTMESS + +# Doit etre en accord avec semeco.f +prefix = '&NOSENSI.MEMO' +nommem = '%-24s' % (prefix + '.CORR') + +def NomCompose(nomsd, nomps, msg='A'): + """Utilitaire analogue à la routine fortran PSRENC. + Retourne le nom composé à partir du couple (SD de base, paramètre sensible). + `msg` : 'A', 'F' ou 'silence' (pas de message) + """ + nomcomp = None + vect = aster.getvectjev(nommem) + if not type(nomsd) in StringTypes: + nomsd = nomsd.get_name() + if not type(nomps) in StringTypes: + nomps = nomps.get_name() + if vect: + trouv = False + for ch in vect[0:len(vect):2]: + if ch[0:8].strip() == nomsd and ch[8:16].strip() == nomps: + trouv=True + nomcomp = ch[16:24].strip() + if not trouv and msg != 'silence': + UTMESS(msg, 'NomCompose', 'Dérivée de %s par rapport à %s non disponible'\ + % (nomsd, nomps)) + elif msg != 'silence': + UTMESS(msg, 'NomCompose', 'Pas de calcul de sensibilité accessible.') + return nomcomp + +def SdPara(nomcomp, msg='A'): + """Retourne le couple (SD de base, paramètre sensible) correspondant au nom + composé `nomcomp`. + `msg` : 'A', 'F' ou 'silence' (pas de message) + """ + nomsd = None + nomps = None + vect = aster.getvectjev(nommem) + if not type(nomcomp) in StringTypes: + UTMESS('F', 'SdPara', "Argument de type '%s' invalide" % type(nomcomp).__name__) + if vect: + trouv = False + for ch in vect[0:len(vect):2]: + if ch[16:24].strip() == nomcomp: + trouv = True + nomsd = ch[0:8].strip() + nomps = ch[8:16].strip() + if not trouv and msg != 'silence': + UTMESS(msg, 'SdPara', 'Dérivée de %s par rapport à %s non disponible'\ + % (nomsd, nomps)) + elif msg != 'silence': + UTMESS(msg, 'SdPara', 'Pas de calcul de sensibilité accessible.') + return nomsd, nomps diff --git a/Aster/Cata/Utilitai/System.py b/Aster/Cata/Utilitai/System.py new file mode 100644 index 00000000..9ecb9b66 --- /dev/null +++ b/Aster/Cata/Utilitai/System.py @@ -0,0 +1,217 @@ +#@ MODIF System Utilitai DATE 29/08/2006 AUTEUR MCOURTOI M.COURTOIS +# -*- coding: iso-8859-1 -*- +# CONFIGURATION MANAGEMENT OF EDF VERSION +# ====================================================================== +# COPYRIGHT (C) 1991 - 2006 EDF R&D WWW.CODE-ASTER.ORG +# THIS PROGRAM IS FREE SOFTWARE; YOU CAN REDISTRIBUTE IT AND/OR MODIFY +# IT UNDER THE TERMS OF THE GNU GENERAL PUBLIC LICENSE AS PUBLISHED BY +# THE FREE SOFTWARE FOUNDATION; EITHER VERSION 2 OF THE LICENSE, OR +# (AT YOUR OPTION) ANY LATER VERSION. +# +# THIS PROGRAM IS DISTRIBUTED IN THE HOPE THAT IT WILL BE USEFUL, BUT +# WITHOUT ANY WARRANTY; WITHOUT EVEN THE IMPLIED WARRANTY OF +# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. SEE THE GNU +# GENERAL PUBLIC LICENSE FOR MORE DETAILS. +# +# YOU SHOULD HAVE RECEIVED A COPY OF THE GNU GENERAL PUBLIC LICENSE +# ALONG WITH THIS PROGRAM; IF NOT, WRITE TO EDF R&D CODE_ASTER, +# 1 AVENUE DU GENERAL DE GAULLE, 92141 CLAMART CEDEX, FRANCE. +# ====================================================================== + + +"""Ce module définit la classe `SYSTEM` et la fonction `ExecCommand` +qui est présente uniquement pour commodité pour les Macros. + +La classe SYSTEM est semblable à celle utilisée dans ASTK_SERV. +""" + +__all__ = ["SYSTEM", "ExecCommand"] + +import sys +import os +import popen2 +import re +from sets import Set +from types import FileType + +# ----- differ messages translation +def _(mesg): + return mesg + +#------------------------------------------------------------------------------- +def _exitcode(status, default=0): + """Extrait le code retour du status. Retourne `default` si le process + n'a pas fini pas exit. + """ + iret = default + if os.WIFEXITED(status): + iret = os.WEXITSTATUS(status) + return iret + +#------------------------------------------------------------------------------- +#------------------------------------------------------------------------------- +#------------------------------------------------------------------------------- +class SYSTEM: + """Class to encapsultate "system" commands (this a simplified version of + ASTER_SYSTEM class defined in ASTK_SERV part). + """ + # this value should be set during installation step. + MaxCmdLen = 1024 + # line length -9 + _LineLen = 80-9 + +#------------------------------------------------------------------------------- + def __init__(self, **kargs): + """Initialization. + Optionnal arguments : silent, verbose, debug, cc_files, maxcmdlen. + """ + self.verbose = kargs.get('verbose', True) + self.debug = kargs.get('debug', False) + self.cc_files = kargs.get('cc_files', None) + if kargs.has_key('maxcmdlen'): + self.MaxCmdLen = kargs['maxcmdlen'] + +#------------------------------------------------------------------------------- + def _mess(self, msg, cod=''): + """Just print a message + """ + self._print('%-18s %s' % (cod, msg)) + +#------------------------------------------------------------------------------- + def _print(self, *args, **kargs): + """print replacement. + Optionnal argument : + term : line terminator (default to os.linesep). + """ + term = kargs.get('term', os.linesep) + files = Set([sys.stdout]) + if self.cc_files: + files.add(self.cc_files) + for f in files: + if type(f) is FileType: + txt = ' '.join(['%s'%a for a in args]) + f.write(txt.replace(os.linesep+' ', os.linesep)+term) + f.flush() + else: + print _('FileType object expected : %s / %s') % (type(f), repr(f)) + +#------------------------------------------------------------------------------- + def VerbStart(self, cmd, verbose=None): + """Start message in verbose mode + """ + Lm = self._LineLen + if verbose == None: + verbose = self.verbose + if verbose: + pcmd = cmd + if len(cmd) > Lm-2 or cmd.count('\n') > 0: + pcmd = pcmd+'\n'+' '*Lm + self._print(('%-'+str(Lm)+'s') % (pcmd,), term='') + +#------------------------------------------------------------------------------- + def VerbEnd(self, iret, output='', verbose=None): + """Ends message in verbose mode + """ + if verbose == None: + verbose = self.verbose + if verbose: + if iret == 0: + self._print('[ OK ]') + else: + self._print(_('[FAILED]')) + self._print(_('Exit code : %d') % iret) + if (iret != 0 or self.debug) and output: + self._print(output) + +#------------------------------------------------------------------------------- + def VerbIgnore(self, verbose=None): + """Ends message in verbose mode + """ + if verbose == None: + verbose = self.verbose + if verbose: + self._print(_('[ SKIP ]')) + +#------------------------------------------------------------------------------- + def Shell(self, cmd, bg=False, verbose=None, follow_output=False, + alt_comment=None, interact=False): + """Execute a command shell + cmd : command + bg : put command in background if True + verbose : print status messages during execution if True + follow_output : follow interactively output of command + alt_comment : print this "alternative comment" instead of "cmd" + interact : allow the user to interact with the process + (don't close stdin). bg=True implies interact=False. + Return : + iret : exit code if bg = False, + process id if bg = True + output : output lines (as string) + """ + if not alt_comment: + alt_comment = cmd + if verbose == None: + verbose = self.verbose + if bg: + interact = False + if len(cmd) > self.MaxCmdLen: + self._mess((_('length of command shell greater '\ + 'than %d characters.') % self.MaxCmdLen), _('_ALARM')) + if self.debug: + self._print(' ', cmd) + self._print(' background mode : ', bg) + # exec + self.VerbStart(alt_comment, verbose=verbose) + if follow_output and verbose: + self._print(_('\nCommand output :')) + # run interactive command + if interact: + iret = os.system(cmd) + return _exitcode(iret), '' + # use popen to manipulate stdout/stderr + output = [] + p = popen2.Popen4(cmd) + p.tochild.close() + if not bg: + if not follow_output: + output = p.fromchild.readlines() + else: + while p.poll() == -1: + output.append(p.fromchild.readline()) + # \n already here... + self._print(output[-1], term='') + # to be sure to empty the buffer + end = p.fromchild.readlines() + self._print(''.join(end)) + output.extend(end) + iret = _exitcode(p.wait()) + else: + iret = 0 + p.fromchild.close() + output = ''.join(output) + + # repeat header message + if follow_output: + self.VerbStart(alt_comment, verbose=verbose) + mat = re.search('EXIT_CODE=([0-9]+)', output) + if mat: + iret = int(mat.group(1)) + self.VerbEnd(iret, output, verbose=verbose) + if bg: + iret = p.pid + if verbose: + self._print(_('Process ID : '), iret) + return iret, output + +#------------------------------------------------------------------------------- +# Juste par commodité. +system = SYSTEM() +ExecCommand = system.Shell + + +#------------------------------------------------------------------------------- +#------------------------------------------------------------------------------- +#------------------------------------------------------------------------------- +if __name__ == '__main__': + iret, output = ExecCommand('ls', alt_comment='Lancement de la commande...') + diff --git a/Aster/Cata/Utilitai/liss_enveloppe.py b/Aster/Cata/Utilitai/liss_enveloppe.py new file mode 100644 index 00000000..0545b1d5 --- /dev/null +++ b/Aster/Cata/Utilitai/liss_enveloppe.py @@ -0,0 +1,931 @@ +#@ MODIF liss_enveloppe Utilitai DATE 29/08/2005 AUTEUR THOMASSO D.THOMASSON +# -*- coding: iso-8859-1 -*- +# CONFIGURATION MANAGEMENT OF EDF VERSION +# ====================================================================== +# COPYRIGHT (C) 1991 - 2005 EDF R&D WWW.CODE-ASTER.ORG +# THIS PROGRAM IS FREE SOFTWARE; YOU CAN REDISTRIBUTE IT AND/OR MODIFY +# IT UNDER THE TERMS OF THE GNU GENERAL PUBLIC LICENSE AS PUBLISHED BY +# THE FREE SOFTWARE FOUNDATION; EITHER VERSION 2 OF THE LICENSE, OR +# (AT YOUR OPTION) ANY LATER VERSION. +# +# THIS PROGRAM IS DISTRIBUTED IN THE HOPE THAT IT WILL BE USEFUL, BUT +# WITHOUT ANY WARRANTY; WITHOUT EVEN THE IMPLIED WARRANTY OF +# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. SEE THE GNU +# GENERAL PUBLIC LICENSE FOR MORE DETAILS. +# +# YOU SHOULD HAVE RECEIVED A COPY OF THE GNU GENERAL PUBLIC LICENSE +# ALONG WITH THIS PROGRAM; IF NOT, WRITE TO EDF R&D CODE_ASTER, +# 1 AVENUE DU GENERAL DE GAULLE, 92141 CLAMART CEDEX, FRANCE. +# ====================================================================== +""" + Maquette demande SEPTEN fonction de lissage enveloppe + Les données se présentent sous la forme d'un fichier texte comportant + un ensemble de groupe de lignes organisé comme suit : + - ligne 1 : Informations générales + - ligne 2 : une liste de valeur d'amortissement + - lignes 3...n : une liste de valeur commencant par une frequence suivit + des amplitudes du spectre pour chacun des amortissements + liste en ligne 2 + Points importants : + - Le nombre de lignes définissant le spectre peut varier + - Le nombre de valeur d'amortissement peut varier ? + + ==> On propose d'opérer ligne par ligne + ==> L'absence d'informations sur la variabilité du nombre d'éléments oblige à traiter le cas général + + + + Etapes du développement : + 24/05/2005 : Test de lecture du fichier, choix d'une stratégie de gestion + 25/05/2005 : Objet itérable pour la lecture du fichier + 29/05/2005 : Créations de filtres pour les spectres +""" + +import math + +def nearestKeys(k1, dct) : + """ + retourne les clés (doublet) les plus proches de 'key' dans le dictionnaire dct + par valeur inférieure et supérieures + """ + kr = min(dct.keys()) + for k2 in dct.keys() : + if (k2kr) : kr = k2 + kinf = kr + + kr = max(dct.keys()) + for k2 in dct.keys() : + if (k2>k1) and (k2 self.upperBound : + toDel = toDel + [i] + + # Nettoyage des fréquences à suppimer (on commence par les plus hautes) + for i in toDel[::-1] : + del spr.listFreq[i] + del spr.dataVal[i] + + toDel = [] + for i in range(0, len(spr.listFreq)) : + if spr.listFreq[i] < self.lowerBound : + toDel = toDel + [i] + else : + break + + # Nettoyage des fréquences à suppimer (on finit par les plus basses) + for i in toDel[::-1] : + del spr.listFreq[i] + del spr.dataVal[i] + + return spr + +class filtreCrible(filtre): + """ + Criblage du spectre selon specif SEPTEN §C-5 (ce que j'en comprend) + """ + def __init__(self, **listOpt): + try : + self.tolerance = listOpt['tolerance'] + except KeyError : + self.tolerance = 0.25 + + self.listEtats = [] + + def _filtre(self, sp) : + self._initListeEtats(sp) # Création de la table des étsts des valeurs du spectre + coef = 1 + + # Parcours de la liste des fréquences + i1, i2, i3 = 0, 2, 1 + bTest = True + while True : + try : + bTest = self._amplitude(sp, i1, i2, i3, coef) + if not(bTest) and ((i2-i1) > 2) : + # Le point a été éliminé, on réexamine le point précédent sauf si c'est le premier examiné + i3 -= 1 + if self._amplitude(sp, i1, i2, i3, coef) : + # Le point a été "récupéré", il devient la nouvelle origine + i1 = i3 + i2 = i2 # écrit quand meme pour la compréhension + i3 += 1 + else : + # Le point reste désactivé, on avance au point suivant, le point d'origine est conservé + i1 = i1 + i2 += 1 + i3 += 2 + elif not(bTest) and not((i2-i1) > 2) : + i1 = i1 + i2 += 1 + i3 += 1 + else : # Le point est conservé, il devient la nouvelle origine + i1 = i3 + i2 += 1 + i3 += 1 + except IndexError : + break + + return self._crible(sp) + + def _initListeEtats(self, sp) : + """ + Crée une liste associant à chaque fréquence du spectre passé en paramètre, un état booléen + qui spécifie si ce couple fréquence-valeur est supprimé ou pas + NB : au départ toutes les valeur sont "True" car aucune valeur n'a été supprimée + """ + self.listEtats = [True for x in sp.listFreq] + + def _crible(self, sp) : + """ + Supprime les points de fréquence qui sont marqué False dans listEtats + """ + sp2 = spectre([], []) # On force car il y a un problème de persistance su spectre précédent + for x,y,z in zip(self.listEtats, sp.listFreq, sp.dataVal) : + if x : + sp2.listFreq.append(y) + sp2.dataVal.append(z) + + return sp2 + + def _amplitude(self, sp, id1, id2, id3, coef=1) : + """ + teste le point d'indice id3 par rapport aux points à sa gauche(p1 d'indice id1) et + à sa droite (p2 d'indice id2). + Le point est éliminé si sa valeur est en dessous de la droite reliant les points + d'indice id1 et id2 sauf si sa distance à cette droite est supérieure à : + tolerance*ordonnée + Le critère est purement sur l'amplitude du point indépendemment de l'intervalle + sur lequel il s'applique + """ + x0 = sp.listFreq[id1] + y0 = sp.dataVal[id1] + x1 = sp.listFreq[id2] + y1 = sp.dataVal[id2] + x2 = sp.listFreq[id3] + y2 = sp.dataVal[id3] + + yp2 = interpole(x2, x0, y0, x1, y1) + + # Le point est il susceptible d'etre supprimé (est il en dessous de la droite p1-p2 ?) + # Faut-il le supprimer pour autant (distance y2 à yp2 > tolerance% de y2) + bSup = not((y2 < yp2) and (abs(yp2-y2)/y2 < self.tolerance)) + + # Changement de l'état du point + self.listEtats[id3] = bSup + + return bSup + +class filtreChevauchement(filtre): + """ + Compare un spectre à un spectre de référence fréquence par fréquence. + Si une fréquence n'existe pas, on cherche la valeur équivalent par interpolation + Pour éviter tout recouvrement, il est éventuellement nécessaire de rajouter + des informations à certaines fréquences + """ + def __init__(self, **listOpt) : + try : + self.spRef = listOpt['ref'] + except KeyError : + self.spRef = spectre() + + try : + signe = listOpt['ordre'] + self.ordre = signe/abs(signe) #coefficient +1 ou -1 + except KeyError : + self.ordre = +1 + except ZeroDivisionError : + self.ordre = +1 + + def _filtre(self, sp) : + spDict = sp.buildMap() + spRefDict = self.spRef.buildMap() + spTestDict = {} + + # On commence par construire un dictionnaire des valeurs à tester comportant toutes les clés contenues + for k in spDict.keys() : spTestDict[k] = True + for k in spRefDict.keys() : spTestDict[k] = True + + # On teste ensuite toutes les valeurs du dictionnaire + for k in spTestDict.keys() : + # Test d'existence dans le dictionnaire du spectre de référence + try : + vr = spRefDict[k] + except KeyError : + ki = nearestKeys(k, spRefDict) + vr = interpole(k, ki[0], spRefDict[ki[0]], ki[1], spRefDict[ki[1]]) + # Test d'existence dans le dictionnaire du spectre à tester + try : + vt = spDict[k] + except KeyError : + ki = nearestKeys(k, spDict) + vt = interpole(k, ki[0], spDict[ki[0]], ki[1], spDict[ki[1]]) + + # Comparaison des deux valeurs. La clé est ajoutée si elle n'existe pas + if vt*self.ordre < vr*self.ordre : spDict[k] = vr + + return spectre.sortSpectre(spDict) + +class spectre : + """ + décrit un spectre composé d'un ensemble de résultat associé à un ensemble de fréquence + """ + def __init__(self, listFreq = [], dataVal = []) : + self.listFreq = [v for v in listFreq] + self.dataVal = [v for v in dataVal] + + def filtre(self, fi) : + """ + Applique le filtre passé en paramètre au spectre et retourne un nouveau spectre + """ + return fi(self) + + def __staticSortSpectre(dict) : + """ + Convertit un spectre présenté sous forme d'un dictionnaire en un spectre normal + Fonction créé parceque les clés du dictionnaire ne sont pas ordonnées + """ + lstFrq = dict.keys() + lstFrq.sort() + lstVal = [] + for fr in lstFrq : + try : + lstVal.append(dict[fr]) + except KeyError : # Ne devrait jamais arriver + lstVal.append(-1E15) + + return spectre(lstFrq, lstVal) + + sortSpectre = staticmethod(__staticSortSpectre) # définition en tant que méthode statique + + def getCoupleVal(self,indice) : + return (self.listFreq[indice], self.dataVal[indice]) + + def moyenne(self) : + """ + Calcule la moyenne pondéré : somme(An* dfn) /F + """ + somme = 0.0 + X0 = self.listFreq[0] + X1 = self.listFreq[len(self.listFreq)-1] + for i in range(0,len(self.listFreq)-1) : + x0 = self.listFreq[i] + y0 = self.dataVal[i] + x1 = self.listFreq[i+1] + y1 = self.dataVal[i+1] + + somme = somme + (y0+y1) * abs(x1-x0) / 2 + + return somme/abs(X1-X0) + + def seuil(self, limit=75) : + """ + retourne un couple d'index délimitant l'ensemble des valeurs du spectre + définissant "limit" pourcent du total cumulé des valeurs + [borne à gauche inclue, borne à droite exclue[ + ATTENTION on fait l'hypothèse que le spectre a une forme en cloche. + """ + resu = [0 for v in self.dataVal] # initialisation du tableau resultat + + maxi = max(self.dataVal) # Valeur maxu du spectre + iMaxi = self.dataVal.index(maxi) # Index de la valeur max du spectre + + # ETAPE 1 : SOMMATION + somme = 0.0 + for v, i in zip(self.dataVal[0:iMaxi], range(0, iMaxi)) : + somme = somme + v + resu[i] = somme + + somme = 0.0 + for v, i in zip(self.dataVal[:iMaxi:-1], range(len(self.dataVal)-1, iMaxi, -1)) : + somme = somme + v + resu[i] = somme + + resu[iMaxi] = resu[iMaxi-1] + self.dataVal[iMaxi] + resu[iMaxi+1] + + #ETAPE 2 : POURCENTAGE (PAS NECESSAIRE MAIS PLUS LISIBLE) + for v, i in zip(self.dataVal[0:iMaxi], range(0, iMaxi)) : + resu[i] = (resu[i] + maxi/2) / resu[iMaxi] * 100 + + for v, i in zip(self.dataVal[iMaxi+1:], range(iMaxi+1, len(self.dataVal))) : + resu[i] = (resu[i] + maxi/2) / resu[iMaxi] * 100 + + resu[iMaxi] = resu[iMaxi-1] + resu[iMaxi+1] + + # ETAPE 3 : RECHERCHE DES BORNES + limit = (100.0 - limit) / 2.0 + b1 = b2 = True + for v1, v2 in zip(resu[:], resu[::-1]): # Parcours simultané dans les deux sens + if b1 and v1 >= limit : # Borne à gauche trouvée + i1 = resu.index(v1) + b1 = False + if b2 and v2 >= limit : # Borne à droite trouvée + i2 = resu.index(v2) + 1 # Borne à droit exclue de l'intervalle + b2 = False + + return (i1, i2) + + def cut(self, nuplet) : + """ + Découpe un spectre en sous-spectres qui sont retournés en sortie de la fonction + sous la forme d'un tableau de spectres + """ + # transformation du nuplet en tableau (permet de lui ajouter un élément) + tabNuplet = [v for v in nuplet] + tabNuplet.append(len(self.listFreq)) + + # Traitement + tableRes = list() + bGauche = 0 + for borne in tabNuplet : + bDroite = borne + sp = spectre() + for i in range(bGauche, bDroite) : + sp.listFreq.append(self.listFreq[i]) + sp.dataVal.append(self.dataVal[i]) + + tableRes.append(sp) + bGauche = bDroite + + return tableRes + + def __staticMerge(tabSpectre) : + """ + A l'inverse de la fonction cut, construit un seul spectre à partir d'un ensemble de spectres + """ + # On vérifie d'abord que les spectres ne partagent pas la meme bande de fréquence (fut ce partiellement) + for i in range(0, len(tabSpectre)-1) : + if exclus(tabSpectre[i].listFreq, tabSpectre[i+1].listFreq) : raise SpectreError + if exclus(tabSpectre[0].listFreq, tabSpectre[len(tabSpectre)-1].listFreq) : raise SpectreError + + spRes = spectre() + #cumul des spectres + for sp in tabSpectre : + for f, v in zip(sp.listFreq, sp.dataVal) : + spRes.listFreq.append(f) + spRes.dataVal.append(v) + + return spRes + + merge = staticmethod(__staticMerge) # définition en tant que méthode statique + + def buildMap(self) : + """ + Construit un dictionnaire à partir d'un spectre + """ + dict = {} + for i, j in zip(self.listFreq, self.dataVal) : + dict[i] = j + + return dict + +class nappe : + """ + décrit un objet nappe qui associe à un ensemble de fréquence à une enesmble de résultats + """ + def __init__(self, listFreq = [], listeTable = [], listAmor = [], entete = ""): + self.listFreq = [v for v in listFreq] # recopie physique ! + self.listTable = [list() for t in listeTable] + for t, st in zip(listeTable, self.listTable) : + for v in t : st.append(v) + + self.listAmor = [l for l in listAmor] + self.entete = entete + + def __staticBuildFromListSpectre(lsp) : + """ + Construction d'une nappe à partir d'une liste de spectres + """ + # On commence par vérifier que toutes les nappes on la meme base de fréquences + # A inclus dans B inclus dans C inclus dans .... et DERNIER inclus dans PREMIER ==> tous égaux + for i in range(0,len(lsp.listSp)-1) : + if inclus(lsp.listSp[i].listFreq, lsp.listSp[i+1].listFreq) : raise NappeCreationError + if inclus(lsp.listSp[i+1].listFreq, lsp.listSp[0].listFreq) : raise NappeCreationError + + # Construction de la nappe à proprement parler + listeFreq = [fr for fr in lsp.listSp[0].listFreq] + listeTable = [list() for sp in lsp.listSp] + for sp, lv in zip(lsp.listSp, listeTable) : + for v in sp.dataVal : + lv.append(v) + return nappe(listeFreq, listeTable, [], 'toto') + + buildFromListSpectre = staticmethod(__staticBuildFromListSpectre) # définition en tant que méthode statique + + def getNbSpectres(self) : + """ Retourne le nombre d'éléments dans la nappe """ + return len(self.listAmor) + + def getNbFreq(self) : + """ Retourne le nombre d'éléments dans chaque spectre """ + return len(self.listFreq) + + def getSpectre(self, index) : + """ + Retourne le spectre d'indice 'index' dans la nappe + """ + return spectre(self.listFreq, self.listTable[index]) + + def filtreDoublons(self): + """ + Supprime bandes de fréquences constantes + """ + prevCpl = None + bCount = False + i=0 + # Recherche des doublons + lstBor = list() # Liste de listes de bornes + lst = list() + for cpl in self.__getListFreq() : + if not(prevCpl) : + prevCpl = cpl + continue + bTest = True + for v1, v2 in zip(cpl[1], prevCpl[1]) : + bTest &= (v1==v2) + if bTest and not bCount : # Début d'une suite de valeurs égales + bCount = True + lst.append(i) + elif not bTest and bCount : # Fin d'une suite de valeurs égales + bCount = False + lst.append(i) + lstBor.append(lst) + lst = list() # Nouvelle liste + + prevCpl = cpl + i += 1 + + # Suppression des doublons si plus de deux valeurs + for cpl in lstBor : + if (cpl[1]-cpl[0]) < 2 : continue + for i in range(cpl[1]-1, cpl[0], -1) : + del self.listFreq[i] + for j in range(0, len(self.listTable)) : + del self.listTable[j][i] + + + + + def __getListFreq(self) : + """ + Fonction privé qui parcours la matrice ligne par ligne + Retourne à chaque itération un couple frequence, liste de valeurs + """ + fr = 0.0 + + for i in range(0, self.getNbFreq()) : + fr = self.listFreq[i] + listVal = [] + for j in range(0, len(self.listTable)): + listVal.append(self.listTable[j][i]) + yield (fr, listVal) + + raise StopIteration + +class listSpectre : + """ + classe container d'une liste de spectre ne partageant pas la meme base de fréquence + cas des spectres à l'issue de la première passe de l'opération de filtrage d'enveloppe + """ + def __init__(self, *listSp) : + self.listSp = [] + for sp in listSp : + self.listSp = sp + + def append(self, spectre) : + """ Ajoute un spectre à la liste """ + self.listSp.append(spectre) + + def __staticBuildFromNappe(uneNappe) : + """ + Construit une liste de spectres (indépendants) à partir d'une nappe + """ + res = listSpectre() + for i in range(0, len(uneNappe.listAmor)) : + res.append(uneNappe.getSpectre(i)) + + return res + + buildFromNappe = staticmethod(__staticBuildFromNappe) #Définition en tant que méthode statique + + def testChevauchement(self) : + """ + Supprime les effets de chevauchement entre les spectres + """ + for i in range(0, len(self.listSp)-1) : + filter = filtreChevauchement(ref=self.listSp[i+1]) + self.listSp[i] = self.listSp[i].filtre(filter) + + def createBase(self, lspRef = None) : + """ + Crée une base de fréquence commune pour l'ensemble des spectres + En s'assurant que le l'on reste enveloppe des spectre de la liste lspRef + """ + lspRes = listSpectre([spectre() for sp in self.listSp]) # Liste résultante + + # Recherche des fréquences attribuées à 5 spectres, 4 spectres, ... classées dans un dictionnaire + dctOc = self.__sortByOccurence() + + iOcc = max(dctOc.keys()) + lst = dctOc[iOcc] # On comence par mettre les frequences communes à tous les spectres + lst.sort() + iOcc -= 1 + test = 0 + while True : + lspRes.__addFreqFromList(self, lst) + # On vérifie si on reste enveloppe du spectre initial + spTest = spectre() + lstComp = list() + for sp0, sp1 in zip(lspRes.listSp, self.listSp) : + filter = filtreChevauchement(ref=sp1) + spTest = sp0.filtre(filter) + # Crée une liste des fréquences ajoutées (s'il y en a...) + for fr in spTest.listFreq : + try : + idx = sp0.listFreq.index(fr) + except ValueError : # Valeur non trouvée dans le tableau + lstComp.append(fr) + + if len(lstComp) > 0 : # Il est nécessaire de compléter les spectres + # on prend de préférence les fréquences définies sur le plus de spectre possible + while True : + lstFreq = dctOc[iOcc] + prevLst = lst # On sauvegarde la liste précédente pour comparaison + lst = self.__buildList(lstComp, lstFreq) + if not(inclus(lst, prevLst)) : + iOcc -= 1 + else : + break + continue + else : + break # On peut sortir les spectres sont complets + + self.listSp = lspRes.listSp # Remplacement de la liste des spectres + + # On s'assure que le spectre reste enveloppe du spectre de référence rajoute des fréquences si nécessaire + # 1. filtre chevauchement + if lspRef : # Si une liste de spectre de référence a été définie, on vérifie le caractère enveloppe du résultat + listComp = list() + + for sp1, sp2 in zip(self.listSp, lspRef.listSp) : + filter = filtreChevauchement(ref=sp2) + spTest = sp1.filtre(filter) + test = inclus(spTest.listFreq, sp1.listFreq) + if test : listComp.append(test) + # 3. Complément éventuel de l'ensemble des spectres + if listComp : lspRes.__addFreqFromList(self, listComp) + + self.listSp = lspRes.listSp # Remplacement de la liste des spectres + + def filtre(self, filter): + """ + Applique un filtre à l'ensemble des spectres de la liste + """ + self.listSp = [sp.filtre(filter) for sp in self.listSp] + + + def __sortByOccurence(self) : + """ + Fonction qui trie les fréquences par leur occurence d'apparition dans la liste de spectre + """ + dct = {} + for sp in self.listSp : # Boucle sur tous les spectres + for fr in sp.listFreq : # Boucle sur toutes les fréquences de chaque spectre + try : + dct[fr] += 1 + except KeyError : + dct[fr] = 1 + + # "Inversion" du dictionnaire + dctOc = {} # Dictionnaire des occurences + for k in dct.keys() : + try : + dctOc[dct[k]].append(k) + except KeyError : + dctOc[dct[k]]=[k] + + + return dctOc + + def __addFreqFromList(self, lstSp, lstFreq) : + """ + Rajoute les fréquences contenues dans lstFreq aux spectres d'un listeSpectre + à partir des spectres fournis par le listeSpectre (lstSp) passé en paramètre + en procédant éventuellement à une interpolation linéaire + """ + # Suppression des doublons de la liste des fréquences + lstFreq = listToDict(lstFreq).keys() # lst est la liste des points qu'il faudrait ajouter pour rester enveloppe + lstFreq.sort() + + for i in range(0, len(self.listSp)) : + # Conversion des spectres en dictionnaire pour pouvoir les traiter + spDctSelf = self.listSp[i].buildMap() + spDctRef = lstSp.listSp[i].buildMap() + for fr in lstFreq : + # On cherche la valeur dans le spectre de référence + try : + vr = spDctRef[fr] + except KeyError : + ki = nearestKeys(fr, spDctRef) + vr = interpole(fr, ki[0], spDctRef[ki[0]], ki[1], spDctRef[ki[1]]) + + # On rajoute la valeur dans le spectre résultat + spDctSelf[fr] = vr + + # Conversion du dictionnaire en spectre réel + self.listSp[i] = spectre.sortSpectre(spDctSelf) + + def __buildList(self, lstComp, lstFreq) : + """ + Construit une liste de fréquences à ajouter à partir d'une liste de fréquences + à ajouter (listComp) et d'une liste de référence, (listFreq) + retourne une liste + """ + lst = list() + for fr in lstComp : + try : + idx = lstFreq.index(fr) + lst.append(fr) + except ValueError : # Fréquence non présente, recherche de la plus proche + couple = nearestKeys(fr, listToDict(lstFreq)) + if abs(couple[0]-fr) > abs(couple[1]-fr) : + lst.append(couple[1]) + else : + lst.append(couple[0]) + + lst = listToDict(lst).keys() # Suppression des doublons + lst.sort() + return lst + + +def lissage(nappe=nappe,fmin=0.2,fmax=35.5,elarg=0.1,tole_liss=0.25) : + resultat = listSpectre() # Le résultat sera contenu dans une liste de spectre + lspBrut = listSpectre.buildFromNappe(nappe) + # Passage en LogLog + lspBrut.filtre(filtreLog()) + for j in range(0,nappe.getNbSpectres()) : + # Spectre brut + sp = nappe.getSpectre(j) + # Limitation de la bande de fréquence + filter = filtreBandWidth(lower=fmin, upper=fmax) + sp = sp.filtre(filter) + # Expansion du spectre + filter = filtreExpand(coef=elarg) + sp = sp.filtre(filter) + # Passage en LogLin + filter = filtreLog(logOrd=False) + sp = sp.filtre(filter) + # éclatement du spectre en 3 sous-parties + tabSpectre = sp.cut(sp.seuil()) + # traitement individuel des sous parties + filter = filtreCrible(tolerance=2.*tole_liss) + tabSpectre[0] = tabSpectre[0].filtre(filter) + tabSpectre[2] = tabSpectre[2].filtre(filter) + filter.tolerance = tole_liss + tabSpectre[1] = tabSpectre[1].filtre(filter) + # Fusion des sous-spectres + sp = spectre.merge(tabSpectre) + + # Seconde passe de filtrage + sp = sp.filtre(filter) + + # On passe en log-log pour les tests de chevauchement + filter = filtreLog(logAbc=False) + sp = sp.filtre(filter) + # Ecriture dans la liste de spectres résultat + resultat.append(sp) # Ajoute la spectre lissé à la liste des spectres + + resultat.testChevauchement() # Test de chevauchement entre les spectre de la liste + resultat.createBase(lspBrut) # construction d'une base commune de fréquence + + # Passage en lin + resultat.filtre(filtreLin()) + + # Construction de la nappe résultat + nappeRes = nappe.buildFromListSpectre(resultat) + nappeRes.listAmor=nappe.listAmor + nappeRes.filtreDoublons() # Suppression des valeurs identiques accolées + + return nappeRes diff --git a/Aster/Cata/Utilitai/optimize.py b/Aster/Cata/Utilitai/optimize.py new file mode 100644 index 00000000..7735eb75 --- /dev/null +++ b/Aster/Cata/Utilitai/optimize.py @@ -0,0 +1,593 @@ +#@ MODIF optimize Utilitai DATE 31/10/2006 AUTEUR ASSIRE A.ASSIRE +# -*- coding: iso-8859-1 -*- +# CONFIGURATION MANAGEMENT OF EDF VERSION +# ====================================================================== +# COPYRIGHT (C) 1991 - 2006 EDF R&D WWW.CODE-ASTER.ORG +# THIS PROGRAM IS FREE SOFTWARE; YOU CAN REDISTRIBUTE IT AND/OR MODIFY +# IT UNDER THE TERMS OF THE GNU GENERAL PUBLIC LICENSE AS PUBLISHED BY +# THE FREE SOFTWARE FOUNDATION; EITHER VERSION 2 OF THE LICENSE, OR +# (AT YOUR OPTION) ANY LATER VERSION. +# +# THIS PROGRAM IS DISTRIBUTED IN THE HOPE THAT IT WILL BE USEFUL, BUT +# WITHOUT ANY WARRANTY; WITHOUT EVEN THE IMPLIED WARRANTY OF +# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. SEE THE GNU +# GENERAL PUBLIC LICENSE FOR MORE DETAILS. +# +# YOU SHOULD HAVE RECEIVED A COPY OF THE GNU GENERAL PUBLIC LICENSE +# ALONG WITH THIS PROGRAM; IF NOT, WRITE TO EDF R&D CODE_ASTER, +# 1 AVENUE DU GENERAL DE GAULLE, 92141 CLAMART CEDEX, FRANCE. +# ====================================================================== +# RESPONSABLE ASSIRE A.ASSIRE +# + +# ******NOTICE*************** +# optimize.py module by Travis E. Oliphant +# +# You may copy and use this module as you see fit with no +# guarantee implied provided you keep this notice in all copies. +# *****END NOTICE************ + +# A collection of optimization algorithms. Version 0.3.1 + +# Minimization routines +"""optimize.py + +A collection of general-purpose optimization routines using Numeric + +fmin --- Nelder-Mead Simplex algorithm (uses only function calls) +fminBFGS --- Quasi-Newton method (uses function and gradient) +fminNCG --- Line-search Newton Conjugate Gradient (uses function, gradient + and hessian (if it's provided)) + +""" +import Numeric +import MLab +Num = Numeric +max = MLab.max +min = MLab.min +abs = Num.absolute +__version__="0.3.1" + +def rosen(x): # The Rosenbrock function + return MLab.sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0) + +def rosen_der(x): + xm = x[1:-1] + xm_m1 = x[:-2] + xm_p1 = x[2:] + der = MLab.zeros(x.shape,x.typecode()) + der[1:-1] = 200*(xm-xm_m1**2) - 400*(xm_p1 - xm**2)*xm - 2*(1-xm) + der[0] = -400*x[0]*(x[1]-x[0]**2) - 2*(1-x[0]) + der[-1] = 200*(x[-1]-x[-2]**2) + return der + +def rosen3_hess_p(x,p): + assert(len(x)==3) + assert(len(p)==3) + hessp = Num.zeros((3,),x.typecode()) + hessp[0] = (2 + 800*x[0]**2 - 400*(-x[0]**2 + x[1])) * p[0] \ + - 400*x[0]*p[1] \ + + 0 + hessp[1] = - 400*x[0]*p[0] \ + + (202 + 800*x[1]**2 - 400*(-x[1]**2 + x[2]))*p[1] \ + - 400*x[1] * p[2] + hessp[2] = 0 \ + - 400*x[1] * p[1] \ + + 200 * p[2] + + return hessp + +def rosen3_hess(x): + assert(len(x)==3) + hessp = Num.zeros((3,3),x.typecode()) + hessp[0,:] = [2 + 800*x[0]**2 -400*(-x[0]**2 + x[1]), -400*x[0], 0] + hessp[1,:] = [-400*x[0], 202+800*x[1]**2 -400*(-x[1]**2 + x[2]), -400*x[1]] + hessp[2,:] = [0,-400*x[1], 200] + return hessp + + +def fmin(func, x0, args=(), xtol=1e-4, ftol=1e-4, maxiter=None, maxfun=None, fulloutput=0, printmessg=1): + """xopt,{fval,warnflag} = fmin(function, x0, args=(), xtol=1e-4, ftol=1e-4, + maxiter=200*len(x0), maxfun=200*len(x0), fulloutput=0, printmessg=0) + + Uses a Nelder-Mead Simplex algorithm to find the minimum of function + of one or more variables. + """ + x0 = Num.asarray(x0) + assert (len(x0.shape)==1) + N = len(x0) + if maxiter is None: + maxiter = N * 200 + if maxfun is None: + maxfun = N * 200 + + rho = 1; chi = 2; psi = 0.5; sigma = 0.5; + one2np1 = range(1,N+1) + + sim = Num.zeros((N+1,N),x0.typecode()) + fsim = Num.zeros((N+1,),'d') + sim[0] = x0 + fsim[0] = apply(func,(x0,)+args) + nonzdelt = 0.05 + zdelt = 0.00025 + for k in range(0,N): + y = Num.array(x0,copy=1) + if y[k] != 0: + y[k] = (1+nonzdelt)*y[k] + else: + y[k] = zdelt + + sim[k+1] = y + f = apply(func,(y,)+args) + fsim[k+1] = f + + ind = Num.argsort(fsim) + fsim = Num.take(fsim,ind) # sort so sim[0,:] has the lowest function value + sim = Num.take(sim,ind,0) + + iterations = 1 + funcalls = N+1 + + while (funcalls < maxfun and iterations < maxiter): + if (max(Num.ravel(abs(sim[1:]-sim[0]))) <= xtol \ + and max(abs(fsim[0]-fsim[1:])) <= ftol): + break + + xbar = Num.add.reduce(sim[:-1],0) / N + xr = (1+rho)*xbar - rho*sim[-1] + fxr = apply(func,(xr,)+args) + funcalls = funcalls + 1 + doshrink = 0 + + if fxr < fsim[0]: + xe = (1+rho*chi)*xbar - rho*chi*sim[-1] + fxe = apply(func,(xe,)+args) + funcalls = funcalls + 1 + + if fxe < fxr: + sim[-1] = xe + fsim[-1] = fxe + else: + sim[-1] = xr + fsim[-1] = fxr + else: # fsim[0] <= fxr + if fxr < fsim[-2]: + sim[-1] = xr + fsim[-1] = fxr + else: # fxr >= fsim[-2] + # Perform contraction + if fxr < fsim[-1]: + xc = (1+psi*rho)*xbar - psi*rho*sim[-1] + fxc = apply(func,(xc,)+args) + funcalls = funcalls + 1 + + if fxc <= fxr: + sim[-1] = xc + fsim[-1] = fxc + else: + doshrink=1 + else: + # Perform an inside contraction + xcc = (1-psi)*xbar + psi*sim[-1] + fxcc = apply(func,(xcc,)+args) + funcalls = funcalls + 1 + + if fxcc < fsim[-1]: + sim[-1] = xcc + fsim[-1] = fxcc + else: + doshrink = 1 + + if doshrink: + for j in one2np1: + sim[j] = sim[0] + sigma*(sim[j] - sim[0]) + fsim[j] = apply(func,(sim[j],)+args) + funcalls = funcalls + N + + ind = Num.argsort(fsim) + sim = Num.take(sim,ind,0) + fsim = Num.take(fsim,ind) + iterations = iterations + 1 + + x = sim[0] + fval = min(fsim) + warnflag = 0 + + if funcalls >= maxfun: + warnflag = 1 + if printmessg: + print "Warning: Maximum number of function evaluations has been exceeded." + elif iterations >= maxiter: + warnflag = 2 + if printmessg: + print "Warning: Maximum number of iterations has been exceeded" + else: + if printmessg: + print "Optimization terminated successfully." + print " Current function value: %f" % fval + print " Iterations: %d" % iterations + print " Function evaluations: %d" % funcalls + + if fulloutput: + return x, fval, warnflag + else: + return x + + +def zoom(a_lo, a_hi): + pass + + + +def line_search(f, fprime, xk, pk, gfk, args=(), c1=1e-4, c2=0.9, amax=50): + """alpha, fc, gc = line_search(f, xk, pk, gfk, + args=(), c1=1e-4, c2=0.9, amax=1) + + minimize the function f(xk+alpha pk) using the line search algorithm of + Wright and Nocedal in 'Numerical Optimization', 1999, pg. 59-60 + """ + + fc = 0 + gc = 0 + alpha0 = 1.0 + phi0 = apply(f,(xk,)+args) + phi_a0 = apply(f,(xk+alpha0*pk,)+args) + fc = fc + 2 + derphi0 = Num.dot(gfk,pk) + derphi_a0 = Num.dot(apply(fprime,(xk+alpha0*pk,)+args),pk) + gc = gc + 1 + + # check to see if alpha0 = 1 satisfies Strong Wolfe conditions. + if (phi_a0 <= phi0 + c1*alpha0*derphi0) \ + and (abs(derphi_a0) <= c2*abs(derphi0)): + return alpha0, fc, gc + + alpha0 = 0 + alpha1 = 1 + phi_a1 = phi_a0 + phi_a0 = phi0 + + i = 1 + while 1: + if (phi_a1 > phi0 + c1*alpha1*derphi0) or \ + ((phi_a1 >= phi_a0) and (i > 1)): + return zoom(alpha0, alpha1) + + derphi_a1 = Num.dot(apply(fprime,(xk+alpha1*pk,)+args),pk) + gc = gc + 1 + if (abs(derphi_a1) <= -c2*derphi0): + return alpha1 + + if (derphi_a1 >= 0): + return zoom(alpha1, alpha0) + + alpha2 = (amax-alpha1)*0.25 + alpha1 + i = i + 1 + alpha0 = alpha1 + alpha1 = alpha2 + phi_a0 = phi_a1 + phi_a1 = apply(f,(xk+alpha1*pk,)+args) + + + +def line_search_BFGS(f, xk, pk, gfk, args=(), c1=1e-4, alpha0=1): + """alpha, fc, gc = line_search(f, xk, pk, gfk, + args=(), c1=1e-4, alpha0=1) + + minimize over alpha, the function f(xk+alpha pk) using the interpolation + algorithm (Armiijo backtracking) as suggested by + Wright and Nocedal in 'Numerical Optimization', 1999, pg. 56-57 + """ + + fc = 0 + phi0 = apply(f,(xk,)+args) # compute f(xk) + phi_a0 = apply(f,(xk+alpha0*pk,)+args) # compute f + fc = fc + 2 + derphi0 = Num.dot(gfk,pk) + + if (phi_a0 <= phi0 + c1*alpha0*derphi0): + return alpha0, fc, 0 + + # Otherwise compute the minimizer of a quadratic interpolant: + + alpha1 = -(derphi0) * alpha0**2 / 2.0 / (phi_a0 - phi0 - derphi0 * alpha0) + phi_a1 = apply(f,(xk+alpha1*pk,)+args) + fc = fc + 1 + + if (phi_a1 <= phi0 + c1*alpha1*derphi0): + return alpha1, fc, 0 + + # Otherwise loop with cubic interpolation until we find an alpha which satifies + # the first Wolfe condition (since we are backtracking, we will assume that + # the value of alpha is not too small and satisfies the second condition. + + while 1: # we are assuming pk is a descent direction + factor = alpha0**2 * alpha1**2 * (alpha1-alpha0) + a = alpha0**2 * (phi_a1 - phi0 - derphi0*alpha1) - \ + alpha1**2 * (phi_a0 - phi0 - derphi0*alpha0) + a = a / factor + b = -alpha0**3 * (phi_a1 - phi0 - derphi0*alpha1) + \ + alpha1**3 * (phi_a0 - phi0 - derphi0*alpha0) + b = b / factor + + alpha2 = (-b + Num.sqrt(abs(b**2 - 3 * a * derphi0))) / (3.0*a) + phi_a2 = apply(f,(xk+alpha2*pk,)+args) + fc = fc + 1 + + if (phi_a2 <= phi0 + c1*alpha2*derphi0): + return alpha2, fc, 0 + + if (alpha1 - alpha2) > alpha1 / 2.0 or (1 - alpha2/alpha1) < 0.96: + alpha2 = alpha1 / 2.0 + + alpha0 = alpha1 + alpha1 = alpha2 + phi_a0 = phi_a1 + phi_a1 = phi_a2 + +epsilon = 1e-8 + +def approx_fprime(xk,f,*args): + f0 = apply(f,(xk,)+args) + grad = Num.zeros((len(xk),),'d') + ei = Num.zeros((len(xk),),'d') + for k in range(len(xk)): + ei[k] = 1.0 + grad[k] = (apply(f,(xk+epsilon*ei,)+args) - f0)/epsilon + ei[k] = 0.0 + return grad + +def approx_fhess_p(x0,p,fprime,*args): + f2 = apply(fprime,(x0+epsilon*p,)+args) + f1 = apply(fprime,(x0,)+args) + return (f2 - f1)/epsilon + + +def fminBFGS(f, x0, fprime=None, args=(), avegtol=1e-5, maxiter=None, fulloutput=0, printmessg=1): + """xopt = fminBFGS(f, x0, fprime=None, args=(), avegtol=1e-5, + maxiter=None, fulloutput=0, printmessg=1) + + Optimize the function, f, whose gradient is given by fprime using the + quasi-Newton method of Broyden, Fletcher, Goldfarb, and Shanno (BFGS) + See Wright, and Nocedal 'Numerical Optimization', 1999, pg. 198. + """ + + app_fprime = 0 + if fprime is None: + app_fprime = 1 + + x0 = Num.asarray(x0) + if maxiter is None: + maxiter = len(x0)*200 + func_calls = 0 + grad_calls = 0 + k = 0 + N = len(x0) + gtol = N*avegtol + I = MLab.eye(N) + Hk = I + + if app_fprime: + gfk = apply(approx_fprime,(x0,f)+args) + func_calls = func_calls + len(x0) + 1 + else: + gfk = apply(fprime,(x0,)+args) + grad_calls = grad_calls + 1 + xk = x0 + sk = [2*gtol] + while (Num.add.reduce(abs(gfk)) > gtol) and (k < maxiter): + pk = -Num.dot(Hk,gfk) + alpha_k, fc, gc = line_search_BFGS(f,xk,pk,gfk,args) + func_calls = func_calls + fc + xkp1 = xk + alpha_k * pk + sk = xkp1 - xk + xk = xkp1 + if app_fprime: + gfkp1 = apply(approx_fprime,(xkp1,f)+args) + func_calls = func_calls + gc + len(x0) + 1 + else: + gfkp1 = apply(fprime,(xkp1,)+args) + grad_calls = grad_calls + gc + 1 + + yk = gfkp1 - gfk + k = k + 1 + + rhok = 1 / Num.dot(yk,sk) + A1 = I - sk[:,Num.NewAxis] * yk[Num.NewAxis,:] * rhok + A2 = I - yk[:,Num.NewAxis] * sk[Num.NewAxis,:] * rhok + Hk = Num.dot(A1,Num.dot(Hk,A2)) + rhok * sk[:,Num.NewAxis] * sk[Num.NewAxis,:] + gfk = gfkp1 + + + if printmessg or fulloutput: + fval = apply(f,(xk,)+args) + if k >= maxiter: + warnflag = 1 + if printmessg: + print "Warning: Maximum number of iterations has been exceeded" + print " Current function value: %f" % fval + print " Iterations: %d" % k + print " Function evaluations: %d" % func_calls + print " Gradient evaluations: %d" % grad_calls + else: + warnflag = 0 + if printmessg: + print "Optimization terminated successfully." + print " Current function value: %f" % fval + print " Iterations: %d" % k + print " Function evaluations: %d" % func_calls + print " Gradient evaluations: %d" % grad_calls + + if fulloutput: + return xk, fval, func_calls, grad_calls, warnflag + else: + return xk + + +def fminNCG(f, x0, fprime, fhess_p=None, fhess=None, args=(), avextol=1e-5, maxiter=None, fulloutput=0, printmessg=1): + """xopt = fminNCG(f, x0, fprime, fhess_p=None, fhess=None, args=(), avextol=1e-5, + maxiter=None, fulloutput=0, printmessg=1) + + Optimize the function, f, whose gradient is given by fprime using the + Newton-CG method. fhess_p must compute the hessian times an arbitrary + vector. If it is not given, finite-differences on fprime are used to + compute it. See Wright, and Nocedal 'Numerical Optimization', 1999, + pg. 140. + """ + + x0 = Num.asarray(x0) + fcalls = 0 + gcalls = 0 + hcalls = 0 + approx_hessp = 0 + if fhess_p is None and fhess is None: # Define hessian product + approx_hessp = 1 + + xtol = len(x0)*avextol + update = [2*xtol] + xk = x0 + k = 0 + while (Num.add.reduce(abs(update)) > xtol) and (k < maxiter): + # Compute a search direction pk by applying the CG method to + # del2 f(xk) p = - grad f(xk) starting from 0. + b = -apply(fprime,(xk,)+args) + gcalls = gcalls + 1 + maggrad = Num.add.reduce(abs(b)) + eta = min([0.5,Num.sqrt(maggrad)]) + termcond = eta * maggrad + xsupi = 0 + ri = -b + psupi = -ri + i = 0 + dri0 = Num.dot(ri,ri) + + if fhess is not None: # you want to compute hessian once. + A = apply(fhess,(xk,)+args) + hcalls = hcalls + 1 + + while Num.add.reduce(abs(ri)) > termcond: + if fhess is None: + if approx_hessp: + Ap = apply(approx_fhess_p,(xk,psupi,fprime)+args) + gcalls = gcalls + 2 + else: + Ap = apply(fhess_p,(xk,psupi)+args) + hcalls = hcalls + 1 + else: + Ap = Num.dot(A,psupi) + # check curvature + curv = Num.dot(psupi,Ap) + if (curv <= 0): + if (i > 0): + break + else: + xsupi = xsupi + dri0/curv * psupi + break + alphai = dri0 / curv + xsupi = xsupi + alphai * psupi + ri = ri + alphai * Ap + dri1 = Num.dot(ri,ri) + betai = dri1 / dri0 + psupi = -ri + betai * psupi + i = i + 1 + dri0 = dri1 # update Num.dot(ri,ri) for next time. + + pk = xsupi # search direction is solution to system. + gfk = -b # gradient at xk + alphak, fc, gc = line_search_BFGS(f,xk,pk,gfk,args) + fcalls = fcalls + fc + gcalls = gcalls + gc + + update = alphak * pk + xk = xk + update + k = k + 1 + + if printmessg or fulloutput: + fval = apply(f,(xk,)+args) + if k >= maxiter: + warnflag = 1 + if printmessg: + print "Warning: Maximum number of iterations has been exceeded" + print " Current function value: %f" % fval + print " Iterations: %d" % k + print " Function evaluations: %d" % fcalls + print " Gradient evaluations: %d" % gcalls + print " Hessian evaluations: %d" % hcalls + else: + warnflag = 0 + if printmessg: + print "Optimization terminated successfully." + print " Current function value: %f" % fval + print " Iterations: %d" % k + print " Function evaluations: %d" % fcalls + print " Gradient evaluations: %d" % gcalls + print " Hessian evaluations: %d" % hcalls + + if fulloutput: + return xk, fval, fcalls, gcalls, hcalls, warnflag + else: + return xk + + + +if __name__ == "__main__": + import string + import time + + + times = [] + algor = [] + x0 = [0.8,1.2,0.7] + start = time.time() + x = fmin(rosen,x0) + print x + times.append(time.time() - start) + algor.append('Nelder-Mead Simplex\t') + + start = time.time() + x = fminBFGS(rosen, x0, fprime=rosen_der, maxiter=80) + print x + times.append(time.time() - start) + algor.append('BFGS Quasi-Newton\t') + + start = time.time() + x = fminBFGS(rosen, x0, avegtol=1e-4, maxiter=100) + print x + times.append(time.time() - start) + algor.append('BFGS without gradient\t') + + + start = time.time() + x = fminNCG(rosen, x0, rosen_der, fhess_p=rosen3_hess_p, maxiter=80) + print x + times.append(time.time() - start) + algor.append('Newton-CG with hessian product') + + + start = time.time() + x = fminNCG(rosen, x0, rosen_der, fhess=rosen3_hess, maxiter=80) + print x + times.append(time.time() - start) + algor.append('Newton-CG with full hessian') + + print "\nMinimizing the Rosenbrock function of order 3\n" + print " Algorithm \t\t\t Seconds" + print "===========\t\t\t =========" + for k in range(len(algor)): + print algor[k], "\t -- ", times[k] + + + + + + + + + + + + + + + + -- 2.39.2