L-BFGS: New input options

Now it is possible to specify WOLFE_C1_SL and WOLFE_C2_SL in the JSON input file.
Documentation is updated as well.
Smaller output improvements.

doc/latex/6_Parameter_Definition.tex

@@ -535,6 +535,8 @@ The corresponding functions are copied out of DENISE Black Edition, which is mai
 			"WOLFE_CONDITION" : "1",
 			"WOLFE_TRY_OLD_STEPLENGTH" : "1",
 			"WOLFE_NUM_TEST" : "5", 
+			"WOLFE_C1_SL" : "1e-4",
+			"WOLFE_C2_SL" : "0.9",
 			
 			"LBFGS_STEP_LENGTH" : "0",
 \end{verbatim}}}
@@ -561,7 +563,7 @@ The ensure a stable convergence when L-BFGS is used the step length has to satis
  	\bigtriangledown f(x+\alpha \cdot p)^{T} \cdot p \geq c_2 \cdot\bigtriangledown f(x)^{T} \cdot p
  \end{align}
  To verify this condition the full gradient (with all shots!) of the updated model have to be calculated, so the verification of this condition is very expensive in terms of computation time. In theory always a step length of one should be tried first, however to save time it is possible to chose the step length from the iteration before which satisfied the wolfe condition. This behavior can be controlled with WOLFE\_TRY\_OLD\_STEPLENGTH, if this is set to 1 always the old step length will be tried first. In practice the experience has shown that at the beginning of a FWI stage a couple of step lengths will be tested and then the algorithm converges to a step length which will be accepted by the wolfe condition until the misfit cannot be reduced further. With the variable WOLFE\_NUM\_TEST a maximum number of test calculations can be defined. If after WOLFE\_NUM\_TEST tests no step length could be found, the algorithm takes the step length of the tests which reduced the misfit most, even if this step length does not satisfy the wolfe condition. If no step length was found which reduces the misfit, the algorithm will switch to the next FWI stage in the workflow file or abort the inversion. 
