Context navigation

Back to Ticket #646

Ticket #646: ropp_1dvar_levmarq_tr1.2.f90

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

Download in other formats:

Original format