/*
 * The MIT License

   Conograph (powder auto-indexing program)

Copyright (c) <2012> <Ryoko Oishi-Tomiyasu, KEK>

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
 *
 */
#ifdef _OPENMP
 # include <omp.h>
#endif
#include "utility_data_structure/index_set.hh"
#include "utility_func/chToDouble.hh"
#include "utility_func/transform_sym_matrix.hh"
#include "utility_lattice_reduction/super_basis3.hh"
#include "lattice_symmetry/LatticeFigureOfMeritToCheckSymmetry.hh"
#include "lattice_symmetry/ReducedLatticeToCheckBravais.hh"
#include "lattice_symmetry/ReducedLatticeToCheckEquiv.hh"
#include "p_out_indexing.hh"
#include "zerror_type/error_out.hh"
#include "zlog/zlog.hh"
#include "ControlParam.hh"
#include "SortingLattice.hh"
#include "utility_func/stopx.hh"

const bool SortingLattice::m_DoesPrudentSymSearch = false;
const Double SortingLattice::m_cv2 = 0.5;

SortingLattice::SortingLattice()
{
	for(ArrayIndex i=0; i<NUM_LS; i++)
	{
		OutputSymmetry[i] = false;
		JudgeSymmetry[i] = false;
	}
	
   	m_resol2 = 0.0;
	m_num_ref_figure_of_merit = 20;
	m_etype_peak_shift = kPeakShiftFunction_Type0;
	m_WlengthX = 1.54056;
}


SortingLattice::~SortingLattice()
{
}


// Set the member variables.
void SortingLattice::setParam(const ControlParam& cont) 
{
	OutputSymmetry[(Int4)Triclinic] = cont.putOutputSymmetry(Triclinic);
	JudgeSymmetry[(Int4)Triclinic] = false;
	for(ArrayIndex i=1; i<NUM_LS; i++)
	{
		OutputSymmetry[i] = cont.putOutputSymmetry(eCrystalSystem(i));
		JudgeSymmetry[i] = cont.putOutputSymmetry(eCrystalSystem(i));
	}

	if( JudgeSymmetry[(Int4)Cubic_P] )
	{
		JudgeSymmetry[(Int4)Tetragonal_P] = true;
	}
	if( JudgeSymmetry[(Int4)Hexagonal] )
	{
		JudgeSymmetry[(Int4)Monoclinic_P] = true;
	}
	if( JudgeSymmetry[(Int4)Tetragonal_P] )
	{
		JudgeSymmetry[(Int4)Orthorhombic_P] = true;
	}
	if( JudgeSymmetry[(Int4)Orthorhombic_P] )
	{
		JudgeSymmetry[(Int4)Monoclinic_P] = true;
	}
	
	if( JudgeSymmetry[(Int4)Orthorhombic_C] )
	{
		JudgeSymmetry[(Int4)Monoclinic_B] = true;
	}

	if( JudgeSymmetry[(Int4)Cubic_I] )
	{
		JudgeSymmetry[(Int4)Tetragonal_I] = true;
	}
	if( JudgeSymmetry[(Int4)Tetragonal_I] )
	{
		JudgeSymmetry[(Int4)Orthorhombic_I] = true;
	}

	if( JudgeSymmetry[(Int4)Cubic_F] )
	{
		JudgeSymmetry[(Int4)Orthorhombic_F] = true;
	}

	m_resol2 = cont.putResolution() * 2.0;
	m_num_ref_figure_of_merit = cont.putNumberOfReflectionsForFigureOfMerit();
	m_etype_peak_shift = cont.putPeakShiftFunctionType();
	m_WlengthX = cont.putWaveLength();

	const vector<Double>& peak_shift_param_rad = cont.putPeakShiftParamRadian();
	const Int4 param_num = peak_shift_param_rad.size();
	assert( m_etype_peak_shift != kPeakShiftFunction_Type0 || param_num == 0 );
	assert( m_etype_peak_shift != kPeakShiftFunction_Type1 || param_num == 1 );

	m_peak_shift_param_rad.resize(param_num);
	for(Int4 i=0; i<param_num; i++)
	{
		m_peak_shift_param_rad[i] = peak_shift_param_rad[i];
	}
}


void SortingLattice::putLatticeFigureOfMerit(const LatticeFigureOfMeritToCheckSymmetry& lattice_original,
		const ePointGroup& epg, const Double& cv2,
		vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result) const
{
	lattice_result.clear();
	map< SymMat<VCData>, NRMat<Int4> > S_red_tray;
	if( !lattice_original.checkLatticeSymmetry(epg, cv2, S_red_tray) ) return;

	const BravaisType& ebrat_original = lattice_original.putLatticeFigureOfMerit().putBravaisType();
	const eBravaisLattice eblat = (ebrat_original.enumCrystalSystem()==Monoclinic_B?
									(epg==D31d_rho?Prim:(epg==D3d_1_hex?Rhom_hex:BaseZ)):ebrat_original.enumBravaisLattice());

	const NRMat<Int4> matrix_min_to_sell = lattice_original.putInitialForm().second;

	SymMat<Double> S_super(4);
	Int4 itnum;
	NRMat<Int4> trans_mat(4,3);
	
	for(map< SymMat<VCData>, NRMat<Int4> >::const_iterator it=S_red_tray.begin(); it!=S_red_tray.end(); it++)
	{
		// S_super = it->second * it->first * Transpose(it->second) is close to
		// Delone-reduced form of the original lattice.
		S_super = transform_sym_matrix(it->second, chToDouble(it->first) );

		trans_mat = identity_matrix<Int4>(4);

		// S_super = trans_mat * it->second * it->first * Transpose(trans_mat * it->second).
		put_super_Gram_matrix_obtuse_angle< Double, SymMat<Double> >(trans_mat, S_super, itnum);
		moveSmallerDiagonalLeftUpper< Double, SymMat<Double> >(S_super, trans_mat);
	
		lattice_result.push_back( LatticeFigureOfMeritToCheckSymmetry( BravaisType( pair<eBravaisLattice, ePointGroup>(eblat, epg) ),
										SymMat43_VCData(it->first, mprod(trans_mat, it->second) ),
										lattice_original.putLatticeFigureOfMerit().putPeakShiftFunctionType(),
										lattice_original.putLatticeFigureOfMerit().putWaveLength(),
										lattice_original.putLatticeFigureOfMerit().putPeakShiftParamRadian() ) );
	}
}


void SortingLattice::putBravaisLatticeFigureOfMerit(const ReducedLatticeToCheckBravais& RLCB,
		const LatticeFigureOfMeritToCheckSymmetry& lattice_original,
		const BravaisType& brat,
		vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result) const
{
	lattice_result.clear();
	
	const map< SymMat<VCData>, NRMat<Int4> >& S_red_tray = RLCB.checkBravaisLatticeType(brat);
	if( S_red_tray.empty() ) return;

	// The lattice of RLCB has at least the symmetry given by eblat.
	SymMat<VCData> S_super_obtuse(4);
	Int4 itnum;
	NRMat<Int4> trans_mat(4,3);

	for(map< SymMat<VCData>, NRMat<Int4> >::const_iterator it=S_red_tray.begin(); it!=S_red_tray.end(); it++)
	{
		S_super_obtuse = transform_sym_matrix(it->second, it->first);
		trans_mat = identity_matrix<Int4>(4);

		// S_super = trans_mat * it->second * it->first * Transpose(trans_mat * it->second) is Delone reduced.
		if( !put_super_Gram_matrix_obtuse_angle< VCData, SymMat<VCData> >(trans_mat, S_super_obtuse, itnum) )
		{
			assert( false );
		}
		moveSmallerDiagonalLeftUpper< VCData, SymMat<VCData> >(S_super_obtuse, trans_mat);
		lattice_result.push_back( LatticeFigureOfMeritToCheckSymmetry( brat, SymMat43_VCData(it->first, mprod(trans_mat, it->second) ),
									lattice_original.putLatticeFigureOfMerit().putPeakShiftFunctionType(),
									lattice_original.putLatticeFigureOfMerit().putWaveLength(),
									lattice_original.putLatticeFigureOfMerit().putPeakShiftParamRadian() ) );
	}
}






void SortingLattice::putLatticeCandidatesForTriclinic(const vector<SymMat43_VCData>& S_super,
		const Double& MIN_NormM,
		const Double& MIN_RevM,
		vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_tri) const
{
	const Int4 num_topo = S_super.size();
	lattice_result_tri.clear();

/*   2011.10.19 VIC Tamura */
Int4 LOOP_COUNTER = 0;
	
#ifdef _OPENMP
	#pragma omp parallel
#endif
	{
		vector< VecDat3<Int4> > closest_hkl_tray;
		vector<bool> is_cal_Q_observed_tray;
		vector<LatticeFigureOfMeritToCheckSymmetry> latFOM_tray;

#ifdef _OPENMP
		#pragma omp for
#endif
		for(Int4 n=0; n<num_topo; n++)
	   	{
/*   2011.10.19 VIC Tamura */
SET_PROGRESS(100*(LOOP_COUNTER++)/num_topo, 65, 1); // critical, but works
if(IS_CANSELED()) continue;

			LatticeFigureOfMeritToCheckSymmetry latFOM(BravaisType( pair<eBravaisLattice, ePointGroup>(Prim, Ci) ), S_super[n],
											m_etype_peak_shift, m_WlengthX, m_peak_shift_param_rad);

			latFOM.setFigureOfMerit(m_num_ref_figure_of_merit,
										VCData::putPeakQData(),
										closest_hkl_tray, is_cal_Q_observed_tray);

			LatticeFigureOfMeritZeroShift latFOM2 = latFOM.putLatticeFigureOfMerit();
			pair<bool, ZErrorMessage> ans = latFOM2.fitLatticeParameterLinear(VCData::putPeakQData(),
												closest_hkl_tray, is_cal_Q_observed_tray, false);

			if( ans.first )
			{
				assert( latFOM.putLatticeFigureOfMerit().putFiguresOfMerit().putNumberOfReflectionsForFigureOfMerit() > 0 );
				latFOM2.setFigureOfMerit(latFOM.putLatticeFigureOfMerit().putFiguresOfMerit().putNumberOfReflectionsForFigureOfMerit(),
											VCData::putPeakQData(),
											closest_hkl_tray, is_cal_Q_observed_tray);
				if( LatticeFigureOfMerit::cmpFOMdeWolff( latFOM2, latFOM.putLatticeFigureOfMerit() ) )
				{
					latFOM.setLatticeFigureOfMerit(latFOM2);
				}
			}
			const SetOfFigureOfMerit& setFOM = latFOM.putLatticeFigureOfMerit().putFiguresOfMerit();
			if( setFOM.putFigureOfMeritWu() < MIN_NormM ) continue;
			if( setFOM.putReversedFigureOfMerit() < MIN_RevM ) continue;

#ifdef _OPENMP
		#pragma omp critical
#endif
			{
				lattice_result_tri.push_back( latFOM );
			}
	   	}
	}

/*   2011.10.19 VIC Tamura */
CHECK_INTERRUPTION();

	sort( lattice_result_tri.begin(), lattice_result_tri.end() );
}


void SortingLattice::putLatticeCandidatesForEachBravaisTypes(
					const Double& MIN_NormM,
					const Double& MIN_RevM,
					const eABCaxis& abc_axis,
					const eRHaxis& rh_axis,
					vector<LatticeFigureOfMeritToCheckSymmetry> lattice_result[NUM_LS]) const
{
	try{

		for(ArrayIndex i=1; i<NUM_LS; i++)
	{
		lattice_result[i].clear();
	}
	vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_tri = lattice_result[(ArrayIndex)Triclinic];
	vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_mono_P = lattice_result[(ArrayIndex)Monoclinic_P];
	vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_mono_B = lattice_result[(ArrayIndex)Monoclinic_B];
	vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_ortho_P = lattice_result[(ArrayIndex)Orthorhombic_P];
	vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_ortho_B = lattice_result[(ArrayIndex)Orthorhombic_C];
	vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_ortho_I = lattice_result[(ArrayIndex)Orthorhombic_I];
	vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_ortho_F = lattice_result[(ArrayIndex)Orthorhombic_F];
	vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_tetra_P = lattice_result[(ArrayIndex)Tetragonal_P];
	vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_tetra_I = lattice_result[(ArrayIndex)Tetragonal_I];
	vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_rhom = lattice_result[(ArrayIndex)Rhombohedral];
	vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_hex = lattice_result[(ArrayIndex)Hexagonal];
	vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_cubic_P = lattice_result[(ArrayIndex)Cubic_P];
	vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_cubic_I = lattice_result[(ArrayIndex)Cubic_I];
	vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_cubic_F = lattice_result[(ArrayIndex)Cubic_F];

	const Int4 num_tri = lattice_result_tri.size();

/*   2011.10.19 VIC Tamura */
Int4 LOOP_COUNTER = 0;

#ifdef _OPENMP
	#pragma omp parallel
#endif
	{
		vector<LatticeFigureOfMeritToCheckSymmetry> latFOM_tray;

#ifdef _OPENMP
		#pragma omp for
#endif
		for(Int4 n=0; n<num_tri; n++)
	   	{
/*   2011.10.19 VIC Tamura */
SET_PROGRESS(99*(LOOP_COUNTER++)/num_tri, 66, 30); // critical, but works
if(IS_CANSELED()) continue;

			LatticeFigureOfMeritToCheckSymmetry& latFOM = lattice_result_tri[n];
			latFOM.setLabel(n+1);
			const SymMat43_VCData& S_red = latFOM.putInitialForm();

			// Calculate figures of merit as triclinic
			const ReducedLatticeToCheckBravais RLCB(abc_axis, rh_axis, m_DoesPrudentSymSearch, m_cv2, S_red);

			if( JudgeSymmetry[Monoclinic_B] )
			{
				putBravaisLatticeFigureOfMerit(RLCB, latFOM, BravaisType( put_monoclinic_b_type(abc_axis) ), latFOM_tray);
#ifdef _OPENMP
				#pragma omp critical(monoB)
#endif
				{
					lattice_result_mono_B.insert(lattice_result_mono_B.end(), latFOM_tray.begin(), latFOM_tray.end());
				}
			}
			if( JudgeSymmetry[Orthorhombic_I] )
			{
				putBravaisLatticeFigureOfMerit(RLCB, latFOM, BravaisType( pair<eBravaisLattice, ePointGroup>(Inner, D2h) ), latFOM_tray);
#ifdef _OPENMP
				#pragma omp critical(orthoI)
#endif
				{
					lattice_result_ortho_I.insert(lattice_result_ortho_I.end(), latFOM_tray.begin(), latFOM_tray.end());
				}
			}
			if( JudgeSymmetry[Orthorhombic_F] )
			{
				putBravaisLatticeFigureOfMerit(RLCB, latFOM, BravaisType( pair<eBravaisLattice, ePointGroup>(Face, D2h) ), latFOM_tray);
#ifdef _OPENMP
				#pragma omp critical(orthoF)
#endif
				{
					lattice_result_ortho_F.insert(lattice_result_ortho_F.end(), latFOM_tray.begin(), latFOM_tray.end());
				}
			}
			if( JudgeSymmetry[Rhombohedral] )
			{
				putBravaisLatticeFigureOfMerit(RLCB, latFOM, BravaisType( put_rhombohedral_type(rh_axis) ), latFOM_tray);
#ifdef _OPENMP
				#pragma omp critical(rhom)
#endif
				{
					lattice_result_rhom.insert(lattice_result_rhom.end(), latFOM_tray.begin(), latFOM_tray.end());
				}
			}

			if( JudgeSymmetry[Monoclinic_P] )
			{
				putLatticeFigureOfMerit(latFOM, put_monoclinic_p_type(abc_axis), m_cv2, latFOM_tray);
#ifdef _OPENMP
				#pragma omp critical(monoP)
#endif
				{
					lattice_result_mono_P.insert(lattice_result_mono_P.end(), latFOM_tray.begin(), latFOM_tray.end());
				}
			}
			if( JudgeSymmetry[Orthorhombic_P] )
			{
				putLatticeFigureOfMerit(latFOM, D2h, m_cv2, latFOM_tray);
#ifdef _OPENMP
				#pragma omp critical (ortho_P)
#endif
				{
					lattice_result_ortho_P.insert(lattice_result_ortho_P.end(), latFOM_tray.begin(), latFOM_tray.end());
				}
			}
	   	}
	}

/*   2011.10.19 VIC Tamura */
CHECK_INTERRUPTION();

ZLOG_INFO( "All the lattice constants are being optimized by linear least squares...\n" );

/*   2011.10.19 VIC Tamura */
LOOP_COUNTER = 0;
Int4 SUM = 0;
for(Int4 i=0; i<NUM_LS; i++) { SUM += lattice_result[i].size(); }

	for(ArrayIndex i=1; i<NUM_LS; i++)
	{
		if( !JudgeSymmetry[i] ) continue;
		sort( lattice_result[i].begin(), lattice_result[i].end() );

		const Int4 num_lattice = lattice_result[i].size();

#ifdef _OPENMP
		#pragma omp parallel
#endif
		{
			vector< VecDat3<Int4> > closest_hkl_tray;
			vector<bool> is_cal_Q_observed_tray;
			vector<LatticeFigureOfMeritToCheckSymmetry> latFOM_tray;

#ifdef _OPENMP
			#pragma omp for
#endif
			for(Int4 index=0; index<num_lattice; index++)
			{
/*   2011.10.19 VIC Tamura */
SET_PROGRESS(99+1*(LOOP_COUNTER++)/SUM, 66, 30); // critical, but works
if(IS_CANSELED()) continue;

				LatticeFigureOfMeritToCheckSymmetry& latFOM0 = lattice_result[i][index];
				latFOM0.setLabel(index+1);

				latFOM0.setFigureOfMerit(m_num_ref_figure_of_merit,
											VCData::putPeakQData(),
											closest_hkl_tray, is_cal_Q_observed_tray);

				LatticeFigureOfMeritZeroShift latFOM2 = latFOM0.putLatticeFigureOfMerit();
				pair<bool, ZErrorMessage> ans = latFOM2.fitLatticeParameterLinear(VCData::putPeakQData(),
															closest_hkl_tray, is_cal_Q_observed_tray, false);
				if( ans.first )
				{
					assert( latFOM0.putLatticeFigureOfMerit().putFiguresOfMerit().putNumberOfReflectionsForFigureOfMerit() > 0 );
					latFOM2.setFigureOfMerit(latFOM0.putLatticeFigureOfMerit().putFiguresOfMerit().putNumberOfReflectionsForFigureOfMerit(),
												VCData::putPeakQData(),
												closest_hkl_tray, is_cal_Q_observed_tray);
					if( LatticeFigureOfMerit::cmpFOMdeWolff( latFOM2, latFOM0.putLatticeFigureOfMerit() ) )
					{
						latFOM0.setLatticeFigureOfMerit(latFOM2);
					}
				}

				const SetOfFigureOfMerit& setFOM = latFOM0.putLatticeFigureOfMerit().putFiguresOfMerit();
				if( setFOM.putFigureOfMeritWu() < MIN_NormM ) continue;
				if( setFOM.putReversedFigureOfMerit() < MIN_RevM ) continue;

				if( eCrystalSystem(i) == Monoclinic_P )
				{
					if( JudgeSymmetry[Hexagonal] )
					{
						putLatticeFigureOfMerit(latFOM0, D6h, m_cv2, latFOM_tray);
#ifdef _OPENMP
						#pragma omp critical (hex)
#endif
						{
							lattice_result_hex.insert(lattice_result_hex.end(), latFOM_tray.begin(), latFOM_tray.end());
						}
					}
				}
				else if( eCrystalSystem(i) == Monoclinic_B )
				{
					if( JudgeSymmetry[Orthorhombic_C] )
					{
						putLatticeFigureOfMerit(latFOM0, D2h, m_cv2, latFOM_tray);
#ifdef _OPENMP
						#pragma omp critical (ortho_B)
#endif
						{
							lattice_result_ortho_B.insert(lattice_result_ortho_B.end(), latFOM_tray.begin(), latFOM_tray.end());
						}
					}
				}
				else if( eCrystalSystem(i) == Orthorhombic_P )
				{
					if( JudgeSymmetry[Tetragonal_P] )
					{
						putLatticeFigureOfMerit(latFOM0, D4h_Z, m_cv2, latFOM_tray);
#ifdef _OPENMP
						#pragma omp critical (tetra_P)
#endif
						{
							lattice_result_tetra_P.insert(lattice_result_tetra_P.end(), latFOM_tray.begin(), latFOM_tray.end());
						}
					}
					if( JudgeSymmetry[Cubic_P] )
					{
						putLatticeFigureOfMerit(latFOM0, Oh, m_cv2, latFOM_tray);
#ifdef _OPENMP
						#pragma omp critical (cubic_P)
#endif
						{
							lattice_result_cubic_P.insert(lattice_result_cubic_P.end(), latFOM_tray.begin(), latFOM_tray.end());
						}
					}
				}
				else if( eCrystalSystem(i) == Orthorhombic_I )
				{
					if( JudgeSymmetry[Tetragonal_I] )
					{
						putLatticeFigureOfMerit(latFOM0, D4h_Z, m_cv2, latFOM_tray);
#ifdef _OPENMP
						#pragma omp critical (tetra_I)
#endif
						{
							lattice_result_tetra_I.insert(lattice_result_tetra_I.end(), latFOM_tray.begin(), latFOM_tray.end());
						}
					}
					if( JudgeSymmetry[Cubic_I] )
					{
						putLatticeFigureOfMerit(latFOM0, Oh, m_cv2, latFOM_tray);
#ifdef _OPENMP
						#pragma omp critical (cubic_I)
#endif
						{
							lattice_result_cubic_I.insert(lattice_result_cubic_I.end(), latFOM_tray.begin(), latFOM_tray.end());
						}
					}
				}
				else if( eCrystalSystem(i) == Orthorhombic_F )
				{
					if( JudgeSymmetry[Cubic_F] )
					{
						putLatticeFigureOfMerit(latFOM0, Oh, m_cv2, latFOM_tray);
#ifdef _OPENMP
						#pragma omp critical (cubic_F)
#endif
						{
							lattice_result_cubic_F.insert(lattice_result_cubic_F.end(), latFOM_tray.begin(), latFOM_tray.end());
						}
					}
				}
			}
		}
/*   2011.10.19 VIC Tamura */
CHECK_INTERRUPTION();

ZLOG_INFO( "(" + num2str( i+1 ) + ") The number of candidates for " + put_cs_name(eCrystalSystem(i), abc_axis)
			+ " : " + num2str<Int4>( lattice_result[i].size() ) + "\n" );
	}
ZLOG_INFO( "\n" );
    }
    catch(bad_alloc& ball)
    {
    	throw nerror(ball, __FILE__, __LINE__, __FUNCTION__);
    }
}


void SortingLattice::putLatticeCandidatesForEachBravaisTypes(const vector<SymMat43_VCData>& S_super,
					const Double& MIN_NormM,
					const Double& MIN_RevM,
					const eABCaxis& abc_axis,
					const eRHaxis& rh_axis,
					vector<LatticeFigureOfMeritToCheckSymmetry> lattice_result[NUM_LS]) const
{
	vector<LatticeFigureOfMeritToCheckSymmetry>& lattice_result_tri = lattice_result[(Int4)Triclinic];
	putLatticeCandidatesForTriclinic(S_super, MIN_NormM, MIN_RevM, lattice_result_tri);

ZLOG_INFO( "Determination of the Bravais type is being carried out...\n(Solutions with " + putLabel(SCWuM) + " less than " + num2str(MIN_NormM)
		+ " or " + putLabel(SCRevM) + " less than " + num2str(MIN_RevM)
		+ " are automatically removed).\n\n" );

ZLOG_INFO( "The program has removed "	+ num2str<Int4>( S_super.size() - lattice_result_tri.size() )
			+ " solutions because their " + putLabel(SCWuM) + " is less than " + num2str(MIN_NormM)
			+ " or their " + putLabel(SCRevM) + " is less than " + num2str(MIN_RevM) + ".\n\n" );
ZLOG_INFO( "Determination of the Bravais type is being carried out...\n" );
	putLatticeCandidatesForEachBravaisTypes(MIN_NormM, MIN_RevM, abc_axis, rh_axis, lattice_result);
}


void SortingLattice::setNumberOfNeighbors(const eABCaxis& baxis_flag,
					vector<LatticeFigureOfMeritToCheckSymmetry> lattice_result[NUM_LS]) const
{

#ifdef _OPENMP
		#pragma omp for
#endif
	for(ArrayIndex i=0; i<NUM_LS; i++)
	{
		if( !OutputSymmetry[i] ) continue;

		stable_sort( lattice_result[i].begin(), lattice_result[i].end() ); // Sort by the unit-cell volume.
		for(vector<LatticeFigureOfMeritToCheckSymmetry>::iterator it=lattice_result[i].begin(); it<lattice_result[i].end(); it++)
		{
			it->putNumberOfLatticesInNeighborhood() = 0;
		}
	}

	const Double coef_lower = 1.0 - sqrt(3.0)*m_resol2;
	const Double coef_upper = 1.0 + sqrt(3.0)*m_resol2;

	Vec_INT index_tray(put_cs_num());

/*   2011.10.19 VIC Tamura */
Int4 SUM=0, LOOP_COUNTER=0;
for(Int4 i=0; i<NUM_LS; i++ ) { SUM += lattice_result[i].size(); }

#ifdef _OPENMP
			#pragma omp for
#endif
	for(ArrayIndex i=0; i<NUM_LS; i++)
	{
		if( !OutputSymmetry[i] ) continue;

		const Int4 num_lattice = lattice_result[i].size();

		for(Int4 index=0; index<num_lattice; index++)
		{
/*   2011.10.19 VIC Tamura */
SET_PROGRESS(100*(LOOP_COUNTER++)/SUM, 97, 1); // critical, but works
if(IS_CANSELED()) continue;

			LatticeFigureOfMeritToCheckSymmetry& latFOM0 = lattice_result[i][index];
			const LatticeFigureOfMerit& latFOM0_prim = latFOM0.putLatticeFigureOfMerit();
			if( latFOM0.putNumberOfLatticesInNeighborhood() < 0 ) continue;

			const Double& detS = latFOM0_prim.putDeterminantOfGramMatrix();
			const Int4 ibegin = distance( lattice_result[i].begin(), lower_bound( lattice_result[i].begin(), lattice_result[i].end(), detS*coef_lower ) );
			const Int4 iend = distance( lattice_result[i].begin(), upper_bound( lattice_result[i].begin(), lattice_result[i].end(), detS*coef_upper ) );

			Int4 count=0;
			if( i == (Int4)Triclinic )
			{
				const ReducedLatticeToCheckEquiv RLCS(m_resol2, latFOM0_prim.putOptimizedForm());
				for(Int4 index2=ibegin; index2<iend; index2++)
				{
					if( index2 == index ) continue;

					LatticeFigureOfMeritToCheckSymmetry& latFOM2 = lattice_result[i][index2];
					const LatticeFigureOfMerit& latFOM2_prim = latFOM2.putLatticeFigureOfMerit();

					// lattice_result_tri[index2] equals trans_mat * RLCB.m_S_super_obtuse * Transpose(trans_mat)
					if( RLCS.equiv( latFOM2_prim.putSellingReducedForm() ) )
					{
						// Compare the figures of merit.
						if( latFOM2.putNumberOfLatticesInNeighborhood() >= 0
							&& LatticeFigureOfMeritToCheckSymmetry::cmpFOMdeWolff( latFOM2, latFOM0 ) )
						{
							count = -1;
							break;
						}
						else
						{
							count++;
							latFOM2.putNumberOfLatticesInNeighborhood() = -1;
						}
					}
				}
			}
			else
			{
				for(Int4 index2=ibegin; index2<iend; index2++)
				{
					if( index2 == index ) continue;

					LatticeFigureOfMeritToCheckSymmetry& latFOM2 = lattice_result[i][index2];
					const LatticeFigureOfMerit& latFOM2_prim = latFOM2.putLatticeFigureOfMerit();

					// *it2 equals trans_mat * RLCS.m_S_super_obtuse * Transpose(trans_mat)
					if( check_equiv_m( latFOM0_prim.putInverseOfMinkowskiReducedForm(), latFOM2_prim.putInverseOfMinkowskiReducedForm(), m_resol2 ) )
					{
						// Compare the figures of merit.
						if( latFOM2.putNumberOfLatticesInNeighborhood() >= 0
							&& LatticeFigureOfMeritToCheckSymmetry::cmpFOMdeWolff( latFOM2, latFOM0 ) )
						{
							count = -1;
							break;
						}
						else
						{
							count++;
							latFOM2.putNumberOfLatticesInNeighborhood() = -1;
						}
					}
				}
			}

			latFOM0.putNumberOfLatticesInNeighborhood() = count;
		}

		Int4& index = index_tray[i];
		for(vector<LatticeFigureOfMeritToCheckSymmetry>::const_iterator it=lattice_result[i].begin(); it<lattice_result[i].end(); it++)
	   	{
			if( it->putNumberOfLatticesInNeighborhood() >= 0 ) index++;
	   	}
	}

/*   2011.10.19 VIC Tamura */
CHECK_INTERRUPTION();

for(ArrayIndex i=0; i<NUM_LS; i++)
{
	if( !OutputSymmetry[i] ) continue;
	ZLOG_INFO( "(" + num2str( i+1 ) + ") The number of candidates for " + put_cs_name(eCrystalSystem(i), baxis_flag)
					+ " : " + num2str( lattice_result[i].size() ) + " -> " + num2str( index_tray[i] ) + "\n" );
}
ZLOG_INFO( "\n" );
}
