Source
xxxxxxxxxx
1455
1455
}
1456
1456
}
1457
1457
for (auto it=aggregateTime.begin(); it!=aggregateTime.end(); ++it) {
1458
1458
Int a = it->first;
1459
1459
it->second /= aggregateWeight[a];
1460
1460
}
1461
1461
1462
1462
}
1463
1463
1464
1464
1465
+
void
1466
+
print_gsl_vector(gsl_vector *v)
1467
+
{
1468
+
const size_t n = v->size;
1469
+
for (int i=0; i!=n; i++) {
1470
+
cerr << gsl_vector_get(v, i) << " ";
1471
+
if (i>0 && (i % 4)==0) cerr << endl;
1472
+
}
1473
+
cerr << endl;
1474
+
}
1465
1475
1466
1476
void
1467
-
least_squares_driver(SDBList& sdbs, Matrix<Float>& casa_param, Matrix<Bool>& casa_flags, Matrix<Float>& casa_snr,
1468
-
Int refant, const std::map< Int, std::set<Int> >& activeAntennas, LogIO& logSink) {
1477
+
print_max_gsl3(gsl_vector *v)
1478
+
{
1479
+
double phi_max = 0.0;
1480
+
double del_max = 0.0;
1481
+
double rat_max = 0.0;
1482
+
1483
+
const size_t n = v->size;
1484
+
for (int i=0; i!=n/3; i++) {
1485
+
if (fabs(gsl_vector_get(v, 3*i+0)) > fabs(phi_max)) phi_max = gsl_vector_get(v, 3*i+0);
1486
+
if (fabs(gsl_vector_get(v, 3*i+1)) > fabs(del_max)) del_max = gsl_vector_get(v, 3*i+1);
1487
+
if (fabs(gsl_vector_get(v, 3*i+2)) > fabs(rat_max)) rat_max = gsl_vector_get(v, 3*i+2);
1488
+
}
1489
+
cerr << "phi_max " << phi_max << " del_max " << del_max << " rat_max " << rat_max << endl;
1490
+
}
1491
+
1492
+
1493
+
1494
+
/*
1495
+
gsl-2.4/multilarge_nlinear/fdf.c defines gsl_multilarge_nlinear_driver,
1496
+
which I have butchered for my purposes here into
1497
+
not_gsl_multilarge_nlinear_driver(). We still iterate the nonlinear
1498
+
least squares solver until completion, but we adopt a convergence
1499
+
criterion copied from AIPS.
1500
+
1501
+
Inputs: maxiter - maximum iterations to allow
1502
+
w - workspace
1503
+
1504
+
Additionally I've removed the info parameter, and I may yet regret it.
1505
+
1506
+
Originally:
1507
+
info - (output) info flag on why iteration terminated
1508
+
1 = stopped due to small step size ||dx|
1509
+
2 = stopped due to small gradient
1510
+
3 = stopped due to small change in f
1511
+
GSL_ETOLX = ||dx|| has converged to within machine
1512
+
precision (and xtol is too small)
1513
+
GSL_ETOLG = ||g||_inf is smaller than machine
1514
+
precision (gtol is too small)
1515
+
GSL_ETOLF = change in ||f|| is smaller than machine
1516
+
precision (ftol is too small)
1517
+
1518
+
Return:
1519
+
GSL_SUCCESS if converged
1520
+
GSL_MAXITER if maxiter exceeded without converging
1521
+
GSL_ENOPROG if no accepted step found on first iteration
1522
+
*/
1523
+
1524
+
int
1525
+
least_squares_inner_driver (const size_t maxiter,
1526
+
gsl_multilarge_nlinear_workspace * w)
1527
+
{
1528
+
int status;
1529
+
size_t iter = 0;
1530
+
/* call user callback function prior to any iterations
1531
+
* with initial system state */
1532
+
Double s;
1533
+
Double last_s = 1.0e30;
1534
+
Bool converged = false;
1535
+
do {
1536
+
status = gsl_multilarge_nlinear_iterate (w);
1537
+
/*
1538
+
* If the solver reports no progress on the first iteration,
1539
+
* then it didn't find a single step to reduce the
1540
+
* cost function and more iterations won't help so return.
1541
+
*
1542
+
* If we get a no progress flag on subsequent iterations,
1543
+
* it means we did find a good step in a previous iteration,
1544
+
* so continue iterating since the solver has now reset
1545
+
* mu to its initial value.
1546
+
*/
1547
+
if (status == GSL_ENOPROG && iter == 0) {
1548
+
return GSL_EMAXITER;
1549
+
}
1550
+
1551
+
Double fnorm = gsl_blas_dnrm2(w->f);
1552
+
s = 0.5 * fnorm * fnorm;
1553
+
if ((iter > 0) && DEVDEBUG) {
1554
+
// cerr << "Parameters: " << endl;
1555
+
// print_gsl_vector(w->x);
1556
+
cerr << "Iter: " << iter << " ";
1557
+
print_max_gsl3(w->dx);
1558
+
cerr << "1 - s/last_s=" << 1 - s/last_s << endl;
1559
+
}
1560
+
++iter;
1561
+
if ((1 - s/last_s < 5e-6) && (iter > 1)) converged = true;
1562
+
last_s = s;
1563
+
/* old test for convergence:
1564
+
status = not_gsl_multilarge_nlinear_test(xtol, gtol, ftol, info, w); */
1565
+
} while (!converged && iter < maxiter);
1566
+
/*
1567
+
* the following error codes mean that the solution has converged
1568
+
* to within machine precision, so record the error code in info
1569
+
* and return success
1570
+
*/
1571
+
if (status == GSL_ETOLF || status == GSL_ETOLX || status == GSL_ETOLG)
1572
+
{
1573
+
status = GSL_SUCCESS;
1574
+
}
1575
+
/* check if max iterations reached */
1576
+
if (iter >= maxiter && status != GSL_SUCCESS)
1577
+
status = GSL_EMAXITER;
1578
+
return status;
1579
+
} /* gsl_multilarge_nlinear_driver() */
1580
+
1581
+
1582
+
1583
+
1584
+
void
1585
+
least_squares_driver(SDBList& sdbs, Matrix<Float>& casa_param, Matrix<Bool>& casa_flags, Matrix<Float>& casa_snr, Int refant,
1586
+
const std::map< Int, std::set<Int> >& activeAntennas, LogIO& logSink) {
1469
1587
// The variable casa_param is the Casa calibration framework's RParam matrix; we transcribe our results into it only at the end.
1470
1588
// n below is number of variables,
1471
1589
// p is number of parameters
1472
1590
1473
1591
AuxParamBundle bundle(sdbs, refant, activeAntennas);
1592
+
1474
1593
for (size_t icor=0; icor != bundle.get_num_corrs(); icor++) {
1475
1594
bundle.set_active_corr(icor);
1476
1595
if (bundle.get_num_antennas() == 0) {
1477
1596
logSink << "No antennas for correlation " << icor << "; not running least-squares solver." << LogIO::POST;
1478
1597
continue;
1479
1598
}
1480
1599
if (bundle.get_num_antennas() == 1) {
1481
1600
logSink << "No baselines for correlation " << icor << "; not running least-squares solver." << LogIO::POST;
1482
1601
continue;
1483
1602
}
1535
1654
gsl_vector *gp_orig = gsl_vector_alloc(p);
1536
1655
// Keep a copy of original parameters
1537
1656
gsl_vector_memcpy (gp_orig, gp);
1538
1657
// initialise workspace
1539
1658
gsl_multilarge_nlinear_init(gp, &f, w);
1540
1659
1541
1660
// compute initial residual norm */
1542
1661
gsl_vector *res_f = gsl_multilarge_nlinear_residual(w);
1543
1662
1544
1663
int info;
1545
-
int status = gsl_multilarge_nlinear_driver(max_iter, param_tol, gtol, ftol,
1546
-
NULL, NULL, &info, w);
1664
+
int status = least_squares_inner_driver(max_iter, w);
1547
1665
double chi1 = gsl_blas_dnrm2(res_f);
1548
1666
1549
1667
gsl_vector_sub(gp_orig, w->x);
1550
1668
gsl_vector *diff = gp_orig;
1551
-
// double diffsize =
1552
-
gsl_blas_dnrm2(diff);
1669
+
double diffsize = gsl_blas_dnrm2(diff);
1553
1670
1554
1671
gsl_vector *res = gsl_multilarge_nlinear_position(w);
1555
1672
1556
1673
// We transcribe values back from gsl_vector to the param matrix
1557
1674
1558
1675
gsl_matrix *hess = gsl_matrix_alloc(p,p);
1559
1676
gsl_vector *snr_vector = gsl_vector_alloc(p);
1560
1677
gsl_matrix_set_zero(hess);
1561
1678
gsl_vector_set_zero(snr_vector);
1562
1679
expb_hess(gp, &bundle, hess, chi1*chi1, snr_vector, logSink);