Carry out optimal estimation retrieval.
454 {
455
456 static int ipa[
N], iqa[
N];
457
458 FILE *out;
459
460 char filename[2 *
LEN];
461
462 double chisq, disq = 0, lmpar = 0.001;
463
464 int it = 0;
465
466
467
468
469
470
471 const size_t m =
obs2y(ctl, obs_meas, NULL, NULL, NULL);
472 const size_t n =
atm2x(ctl, atm_apr, NULL, iqa, ipa);
473 if (m == 0 || n == 0)
474 ERRMSG(
"Check problem definition!");
475
476
477 LOG(1,
"Problem size: m= %d / n= %d "
478 "(alloc= %.4g MB / stat= %.4g MB)",
479 (int) m, (int) n,
480 (double) (3 * m * n + 4 * n * n + 8 * m +
481 8 * n) * sizeof(double) / 1024. / 1024.,
482 (
double) (5 *
sizeof(
atm_t) + 3 *
sizeof(
obs_t)
483 + 2 *
N *
sizeof(
int)) / 1024. / 1024.);
484
485
486 gsl_matrix *a = gsl_matrix_alloc(n, n);
487 gsl_matrix *cov = gsl_matrix_alloc(n, n);
488 gsl_matrix *k_i = gsl_matrix_alloc(m, n);
489 gsl_matrix *s_a_inv = gsl_matrix_alloc(n, n);
490
491 gsl_vector *b = gsl_vector_alloc(n);
492 gsl_vector *dx = gsl_vector_alloc(n);
493 gsl_vector *dy = gsl_vector_alloc(m);
494 gsl_vector *sig_eps_inv = gsl_vector_alloc(m);
495 gsl_vector *sig_formod = gsl_vector_alloc(m);
496 gsl_vector *sig_noise = gsl_vector_alloc(m);
497 gsl_vector *x_a = gsl_vector_alloc(n);
498 gsl_vector *x_i = gsl_vector_alloc(n);
499 gsl_vector *x_step = gsl_vector_alloc(n);
500 gsl_vector *y_aux = gsl_vector_alloc(m);
501 gsl_vector *y_i = gsl_vector_alloc(m);
502 gsl_vector *y_m = gsl_vector_alloc(m);
503
504
507 formod(ctl, atm_i, obs_i);
508
509
510 atm2x(ctl, atm_apr, x_a, NULL, NULL);
511 atm2x(ctl, atm_i, x_i, NULL, NULL);
512 obs2y(ctl, obs_meas, y_m, NULL, NULL);
513 obs2y(ctl, obs_i, y_i, NULL, NULL);
514
515
518 atm_i, obs_i, "x", "x", "r");
520
521
522 set_cov_meas(ret, ctl, obs_meas, sig_noise, sig_formod, sig_eps_inv);
523
524
525 sprintf(filename,
"%s/costs.tab", ret->
dir);
526 if (!(out = fopen(filename, "w")))
527 ERRMSG(
"Cannot create cost function file!");
528
529
530 fprintf(out,
531 "# $1 = iteration number\n"
532 "# $2 = normalized cost function\n"
533 "# $3 = number of measurements\n"
534 "# $4 = number of state vector elements\n\n");
535
536
537 gsl_vector_memcpy(dx, x_i);
538 gsl_vector_sub(dx, x_a);
539 gsl_vector_memcpy(dy, y_m);
540 gsl_vector_sub(dy, y_i);
541
542
544
545
546 LOG(1,
"it= %d / chi^2/m= %g", it, chisq);
547
548
549 fprintf(out, "%d %g %d %d\n", it, chisq, (int) m, (int) n);
550
551
552 kernel(ctl, atm_i, obs_i, k_i);
553
554
555
556
557
558
560
561
562 double chisq_old = chisq;
563
564
566 kernel(ctl, atm_i, obs_i, k_i);
567
568
571
572
573 for (size_t i = 0; i < m; i++)
574 gsl_vector_set(y_aux, i, gsl_vector_get(dy, i)
575 *
POW2(gsl_vector_get(sig_eps_inv, i)));
576 gsl_blas_dgemv(CblasTrans, 1.0, k_i, y_aux, 0.0, b);
577 gsl_blas_dgemv(CblasNoTrans, -1.0, s_a_inv, dx, 1.0, b);
578
579
580 for (int it2 = 0; it2 < 20; it2++) {
581
582
583 gsl_matrix_memcpy(a, s_a_inv);
584 gsl_matrix_scale(a, 1 + lmpar);
585 gsl_matrix_add(a, cov);
586
587
588 gsl_linalg_cholesky_decomp(a);
589 gsl_linalg_cholesky_solve(a, b, x_step);
590
591
592 gsl_vector_add(x_i, x_step);
595 x2atm(ctl, x_i, atm_i);
596
597
598 for (
int ip = 0; ip < atm_i->
np; ip++) {
599 atm_i->
p[ip] =
MIN(
MAX(atm_i->
p[ip], 5e-7), 5e4);
600 atm_i->
t[ip] =
MIN(
MAX(atm_i->
t[ip], 100), 400);
601 for (
int ig = 0; ig < ctl->
ng; ig++)
602 atm_i->
q[ig][ip] =
MIN(
MAX(atm_i->
q[ig][ip], 0), 1);
603 for (
int iw = 0; iw < ctl->
nw; iw++)
604 atm_i->
k[iw][ip] =
MAX(atm_i->
k[iw][ip], 0);
605 }
608 for (
int icl = 0; icl < ctl->
ncl; icl++)
609 atm_i->
clk[icl] =
MAX(atm_i->
clk[icl], 0);
613 for (
int isf = 0; isf < ctl->
nsf; isf++)
615
616
617 formod(ctl, atm_i, obs_i);
618 obs2y(ctl, obs_i, y_i, NULL, NULL);
619
620
621 gsl_vector_memcpy(dx, x_i);
622 gsl_vector_sub(dx, x_a);
623 gsl_vector_memcpy(dy, y_m);
624 gsl_vector_sub(dy, y_i);
625
626
628
629
630 if (chisq > chisq_old) {
631 lmpar *= 10;
632 gsl_vector_sub(x_i, x_step);
633 } else {
634 lmpar /= 10;
635 break;
636 }
637 }
638
639
640 LOG(1,
"it= %d / chi^2/m= %g", it, chisq);
641
642
643 fprintf(out, "%d %g %d %d\n", it, chisq, (int) m, (int) n);
644
645
646 gsl_blas_ddot(x_step, b, &disq);
647 disq /= (double) n;
648
649
651 break;
652 }
653
654
655 fclose(out);
656
657
661 atm_i, obs_i, "y", "x", "r");
662
663
664
665
666
667
669
670
671 gsl_matrix *auxnm = gsl_matrix_alloc(n, m);
672 gsl_matrix *corr = gsl_matrix_alloc(n, n);
673 gsl_matrix *gain = gsl_matrix_alloc(n, m);
674
675
676
678 gsl_matrix_add(cov, s_a_inv);
679
680
683 atm_i, obs_i, "x", "x", "r");
685
686
687 for (size_t i = 0; i < n; i++)
688 for (size_t j = 0; j < n; j++)
689 gsl_matrix_set(corr, i, j, gsl_matrix_get(cov, i, j)
690 / sqrt(gsl_matrix_get(cov, i, i))
691 / sqrt(gsl_matrix_get(cov, j, j)));
693 atm_i, obs_i, "x", "x", "r");
694
695
696
697 for (size_t i = 0; i < n; i++)
698 for (size_t j = 0; j < m; j++)
699 gsl_matrix_set(auxnm, i, j, gsl_matrix_get(k_i, j, i)
700 *
POW2(gsl_vector_get(sig_eps_inv, j)));
701 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, cov, auxnm, 0.0, gain);
703 atm_i, obs_i, "x", "y", "c");
704
705
708
709
712
713
714
715 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, gain, k_i, 0.0, a);
717 atm_i, obs_i, "x", "x", "r");
718
719
721
722
723 gsl_matrix_free(auxnm);
724 gsl_matrix_free(corr);
725 gsl_matrix_free(gain);
726 }
727
728
729
730
731
732 gsl_matrix_free(a);
733 gsl_matrix_free(cov);
734 gsl_matrix_free(k_i);
735 gsl_matrix_free(s_a_inv);
736
737 gsl_vector_free(b);
738 gsl_vector_free(dx);
739 gsl_vector_free(dy);
740 gsl_vector_free(sig_eps_inv);
741 gsl_vector_free(sig_formod);
742 gsl_vector_free(sig_noise);
743 gsl_vector_free(x_a);
744 gsl_vector_free(x_i);
745 gsl_vector_free(x_step);
746 gsl_vector_free(y_aux);
747 gsl_vector_free(y_i);
748 gsl_vector_free(y_m);
749}
void x2atm(const ctl_t *ctl, const gsl_vector *x, atm_t *atm)
Decompose parameter vector or state vector.
void formod(const ctl_t *ctl, atm_t *atm, obs_t *obs)
Determine ray paths and compute radiative transfer.
void copy_obs(const ctl_t *ctl, obs_t *obs_dest, const obs_t *obs_src, const int init)
Copy and initialize observation data.
void write_obs(const char *dirname, const char *filename, const ctl_t *ctl, const obs_t *obs)
Write observation data.
size_t obs2y(const ctl_t *ctl, const obs_t *obs, gsl_vector *y, int *ida, int *ira)
Compose measurement vector.
size_t atm2x(const ctl_t *ctl, const atm_t *atm, gsl_vector *x, int *iqa, int *ipa)
Compose state vector or parameter vector.
void write_matrix(const char *dirname, const char *filename, const ctl_t *ctl, const gsl_matrix *matrix, const atm_t *atm, const obs_t *obs, const char *rowspace, const char *colspace, const char *sort)
Write matrix.
void kernel(ctl_t *ctl, atm_t *atm, obs_t *obs, gsl_matrix *k)
Compute Jacobians.
#define LEN
Maximum length of ASCII data lines.
#define MIN(a, b)
Macro to determine the minimum of two values.
#define ERRMSG(...)
Print error message and quit program.
#define LOG(level,...)
Print log message.
#define MAX(a, b)
Macro to determine the maximum of two values.
void matrix_product(gsl_matrix *a, gsl_vector *b, int transpose, gsl_matrix *c)
Compute matrix product A^TBA or ABA^T for diagonal matrix B.
void set_cov_apr(ret_t *ret, ctl_t *ctl, atm_t *atm, int *iqa, int *ipa, gsl_matrix *s_a)
Set a priori covariance.
void analyze_avk(ret_t *ret, ctl_t *ctl, atm_t *atm, int *iqa, int *ipa, gsl_matrix *avk)
Compute information content and resolution.
void set_cov_meas(ret_t *ret, ctl_t *ctl, obs_t *obs, gsl_vector *sig_noise, gsl_vector *sig_formod, gsl_vector *sig_eps_inv)
Set measurement errors.
void write_stddev(const char *quantity, ret_t *ret, ctl_t *ctl, atm_t *atm, gsl_matrix *s)
Write retrieval error to file.
void matrix_invert(gsl_matrix *a)
Invert symmetric matrix.
double cost_function(gsl_vector *dx, gsl_vector *dy, gsl_matrix *s_a_inv, gsl_vector *sig_eps_inv)
Compute cost function.
int np
Number of data points.
Observation geometry and radiance data.
int err_ana
Carry out error analysis (0=no, 1=yes).
double conv_dmin
Minimum normalized step size in state space.
int kernel_recomp
Re-computation of kernel matrix (number of iterations).
int conv_itmax
Maximum number of iterations.