Context navigation

Back to Ticket #646

Ticket #646: ropp_1dvar_levmarq_tr2.2.f90

File ropp_1dvar_levmarq_tr2.2.f90, 26.5 KB (added by marquardt, 5 years ago)
ropp_1dvar_levmarq_tr2.f90 v2

Line
1	! $Id: ropp_1dvar_levmarq.f90 4452 2015-01-29 14:42:02Z idculv $
2
3	!***s 1DVar/ropp_1dvar_levmarq
4	!
5	! NAME
6	! ropp_1dvar_levmarq - Solve the 1DVar for background data using the
7	! Levenberg-Marquardt minimiser
8	!
9	! SYNOPSIS
10	! CALL ropp_1dvar_levmarq(obs, bg, state, config, diag)
11	!
12	!
13	! DESCRIPTION
14	! This subroutine evaluates a quadratic cost function for a
15	! variational data assimilation procedure.
16	!
17	! More specifically, this routine calculates a cost function
18	!
19	! 1 / \| -1 \| \
20	! J = - < y - H(x) \| O \| y - H(x) > +
21	! 2 \ \| \| /
22	!
23	! 1 / \| -1 \| \
24	! - < x - x \| B \| x - x >
25	! 2 \ b \| \| b/
26	!
27	! where the background state x_b is given by the state vector
28	! state.
29	!
30	! A solution for x is obtained by minimising J using the Levenberg-Marquardt
31	! minimisation method.
32	!
33	! INPUTS
34	! obs Observation data structure.
35	! bg Background data structure.
36	! state State vector structure.
37	! config Configuration structure.
38	! diag Diagnostics structure.
39	!
40	! OUTPUT
41	! state
42	! diag
43	!
44	! REFERENCES
45	! Madsen, K., Nielsen, H.B. & Tingleff, O., Methods for Non-Linear Least
46	! Squares Problems, Lyngby, Denmark, 2004. Available at:
47	! http://www2.imm.dtu.dk/pubdb/views/publication_details.php?id=3215
48	! Nielsen, H.B., Damping Parameter in Marquardtâs Method, Lyngby, Denmark,
49	! 1999. Available at:
50	! http://www2.imm.dtu.dk/pubdb/views/publication_details.php?id=648
51	!
52	! AUTHOR
53	! Met Office, Exeter, UK.
54	! Any comments on this software should be given via the ROM SAF
55	! Helpdesk at http://www.romsaf.org
56	!
57	! COPYRIGHT
58	! (c) EUMETSAT. All rights reserved.
59	! For further details please refer to the file COPYRIGHT
60	! which you should have received as part of this distribution.
61	!
62	!****
63
64	!-------------------------------------------------------------------------------
65	! 1. Bending angle
66	!-------------------------------------------------------------------------------
67
68	SUBROUTINE ropp_1dvar_levmarq_bangle(obs, bg, state, config, diag)
69
70	! 1.1 Declarations
71	! ----------------
72
73	USE typesizes, ONLY: wp => EightByteReal
74	USE ropp_utils
75	USE ropp_fm
76	USE ropp_1dvar, not_this => ropp_1dvar_levmarq_bangle
77	USE matrix
78
79	IMPLICIT NONE
80
81	TYPE(Obs1dBangle), INTENT(inout) :: obs ! Observation data
82	TYPE(State1dFM), INTENT(inout) :: bg ! Background data
83	TYPE(State1dFM), INTENT(inout) :: state ! State vector
84	TYPE(VarConfig), INTENT(in) :: config ! Configuration options
85	TYPE(VarDiag), INTENT(inout) :: diag ! Diagnostic output
86
87	REAL(wp) :: J ! Cost function value
88	TYPE(State1dFM) :: x ! Control vector
89	TYPE(State1dFM) :: x_old ! Control for LM parameter testing
90	TYPE(Obs1dBangle) :: y ! Forward model obs
91	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: K ! K-matrix
92	REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_x ! Change of state
93	REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_y ! Change of observation
94	REAL(wp), DIMENSION(:), ALLOCATABLE :: x_incr ! Increment of the state vector
95	REAL(wp), DIMENSION(:), ALLOCATABLE :: dJ_dx ! Cost function gradient
96	REAL(wp), DIMENSION(:), ALLOCATABLE :: diag_d2J ! Diagonal d2J/dx2
97	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: d2J_dx2 ! 2nd derivative cost fn
98	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: KO ! K O^-1
99	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: Bm1 ! B^-1
100	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: Om1 ! O^-1
101	REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_J ! Change of cost fn
102	REAL(wp), DIMENSION(:), ALLOCATABLE :: state_last ! Previous state vector
103	REAL(wp), DIMENSION(:), ALLOCATABLE :: state_sigma ! State std deviation
104	REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_state ! Change of state vector
105	REAL(wp) :: J_last ! Previous cost function
106	REAL(wp) :: J_min ! Minimum cost function
107	INTEGER :: i_pointer ! Value index
108	INTEGER :: i ! Counter
109	INTEGER :: n_iter ! Number of iterations
110	INTEGER :: n_iter_max ! Maximum number of iterations
111	INTEGER :: n_grad ! Number of gradient updates
112	INTEGER :: nobs ! Number of observations
113	INTEGER :: nstate ! No. of state elements
114	REAL(wp) :: mu ! Levenberg-Marquardt factor
115	REAL(wp) :: mu_max ! Maximum Levenberg-Marquardt factor
116	REAL(wp) :: rho ! Trust region radius
117	REAL(wp) :: beta ! Increment factor for mu
118	REAL(wp) :: gamma ! Decrement factor for mu
119	REAL(wp) :: nu ! Variation factor for mu
120
121	CHARACTER(len = 4) :: it_str
122	CHARACTER(len = 4) :: gr_str
123	CHARACTER(len = 15) :: ch_str, co_str
124	CHARACTER(len = 256) :: routine
125
126	! 1.2 Message handling
127	! --------------------
128
129	CALL message_get_routine(routine)
130	CALL message_set_routine('ropp_1dvar_levmarq')
131
132	! 1.3 Initialise variables
133	! ------------------------
134
135	i_pointer = 0
136	n_grad = 0
137	n_iter = 0
138	n_iter_max = 50
139
140	mu = 1.0E-2_wp
141	mu_max = 1.0E10_wp
142	beta = 2.0_wp
143	gamma = 3.0_wp
144	nu = beta
145
146	J = 0.0_wp
147	J_last = 1.0E30_wp
148	J_min = J_last
149
150	nstate = SIZE(bg%state)
151	nobs = SIZE(obs%bangle)
152
153	ALLOCATE(K(nobs, nstate))
154	ALLOCATE(delta_x(nstate))
155	ALLOCATE(delta_y(nobs))
156	ALLOCATE(x_incr(nstate))
157	ALLOCATE(dJ_dx(nstate))
158	ALLOCATE(diag_d2J(nstate))
159	ALLOCATE(d2J_dx2(nstate,nstate))
160	ALLOCATE(KO(nstate,nobs))
161	ALLOCATE(Bm1(nstate,nstate))
162	ALLOCATE(Om1(nobs,nobs))
163
164	IF (ALLOCATED(delta_J)) DEALLOCATE(delta_J)
165	ALLOCATE(delta_J(config%conv_check_n_previous))
166	delta_J(:) = 0.0_wp
167
168	IF (ALLOCATED(state_last)) DEALLOCATE(state_last)
169	ALLOCATE(state_last(SIZE(bg%state)))
170	state_last(:) = 0.0_wp
171
172	IF (ALLOCATED(state_sigma)) DEALLOCATE(state_sigma)
173	ALLOCATE(state_sigma(SIZE(bg%state)))
174	DO i = 1, SIZE(state_sigma)
175	state_sigma(i) = SQRT(bg%cov%d(i + i*(i-1)/2)) ! Direct read from matrix_pp
176	ENDDO
177
178	IF (ALLOCATED(delta_state)) DEALLOCATE(delta_state)
179	ALLOCATE(delta_state(config%conv_check_n_previous))
180	delta_state(:) = 0.0_wp
181
182	! 1.4 Inverse error covariances
183	! -----------------------------
184
185	Bm1 = matrix_invert(bg%cov)
186	Om1 = matrix_invert(obs%cov)
187
188	! 1.5 First guess and pseudo-observations
189	! ---------------------------------------
190
191	x = bg ; y = obs
192
193	! 1.6 Initial cost function
194	! -------------------------
195
196	IF (ASSOCIATED(x%ak)) THEN
197	CALL ropp_fm_state2state_ecmwf(x)
198	ELSE
199	CALL ropp_fm_state2state_meto(x)
200	END IF
201
202	CALL ropp_fm_bangle_1d(x, y)
203
204	delta_x = (x%state - bg%state)
205	delta_y = (y%bangle - obs%bangle) * y%weights
206
207	J = 0.5_wp * DOT_PRODUCT(delta_y, MATMUL(Om1, delta_y)) &
208	+ 0.5_wp * DOT_PRODUCT(delta_x, MATMUL(Bm1, delta_x))
209
210	diag%J_init = J
211	J_last = J
212	J_min = J
213
214	IF (config%minropp%impres == 0) THEN
215	WRITE(it_str, '(i4)') n_iter
216	WRITE(ch_str, '(g15.5)') J
217	ch_str = ADJUSTL(ch_str)
218	co_str = ' - '
219
220	CALL message(msg_cont, &
221	' n_iter = ' // it_str // ' J = ' // ch_str // &
222	' max(relative change in state) = ' // co_str)
223	END IF
224
225	! 1.7 Main minimisation iteration loop
226	! ------------------------------------
227
228	DO WHILE(n_iter < n_iter_max)
229
230	! 1.7.1 Bookkeeping
231
232	n_grad = n_grad + 1
233
234	! 1.7.2 Compute gradient and Hessian of cost function
235
236	CALL ropp_fm_bangle_1d_grad(x, y, K)
237
238	!TODO: (a) Doe we need this? (b) Does it change the cost function?
239	WHERE(ABS(delta_y) > 50.0_wp)
240	delta_y = 0.0_wp
241	END WHERE
242
243	KO = MATMUL(TRANSPOSE(K), Om1)
244	dJ_dx = MATMUL(KO, delta_y) + MATMUL(Bm1, delta_x)
245	d2J_dx2 = Bm1 + MATMUL(KO, K)
246
247	! 1.7.3 Levenberg-Marquardt adjustment of Hessian diagonal
248
249	DO i = 1, nstate
250	d2J_dx2(i,i) = d2J_dx2(i,i) * (1.0_wp + mu)
251	END DO
252
253	! 1.7.4 Solve matrix equation d2J_dx2 . dx = -dJ_dx and update state
254
255	x_incr = matrix_solve(d2J_dx2, - dJ_dx)
256
257	! 1.7.5 Update test state vector
258
259	x_old = x
260
261	IF (x%use_logq) THEN
262	x%state = x%state + SIGN(MIN(ABS(x_incr), ABS(x%state/2.0_wp)), x_incr)
263	ELSE
264	x%state = x%state + x_incr
265	END IF
266
267	! 1.7.6 Special handling of ionospheric state vector elements
268
269	IF (bg%direct_ion) THEN
270
271	i = nstate - 2
272	IF (x%state(i) < ropp_ZERO) THEN
273	CALL message(msg_warn, "Levenberg-Marquardt solver returns " // &
274	"Ne_max < 0 ... suggest examining final value. \n")
275	END IF
276
277	i = nstate - 1
278	IF (x%state(i) < 0.01_wp*bg%state(i)) THEN
279	CALL message(msg_warn, "Levenberg-Marquardt solver returns " // &
280	"H_peak < 1% of background ... resetting to background value. \n")
281	x%state(i) = bg%state(i)
282	END IF
283
284	i = nstate
285	IF (x%state(i) < 0.01_wp*bg%state(i)) THEN
286	CALL message(msg_warn, "Levenberg-Marquardt solver returns " // &
287	"H_width < 1% of background ... resetting to background value. \n")
288	x%state(i) = bg%state(i)
289	END IF
290
291	END IF
292
293	! 1.7.3 Compute cost function
294
295	IF (ASSOCIATED(x%ak)) THEN
296	CALL ropp_fm_state2state_ecmwf(x)
297	ELSE
298	CALL ropp_fm_state2state_meto(x)
299	END IF
300
301	CALL ropp_fm_bangle_1d(x, y)
302
303	delta_x = (x%state - bg%state)
304	delta_y = (y%bangle - obs%bangle) * y%weights
305
306	J = 0.5_wp * DOT_PRODUCT(delta_y, MATMUL(Om1, delta_y)) &
307	+ 0.5_wp * DOT_PRODUCT(delta_x, MATMUL(Bm1, delta_x))
308
309	! 1.8 Levenberg-Marquardt update and convergence check
310	! ----------------------------------------------------
311
312	! Trust region / gain parameter
313
314	rho = 2.0_wp * (J_last - J) / DOT_PRODUCT(x_incr, mu * x_incr - dJ_dx)
315
316	IF (rho > 0.0_wp) THEN ! Check convergence
317
318	! 1.8.4 Update mu
319
320	mu = mu * max(1.0_wp/gamma, 1.0_wp - (beta - 1.0_wp)(2.0_wprho - 1.0_wp)**3)
321	nu = beta
322
323	! 1.8.5 Keep track of current state
324
325	n_iter = n_iter + 1
326
327	IF (config % conv_check_apply) THEN
328	i_pointer = MOD(i_pointer + 1, config % conv_check_n_previous)
329	IF (i_pointer == 0) i_pointer = config % conv_check_n_previous
330
331	delta_J(i_pointer) = J_last - J
332	delta_state(i_pointer) = MAXVAL(ABS(state_last - x%state)/state_sigma)
333
334	state_last = x%state
335	J_last = J
336	END IF
337
338	IF (J < J_min) THEN
339	J_min = J
340	END IF
341
342	! 1.8.6 Check for convergence
343
344	IF (config % conv_check_apply) THEN
345
346	IF (config%minropp%impres == 0) THEN
347	WRITE(it_str, '(i4)') n_iter
348	WRITE(ch_str, '(g15.5)') J
349	ch_str = ADJUSTL(ch_str)
350	IF (n_iter > 0) THEN
351	WRITE(co_str, '(g15.5)') delta_state(i_pointer)
352	co_str = ADJUSTL(co_str)
353	ELSE
354	co_str = ' - '
355	END IF
356
357	CALL message(msg_cont, &
358	' n_iter = ' // it_str // ' J = ' // ch_str // &
359	' max(relative change in state) = ' // co_str)
360	END IF
361
362	IF (MAXVAL(delta_state) < config%conv_check_max_delta_state &
363	.AND. n_iter > config % conv_check_n_previous) THEN
364	WRITE(ch_str, '(g15.5)') config%conv_check_max_delta_state
365	ch_str = ADJUSTL(ch_str)
366	WRITE(it_str, '(i2)') config%conv_check_n_previous
367	it_str = ADJUSTL(it_str)
368	CALL message(msg_cont, '')
369	CALL message(msg_info, &
370	'Convergence assumed to be achieved as the state vector did ' // &
371	'not change by more\n ' // 'than ' // TRIM(ch_str) // ' ' // &
372	'relative to the assumed background errors for the last ' // &
373	TRIM(it_str) // ' iterations.\n')
374	EXIT
375	ELSE IF (MAXVAL(ABS(delta_J)) < config%conv_check_max_delta_J &
376	.AND. n_iter > config % conv_check_n_previous) THEN
377	WRITE(ch_str, '(g15.5)') config%conv_check_max_delta_J
378	ch_str = ADJUSTL(ch_str)
379	WRITE(it_str, '(i2)') config%conv_check_n_previous
380	it_str = ADJUSTL(it_str)
381	CALL message(msg_cont, '')
382	CALL message(msg_info, &
383	'Convergence assumed to be achieved as the cost function did ' // &
384	'not change by more\n ' // 'than ' // TRIM(ch_str) // &
385	' for the last ' // TRIM(it_str) // ' iterations.\n')
386	EXIT
387	ELSE IF (ABS(mu) > mu_max) THEN
388	WRITE(ch_str, '(g15.5)') mu
389	ch_str = ADJUSTL(ch_str)
390	CALL message(msg_cont, '')
391	CALL message(msg_info, &
392	'Levenberg-Marquardt parameter grew beyond ' // TRIM(ch_str) // &
393	'; convergence will not get much better any more.\n ' // &
394	'Check convergence parameters; conv_check_max_delta_J might ' // &
395	'be too demanding.\n')
396	EXIT
397	END IF
398
399	END IF ! Convergence check
400
401	ELSE ! rho < 0; retry with old state and updated mu
402
403	x = x_old
404	mu = mu * nu ; nu = 2.0_wp * nu
405
406	!TODO: Add a configuration option for mu update outputs
407	IF (config%minropp%impres == 0) THEN
408	WRITE(it_str, '(i4)') n_iter
409	WRITE(ch_str, '(g15.5)') J
410	ch_str = ADJUSTL(ch_str)
411	WRITE(co_str, '(g15.5)') mu
412	co_str = ADJUSTL(co_str)
413
414	CALL message(msg_cont, &
415	' n_iter = ' // it_str // ' J = ' // ch_str // &
416	' mu -> ' // co_str)
417	END IF
418
419	END IF ! Levenberg-Marquardt update
420	END DO ! end main iteration loop
421
422	IF (n_iter < n_iter_max) THEN
423	WRITE(it_str, '(i4)') n_iter ; it_str = ADJUSTL(it_str)
424	WRITE(gr_str, '(i4)') n_grad ; gr_str = ADJUSTL(gr_str)
425	CALL message(msg_info, 'Finished after ' // TRIM(it_str) // ' iterations (' // &
426	TRIM(gr_str) // ' forward model / gradient evaluations).\n')
427	ELSE
428	WRITE(it_str, '(i4)') n_iter ; it_str = ADJUSTL(it_str)
429	CALL message(msg_warn, &
430	'Iteration ended after ' // TRIM(it_str) // &
431	' iterations without convergence achieved.\n')
432	END IF
433
434	! 1.9 Copy solution back to state
435	! -------------------------------
436
437	state = x
438
439	! 1.10 Diagnostic data
440	! --------------------
441
442	diag%J = J_min
443
444	IF (COUNT(obs%weights > 0.0_wp) > 0) THEN
445	diag%J_scaled = 2.0_wp * J_min / REAL(COUNT(obs%weights > 0.0_wp), wp)
446	ENDIF
447
448	diag%n_iter = n_iter
449
450	ALLOCATE (diag%J_bgr(SIZE(state%state)))
451	delta_x = state%state - bg%state
452	diag%J_bgr = 0.5_wp * delta_x * matrix_solve(bg%cov, delta_x)
453
454	! 1.11 Clean up
455	! -------------
456
457	DEALLOCATE(K)
458	DEALLOCATE(delta_x)
459	DEALLOCATE(delta_y)
460	DEALLOCATE(dJ_dx)
461	DEALLOCATE(x_incr)
462	DEALLOCATE(diag_d2J)
463	DEALLOCATE(d2J_dx2)
464	DEALLOCATE(KO)
465	DEALLOCATE(Bm1)
466	DEALLOCATE(Om1)
467
468	CALL message_set_routine(routine)
469
470	END SUBROUTINE ropp_1dvar_levmarq_bangle
471
472	!***
473
474	!-------------------------------------------------------------------------------
475	! 2. Refractivity
476	!-------------------------------------------------------------------------------
477
478	SUBROUTINE ropp_1dvar_levmarq_refrac(obs, bg, state, config, diag)
479
480	! 2.1 Declarations
481	! ----------------
482
483	USE typesizes, ONLY: wp => EightByteReal
484	USE ropp_utils
485	USE ropp_fm
486	USE ropp_1dvar, not_this => ropp_1dvar_levmarq_refrac
487	USE matrix
488
489	IMPLICIT NONE
490
491	TYPE(Obs1dRefrac), INTENT(inout) :: obs ! Observation data
492	TYPE(State1dFM), INTENT(inout) :: bg ! Background data
493	TYPE(State1dFM), INTENT(inout) :: state ! State vector
494	TYPE(VarConfig), INTENT(in) :: config ! Configuration options
495	TYPE(VarDiag), INTENT(inout) :: diag ! Diagnostic output
496
497	REAL(wp) :: J ! Cost function value
498	TYPE(State1dFM) :: x ! Control vector
499	TYPE(State1dFM) :: xmin ! Control at minimum J
500	TYPE(Obs1dRefrac) :: y ! Forward model obs
501	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: K ! K-matrix
502	REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_x ! Change of state
503	REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_y ! Change of observation
504	REAL(wp), DIMENSION(:), ALLOCATABLE :: dJ_dx ! Cost function gradient
505	REAL(wp), DIMENSION(:), ALLOCATABLE :: diag_d2J ! Diagonal d2J/dx2
506	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: d2J_dx2 ! 2nd derivative cost
507	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: KO ! K O^-1
508	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: Bm1 ! B^-1
509	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: Om1 ! O^-1
510	REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_J ! Change of cost fn
511	REAL(wp), DIMENSION(:), ALLOCATABLE :: state_last ! Previous state vector
512	REAL(wp), DIMENSION(:), ALLOCATABLE :: state_sigma ! State std deviation
513	REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_state ! Change of state vector
514	REAL(wp) :: J_last ! Previous cost function
515	REAL(wp) :: J_min ! Minimum cost function
516	INTEGER :: i_pointer ! Value index
517	INTEGER :: i ! Counter
518	INTEGER :: n_iter ! Number of iterations
519	INTEGER :: nobs ! Number of obserations
520	INTEGER :: nstate ! No. of state elements
521	LOGICAL :: marq ! Flag to L-M minimise
522	REAL(wp) :: lambda ! Iteration factor
523
524	CHARACTER(len = 4) :: it_str
525	CHARACTER(len = 15) :: ch_str, co_str
526	CHARACTER(len = 256) :: routine
527
528	! 2.2 Message handling
529	! --------------------
530
531	CALL message_get_routine(routine)
532	CALL message_set_routine('ropp_1dvar_levmarq')
533
534	! 2.3 Initialise rolling buffers and pointer for convergence checks
535	! -----------------------------------------------------------------
536
537	i_pointer = 0
538	n_iter = 0
539	lambda = 1.0E-4_wp
540	marq = .FALSE.
541	J = 0.0_wp
542	J_last = 1.0E30_wp
543	J_min = J_last
544
545	nstate = SIZE(bg%state)
546	nobs = SIZE(obs%refrac)
547
548	ALLOCATE(K(nobs, nstate))
549	ALLOCATE(delta_x(nstate))
550	ALLOCATE(delta_y(nobs))
551	ALLOCATE(dJ_dx(nstate))
552	ALLOCATE(diag_d2J(nstate))
553	ALLOCATE(d2J_dx2(nstate,nstate))
554	ALLOCATE(KO(nstate,nobs))
555	ALLOCATE(Bm1(nstate,nstate))
556	ALLOCATE(Om1(nobs,nobs))
557
558	Bm1 = matrix_invert(bg%cov)
559	Om1 = matrix_invert(obs%cov)
560
561	IF (ALLOCATED(delta_J)) DEALLOCATE(delta_J)
562	ALLOCATE(delta_J(config%conv_check_n_previous))
563	delta_J(:) = 0.0_wp
564
565	IF (ALLOCATED(state_last)) DEALLOCATE(state_last)
566	ALLOCATE(state_last(SIZE(bg%state)))
567	state_last(:) = 0.0_wp
568
569	IF (ALLOCATED(state_sigma)) DEALLOCATE(state_sigma)
570	ALLOCATE(state_sigma(SIZE(bg%state)))
571	DO i = 1, SIZE(state_sigma)
572	state_sigma(i) = SQRT(bg%cov%d(i + i*(i-1)/2)) ! Direct read from matrix_pp
573	ENDDO
574
575	IF (ALLOCATED(delta_state)) DEALLOCATE(delta_state)
576	ALLOCATE(delta_state(config%conv_check_n_previous))
577	delta_state(:) = 0.0_wp
578
579	! 2.4 First guess state
580	! ---------------------
581
582	x = bg
583	xmin = bg
584
585	! 2.5 Main minimisation iteration loop
586	! ------------------------------------
587
588	DO WHILE (J <= J_min)
589
590	! 2.5.1 Update iteration counter
591
592	n_iter = n_iter + 1
593
594	! 2.6 Compute cost function
595	! -------------------------
596
597	! 2.6.1 Calculate pseudo observations
598	! -----------------------------------
599
600	ALLOCATE(y%refrac(SIZE(obs%refrac)))
601
602	CALL copy_alloc(obs % geop, y % geop)
603	CALL copy_alloc(obs % weights, y % weights)
604
605	! 2.6.2 Forward model
606	! -------------------
607
608	IF(ASSOCIATED(x%ak))THEN
609	CALL ropp_fm_state2state_ecmwf(x)
610	ELSE
611	CALL ropp_fm_state2state_meto(x)
612	ENDIF
613
614	IF (x%new_ref_op) THEN
615	CALL ropp_fm_refrac_1d_new(x, y)
616	ELSE
617	CALL ropp_fm_refrac_1d(x, y)
618	END IF
619
620	! 2.6.3 Compute cost function
621	! ---------------------------
622
623	delta_x = (x%state - bg%state)
624	delta_y = (y%refrac - obs%refrac) * y%weights
625
626	J = 0.5_wp * DOT_PRODUCT(delta_y, MATMUL(Om1, delta_y)) &
627	+ 0.5_wp * DOT_PRODUCT(delta_x, MATMUL(Bm1, delta_x))
628
629	IF (n_iter == 1) diag % J_init = J
630
631	IF (J > J_last) marq = .TRUE.
632
633	! 2.7 Levenberg-Marquardt minimisation
634	! ------------------------------------
635
636	IF(.NOT. marq)THEN
637
638	! 2.7.1 Normal Newtonian iteration
639	! --------------------------------
640
641	lambda = 0.1_wp * lambda
642	marq = .FALSE.
643
644	! 2.7.2 Evaluate K gradient matrix for current x
645	! ----------------------------------------------
646
647	CALL ropp_fm_refrac_1d_grad(x, y, K)
648
649	! 2.7.3 Calculate -dJ_dx vector and d2J_dx2 matrix at x
650	! -----------------------------------------------------
651
652	delta_x = x%state - bg%state
653	delta_y = (obs%refrac - y%refrac) * y%weights
654
655	WHERE(ABS(delta_y) > 50.0_wp)
656	delta_y = 0.0_wp
657	END WHERE
658
659	KO = MATMUL(TRANSPOSE(K), Om1)
660	dJ_dx = MATMUL(KO, delta_y) - MATMUL(Bm1, delta_x)
661	d2J_dx2 = Bm1 + MATMUL(KO, K)
662
663	! 2.7.4 Store inverse of solution covariance matrix
664	! -------------------------------------------------
665
666	DO i=1,SIZE(bg%state)
667	diag_d2J(i) = d2J_dx2(i,i)
668	ENDDO
669
670	ELSE
671
672	! 2.7.5 Levenberg-Marquardt iteration
673	! -----------------------------------
674	! The previous increment increased the value of the penalty function.
675	! Use previous values of -dJ_dx, d2J_dx2 and adjust value of lambda
676
677	marq = .TRUE.
678	lambda = 100.0_wp * lambda
679
680	ENDIF
681
682	! 2.7.6 Levenberg-Marquardt adjustment to diagonal terms
683	! ------------------------------------------------------
684
685	DO i=1, SIZE(bg%state)
686	d2J_dx2(i,i) = diag_d2J(i) * (lambda+1.0_wp)
687	ENDDO
688
689	! 2.7.7 Solve matrix equation d2J_dx2 . dx = -dJ_dx and update state
690	! ------------------------------------------------------------------
691
692	if (x%use_logq) then
693	delta_x = matrix_solve(d2J_dx2, dJ_dx)
694	x%state = x%state + SIGN( MIN(ABS(delta_x),ABS(x%state/2.0_wp)), &
695	delta_x)
696	else
697	x%state = x%state + matrix_solve(d2J_dx2, dJ_dx)
698	endif
699
700	! 2.8 Update state vector variables - model dependent
701	! ---------------------------------
702
703	IF(ASSOCIATED(x%ak))THEN
704	CALL ropp_fm_state2state_ecmwf(x)
705	ELSE
706	CALL ropp_fm_state2state_meto(x)
707	ENDIF
708
709	! 2.9 Keep track of current state
710	! -------------------------------
711
712	IF( J_min > J ) THEN
713	J_min = J
714	xmin = x
715	ENDIF
716
717	IF (config % conv_check_apply) THEN
718
719	i_pointer = MOD(i_pointer + 1, config % conv_check_n_previous)
720	IF (i_pointer == 0) i_pointer = config % conv_check_n_previous
721
722	delta_J(i_pointer) = J_last - J
723	J_last = J
724	delta_state(i_pointer) = MAXVAL(ABS(state_last - x%state)/state_sigma)
725	state_last = x%state
726
727	ENDIF
728
729	! 2.10 Check for convergence
730	! --------------------------
731
732	IF (config % conv_check_apply) THEN
733
734	IF (config%minropp%impres == 0)THEN
735	WRITE(it_str, '(i4)') n_iter
736	WRITE(ch_str, '(g15.5)') J
737	ch_str = ADJUSTL(ch_str)
738	IF (n_iter > 1) THEN
739	WRITE(co_str, '(g15.5)') delta_state(i_pointer)
740	co_str = ADJUSTL(co_str)
741	ELSE
742	co_str = ' - '
743	ENDIF
744
745	CALL message(msg_cont, &
746	' n_iter = ' // it_str // ' J = ' // ch_str // &
747	' max(relative change in state) = ' // co_str)
748	ENDIF
749
750
751	IF (MAXVAL(delta_state) < config%conv_check_max_delta_state &
752	.AND. n_iter > config % conv_check_n_previous) THEN
753	WRITE(ch_str, '(g15.5)') config%conv_check_max_delta_state
754	ch_str = ADJUSTL(ch_str)
755	WRITE(it_str, '(i2)') config%conv_check_n_previous
756	it_str = ADJUSTL(it_str)
757	CALL message(msg_cont, '')
758	CALL message(msg_info, &
759	'Convergence assumed to be achieved as the state vector did ' // &
760	'not change by more\n ' // 'than ' // TRIM(ch_str) // ' ' // &
761	'relative to the assumed background errors for the last ' // &
762	TRIM(it_str) // ' iterations.\n')
763	EXIT
764	ELSE IF (MAXVAL(ABS(delta_J)) < config%conv_check_max_delta_J &
765	.AND. n_iter > config % conv_check_n_previous) THEN
766	WRITE(ch_str, '(g15.5)') config%conv_check_max_delta_J
767	ch_str = ADJUSTL(ch_str)
768	WRITE(it_str, '(i2)') config%conv_check_n_previous
769	it_str = ADJUSTL(it_str)
770	CALL message(msg_cont, '')
771	CALL message(msg_info, &
772	'Convergence assumed to be achieved as the cost function did ' // &
773	'not change by more\n ' // 'than ' // TRIM(ch_str) // &
774	' for the last ' // TRIM(it_str) // ' iterations.\n')
775	EXIT
776	ENDIF
777
778	ENDIF
779
780	ENDDO ! end main iteration loop
781
782	! 2.11 Copy solution back to state
783	! --------------------------------
784
785	state = xmin
786
787	! 2.12 Diagnostic data
788	! --------------------
789
790	diag % J = J_min
791
792	IF (COUNT(obs % weights > 0.0_wp) > 0) THEN
793	diag % J_scaled = 2.0_wp * J_min / REAL(COUNT(obs % weights > 0.0_wp), wp)
794	ENDIF
795
796	diag % n_iter = n_iter
797
798	ALLOCATE (diag%J_bgr(SIZE(state%state)))
799	delta_x = state%state - bg%state
800	diag % J_bgr = 0.5_wp * delta_x * matrix_solve(bg%cov, delta_x)
801
802	! 2.13 Clean up
803	! -------------
804
805	DEALLOCATE(K)
806	DEALLOCATE(delta_x)
807	DEALLOCATE(delta_y)
808	DEALLOCATE(dJ_dx)
809	DEALLOCATE(diag_d2J)
810	DEALLOCATE(d2J_dx2)
811	DEALLOCATE(KO)
812	DEALLOCATE(Bm1)
813	DEALLOCATE(Om1)
814
815	CALL message_set_routine(routine)
816
817	END SUBROUTINE ropp_1dvar_levmarq_refrac

Download in other formats:

Original format