Skip to content

Commit

Permalink
Refs #5772 Tweaked cutoffs, etc., to better match ISAW results
Browse files Browse the repository at this point in the history
  • Loading branch information
@RuthFromDuluth
RuthFromDuluth committed Oct 14, 2012
1 parent 9192cec commit cd73759
Show file tree
Hide file tree
Showing 4 changed files with 92 additions and 64 deletions.
90 changes: 51 additions & 39 deletions Code/Mantid/Framework/Crystal/src/IntegratePeakTimeSlices.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -251,6 +251,7 @@ namespace Mantid
double &Varxx, double &Varxy, double &Varyy,
std::vector<double>&ParameterValues)
{
UNUSED_ARG(ParameterValues);
double Den = StatBase[IIntensities]-background*StatBase[ISS1];
Varxx =(StatBase[ISSIxx]-2*meanx*StatBase[ISSIx]+meanx*meanx*StatBase[IIntensities]-
background*(StatBase[ISSxx]-2*meanx*StatBase[ISSx]+meanx*meanx*StatBase[ISS1]))
Expand All @@ -264,12 +265,14 @@ namespace Mantid
meanx*meany*StatBase[IIntensities]-
background*(StatBase[ISSxy]-meanx*StatBase[ISSy]-meany*StatBase[ISSx]+meanx*meany*StatBase[ISS1]))
/Den;
if( ParameterValues.size() >IVXX)

if( CalcVariances())
{
Varxx = std::min<double>(Varxx, 1.21*ParameterValues[IVXX]);//copied from BiVariateNormal
Varxx = std::max<double>(Varxx, .79*ParameterValues[IVXX]);
Varyy = std::min<double>(Varyy, 1.21*ParameterValues[IVYY]);
Varyy = std::max<double>(Varyy, .79*ParameterValues[IVYY]);

Varxx = std::min<double>(Varxx, 1.21*getInitVarx());//copied from BiVariateNormal
Varxx = std::max<double>(Varxx, .79*getInitVarx());
Varyy = std::min<double>(Varyy, 1.21*getInitVary());
Varyy = std::max<double>(Varyy, .79*getInitVary());
}

};
Expand All @@ -288,25 +291,28 @@ namespace Mantid
double backVar= TotBoundaryVariances/nBoundaryCells/nBoundaryCells;
double IntensVar = Variance +ncells*ncells*backVar;

double relError = .15;
double relError = .25;

if( back_calc != back )
relError = .35;
relError = .45;

int N = NParameters;
if( CalcVariances)
N = N-3;


double NSigs=2;
if( back_calc >0)
NSigs = std::max<double>( 2,7-5*back_calc/back);//background too high
ostringstream str;
str<< max<double>(0.0, back_calc-(1+2*relError)*sqrt(backVar))<<"<Background<"<<(back+(1+2*relError)*sqrt(backVar))
<<","<< max<double>(0.0, Intensity_calc-(1+2*relError)*sqrt(IntensVar))<<
"<Intensity<"<<TotIntensity-ncells*back+(1+2*relError)*sqrt(IntensVar);
str<< max<double>(0.0, back_calc-NSigs*(1+relError)*sqrt(backVar))<<"<Background<"<<(back+NSigs*(1+relError)*sqrt(backVar))
<<","<< max<double>(0.0, Intensity_calc-NSigs*(1+relError)*sqrt(IntensVar))<<
"<Intensity<"<<TotIntensity-ncells*back+NSigs*(1+relError)*sqrt(IntensVar);

double min =max<double>(0.0, back_calc-(1+2*relError)*sqrt(backVar));
double maxx =back+(1+2*relError)*sqrt(backVar);
double min =max<double>(0.0, back_calc-NSigs*(1+relError)*sqrt(backVar));
double maxx =back+NSigs*(1+relError)*sqrt(backVar);
Bounds.push_back( pair<double,double>(min,maxx ));
Bounds.push_back( pair<double,double>(max<double>(0.0, Intensity_calc-(1+2*relError)*sqrt(IntensVar)),TotIntensity-ncells*back+(1+2*relError)*sqrt(IntensVar)));
Bounds.push_back( pair<double,double>(max<double>(0.0, Intensity_calc-NSigs*(1+relError)*sqrt(IntensVar)),
TotIntensity-ncells*back+NSigs*(1+relError)*sqrt(IntensVar)));
/* for( int i=2; i<N; i++ )
{
if( i>0)
Expand All @@ -324,7 +330,7 @@ namespace Mantid
*/
double val =col_calc;

double relErr1=.75*relError;
double relErr1= relError*.75;


str<<","<<(1-relErr1)*val <<"<"<<"Mcol"<<"<"<<(1+relErr1)*val;
Expand All @@ -344,9 +350,9 @@ namespace Mantid
Bounds.push_back(pair<double,double>((1-relErr1)*val,(1+relErr1)*val));


