Commit 5042d349 authored by Nicolas Fressengeas's avatar Nicolas Fressengeas
Browse files

Merge branch 'master' of https://github.com/pierrelouislee/Escapade into supelec

Conflicts:
	Escapade.sage
	Makefile
	Mesh.spyx
	vecteurs.so
parents 3569286d 9312177e
......@@ -33,14 +33,20 @@ load(os.path.join(os.path.expanduser(EscapadeInstallDir),'dd.sage'))
# La classe principale
class Escapade():
def __init__(self,system,mesh,unknowns,constants=[],observers=[],dirichlet=[]):
"""The constructor constructs the three members : system, unknowns and mesh.
The three arguments are respectively a list of equations, a list of functions (the unknowns) and the mesh, of the Mesh class"
It also sets to 0 the state variables that indicate if the computations have been done and the fils written.
The constants argument states which variables are supposed to be known.
The observers are a list of expressions of the variables and constants which will computed at each automaton step, allowing one to "observe" a given value which is not directly computed
dirichlet is a list of lists, each of one corresponds to an equation in the system (first arg). This list contains those of the variables which are assumed to be known on that given system (in other words, those variables for which you give a Dirichlet condition)
def __init__(self,system,mesh,unknowns,constants=[],observers=[],dirichlet=[],method=method):
"""
Le constructeur du fichier ESCAPADE
* définit :
- le système : une liste d'équations
- les inconnus : une liste de fonctions
- les mailles : instance de la classe « Mesh »
* met à 0 les variables d'états qui indiquent si les calculs ont été faits et si les fichiers .cell sont écrits
* pour ses arguments :
- les constantes : fixent les valeurs des variables supposées être connues
- les observateurs : une liste d'expressions (contenant des variables et des constantes) qui seront calculés à chaque pas de l'automate ce
qui permet « d'observer » une valeur qui n'est pas directement calculée.
- « dirichlet » : une liste de listes, chacune correspondant à une équation du système. La liste contient les variables connus sur ce
système donné (des variables qui ont des conditions de Dirichlet)
"""
self.systemlist=FiniteDifferenceSystemList(system)
self.unknowns=unknowns
......@@ -76,7 +82,7 @@ class Escapade():
# 0 : the false time method
# 1 : the Newton method
# ...
self.method=1
self.method=method=4
#Now the value of state variables
self.patterns_done=0
self.field_done=0
......@@ -85,8 +91,6 @@ class Escapade():
self.zerovector=vector(self.local_args)
#The differentiation vector
self.diff_vector=vector(apply(f,self.local_args) for f in unknowns)
#Choice of the method to compute the rules
self.rulechoice=1
#The number of pattern points computed so far : usefull to recover from a maxima crash
self.points_computed_so_far=-1
#Time step used in the false time method
......@@ -118,60 +122,7 @@ class Escapade():
return("The systems that are to be solved (==0 are implicit):\n"+self.systemlist.__repr__()+"\nThe sought functions:\n"+self.unknowns.__repr__()+cs+os+"\nThe mesh :\n"+self.mesh.__repr__())
def rules(self,pos):
"""Computes the update rules at a given point : it is a list of rules.
Constants and observers are not included in it. They are added later.
Should be computed by summing the squared error on all the points pointed by dependency.
"""
# The system on the particular point we are on
localsystem=(self.systemlist[self.mesh[tuple(pos)]])
# Whatever the method, if the localsystem is null, return the no move rule
if localsystem==[0]:
return list(self.diff_vector)
# if localsystem.squared_norm()==0:
# print("squarednorm simplification")
# return list(self.diff_vector)
# For the False Time method, it is really easy, with a tricky part if the null system has only one component
if self.method==0:
rule=[]
i=0
for equation in localsystem:
rule.append(self.diff_vector[i]+self.timestep*equation)
i+=1
return rule
# somme is the sum to be minimized
diffsomme=0
# Sum of squared error is done over all the points designated by self.dependency when centered on pos
# Dp not forget to translate the arguments of all unknowns and constants
for point in self.dependency.points(pos):
# If point is outside the mesh, ignore it
if not self.mesh.is_in(point):
continue
# What system do we want to solve here
system=self.systemlist[self.mesh[tuple(point)]]
translatedsystem=(system.translation(vecteurs.diffv(point,pos),self.unknowns+self.constants))
# print (cputime())
somme_item=FiniteDifference(translatedsystem.squared_norm())
tmpgrad=self.Dirichlet[self.mesh[tuple(point)]].pairwise_product(grad(somme_item,self.diff_vector))
# print (cputime())
if tmpgrad!=self.zerovector:
if self.method==1:
diffsomme_item=(hessian(somme_item,self.diff_vector)).solve_right(tmpgrad)
elif self.method==2:
diffsomme_item=(self.descent*tmpgrad)
elif self.method==3:
diffsomme_item=((hessian(somme_item,self.diff_vector)-self.lambdanewton*identity_matrix(len(self.diff_vector))).solve_right(tmpgrad))
elif self.method==4:
diffsomme_item=(diaghessian(somme_item,self.diff_vector)*(tmpgrad))
else :
print("method==0: False Time method\nmethod==1: Newton method\nmethod=2: steepest gradient\nmethod=3: Modified Newton\nmethod=4: Simplified Newton")
return(-1);
diffsomme+=self.Dirichlet[self.mesh[tuple(point)]].pairwise_product(diffsomme_item)
# print (cputime())
return(list(self.diff_vector-diffsomme))
def rules1(self,pos):
......@@ -184,7 +135,7 @@ class Escapade():
localsystem=(self.systemlist[self.mesh[tuple(pos)]])
# Whatever the method, if the localsystem is null, return the no move rule
if localsystem==[0]:
#print("nomove rule")
return list(self.diff_vector)
# if localsystem.squared_norm()==0:
# print("squarednorm simplification")
......@@ -223,13 +174,12 @@ class Escapade():
elif self.method==3:
diffsomme=((hessian(somme,self.diff_vector)-self.lambdanewton*identity_matrix(len(self.diff_vector))).solve_right(tmpgrad))
elif self.method==4:
diffsomme=(diaghessian(somme,self.diff_vector)*(tmpgrad))
else :
print("method==0: False Time method\nmethod==1: Newton method\nmethod=2: steepest gradient\nmethod=3: Modified Newton\nmethod=4: Simplified Newton")
return(-1);
#print("point :",point)
#print("self.mesh[tuple(point)]:",self.mesh[tuple(point)])
#print("self.Dirichlet[self.mesh[tuple(point)]]:",self.Dirichlet[self.mesh[tuple(point)]])
return(list(self.diff_vector-self.Dirichlet[self.mesh[tuple(point)]].pairwise_product(diffsomme)))
......@@ -245,6 +195,7 @@ class Escapade():
print("Method 3: Modified Newton. the Hessian is replaced by H-lambda I. lambda is the attribute self.lambdanewton. For low lambdas, this is Newton Method while for large lambda, it ressembles the steepest gradient")
print("Method 4: Simplified Newton: the Hessian is replaced by its diagonal. Easier to compute.")
def makefield(self):
"""Computes the field attribute.
It is mesh of the same dimensions as the mesh attribute.
......@@ -269,6 +220,7 @@ class Escapade():
self.dependency.save("dependency.escapade.dump");
self.field_done=1;
def loadfield(self):
"""Loads the result of makefield instead of computing them.
Needs the following files in the current directory :
......@@ -332,10 +284,7 @@ class Escapade():
currentpoint+=1
if currentpoint<=self.points_computed_so_far:
continue
if self.rulechoice==0:
rule=self.rules(point)
else:
rule=self.rules1(point)
# print("Done point",point,"#",currentpoint, "time :",walltime()-debut)
#Simplify if required
if self.simplifyQ==1:
......@@ -343,18 +292,17 @@ class Escapade():
#Pickle and compress the rule before listing for memory saving
rule=zlib.compress(dumps(rule))
#Skip, or NOT, the lengthy test to limit the number of rules
#BUG : only one point is set to the new rule, leaving rule numbers in field which do not mean anything
#First work around : do not limit the number of rules
#inlist=False
#for i in range(self.rulenumber):
# if rule==self.rulelist[i]:
# inlist=True
# self.ruledict[self.field[tuple(point)]]=i #<---BUG
# break;
#if not inlist:
self.rulelist.append(rule)
self.ruledict[self.field[tuple(point)]]=self.rulenumber
self.rulenumber+=1
inlist=False
for i in range(self.rulenumber):
if rule==self.rulelist[i]:
inlist=True
self.ruledict[self.field[tuple(point)]]=i
break;
if not inlist:
self.rulelist.append(rule)
self.ruledict[self.field[tuple(point)]]=self.rulenumber
self.rulenumber+=1
#end of list adding
self.points_computed_so_far+=1
f=open("rulelist.escapade.dump","w")
......@@ -412,24 +360,19 @@ class Escapade():
if self.field_done==0:
self.makefield()
print("Computing update rules")
if (self.workers>1) and (self.method!=0):
self.makepatternsparallel()
else:
self.makepatterns()
self.makepatterns()
# Now let us add the rules for the constants and the obervers
constantrules=list(apply(f,self.local_args) for f in self.constants)
newrulelist=[]
while True:
try:
rule=loads(zlib.decompress(self.rulelist.pop(0)))
rule.extend(constantrules)
rule.extend(self.observers)
newrulelist.append(zlib.compress(dumps(rule)))
except IndexError:
self.rulelist=newrulelist
newrulelist=[]
break
newrulelist=[]
while self.rulelist :
rule=loads(zlib.decompress(self.rulelist.pop(0)))
rule.extend(constantrules)
rule.extend(self.observers)
newrulelist.append(zlib.compress(dumps(rule)))
self.rulelist=newrulelist
newrulelist=[]
self.patterns_done=1
......@@ -546,126 +489,6 @@ class Escapade():
self.points_computed_so_far=0
def write_simpffiles(self):
"""Writes the files for the old version of simpf.
"""
if self.patterns_done==0:
#print("Computations not done: run self.makeCNN first")
self.makeCNN()
# Apply a series of rules for compliance with booz
Crulelist=[]
for rules in self.rulelist:
Crules=[]
for rule in rules:
Crules.append(FiniteDifference(loads(zlib.decompress(rule))).cform(self.max_args))
Crulelist.append(Crules)
#Patterns.simpf
stringpattern=Crulelist.__repr__()
#replace each unknown and constant by ?(i) where i is its position in the list
i=0
stringpattern=stringpattern.replace("simpfarobas","@")
stringpattern=stringpattern.replace("[","List(")
stringpattern=stringpattern.replace("]",")")
for f in (self.unknowns+self.constants):
stringpattern=stringpattern.replace(f.__repr__(),"?("+i.__repr__()+")")
i+=1
stringpattern=simpfpattern(stringpattern)
f=open("Patterns.simpf","w")
f.write(stringpattern)
f.close()
#Field.simpf
f=open("Field.simpf","w")
f.write("List(")
for i,val in numpy.ndenumerate(self.field):
f.write(str(val))
f.write(",")
f.write(")")
f.close()
#Layout.simpf
f=open("Layout.simpf","w")
f.write(str(len(self.unknowns+self.constants+self.observers))+"\n")
f.write(str(len(self.rulelist))+"\n")
l=len(self.mesh.shape)
if l==4:
simpfshape=[self.mesh.shape[0],self.mesh.shape[1],self.mesh.shape[2],self.mesh.shape[3]]
elif l==3:
simpfshape=[1,self.mesh.shape[0],self.mesh.shape[1],self.mesh.shape[2]]
elif l==2:
simpfshape=[1,1,self.mesh.shape[0],self.mesh.shape[1]]
elif l==1:
simpfshape=[1,1,1,self.mesh.shape[0]]
else:
print("For simpf, Mesh dimensions must be between 1 and 4")
for val in simpfshape:
f.write(str(val)+" ")
f.write("\n")
f.close()
#A few parameters for Params.simpf
f=open("Params.simpf","w")
f.write(str(self.timestep))
#f.write("\t1\t#Time step used in the false time method\n")
f.write("\t1\n")
f.write(str(self.descent))
#f.write(str("\t1\t#Gradient coefficient for the steepest gradient descent\n"))
f.write(str("\t1\n"))
for params in self.params:
for param in params:
f.write(str(param)+"\t")
f.write("\n")
f.close()
#A minimal Functions.simpf
f=open("Functions.simpf","w")
f.write("pow\tpow\n")
f.close()
def simpfinit(self,initlist):
"""The method simpfinit takes as arguments a list of ndarrays corresponding to the initial values for unknowns and constants.
This is output to Init.simpf file.
"""
f=open("Init.simpf","w")
f.write('List(')
l=len(initlist)
for i,val in numpy.ndenumerate(initlist[0]):
for j in range(l):
f.write(str(initlist[j][i]))
f.write(",")
f.write("\n")
f.write(")")
f.close()
def read_simpf_result(self,snapshot):
"""The method read results reads the results from the snapshot indicated file.
And outputs it in a list of arrays.
"""
# The results are all output, including constants and obersvers
n_results=len(self.unknowns+self.constants+self.observers)
# start by a bunch of arrays filled with zeros
result_list=list(numpy.zeros(self.mesh.shape) for i in range(n_results))
# Open the file and read all the data in
# This gulps a lot of memory: twice as needed
# Can be corrected by either reading numbers one by one
# or line by line
reader=csv.reader(open(snapshot),delimiter=" ")
#skip 2 lines
# The first line could be used to check the file dimensionality agains self.shape
_=reader.next()
_=reader.next()
data=reader.next()
#remove last empty string
data.pop()
#convert to float
data=(map(float,data))
#fill the result_list
index=0 # index in the data list
for i,val in numpy.ndenumerate(result_list[0]):
for j in range(n_results):
result_list[j][i]=data[index]
index+=1
return (result_list)
#TODO :
......
This diff is collapsed.
###############################################################################
# ESCAPADE
# Ergonomic Solver using Cellular Automata for PArtial Differential Equation
# Copyright (C) 2009 Nicolas Fressengeas <nicolas@fressengeas.net>
# Copyright (C) 2009-2015 Nicolas Fressengeas <nicolas@fressengeas.net>
# Distributed under the terms of the GNU General Public License (GPL),
# version 2 or any later version. The full text of the GPL is available at:
# http://www.gnu.org/licenses/
......@@ -9,15 +9,14 @@
#!/bin/bash
VERSION = 1.04
DIST_NAME = escabooz.sage
PYTHON = sage
SETUP = setup.py
SETUPOPTIONS = build_ext --inplace
DISTDIR = /home/fresseng/Documents/HTML/Escapade
RELEASEDIR= releases
#DISTDIR = /home/fresseng/Documents/HTML/Escapade
RELEASEDIR= ../releases
INSTALLDIR= ~/.escapade
SAGEDIR=~/.sage
......@@ -61,7 +60,8 @@ dist :
@cp $(FILES) $(DIST_NAME).$(VERSION)
@tar zcvf $(DIST_NAME).$(VERSION).tar.gz $(DIST_NAME).$(VERSION)
@rm -rf $(DIST_NAME).$(VERSION)
@cp $(DIST_NAME).$(VERSION).tar.gz $(DISTDIR)/
#@cp $(DIST_NAME).$(VERSION).tar.gz $(DISTDIR)/
@mv $(DIST_NAME).$(VERSION).tar.gz $(RELEASEDIR)/
......@@ -79,7 +79,8 @@ install :
@cp *.sage $(INSTALLDIR)/
@cp *.spyx $(INSTALLDIR)/
@cp *.so $(INSTALLDIR)/
@cat EscapadeSageInit.sage >> $(SAGEDIR)/init.sage
@cat EscapadeSageInit.sage > $(SAGEDIR)/init.sage
uninstall :
@if [ ! -d $(INSTALLDIR) ]; then \
......
......@@ -30,10 +30,13 @@ import vecteurs
class Mesh():
#It does not inherit.
#It will have a ndarry object for the initial mesh
#It will have a ndarray object for the initial mesh
#And another for the pattern mesh
# The constructor builds a matrix dor the object mesh out of the list argument
def __init__(self,l):
""" c'est une fonction qui initialise le maillage à "l" (ayant la même forme et ne contenant que des entiers positifs).
"""
if (l.dtype!=int):
print("For the sake of speed : dtype type must be int")
# #Integers only are allowed for performance issues
......@@ -41,16 +44,25 @@ class Mesh():
self.mesh=l
self.shape=l.shape
if l.min()<0:
print("All array elements should be positive integers")
print("All array elements should be positive integers")
def __repr__(self):
""" Ceci donne la représentation du maillage
"""
return (self.mesh).__repr__()
# The printer prints the mesh object
#The printer prints the mesh object
# __getitem__ mimmicks the getitem from ndarray
# but returns -1 if oustside the mesh
def __getitem__(self,pos):
""" x.__getitem__(y) <==> x[y] si y appartient au maillage
<==> -1 sinon
len() retourne la dimension de l'objet
"pos" est la position à laquelle on veut regarder la valeur de la grandeur physique étudiée (y selon précédemment)
"""
if len(pos)>len(self.shape):
print("dimension of argument ("+str(len(pos))+") must be less or equal to dimension of mesh ("+str(len(self.shape))+")")
n=len(self.shape)
......@@ -64,6 +76,12 @@ class Mesh():
def __setitem__(self,pos,x):
self.mesh[pos]=x
""" assigner à la position "pos" de la grille la valeur "x"
x correspond au numéro de l'équation à résoudre en ce point de la grille
"""
# The method local_pattern takes as arguments:
......@@ -71,9 +89,14 @@ class Mesh():
# a position in the mesh as a vector
# and returns
# a list containing the elements of mesh referenced by the patter (-1 if outside the mesh)
def local_pattern(self,pattern, pos):
""" Cette méthode prend en arguments un pattern (un modèle), et une position dans la grille (un vecteur) et retourne une matrice de mêmes dimensions (de classe Pattern) contenant des entiers dont chacun correspond à la valeur contenue dans la maille.
--> C'est une sorte de translation du pattern ou une "extraction" du voisinage du point qui nous intéresse (pos).
la méthode prédéfinie "append" et utilisée dans cette fonction permet de rajouter un élèment à la fin d'une liste (list.append(x))
"""
v=Pattern([])
for i in pattern:
v.append(self[vecteurs.sumv(i,pos)])
return (v)
......@@ -85,29 +108,47 @@ class Mesh():
# each integer stands for one possibly different update rule in the cellular network
# - a list of points, each of which corresponds to one different integer
# this list will be used to actually compute the update rules
def field(self,pattern):
""" Cette méthode prend un pattern en argument et renvoie
- Un tableau d'entiers ayant la même dimension que la grille dont chaque case contient le numéro de la règle de mise à jour
correspondante à cette position là dans le réseau cellulaire.
- Une liste de points dont chacun correspond à un entier différent. Cette liste sera utilisée pour calculer les règles de mise à jour.
"""
# Let's start with an empty pattern list
# Commençons par une liste de pattern vide
patt_list=[]
# Et une liste de points vide
# And an empty point list
point_list=[]
# The return array will have the same dimensions as self
# Le tableau à renvoyer doit avoir les mêmes dimensions que la grille et ne contiendra que des entiers
pattnum_array=numpy.zeros(self.shape,dtype=int)
# Happily, numpy provides an N dimensionnal iterator even if N is unknown
# Happily, numpy provides an N dimensional iterator even if N is unknown
# It yields the coordinates i and the value val (only i here is used)
# numpy possède un itérateur d'index multidimensionnel "numpy.ndenumerate"
# Il rend un itérateur rapportant les paires de cordonnées d'un tableau et des valeurs
for (i,val) in numpy.ndenumerate(self.mesh):
# the the Pattern pattern is mapped onto the mesh on point i
# the Pattern pattern is mapped onto the mesh on point i
# le pattern est reporté sur la grille et centré au point i.
# it yields a given local pattern
# il rapporte un pattern local donné
loc_patt=self.local_pattern(pattern,i)
# if it is not in the list, then add it
# and add the point to the point list
if patt_list.count(loc_patt)==0:
try:
pattnum_array[i]=patt_list.index(loc_patt)
except ValueError:
pattnum_array[i]=len(patt_list)
patt_list.append(loc_patt)
point_list.append(i)
# give the local point its position number in the list
pattnum_array[i]=patt_list.index(loc_patt)
# The return by value of the big array pattnum_array is not a good idea
# I do not know how to do better in python
# C could handle it
return([point_list,pattnum_array])
......@@ -115,7 +156,12 @@ class Mesh():
def is_in(self,point):
"""This method returns True if point is within self.
""" Cette méthode prend pour argument "point" de l'univers et vérifier s'il appartient à notre grille !
Ainsi, "point" est une liste de même dimension que "la dimension de la grille", c'est à dire si notre grille est en 3D "point" sera une liste de 3 points.
La méthode renvoie True si "point" appartient à la grille, sinon elle renvoie False.
Elle utilise la méthode prédéfinie "shape" qui renvoie les dimensions de la grille.
"""
n=len(self.shape)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment