67 proshade_double ret = 0.0;
68 proshade_unsign vecSize =
static_cast< proshade_unsign
> ( inArr.size() );
69 proshade_double* meadianAndIQR =
new proshade_double[2];
72 if ( vecSize == 0 ) {
delete[] meadianAndIQR;
return ( ret ); }
73 if ( vecSize <= 4 ) { ret = std::accumulate ( inArr.begin(), inArr.end(), 0.0 ) /
static_cast< proshade_double
> ( vecSize ); }
85 ret = meadianAndIQR[0] + ( meadianAndIQR[1] * noIQRsFromMedian );
89 if ( ret > *( std::max_element ( inArr.begin(), inArr.end() ) ) )
91 ret = *( std::max_element ( inArr.begin(), inArr.end() ) );
95 delete[] meadianAndIQR;
119 proshade_double shellSpacing = ( 2.0 * M_PI ) /
static_cast<proshade_double
> ( this->maxShellBand * 2.0 );
120 std::vector< proshade_double > allPeakHeights;
123 for ( proshade_unsign spIt = 1; spIt < ( this->maxShellBand * 2 ); spIt++ )
127 this->maxShellBand * 2.0,
128 static_cast<proshade_double
> ( spIt ) * shellSpacing,
133 for ( proshade_unsign shIt = 0; shIt < static_cast<proshade_unsign> ( sphereMappedRotFun.size() ); shIt++ )
136 std::stringstream hlpSS;
137 hlpSS <<
"Interpolating sphere " << shIt <<
" ( radius: " << this->sphereMappedRotFun.at(shIt)->getRadius() <<
" ).";
141 this->sphereMappedRotFun.at(shIt)->interpolateSphereValues ( this->getInvSO3Coeffs ( ) );
149 for ( proshade_unsign shIt = 0; shIt < static_cast<proshade_unsign> ( this->sphereMappedRotFun.size() ); shIt++ )
151 this->sphereMappedRotFun.at(shIt)->findAllPeaks ( settings->
peakNeighbours, &allPeakHeights );
155 std::stringstream hlpSS;
156 hlpSS <<
"Detected " << allPeakHeights.size() <<
" peaks with any height.";
163 std::stringstream hlpSS2;
164 hlpSS2 <<
"From these peaks, decided the threshold will be " << peakThres <<
" peak height.";
168 for ( proshade_unsign shIt = 0; shIt < static_cast<proshade_unsign> ( this->sphereMappedRotFun.size() ); shIt++ )
170 this->sphereMappedRotFun.at(shIt)->removeSmallPeaks ( peakThres );
194 std::vector< proshade_double* > ret;
211 std::vector< std::vector< proshade_unsign > > detectedCSymmetries;
212 for ( proshade_signed iter =
static_cast<proshade_unsign
> ( peakGroupsBoundaries.size() - 1 ); iter >= 0; iter-- )
215 std::vector< proshade_double* > symPeaks;
216 for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( peaksAA.size() ); it++ )
226 settings->axisErrToleranceDefault,
237 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( peaksAA.size() ); iter++ ) {
delete[] allPeaks.at(iter);
delete[] peaksAA.at(iter); }
259 std::vector< proshade_double* > ret;
260 proshade_double* hlpP = NULL;
261 proshade_double* rotMat =
new proshade_double [9];
265 for ( proshade_unsign peakIter = 0; peakIter < static_cast<proshade_unsign> ( allPeaks.size() ); peakIter++ )
271 hlpP =
new proshade_double [5];
276 hlpP[4] = allPeaks.at(peakIter)[3];
286 std::stringstream hlpSSP;
287 hlpSSP <<
"Found " << ret.size() <<
" possible peaks.";
291 if ( ret.size() < 1 )
293 ProSHADE_internal_messages::printWarningMessage ( verbose,
"!!! ProSHADE WARNING !!! Failed to detect any symmetries. There are no reasonable peaks in the self-rotation map. If you believe there should be some symmetry, you can try decreasing the resolution or changing the peak IQR threshold.",
"WS00029" );
313 std::vector< proshade_double > boundaries;
315 proshade_double peakContribution = 0.0;
318 std::vector< proshade_double > pdf;
319 for ( proshade_double iter = 0.0; iter <= 1.0; iter += 0.01 )
322 peakContribution = 0.0;
325 for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( allPeaks.size() ); it++ )
335 proshade_double prev = pdf.at(0);
336 for ( proshade_unsign iter = 1; iter < static_cast<proshade_unsign> ( pdf.size() - 1 ); iter ++ )
339 if ( ( prev > pdf.at(iter) ) && ( pdf.at(iter+1) > pdf.at(iter) ) )
349 return ( boundaries );
371 std::vector< std::vector< proshade_unsign > > ret;
372 std::vector< proshade_double > triedAlready;
373 std::vector< proshade_unsign > angsToTry, testedAlready;
374 proshade_double angDist, angDivisionRemainder, angDivisionBasis, nextSymmetryError, nextPeakError = ( M_PI * 2.0 ) / (
static_cast<proshade_double
> ( band ) * 2.0 );
377 if ( peaks->size() < 1 ) {
return ( ret ); }
383 for ( proshade_unsign grpIt = 0; grpIt < static_cast<proshade_unsign> ( sameAxesGroups.size() ); grpIt++ )
389 triedAlready.clear ( );
390 testedAlready.clear ( );
422 std::vector< std::vector< proshade_unsign > > ret;
423 bool sameAxisFound =
false;
424 proshade_double angTolerance = std::acos ( 1.0 - errTolerance );
430 for ( proshade_unsign peakIter = 0; peakIter < static_cast<proshade_unsign> ( peaks.size() ); peakIter++ )
433 sameAxisFound =
false;
436 if ( ( peaks.at(peakIter)[3] - angTolerance <= 0.0 ) && ( peaks.at(peakIter)[3] + angTolerance > 0.0 ) ) {
continue; }
440 if ( ( ( peaks.at(peakIter)[0] - 0.1 <= 0.0 ) && ( peaks.at(peakIter)[0] + 0.1 > 0.0 ) ) &&
441 ( ( peaks.at(peakIter)[1] - 0.1 <= 0.0 ) && ( peaks.at(peakIter)[1] + 0.1 > 0.0 ) ) &&
442 ( ( peaks.at(peakIter)[2] - 0.1 <= 0.0 ) && ( peaks.at(peakIter)[2] + 0.1 > 0.0 ) ) ) {
continue; }
445 for ( proshade_unsign sameAxisGrp = 0; sameAxisGrp < static_cast<proshade_unsign> ( ret.size() ); sameAxisGrp++ )
448 for ( proshade_unsign sameAxis = 0; sameAxis < static_cast<proshade_unsign> ( ret.at(sameAxisGrp).size() ); sameAxis++ )
452 peaks.at(ret.at(sameAxisGrp).at(sameAxis))[1],
453 peaks.at(ret.at(sameAxisGrp).at(sameAxis))[2],
454 peaks.at(peakIter)[0],
455 peaks.at(peakIter)[1],
456 peaks.at(peakIter)[2],
459 sameAxisFound =
true;
467 if ( sameAxisFound ) {
continue; }
470 std::vector<proshade_unsign> hlpVec;
494 for ( proshade_unsign i = 0; i < static_cast<proshade_unsign> ( peaks.size() ); i++ )
496 if ( ( ( std::max ( std::abs ( peaks.at(i)[0] ), std::max( std::abs ( peaks.at(i)[1] ), std::abs ( peaks.at(i)[2] ) ) ) == std::abs ( peaks.at(i)[0] ) ) && ( peaks.at(i)[0] < 0.0 ) ) ||
497 ( ( std::max ( std::abs ( peaks.at(i)[0] ), std::max( std::abs ( peaks.at(i)[1] ), std::abs ( peaks.at(i)[2] ) ) ) == std::abs ( peaks.at(i)[1] ) ) && ( peaks.at(i)[1] < 0.0 ) ) ||
498 ( ( std::max ( std::abs ( peaks.at(i)[0] ), std::max( std::abs ( peaks.at(i)[1] ), std::abs ( peaks.at(i)[2] ) ) ) == std::abs ( peaks.at(i)[2] ) ) && ( peaks.at(i)[2] < 0.0 ) ) )
500 peaks.at(i)[0] *= -1.0;
501 peaks.at(i)[1] *= -1.0;
502 peaks.at(i)[2] *= -1.0;
503 peaks.at(i)[3] *= -1.0;
521 std::stringstream hlpSS;
522 hlpSS <<
"Symmetry axis group " << groupNo;
527 for ( proshade_unsign axIt = 0; axIt < static_cast<proshade_unsign> ( grp.size() ); axIt++ )
529 std::stringstream SS;
530 SS <<
" " << axIt <<
"\t " <<
static_cast<int>( peaks.at(grp.at(axIt))[0] * 100.0 ) / 100.0 <<
"\t" <<
static_cast<int>( peaks.at(grp.at(axIt))[1] * 100.0 ) / 100.0 <<
"\t" <<
static_cast<int>( peaks.at(grp.at(axIt))[2] * 100.0 ) / 100.0 <<
"\t" <<
static_cast<int>( peaks.at(grp.at(axIt))[3] * 100.0 ) / 100.0 <<
"\t" <<
static_cast<int>( peaks.at(grp.at(axIt))[4] * 100.0 ) / 100.0;
556 proshade_unsign g1 = 0, g2 = 0;
560 for ( proshade_unsign gr1It = 0; gr1It < static_cast<proshade_unsign> ( grp.size() ); gr1It++ )
562 for ( proshade_unsign gr2It = 1; gr2It < static_cast<proshade_unsign> ( grp.size() ); gr2It++ )
565 if ( gr1It >= gr2It ) {
continue; }
569 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( tried->size() ); iter += 3 )
572 if ( ( gr2It == tried->at( iter + 1 ) ) && ( gr1It == tried->at( iter ) ) ) { skip =
true; }
576 ( ( std::abs( std::abs ( peaks.at(grp.at(gr1It))[3] ) - std::abs ( peaks.at(grp.at(gr2It))[3] ) ) - 0.01 ) < tried->at( iter + 2 ) ) &&
577 ( ( std::abs( std::abs ( peaks.at(grp.at(gr1It))[3] ) - std::abs ( peaks.at(grp.at(gr2It))[3] ) ) + 0.01 ) > tried->at( iter + 2 ) ) )
582 if ( skip ) {
continue; }
585 if ( std::abs( std::abs ( peaks.at(grp.at(gr1It))[3] ) - std::abs ( peaks.at(grp.at(gr2It))[3] ) ) < (*dist) )
588 if ( std::abs( std::abs ( peaks.at(grp.at(gr1It))[3] ) - std::abs ( peaks.at(grp.at(gr2It))[3] ) ) > 0.01 )
592 *dist = std::abs( std::abs ( peaks.at(grp.at(gr1It))[3] ) - std::abs ( peaks.at(grp.at(gr2It))[3] ) );
599 if ( *dist != 999.9 )
624 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( grpsVec.size() ); iter++ )
626 proshade_double* hlpP =
new proshade_double [5];
628 hlpP[0] = peaks.at(grpsVec.at(iter).at(0))[0];
629 hlpP[1] = peaks.at(grpsVec.at(iter).at(0))[1];
630 hlpP[2] = peaks.at(grpsVec.at(iter).at(0))[2];
632 hlpP[4] = peaks.at(grpsVec.at(iter).at(0))[4];
664 *divRem = std::modf (
static_cast<proshade_double
> ( ( 2.0 * M_PI ) / std::abs ( dist ) ), divBasis );
674 *symmErr = ( M_PI * 2.0 / *divBasis ) - ( M_PI * 2.0 / ( *divBasis + 1.0 ) );
675 proshade_double angTolerance = ( peakErr / *symmErr );
678 if ( ( *divRem < ( 0.0 + angTolerance ) ) && ( *divRem > ( 0.0 - angTolerance ) ) )
682 for ( proshade_signed iter = -angTolRound; iter <= angTolRound; iter++ )
689 if ( angsToTry->size() == 0 ) { ret =
false; }
708 proshade_double groupAngle = ( 2.0 * M_PI ) /
static_cast<proshade_double
> ( fold );
711 for ( proshade_signed iter =
static_cast<proshade_signed
> ( -(
static_cast<proshade_double
> ( fold ) / 2.0 + 1.0) ); iter <= static_cast<proshade_signed> (
static_cast<proshade_double
> ( fold )/2.0 + 1.0 ); iter++ )
734 proshade_unsign
ProSHADE_internal_symmetry::checkExpectedAgainstFound ( std::vector< proshade_unsign > grp, std::vector< proshade_double* > peaks, std::vector< proshade_double >* expAngs, std::vector< proshade_unsign >* matchedAngs, std::vector< proshade_unsign >* missingAngs, proshade_double angTol )
737 proshade_unsign ret = 0;
738 proshade_unsign retHlp = 0;
739 proshade_double groupAngle = expAngs->at(1) - expAngs->at(0);
740 bool matchedThisPeak =
false;
741 bool noDoubleMatches =
false;
742 std::vector < proshade_unsign > matchedAlready;
745 for ( proshade_unsign expAngIt = 0; expAngIt < static_cast<proshade_unsign> ( expAngs->size() ); expAngIt++ )
748 matchedThisPeak =
false;
749 for ( proshade_unsign peakIt = 0; peakIt < static_cast<proshade_unsign> ( grp.size() ); peakIt++ )
751 if ( ( expAngs->at(expAngIt) < ( peaks.at(grp.at(peakIt))[3] + angTol ) ) &&
752 ( expAngs->at(expAngIt) > ( peaks.at(grp.at(peakIt))[3] - angTol ) ) )
754 noDoubleMatches =
false;
755 for ( proshade_unsign ndm = 0; ndm < static_cast<proshade_unsign> ( matchedAlready.size() ); ndm++ )
757 if ( matchedAlready.at(ndm) == grp.at(peakIt) ) { noDoubleMatches =
true;
break; }
760 if ( !noDoubleMatches )
764 matchedThisPeak =
true;
771 if ( !matchedThisPeak )
778 if ( matchedAngs->size () > 1 )
780 for ( proshade_unsign iter = 1; iter < static_cast<unsigned int> ( matchedAngs->size () ); iter++ )
782 if ( ( ( peaks.at(matchedAngs->at(iter-1))[3] + groupAngle ) < ( peaks.at(matchedAngs->at(iter))[3] + angTol ) ) &&
783 ( ( peaks.at(matchedAngs->at(iter-1))[3] + groupAngle ) > ( peaks.at(matchedAngs->at(iter))[3] - angTol ) ) )
791 if ( retHlp > ret ) { ret = retHlp; }
819 proshade_double ret = 0.0;
820 proshade_unsign arrIndex = 0;
821 proshade_double* rotMat =
new proshade_double [9];
823 proshade_double pointHeight, euA, euB, euG, xPk, yPk, zPk, anglPk;
824 proshade_double angTol = std::acos ( 1.0 - axTol );
827 for ( proshade_unsign xIt = 0; xIt < ( dataObj->
getMaxBand() * 2 ); xIt++ )
829 for ( proshade_unsign yIt = 0; yIt < ( dataObj->
getMaxBand() * 2 ); yIt++ )
831 for ( proshade_unsign zIt = 0; zIt < ( dataObj->
getMaxBand() * 2 ); zIt++ )
836 if ( pointHeight < heightThres ) {
continue; }
840 static_cast<proshade_signed
> ( yIt ),
static_cast<proshade_signed
> ( zIt ),
846 if ( ( ( std::abs( anglPk ) - angTol ) < std::abs ( angle ) ) && ( ( std::abs( anglPk ) + angTol ) > std::abs ( angle ) ) )
849 if ( ( ( std::max( std::abs( xPk ), std::max( std::abs( yPk ), std::abs( zPk ) ) ) == std::abs( xPk ) ) && ( xPk < 0.0 ) ) ||
850 ( ( std::max( std::abs( xPk ), std::max( std::abs( yPk ), std::abs( zPk ) ) ) == std::abs( yPk ) ) && ( yPk < 0.0 ) ) ||
851 ( ( std::max( std::abs( xPk ), std::max( std::abs( yPk ), std::abs( zPk ) ) ) == std::abs( zPk ) ) && ( zPk < 0.0 ) ) )
862 if ( ret < pointHeight ) { ret = pointHeight; }
887 std::vector< proshade_unsign > hlpVec;
891 for ( proshade_unsign pIt = 0; pIt < static_cast<proshade_unsign> ( matchedPeaks->size() ); pIt++ )
898 std::stringstream hlpS;
899 hlpS <<
"Found symmetry C" << fold;
930 std::stringstream hlpSSP;
931 hlpSSP <<
"Searching for missing peaks for symmetry C" << fold;
935 proshade_double heightThreshold = 0.0;
936 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( grp->size() ); grIt++ ) { heightThreshold += peaks->at(grp->at(grIt))[4]; }
937 heightThreshold /=
static_cast<proshade_double
> ( grp->size() );
938 heightThreshold *= 0.5;
941 for ( proshade_unsign misPkIt = 0; misPkIt < static_cast<proshade_unsign> ( missingPeaks->size() ); misPkIt++ )
944 if ( expectedAngles->at(missingPeaks->at(misPkIt)) > M_PI ) {
continue; }
945 if ( expectedAngles->at(missingPeaks->at(misPkIt)) < -M_PI ) {
continue; }
948 proshade_double misHeight =
ProSHADE_internal_symmetry::checkForMissingPeak ( dataObj, peaks->at(grp->at(0))[0], peaks->at(grp->at(0))[1], peaks->at(grp->at(0))[2], expectedAngles->at(missingPeaks->at(misPkIt)), heightThreshold, axErrTolerance );
949 if ( misHeight != 0.0 )
952 proshade_double* hlpP =
new proshade_double [5];
954 hlpP[0] = peaks->at(grp->at(0))[0];
955 hlpP[1] = peaks->at(grp->at(0))[1];
956 hlpP[2] = peaks->at(grp->at(0))[2];
957 hlpP[3] = expectedAngles->at(missingPeaks->at(misPkIt));
992 bool skipFold =
false;
993 std::vector< proshade_unsign > matchedPeaks, missingPeaks;
994 std::vector< proshade_double > expectedAngles;
995 proshade_double angTolerance = std::acos ( 1.0 - axErrTolerance );
998 for ( proshade_unsign fIt = 0; fIt < static_cast<proshade_unsign> ( angsToTry->size() ); fIt++ )
1002 for ( proshade_unsign ftIt = 0; ftIt < static_cast<proshade_unsign> ( testedAlready->size() ); ftIt++ ) {
if ( testedAlready->at(ftIt) == angsToTry->at(fIt) ) { skipFold =
true; } }
1003 if ( skipFold ) {
continue; }
1007 if ( axErrToleranceDefault )
1009 angTolerance = std::max ( std::min ( angTolerance, ( ( (M_PI * 2.0) /
static_cast<double> ( angsToTry->at(fIt) ) ) -
1010 ( (M_PI * 2.0) /
static_cast<double> ( angsToTry->at(fIt) + 1 ) ) ) * 2.0 ), 0.02 );
1011 axErrTolerance = std::max ( 1.0 - std::cos ( angTolerance ), 0.0008 );
1015 expectedAngles.clear ( );
1019 matchedPeaks.clear ( );
1020 missingPeaks.clear ( );
1022 &matchedPeaks, &missingPeaks, angTolerance );
1025 if ( consecMatches >= angsToTry->at(fIt) )
1031 if ( (
static_cast<proshade_double
> ( matchedPeaks.size() ) /
static_cast<proshade_double
> ( angsToTry->at(fIt) ) ) >= ( 1.0 - missPeakThres ) )
1035 &expectedAngles, &matchedPeaks, axErrTolerance, verbose ) )
1062 std::stringstream ss;
1063 ss <<
"Detected C" << grp.at(0) <<
" symmetry with following peaks:";
1068 for ( proshade_unsign pkIt = 1; pkIt < static_cast<proshade_unsign> ( grp.size() ); pkIt++ )
1070 std::stringstream SS;
1071 SS <<
" " <<
static_cast<int>( peaks.at(grp.at(pkIt))[0] * 100.0 ) / 100.0 <<
"\t" <<
static_cast<int>( peaks.at(grp.at(pkIt))[1] * 100.0 ) / 100.0 <<
"\t" <<
static_cast<int>( peaks.at(grp.at(pkIt))[2] * 100.0 ) / 100.0 <<
"\t" <<
static_cast<int>( peaks.at(grp.at(pkIt))[3] * 100.0 ) / 100.0 <<
"\t" <<
static_cast<int>( peaks.at(grp.at(pkIt))[4] * 100.0 ) / 100.0;
1088 std::stringstream ss;
1089 ss <<
"Detected " << noSyms <<
" Cyclic symmetries.";
1095 ProSHADE_internal_messages::printWarningMessage ( verbose,
"!!! ProSHADE WARNING !!! Failed to detect any symmetries. If you believe there should be one, you can try decreasing the resolution or checking that the map is centred on the centry of symmetry (or use map centering option in ProSHADE).",
"WS00030" );
1118 proshade_double sumX, sumY, sumZ, sumH;
1121 for ( proshade_unsign symIt = 0; symIt < static_cast<proshade_unsign> ( detected.size() ); symIt++ )
1124 proshade_double* hlpP =
new proshade_double [6];
1128 hlpP[0] =
static_cast<proshade_double
> ( detected.at(symIt).at(0) );
1129 hlpP[4] =
static_cast<proshade_double
> ( ( 2.0 * M_PI ) / hlpP[0] );
1132 sumX = 0.0; sumY = 0.0; sumZ = 0.0; sumH = 0.0;
1133 for ( proshade_unsign pkIt = 1; pkIt < static_cast<proshade_unsign> ( detected.at(symIt).size() ); pkIt++ )
1135 sumX += peaks.at(detected.at(symIt).at(pkIt))[0];
1136 sumY += peaks.at(detected.at(symIt).at(pkIt))[1];
1137 sumZ += peaks.at(detected.at(symIt).at(pkIt))[2];
1138 sumH += peaks.at(detected.at(symIt).at(pkIt))[4];
1140 sumX /=
static_cast<proshade_double
> ( detected.at(symIt).size() - 1 );
1141 sumY /=
static_cast<proshade_double
> ( detected.at(symIt).size() - 1 );
1142 sumZ /=
static_cast<proshade_double
> ( detected.at(symIt).size() - 1 );
1143 sumH /=
static_cast<proshade_double
> ( detected.at(symIt).size() - 1 );
1180 proshade_double dotProduct = 0.0;
1183 for ( proshade_unsign symIt = 0; symIt < static_cast<proshade_unsign> ( ret->size() ); symIt++ )
1186 if ( ret->at(symIt)[0] == sym[0] )
1190 &ret->at(symIt)[3], &sym[1], &sym[2], &sym[3] );
1191 if ( ( ( 1.0 > ( dotProduct - simThres ) ) && ( 1.0 < ( dotProduct + simThres ) ) ) || ( ( -1.0 > ( dotProduct - simThres ) ) && ( -1.0 < ( dotProduct + simThres ) ) ) )
1194 if ( ret->at(symIt)[5] >= sym[5] ) {
return (
true ); }
1197 ret->at(symIt)[1] = sym[1];
1198 ret->at(symIt)[2] = sym[2];
1199 ret->at(symIt)[3] = sym[3];
1200 ret->at(symIt)[5] = sym[5];
1225 std::vector< proshade_double* > ret;
1226 proshade_double dotProduct;
1232 if ( CSymList->size() < 2 ) {
return ( ret ); }
1235 for ( proshade_unsign ax1 = 0; ax1 < static_cast<proshade_unsign> ( CSymList->size() ); ax1++ )
1238 if ( CSymList->at(ax1)[5] < settings->
minSymPeak ) {
continue; }
1240 for ( proshade_unsign ax2 = 1; ax2 < static_cast<proshade_unsign> ( CSymList->size() ); ax2++ )
1243 if ( ax1 >= ax2 ) {
continue; }
1246 if ( CSymList->at(ax2)[5] < settings->
minSymPeak ) {
continue; }
1250 &CSymList->at(ax1)[3], &CSymList->at(ax2)[1],
1251 &CSymList->at(ax2)[2], &CSymList->at(ax2)[3] );
1257 if ( CSymList->at(ax1)[0] >= CSymList->at(ax2)[0] )
1261 std::vector< proshade_unsign > DSymInd;
1271 std::vector< proshade_unsign > DSymInd;
1281 std::stringstream hlpSS;
1282 hlpSS <<
"Detected " << ret.size() <<
" D symmetries.";
1304 proshade_double* hlpP =
new proshade_double [12];
1308 hlpP[0] =
static_cast<proshade_double
> ( CSymList->at(axisOne)[0] );
1309 hlpP[4] =
static_cast<proshade_double
> ( ( 2.0 * M_PI ) / hlpP[0] );
1310 hlpP[6] =
static_cast<proshade_double
> ( CSymList->at(axisTwo)[0] );
1311 hlpP[10] =
static_cast<proshade_double
> ( ( 2.0 * M_PI ) / hlpP[6] );
1314 hlpP[1] = CSymList->at(axisOne)[1];
1315 hlpP[2] = CSymList->at(axisOne)[2];
1316 hlpP[3] = CSymList->at(axisOne)[3];
1317 hlpP[5] = CSymList->at(axisOne)[5];
1318 hlpP[7] = CSymList->at(axisTwo)[1];
1319 hlpP[8] = CSymList->at(axisTwo)[2];
1320 hlpP[9] = CSymList->at(axisTwo)[3];
1321 hlpP[11] = CSymList->at(axisTwo)[5];
1343 std::vector< proshade_double* > ret;
1359 for ( proshade_unsign csIt = 0; csIt < static_cast<proshade_unsign> ( CSymList->size() ); csIt++ )
1361 for ( proshade_unsign retIt = 0; retIt < static_cast<proshade_unsign> ( ret.size() ); retIt++ )
1363 if ( ( CSymList->at(csIt)[0] == ret.at(retIt)[0] ) &&
1364 ( CSymList->at(csIt)[1] == ret.at(retIt)[1] ) &&
1365 ( CSymList->at(csIt)[2] == ret.at(retIt)[2] ) &&
1366 ( CSymList->at(csIt)[3] == ret.at(retIt)[3] ) &&
1367 ( CSymList->at(csIt)[4] == ret.at(retIt)[4] ) &&
1368 ( CSymList->at(csIt)[5] == ret.at(retIt)[5] ) )
1399 std::vector< proshade_unsign > C3List;
1400 proshade_double dotProduct;
1403 for ( proshade_unsign cSym = 0; cSym < static_cast<proshade_unsign> ( CSymList->size() ); cSym++ )
1409 for ( proshade_unsign c31 = 0; c31 < static_cast<proshade_unsign> ( C3List.size() ); c31++ )
1411 for ( proshade_unsign c32 = 1; c32 < static_cast<proshade_unsign> ( C3List.size() ); c32++ )
1414 if ( c31 >= c32 ) {
continue; }
1417 dotProduct =
ProSHADE_internal_maths::computeDotProduct ( &CSymList->at(C3List.at(c31))[1], &CSymList->at(C3List.at(c31))[2], &CSymList->at(C3List.at(c31))[3], &CSymList->at(C3List.at(c32))[1], &CSymList->at(C3List.at(c32))[2], &CSymList->at(C3List.at(c32))[3] );
1420 if ( ( ( 1.0 / 3.0 ) > ( dotProduct - axErr ) ) && ( ( 1.0 / 3.0 ) < ( dotProduct + axErr ) ) )
1448 std::vector< proshade_unsign > C3PossibilitiesHlp;
1449 std::vector< std::vector< proshade_unsign > > C3Possibilities;
1456 for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( CSymList->size() ); cIt++ )
1459 if ( CSymList->at(cIt)[0] != 3.0 || CSymList->at(cIt)[0] < minPeakHeight ) {
continue; }
1465 groupMatched =
false;
1466 for ( proshade_unsign gIt = 0; gIt < static_cast<proshade_unsign> ( C3Possibilities.size() ); gIt++ )
1479 proshade_double maxHeight = 0.0; proshade_unsign maxGrp = 0;
1480 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( C3Possibilities.size() ); iter++ ) {
if ( C3Possibilities.at(iter).size() == 4 ) {
if ( ( ( CSymList->at(C3Possibilities.at(iter).at(0))[5] + CSymList->at(C3Possibilities.at(iter).at(1))[5] + CSymList->at(C3Possibilities.at(iter).at(2))[5] + CSymList->at(C3Possibilities.at(iter).at(3))[5] ) / 4.0 ) > maxHeight ) { maxHeight = ( ( CSymList->at(C3Possibilities.at(iter).at(0))[5] + CSymList->at(C3Possibilities.at(iter).at(1))[5] + CSymList->at(C3Possibilities.at(iter).at(2))[5] + CSymList->at(C3Possibilities.at(iter).at(3))[5] ) / 4.0 ); maxGrp = iter; } } }
1482 if ( C3Possibilities.at(maxGrp).size() == 4 )
1485 for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( C3Possibilities.at(maxGrp).size() ); it++ ) {
ProSHADE_internal_misc::addToDblPtrVector ( ret, CSymList->at(C3Possibilities.at(maxGrp).at(it)) ); }
1518 bool allAnglesMet =
true;
1519 proshade_double dotProduct;
1524 for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( grp->size() ); mIt++ )
1528 if ( ( ( 1.0 > ( dotProduct - axErr ) ) && ( 1.0 < ( dotProduct + axErr ) ) ) || ( ( -1.0 > ( dotProduct - axErr ) ) && ( -1.0 < ( dotProduct + axErr ) ) ) )
1530 if ( sym[5] > CSymList->at(grp->at(mIt))[5] )
1536 allAnglesMet =
false;
1537 return ( allAnglesMet );
1544 for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( grp->size() ); mIt++ )
1548 if ( ( angle > ( std::abs ( dotProduct ) - axErr ) ) &&
1549 ( angle < ( std::abs ( dotProduct ) + axErr ) ) )
1556 allAnglesMet =
false;
1562 return ( allAnglesMet );
1585 std::vector< proshade_double* > hlpVec;
1586 bool atLeastOne =
false;
1589 for ( proshade_unsign gIt = 0; gIt < static_cast<proshade_unsign> ( possibilities->size() ); gIt++ )
1591 if (
static_cast<proshade_unsign
> ( possibilities->at(gIt).size() ) == requiredNoAxes ) { atLeastOne =
true;
return ( atLeastOne ); }
1595 for ( proshade_unsign gIt = 0; gIt < static_cast<proshade_unsign> ( possibilities->size() ); gIt++ )
1598 if ( possibilities->at(gIt).size() < 2 ) {
continue; }
1607 if ( hlpVec.size() > 0 )
1613 for ( proshade_unsign axIt = 0; axIt < static_cast<proshade_unsign> ( hlpVec.size() ); axIt++ )
1628 if ( possibilities->at(gIt).size() == requiredNoAxes ) { atLeastOne =
true; }
1632 return ( atLeastOne );
1645 return ( a[0] < b[0] );
1666 proshade_double ret = 0.0;
1667 proshade_double curSum = 0.0;
1668 proshade_double maxVal = 0.0;
1669 proshade_double angStep = std::acos ( 1.0 - axErr ) / 2;
1670 std::vector< proshade_double* > angVec;
1679 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( std::floor ( ( 2.0 * M_PI / angStep ) / fold ) ); iter++ )
1685 for ( proshade_unsign angCmb = 0; angCmb < static_cast<proshade_unsign> ( fold ); angCmb++ )
1691 for ( proshade_unsign angIt = 0; angIt < static_cast<proshade_unsign> ( angVec.size() ); angIt++ )
1693 if ( angVec.at(angIt)[0] < ( ( iter*angStep ) + ( ( 2.0 * M_PI / fold ) * angCmb ) ) ) {
continue; }
1694 if ( angVec.at(angIt)[0] > ( ( (iter+1)*angStep ) + ( ( 2.0 * M_PI / fold ) * angCmb ) ) ) {
break; }
1696 if ( angVec.at(angIt)[1] > maxVal ) { maxVal = angVec.at(angIt)[1]; }
1700 curSum /=
static_cast<proshade_double
> ( fold );
1701 if ( ret < curSum ) { ret = curSum; }
1705 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( angVec.size() ); iter++ ) {
delete[] angVec.at(iter); }
1728 proshade_double euA, euB, euG, xPk, yPk, zPk, anglPk;
1729 proshade_double* rotMat =
new proshade_double [9];
1731 proshade_unsign arrIndex;
1732 std::vector< proshade_double* > angVec;
1735 for ( proshade_unsign xIt = 0; xIt < ( dataObj->
getMaxBand() * 2 ); xIt++ )
1737 for ( proshade_unsign yIt = 0; yIt < ( dataObj->
getMaxBand() * 2 ); yIt++ )
1739 for ( proshade_unsign zIt = 0; zIt < ( dataObj->
getMaxBand() * 2 ); zIt++ )
1742 arrIndex = zIt + ( dataObj->
getMaxBand() * 2 ) * ( yIt + ( dataObj->
getMaxBand() * 2 ) * xIt );
1746 static_cast<proshade_signed
> ( yIt ),
static_cast<proshade_signed
> ( zIt ),
1752 if ( ( ( std::max ( std::abs ( xPk ), std::max( std::abs ( yPk ), std::abs ( zPk ) ) ) == std::abs ( xPk ) ) && ( xPk < 0.0 ) ) ||
1753 ( ( std::max ( std::abs ( xPk ), std::max( std::abs ( yPk ), std::abs ( zPk ) ) ) == std::abs ( yPk ) ) && ( yPk < 0.0 ) ) ||
1754 ( ( std::max ( std::abs ( xPk ), std::max( std::abs ( yPk ), std::abs ( zPk ) ) ) == std::abs ( zPk ) ) && ( zPk < 0.0 ) ) )
1766 proshade_double* hlpArr =
new proshade_double [2];
1768 hlpArr[0] = anglPk + M_PI;
1802 proshade_double* hlpSym =
new proshade_double [6];
1806 hlpSym[0] =
static_cast<proshade_double
> ( fold );
1810 hlpSym[4] = ( 2.0 * M_PI ) /
static_cast<proshade_double
> ( fold );
1814 for ( proshade_unsign symIt = 0; symIt < static_cast<proshade_unsign> ( axVec->size() ); symIt++ )
1817 if ( axVec->at(symIt)[0] == hlpSym[0] )
1820 axVec->at(symIt)[2],
1821 axVec->at(symIt)[3],
1828 if ( axVec->at(symIt)[5] < hlpSym[5] )
1831 delete[] axVec->at(symIt);
1832 axVec->at(symIt) = hlpSym;
1870 if ( grp->size() < 2 ) {
return; }
1873 proshade_double axHeight = 0.0;
1874 proshade_double* symHlp =
new proshade_double[6];
1878 for ( proshade_unsign fAx = 0; fAx < static_cast<proshade_unsign> ( grp->size() ); fAx++ )
1880 for ( proshade_unsign sAx = 1; sAx < static_cast<proshade_unsign> ( grp->size() ); sAx++ )
1883 if ( fAx >= sAx ) {
continue; }
1887 CSymList->at(grp->at(fAx))[2],
1888 CSymList->at(grp->at(fAx))[3],
1889 CSymList->at(grp->at(sAx))[1],
1890 CSymList->at(grp->at(sAx))[2],
1891 CSymList->at(grp->at(sAx))[3], angle, angle );
1894 if ( ( ( std::max ( std::abs ( solVec.at(0) ), std::max( std::abs ( solVec.at(1) ), std::abs ( solVec.at(2) ) ) ) == std::abs ( solVec.at(0) ) ) && ( solVec.at(0) < 0.0 ) ) ||
1895 ( ( std::max ( std::abs ( solVec.at(0) ), std::max( std::abs ( solVec.at(1) ), std::abs ( solVec.at(2) ) ) ) == std::abs ( solVec.at(1) ) ) && ( solVec.at(1) < 0.0 ) ) ||
1896 ( ( std::max ( std::abs ( solVec.at(0) ), std::max( std::abs ( solVec.at(1) ), std::abs ( solVec.at(2) ) ) ) == std::abs ( solVec.at(2) ) ) && ( solVec.at(2) < 0.0 ) ) )
1898 solVec.at(0) *= -1.0;
1899 solVec.at(1) *= -1.0;
1900 solVec.at(2) *= -1.0;
1904 symHlp[1] = solVec.at(0); symHlp[2] = solVec.at(1); symHlp[3] = solVec.at(2);
1916 CSymList->at(grp->at(fAx))[2],
1917 CSymList->at(grp->at(fAx))[3],
1918 CSymList->at(grp->at(sAx))[1],
1919 CSymList->at(grp->at(sAx))[2],
1920 CSymList->at(grp->at(sAx))[3], -angle, -angle );
1923 if ( ( ( std::max ( std::abs ( solVec.at(0) ), std::max( std::abs ( solVec.at(1) ), std::abs ( solVec.at(2) ) ) ) == std::abs ( solVec.at(0) ) ) && ( solVec.at(0) < 0.0 ) ) ||
1924 ( ( std::max ( std::abs ( solVec.at(0) ), std::max( std::abs ( solVec.at(1) ), std::abs ( solVec.at(2) ) ) ) == std::abs ( solVec.at(1) ) ) && ( solVec.at(1) < 0.0 ) ) ||
1925 ( ( std::max ( std::abs ( solVec.at(0) ), std::max( std::abs ( solVec.at(1) ), std::abs ( solVec.at(2) ) ) ) == std::abs ( solVec.at(2) ) ) && ( solVec.at(2) < 0.0 ) ) )
1927 solVec.at(0) *= -1.0;
1928 solVec.at(1) *= -1.0;
1929 solVec.at(2) *= -1.0;
1933 symHlp[1] = solVec.at(0); symHlp[2] = solVec.at(1); symHlp[3] = solVec.at(2);
1969 std::vector< proshade_unsign > C3s, prospectiveC2s, C2PossibilitiesHlp;
1970 std::vector< std::vector< proshade_unsign > > C2Possibilities;
1971 proshade_double dotProd;
1979 for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( ret->size() ); rIt++ )
1982 for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( CSymList->size() ); cIt++ )
1985 if ( CSymList->at(cIt)[0] != 2.0 || CSymList->at(cIt)[5] < minPeakHeight ) {
continue; }
1989 &CSymList->at(cIt)[1], &CSymList->at(cIt)[2], &CSymList->at(cIt)[3] );
1996 C2Possibilities.clear(); C2PossibilitiesHlp.clear();
1997 for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( prospectiveC2s.size() ); cIt++ )
2000 groupMatched =
false;
2001 for ( proshade_unsign gIt = 0; gIt < static_cast<proshade_unsign> ( C2Possibilities.size() ); gIt++ )
2011 while ( C2Possibilities.size() != 0 )
2017 if ( C2Possibilities.at(0).size() == 3 )
2028 else { C2Possibilities.erase ( C2Possibilities.begin() ); }
2050 bool ProSHADE_internal_symmetry::testGroupAgainstGroup ( std::vector< proshade_double* >* GrList1, std::vector< proshade_unsign >* grp1, std::vector< proshade_double* >* GrList2, std::vector< proshade_unsign >* grp2, proshade_double angle, proshade_double axErr )
2054 proshade_double dotProduct;
2057 for ( proshade_unsign g1It = 0; g1It < static_cast<proshade_unsign> ( grp1->size() ); g1It++ )
2059 for ( proshade_unsign g2It = 0; g2It < static_cast<proshade_unsign> ( grp2->size() ); g2It++ )
2063 &GrList1->at(grp1->at(g1It))[2],
2064 &GrList1->at(grp1->at(g1It))[3],
2065 &GrList2->at(grp2->at(g2It))[1],
2066 &GrList2->at(grp2->at(g2It))[2],
2067 &GrList2->at(grp2->at(g2It))[3] );
2070 if ( ( angle > ( dotProduct - axErr ) ) && ( angle < ( dotProduct + axErr ) ) )
2096 std::vector< proshade_double* > ret;
2115 for ( proshade_unsign csIt = 0; csIt < static_cast<proshade_unsign> ( CSymList->size() ); csIt++ )
2117 for ( proshade_unsign retIt = 0; retIt < static_cast<proshade_unsign> ( ret.size() ); retIt++ )
2119 if ( ( CSymList->at(csIt)[0] == ret.at(retIt)[0] ) &&
2120 ( CSymList->at(csIt)[1] == ret.at(retIt)[1] ) &&
2121 ( CSymList->at(csIt)[2] == ret.at(retIt)[2] ) &&
2122 ( CSymList->at(csIt)[3] == ret.at(retIt)[3] ) &&
2123 ( CSymList->at(csIt)[4] == ret.at(retIt)[4] ) &&
2124 ( CSymList->at(csIt)[5] == ret.at(retIt)[5] ) )
2155 std::vector< proshade_unsign > C4List;
2156 proshade_double dotProduct;
2159 for ( proshade_unsign cSym = 0; cSym < static_cast<proshade_unsign> ( CSymList->size() ); cSym++ )
2165 for ( proshade_unsign c4 = 0; c4 < static_cast<proshade_unsign> ( C4List.size() ); c4++ )
2167 for ( proshade_unsign cSym = 0; cSym < static_cast<proshade_unsign> ( CSymList->size() ); cSym++ )
2170 if ( CSymList->at(cSym)[0] != 3 ) {
continue; }
2174 &CSymList->at(C4List.at(c4))[2],
2175 &CSymList->at(C4List.at(c4))[3],
2176 &CSymList->at(cSym)[1],
2177 &CSymList->at(cSym)[2],
2178 &CSymList->at(cSym)[3] );
2181 if ( ( ( 1.0 / sqrt ( 3.0 ) ) > ( dotProduct - axErr ) ) && ( ( 1.0 / sqrt ( 3.0 ) ) < ( dotProduct + axErr ) ) )
2209 std::vector< proshade_unsign > C4PossibilitiesHlp;
2210 std::vector< std::vector< proshade_unsign > > C4Possibilities;
2217 for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( CSymList->size() ); cIt++ )
2220 if ( CSymList->at(cIt)[0] != 4.0 || CSymList->at(cIt)[5] < minPeakHeight ) {
continue; }
2223 groupMatched =
false;
2224 for ( proshade_unsign gIt = 0; gIt < static_cast<proshade_unsign> ( C4Possibilities.size() ); gIt++ )
2237 proshade_double maxHeight = 0.0; proshade_unsign maxGrp = 0;
2238 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( C4Possibilities.size() ); iter++ ) {
if ( C4Possibilities.at(iter).size() == 3 ) {
if ( ( ( CSymList->at(C4Possibilities.at(iter).at(0))[5] + CSymList->at(C4Possibilities.at(iter).at(1))[5] + CSymList->at(C4Possibilities.at(iter).at(2))[5] ) / 3.0 ) > maxHeight ) { maxHeight = ( ( CSymList->at(C4Possibilities.at(iter).at(0))[5] + CSymList->at(C4Possibilities.at(iter).at(1))[5] + CSymList->at(C4Possibilities.at(iter).at(2))[5] ) / 3.0 ); maxGrp = iter; } } }
2240 if ( C4Possibilities.at(maxGrp).size() == 3 )
2243 for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( C4Possibilities.at(maxGrp).size() ); it++ ) {
ProSHADE_internal_misc::addToDblPtrVector ( ret, CSymList->at(C4Possibilities.at(maxGrp).at(it)) ); }
2275 std::vector< proshade_unsign > C4s, prospectiveC3s, C3PossibilitiesHlp;
2276 std::vector< std::vector< proshade_unsign > > C3Possibilities;
2277 proshade_double dotProd;
2285 for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( ret->size() ); rIt++ )
2288 for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( CSymList->size() ); cIt++ )
2291 if ( CSymList->at(cIt)[0] != 3.0 || CSymList->at(cIt)[5] < minPeakHeight ) {
continue; }
2301 C3Possibilities.clear(); C3PossibilitiesHlp.clear();
2302 for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( prospectiveC3s.size() ); cIt++ )
2305 groupMatched =
false;
2306 for ( proshade_unsign gIt = 0; gIt < static_cast<proshade_unsign> ( C3Possibilities.size() ); gIt++ )
2316 while ( C3Possibilities.size() != 0 )
2322 if ( C3Possibilities.at(0).size() == 4 )
2333 else { C3Possibilities.erase ( C3Possibilities.begin() ); }
2357 std::vector< proshade_unsign > prospectiveC2s, retGrp;
2358 proshade_double dotProd;
2359 proshade_unsign noPerpendicular, noSqrtTwo;
2365 for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( CSymList->size() ); cIt++ )
2368 if ( CSymList->at(cIt)[0] != 2.0 || CSymList->at(cIt)[5] < minPeakHeight ) {
continue; }
2371 noPerpendicular = 0; noSqrtTwo = 0;
2372 for ( proshade_unsign rIt = 0; rIt < 3; rIt++ )
2377 &CSymList->at(cIt)[1],
2378 &CSymList->at(cIt)[2],
2379 &CSymList->at(cIt)[3] );
2381 if ( ( std::abs ( dotProd ) > ( ( 1.0 / sqrt(2.0) ) - axErr ) ) && ( std::abs ( dotProd ) < ( ( 1.0 / sqrt(2.0) ) + axErr ) ) ) { noSqrtTwo += 1;
continue; }
2382 if ( ( std::abs ( dotProd ) > ( 0.0 - axErr ) ) && ( std::abs ( dotProd ) < ( 0.0 + axErr ) ) ) { noPerpendicular += 1;
continue; }
2386 if ( ( noSqrtTwo == 2 ) && ( noPerpendicular == 1 ) )
2394 if ( !
ProSHADE_internal_symmetry::findMissingAxesDual ( &prospectiveC2s, CSymList, ret, &retGrp, 6, axErr, 1, 0.0, 2, 1/sqrt(2.0), 2, dataObj ) )
2400 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( prospectiveC2s.size() ); iter++ )
2434 bool ProSHADE_internal_symmetry::findMissingAxesDual ( std::vector< proshade_unsign >* possibilities, std::vector< proshade_double* >* CSymList, std::vector< proshade_double* >* ret, std::vector< proshade_unsign >* retGroup, proshade_unsign requiredNoAxes, proshade_double axErr, proshade_unsign noMatchesG1, proshade_double angle1, proshade_unsign noMatchesG2, proshade_double angle2, proshade_unsign fold,
ProSHADE_internal_data::ProSHADE_data* dataObj )
2437 bool atLeastOne =
false;
2438 std::vector< proshade_double* > prosp;
2439 std::vector< proshade_double > sol;
2442 if (
static_cast<proshade_unsign
> ( possibilities->size() ) == requiredNoAxes ) { atLeastOne =
true;
return ( atLeastOne ); }
2445 for ( proshade_unsign prIt = 0; prIt < static_cast<proshade_unsign> ( possibilities->size() ); prIt++ )
2448 CSymList->at(possibilities->at(prIt))[1],
2449 CSymList->at(possibilities->at(prIt))[2],
2450 CSymList->at(possibilities->at(prIt))[3],
2451 CSymList->at(possibilities->at(prIt))[5], &prosp, axErr );
2455 for ( proshade_unsign rgIt1 = 0; rgIt1 < static_cast<proshade_unsign> ( retGroup->size() ); rgIt1++ )
2457 for ( proshade_unsign rgIt2 = 0; rgIt2 < static_cast<proshade_unsign> ( retGroup->size() ); rgIt2++ )
2460 if ( rgIt1 == rgIt2 ) {
continue; }
2464 ret->at(rgIt2)[1], ret->at(rgIt2)[2], ret->at(rgIt2)[3], angle1, angle2 );
2467 ProSHADE_internal_symmetry::checkFittingAxisDualAndSave ( retGroup, ret, fold, sol.at(0), sol.at(1), sol.at(2), &prosp, axErr, noMatchesG1, angle1, noMatchesG2, angle2, dataObj );
2468 if ( prosp.size() == requiredNoAxes ) {
break; }
2472 ret->at(rgIt2)[1], ret->at(rgIt2)[2], ret->at(rgIt2)[3], -angle1, -angle2 );
2475 ProSHADE_internal_symmetry::checkFittingAxisDualAndSave ( retGroup, ret, fold, sol.at(0), sol.at(1), sol.at(2), &prosp, axErr, noMatchesG1, angle1, noMatchesG2, angle2, dataObj );
2476 if ( prosp.size() == requiredNoAxes ) {
break; }
2479 if ( prosp.size() == requiredNoAxes ) {
break; }
2483 if (
static_cast<proshade_unsign
> ( prosp.size() ) == requiredNoAxes )
2486 for ( proshade_unsign iter =
static_cast<proshade_unsign
> ( possibilities->size() ); iter <
static_cast<proshade_unsign
> ( prosp.size() ); iter++ )
2498 return ( atLeastOne );
2503 for ( proshade_unsign iter =
static_cast<proshade_unsign
> ( possibilities->size() ); iter <
static_cast<proshade_unsign
> ( prosp.size() ); iter++ )
2505 delete[] prosp.at(iter);
2510 return ( atLeastOne );
2530 proshade_double* symHlp =
new proshade_double[6];
2534 symHlp[0] =
static_cast<proshade_double
> ( fold );
2538 symHlp[4] = 2.0 * M_PI / symHlp[0];
2539 symHlp[5] = axHeight;
2577 bool ProSHADE_internal_symmetry::checkFittingAxisDualAndSave ( std::vector< proshade_unsign >* retGroup, std::vector< proshade_double* >* ret, proshade_unsign fold, proshade_double axX, proshade_double axY, proshade_double axZ, std::vector< proshade_double* >* prosp, proshade_double axErr, proshade_unsign noMatchesG1, proshade_double angle1, proshade_unsign noMatchesG2, proshade_double angle2,
ProSHADE_internal_data::ProSHADE_data* dataObj )
2580 proshade_unsign noG1 = 0;
2581 proshade_unsign noG2 = 0;
2582 proshade_double dotProd = 0.0;
2583 proshade_double axHeight = 0.0;
2586 for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( retGroup->size() ); rIt++ )
2589 &ret->at(retGroup->at(rIt))[2],
2590 &ret->at(retGroup->at(rIt))[3],
2593 if ( ( std::abs ( dotProd ) > ( angle1 - axErr ) ) && ( std::abs ( dotProd ) < ( angle1 + axErr ) ) ) { noG1 += 1;
continue; }
2594 if ( ( std::abs ( dotProd ) > ( angle2 - axErr ) ) && ( std::abs ( dotProd ) < ( angle2 + axErr ) ) ) { noG2 += 1;
continue; }
2598 if ( ( noG1 == noMatchesG1 ) && ( noG2 == noMatchesG2 ) )
2604 if ( axHeight > 0.1 )
2606 proshade_unsign prevProsp =
static_cast<proshade_unsign
> ( prosp->size() );
2609 if (
static_cast<proshade_unsign
> ( prosp->size() ) > prevProsp ) {
return (
true ); }
2610 else {
return (
false ); }
2632 std::vector< proshade_double* > ret;
2651 for ( proshade_unsign csIt = 0; csIt < static_cast<proshade_unsign> ( CSymList->size() ); csIt++ )
2653 for ( proshade_unsign retIt = 0; retIt < static_cast<proshade_unsign> ( ret.size() ); retIt++ )
2655 if ( ( CSymList->at(csIt)[0] == ret.at(retIt)[0] ) &&
2656 ( CSymList->at(csIt)[1] == ret.at(retIt)[1] ) &&
2657 ( CSymList->at(csIt)[2] == ret.at(retIt)[2] ) &&
2658 ( CSymList->at(csIt)[3] == ret.at(retIt)[3] ) &&
2659 ( CSymList->at(csIt)[4] == ret.at(retIt)[4] ) &&
2660 ( CSymList->at(csIt)[5] == ret.at(retIt)[5] ) )
2694 std::vector< proshade_double* > ret;
2709 for ( proshade_unsign retIt = 0; retIt < static_cast < proshade_unsign > ( ret.size() ); retIt++ )
2742 std::vector< proshade_double* > ret;
2757 for ( proshade_unsign retIt = 0; retIt < static_cast < proshade_unsign > ( ret.size() ); retIt++ )
2786 std::vector< proshade_unsign > C5List;
2787 proshade_double dotProduct;
2790 for ( proshade_unsign cSym = 0; cSym < static_cast<proshade_unsign> ( CSymList->size() ); cSym++ )
2796 for ( proshade_unsign c5 = 0; c5 < static_cast<proshade_unsign> ( C5List.size() ); c5++ )
2798 for ( proshade_unsign cSym = 0; cSym < static_cast<proshade_unsign> ( CSymList->size() ); cSym++ )
2801 if ( CSymList->at(cSym)[0] != 3 ) {
continue; }
2805 &CSymList->at(C5List.at(c5))[2],
2806 &CSymList->at(C5List.at(c5))[3],
2807 &CSymList->at(cSym)[1],
2808 &CSymList->at(cSym)[2],
2809 &CSymList->at(cSym)[3] );
2812 if ( std::abs ( std::abs( std::sqrt ( ( 1.0 + 2.0 / std::sqrt ( 5.0 ) ) / 3.0 ) ) - std::abs( dotProduct ) ) < axErr )
2844 std::vector< proshade_unsign > C5PossibilitiesHlp;
2845 std::vector< std::vector< proshade_unsign > > C5Possibilities;
2852 for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( CSymList->size() ); cIt++ )
2855 if ( CSymList->at(cIt)[0] != 5.0 || CSymList->at(cIt)[5] < minPeakHeight ) {
continue; }
2858 groupMatched =
false;
2859 for ( proshade_unsign gIt = 0; gIt < static_cast<proshade_unsign> ( C5Possibilities.size() ); gIt++ )
2872 proshade_double maxHeight = 0.0; proshade_unsign maxGrp = 0;
2873 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( C5Possibilities.size() ); iter++ ) {
if ( C5Possibilities.at(iter).size() == 6 ) {
if ( ( ( CSymList->at(C5Possibilities.at(iter).at(0))[5] + CSymList->at(C5Possibilities.at(iter).at(1))[5] + CSymList->at(C5Possibilities.at(iter).at(2))[5] + CSymList->at(C5Possibilities.at(iter).at(3))[5] + CSymList->at(C5Possibilities.at(iter).at(4))[5] + CSymList->at(C5Possibilities.at(iter).at(5))[5] ) / 6.0 ) > maxHeight ) { maxHeight = ( ( CSymList->at(C5Possibilities.at(iter).at(0))[5] + CSymList->at(C5Possibilities.at(iter).at(1))[5] + CSymList->at(C5Possibilities.at(iter).at(2))[5] + CSymList->at(C5Possibilities.at(iter).at(3))[5] + CSymList->at(C5Possibilities.at(iter).at(4))[5] + CSymList->at(C5Possibilities.at(iter).at(5))[5] ) / 6.0 ); maxGrp = iter; } } }
2875 if ( C5Possibilities.at(maxGrp).size() == 6 )
2878 for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( C5Possibilities.at(maxGrp).size() ); it++ ) {
ProSHADE_internal_misc::addToDblPtrVector ( ret, CSymList->at(C5Possibilities.at(maxGrp).at(it)) ); }
2899 std::pair< proshade_unsign, proshade_unsign >
findBestIcosDihedralPair ( std::vector< proshade_double* >* CSymList, proshade_double minPeakHeight, proshade_double axErr )
2902 std::pair< proshade_unsign, proshade_unsign > ret;
2903 std::vector< proshade_unsign > C5List;
2904 proshade_double bestDihedralAngle = 999.9;
2905 proshade_double dotProduct;
2908 for ( proshade_unsign cSym = 0; cSym < static_cast<proshade_unsign> ( CSymList->size() ); cSym++ ) {
if ( CSymList->at(cSym)[0] == 5 && CSymList->at(cSym)[5] >= minPeakHeight ) {
ProSHADE_internal_misc::addToUnsignVector ( &C5List, cSym ); } }
2911 for ( proshade_unsign c5 = 0; c5 < static_cast<proshade_unsign> ( C5List.size() ); c5++ )
2913 for ( proshade_unsign cSym = 0; cSym < static_cast<proshade_unsign> ( CSymList->size() ); cSym++ )
2916 if ( CSymList->at(cSym)[0] != 3 ) {
continue; }
2917 if ( CSymList->at(cSym)[5] < minPeakHeight ) {
continue; }
2921 &CSymList->at(C5List.at(c5))[2],
2922 &CSymList->at(C5List.at(c5))[3],
2923 &CSymList->at(cSym)[1],
2924 &CSymList->at(cSym)[2],
2925 &CSymList->at(cSym)[3] );
2928 if ( std::abs ( std::abs( std::sqrt ( ( 1.0 + 2.0 / std::sqrt ( 5.0 ) ) / 3.0 ) ) - std::abs( dotProduct ) ) < axErr )
2930 if ( bestDihedralAngle > std::abs ( std::abs( std::sqrt ( ( 1.0 + 2.0 / std::sqrt ( 5.0 ) ) / 3.0 ) ) - std::abs( dotProduct ) ) )
2932 bestDihedralAngle = std::abs ( std::abs( std::sqrt ( ( 1.0 + 2.0 / std::sqrt ( 5.0 ) ) / 3.0 ) ) - std::abs( dotProduct ) );
2933 ret.first = C5List.at(c5);
2967 std::pair< proshade_unsign, proshade_unsign > initAxes =
findBestIcosDihedralPair ( CSymList, minPeakHeight, axErr );
2971 ProSHADE_internal_precomputedVals::icosahedronAxes.at(0).at(2),
2972 ProSHADE_internal_precomputedVals::icosahedronAxes.at(0).at(3),
2973 CSymList->at(initAxes.first)[1],
2974 CSymList->at(initAxes.first)[2],
2975 CSymList->at(initAxes.first)[3] );
2979 ProSHADE_internal_precomputedVals::icosahedronAxes.at(6).at(1),
2980 ProSHADE_internal_precomputedVals::icosahedronAxes.at(6).at(2),
2981 ProSHADE_internal_precomputedVals::icosahedronAxes.at(6).at(3) );
2984 proshade_double bestAng = 0.0, curAngDist, bestAngDist = 999.9;
2985 for ( proshade_double ang = 0.0; ang < ( M_PI * 2.0 ); ang += 0.01 )
2988 proshade_double* rotMatHlp =
new proshade_double[9];
2999 curAngDist = std::sqrt ( std::pow ( rotRotModelC3[0] - CSymList->at(initAxes.second)[1], 2.0 ) +
3000 std::pow ( rotRotModelC3[1] - CSymList->at(initAxes.second)[2], 2.0 ) +
3001 std::pow ( rotRotModelC3[2] - CSymList->at(initAxes.second)[3], 2.0 ) );
3004 if ( curAngDist < bestAngDist ) { bestAngDist = curAngDist; bestAng = ang; }
3008 delete[] rotRotModelC3;
3012 proshade_double* rotMat2 =
new proshade_double[9];
3020 for ( proshade_unsign iter = 0; iter < static_cast < proshade_unsign > ( ProSHADE_internal_precomputedVals::icosahedronAxes.size() ); iter++ )
3024 ProSHADE_internal_precomputedVals::icosahedronAxes.at(iter).at(1),
3025 ProSHADE_internal_precomputedVals::icosahedronAxes.at(iter).at(2),
3026 ProSHADE_internal_precomputedVals::icosahedronAxes.at(iter).at(3) );
3029 proshade_double* axis =
new proshade_double[6];
3032 axis[0] = ProSHADE_internal_precomputedVals::icosahedronAxes.at(iter).at(0);
3033 axis[1] = rotAxis[0];
3034 axis[2] = rotAxis[1];
3035 axis[3] = rotAxis[2];
3036 axis[4] = ( 2.0 * M_PI ) / axis[0];
3050 delete[] rotModelC3;
3064 std::pair< proshade_unsign, proshade_unsign >
findBestOctaDihedralPair ( std::vector< proshade_double* >* CSymList, proshade_double minPeakHeight, proshade_double axErr )
3067 std::pair< proshade_unsign, proshade_unsign > ret;
3068 std::vector< proshade_unsign > C4List;
3069 proshade_double bestDihedralAngle = 999.9;
3070 proshade_double dotProduct;
3073 for ( proshade_unsign cSym = 0; cSym < static_cast<proshade_unsign> ( CSymList->size() ); cSym++ ) {
if ( CSymList->at(cSym)[0] == 4 && CSymList->at(cSym)[5] >= minPeakHeight ) {
ProSHADE_internal_misc::addToUnsignVector ( &C4List, cSym ); } }
3076 for ( proshade_unsign c4 = 0; c4 < static_cast<proshade_unsign> ( C4List.size() ); c4++ )
3078 for ( proshade_unsign cSym = 0; cSym < static_cast<proshade_unsign> ( CSymList->size() ); cSym++ )
3081 if ( CSymList->at(cSym)[0] != 3 ) {
continue; }
3082 if ( CSymList->at(cSym)[5] < minPeakHeight ) {
continue; }
3086 &CSymList->at(C4List.at(c4))[2],
3087 &CSymList->at(C4List.at(c4))[3],
3088 &CSymList->at(cSym)[1],
3089 &CSymList->at(cSym)[2],
3090 &CSymList->at(cSym)[3] );
3093 if ( ( ( 1.0 / sqrt ( 3.0 ) ) > ( std::abs( dotProduct ) - axErr ) ) && ( ( 1.0 / sqrt ( 3.0 ) ) < ( std::abs( dotProduct ) + axErr ) ) )
3095 if ( bestDihedralAngle > std::abs( ( 1.0 / sqrt ( 3.0 ) ) - std::abs( dotProduct ) ) )
3097 bestDihedralAngle = std::abs( ( 1.0 / sqrt ( 3.0 ) ) - std::abs( dotProduct ) );
3098 ret.first = C4List.at(c4);
3133 std::pair< proshade_unsign, proshade_unsign > initAxes =
findBestOctaDihedralPair ( CSymList, minPeakHeight, axErr );
3137 ProSHADE_internal_precomputedVals::octahedronAxes.at(0).at(2),
3138 ProSHADE_internal_precomputedVals::octahedronAxes.at(0).at(3),
3139 CSymList->at(initAxes.first)[1],
3140 CSymList->at(initAxes.first)[2],
3141 CSymList->at(initAxes.first)[3] );
3145 ProSHADE_internal_precomputedVals::octahedronAxes.at(3).at(1),
3146 ProSHADE_internal_precomputedVals::octahedronAxes.at(3).at(2),
3147 ProSHADE_internal_precomputedVals::octahedronAxes.at(3).at(3) );
3150 proshade_double bestAng, curAngDist, bestAngDist = 999.9;
3151 for ( proshade_double ang = 0.0; ang < ( M_PI * 2.0 ); ang += 0.01 )
3154 proshade_double* rotMatHlp =
new proshade_double[9];
3165 curAngDist = std::sqrt ( std::pow ( rotRotModelC3[0] - CSymList->at(initAxes.second)[1], 2.0 ) +
3166 std::pow ( rotRotModelC3[1] - CSymList->at(initAxes.second)[2], 2.0 ) +
3167 std::pow ( rotRotModelC3[2] - CSymList->at(initAxes.second)[3], 2.0 ) );
3170 if ( curAngDist < bestAngDist ) { bestAngDist = curAngDist; bestAng = ang; }
3174 delete[] rotRotModelC3;
3178 proshade_double* rotMat2 =
new proshade_double[9];
3186 for ( proshade_unsign iter = 0; iter < static_cast < proshade_unsign > ( ProSHADE_internal_precomputedVals::octahedronAxes.size() ); iter++ )
3190 ProSHADE_internal_precomputedVals::octahedronAxes.at(iter).at(1),
3191 ProSHADE_internal_precomputedVals::octahedronAxes.at(iter).at(2),
3192 ProSHADE_internal_precomputedVals::octahedronAxes.at(iter).at(3) );
3195 proshade_double* axis =
new proshade_double[6];
3198 axis[0] = ProSHADE_internal_precomputedVals::octahedronAxes.at(iter).at(0);
3199 axis[1] = rotAxis[0];
3200 axis[2] = rotAxis[1];
3201 axis[3] = rotAxis[2];
3202 axis[4] = ( 2.0 * M_PI ) / axis[0];
3216 delete[] rotModelC3;
3239 std::vector< proshade_unsign > prospectiveC3s, retGrp;
3240 proshade_double dotProd;
3241 proshade_unsign noClose, noAway;
3247 for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( CSymList->size() ); cIt++ )
3250 if ( CSymList->at(cIt)[0] != 3.0 || CSymList->at(cIt)[0] < minPeakHeight ) {
continue; }
3253 noClose = 0; noAway = 0;
3254 for ( proshade_unsign rIt = 0; rIt < 6; rIt++ )
3259 &CSymList->at(cIt)[1],
3260 &CSymList->at(cIt)[2],
3261 &CSymList->at(cIt)[3] );
3263 if ( ( std::abs ( dotProd ) > ( ( std::sqrt ( ( 1.0 + 2.0 / std::sqrt ( 5.0 ) ) / 3.0 ) ) - axErr ) ) && ( std::abs ( dotProd ) < ( ( std::sqrt ( ( 1.0 + 2.0 / std::sqrt ( 5.0 ) ) / 3.0 ) ) + axErr ) ) ) { noClose += 1;
continue; }
3264 if ( ( std::abs ( dotProd ) > ( 1.0 - ( std::sqrt ( ( 1.0 + 2.0 / std::sqrt ( 5.0 ) ) / 3.0 ) ) - axErr ) ) && ( std::abs ( dotProd ) < ( 1.0 - ( std::sqrt ( ( 1.0 + 2.0 / std::sqrt ( 5.0 ) ) / 3.0 ) ) + axErr ) ) ) { noAway += 1;
continue; }
3268 if ( ( noClose == 3 ) && ( noAway == 3 ) )
3276 if ( !
ProSHADE_internal_symmetry::findMissingAxesDual ( &prospectiveC3s, CSymList, ret, &retGrp, 10, axErr, 3, std::sqrt ( ( 1.0 + 2.0 / std::sqrt ( 5.0 ) ) / 3.0 ), 3, 1.0 - ( std::sqrt ( ( 1.0 + 2.0 / std::sqrt ( 5.0 ) ) / 3.0 ) ), 3, dataObj ) )
3282 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( prospectiveC3s.size() ); iter++ )
3311 std::vector< proshade_unsign > prospectiveC2s, retGrp;
3312 proshade_double dotProd;
3313 proshade_unsign noClose, noMidway, noAway;
3319 for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( CSymList->size() ); cIt++ )
3322 if ( CSymList->at(cIt)[0] != 2.0 || CSymList->at(cIt)[0] < minPeakHeight ) {
continue; }
3325 noClose = 0; noMidway = 0; noAway = 0;
3326 for ( proshade_unsign rIt = 0; rIt < 6; rIt++ )
3331 &CSymList->at(cIt)[1],
3332 &CSymList->at(cIt)[2],
3333 &CSymList->at(cIt)[3] );
3335 if ( ( std::abs ( dotProd ) > ( ( sqrt ( 3.0 ) / 2.0 ) - axErr ) ) && ( std::abs ( dotProd ) < ( ( sqrt ( 3.0 ) / 2.0 ) + axErr ) ) ) { noAway += 1;
continue; }
3336 if ( ( std::abs ( dotProd ) > ( ( 1.0 / 2.0 ) - axErr ) ) && ( std::abs ( dotProd ) < ( ( 1.0 / 2.0 ) + axErr ) ) ) { noMidway += 1;
continue; }
3337 if ( ( std::abs ( dotProd ) > ( ( 0.0 ) - axErr ) ) && ( std::abs ( dotProd ) < ( ( 0.0 ) + axErr ) ) ) { noClose += 1;
continue; }
3341 if ( ( noClose == 2 ) && ( noMidway == 2 ) && ( noAway == 2 ) )
3349 if ( !
ProSHADE_internal_symmetry::findMissingAxesTriple ( &prospectiveC2s, CSymList, ret, &retGrp, 15, axErr, 2, 0.0, 2, 1.0/2.0, 2, sqrt ( 3.0 ) / 2.0, 2, dataObj ) )
3355 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( prospectiveC2s.size() ); iter++ )
3390 bool ProSHADE_internal_symmetry::findMissingAxesTriple ( std::vector< proshade_unsign >* possibilities, std::vector< proshade_double* >* CSymList, std::vector< proshade_double* >* ret, std::vector< proshade_unsign >* retGroup, proshade_unsign requiredNoAxes, proshade_double axErr, proshade_unsign noMatchesG1, proshade_double angle1, proshade_unsign noMatchesG2, proshade_double angle2, proshade_unsign noMatchesG3, proshade_double angle3, proshade_unsign fold,
ProSHADE_internal_data::ProSHADE_data* dataObj )
3393 bool atLeastOne =
false;
3394 std::vector< proshade_double* > prosp;
3395 std::vector< proshade_double > sol;
3398 if (
static_cast<proshade_unsign
> ( possibilities->size() ) == requiredNoAxes ) { atLeastOne =
true;
return ( atLeastOne ); }
3401 for ( proshade_unsign prIt = 0; prIt < static_cast<proshade_unsign> ( possibilities->size() ); prIt++ )
3404 CSymList->at(possibilities->at(prIt))[1],
3405 CSymList->at(possibilities->at(prIt))[2],
3406 CSymList->at(possibilities->at(prIt))[3],
3407 CSymList->at(possibilities->at(prIt))[5], &prosp, axErr );
3411 for ( proshade_unsign rgIt1 = 0; rgIt1 < static_cast<proshade_unsign> ( retGroup->size() ); rgIt1++ )
3413 for ( proshade_unsign rgIt2 = 0; rgIt2 < static_cast<proshade_unsign> ( retGroup->size() ); rgIt2++ )
3416 if ( rgIt1 == rgIt2 ) {
continue; }
3418 for ( proshade_unsign rgIt3 = 0; rgIt3 < static_cast<proshade_unsign> ( retGroup->size() ); rgIt3++ )
3421 if ( ( rgIt1 == rgIt3 ) || ( rgIt2 == rgIt3 ) ) {
continue; }
3425 ret->at(rgIt2)[1], ret->at(rgIt2)[2], ret->at(rgIt2)[3],
3426 ret->at(rgIt3)[1], ret->at(rgIt3)[2], ret->at(rgIt3)[3], angle1, angle2, angle3 );
3429 ProSHADE_internal_symmetry::checkFittingAxisTripleAndSave ( retGroup, ret, fold, sol.at(0), sol.at(1), sol.at(2), &prosp, axErr, noMatchesG1, angle1, noMatchesG2, angle2, noMatchesG3, angle3, dataObj );
3430 if ( prosp.size() == requiredNoAxes ) {
break; }
3434 ret->at(rgIt2)[1], ret->at(rgIt2)[2], ret->at(rgIt2)[3],
3435 ret->at(rgIt3)[1], ret->at(rgIt3)[2], ret->at(rgIt3)[3], -angle1, -angle2, -angle3 );
3438 ProSHADE_internal_symmetry::checkFittingAxisTripleAndSave ( retGroup, ret, fold, sol.at(0), sol.at(1), sol.at(2), &prosp, axErr, noMatchesG1, angle1, noMatchesG2, angle2, noMatchesG3, angle3, dataObj );
3439 if ( prosp.size() == requiredNoAxes ) {
break; }
3442 if ( prosp.size() == requiredNoAxes ) {
break; }
3445 if ( prosp.size() == requiredNoAxes ) {
break; }
3449 if ( prosp.size() == requiredNoAxes )
3452 for ( proshade_unsign axIt =
static_cast<proshade_unsign
> ( possibilities->size() ); axIt <
static_cast<proshade_unsign
> ( prosp.size() ); axIt++ )
3463 return ( atLeastOne );
3468 for ( proshade_unsign axIt =
static_cast<proshade_unsign
> ( possibilities->size() ); axIt <
static_cast<proshade_unsign
> ( prosp.size() ); axIt++ )
3470 delete[] prosp.at(axIt);
3475 return ( atLeastOne );
3501 void ProSHADE_internal_symmetry::checkFittingAxisTripleAndSave ( std::vector< proshade_unsign >* retGroup, std::vector< proshade_double* >* ret, proshade_unsign fold, proshade_double axX, proshade_double axY, proshade_double axZ, std::vector< proshade_double* >* prosp, proshade_double axErr, proshade_unsign noMatchesG1, proshade_double angle1, proshade_unsign noMatchesG2, proshade_double angle2, proshade_unsign noMatchesG3, proshade_double angle3,
ProSHADE_internal_data::ProSHADE_data* dataObj )
3504 proshade_unsign noG1 = 0;
3505 proshade_unsign noG2 = 0;
3506 proshade_unsign noG3 = 0;
3507 proshade_double dotProd = 0.0;
3508 proshade_double axHeight = 0.0;
3511 for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( retGroup->size() ); rIt++ )
3514 &ret->at(retGroup->at(rIt))[2],
3515 &ret->at(retGroup->at(rIt))[3],
3518 if ( ( std::abs ( dotProd ) > ( angle1 - axErr ) ) && ( std::abs ( dotProd ) < ( angle1 + axErr ) ) ) { noG1 += 1;
continue; }
3519 if ( ( std::abs ( dotProd ) > ( angle2 - axErr ) ) && ( std::abs ( dotProd ) < ( angle2 + axErr ) ) ) { noG2 += 1;
continue; }
3520 if ( ( std::abs ( dotProd ) > ( angle3 - axErr ) ) && ( std::abs ( dotProd ) < ( angle3 + axErr ) ) ) { noG3 += 1;
continue; }
3524 if ( ( noG1 == noMatchesG1 ) && ( noG2 == noMatchesG2 ) && ( noG3 == noMatchesG3 ) )
3530 if ( axHeight > 0.1 )
3557 std::vector< proshade_double* > ret, tmpHolder;
3558 std::vector< proshade_unsign > testedFolds;
3559 proshade_double symThres;
3560 proshade_unsign foldToTest;
3561 bool foldDone, anyNewSyms =
true;
3564 for ( proshade_unsign prIt = 0; prIt < static_cast< proshade_unsign > ( primes.size() ); prIt++ )
3567 std::stringstream hlpSS;
3568 hlpSS <<
"Searching for prime fold symmetry C" << primes.at(prIt) <<
".";
3572 std::vector< proshade_double* > prSyms = this->findRequestedCSymmetryFromAngleAxis ( settings, primes.at(prIt), &symThres );
3575 for ( proshade_unsign axIt = 0; axIt < static_cast< proshade_unsign > ( prSyms.size() ); axIt++ )
3578 if ( prSyms.at(axIt)[5] >= symThres )
3588 delete[] prSyms.at(axIt);
3593 if ( ret.size() < 1 ) {
return ( ret ); }
3596 while ( anyNewSyms )
3602 for ( proshade_unsign axIt1 = 0; axIt1 < static_cast< proshade_unsign > ( ret.size() ); axIt1++ )
3604 for ( proshade_unsign axIt2 = 0; axIt2 < static_cast< proshade_unsign > ( ret.size() ); axIt2++ )
3607 foldToTest = ret.at(axIt1)[0] * ret.at(axIt2)[0];
3612 for ( proshade_unsign fIt = 0; fIt < static_cast< proshade_unsign > ( testedFolds.size() ); fIt++ ) {
if ( testedFolds.at(fIt) == foldToTest ) { foldDone =
true;
break; } }
3613 if ( foldDone ) {
continue; }
3617 std::stringstream hlpSS2;
3618 hlpSS2 <<
"Searching for fold combination of detected folds " << ret.at(axIt1)[0] <<
" and " << ret.at(axIt2)[0] <<
", i.e. C" << foldToTest <<
".";
3622 std::vector< proshade_double* > prSyms = this->findRequestedCSymmetryFromAngleAxis ( settings, foldToTest, &symThres );
3625 for ( proshade_unsign newAxIt = 0; newAxIt < static_cast< proshade_unsign > ( prSyms.size() ); newAxIt++ )
3627 if ( prSyms.at(newAxIt)[5] >= symThres )
3637 delete[] prSyms.at(newAxIt);
3643 if ( tmpHolder.size() > 0 )
3645 for ( proshade_unsign tmpIt = 0; tmpIt < static_cast< proshade_unsign > ( tmpHolder.size() ); tmpIt++ )
3648 delete[] tmpHolder.at(tmpIt);
3652 tmpHolder.clear ( );
3672 return ( a[5] > b[5] );
3695 proshade_double soughtAngle;
3696 std::vector< proshade_double > allPeakHeights;
3697 std::vector< ProSHADE_internal_spheres::ProSHADE_rotFun_spherePeakGroup* > peakGroups;
3698 std::vector< proshade_double* > ret;
3702 this->sphereMappedRotFun.clear();
3705 for ( proshade_double angIt = 1.0; angIt < static_cast<proshade_double> ( fold ); angIt += 1.0 )
3708 soughtAngle = angIt * ( 2.0 * M_PI /
static_cast<proshade_double
> ( fold ) );
3712 M_PI /
static_cast<proshade_double
> ( fold ),
3713 this->maxShellBand * 2.0,
3715 static_cast<proshade_unsign
> ( angIt - 1.0 ) ) );
3718 this->sphereMappedRotFun.at(
static_cast<proshade_unsign
> ( angIt - 1.0 ))->interpolateSphereValues ( this->getInvSO3Coeffs ( ) );
3721 this->sphereMappedRotFun.at(
static_cast<proshade_unsign
> ( angIt - 1.0 ))->findAllPeaks ( settings->
peakNeighbours, &allPeakHeights );
3725 std::stringstream hlpSS;
3726 hlpSS <<
"Found a total of " << std::pow ( this->maxShellBand * 2.0 * (fold-1), 2.0 ) - allPeakHeights.size() <<
" non-peaks for thresholding.";
3733 std::stringstream hlpSS2;
3734 hlpSS2 <<
"Determined peak threshold " << *peakThres <<
".";
3738 for ( proshade_unsign shIt = 0; shIt < static_cast<proshade_unsign> ( this->sphereMappedRotFun.size() ); shIt++ )
3740 this->sphereMappedRotFun.at(shIt)->removeSmallPeaks ( *peakThres );
3744 for ( proshade_unsign sphIt = 0; sphIt < static_cast<proshade_unsign> ( this->sphereMappedRotFun.size() ); sphIt++ )
3747 for ( proshade_unsign pkIt = 0; pkIt < static_cast<proshade_unsign> ( this->sphereMappedRotFun.at(sphIt)->getPeaks().size() ); pkIt++ )
3751 for ( proshade_unsign pkGrpIt = 0; pkGrpIt < static_cast<proshade_unsign> ( peakGroups.size() ); pkGrpIt++ )
3753 if ( peakGroups.at(pkGrpIt)->checkIfPeakBelongs ( this->sphereMappedRotFun.at(sphIt)->getPeaks().at(pkIt).first, this->sphereMappedRotFun.at(sphIt)->getPeaks().at(pkIt).second, sphIt, settings->
axisErrTolerance, settings->
verbose ) ) { newPeak =
false;
break; }
3757 if ( !newPeak ) {
continue; }
3761 this->sphereMappedRotFun.at(sphIt)->getPeaks().at(pkIt).second,
3763 this->sphereMappedRotFun.at(sphIt)->getAngularDim() ) );
3768 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( peakGroups.size() ); grIt++ )
3771 std::stringstream hlpSS3;
3772 hlpSS3 <<
"Now considering group with LAT " << peakGroups.at(grIt)->getLatFromIndices() <<
" - " << peakGroups.at(grIt)->getLatToIndices() <<
" and LON " << peakGroups.at(grIt)->getLonFromIndices() <<
" - " << peakGroups.at(grIt)->getLonToIndices() <<
" spanning spheres ";
3773 for ( proshade_unsign sphIt = 0; sphIt < static_cast<proshade_unsign> ( peakGroups.at(grIt)->getSpherePositions().size() ); sphIt++ ) { hlpSS3 << peakGroups.at(grIt)->getSpherePositions().at(sphIt) <<
" ; "; }
3780 delete peakGroups.at(grIt);
3806 std::vector < proshade_unsign > folds;
3808 proshade_double lat, lon;
3809 proshade_double latSamlUnit = ( 2.0 * M_PI ) / ( dataObj->
maxShellBand * 2.0 );
3810 proshade_double lonSamlUnit = ( 1.0 * M_PI ) / ( dataObj->
maxShellBand * 2.0 );
3813 for ( proshade_unsign iter = 0; iter < static_cast < proshade_unsign > ( ret->size() ); iter++ )
3815 alreadyFound =
false;
3816 for ( proshade_unsign it = 0; it < static_cast < proshade_unsign > ( folds.size() ); it++ ) {
if ( folds.at(it) == ret->at(iter)[0] ) { alreadyFound =
true;
break; } }
3822 for ( proshade_unsign foldIt = 0; foldIt < static_cast < proshade_unsign > ( folds.size() ); foldIt++ )
3825 dataObj->sphereMappedRotFun.clear();
3828 for ( proshade_double angIt = 1.0; angIt < static_cast<proshade_double> ( folds.at(foldIt) ); angIt += 1.0 )
3832 M_PI /
static_cast < proshade_double
> ( folds.at(foldIt) ),
3834 ( 2.0 * M_PI ) /
static_cast < proshade_double
> ( folds.at(foldIt) ),
3835 static_cast<proshade_unsign
> ( angIt - 1.0 ) ) );
3838 dataObj->sphereMappedRotFun.at(
static_cast < proshade_unsign
> ( angIt - 1.0 ))->interpolateSphereValues ( dataObj->
getInvSO3Coeffs ( ) );
3842 for ( proshade_unsign axIt = 0; axIt < static_cast< proshade_unsign > ( ret->size() ); axIt++ )
3845 if ( ret->at(axIt)[0] != folds.at(foldIt) ) {
continue; }
3848 lat = std::atan2( ret->at(axIt)[2], ret->at(axIt)[1] ) / latSamlUnit;
3849 lon = std::acos ( ret->at(axIt)[3] ) / lonSamlUnit;
3851 if ( lat < 0.0 ) { lat += ( dataObj->
maxShellBand * 2.0 ); }
3852 if ( lon < 0.0 ) { lon += ( dataObj->
maxShellBand * 2.0 ); }
3854 lat = std::round ( lat );
3855 lon = std::round ( lon );
3861 for ( proshade_unsign sphIt = 0; sphIt < static_cast<proshade_unsign> ( dataObj->sphereMappedRotFun.size() ); sphIt++ )
3876 std::vector < proshade_double* > detectedAxis;
3880 ret->at(axIt)[5] = detectedAxis.at(0)[5];
3883 for ( proshade_unsign i = 0; i < static_cast < proshade_unsign > ( detectedAxis.size() ); i++ ) {
delete detectedAxis.at(i); }
3900 std::pair< proshade_unsign, proshade_unsign >
findBestTetraDihedralPair ( std::vector< proshade_double* >* CSymList, proshade_double minPeakHeight, proshade_double axErr )
3903 std::pair< proshade_unsign, proshade_unsign > ret;
3904 std::vector< proshade_unsign > C3List;
3905 proshade_double bestDihedralAngle = 999.9;
3906 proshade_double dotProduct;
3909 for ( proshade_unsign cSym = 0; cSym < static_cast<proshade_unsign> ( CSymList->size() ); cSym++ ) {
if ( CSymList->at(cSym)[0] == 3 && CSymList->at(cSym)[5] >= minPeakHeight ) {
ProSHADE_internal_misc::addToUnsignVector ( &C3List, cSym ); } }
3912 for ( proshade_unsign c3 = 0; c3 < static_cast<proshade_unsign> ( C3List.size() ); c3++ )
3914 for ( proshade_unsign cSym = 0; cSym < static_cast<proshade_unsign> ( CSymList->size() ); cSym++ )
3917 if ( CSymList->at(cSym)[0] != 2 ) {
continue; }
3918 if ( CSymList->at(cSym)[5] < minPeakHeight ) {
continue; }
3922 &CSymList->at(C3List.at(c3))[2],
3923 &CSymList->at(C3List.at(c3))[3],
3924 &CSymList->at(cSym)[1],
3925 &CSymList->at(cSym)[2],
3926 &CSymList->at(cSym)[3] );
3929 if ( ( ( 1.0 / sqrt ( 3.0 ) ) > ( std::abs( dotProduct ) - axErr ) ) && ( ( 1.0 / sqrt ( 3.0 ) ) < ( std::abs( dotProduct ) + axErr ) ) )
3931 if ( bestDihedralAngle > std::abs( ( 1.0 / sqrt ( 3.0 ) ) - std::abs( dotProduct ) ) )
3933 bestDihedralAngle = std::abs( ( 1.0 / sqrt ( 3.0 ) ) - std::abs( dotProduct ) );
3934 ret.first = C3List.at(c3);
3963 std::vector< proshade_double* > ret;
3978 for ( proshade_unsign retIt = 0; retIt < static_cast < proshade_unsign > ( ret.size() ); retIt++ )
4020 ProSHADE_internal_precomputedVals::tetrahedronAxes.at(0).at(2),
4021 ProSHADE_internal_precomputedVals::tetrahedronAxes.at(0).at(3),
4022 CSymList->at(initAxes.first)[1],
4023 CSymList->at(initAxes.first)[2],
4024 CSymList->at(initAxes.first)[3] );
4028 ProSHADE_internal_precomputedVals::tetrahedronAxes.at(4).at(1),
4029 ProSHADE_internal_precomputedVals::tetrahedronAxes.at(4).at(2),
4030 ProSHADE_internal_precomputedVals::tetrahedronAxes.at(4).at(3) );
4033 proshade_double bestAng, curAngDist, bestAngDist = 999.9;
4034 for ( proshade_double ang = 0.0; ang < ( M_PI * 2.0 ); ang += 0.01 )
4037 proshade_double* rotMatHlp =
new proshade_double[9];
4048 curAngDist = std::sqrt ( std::pow ( rotRotModelC2[0] - CSymList->at(initAxes.second)[1], 2.0 ) +
4049 std::pow ( rotRotModelC2[1] - CSymList->at(initAxes.second)[2], 2.0 ) +
4050 std::pow ( rotRotModelC2[2] - CSymList->at(initAxes.second)[3], 2.0 ) );
4053 if ( curAngDist < bestAngDist ) { bestAngDist = curAngDist; bestAng = ang; }
4057 delete[] rotRotModelC2;
4061 proshade_double* rotMat2 =
new proshade_double[9];
4069 for ( proshade_unsign iter = 0; iter < static_cast < proshade_unsign > ( ProSHADE_internal_precomputedVals::tetrahedronAxes.size() ); iter++ )
4073 ProSHADE_internal_precomputedVals::tetrahedronAxes.at(iter).at(1),
4074 ProSHADE_internal_precomputedVals::tetrahedronAxes.at(iter).at(2),
4075 ProSHADE_internal_precomputedVals::tetrahedronAxes.at(iter).at(3) );
4078 proshade_double* axis =
new proshade_double[6];
4081 axis[0] = ProSHADE_internal_precomputedVals::tetrahedronAxes.at(iter).at(0);
4082 axis[1] = rotAxis[0];
4083 axis[2] = rotAxis[1];
4084 axis[3] = rotAxis[2];
4085 axis[4] = ( 2.0 * M_PI ) / axis[0];
4099 delete[] rotModelC2;