- The variable $c_1$ will be $10^{-7}\cdot max(f(x)^{-1}$ and $c_2$ is set to 0.9 (they are hard coded). 
+ The variable $c_1$ will be $10^{-7}\cdot max(f(x)^{-1}$ and $c_2$ is set to 0.9. With WOLFE\_C1\_SL and WOLFE\_C2\_SL it is possible to manipulate $c_1$ and $c_2$ with the JSON input file. If the parameters are not set, they will be set to default values, so in most cases it should be the best to do not specify theses parameters in the JSON file.
  A step length search which is based on the wolfe condition can be used which the switch WOLFE\_CONDITION. Please note that this step length search is only available if GRAD\_METHOD==2. 
   
   

par/in_and_out/DENISE_INV_all_parameters.json

@@ -111,12 +111,15 @@
 			"WOLFE_CONDITION" : "1",
 			"WOLFE_NUM_TEST" : "5",
 			"WOLFE_TRY_OLD_STEPLENGTH" : "0",
+			"WOLFE_C1_SL" : "1e-4",
+			"WOLFE_C2_SL" : "0.9",
 
 "Approx. Hessian" : "comment",
 			"EPRECOND" : "3",
 			"EPSILON_WE" : "0.005",
 
 			"EPRECOND_ITER" : "0",
+			"EPRECOND_PER_SHOT" : "1",
 
 "Workflow" : "comment",
 			"USE_WORKFLOW" : "0",

src/LBFGS.c

@@ -245,7 +245,7 @@ void lbfgs_core(int iteration, int N_LBFGS, int NPAR_LBFGS,float ** s_LBFGS, flo
     /* give output so stdout */
     if(VERBOSE || 1) {
         for(w=1;w<=N_LBFGS;w++) {
-            fprintf(FP,"\n rho_LBFGS(%d)=%e\n",w,rho_LBFGS[w]);
+            fprintf(FP,"\n rho_LBFGS(%2d)=%e",w,rho_LBFGS[w]);
         }
         fprintf(FP,"\n h0=%e\n",h0);
     }

src/denise.c

@@ -111,7 +111,7 @@ int main(int argc, char **argv){
     int LBFGS=0,LBFGS_NPAR=3;
     int LBFGS_iter_start=1;
     float LBFGS_L2_temp;
-    float **s_LBFGS,**bfgsmod2,**bfgsmod3,**y_LBFGS, **bfgsgrad2, **bfgsgrad3, *rho_LBFGS,*bfgsscale2,*bfgsscale3;
+    float **s_LBFGS,**y_LBFGS, *rho_LBFGS;
     int l=0;
     int w=0;
     int m=0;
@@ -581,6 +581,10 @@ int main(int argc, char **argv){
         waveconv_rho_shot = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
         
         if(WOLFE_CONDITION){
+            
+            c1_SL=WOLFE_C1_SL;
+            c2_SL=WOLFE_C2_SL;
+            
             waveconv_old= matrix(-nd+1,NY+nd,-nd+1,NX+nd);
             if(!ACOUSTIC) waveconv_u_old= matrix(-nd+1,NY+nd,-nd+1,NX+nd);
             waveconv_rho_old= matrix(-nd+1,NY+nd,-nd+1,NX+nd);
@@ -1095,6 +1099,14 @@ int main(int argc, char **argv){
                     if(MYID==0){
                         if(iter==1){
                             FPL2=fopen(MISFIT_LOG_FILE,"w");
+                            /* Write header for misfit log file */
+                            if(GRAD_METHOD==1) {
+                                if (TIME_FILT==0){
+                                    fprintf(FPL2,"opteps_vp \t epst1[1] \t epst1[2] \t epst1[3] \t L2t[1] \t L2t[2] \t L2t[3] \t L2t[4] \n");}
+                                else{
+                                    fprintf(FPL2,"opteps_vp \t epst1[1] \t epst1[2] \t epst1[3] \t L2t[1] \t L2t[2] \t L2t[3] \t L2t[4] \t FC \n");
+                                }
+                            }
                         }
                         if(iter>1){FPL2=fopen(MISFIT_LOG_FILE,"a");}
                     }
@@ -1518,7 +1530,7 @@ int main(int argc, char **argv){
                         
                         fprintf(FP,"\n==================================================================================\n");
                         fprintf(FP,"\n MYID=%d * Starting simulation (forward model) for shot %d of %d. Iteration %d ** \n",MYID,ishot,nshots,iter);
-                        fprintf(FP,"\n==================================================================================\n\n");
+                        fprintf(FP,"\n==================================================================================\n");
                         
                         for (nt=1;nt<=8;nt++) srcpos1[nt][1]=srcpos[nt][ishot];
                         
@@ -1561,13 +1573,12 @@ int main(int argc, char **argv){
                         /*------------------------------------------------------------------------------*/
                         
                         if (((TIME_FILT==1) || (TIME_FILT==2)) && (QUELLART!=6) && (INV_STF==0)){
-                            fprintf(FP," Time Domain Filter applied: lowpass filter, corner frequency of %.2f Hz, order %d\n",FC,ORDER);
+                            fprintf(FP,"\n Time Domain Filter applied: Lowpass with corner frequency of %.2f Hz, order %d\n",FC,ORDER);
                             
                             /*time domain filtering of the source signal */
                             if(WAVETYPE==1||WAVETYPE==3) timedomain_filt(signals,FC,ORDER,nsrc_loc,ns,1);
                             if(WAVETYPE==2||WAVETYPE==3) timedomain_filt(signals_SH,FC,ORDER,nsrc_loc,ns,1);
                             
-                            
                             if(WAVETYPE==1||WAVETYPE==3){
                                 if ((QUELLTYPB==1)|| (QUELLTYPB==2)){
                                     /*time domain filtering of the observed data sectionvy_obs */
@@ -2122,7 +2133,7 @@ int main(int argc, char **argv){
                                 if(MYID==0){
                                     printf("\n==================================================================================\n");
                                     printf("\n MYID=%d *****  Starting simulation (backward model) for shot %d of %d  ********** \n",MYID,irec,nshots1);
-                                    printf("\n==================================================================================\n\n");
+                                    printf("\n==================================================================================\n");
                                 }
                                 
                                 
@@ -2972,6 +2983,8 @@ int main(int argc, char **argv){
                         use_wolfe_failsafe=1;
                         break;
                     } else {
+                        fprintf(FP,"\n After %d simulations no step length could be found which reduces the misfit and satisfy the wolfe condition.",countstep);
+                        fprintf(FP,"\n Will continue without model update.");
                         wolfe_SLS_failed=1;
                         break;
                     }
@@ -3028,6 +3041,21 @@ int main(int argc, char **argv){
             if(INVMAT!=0) break;
         }
         
+        if(wolfe_SLS_failed) {
+            
+            if (TIME_FILT==0){
+                if(MYID==0) fprintf(FPL2,"%e \t %d \t %d \t %f \t 0 \t %d \t %e \t %e \n",0.0,iter,wolfe_sum_FWI,0.0,countstep-1,L2_SL_old,L2_SL_old);}
+            else{
+                if(MYID==0) fprintf(FPL2,"%e \t %d \t %d \t %f \t 0 \t %d \t %e \t %e \t %f\n",0.0,iter,wolfe_sum_FWI,0.0,countstep-1,L2_SL_old,L2_SL_old,FC);
+            }
+            
+            /* No update is done here, however model fils are written to disk for easy post processing */
+            alpha_SL=0.0;
+            calc_mat_change_test(waveconv_up,waveconv_rho_up,waveconv_u_up,prhonp1,prho,ppinp1,ppi,punp1,pu,iter,1,INVMAT,alpha_SL,0,nfstart,Vs0,Vp0,Rho0,wavetype_start,s_LBFGS,N_LBFGS,LBFGS_NPAR,Vs_avg,Vp_avg,rho_avg,LBFGS_iter_start);
+            
+            alpha_SL_old=1;
+        }
+        
         /*-----------------------------------------------------*/
         /*       Wolfe condition: Model update                 */
         /*-----------------------------------------------------*/
@@ -3082,6 +3110,7 @@ int main(int argc, char **argv){
             gradient_optimization=1;
         }
         
+        
         opteps_vp=0.0;
         opteps_vs=0.0;
         opteps_rho=0.0;
@@ -3828,11 +3857,12 @@ int main(int argc, char **argv){
                     break;
                 }
                 
-                if(MYID==0) printf("Changing to corner frequency of %4.2f Hz \n",FC);
+                if(MYID==0) printf("\n Changing to corner frequency of %4.2f Hz \n",FC);
                 
                 /* Restart L-BFGS at next iteration */
                 LBFGS_iter_start=iter+1;
                 wolfe_SLS_failed=0;
+                alpha_SL_old=1;
             }
             
             /* ------------------------------------------------- */
@@ -3858,11 +3888,12 @@ int main(int argc, char **argv){
                 s=1;
                 min_iter_help=0;
                 min_iter_help=iter+MIN_ITER;
-                if(MYID==0) printf("Changing to corner frequency of %4.2f Hz \n",FC);
+                if(MYID==0) printf("\n Changing to corner frequency of %4.2f Hz \n",FC);
                 
                 /* Restart L-BFGS at next iteration */
                 LBFGS_iter_start=iter+1;
                 wolfe_SLS_failed=0;
+                alpha_SL_old=1;
             }
             
         }

src/exchange_par.c

@@ -110,6 +110,8 @@ void exchange_par(void){
     extern int WOLFE_CONDITION;
     extern int WOLFE_NUM_TEST;
     extern int WOLFE_TRY_OLD_STEPLENGTH;
+    extern float WOLFE_C1_SL;
+    extern float WOLFE_C2_SL;
     
     /* definition of local variables */
     
@@ -200,6 +202,9 @@ void exchange_par(void){
         
         fdum[62]=EPSILON_WE;
         
+        fdum[63]=WOLFE_C1_SL;
+        fdum[64]=WOLFE_C2_SL;
+        
         idum[1]  = NPROCX;
         idum[2]  = NPROCY;
         
@@ -465,6 +470,9 @@ void exchange_par(void){
     
     EPSILON_WE=fdum[62];
     
+    WOLFE_C1_SL=fdum[63];
+    WOLFE_C2_SL=fdum[64];
+    
     NPROCX = idum[1];
     NPROCY = idum[2];
     NPROC  = idum[4];

src/globvar.h

@@ -96,6 +96,8 @@ int LBFGS_STEP_LENGTH;
 int WOLFE_CONDITION=0;
 int WOLFE_NUM_TEST=5;
 int WOLFE_TRY_OLD_STEPLENGTH=1;
+float WOLFE_C1_SL;
+float WOLFE_C2_SL;
 
 /* trace kill variables */
 int TRKILL;

src/read_par_json.c

@@ -115,6 +115,9 @@ void read_par_json(FILE *fp, char *fileinp){
     extern int WOLFE_CONDITION;
     extern int WOLFE_NUM_TEST;
     extern int WOLFE_TRY_OLD_STEPLENGTH;
+    extern float WOLFE_C1_SL;
+    extern float WOLFE_C2_SL;
+    
     /* definition of local variables */
     
     int number_readobjects=0,fserr=0;
@@ -731,6 +734,14 @@ void read_par_json(FILE *fp, char *fileinp){
                                 WOLFE_TRY_OLD_STEPLENGTH=1;
                                 fprintf(fp,"Variable WOLFE_TRY_OLD_STEPLENGTH is set to default value %d.\n",WOLFE_TRY_OLD_STEPLENGTH);
                             }
+                            if (get_float_from_objectlist("WOLFE_C1_SL",number_readobjects,&WOLFE_C1_SL,varname_list, value_list)){
+                                WOLFE_C1_SL=1e-4;
+                                fprintf(fp,"Variable WOLFE_C1_SL is set to default value %f.\n",WOLFE_C1_SL);
+                            }
+                            if (get_float_from_objectlist("WOLFE_C2_SL",number_readobjects,&WOLFE_C2_SL,varname_list, value_list)){
+                                WOLFE_C2_SL=0.9;
+                                fprintf(fp,"Variable WOLFE_C2_SL is set to default value %f.\n",WOLFE_C2_SL);
+                            }
                         }
                     }
                 }

src/wolfe_condition.c

@@ -59,7 +59,7 @@ int check_wolfe(float steplength, float misfit_old, float misfit_new, float ** g
         grad_dot_p_new+=matrix_product(grad_new_vp, update_vp);
     }
     
-    fprintf(FP,"\n\n--------- Wolfe condition ----------- ");
+    fprintf(FP,"\n\n ------- Wolfe condition -------- ");
     
     /*---------------------------------*/
     /*       Check condition 1         */