/*
 * 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.
 *
 */
#include "../utility_func/chToDouble.hh"
#include "../utility_lattice_reduction/put_Minkowski_reduced_lattice.hh"
#include "OutputInfo.hh"
#include "ReducedLatticeToCheckBravais.hh"
#include "LatticeFigureOfMeritToCheckSymmetry.hh"

const string LatticeFigureOfMeritToCheckSymmetry::CS_LABEL[NUM_LS] =
		{	"01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12", "13", "14" };

LatticeFigureOfMeritToCheckSymmetry::LatticeFigureOfMeritToCheckSymmetry(const Double& rhs)
	: m_label(-1), m_latfom(rhs),
	  m_S_red( SymMat43_VCData( SymMat<VCData>(3), NRMat<Int4>(4,3) ) )
{
}


LatticeFigureOfMeritToCheckSymmetry::LatticeFigureOfMeritToCheckSymmetry(const BravaisType& brat,
		const SymMat43_VCData& S,
		const ePeakShiftFunctionType& type,
		const Double& wave_length,
		const vector<ZParawError>& peak_shift_param_rad)
	: m_label(-1),
	  m_S_red( SymMat43_VCData( SymMat<VCData>(3), NRMat<Int4>(4,3) ) )
{
	this->setLatticeConstants43(brat, S);
	m_latfom.setPeakShiftParamRadian(type, wave_length, peak_shift_param_rad);
	m_num_lattice_found = 0;
}


#ifdef DEBUG
static bool checkInitialLatticeParameters(
		const BravaisType& brat,
		const SymMat43_VCData& S_red)
{
	const SymMat<Double> dbl_S_red( chToDouble(S_red.first) );

	if( brat.enumPointGroup() == Ci && brat.enumBravaisLattice() == Prim )
	{
		assert( dbl_S_red(2,2)*0.9999 < dbl_S_red(1,1) && dbl_S_red(1,1)*0.9999 < dbl_S_red(0,0)
				&& fabs( dbl_S_red(0,1) ) * 1.9999 < dbl_S_red(1,1)
				&& fabs( dbl_S_red(0,2) ) * 1.9999 < dbl_S_red(2,2)
				&& fabs( dbl_S_red(1,2) ) * 1.9999 < dbl_S_red(2,2) );
	}
	else if( brat.enumPointGroup() == C2h_Y && brat.enumBravaisLattice() == Prim )
	{
		assert( 0.0 <= dbl_S_red(0,2) && dbl_S_red(2,2)*0.9999 < dbl_S_red(0,0)
				&& fabs( dbl_S_red(0,2) ) * 1.9999 < dbl_S_red(2,2) && fabs( dbl_S_red(0,2) ) * 1.9999 < dbl_S_red(0,0) );
	}
	else if( brat.enumPointGroup() == C2h_Z && brat.enumBravaisLattice() == Prim )
	{
		assert( 0.0 <= dbl_S_red(0,1) && dbl_S_red(1,1)*0.9999 < dbl_S_red(0,0)
				&& fabs( dbl_S_red(0,1) ) * 1.9999 < dbl_S_red(0,0) && fabs( dbl_S_red(0,1) ) * 1.9999 < dbl_S_red(1,1) );
	}
	else if( brat.enumPointGroup() == C2h_X && brat.enumBravaisLattice() == Prim )
	{
		assert( 0.0 <= dbl_S_red(1,2) && dbl_S_red(2,2)*0.9999 < dbl_S_red(1,1)
				&& fabs( dbl_S_red(1,2) ) * 1.9999 < dbl_S_red(1,1) && fabs( dbl_S_red(1,2) ) * 1.9999 < dbl_S_red(2,2) );
	}
	else if( brat.enumPointGroup() == C2h_Y && brat.enumBravaisLattice() == BaseZ )
	{
		assert( 0.0 <= dbl_S_red(0,2) && fabs( dbl_S_red(0,2) ) * 1.9999 < dbl_S_red(2,2) && fabs( dbl_S_red(0,2) ) * 0.9999 < dbl_S_red(0,0) );
	}
	else if( brat.enumPointGroup() == C2h_Z && brat.enumBravaisLattice() == BaseX )
	{
		assert( 0.0 <= dbl_S_red(0,1) && fabs( dbl_S_red(0,1) ) * 1.9999 < dbl_S_red(0,0) && fabs( dbl_S_red(0,1) ) * 0.9999 < dbl_S_red(1,1) );
	}
	else if( brat.enumPointGroup() == C2h_X && brat.enumBravaisLattice() == BaseY )
	{
		assert( 0.0 <= dbl_S_red(1,2) && fabs( dbl_S_red(1,2) ) * 1.9999 < dbl_S_red(1,1) && fabs( dbl_S_red(1,2) ) * 0.9999  < dbl_S_red(2,2) );
	}
	else if( brat.enumPointGroup() == D2h
			&& brat.enumBravaisLattice() != BaseX
			&& brat.enumBravaisLattice() != BaseY
			&& brat.enumBravaisLattice() != BaseZ )
	{
		assert( dbl_S_red(2,2)*0.9999 < dbl_S_red(1,1) && dbl_S_red(1,1)*0.9999 < dbl_S_red(0,0) );
	}

	const SymMat<VCData> S_super = transform_sym_matrix(S_red.second, S_red.first);
	assert( S_super(0,1) <= VCData()
			&& S_super(0,2) <= VCData()
			&& S_super(0,3) <= VCData()
			&& S_super(1,2) <= VCData()
			&& S_super(1,3) <= VCData()
			&& S_super(2,3) <= VCData() );

	SymMat<VCData> S_red_cp = S_red.first;
	cal_average_crystal_system(brat.enumPointGroup(), S_red_cp);
	assert( S_red.first(0,0).putVecCoef() == S_red_cp(0,0).putVecCoef() );
	assert( S_red.first(1,1).putVecCoef() == S_red_cp(1,1).putVecCoef() );
	assert( S_red.first(2,2).putVecCoef() == S_red_cp(2,2).putVecCoef() );
	assert( S_red.first(0,1).putVecCoef() == S_red_cp(0,1).putVecCoef() );
	assert( S_red.first(0,2).putVecCoef() == S_red_cp(0,2).putVecCoef() );
	assert( S_red.first(1,2).putVecCoef() == S_red_cp(1,2).putVecCoef() );

	return true;
}
#endif


void LatticeFigureOfMeritToCheckSymmetry::setLatticeConstants43(const BravaisType& brat, const SymMat43_VCData& S)
{
	m_S_red = S;
	assert( checkInitialLatticeParameters(brat, m_S_red) );

	m_latfom.setLatticeConstants43(brat, SymMat43_Double(chToDouble(m_S_red.first), m_S_red.second));
	m_num_lattice_found = 0;
//	for(Int4 i=0; i<NUM_LS; i++) m_lattice_equiv[i].clear();
}




bool LatticeFigureOfMeritToCheckSymmetry::checkIfLatticeIsMonoclinic(const ePointGroup& epg_new,
		const Double& cv2,
		map< SymMat<VCData>, NRMat<Int4> >& ans) const
{
	ans.clear();

	SymMat<VCData> ans0 = m_S_red.first;
	cal_average_crystal_system(C2h_X, ans0);

	SymMat<VCData> S_red(3);
	NRMat<Int4> trans_mat2;
	if( check_equiv_m(ans0, m_S_red.first, cv2 ) )
	{
		if( epg_new == C2h_X )
		{
			S_red = ans0;
			trans_mat2 = m_S_red.second;
			putMinkowskiReducedMonoclinicP(1, 2, S_red, trans_mat2);
		}
		else if( epg_new == C2h_Y )
		{
			S_red = transform_sym_matrix(put_matrix_YXZ(), ans0);
			trans_mat2 = mprod(m_S_red.second, put_matrix_YXZ());
			putMinkowskiReducedMonoclinicP(0, 2, S_red, trans_mat2);
		}
		else // if( epg_new == C2h_Z )
		{
			S_red = transform_sym_matrix(put_matrix_YZX(), ans0);
			trans_mat2 = mprod(m_S_red.second, put_matrix_ZXY());
			putMinkowskiReducedMonoclinicP(0, 1, S_red, trans_mat2);
		}
		ans.insert( SymMat43_VCData( S_red, trans_mat2) );
	}

	ans0 = m_S_red.first;
	cal_average_crystal_system(C2h_Y, ans0);
	if( check_equiv_m(ans0, m_S_red.first, cv2 ) )
	{
		if( epg_new == C2h_X )
		{
			S_red = transform_sym_matrix(put_matrix_YXZ(), ans0);
			trans_mat2 = mprod(m_S_red.second, put_matrix_YXZ());
			putMinkowskiReducedMonoclinicP(1, 2, S_red, trans_mat2);
		}
		else if( epg_new == C2h_Y )
		{
			S_red = ans0;
			trans_mat2 = m_S_red.second;
			putMinkowskiReducedMonoclinicP(0, 2, S_red, trans_mat2);
		}
		else // if( epg_new == C2h_Z )
		{
			S_red = transform_sym_matrix(put_matrix_XZY(), ans0);
			trans_mat2 = mprod(m_S_red.second, put_matrix_XZY());
			putMinkowskiReducedMonoclinicP(0, 1, S_red, trans_mat2);
		}
		ans.insert( SymMat43_VCData( S_red, trans_mat2) );
	}

	ans0 = m_S_red.first;
	cal_average_crystal_system(C2h_Z, ans0);
	if( check_equiv_m(ans0, m_S_red.first, cv2 ) )
	{
		if( epg_new == C2h_X )
		{
			S_red = transform_sym_matrix(put_matrix_ZXY(), ans0);
			trans_mat2 = mprod(m_S_red.second, put_matrix_YZX());
			putMinkowskiReducedMonoclinicP(1, 2, S_red, trans_mat2);
		}
		else if( epg_new == C2h_Y )
		{
			S_red = transform_sym_matrix(put_matrix_XZY(), ans0);
			trans_mat2 = mprod(m_S_red.second, put_matrix_XZY());
			putMinkowskiReducedMonoclinicP(0, 2, S_red, trans_mat2);
		}
		else // if( epg_new == C2h_Z )
		{
			S_red = ans0;
			trans_mat2 = m_S_red.second;
			putMinkowskiReducedMonoclinicP(0, 1, S_red, trans_mat2);
		}
		ans.insert( SymMat43_VCData( S_red, trans_mat2) );
	}

	return !( ans.empty() );
}


bool LatticeFigureOfMeritToCheckSymmetry::checkIfLatticeIsOrthorhombic(const Double& cv2,
							map< SymMat<VCData>, NRMat<Int4> >& ans) const
{
	ans.clear();

	const BravaisType& brat = m_latfom.putBravaisType();

	SymMat<VCData> ans0 = m_S_red.first;
	cal_average_crystal_system(D2h, ans0);
	if( check_equiv_m(ans0, m_S_red.first, cv2 ) )
	{
		if( brat.enumBravaisLattice() == BaseX )
		{
			if( ans0(1,1) < ans0(2,2) )
			{
				ans.insert( SymMat43_VCData( transform_sym_matrix(put_matrix_ZYX(), ans0), mprod( m_S_red.second, put_matrix_ZYX() ) ) );
			}
			else
			{
				ans.insert( SymMat43_VCData( transform_sym_matrix(put_matrix_YZX(), ans0), mprod( m_S_red.second, put_matrix_ZXY() ) ) );
			}
		}
		else if( brat.enumBravaisLattice() == BaseY )
		{
			if( ans0(0,0) < ans0(2,2) )
			{
				ans.insert( SymMat43_VCData( transform_sym_matrix(put_matrix_ZXY(), ans0), mprod( m_S_red.second, put_matrix_YZX() ) ) );
			}
			else
			{
				ans.insert( SymMat43_VCData( transform_sym_matrix(put_matrix_XZY(), ans0), mprod( m_S_red.second, put_matrix_XZY() ) ) );
			}
		}
		else if( brat.enumBravaisLattice() == BaseZ )
		{
			if( ans0(0,0) < ans0(1,1) )
			{
				ans.insert( SymMat43_VCData( transform_sym_matrix(put_matrix_YXZ(), ans0), mprod( m_S_red.second, put_matrix_YXZ() ) ) );
			}
			else
			{
				ans.insert( SymMat43_VCData( transform_sym_matrix(put_matrix_XYZ(), ans0), m_S_red.second ) );
			}
		}
		else
		{
			NRMat<Int4> trans_mat = m_S_red.second;
			putMinkowskiReducedOrthorhombic(brat.enumBravaisLattice(), ans0, trans_mat);
			ans.insert( SymMat43_VCData(ans0, trans_mat ) );
		}
		return true;
	}
	return false;
}


bool LatticeFigureOfMeritToCheckSymmetry::checkIfLatticeIsTetragonal(const Double& cv2,
		map< SymMat<VCData>, NRMat<Int4> >& ans) const
{
	ans.clear();

	SymMat<VCData> ans0 = m_S_red.first;
	cal_average_crystal_system(D4h_X, ans0);
	if( check_equiv_m(ans0, m_S_red.first, cv2 ) )
	{
		ans.insert( SymMat43_VCData(
				transform_sym_matrix(put_matrix_YZX(), ans0), mprod( m_S_red.second, put_matrix_ZXY() ) ) );
	}

	ans0 = m_S_red.first;
	cal_average_crystal_system(D4h_Y, ans0);
	if( check_equiv_m(ans0, m_S_red.first, cv2 ) )
	{
		ans.insert( SymMat43_VCData(
				transform_sym_matrix(put_matrix_XZY(), ans0), mprod( m_S_red.second, put_matrix_XZY() ) ) );
	}

	ans0 = m_S_red.first;
	cal_average_crystal_system(D4h_Z, ans0);
	if( check_equiv_m(ans0, m_S_red.first, cv2 ) )
	{
		ans.insert( SymMat43_VCData(ans0, m_S_red.second ) );
	}

	return !( ans.empty() );
}




bool LatticeFigureOfMeritToCheckSymmetry::checkIfLatticeIsHexagonal(const ePointGroup& epg_new, const Double& cv2,
		map< SymMat<VCData>, NRMat<Int4> >& ans) const
{
	ans.clear();
	const BravaisType& brat = m_latfom.putBravaisType();

	SymMat43_VCData ans2(SymMat<VCData>(3), NRMat<Int4>(3,3));

	if( brat.enumPointGroup() == C2h_X )
	{
		ans2.first = transform_sym_matrix(put_matrix_YZX(), m_S_red.first);
		ans2.second = mprod( m_S_red.second, put_matrix_ZXY() );
	}
	else if( brat.enumPointGroup() == C2h_Y )
	{
		ans2.first = transform_sym_matrix(put_matrix_ZXY(), m_S_red.first);
		ans2.second = mprod( m_S_red.second, put_matrix_YZX() );
	}
	else // if( brat.enumPointGroup() == C2h_Z )
	{
		ans2.first = transform_sym_matrix(put_matrix_XYZ(), m_S_red.first);
		ans2.second = m_S_red.second;
	}

	if( ans2.first(0,1) < VCData() )
	{
		ans2.first(0,1) *= -1;
		ans2.second[0][0] *= -1;
		ans2.second[1][0] *= -1;
		ans2.second[2][0] *= -1;
	}

	SymMat<VCData> ans0 = ans2.first;
	cal_average_crystal_system(epg_new, ans2.first);
	if( check_equiv_m(ans2.first, ans0, cv2 ) )
	{
		ans.insert( ans2 );
		return true;
	}
	else return false;
}


bool LatticeFigureOfMeritToCheckSymmetry::checkLatticeSymmetry(const ePointGroup& epg_new, const Double& cv2,
		map< SymMat<VCData>, NRMat<Int4> >& ans) const
{
	ans.clear();
	const BravaisType& brat = m_latfom.putBravaisType();
	if( epg_new == Ci || epg_new == brat.enumPointGroup() )
	{
		ans.insert( m_S_red );
		return true;
	}

	if( epg_new == C2h_X || epg_new == C2h_Y ||  epg_new == C2h_Z )
	{
		assert( brat.enumPointGroup() == Ci );
		assert( brat.enumBravaisLattice() == Prim );

		return checkIfLatticeIsMonoclinic(epg_new, cv2, ans);
	}
	else if( epg_new == D4h_Z )
	{
		assert( brat.enumPointGroup() == D2h );
		assert( brat.enumBravaisLattice() == Prim
				|| brat.enumBravaisLattice() == Inner );

		return checkIfLatticeIsTetragonal(cv2, ans);
	}
	else if( epg_new == D2h )
	{
		assert( brat.enumPointGroup() != Ci || brat.enumBravaisLattice() == Prim );
		assert( brat.enumPointGroup() != C2h_Z || brat.enumBravaisLattice() == BaseX );
		assert( brat.enumPointGroup() != C2h_X || brat.enumBravaisLattice() == BaseY );
		assert( brat.enumPointGroup() != C2h_Y || brat.enumBravaisLattice() == BaseZ );
		assert( brat.enumBravaisLattice() != Rhom_hex );

		return checkIfLatticeIsOrthorhombic(cv2, ans);
	}
	else if( epg_new == D6h )
	{
		assert( brat.enumBravaisLattice() == Prim );
		assert( brat.enumPointGroup() == C2h_X
				|| brat.enumPointGroup() == C2h_Y
				|| brat.enumPointGroup() == C2h_Z );
		return checkIfLatticeIsHexagonal(epg_new, cv2, ans);
	}
	else
	{
		assert( epg_new == Oh );
		assert( brat.enumBravaisLattice() == Prim
				|| brat.enumBravaisLattice() == Inner
				|| brat.enumBravaisLattice() == Face );

		SymMat43_VCData ans2 = m_S_red;
		cal_average_crystal_system(epg_new, ans2.first);
		if( check_equiv_m(ans2.first, m_S_red.first, cv2 ) )
		{
			ans.insert( ans2 );
			return true;
		}
	}
	return !(ans.empty());
}


void LatticeFigureOfMeritToCheckSymmetry::putEquivalentLatticeConstantsDegreeWithOtherCentring(
		const eABCaxis& abc_axis, const eRHaxis& rh_axis,
		vector< pair< eCrystalSystem, SymMat<Double> > >& ans) const
{
	ans.clear();

	static const Double cv2 = 0.04;
	static const Double resol2 = 0.06;

	// Calculate figures of merit as triclinic
	const ReducedLatticeToCheckBravais RLCB(abc_axis, rh_axis, false, cv2, this->putInitialForm());
	const SymMat<Double> S_obtuse = this->putInitialSellingReducedForm();

	if( m_latfom.enumCrystalSystem() != Rhombohedral )
	{
		const map< SymMat<VCData>, NRMat<Int4> >& S_red_tray = RLCB.checkBravaisLatticeType(BravaisType(Rhombohedral, abc_axis, rh_axis));
		for(map< SymMat<VCData>, NRMat<Int4> >::const_iterator it=S_red_tray.begin(); it!=S_red_tray.end(); it++)
		{
			const SymMat<Double> S_red = chToDouble(it->first);
			if( !check_equiv_s(S_obtuse, transform_sym_matrix(it->second, S_red), resol2) ) continue;
			ans.push_back( pair< eCrystalSystem, SymMat<Double> >(Rhombohedral, S_red) );
		}
	}
	if( m_latfom.enumBravaisLattice() != Face )
	{
		const map< SymMat<VCData>, NRMat<Int4> >& S_red_tray = RLCB.checkBravaisLatticeType(BravaisType(Orthorhombic_F, abc_axis, rh_axis));
		for(map< SymMat<VCData>, NRMat<Int4> >::const_iterator it=S_red_tray.begin(); it!=S_red_tray.end(); it++)
		{
			const SymMat<Double> S_red = chToDouble(it->first);
			if( !check_equiv_s(S_obtuse, transform_sym_matrix(it->second, S_red), resol2) ) continue;
			ans.push_back( pair< eCrystalSystem, SymMat<Double> >(Orthorhombic_F, S_red ) );
		}
	}
	if( m_latfom.enumBravaisLattice() != Inner )
	{
		const map< SymMat<VCData>, NRMat<Int4> >& S_red_tray = RLCB.checkBravaisLatticeType(BravaisType(Orthorhombic_I, abc_axis, rh_axis));
		for(map< SymMat<VCData>, NRMat<Int4> >::const_iterator it=S_red_tray.begin(); it!=S_red_tray.end(); it++)
		{
			const SymMat<Double> S_red = chToDouble(it->first);
			if( !check_equiv_s(S_obtuse, transform_sym_matrix(it->second, S_red), resol2) ) continue;
			ans.push_back( pair< eCrystalSystem, SymMat<Double> >(Orthorhombic_I, S_red ) );
		}
	}
	if( m_latfom.enumBravaisLattice() != BaseX && m_latfom.enumBravaisLattice() != BaseY && m_latfom.enumBravaisLattice() != BaseZ )
	{
		const map< SymMat<VCData>, NRMat<Int4> >& S_red_tray = RLCB.checkBravaisLatticeType(BravaisType(Monoclinic_B, abc_axis, rh_axis));
		for(map< SymMat<VCData>, NRMat<Int4> >::const_iterator it=S_red_tray.begin(); it!=S_red_tray.end(); it++)
		{
			const SymMat<Double> S_red = chToDouble(it->first);
			if( !check_equiv_s(S_obtuse, transform_sym_matrix(it->second, S_red), resol2) ) continue;
			ans.push_back( pair< eCrystalSystem, SymMat<Double> >(Monoclinic_B, S_red ) );
		}
	}

	NRMat<Int4> trans_mat;
	SymMat<Double> S_red(3);
	if( m_latfom.enumCrystalSystem() == Rhombohedral || m_latfom.enumBravaisLattice() != Prim )
	{
		const SymMat<Double> S_super = put_sym_matrix_size4to3(this->putInitialSellingReducedForm());
		putTransformMatrixToMinkowskiReduced(Inverse3(S_super), trans_mat);
		transpose_square_matrix(trans_mat);
		ans.push_back( pair< eCrystalSystem, SymMat<Double> >(Triclinic, transform_sym_matrix(Inverse3(trans_mat), S_super) ) );
	}
}


void LatticeFigureOfMeritToCheckSymmetry::putEquivalentLatticeConstantsDegreeWithOtherCentring(
		const eABCaxis& abc_axis, const eRHaxis& rh_axis,
		vector< pair< eCrystalSystem, pair< VecDat3<Double>, VecDat3<Double> > > >& ans) const
{
	vector< pair< eCrystalSystem, SymMat<Double> > > ans0;
	putEquivalentLatticeConstantsDegreeWithOtherCentring(abc_axis, rh_axis, ans0);

	ans.clear();
	ans.resize( ans0.size() );
	vector< pair< eCrystalSystem, pair< VecDat3<Double>, VecDat3<Double> > > >::iterator it2 = ans.begin();
	for(vector< pair< eCrystalSystem, SymMat<Double> > >::const_iterator it=ans0.begin(); it<ans0.end(); it++, it2++)
	{
		it2->first = it->first;
		LatticeFigureOfMerit::putLatticeConstantsDegree( BravaisType(it->first, abc_axis, rh_axis), it->second, abc_axis, rh_axis, it2->second.first, it2->second.second );
	}
}


void LatticeFigureOfMeritToCheckSymmetry::printLatticeInformation(
					const vector<LatticeFigureOfMeritToCheckSymmetry> lattice_result[],
					const OutputInfo outinfo[],
					const eABCaxis& abc_axis,
					const eRHaxis& rh_axis,
					const Int4& label_start0,
					ostream* os) const
{
	m_latfom.printLatticeInformation(abc_axis, rh_axis, label_start0, os);

//	const FracMat mat_sell_to_min = FInverse3( put_transform_matrix_row4to3( putLatticeFigureOfMerit().putOptimizedFormWithTransformMatrixToSellingReduced().second) );
//  	const SymMat<Double> dbl_S_init = transform_sym_matrix(mat_sell_to_min.mat,
//  										put_sym_matrix_size4to3( this->putInitialSellingReducedForm() ) ) / (mat_sell_to_min.denom*mat_sell_to_min.denom);
//  	const SymMat<Double>& dbl_S_init = this->putInitialSellingReducedForm();

//	os->width(label_start); *os << "";
//  	*os << "<!-- A*, B*, C*, D*, E*, F*(angstrom^(-2)) first given by peak-positions.-->\n";
//	os->width(label_start); *os << "";
//  	*os << "<InitialReciprocalLatticeParameters>";
//  	os->width(14);
//   	*os << dbl_S_init(0,0);
//  	os->width(14);
//   	*os << dbl_S_init(1,1);
//  	os->width(14);
//   	*os << dbl_S_init(2,2);
//  	os->width(14);
//   	*os << dbl_S_init(1,2);
//  	os->width(14);
//   	*os << dbl_S_init(0,2);
//  	os->width(14);
//   	*os << dbl_S_init(0,1);
//  	*os << " </InitialReciprocalLatticeParameters>\n";

  	Int4 label_start = label_start0;
	os->width(label_start);
   	*os << "" << "<NumberOfLatticesInNeighborhood>";
  	os->width(14);
   	*os << this->putNumberOfLatticesInNeighborhood();
   	*os << " </NumberOfLatticesInNeighborhood>\n\n";

	os->width(label_start);
	*os << "" << "<EquivalentLatticeCandidates>\n";
	label_start++;

	vector< pair< eCrystalSystem, pair< VecDat3<Double>, VecDat3<Double> > > > lattice_equiv;
	this->putEquivalentLatticeConstantsDegreeWithOtherCentring(abc_axis, rh_axis, lattice_equiv);

		for(vector< pair< eCrystalSystem, pair< VecDat3<Double>, VecDat3<Double> > > >::const_iterator it=lattice_equiv.begin(); it<lattice_equiv.end(); it++)
		{
			os->width(label_start); *os << "";
	      	*os << "<LatticeCandidate>\n";
			label_start++;

			os->width(label_start); *os << "";
	      	*os << "<CrystalSystem>";
			os->width(17);
			*os << put_cs_name(it->first, abc_axis);
	      	*os << " </CrystalSystem>\n";

			os->width(label_start);
	      	*os << "" << "<LatticeParameters>";
	      	os->width(14);
	      	*os << it->second.first[0];
	      	os->width(14);
	   	   	*os << it->second.first[1];
	      	os->width(14);
	   	   	*os << it->second.first[2];
	      	os->width(14);
	   	   	*os << it->second.second[0];
	      	os->width(14);
	   	   	*os << it->second.second[1];
	      	os->width(14);
	   	   	*os << it->second.second[2];
	      	*os << " </LatticeParameters>\n";

	   	   	label_start--;
			os->width(label_start); *os << "";
	      	*os << "</LatticeCandidate>\n";
		}

	label_start--;
	os->width(label_start); *os << "";
	*os << "</EquivalentLatticeCandidates>\n\n";
}
