Skip to content

GitLab

Commit 1393d52d authored by Zhang, Chen's avatar Zhang, Chen
Browse files

cleanup PR to speed up review process

parent 73a6c68c
Hide whitespace changes
Inline Side-by-side
Framework/Crystal/inc/MantidCrystal/SCDCalibratePanels2.h
View file @ 1393d52d
......@@ -87,18 +87,6 @@ private:
/// Private function for calibrating banks
void optimizeBanks(Mantid::API::IPeaksWorkspace_sptr pws);
/// Twiddle search for getting L1
double twiddle_search_L1(Mantid::API::IPeaksWorkspace_sptr pws,
double deltaL1, double threshold);
double objfunc_L1(Mantid::API::IPeaksWorkspace_sptr pws, double source_z,
double init_z);
/// Twiddle search for calibrating bank
void twiddle_search_banks(Mantid::API::IPeaksWorkspace_sptr pws,
double searchStep[6], double threshold);
double objfunc_bank(Mantid::API::IPeaksWorkspace_sptr pwsBank,
double params[6], std::string bankname);
/// Helper function for selecting peaks based on given bank name
Mantid::API::IPeaksWorkspace_sptr
selectPeaksByBankName(Mantid::API::IPeaksWorkspace_sptr pws,
......
Framework/Crystal/src/SCDCalibratePanels2.cpp
View file @ 1393d52d
......@@ -266,38 +266,11 @@ void SCDCalibratePanels2::exec() {
if (calibrateL1) {
g_log.notice() << "** Calibrating L1 (moderator) as requested\n";
optimizeL1(m_pws);
///
/// Twiddle search here is for dev to explore the parameter space
///
// double search_step_L1 = 1e-4; // 100um search step
// double threshold_L1 = 1e-6; // when to stop
// double new_L1 = twiddle_search_L1(m_pws, search_step_L1, threshold_L1);
//
// g_log.notice() << "** -- New L1 = " << new_L1 << "\n";
//
// g_log.notice() << "Profiling objfunc sensitivity (coarse) \n";
// std::ofstream proffile;
// proffile.precision(17);
// proffile.open("/tmp/prof_unittest_nullcase.csv");
// proffile << "dz\terr\n";
// std::cout.precision(17);
// for (double z = -0.02; z <= 0.02; z = z + 1e-5) {
// double err = objfunc_L1(m_pws, -20 + z, -20);
// std::cout << std::scientific << z << "\t" << std::fixed << err << "\n";
// proffile << std::scientific << z << "\t" << std::fixed << err << "\n";
// }
// proffile.close();
}
if (calibrateBanks) {
g_log.notice() << "** Calibrating L2 and orientation (bank) as requested\n";
optimizeBanks(m_pws);
// dx, dy, dz --> 1um
// theta, phi, rotang --> 1e-3 deg
// double search_step_bank[6] = {1e-6, 1e-6, 1e-6, 0.0, 0.0, 0.0};
// double threshold_bank = 1e-6;
// twiddle_search_banks(m_pws, search_step_bank, threshold_bank);
}
// STEP_3: generate a table workspace to save the calibration results
......@@ -344,7 +317,7 @@ void SCDCalibratePanels2::optimizeT0(IPeaksWorkspace_sptr pws) {
// class just to get Fit serving as an optimizer.
// For this particular case, we are constructing an objective function
// based on IFunction1D that outputs a fake histogram consist of
// qSample calucated based on perturbed instrument positions and
// qSample calculated based on perturbed instrument positions and
// orientations.
MatrixWorkspace_sptr t0ws = getIdealQSampleAsHistogram1D(pws);
......@@ -386,7 +359,8 @@ void SCDCalibratePanels2::optimizeT0(IPeaksWorkspace_sptr pws) {
*/
void SCDCalibratePanels2::optimizeL1(IPeaksWorkspace_sptr pws) {
// cache starting L1 position
double original_L1 = -pws->getInstrument()->getSource()->getPos().Z();
double original_L1 =
std::abs(pws->getInstrument()->getSource()->getPos().Z());
MatrixWorkspace_sptr l1ws = getIdealQSampleAsHistogram1D(pws);
......@@ -529,322 +503,6 @@ void SCDCalibratePanels2::optimizeBanks(IPeaksWorkspace_sptr pws) {
PARALLEL_CHECK_INTERUPT_REGION
}
/**
* @brief Use twiddle search optimizer to find the L1 position by minimizing
* the qSample calculated from Miller-indices (int) and qSample from
* purturbed instrument.
*
* @param pws Input peak workspace
* @param deltaL1 Initial search step in meters
* @param threshold Stop searching when deltaL1 is below this value
* @return double
*/
double SCDCalibratePanels2::twiddle_search_L1(IPeaksWorkspace_sptr pws,
double deltaL1,
double threshold) {
// calculate the error
// NOTE: use the current position as the starting point
double init_L1 = pws->getInstrument()->getSource()->getPos().Z();
double L1 = init_L1;
double best_err = objfunc_L1(pws, L1, init_L1);
int niter = 0;
g_log.notice() << "-- twiddle search\n"
<< " err_0 = " << best_err << "\n";
while (deltaL1 > threshold) {
L1 += deltaL1;
double err = objfunc_L1(pws, L1, init_L1);
if (err < best_err) {
// There was some improvement
best_err = err;
deltaL1 *= 1.1;
} else {
// There was no improvement
L1 -= 2.0 * deltaL1;
err = objfunc_L1(pws, L1, init_L1);
if (err < best_err) {
// There was an improvement
best_err = err;
deltaL1 *= 1.05;
} else {
L1 += deltaL1;
// As there was no improvement, the step size in either
// direction, the step size might simply be too big.
deltaL1 *= 0.95;
}
}
niter += 1;
std::ostringstream msg;
msg.precision(17);
msg << "-----@iter_" << niter << "-----\n"
<< " deltaL1 = " << deltaL1 << "\n"
<< " L1 = " << L1 << "\n"
<< " besterr = " << best_err << "\n";
g_log.notice() << msg.str();
}
// update the source position
IAlgorithm_sptr mv_alg =
createChildAlgorithm("MoveInstrumentComponent", -1, -1, false);
mv_alg->setLogging(LOGCHILDALG);
mv_alg->setProperty<Workspace_sptr>("Workspace", pws);
mv_alg->setProperty("ComponentName",
pws->getInstrument()->getSource()->getName());
mv_alg->setProperty("X", 0.0);
mv_alg->setProperty("Y", 0.0);
mv_alg->setProperty("Z", L1 - init_L1);
mv_alg->setProperty("RelativePosition", true);
mv_alg->executeAsChildAlg();
return L1;
}
/**
* @brief The objective function computes the difference between
* two errors (one use tof@n-1 and one use tof@n). The intecept
* of the two errors is the solution
*
* @param pws Input peak workspace
* @param source_z Current source position (L1)
* @param init_z Initial source position (L1 engineering value)
* @return double
*/
double SCDCalibratePanels2::objfunc_L1(IPeaksWorkspace_sptr pws_org,
double source_z, double init_z) {
IPeaksWorkspace_sptr pws = pws_org->clone();
double err = 0.0;
V3D qv_target;
V3D qv_calc;
V3D delta_qv;
Units::Wavelength wl;
// move the source to new location
double z_shift = source_z - init_z;
IAlgorithm_sptr mv_alg =
createChildAlgorithm("MoveInstrumentComponent", -1, -1, false);
mv_alg->setLogging(LOGCHILDALG);
mv_alg->setProperty<Workspace_sptr>("Workspace", pws);
mv_alg->setProperty("ComponentName",
pws->getInstrument()->getSource()->getName());
mv_alg->setProperty("X", 0.0);
mv_alg->setProperty("Y", 0.0);
mv_alg->setProperty("Z", z_shift);
mv_alg->setProperty("RelativePosition", true);
mv_alg->executeAsChildAlg();
// calculate the residual
auto ubmatrix = pws->sample().getOrientedLattice().getUB();
for (int i = 0; i < pws->getNumberPeaks(); ++i) {
// get the reference
qv_target = ubmatrix * pws_org->getPeak(i).getIntHKL();
qv_target *= 2.0 * PI;
// cache time-of-flight
double tof = pws->getPeak(i).getTOF();
// make a peak that has the new instrument
// NOTE:
// Change the instrument in a peak workspace does not
// trigger an auto update of all the peaks inside,
// which is why we need to do it here explicitly
Peak pk = Peak(pws->getPeak(i));
pk.setInstrument(pws->getInstrument());
wl.initialize(pk.getL1(), pk.getL2(), pk.getScattering(), 0,
pk.getInitialEnergy(), 0.0);
pk.setDetectorID(pws->getPeak(i).getDetectorID());
pk.setWavelength(wl.singleFromTOF(tof));
// calculate the residual
qv_calc = pk.getQSampleFrame();
delta_qv = qv_calc - qv_target;
err += delta_qv.norm2();
}
err /= pws->getNumberPeaks();
return err;
}
void SCDCalibratePanels2::twiddle_search_banks(IPeaksWorkspace_sptr pws,
double searchStep[6],
double threshold) {
PARALLEL_FOR_IF(Kernel::threadSafe(*pws))
for (int i = 0; i < static_cast<int>(m_BankNames.size()); ++i) {
PARALLEL_START_INTERUPT_REGION
const std::string bankname = *std::next(m_BankNames.begin(), i);
const std::string pwsBankiName = "_pws_" + bankname;
//-- step 0: extract peaks that lies on the current bank
IPeaksWorkspace_sptr pwsBanki =
selectPeaksByBankName(pws, bankname, pwsBankiName);
// Do not attempt correct panels with less than 6 peaks
int nBankPeaks = pwsBanki->getNumberPeaks();
if (nBankPeaks < MINIMUM_PEAKS_PER_BANK) {
// use ostringstream to prevent OPENMP breaks log info
std::ostringstream msg_npeakCheckFail;
msg_npeakCheckFail << "-- Bank " << bankname << " have only "
<< nBankPeaks << " (<" << MINIMUM_PEAKS_PER_BANK
<< ") Peaks, skipping\n";
g_log.notice() << msg_npeakCheckFail.str();
continue;
}
//-- step 1: twiddle search
double params[6] = {0.0, 0.0, 0.0, // dx, dy, dz
0.0, 0.0, 0.0}; // theta, phi, ang
double best_err = objfunc_bank(pwsBanki, params, bankname);
double err;
int niter = 0;
double stepsSum = 0.0;
for (int i = 0; i < 6; ++i)
stepsSum += searchStep[i];
// step size adjustment is fixed at 5%
while (stepsSum > threshold) {
// iterate through feature vector (params)
// assuming they are independent from each other
for (int i = 0; i < 6; ++i) {
// search positive direction +1
params[i] += searchStep[i];
err = objfunc_bank(pwsBanki, params, bankname);
if (err < best_err) {
// got a hit, increase step size
best_err = err;
searchStep[i] *= 1.1;
} else {
// search negative direction +1-2 => -1
params[i] -= 2 * searchStep[i];
err = objfunc_bank(pwsBanki, params, bankname);
if (err < best_err) {
// got a hit, increase step size
best_err = err;
searchStep[i] *= 1.05;
} else {
// no hit, step size might be too large,
// reset posiiton +1-2+1 => 0
params[i] += searchStep[i];
// reduce step size
searchStep[i] *= 0.95;
}
}
}
niter += 1;
stepsSum = 0.0;
for (int i = 0; i < 6; ++i)
stepsSum += searchStep[i];
// logging
V3D xyz = V3D(params[0], params[1], params[2]);
double _theta = params[3];
double _phi = params[4];
double _ang = params[5];
V3D rotvec =
V3D(sin(_theta) * cos(_phi), sin(_theta) * sin(_phi), cos(_theta));
std::ostringstream logmsg;
logmsg.precision(12);
logmsg << bankname << "@iter_" << niter << "\n"
<< "-- best_err: " << best_err << "\n"
<< "-- (dx,dy,dz): " << xyz << "\n"
<< "-- ang@rotvec: " << _ang << "@" << rotvec << "\n";
g_log.notice() << logmsg.str();
}
//-- step 2: move bank to corresponding location
const double dx = params[0];
const double dy = params[1];
const double dz = params[2];
const double theta = params[3];
const double phi = params[4];
const double rotang = params[5];
double rvx = sin(theta) * cos(phi);
double rvy = sin(theta) * sin(phi);
double rvz = cos(theta);
std::string bn = bankname;
if (pws->getInstrument()->getName().compare("CORELLI") == 0)
bn.append("/sixteenpack");
adjustComponent(dx, dy, dz, rvx, rvy, rvz, rotang, bn, pws);
//-- step 3: log
V3D dxyz = V3D(dx, dy, dz);
V3D rv = V3D(rvx, rvy, rvz);
std::ostringstream msg;
msg.precision(8);
msg << bankname << " calibration results:\n"
<< "(dx, dy, dz) = " << dxyz << "\n"
<< "ang@rot-axis = " << rotang << "@" << rv << "\n";
g_log.notice() << msg.str();
PARALLEL_END_INTERUPT_REGION
}
PARALLEL_CHECK_INTERUPT_REGION
}
double SCDCalibratePanels2::objfunc_bank(IPeaksWorkspace_sptr pwsBank,
double params[6],
std::string bankname) {
// prep
double err = 0.0;
V3D qv_target;
V3D qv_calc;
V3D delta_qv;
Units::Wavelength wl;
// only operate on the clone
IPeaksWorkspace_sptr pws = pwsBank->clone();
// purturb the instrument
const double dx = params[0];
const double dy = params[1];
const double dz = params[2];
const double theta = params[3];
const double phi = params[4];
const double rotang = params[5];
double rvx = sin(theta) * cos(phi);
double rvy = sin(theta) * sin(phi);
double rvz = cos(theta);
std::string bn = bankname;
if (pws->getInstrument()->getName().compare("CORELLI") == 0)
bn.append("/sixteenpack");
adjustComponent(dx, dy, dz, rvx, rvy, rvz, rotang, bn, pws);
// calculate the residual
auto ubmatrix = pws->sample().getOrientedLattice().getUB();
for (int i = 0; i < pws->getNumberPeaks(); ++i) {
qv_target = ubmatrix * pwsBank->getPeak(i).getIntHKL();
qv_target *= 2 * PI;
// cache time-of-flight
double tof = pws->getPeak(i).getTOF();
// get qv with updated instrument
Peak pk = Peak(pws->getPeak(i));
pk.setInstrument(pws->getInstrument());
wl.initialize(pk.getL1(), pk.getL2(), pk.getScattering(), 0,
pk.getInitialEnergy(), 0.0);
pk.setDetectorID(pws->getPeak(i).getDetectorID());
pk.setWavelength(wl.singleFromTOF(tof));
// calculate the residual
qv_calc = pk.getQSampleFrame();
delta_qv = qv_calc - qv_target;
err += delta_qv.norm2();
}
err /= pws->getNumberPeaks();
return err;
}
/// ---------------- ///
/// helper functions ///
/// ---------------- ///
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment