Carry out optimal estimation retrieval.
462 {
463
464 static int ipa[
N], iqa[
N];
465
466 FILE *out;
467
468 char filename[2 *
LEN];
469
470 double chisq, disq = 0, lmpar = 0.001;
471
472 int it = 0;
473
474
475
476
477
478
479 const size_t m =
obs2y(ctl, obs_meas, NULL, NULL, NULL);
480 const size_t n =
atm2x(ctl, atm_apr, NULL, iqa, ipa);
481 if (m == 0 || n == 0)
482 ERRMSG(
"Check problem definition!");
483
484
485 LOG(1,
"Problem size: m= %d / n= %d "
486 "(alloc= %.4g MB / stat= %.4g MB)",
487 (int) m, (int) n,
488 (double) (3 * m * n + 4 * n * n + 8 * m +
489 8 * n) * sizeof(double) / 1024. / 1024.,
490 (
double) (5 *
sizeof(
atm_t) + 3 *
sizeof(
obs_t)
491 + 2 *
N *
sizeof(
int)) / 1024. / 1024.);
492
493
494 gsl_matrix *a = gsl_matrix_alloc(n, n);
495 gsl_matrix *cov = gsl_matrix_alloc(n, n);
496 gsl_matrix *k_i = gsl_matrix_alloc(m, n);
497 gsl_matrix *s_a_inv = gsl_matrix_alloc(n, n);
498
499 gsl_vector *b = gsl_vector_alloc(n);
500 gsl_vector *dx = gsl_vector_alloc(n);
501 gsl_vector *dy = gsl_vector_alloc(m);
502 gsl_vector *sig_eps_inv = gsl_vector_alloc(m);
503 gsl_vector *sig_formod = gsl_vector_alloc(m);
504 gsl_vector *sig_noise = gsl_vector_alloc(m);
505 gsl_vector *x_a = gsl_vector_alloc(n);
506 gsl_vector *x_i = gsl_vector_alloc(n);
507 gsl_vector *x_step = gsl_vector_alloc(n);
508 gsl_vector *y_aux = gsl_vector_alloc(m);
509 gsl_vector *y_i = gsl_vector_alloc(m);
510 gsl_vector *y_m = gsl_vector_alloc(m);
511
512
515 formod(ctl, tbl, atm_i, obs_i);
516
517
518 atm2x(ctl, atm_apr, x_a, NULL, NULL);
519 atm2x(ctl, atm_i, x_i, NULL, NULL);
520 obs2y(ctl, obs_meas, y_m, NULL, NULL);
521 obs2y(ctl, obs_i, y_i, NULL, NULL);
522
523
526 atm_i, obs_i, "x", "x", "r");
528
529
530 set_cov_meas(ret, ctl, obs_meas, sig_noise, sig_formod, sig_eps_inv);
531
532
533 sprintf(filename,
"%s/costs.tab", ret->
dir);
534 if (!(out = fopen(filename, "w")))
535 ERRMSG(
"Cannot create cost function file!");
536
537
538 fprintf(out,
539 "# $1 = iteration number\n"
540 "# $2 = normalized cost function\n"
541 "# $3 = number of measurements\n"
542 "# $4 = number of state vector elements\n\n");
543
544
545 gsl_vector_memcpy(dx, x_i);
546 gsl_vector_sub(dx, x_a);
547 gsl_vector_memcpy(dy, y_m);
548 gsl_vector_sub(dy, y_i);
549
550
552
553
554 LOG(1,
"it= %d / chi^2/m= %g", it, chisq);
555
556
557 fprintf(out, "%d %g %d %d\n", it, chisq, (int) m, (int) n);
558
559
560 kernel(ctl, tbl, atm_i, obs_i, k_i);
561
562
563
564
565
566
568
569
570 double chisq_old = chisq;
571
572
574 kernel(ctl, tbl, atm_i, obs_i, k_i);
575
576
579
580
581 for (size_t i = 0; i < m; i++)
582 gsl_vector_set(y_aux, i, gsl_vector_get(dy, i)
583 *
POW2(gsl_vector_get(sig_eps_inv, i)));
584 gsl_blas_dgemv(CblasTrans, 1.0, k_i, y_aux, 0.0, b);
585 gsl_blas_dgemv(CblasNoTrans, -1.0, s_a_inv, dx, 1.0, b);
586
587
588 for (int it2 = 0; it2 < 20; it2++) {
589
590
591 gsl_matrix_memcpy(a, s_a_inv);
592 gsl_matrix_scale(a, 1 + lmpar);
593 gsl_matrix_add(a, cov);
594
595
596 gsl_linalg_cholesky_decomp(a);
597 gsl_linalg_cholesky_solve(a, b, x_step);
598
599
600 gsl_vector_add(x_i, x_step);
603 x2atm(ctl, x_i, atm_i);
604
605
606 for (
int ip = 0; ip < atm_i->
np; ip++) {
607 atm_i->
p[ip] =
MIN(
MAX(atm_i->
p[ip], 5e-7), 5e4);
608 atm_i->
t[ip] =
MIN(
MAX(atm_i->
t[ip], 100), 400);
609 for (
int ig = 0; ig < ctl->
ng; ig++)
610 atm_i->
q[ig][ip] =
MIN(
MAX(atm_i->
q[ig][ip], 0), 1);
611 for (
int iw = 0; iw < ctl->
nw; iw++)
612 atm_i->
k[iw][ip] =
MAX(atm_i->
k[iw][ip], 0);
613 }
616 for (
int icl = 0; icl < ctl->
ncl; icl++)
617 atm_i->
clk[icl] =
MAX(atm_i->
clk[icl], 0);
621 for (
int isf = 0; isf < ctl->
nsf; isf++)
623
624
625 formod(ctl, tbl, atm_i, obs_i);
626 obs2y(ctl, obs_i, y_i, NULL, NULL);
627
628
629 gsl_vector_memcpy(dx, x_i);
630 gsl_vector_sub(dx, x_a);
631 gsl_vector_memcpy(dy, y_m);
632 gsl_vector_sub(dy, y_i);
633
634
636
637
638 if (chisq > chisq_old) {
639 lmpar *= 10;
640 gsl_vector_sub(x_i, x_step);
641 } else {
642 lmpar /= 10;
643 break;
644 }
645 }
646
647
648 LOG(1,
"it= %d / chi^2/m= %g", it, chisq);
649
650
651 fprintf(out, "%d %g %d %d\n", it, chisq, (int) m, (int) n);
652
653
654 gsl_blas_ddot(x_step, b, &disq);
655 disq /= (double) n;
656
657
659 break;
660 }
661
662
663 fclose(out);
664
665
669 atm_i, obs_i, "y", "x", "r");
670
671
672
673
674
675
677
678
679 gsl_matrix *auxnm = gsl_matrix_alloc(n, m);
680 gsl_matrix *corr = gsl_matrix_alloc(n, n);
681 gsl_matrix *gain = gsl_matrix_alloc(n, m);
682
683
684
686 gsl_matrix_add(cov, s_a_inv);
687
688
691 atm_i, obs_i, "x", "x", "r");
693
694
695 for (size_t i = 0; i < n; i++)
696 for (size_t j = 0; j < n; j++)
697 gsl_matrix_set(corr, i, j, gsl_matrix_get(cov, i, j)
698 / sqrt(gsl_matrix_get(cov, i, i))
699 / sqrt(gsl_matrix_get(cov, j, j)));
701 atm_i, obs_i, "x", "x", "r");
702
703
704
705 for (size_t i = 0; i < n; i++)
706 for (size_t j = 0; j < m; j++)
707 gsl_matrix_set(auxnm, i, j, gsl_matrix_get(k_i, j, i)
708 *
POW2(gsl_vector_get(sig_eps_inv, j)));
709 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, cov, auxnm, 0.0, gain);
711 atm_i, obs_i, "x", "y", "c");
712
713
716
717
720
721
722
723 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, gain, k_i, 0.0, a);
725 atm_i, obs_i, "x", "x", "r");
726
727
729
730
731 gsl_matrix_free(auxnm);
732 gsl_matrix_free(corr);
733 gsl_matrix_free(gain);
734 }
735
736
737
738
739
740 gsl_matrix_free(a);
741 gsl_matrix_free(cov);
742 gsl_matrix_free(k_i);
743 gsl_matrix_free(s_a_inv);
744
745 gsl_vector_free(b);
746 gsl_vector_free(dx);
747 gsl_vector_free(dy);
748 gsl_vector_free(sig_eps_inv);
749 gsl_vector_free(sig_formod);
750 gsl_vector_free(sig_noise);
751 gsl_vector_free(x_a);
752 gsl_vector_free(x_i);
753 gsl_vector_free(x_step);
754 gsl_vector_free(y_aux);
755 gsl_vector_free(y_i);
756 gsl_vector_free(y_m);
757}
void x2atm(const ctl_t *ctl, const gsl_vector *x, atm_t *atm)
Decompose parameter vector or state vector.
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.
void kernel(const ctl_t *ctl, const tbl_t *tbl, atm_t *atm, obs_t *obs, gsl_matrix *k)
Compute Jacobians.
void formod(const ctl_t *ctl, const tbl_t *tbl, atm_t *atm, obs_t *obs)
Determine ray paths and compute radiative transfer.
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.
#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.