Context navigation

Back to Ticket #646

Ticket #646: ropp_1dvar_levmarq.f90.r6123

File ropp_1dvar_levmarq.f90.r6123, 28.7 KB (added by Ian Culverwell, 6 years ago)
ropp_1dvar_levmarq.f90.r6123

Line
1	! $Id: ropp_1dvar_levmarq.f90 6213 2020-03-05 10:02:39Z 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	! W.H. Press, S.A. Teukolsjy, W.T. Vetterling and B.P. Flannery,
46	! Numerical Recipes in C - The Art of Scientific Computing.
47	! 2nd Ed., Cambridge University Press, 1992.
48	!
49	! AUTHOR
50	! Met Office, Exeter, UK.
51	! Any comments on this software should be given via the ROM SAF
52	! Helpdesk at http://www.romsaf.org
53	!
54	! COPYRIGHT
55	! (c) EUMETSAT. All rights reserved.
56	! For further details please refer to the file COPYRIGHT
57	! which you should have received as part of this distribution.
58	!
59	!****
60
61	!-------------------------------------------------------------------------------
62	! 1. Bending angle
63	!-------------------------------------------------------------------------------
64
65	SUBROUTINE ropp_1dvar_levmarq_bangle(obs, bg, state, config, diag)
66
67	! 1.1 Declarations
68	! ----------------
69
70	USE typesizes, ONLY: wp => EightByteReal
71	USE ropp_utils
72	USE ropp_fm
73	USE ropp_1dvar, not_this => ropp_1dvar_levmarq_bangle
74	USE matrix
75
76	IMPLICIT NONE
77
78	TYPE(Obs1dBangle), INTENT(inout) :: obs ! Observation data
79	TYPE(State1dFM), INTENT(inout) :: bg ! Background data
80	TYPE(State1dFM), INTENT(inout) :: state ! State vector
81	TYPE(VarConfig), INTENT(in) :: config ! Configuration options
82	TYPE(VarDiag), INTENT(inout) :: diag ! Diagnostic output
83
84	REAL(wp) :: J ! Cost function value
85	TYPE(State1dFM) :: x ! Control vector
86	TYPE(State1dFM) :: x_old ! Control for LM parameter testing
87	TYPE(Obs1dBangle) :: y ! Forward model obs
88	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: K ! K-matrix
89	REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_x ! Change of state
90	REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_y ! Change of observation
91	REAL(wp), DIMENSION(:), ALLOCATABLE :: dJ_dx ! Cost function gradient
92	REAL(wp), DIMENSION(:), ALLOCATABLE :: diag_d2J ! Diagonal d2J/dx2
93	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: d2J_dx2 ! 2nd derivative cost fn
94	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: KO ! K O^-1
95	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: Bm1 ! B^-1
96	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: Om1 ! O^-1
97	REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_J ! Change of cost fn
98	REAL(wp), DIMENSION(:), ALLOCATABLE :: state_last ! Previous state vector
99	REAL(wp), DIMENSION(:), ALLOCATABLE :: state_sigma ! State std deviation
100	REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_state ! Change of state vector
101	REAL(wp) :: J_last ! Previous cost function
102	REAL(wp) :: J_min ! Minimum cost function
103	INTEGER :: i_pointer ! Value index
104	INTEGER :: i ! Counter
105	INTEGER :: n_iter ! Number of iterations
106	INTEGER :: n_iter_max ! Maximum number of iterations
107	INTEGER :: n_grad ! Number of gradient updates
108	INTEGER :: nobs ! Number of observations
109	INTEGER :: nstate ! No. of state elements
110	REAL(wp) :: lambda ! Levenberg-Marquardt factor
111	REAL(wp) :: lambda_max ! Maximum Levenberg-Marquardt factor
112
113	CHARACTER(len = 4) :: it_str
114	CHARACTER(len = 4) :: gr_str
115	CHARACTER(len = 15) :: ch_str, co_str
116	CHARACTER(len = 256) :: routine
117
118	! 1.2 Message handling
119	! --------------------
120
121	CALL message_get_routine(routine)
122	CALL message_set_routine('ropp_1dvar_levmarq_bangle')
123
124	! 1.3 Initialise variables
125	! ------------------------
126
127	i_pointer = 0
128	n_grad = 0
129	n_iter = 0
130	n_iter_max = 50
131
132	lambda = 1.0E-4_wp
133	IF (bg%direct_ion) lambda = 1.0E-2_wp
134	lambda_max = 1.0E10_wp
135
136	J = 0.0_wp
137	J_last = 1.0E30_wp
138	J_min = J_last
139
140	nstate = SIZE(bg%state)
141	nobs = SIZE(obs%bangle)
142
143	ALLOCATE(K(nobs, nstate))
144	ALLOCATE(delta_x(nstate))
145	ALLOCATE(delta_y(nobs))
146	ALLOCATE(dJ_dx(nstate))
147	ALLOCATE(diag_d2J(nstate))
148	ALLOCATE(d2J_dx2(nstate,nstate))
149	ALLOCATE(KO(nstate,nobs))
150	ALLOCATE(Bm1(nstate,nstate))
151	ALLOCATE(Om1(nobs,nobs))
152
153	IF (ALLOCATED(delta_J)) DEALLOCATE(delta_J)
154	ALLOCATE(delta_J(config%conv_check_n_previous))
155	delta_J(:) = 0.0_wp
156
157	IF (ALLOCATED(state_last)) DEALLOCATE(state_last)
158	ALLOCATE(state_last(SIZE(bg%state)))
159	state_last(:) = 0.0_wp
160
161	IF (ALLOCATED(state_sigma)) DEALLOCATE(state_sigma)
162	ALLOCATE(state_sigma(SIZE(bg%state)))
163	DO i = 1, SIZE(state_sigma)
164	state_sigma(i) = SQRT(bg%cov%d(i + i*(i-1)/2)) ! Direct read from matrix_pp
165	ENDDO
166
167	IF (ALLOCATED(delta_state)) DEALLOCATE(delta_state)
168	ALLOCATE(delta_state(config%conv_check_n_previous))
169	delta_state(:) = 0.0_wp
170
171	! 1.4 Inverse error covariances
172	! -----------------------------
173
174	Bm1 = matrix_invert(bg%cov)
175	Om1 = matrix_invert(obs%cov)
176
177	! 1.5 First guess and pseudo-observations
178	! ---------------------------------------
179
180	x = bg
181	y = obs
182
183	! 1.6 Initial cost function
184	! -------------------------
185
186	IF (ASSOCIATED(x%ak)) THEN
187	CALL ropp_fm_state2state_ecmwf(x)
188	ELSE
189	CALL ropp_fm_state2state_meto(x)
190	END IF
191
192	CALL ropp_fm_bangle_1d(x, y)
193
194	delta_x = (x%state - bg%state)
195	delta_y = (y%bangle - obs%bangle) * y%weights
196
197	J = 0.5_wp * DOT_PRODUCT(delta_y, MATMUL(Om1, delta_y)) + &
198	0.5_wp * DOT_PRODUCT(delta_x, MATMUL(Bm1, delta_x))
199
200	diag%J_init = J
201	J_last = J
202	J_min = J
203
204	IF (config%minropp%impres == 0) THEN
205	WRITE(it_str, '(i4)') n_iter
206	WRITE(ch_str, '(g15.5)') J
207	ch_str = ADJUSTL(ch_str)
208	co_str = ' - '
209
210	CALL message(msg_cont, &
211	' n_iter = ' // it_str // ' J = ' // ch_str // &
212	' max(relative change in state) = ' // co_str)
213	END IF
214
215	! 1.7 Main minimisation iteration loop
216	! ------------------------------------
217
218	DO WHILE (n_iter < n_iter_max)
219
220	! 1.7.1 Bookkeeping
221
222	n_grad = n_grad + 1
223
224	! 1.7.2 Compute gradient and Hessian of cost function
225
226	CALL ropp_fm_bangle_1d_grad(x, y, K)
227
228	!TODO: (a) Do we need this?
229	! (b) Does it change the cost function?
230	! (c) When would a change in bending angle ever exceed 50 rad?
231	! Was this lifted from the refrac code? (delta_N = 50 N-units just about plausible.)
232	WHERE (ABS(delta_y) > 50.0_wp)
233	delta_y = 0.0_wp
234	END WHERE
235
236	KO = MATMUL(TRANSPOSE(K), Om1)
237	dJ_dx = MATMUL(KO, delta_y) + MATMUL(Bm1, delta_x)
238	d2J_dx2 = Bm1 + MATMUL(KO, K)
239
240	! 1.7.3 Levenberg-Marquardt adjustment of Hessian diagonal
241
242	DO i = 1, nstate
243	d2J_dx2(i,i) = d2J_dx2(i,i) * (1.0_wp + lambda)
244	END DO
245
246	! 1.7.4 Solve matrix equation d2J_dx2 . dx = -dJ_dx
247
248	delta_x = matrix_solve(d2J_dx2, - dJ_dx)
249
250	! 1.7.5 Update test state vector
251
252	x_old = x
253
254	IF (x%use_logq) THEN
255	x%state = x%state + SIGN(MIN(ABS(delta_x), ABS(x%state/2.0_wp)), delta_x)
256	ELSE
257	x%state = x%state + delta_x
258	END IF
259
260	! 1.7.6 Special handling of ionospheric state vector elements
261
262	IF (bg%direct_ion) THEN
263
264	i = nstate - 2
265	IF (x%state(i) < ropp_ZERO) THEN
266	CALL message(msg_warn, "Levenberg-Marquardt solver returns " // &
267	"Ne_max < 0 ... suggest examining final value. \n")
268	END IF
269
270	i = nstate - 1
271	IF (x%state(i) < 0.01_wp*bg%state(i)) THEN
272	CALL message(msg_warn, "Levenberg-Marquardt solver returns " // &
273	"H_peak < 1% of background ... resetting to background value. \n")
274	x%state(i) = bg%state(i)
275	END IF
276
277	i = nstate
278	IF (x%state(i) < 0.01_wp*bg%state(i)) THEN
279	CALL message(msg_warn, "Levenberg-Marquardt solver returns " // &
280	"H_width < 1% of background ... resetting to background value. \n")
281	x%state(i) = bg%state(i)
282	END IF
283
284	END IF
285
286	! 1.7.7 Compute cost function
287
288	IF (ASSOCIATED(x%ak)) THEN
289	CALL ropp_fm_state2state_ecmwf(x)
290	ELSE
291	CALL ropp_fm_state2state_meto(x)
292	END IF
293
294	CALL ropp_fm_bangle_1d(x, y)
295
296	delta_x = (x%state - bg%state)
297	delta_y = (y%bangle - obs%bangle) * y%weights
298
299	J = 0.5_wp * DOT_PRODUCT(delta_y, MATMUL(Om1, delta_y)) + &
300	0.5_wp * DOT_PRODUCT(delta_x, MATMUL(Bm1, delta_x))
301
302	! 1.8 Levenberg-Marquardt update and convergence check
303	! ----------------------------------------------------
304
305	IF ( J >= (J_last + config%conv_check_max_delta_state) ) THEN ! Keep previous state vector, increase lambda, repeat
306
307	! 1.8.1 Keep old state, increase lambda, bail out if it's too big
308
309	x = x_old
310	lambda = 10.0_wp * lambda
311
312	IF (ABS(lambda) > lambda_max) THEN
313
314	WRITE(ch_str, '(g15.5)') lambda_max
315	ch_str = ADJUSTL(ch_str)
316	CALL message(msg_cont, '')
317	CALL message(msg_info, &
318	'Levenberg-Marquardt parameter exceeded max (=' // TRIM(ch_str) // &
319	'); convergence is unlikely.\n ' // &
320	'Check convergence parameters; conv_check_max_delta_J might ' // &
321	'be too demanding.\n')
322	EXIT
323
324	END IF
325
326	! 1.8.2 Logging
327
328	!TODO: Add a configuration option for lambda update outputs
329	IF (config%minropp%impres == 0) THEN
330
331	WRITE(it_str, '(i4)') n_iter
332	WRITE(ch_str, '(g15.5)') J
333	ch_str = ADJUSTL(ch_str)
334	IF (n_iter > 0) THEN
335	WRITE(co_str, '(g15.5)') lambda
336	co_str = ADJUSTL(co_str)
337	ELSE
338	co_str = ' - '
339	END IF
340
341	CALL message(msg_cont, &
342	' n_iter = ' // it_str // ' J = ' // ch_str // &
343	' lambda -> ' // co_str)
344
345	END IF
346
347	! 1.8.3 Next iteration
348
349	CYCLE
350
351	ELSE ! Decrease lambda, check convergence, continue iteration
352
353	! 1.8.4 Decrease lambda
354
355	lambda = 0.1_wp *lambda
356
357	! 1.8.5 Keep track of current state
358
359	n_iter = n_iter + 1
360
361	IF (config % conv_check_apply) THEN
362
363	i_pointer = MOD(i_pointer + 1, config % conv_check_n_previous)
364	IF (i_pointer == 0) i_pointer = config % conv_check_n_previous
365
366	delta_J(i_pointer) = J_last - J
367	delta_state(i_pointer) = MAXVAL(ABS(state_last - x%state)/state_sigma)
368
369	state_last = x%state
370
371	J_last = J
372
373	END IF
374
375	IF (J < J_min) THEN
376	J_min = J
377	END IF
378
379	! 1.8.6 Check for convergence
380
381	IF (config % conv_check_apply) THEN
382
383	IF (config%minropp%impres == 0) THEN
384	WRITE(it_str, '(i4)') n_iter
385	WRITE(ch_str, '(g15.5)') J
386	ch_str = ADJUSTL(ch_str)
387	IF (n_iter > 0) THEN
388	WRITE(co_str, '(g15.5)') delta_state(i_pointer)
389	co_str = ADJUSTL(co_str)
390	ELSE
391	co_str = ' - '
392	END IF
393
394	CALL message(msg_cont, &
395	' n_iter = ' // it_str // ' J = ' // ch_str // &
396	' max(relative change in state) = ' // co_str)
397	END IF
398
399	IF (MAXVAL(delta_state) < config%conv_check_max_delta_state &
400	.AND. n_iter > config % conv_check_n_previous) THEN
401
402	WRITE(ch_str, '(g15.5)') config%conv_check_max_delta_state
403	ch_str = ADJUSTL(ch_str)
404	WRITE(it_str, '(i2)') config%conv_check_n_previous
405	it_str = ADJUSTL(it_str)
406	CALL message(msg_cont, '')
407	CALL message(msg_info, &
408	'Convergence assumed to be achieved as the state vector did ' // &
409	'not change by more\n ' // 'than ' // TRIM(ch_str) // ' ' // &
410	'relative to the assumed background errors for the last ' // &
411	TRIM(it_str) // ' iterations.\n')
412	EXIT
413
414	ELSE IF (MAXVAL(ABS(delta_J)) < config%conv_check_max_delta_J &
415	.AND. n_iter > config % conv_check_n_previous) THEN
416
417	WRITE(ch_str, '(g15.5)') config%conv_check_max_delta_J
418	ch_str = ADJUSTL(ch_str)
419	WRITE(it_str, '(i2)') config%conv_check_n_previous
420	it_str = ADJUSTL(it_str)
421	CALL message(msg_cont, '')
422	CALL message(msg_info, &
423	'Convergence assumed to be achieved as the cost function did ' // &
424	'not change by more\n ' // 'than ' // TRIM(ch_str) // &
425	' for the last ' // TRIM(it_str) // ' iterations.\n')
426	EXIT
427
428	ELSE IF (ABS(lambda) > lambda_max) THEN
429
430	WRITE(ch_str, '(g15.5)') lambda_max
431	ch_str = ADJUSTL(ch_str)
432	CALL message(msg_cont, '')
433	CALL message(msg_info, &
434	'Levenberg-Marquardt parameter exceeded max (=' // TRIM(ch_str) // &
435	'); convergence is unlikely.\n ' // &
436	'Check convergence parameters; conv_check_max_delta_J might ' // &
437	'be too demanding.\n')
438	EXIT
439
440	END IF
441
442	END IF ! Convergence check
443
444	END IF ! Levenberg-Marquardt update
445
446	END DO ! End main iteration loop
447
448	IF (n_iter < n_iter_max) THEN
449	WRITE(it_str, '(i4)') n_iter ; it_str = ADJUSTL(it_str)
450	WRITE(gr_str, '(i4)') n_grad ; gr_str = ADJUSTL(gr_str)
451	CALL message(msg_info, 'Finished after ' // TRIM(it_str) // ' iterations (' // &
452	TRIM(gr_str) // ' forward model / gradient evaluations).\n')
453	ELSE
454	WRITE(it_str, '(i4)') n_iter_max ; it_str = ADJUSTL(it_str)
455	CALL message(msg_warn, &
456	'Iteration ended after n_iter_max = ' // TRIM(it_str) // &
457	' iterations without achieving convergence.\n')
458	END IF
459
460	! 1.9 Copy solution back to state
461	! -------------------------------
462
463	state = x
464
465	! 1.10 Diagnostic data
466	! --------------------
467
468	diag%J = J_min
469
470	IF (COUNT(obs%weights > 0.0_wp) > 0) THEN
471	diag%J_scaled = 2.0_wp * J_min / REAL(COUNT(obs%weights > 0.0_wp), wp)
472	ENDIF
473
474	diag%n_iter = n_iter
475
476	ALLOCATE (diag%J_bgr(SIZE(state%state)))
477	delta_x = state%state - bg%state
478	diag%J_bgr = 0.5_wp * delta_x * matrix_solve(bg%cov, delta_x)
479
480	! 1.11 Clean up
481	! -------------
482
483	DEALLOCATE(K)
484	DEALLOCATE(delta_x)
485	DEALLOCATE(delta_y)
486	DEALLOCATE(dJ_dx)
487	DEALLOCATE(diag_d2J)
488	DEALLOCATE(d2J_dx2)
489	DEALLOCATE(KO)
490	DEALLOCATE(Bm1)
491	DEALLOCATE(Om1)
492
493	CALL message_set_routine(routine)
494
495	END SUBROUTINE ropp_1dvar_levmarq_bangle
496
497	!***
498
499	!-------------------------------------------------------------------------------
500	! 2. Refractivity
501	!-------------------------------------------------------------------------------
502
503	SUBROUTINE ropp_1dvar_levmarq_refrac(obs, bg, state, config, diag)
504
505	! 2.1 Declarations
506	! ----------------
507
508	USE typesizes, ONLY: wp => EightByteReal
509	USE ropp_utils
510	USE ropp_fm
511	USE ropp_1dvar, not_this => ropp_1dvar_levmarq_refrac
512	USE matrix
513
514	IMPLICIT NONE
515
516	TYPE(Obs1dRefrac), INTENT(inout) :: obs ! Observation data
517	TYPE(State1dFM), INTENT(inout) :: bg ! Background data
518	TYPE(State1dFM), INTENT(inout) :: state ! State vector
519	TYPE(VarConfig), INTENT(in) :: config ! Configuration options
520	TYPE(VarDiag), INTENT(inout) :: diag ! Diagnostic output
521
522	REAL(wp) :: J ! Cost function value
523	TYPE(State1dFM) :: x ! Control vector
524	TYPE(State1dFM) :: x_old ! Control for LM parameter testing
525	TYPE(Obs1dRefrac) :: y ! Forward model obs
526	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: K ! K-matrix
527	REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_x ! Change of state
528	REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_y ! Change of observation
529	REAL(wp), DIMENSION(:), ALLOCATABLE :: dJ_dx ! Cost function gradient
530	REAL(wp), DIMENSION(:), ALLOCATABLE :: diag_d2J ! Diagonal d2J/dx2
531	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: d2J_dx2 ! 2nd derivative cost fn
532	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: KO ! K O^-1
533	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: Bm1 ! B^-1
534	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: Om1 ! O^-1
535	REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_J ! Change of cost fn
536	REAL(wp), DIMENSION(:), ALLOCATABLE :: state_last ! Previous state vector
537	REAL(wp), DIMENSION(:), ALLOCATABLE :: state_sigma ! State std deviation
538	REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_state ! Change of state vector
539	REAL(wp) :: J_last ! Previous cost function
540	REAL(wp) :: J_min ! Minimum cost function
541	INTEGER :: i_pointer ! Value index
542	INTEGER :: i ! Counter
543	INTEGER :: n_iter ! Number of iterations
544	INTEGER :: n_iter_max ! Maximum number of iterations
545	INTEGER :: n_grad ! Number of gradient updates
546	INTEGER :: nobs ! Number of observations
547	INTEGER :: nstate ! No. of state elements
548	REAL(wp) :: lambda ! Levenberg-Marquardt factor
549	REAL(wp) :: lambda_max ! Maximum Levenberg-Marquardt factor
550
551	CHARACTER(len = 4) :: it_str
552	CHARACTER(len = 4) :: gr_str
553	CHARACTER(len = 15) :: ch_str, co_str
554	CHARACTER(len = 256) :: routine
555
556	! 2.2 Message handling
557	! --------------------
558
559	CALL message_get_routine(routine)
560	CALL message_set_routine('ropp_1dvar_levmarq_refrac')
561
562	! 2.3 Initialise variables
563	! ------------------------
564
565	i_pointer = 0
566	n_grad = 0
567	n_iter = 0
568	n_iter_max = 50
569
570	lambda = 1.0E-4_wp
571	lambda_max = 1.0E10_wp
572
573	J = 0.0_wp
574	J_last = 1.0E30_wp
575	J_min = J_last
576
577	nstate = SIZE(bg%state)
578	nobs = SIZE(obs%refrac)
579
580	ALLOCATE(K(nobs, nstate))
581	ALLOCATE(delta_x(nstate))
582	ALLOCATE(delta_y(nobs))
583	ALLOCATE(dJ_dx(nstate))
584	ALLOCATE(diag_d2J(nstate))
585	ALLOCATE(d2J_dx2(nstate,nstate))
586	ALLOCATE(KO(nstate,nobs))
587	ALLOCATE(Bm1(nstate,nstate))
588	ALLOCATE(Om1(nobs,nobs))
589
590	IF (ALLOCATED(delta_J)) DEALLOCATE(delta_J)
591	ALLOCATE(delta_J(config%conv_check_n_previous))
592	delta_J(:) = 0.0_wp
593
594	IF (ALLOCATED(state_last)) DEALLOCATE(state_last)
595	ALLOCATE(state_last(SIZE(bg%state)))
596	state_last(:) = 0.0_wp
597
598	IF (ALLOCATED(state_sigma)) DEALLOCATE(state_sigma)
599	ALLOCATE(state_sigma(SIZE(bg%state)))
600	DO i = 1, SIZE(state_sigma)
601	state_sigma(i) = SQRT(bg%cov%d(i + i*(i-1)/2)) ! Direct read from matrix_pp
602	ENDDO
603
604	IF (ALLOCATED(delta_state)) DEALLOCATE(delta_state)
605	ALLOCATE(delta_state(config%conv_check_n_previous))
606	delta_state(:) = 0.0_wp
607
608	! 2.4 Inverse error covariances
609	! -----------------------------
610
611	Bm1 = matrix_invert(bg%cov)
612	Om1 = matrix_invert(obs%cov)
613
614	! 2.5 First guess and pseudo-observations
615	! ---------------------------------------
616
617	x = bg
618	y = obs
619
620	! 2.6 Initial cost function
621	! -------------------------
622
623	IF (ASSOCIATED(x%ak)) THEN
624	CALL ropp_fm_state2state_ecmwf(x)
625	ELSE
626	CALL ropp_fm_state2state_meto(x)
627	END IF
628
629	IF (x%new_ref_op) THEN
630	CALL ropp_fm_refrac_1d_new(x, y)
631	ELSE
632	CALL ropp_fm_refrac_1d(x, y)
633	END IF
634
635	delta_x = (x%state - bg%state)
636	delta_y = (y%refrac - obs%refrac) * y%weights
637
638	J = 0.5_wp * DOT_PRODUCT(delta_y, MATMUL(Om1, delta_y)) + &
639	0.5_wp * DOT_PRODUCT(delta_x, MATMUL(Bm1, delta_x))
640
641	diag%J_init = J
642	J_last = J
643	J_min = J
644
645	IF (config%minropp%impres == 0) THEN
646	WRITE(it_str, '(i4)') n_iter
647	WRITE(ch_str, '(g15.5)') J
648	ch_str = ADJUSTL(ch_str)
649	co_str = ' - '
650
651	CALL message(msg_cont, &
652	' n_iter = ' // it_str // ' J = ' // ch_str // &
653	' max(relative change in state) = ' // co_str)
654	END IF
655
656	! 2.7 Main minimisation iteration loop
657	! ------------------------------------
658
659	DO WHILE (n_iter < n_iter_max)
660
661	! 2.7.1 Bookkeeping
662
663	n_grad = n_grad + 1
664
665	! 2.7.2 Compute gradient and Hessian of cost function
666
667	CALL ropp_fm_refrac_1d_grad(x, y, K)
668
669	!TODO: (a) Do we need this?
670	! (b) Does it change the cost function?
671	WHERE (ABS(delta_y) > 50.0_wp)
672	delta_y = 0.0_wp
673	END WHERE
674
675	KO = MATMUL(TRANSPOSE(K), Om1)
676	dJ_dx = MATMUL(KO, delta_y) + MATMUL(Bm1, delta_x)
677	d2J_dx2 = Bm1 + MATMUL(KO, K)
678
679	! 2.7.3 Levenberg-Marquardt adjustment of Hessian diagonal
680
681	DO i = 1, nstate
682	d2J_dx2(i,i) = d2J_dx2(i,i) * (1.0_wp + lambda)
683	END DO
684
685	! 2.7.4 Solve matrix equation d2J_dx2 . dx = -dJ_dx
686
687	delta_x = matrix_solve(d2J_dx2, - dJ_dx)
688
689	! 2.7.5 Update test state vector
690
691	x_old = x
692
693	IF (x%use_logq) THEN
694	x%state = x%state + SIGN(MIN(ABS(delta_x), ABS(x%state/2.0_wp)), delta_x)
695	ELSE
696	x%state = x%state + delta_x
697	END IF
698
699	! 2.7.6 Compute cost function
700
701	IF (ASSOCIATED(x%ak)) THEN
702	CALL ropp_fm_state2state_ecmwf(x)
703	ELSE
704	CALL ropp_fm_state2state_meto(x)
705	END IF
706
707	IF (x%new_ref_op) THEN
708	CALL ropp_fm_refrac_1d_new(x, y)
709	ELSE
710	CALL ropp_fm_refrac_1d(x, y)
711	END IF
712
713	delta_x = (x%state - bg%state)
714	delta_y = (y%refrac - obs%refrac) * y%weights
715
716	J = 0.5_wp * DOT_PRODUCT(delta_y, MATMUL(Om1, delta_y)) + &
717	0.5_wp * DOT_PRODUCT(delta_x, MATMUL(Bm1, delta_x))
718
719	! 2.8 Levenberg-Marquardt update and convergence check
720	! ----------------------------------------------------
721
722	IF ( J >= (J_last + config%conv_check_max_delta_state) ) THEN ! Keep previous state vector, increase lambda, repeat
723
724	! 2.8.1 Keep old state, increase lambda, bail out if it's too big
725
726	x = x_old
727	lambda = 10.0_wp * lambda
728
729	IF (ABS(lambda) > lambda_max) THEN
730
731	WRITE(ch_str, '(g15.5)') lambda_max
732	ch_str = ADJUSTL(ch_str)
733	CALL message(msg_cont, '')
734	CALL message(msg_info, &
735	'Levenberg-Marquardt parameter exceeded max (=' // TRIM(ch_str) // &
736	'); convergence is unlikely.\n ' // &
737	'Check convergence parameters; conv_check_max_delta_J might ' // &
738	'be too demanding.\n')
739	EXIT
740	END IF
741
742	! 2.8.2 Logging
743
744	!TODO: Add a configuration option for lambda update outputs
745	IF (config%minropp%impres == 0) THEN
746
747	WRITE(it_str, '(i4)') n_iter
748	WRITE(ch_str, '(g15.5)') J
749	ch_str = ADJUSTL(ch_str)
750	IF (n_iter > 0) THEN
751	WRITE(co_str, '(g15.5)') lambda
752	co_str = ADJUSTL(co_str)
753	ELSE
754	co_str = ' - '
755	END IF
756
757	CALL message(msg_cont, &
758	' n_iter = ' // it_str // ' J = ' // ch_str // &
759	' lambda -> ' // co_str)
760
761	END IF
762
763	! 2.8.3 Next iteration
764
765	CYCLE
766
767	ELSE ! Decrease lambda, check convergence, continue iteration
768
769	! 2.8.4 Decrease lambda
770
771	lambda = 0.1_wp *lambda
772
773	! 2.8.5 Keep track of current state
774
775	n_iter = n_iter + 1
776
777	IF (config % conv_check_apply) THEN
778
779	i_pointer = MOD(i_pointer + 1, config % conv_check_n_previous)
780	IF (i_pointer == 0) i_pointer = config % conv_check_n_previous
781
782	delta_J(i_pointer) = J_last - J
783	delta_state(i_pointer) = MAXVAL(ABS(state_last - x%state)/state_sigma)
784
785	state_last = x%state
786
787	J_last = J
788
789	END IF
790
791	IF (J < J_min) THEN
792	J_min = J
793	END IF
794
795	! 2.8.6 Check for convergence
796
797	IF (config % conv_check_apply) THEN
798
799	IF (config%minropp%impres == 0) THEN
800	WRITE(it_str, '(i4)') n_iter
801	WRITE(ch_str, '(g15.5)') J
802	ch_str = ADJUSTL(ch_str)
803	IF (n_iter > 0) THEN
804	WRITE(co_str, '(g15.5)') delta_state(i_pointer)
805	co_str = ADJUSTL(co_str)
806	ELSE
807	co_str = ' - '
808	END IF
809
810	CALL message(msg_cont, &
811	' n_iter = ' // it_str // ' J = ' // ch_str // &
812	' max(relative change in state) = ' // co_str)
813	END IF
814
815	IF (MAXVAL(delta_state) < config%conv_check_max_delta_state &
816	.AND. n_iter > config % conv_check_n_previous) THEN
817
818	WRITE(ch_str, '(g15.5)') config%conv_check_max_delta_state
819	ch_str = ADJUSTL(ch_str)
820	WRITE(it_str, '(i2)') config%conv_check_n_previous
821	it_str = ADJUSTL(it_str)
822	CALL message(msg_cont, '')
823	CALL message(msg_info, &
824	'Convergence assumed to be achieved as the state vector did ' // &
825	'not change by more\n ' // 'than ' // TRIM(ch_str) // ' ' // &
826	'relative to the assumed background errors for the last ' // &
827	TRIM(it_str) // ' iterations.\n')
828	EXIT
829
830	ELSE IF (MAXVAL(ABS(delta_J)) < config%conv_check_max_delta_J &
831	.AND. n_iter > config % conv_check_n_previous) THEN
832
833	WRITE(ch_str, '(g15.5)') config%conv_check_max_delta_J
834	ch_str = ADJUSTL(ch_str)
835	WRITE(it_str, '(i2)') config%conv_check_n_previous
836	it_str = ADJUSTL(it_str)
837	CALL message(msg_cont, '')
838	CALL message(msg_info, &
839	'Convergence assumed to be achieved as the cost function did ' // &
840	'not change by more\n ' // 'than ' // TRIM(ch_str) // &
841	' for the last ' // TRIM(it_str) // ' iterations.\n')
842	EXIT
843
844	ELSE IF (ABS(lambda) > lambda_max) THEN
845
846	WRITE(ch_str, '(g15.5)') lambda_max
847	ch_str = ADJUSTL(ch_str)
848	CALL message(msg_cont, '')
849	CALL message(msg_info, &
850	'Levenberg-Marquardt parameter exceeded max (=' // TRIM(ch_str) // &
851	'); convergence is unlikely.\n ' // &
852	'Check convergence parameters; conv_check_max_delta_J might ' // &
853	'be too demanding.\n')
854	EXIT
855
856	END IF
857
858	END IF ! Convergence check
859
860	END IF ! Levenberg-Marquardt update
861
862	END DO ! End main iteration loop
863
864	IF (n_iter < n_iter_max) THEN
865	WRITE(it_str, '(i4)') n_iter ; it_str = ADJUSTL(it_str)
866	WRITE(gr_str, '(i4)') n_grad ; gr_str = ADJUSTL(gr_str)
867	CALL message(msg_info, 'Finished after ' // TRIM(it_str) // ' iterations (' // &
868	TRIM(gr_str) // ' forward model / gradient evaluations).\n')
869	ELSE
870	WRITE(it_str, '(i4)') n_iter_max ; it_str = ADJUSTL(it_str)
871	CALL message(msg_warn, &
872	'Iteration ended after n_iter_max = ' // TRIM(it_str) // &
873	' iterations without achieving convergence.\n')
874	END IF
875
876	! 2.9 Copy solution back to state
877	! -------------------------------
878
879	state = x
880
881	! 2.10 Diagnostic data
882	! --------------------
883
884	diag%J = J_min
885
886	IF (COUNT(obs%weights > 0.0_wp) > 0) THEN
887	diag%J_scaled = 2.0_wp * J_min / REAL(COUNT(obs%weights > 0.0_wp), wp)
888	ENDIF
889
890	diag%n_iter = n_iter
891
892	ALLOCATE (diag%J_bgr(SIZE(state%state)))
893	delta_x = state%state - bg%state
894	diag%J_bgr = 0.5_wp * delta_x * matrix_solve(bg%cov, delta_x)
895
896	! 2.11 Clean up
897	! -------------
898
899	DEALLOCATE(K)
900	DEALLOCATE(delta_x)
901	DEALLOCATE(delta_y)
902	DEALLOCATE(dJ_dx)
903	DEALLOCATE(diag_d2J)
904	DEALLOCATE(d2J_dx2)
905	DEALLOCATE(KO)
906	DEALLOCATE(Bm1)
907	DEALLOCATE(Om1)
908
909	CALL message_set_routine(routine)
910
911	END SUBROUTINE ropp_1dvar_levmarq_refrac

Download in other formats:

Original format