solve.cc

Go to the documentation of this file.

#include <algorithm>
#include <ctime>
#include <thread>

#include <gsl/gsl_multifit_nlin.h>

#include "Graybody.h"
#include "PACS.h"
#include "SPIRE.h"

namespace mu = mutils;

namespace manticore {

const double  mu_g = (2.75*m_H),  
             imu_g = 1.0/mu_g;    

struct GrayData {
  GrayData(Graybody *g, std::vector<mapDataType> *o,
           size_t bands, size_t param) :
    gray(g), outMap(o), nBands(bands), nParam(param)
  {
    Fobs.resize(nBands);
    Ferr.resize(nBands);
    band.resize(nBands);
  }

  Graybody *gray;              

  std::vector<long> axes;      

  std::vector<double> Fobs,    
                      Ferr,    
                      dnu0,    
                      fixed;   

  std::vector<const Detector*>
                      band;    

  std::vector<mapDataType>
                     *outMap;  


  size_t nBands,               
         nParam;               
};


int f_gray(const gsl_vector *x, void *p, gsl_vector *f)
{
  GrayData *gd = static_cast<GrayData *>(p);
  const auto &Fobs = gd->Fobs, &Ferr = gd->Ferr;
  const auto &gray = *gd->gray;

  size_t n = Ferr.size();
  double T = gsl_vector_get(x, 0), S = gsl_vector_get(x, 1);
  for (size_t i=0; i != n; i++) {
    auto Di = gd->band[i];
    double iFerr = 1.0/Ferr[i], F = gray.F(T, S, Di);
    gsl_vector_set(f, i, (F - Fobs[i])*iFerr);
  }

  return GSL_SUCCESS;
}

int df_gray(const gsl_vector *x, void *p, gsl_matrix *J)
{
  GrayData *gd = static_cast<GrayData *>(p);
  const auto &Ferr = gd->Ferr;
  const auto &gray = *gd->gray;

  size_t n = Ferr.size();
  double T = gsl_vector_get(x, 0), S = gsl_vector_get(x, 1);
  for (size_t i=0; i != n; i++) {
    auto Di = gd->band[i];
    double iFerr = 1.0/Ferr[i];
    gsl_matrix_set(J, i, 0, gray.dF_dT(T, S, Di)*iFerr);
    gsl_matrix_set(J, i, 1, gray.dF_dS(T, S, Di)*iFerr);
  }

  return GSL_SUCCESS;
}

int fdf_gray(const gsl_vector *x, void *p, gsl_vector *f, gsl_matrix *J)
{ return f_gray(x, p, f), df_gray(x, p, J); }


int f_gray2(const gsl_vector *x, void *p, gsl_vector *f)
{
  GrayData *gd = static_cast<GrayData *>(p);
  const auto &Fobs = gd->Fobs, &Ferr = gd->Ferr;
  const auto &gray = *gd->gray;

  size_t n = Ferr.size(), nParam = gd->nParam;
  double Sc = gsl_vector_get(x, 0), Tc = gsl_vector_get(x, 1),
         Sh = gsl_vector_get(x, 2),
         Th = (nParam == 4 ? gsl_vector_get(x, 3) : gd->fixed[0]);
  for (size_t i=0; i != n; i++) {
    auto Di = gd->band[i];
    double iFerr = 1.0/Ferr[i];
    gsl_vector_set(f, i, (gray.F(Tc, Sc, Th, Sh, Di) - Fobs[i])*iFerr);
  }

  return GSL_SUCCESS;
}

int df_gray2(const gsl_vector *x, void *p, gsl_matrix *J)
{
  GrayData *gd = static_cast<GrayData *>(p);
  const auto &Ferr = gd->Ferr;
  const auto &gray = *gd->gray;

  size_t n = Ferr.size(), nParam = gd->nParam;
  double Sc = gsl_vector_get(x, 0), Tc = gsl_vector_get(x, 1),
         Sh = gsl_vector_get(x, 2),
         Th = (nParam == 4 ? gsl_vector_get(x, 3) : gd->fixed[0]);
  for (size_t i=0; i != n; i++) {
    auto Di = gd->band[i];
    double iFerr = 1.0/Ferr[i];
    gsl_matrix_set(J, i, 0, gray.dF_dSc(Tc, Sc, Th, Sh, Di)*iFerr);
    gsl_matrix_set(J, i, 1, gray.dF_dTc(Tc, Sc, Th, Sh, Di)*iFerr);
    gsl_matrix_set(J, i, 2, gray.dF_dSh(Tc, Sc, Th, Sh, Di)*iFerr);
    if (nParam == 4) {
      gsl_matrix_set(J, i, 3, gray.dF_dTh(Tc, Sc, Th, Sh, Di)*iFerr);
    }
  }

  return GSL_SUCCESS;
}

int fdf_gray2(const gsl_vector *x, void *p, gsl_vector *f, gsl_matrix *J)
{ return f_gray2(x, p, f), df_gray2(x, p, J); }


std::string asctime()
{
  char buf[64];
  struct tm tmbuf;
  auto t = std::time(nullptr);
  std::string rtn = asctime_r(localtime_r(&t, &tmbuf), buf);
  rtn.pop_back();  // Remove newline.
  return rtn;
}


void makeModelMap(std::vector<mapDataType> &bandData,
                  std::vector<mapDataType> &bandError,
                  const GrayData &gdata)
{
  const auto nParam = gdata.nParam;
  const auto &gray = *gdata.gray;
  const int n = sqrt(1.0*modelMapLen);
  const double in1 = 1.0/(n-1);

  for (unsigned iMap=0; iMap != modelMapLen*modelMapLen; iMap++) {
    auto i = iMap/modelMapLen, j = iMap%modelMapLen;
    double
      Tc = 6.0*(1.0 + (j%n)*in1),           // Tc ~ [6, 12] K
      Sc = mu_g*3e21*pow(100.0, (j/n)*in1), // Nc(H2) ~ [3e21, 3e23] cm^-2
      Th = 15.0 + 3.0*(i%n)*in1,            // Th ~ [15, 18] K
      Sh = mu_g*3e19*pow(100.0, (i/n)*in1); // Nc(H2) ~ [3e19, 3e21] cm^-2
      MU_DEBUG( "makeModelMap",
                  "(%lu, %lu) Tc = %.2f, Sc = %.2e, Th = %.2f, Sh = %.2e",
                  i, j, Tc, Sc, Th, Sh);

    for (unsigned iBand=0; iBand != bandData.size(); iBand++) {
      auto &data = bandData[iBand], &err = bandError[iBand];
      auto Di = gdata.band[iBand];
      data[iMap] = (nParam == 2 ? gray.F(Th, Sh, Di) :
                                  gray.F(Tc, Sc, Th, Sh, Di)) /
                    gdata.dnu0[iBand];
       err[iMap] = 0.001*data[iMap];
      MU_DEBUG("makeModelMap", "  [%lu] %.4e", iBand, data[iMap]);
    }
  }
}


void solveStage1(const std::vector<mapDataType> &bandData,
                 const std::vector<mapDataType> &bandError,
                 const mapDataType &chi2Map,
                 const mu::CommandLine &cli,
                 const GrayData &gdata0, unsigned id, unsigned nThreads)
{
  const char *const fn = "solveStage1";
  size_t nBands = gdata0.nBands, nParam = gdata0.nParam,
         minBands = cli.getInteger('n', 3);
  auto &Tmap  = gdata0.outMap->at(0),
       &dTmap = gdata0.outMap->at(1),
       &Nmap  = gdata0.outMap->at(2),
       &dNmap = gdata0.outMap->at(3),
       &Cmap  = gdata0.outMap->at(4),
       &Imap  = gdata0.outMap->at(5);
  auto *diffMap = &gdata0.outMap->at(6);
  bool modelMap = cli.getOption('M');
  MU_WARN_IF(chi2Map.size(), fn,
             "Fit is 2-param, ignoring --single-temp map");

  // Thread-local storage.
  Graybody gray{*gdata0.gray};
  GrayData gdata{gdata0};
  gdata.gray = &gray;

  // GSL v1.x least-squares fitting setup.
  gsl_multifit_fdfsolver *solver =
    gsl_multifit_fdfsolver_alloc(gsl_multifit_fdfsolver_lmsder, // or lmder
                                 nBands, nParam);
  //
  gsl_multifit_function_fdf fdf;
  fdf.f = f_gray;
  fdf.df = df_gray;
  fdf.fdf = fdf_gray;
  fdf.n = nBands;
  fdf.p = nParam;
  fdf.params = &gdata;
  //
  // Initial guess.
  const double T0 = cli.getDouble('T', 16.5), S0 = 0.1;
  gsl_vector *x0 = gsl_vector_alloc(nParam);
  gsl_vector_set(x0, 0, T0);
  gsl_vector_set(x0, 1, S0);
  //
  // Covariance storage.
  gsl_matrix *cov = gsl_matrix_alloc(nParam, nParam);

  // Loop over map pixels.
  const bool absdiff = cli.getOption('a');
  const double relerr = cli.getDouble('A', 1e-2);
  size_t maxIter = std::max(1L, cli.getInteger('m', 25)),
         px = (modelMap ? modelMapLen*modelMapLen : bandData[0].size()),
         i0 = (id*px)/nThreads, i1 = ((id+1)*px)/nThreads;
  for (unsigned i=i0; i != i1; i++) {
    // Set target data.
    unsigned nOk = 0;
    for (size_t j=0; j != nBands; j++) {
      gdata.Fobs[j] = gdata.dnu0[j]*bandData[j][i];
      gdata.Ferr[j] = gdata.dnu0[j]*bandError[j][i];

      if (std::isfinite(gdata.Fobs[j]) && std::isfinite(gdata.Ferr[j])) {
        nOk++;
      } else {
        gdata.Fobs[j] = 0.0;
        gdata.Ferr[j] = 1e10;  // Weight-out bad data.
      }
    }

    // Progress display.
    MU_INFO_IF(id == 0 && i%std::max(i1/10, 1UL) == 0, fn,
               "Fitting solution %2.0f%% complete @ %s",
               (100.0*i)/i1, asctime());

    // Skip invalid pixels.
    if (nOk < minBands) {
      for (auto &omap: *gdata0.outMap) { omap[i] = NAN; }
      continue;
    }

    // Initialize solver.
    gsl_multifit_fdfsolver_set(solver, &fdf, x0);

    // Iterate on solution.
    size_t iter = 0;
    do {
      iter++;
      gsl_multifit_fdfsolver_iterate(solver);
      MU_DEBUG(fn, "[%2lu] T = %.2f  S = %.4f", iter,
               gsl_vector_get(solver->x, 0), gsl_vector_get(solver->x, 1));
    } while (gsl_multifit_test_delta(solver->dx, solver->x, 0.0, relerr) ==
             GSL_CONTINUE && iter != maxIter);

    // Results.
    gsl_multifit_covar(solver->J, 1e-8, cov);
    double chi2 = 0.0;
    for (size_t j=0; j != nBands; j++) {
      auto fj = gsl_vector_get(solver->f, j);
      diffMap[j][i] = fj*(absdiff ? bandError[j][i] : 1.0);
      chi2 += fj*fj;
    }
    //
    double chi2_dof = chi2/(nBands-nParam), vscale = std::max(1.0, chi2_dof),
           T = gsl_vector_get(solver->x, 0), S = gsl_vector_get(solver->x, 1),
           dT = sqrt(vscale * gsl_matrix_get(cov, 0, 0)),
           dS = sqrt(vscale * gsl_matrix_get(cov, 1, 1));
    MU_DEBUG(fn, "T = %.2f(%.2f), S = %.4f(%.4f), chi2/DOF = %.3f\n",
             T, dT, S, dS, chi2_dof);
    //
    Tmap[i] = T,       dTmap[i] = dT;
    Nmap[i] = S*imu_g, dNmap[i] = dS*imu_g;
    Cmap[i] = chi2_dof; Imap[i] = iter;
    //
    MU_WARN_IF(iter == maxIter, fn,
               "Fitting failed for pixel (%ld,%ld): dT/T = %.4g, dN/N = %.4g",
               1+i%gdata.axes[0], 1+i/gdata.axes[0],
               gsl_vector_get(solver->dx, 0)/T,
               gsl_vector_get(solver->dx, 1)/S);

    // Chain next initial guess to previous solution (if sane).
    // Vary guess for model maps to test convergence robustness.
    if (modelMap) { T += (i%7)-2.5, S *= 1.0 + ((i%11)-5.5)/10.0; }
    gsl_vector_set(x0, 0, T > 1.0  && T < 1e3 ? T : T0);
    gsl_vector_set(x0, 1, S > 1e-3 && S < 1e3 ? S : S0);
  }

  // Free storage.
  gsl_multifit_fdfsolver_free(solver);
  gsl_vector_free(x0);
  gsl_matrix_free(cov);
}


void solveStage2(const std::vector<mapDataType> &bandData,
                 const std::vector<mapDataType> &bandError,
                 const mapDataType &chi2Map,
                 const mu::CommandLine &cli,
                 const GrayData &gdata0, unsigned id, unsigned nThreads)
{
  const char *const fn = "solveStage2";
  size_t nBands = gdata0.nBands, nParam = gdata0.nParam,
         minBands = cli.getInteger('n', 4);
  auto &Tcore   = gdata0.outMap->at(0),
       &dTcore  = gdata0.outMap->at(1),
       &Ncore   = gdata0.outMap->at(2),
       &dNcore  = gdata0.outMap->at(3),
       &Thalo   = gdata0.outMap->at(4),
       &dThalo  = gdata0.outMap->at(5),
       &Nhalo   = gdata0.outMap->at(6),
       &dNhalo  = gdata0.outMap->at(7),
       &Cmap    = gdata0.outMap->at(8),
       &Imap    = gdata0.outMap->at(9);
  auto *diffMap = &gdata0.outMap->at(10);
  bool modelMap = cli.getOption('M');

  // Single-temperature chi^2 bypass mask.
  // Pixels whose chi2Map[i] < chi2Max are skipped (NaN-filled) here.
  double chi2Max = (chi2Map.size() ? cli.getDouble('c', 1.0) : -1.0);
  if (chi2Map.size() && chi2Map.size() != bandData[0].size()) {
    MU_ERROR(fn, "Ignoring --single-temp map (size %lu, require %lu)",
             chi2Map.size(), bandData[0].size());
    chi2Max = -1.0;
  }

  // Thread-local storage.
  Graybody gray{*gdata0.gray};
  GrayData gdata{gdata0};
  gdata.gray = &gray;

  // GSL v1.x least-squares fitting setup.
  gsl_multifit_fdfsolver *solver =
    gsl_multifit_fdfsolver_alloc(gsl_multifit_fdfsolver_lmsder, // or lmder
                                 nBands, nParam);
  //
  gsl_multifit_function_fdf fdf;
  fdf.f = f_gray2;
  fdf.df = df_gray2;
  fdf.fdf = fdf_gray2;
  fdf.n = nBands;
  fdf.p = nParam;
  fdf.params = &gdata;
  //
  // Initial guess.
  const double T0 = cli.getDouble('T', 16.5), S0 = 1.0, S1 = 0.1;
  gsl_vector *x0 = gsl_vector_alloc(nParam);
  gsl_vector_set(x0, 0, S0);
  gsl_vector_set(x0, 1, T0-5.0);
  gsl_vector_set(x0, 2, S1);
  if (nParam == 4) {
    gsl_vector_set(x0, 3, T0+5.0);
  }
  gdata.fixed.resize(1);  // Used when Th held fixed.
  //
  // Covariance storage.
  gsl_matrix *cov = gsl_matrix_alloc(nParam, nParam);

  // Loop over map pixels.
  const double relerr = cli.getDouble('A', 1e-2);
  size_t maxIter = std::max(1L, cli.getInteger('m', 25)),
         px = (modelMap ? modelMapLen*modelMapLen : bandData[0].size()),
         i0 = (id*px)/nThreads, i1 = ((id+1)*px)/nThreads;
  for (unsigned i=i0; i != i1; i++) {
    // Set target data.
    unsigned nOk = 0;
    for (size_t j=0; j != nBands; j++) {
      gdata.Fobs[j] = gdata.dnu0[j]*bandData[j][i];
      gdata.Ferr[j] = gdata.dnu0[j]*bandError[j][i];

      if (std::isfinite(gdata.Fobs[j]) && std::isfinite(gdata.Ferr[j])) {
        nOk++;
      } else {
        gdata.Fobs[j] = 0.0;
        gdata.Ferr[j] = 1e10;  // Weight-out bad data.
      }
    }
    gdata.fixed[0] = T0;

    // Progress display.
    MU_INFO_IF(id == 0 && i%(i1/10) == 0, fn,
               "Fitting solution %2.0f%% complete @ %s",
               (100.0*i)/i1, asctime());

    // Skip invalid/masked pixels.
    if (nOk < minBands || (chi2Max >= 0.0 && chi2Map[i] <= chi2Max)) {
      for (auto &omap: *gdata0.outMap) { omap[i] = NAN; }
      continue;
    }

    // Initialize solver.
    gsl_multifit_fdfsolver_set(solver, &fdf, x0);

    // Iterate on solution.
    size_t iter = 0;
    do {
      MU_DEBUG(fn, "[%2lu] Tc = %5.2f  Sc = %7.4f "
                         " Th = %5.2f  Sh = %7.4f", iter,
               gsl_vector_get(solver->x, 1), gsl_vector_get(solver->x, 0),
               nParam == 4 ? gsl_vector_get(solver->x, 3) : gdata.fixed[0],
               gsl_vector_get(solver->x, 2));
      iter++;
      gsl_multifit_fdfsolver_iterate(solver);
    } while (gsl_multifit_test_delta(solver->dx, solver->x, 0.0, relerr) ==
             GSL_CONTINUE && ++iter != maxIter);

    // Results.
    gsl_multifit_covar(solver->J, 1e-6, cov);
    double chi2 = 0.0;
    for (size_t j=0; j != nBands; j++) {
      auto fj = gsl_vector_get(solver->f, j);
      diffMap[j][i] = fj;
      chi2 += fj*fj;
    }
    //
    double chi2_dof = chi2/std::max(nBands-nParam, 1UL),
           vscale = std::max(1.0, chi2_dof),
           Sc = gsl_vector_get(solver->x, 0),
           Tc = gsl_vector_get(solver->x, 1),
           Sh = gsl_vector_get(solver->x, 2),
           Th = nParam == 4 ? gsl_vector_get(solver->x, 3) : gdata.fixed[0],
           dSc = sqrt(vscale * gsl_matrix_get(cov, 0, 0)),
           dTc = sqrt(vscale * gsl_matrix_get(cov, 1, 1)),
           dSh = sqrt(vscale * gsl_matrix_get(cov, 2, 2)),
           dTh = nParam == 4 ? sqrt(vscale * gsl_matrix_get(cov, 3, 3)) : 0.0;
    MU_DEBUG(fn, "(%ld,%ld) Tc = %5.2f(%.2f), Sc = %7.4f(%.4f), "
                 "Th = %5.2f(%.2f), Sh = %7.4f(%.4f), chi2/DOF = %.3f\n",
             1+i%gdata.axes[0], 1+i/gdata.axes[0],
             Tc, dTc, Sc, dSc, Th, dTh, Sh, dSh, chi2_dof);
    //
    Tcore[i] = Tc,       dTcore[i] = dTc;
    Ncore[i] = Sc*imu_g, dNcore[i] = dSc*imu_g;
    Thalo[i] = Th,       dThalo[i] = dTh;
    Nhalo[i] = Sh*imu_g, dNhalo[i] = dSh*imu_g;
     Cmap[i] = chi2_dof;   Imap[i] = iter;
    //
    MU_WARN_IF(iter == maxIter, fn,
               "Fitting failed for pixel (%ld,%ld): "
               "dTc/Tc = %.4g, dNc/Nc = %.4g, "
               "dTh/Th = %.4g, dNh/Nh = %.4g",
               1+i%gdata.axes[0], 1+i/gdata.axes[0],
               gsl_vector_get(solver->dx, 1)/Tc,
               gsl_vector_get(solver->dx, 0)/Sc,
               nParam == 4 ? gsl_vector_get(solver->dx, 3)/Th : 0.0,
               gsl_vector_get(solver->dx, 2)/Sh);

    // Chain next initial guess to previous solution (if sane).
    // !!!! Disabled pending resolution of GSL integrator convergence.
    gsl_vector_set(x0, 0, S0);//Sc > 1e-2 && Sc < 10.0 ? Sc : S0);
    gsl_vector_set(x0, 1, T0-5.0);//Tc > 3.0 && Tc < 100.0 ? Tc : T0-5.0);
    gsl_vector_set(x0, 2, S1);//Sh > 1e-3 && Sh < 1e3 ? Sh : 0.05*S0);
    if (nParam == 4) {
      gsl_vector_set(x0, 3, T0+5.0);//Th > 10.0 && Th < 100.0 ? Th : T0+5.0);
    }
  }

  // Free storage.
  gsl_multifit_fdfsolver_free(solver);
  gsl_vector_free(x0);
  gsl_matrix_free(cov);
}


mapDataType getChi2(const mu::CommandLine &cli)
{
  mapDataType chi2;
  mu::CommandLine::Arguments files;

  if (cli.getOption('s', &files, 1)) {
    const std::string hduName = "Chi2/DOF";
    std::unique_ptr<CCfits::FITS>
      fchi2{new CCfits::FITS(files.front(), CCfits::Read, hduName, true)};
    fchi2->extension(hduName).read(chi2);
  }

  return chi2;
}


void solve(std::vector<mapDataType> &outMap,
           std::vector<CCfits::ExtHDU *> &outHDU,
           CCfits::FITS *outFITS,
           const std::vector<mapDataType> &bandData,
           const std::vector<mapDataType> &bandError,
           const std::vector<CCfits::PHDU*> &bandHDU,
           const mu::CommandLine &cli)
{
  const char *const fn = "solve";
  bool modelMap = cli.getOption('M');

  // Dust model.
  bool extend = !cli.getOption('p');
  Dust dust(cli.getString('D', "OH5"), cli.getDouble('g', 100.0));
  MU_INFO(fn, "Using %s source detector response",
          extend ? "extended" : "point");
  MU_INFO(fn, "Dust model: %s, gas-to-dust ratio = %g",
          dust.name(), dust.rho());

  // Band detector setup.
  // Currently supports PACS/SPIRE only!
  size_t nBands = bandData.size();
  std::vector<std::unique_ptr<Detector>> detect;
  for (size_t j=0; j != nBands; j++) {
    int band = getBand(*bandHDU[j]);
    auto d = (band < 200 ? (Detector *)new PACS(band, extend)
                         : (Detector *)new SPIRE(band, extend));
    d->init();
    detect.emplace_back(d);
  }

  // Allocate result maps.
  const double relerr = cli.getDouble('A', 1e-2);
  std::vector<long> axes = {bandHDU[0]->axis(0), bandHDU[0]->axis(1)};
  if (modelMap) { axes = {modelMapLen, modelMapLen}; }

  auto px = axes[0]*axes[1];
  char comm[1024];
  snprintf(comm, sizeof(comm),
           "%s sources, %s dust, gas/dust = %g, relerr = %.1e",
           (extend ? "Extended" : "Point"), dust.name().c_str(),
           dust.rho(), relerr);
  //
  size_t nParam = (cli.getOption('4') ? 4 : cli.getOption('3') ? 3 : 2);
  auto solver = &solveStage1;
  std::vector<std::string>
    mapNames = {"T", "dT", "N(H2)", "dN(H2)", "Chi2/DOF", "nIter"};
  if (nParam > 2) {
    solver = &solveStage2;
    mapNames = {"Tc", "dTc", "Nc(H2)", "dNc(H2)",
                "Th", "dTh", "Nh(H2)", "dNh(H2)", "Chi2/DOF", "nIter"};
  }
  for (auto &b: bandHDU) {
    // Band fit residuals.
    mapNames.push_back("diff" + std::to_string(getBand(*b)));
  }
  //
  for (auto name: mapNames) {
    outMap.push_back(mapDataType(0.0, px));
    outHDU.push_back(outFITS->addImage(name, mapDataCode, axes));
    outHDU.back()->writeDate();
    outHDU.back()->writeComment(comm);
  }

  // Emission model.
  constexpr double MJy = 1e-17;  // MJy -> erg/cm^2/s/Hz
  Graybody gray{dust, 0.1*relerr};
  GrayData gdata{&gray, &outMap, nBands, nParam};
  for (size_t j=0; j != nBands; j++) {
    gdata.band[j] = detect[j].get();
    gdata.dnu0.push_back(detect[j]->freqWidth() * MJy);
  }
  gdata.axes = axes;

  // Theoretical model map option.
  if (modelMap) {
    makeModelMap(const_cast<std::vector<mapDataType>&>(bandData),
                 const_cast<std::vector<mapDataType>&>(bandError),
                 gdata);
  }

  // Create threads and process maps.
  unsigned nThreads = std::max(1L, cli.getInteger('t', 1));
  auto chi2Map = getChi2(cli);
  std::vector<std::thread> threads;
  for (unsigned i=1; i != nThreads; i++) {
    threads.emplace_back(solver, bandData, bandError, chi2Map, cli,
                         gdata, i, nThreads);
  }
  solver(bandData, bandError, chi2Map, cli, gdata, 0, nThreads);
  //
  for (auto &t: threads) { t.join(); }

  MU_INFO(fn, "Solution complete @ %s", asctime());
}

}  // namespace manticore