Commit 64fe3230 authored by Alessandro Mirone's avatar Alessandro Mirone

Algorithmes python +

added LIPSHITZFACTOR
parent ab43cbce
......@@ -1185,6 +1185,7 @@ void CCspace_precalculations( CCspace * self, int ncpus ) {
self->gpu_context->CONICITY_MARGIN_DOWN = self->params.CONICITY_MARGIN_DOWN ;
self->gpu_context->ITER_RING_SIZE = self->params.ITER_RING_SIZE ;
self->gpu_context->LIPSCHITZFACTOR = self->params.LIPSCHITZFACTOR ;
self->gpu_context->gpu_mainInit( self->gpu_context,
......@@ -1666,7 +1667,7 @@ void CCspace_Sino_2_Slice( CCspace * self, float * dataA, int nslices, int ns
} else {
int dimslice= self->params.num_x, iter ;
float beta = self->params.BETA_TV;
float gamma = .9/ ( dimslice * self->params.num_projections );
float gamma = 1.0/self->params.LIPSCHITZFACTOR/ ( dimslice * self->params.num_projections );
// x = x0
// u_old = np.zeros((l, l))
......@@ -1678,7 +1679,7 @@ void CCspace_Sino_2_Slice( CCspace * self, float * dataA, int nslices, int ns
if(self->params.DENOISING_TYPE==5 || self->params.DENOISING_TYPE==6 || self->params.DENOISING_TYPE==7) {
// we dont need extra factors
} else if(self->params.ITERATIVE_CORRECTIONS_NOPREC ) {
gamma = .9/ ( dimslice * self->params.num_projections ); // /3 ;
gamma = 1.0/self->params.LIPSCHITZFACTOR/ ( dimslice * self->params.num_projections ); // /3 ;
/* { */
/* int i; */
/* for(i=0; i<100; i++) */
......@@ -1707,7 +1708,7 @@ void CCspace_Sino_2_Slice( CCspace * self, float * dataA, int nslices, int ns
1 ,1, NULL );
}
norma=1.1*norma ;
self->fbp_precalculated.prec_gamma = .9/norma;
self->fbp_precalculated.prec_gamma = 1.0/self->params.LIPSCHITZFACTOR/norma;
self->fbp_precalculated.prec_gamma_is_set=1;
}
......
......@@ -217,7 +217,7 @@ typedef struct {
float ITER_RING_HEIGTH;
float ITER_RING_BETA;
float ITER_RING_SIZE;
float LIPSCHITZFACTOR;
int JOSEPHNOCLIP;
......@@ -617,6 +617,7 @@ struct Gpu_Context_struct {
int * CONICITY_MARGIN_DOWN ;
float ITER_RING_SIZE ;
float LIPSCHITZFACTOR;
DFI_params dfi_params ;
......
......@@ -551,6 +551,7 @@ Cspace_new( PyObject *self, PyObject *args)
res->myCCspace->params.ITER_RING_HEIGTH = cpyutils_o2f (cpyutils_getattribute(res->OParameters ,"ITER_RING_HEIGTH" ) );
res->myCCspace->params.ITER_RING_SIZE = cpyutils_o2f (cpyutils_getattribute(res->OParameters ,"ITER_RING_SIZE" ) );
res->myCCspace->params.ITER_RING_BETA = cpyutils_o2f (cpyutils_getattribute(res->OParameters ,"ITER_RING_BETA" ) );
res->myCCspace->params.LIPSCHITZFACTOR = cpyutils_o2f (cpyutils_getattribute(res->OParameters ,"LIPSCHITZFACTOR" ) );
res->myCCspace->params.LOCAL_TOMO_ROI_RADIUS = cpyutils_o2f (cpyutils_getattribute(res->OParameters ,"LOCAL_TOMO_ROI_RADIUS" ) );
......
......@@ -381,7 +381,7 @@ void read_data_from_edf( char * filename,
FILE *edffile;
int res;
// printf(" apro %s \n", filename );
// printf(" apro %s filedatatype %d sizeofdatatype %d \n", filename, filedatatype , sizeofdatatype );
edffile = fopen( filename ,"r");
if(edffile==NULL) {
fprintf(stderr," error opening %s for reading \n", filename );
......
......@@ -2167,7 +2167,7 @@ void fb_dl(Gpu_Context * self, float *WORK , float * SLICE, int precondition,
}
if(failip) {
*Lipschitz = norma*1.2;
*Lipschitz = norma*self->LIPSCHITZFACTOR;
} else {
memcpy(SLICE, image, dim0*dim1*doppio*sizeof(float));
}
......
......@@ -1192,6 +1192,7 @@ class Parameters:
ITER_RING_HEIGTH=0.0
ITER_RING_SIZE=3.0
ITER_RING_BETA=0.0
LIPSCHITZFACTOR=1.1
ITERATIVE_CORRECTIONS_NOPREC = 1
......
......@@ -104,7 +104,7 @@ class Problem:
nit = 0
while res > 0.001:
x = tensordot(B,x,[[1],[0]])
i = argmax(x)
i = argmax( abs(x) )
x = x/float(x[i])
vp_tmp =vp
vp = tensordot(B,x,[[1],[0]])[i]
......@@ -189,6 +189,9 @@ class Problem:
px = 0
dpx = zeros(n)
ddpx = zeros(n)
## il n' y a pas besoin de faire une boucle.
## On peut utiliser numpy : numpy.abs, numpy.sign, numpy.log
## qui travaillent sur des vecteur
for i in arange(n):
absx = abs(x[i])
px += absx-self.s*log(1+absx/self.s)
......@@ -207,11 +210,19 @@ class Problem:
f : evaluation de la fonction en alpha
df : evaluation de la derivee en alpha
"""
## SideFunc est appelle plusiere fois a l' interieure d' une boucle
## ou x ne change pas mais seulement a
## On calcule a chaque fois Ap, Axmb, normAp , inner(Ap,Axmb)
## qui ne varient pas
Ap = tensordot(self.A,p,[[0],[0]])
normAp = (Ap*Ap).sum()
Axmb = tensordot(self.A,x,[[0],[0]])-self.b
phi = self.eval_phi(x+a*p)
h = 2*inner(Ap,Axmb)+2*a*normAp+self.rho*inner(p,phi[1])
dh = 2*normAp+self.rho*inner(p*p,phi[2])
return(h,dh)
......@@ -233,6 +244,12 @@ class Problem:
alpha : le pas recherche
"""
# localisation de la solution par dichotomie (a1 et a2 en sortie)
## on peu garder alpham a zero et alphap positif :
## c' est parce que on a la direction descente
## et si on se deplace le long de celle ci avec alpha>0
## on va finir par monter
alpham = alpha-1
alphap = alpha+1
hm = self.SideFunc(p,x,alpham)[0]
......@@ -245,6 +262,8 @@ class Problem:
elif hp == 0:
return(alpha+1)
while (hm*h > 0 and hp*h > 0):
## ici on peut utiliser un pas
## d' abord petit et ensuite le doubler a chaque fois
alpham -= 1
alphap += 1
hm = self.SideFunc(p,x,alpham)[0]
......@@ -340,6 +359,9 @@ class Problem:
Axmb = tensordot(self.A,x,[[0],[0]])-self.b
term1 = -inner(Ap,Axmb)/normAp
alpha = term1 # Linesearch initial value
## passe a bisection aussi Ap normAp Axmb , ils dependent seulement de x
alpha = self.Newton_Bisection(100,0.0001,p,x,alpha) # Line search par mixte Newton - Dichotomie
x = x+alpha*p
# Modif du 27 mars : calcul "exact" de beta au lieu de la formule de PR
......@@ -347,9 +369,14 @@ class Problem:
gradQ_tmp = gradQ
gradf_tmp = gradf
gradQ = self.gradient(x)
L = self.eval_phi(x)
gradL = L[1]
hessL = L[2]
## Ajouter ici l'option
## de tester l' original PR
##
gradf = gradQ + self.rho*gradL
var_gradQ = gradQ-gradQ_tmp
term = var_gradQ+self.rho*alpha*p*hessL
......
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