val = Vxy_calc;
str<<","<<(1-relErr1)*val <<"<"<<"SSrc"<<"<"<<(1+relErr1)*val;
Bounds.push_back(pair<double,double>((1-relErr1)*val,(1+relErr1)*val));
// val = Vxy_calc;
// str<<","<<(1-relErr1)*val <<"<"<<"SSrc"<<"<"<<(1+relErr1)*val;
// Bounds.push_back(pair<double,double>((1-relErr1)*val,(1+relErr1)*val));

}

Expand Down Expand Up @@ -496,6 +502,8 @@ namespace Mantid
fit_alg->executeAsSubAlg();

chisqOverDOF = fit_alg->getProperty("OutputChi2overDoF");
std::string outputStatus = fit_alg->getProperty( "OutputStatus" );
g_log.debug()<<"Chisq/OutputStatus="<< chisqOverDOF<< "/"<< outputStatus<< std::endl;

names.clear(); params.clear();errs.clear();
ITableWorkspace_sptr RRes = fit_alg->getProperty( "OutputParameters");
Expand Down Expand Up @@ -678,7 +686,7 @@ namespace Mantid

R = min<double> (36*max<double>(CellWidth,CellHeight), R);
R = max<double> (6*max<double>(CellWidth,CellHeight), R);
R=1.1*R;//Gets a few more background cells.
R=1.4*R;//Gets a few more background cells.
int Chan;

Mantid::MantidVec X = inpWkSpace->dataX(wsIndx);
Expand Down Expand Up @@ -721,12 +729,13 @@ namespace Mantid
double MaxCounts = -1;

for( int dir =1 ; dir >-2; dir -=2)
for( int t= 0; t < dChan && !done; t++)
{
done = false;
for( int t= 0; t < dChan && !done; t++)
if( dir < 0 && t==0 )
{
Centy = Row0;
Centx = Col0;
done = false;
}
else if( Chan+dir*t <0 || Chan+dir*t >= (int)X.size())
done = true;
Expand Down Expand Up @@ -767,7 +776,7 @@ namespace Mantid
done = true;

}

}
if( MaxChan >0)
Chan = MaxChan ;

Expand Down Expand Up @@ -976,8 +985,8 @@ namespace Mantid

if( dir ==1 && (chan ==0||chan==1))
{
lastAttributeList->case4 = true;
origAttributeList =lastAttributeList;
AttributeValues->case4 = true;
origAttributeList =AttributeValues;
}else
LastTableRow = -1;

Expand Down Expand Up @@ -1304,16 +1313,16 @@ namespace Mantid
double Vx,Vy;
Vx = Vy =HalfWidthAtHalfHeightRadius*HalfWidthAtHalfHeightRadius;

if( EdgeX < 0)
if( EdgeX < .5*Vx )
Vx = std::max<double>(HalfWidthAtHalfHeightRadius*HalfWidthAtHalfHeightRadius,
Vx_calc);
Vx_calc);

if( EdgeY < 0)
if( EdgeY < .5*Vy)
Vy =std::max<double>(HalfWidthAtHalfHeightRadius*HalfWidthAtHalfHeightRadius,
Vy_calc);
Vy_calc);

double DD = max<double>( sqrt(Vy)*CellHeight, sqrt(Vx)*CellWidth);
double NewRadius= 1.4* max<double>( 5*max<double>(CellWidth,CellHeight), 4*DD);
double NewRadius= 1.4* max<double>( 5*max<double>(CellWidth,CellHeight), 4.5*DD);
NewRadius = min<double>(baseRCRadius, NewRadius);
//1.4 is needed to get more background cells. In rectangle the corners were background

Expand Down Expand Up @@ -1615,6 +1624,8 @@ namespace Mantid

AttributeValues->updateEdgeYsize(abs(ROW-row));
AttributeValues->updateEdgeXsize(abs(COL-col));
AttributeValues->updateEdgeYsize(abs(row));
AttributeValues->updateEdgeXsize(abs(col));
}
}
}
Expand Down Expand Up @@ -1721,7 +1732,7 @@ namespace Mantid

double x = params[IIntensity] / (AttributeValues->StatBaseVals(IIntensities) - params[Ibk] * ncells);

if ((x < .7 || x > 1.45)&& !EdgePeak)// The fitted intensity should be close to tot intensity - background
if ((x < .25 || x > 2.5)&& !EdgePeak)// The fitted intensity should be close to tot intensity - background
{
g_log.debug()<< " Bad Slice. Fitted Intensity & Observed Intensity(-back) too different. ratio="
<<x<<std::endl;
Expand Down Expand Up @@ -1791,8 +1802,8 @@ namespace Mantid
if (!GoodNums)

{ g_log.debug()<<
" Bad Slice.Background or Intensity params or errs are out of bounds or relative errors are too high";

" Bad Slice.Background or Intensity params or errs are out of bounds or relative errors are too high"
<<std::endl;

return false;
}
Expand All @@ -1813,8 +1824,8 @@ namespace Mantid
return false;
}

double Nrows = AttributeValues->getCurrentRadius();
if( maxPeakHeightTheoretical < 1 && (params[IVXX] > Nrows*Nrows/20 || params[IVYY] > Nrows*Nrows/20))
double Nrows = std::max<double>(AttributeValues->StatBase[INRows], AttributeValues->StatBase[INCol]);
if( maxPeakHeightTheoretical < 1 && (params[IVXX] > Nrows*Nrows/4 || params[IVYY] > Nrows*Nrows/4))
{
g_log.debug()<<"Peak is too flat "<<std::endl;
return false;
Expand Down Expand Up @@ -1945,27 +1956,28 @@ namespace Mantid
UNUSED_ARG( ncells);
// int Ibk = find("Background", names);
double err =0;
double intensity=0;

if( !EdgePeak )
{
err = AttributeValues->CalcISAWIntensityVariance(params.data(),errs.data(), chisqdivDOF);
//CalculateIsawIntegrateError(params[Ibk], errs[Ibk], chisqdivDOF, TotSliceVariance, ncells);

TotIntensity += AttributeValues->CalcISAWIntensity( params.data());
intensity =AttributeValues->CalcISAWIntensity( params.data());
TotIntensity += intensity;
//TotSliceIntensity - params[IBACK] * ncells;

TotVariance += err ;

}else
{
int IInt = find("Intensity", names);

TotIntensity += params[IInt];
intensity = params[IInt];
TotIntensity +=intensity;
err = errs[IInt];
TotVariance += err * err;

}

g_log.debug()<<"Slice/Tot intensity;Var="<<intensity<<";"<<err*err<<"/"<<TotIntensity<<";"<<TotVariance<<std::endl;
}

bool DataModeHandler::CalcVariances( )
Expand Down
10 changes: 5 additions & 5 deletions Code/Mantid/Framework/Crystal/test/IntegratePeakTimeSlicesTest.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -196,7 +196,7 @@ class IntegratePeakTimeSlicesTest: public CxxTest::TestSuite

TS_ASSERT_LESS_THAN(fabs(intensity -60000), 1500.0);
//Not sure why this reduced the error so much in the test
TS_ASSERT_LESS_THAN(fabs(sigma -539), 21.0);
TS_ASSERT_LESS_THAN(fabs(sigma -502), 21.0);


TS_ASSERT_EQUALS( Twk->rowCount(), 7);
Expand All @@ -208,16 +208,16 @@ class IntegratePeakTimeSlicesTest: public CxxTest::TestSuite

TS_ASSERT_LESS_THAN(fabs(Twk->getRef<double> (std::string("Background"), 1) - 1.2619 ), .5);

TS_ASSERT_LESS_THAN(fabs(Twk->getRef<double> ("Intensity", 2) - 11309.8 ), 120);
TS_ASSERT_LESS_THAN(fabs(Twk->getRef<double> ("Intensity", 2) - 11514 ), 120);


TS_ASSERT_LESS_THAN(fabs(Twk->getRef<double> ("NCells", 3) - 553), 5);
TS_ASSERT_LESS_THAN(fabs(Twk->getRef<double> ("NCells", 3) - 885), 5);


TS_ASSERT_LESS_THAN(fabs(Twk->getRef<double> ("ChiSqrOverDOF", 4) - 60.4183), 3.5);
TS_ASSERT_LESS_THAN(fabs(Twk->getRef<double> ("ChiSqrOverDOF", 4) - 35.3), 3.5);


TS_ASSERT_LESS_THAN(fabs(Twk->getRef<double> ("TotIntensity", 0) - 5298.4 ), 10);
TS_ASSERT_LESS_THAN(fabs(Twk->getRef<double> ("TotIntensity", 0) - 6228 ), 10);



Expand Down
4 changes: 2 additions & 2 deletions Code/Mantid/Framework/Crystal/test/PeakIntegrationTest.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -204,10 +204,10 @@ class PeakIntegrationTest : public CxxTest::TestSuite
double sigIntensity =peak.getSigmaIntensity();
//std::cout<<"Peak Intens,sig,slice="<<intensity;
//std::cout<<","<<sigIntensity<<","<<IC<<std::endl;
double intensity0=2509;
double intensity0=4470;
TS_ASSERT_DELTA(intensity,intensity0, 400.0);

double sig0=100;
double sig0=180;


TS_ASSERT_DELTA( sigIntensity,sig0, 25.0);
Expand Down
52 changes: 34 additions & 18 deletions Code/Mantid/Framework/CurveFitting/src/BivariateNormal.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,6 +85,10 @@ BivariateNormal::~BivariateNormal()

initCoeff(D, X, Y, coefNorm, expCoeffx2, expCoeffy2, expCoeffxy, NCells, Varxx,Varxy,Varyy);

bool badParams=false;
if( Varxx<=0 || Varyy <=0 || Varxy*Varxy >.9*Varxx*Varyy)
badParams = true;

NCells = std::min<int>((int)nData, NCells);

double Background = getParameter(IBACK);
Expand All @@ -110,7 +114,7 @@ BivariateNormal::~BivariateNormal()
double pen =0;
if( (i-1)%10==0)
pen = penalty;
if (isNaNs)
if (isNaNs )
out[x] = 10000;
else
{
Expand All @@ -119,11 +123,12 @@ BivariateNormal::~BivariateNormal()
out[x] = Background + coefNorm * Intensity * exp(expCoeffx2 * dx * dx + expCoeffxy * dx * dy
+ expCoeffy2 * dy * dy) + pen;

if (out[x] != out[x])
if (out[x] != out[x] )
{
out[x] = 100000;
isNaNs = true;
}
}else if( badParams)
out[x] = 10000;

}
double diff = out[x]-D[x];
Expand All @@ -133,7 +138,14 @@ BivariateNormal::~BivariateNormal()

x++;
}

inf<<"Constr:";
for (size_t i=0;i< nParams(); i++)
{
IConstraint* constr = getConstraint( (size_t)i);
if( constr)
inf<<i<<"="<<constr->check()<<";";
}
inf << std::endl;
inf << std::endl<<" chiSq =" << chiSq <<" nData "<<nData<<std::endl;
g_log.debug(inf.str());

Expand All @@ -156,7 +168,8 @@ void BivariateNormal::functionDeriv1D(API::Jacobian *out, const double *xValues,
inf << std::endl;
g_log.debug( inf.str());


// std::vector<double>outf(nData,0.0);
// function1D( outf.data(), xValues, nData);
double penDeriv =0;
double uu = LastParams[IVXX] * LastParams[IVYY] - LastParams[IVXY] * LastParams[IVXY];

Expand Down Expand Up @@ -579,15 +592,16 @@ void BivariateNormal::initCommon()
int NCells1;
double Varxx, Varxy,Varyy;


Varxx = Varxy = Varyy=-1;
initCoeff( D, X, Y, coefNorm, expCoeffx2, expCoeffy2, expCoeffxy,
NCells1, Varxx, Varxy,Varyy);

// if( CalcVariances && nParams() < 5 && Varx0 < 0 && Vary0 < 0)
{
Varx0 = Varxx;
Vary0 = Varyy;
}
if( Varx0 < 0)
{
Varx0 = Varxx;
Vary0 = Varyy;
}

LastParams[IVXX] = Varxx;
LastParams[IVXY] = Varxy;
LastParams[IVYY] = Varyy;
Expand Down Expand Up @@ -631,7 +645,7 @@ void BivariateNormal::initCoeff( const MantidVec &D,
TotN) ;

if( Varx0 > 0) Varxx = std::min<double>(Varxx, 1.21*Varx0);
if( Varx0 > 0) Varxx = std::max<double>(Varxx, .79*Varx0);
if( Varx0 > 0 ) Varxx = std::max<double>(Varxx, .79*Varx0);

} else {
Varxx = getParameter( IVXX);
Expand All @@ -642,13 +656,15 @@ void BivariateNormal::initCoeff( const MantidVec &D,

if( CalcVyy||nParams() <6)
{
Varyy = (SIyy +(getParameter("Mrow")- (mIy))*(getParameter("Mrow")- (mIy) )*
TotI -getParameter("Background")* (Syy) -getParameter("Background")*(getParameter("Mrow")- (my) )*
(getParameter("Mrow")-
(my) )* TotN )/( TotI -getParameter("Background")*
double Mrow = getParameter("Mrow");
double Background = getParameter("Background");
Varyy = (SIyy +(Mrow- (mIy))*(Mrow- (mIy) )*
TotI -getParameter("Background")* (Syy) -Background*(Mrow- (my) )*
(Mrow-
(my) )* TotN )/( TotI -Background*
TotN );
if( Vary0 > 0) Varyy = std::min<double>(Varyy, 1.21*Vary0);
if( Vary0 > 0)Varyy = std::max<double>(Varyy, .79*Vary0);
if( Vary0 > 0 ) Varyy = std::min<double>(Varyy, 1.21*Vary0);
if( Vary0 > 0 )Varyy = std::max<double>(Varyy, .79*Vary0);

} else {
Varyy = getParameter( IVYY);
Expand Down

0 comments on commit cd73759

Please sign in to comment.