Skip to content
Snippets Groups Projects
Commit e63ee18d authored by Paul McCarthy's avatar Paul McCarthy :mountain_bicyclist:
Browse files

STY: whitespace

parent 5f1495dc
No related branches found
No related tags found
1 merge request!8Mnt/conda
Hide whitespace changes
Inline Side-by-side
histogram.cc
+ 42
− 76
View file @ e63ee18d
......@@ -14,13 +14,13 @@ using namespace NEWMAT;
namespace MISCMATHS {
float Histogram::getPercentile(float perc)
float Histogram::getPercentile(float perc)
{
if(histogram.Nrows()==0) generate();
generateCDF();
float percentile=getValue(1);
for (int i=2;i<=CDF.Nrows();i++){
if(CDF(i)>perc){
double diff=(CDF(i)-perc)/(CDF(i)-CDF(i-1));
......@@ -35,9 +35,9 @@ namespace MISCMATHS {
void Histogram::generate()
{
Tracer ts("Histogram::generate");
int size = sourceData.Nrows();
if(calcRange)
{
// calculate range automatically
......@@ -50,21 +50,21 @@ namespace MISCMATHS {
histMin=sourceData(i);
}
}
// zero histogram
histogram.ReSize(bins);
histogram=0;
// create histogram; the MIN is so that the maximum value falls in the
// last valid bin, not the (last+1) bin
for(int i=1; i<=size; i++)
{
{
histogram(getBin(sourceData(i)))++;
}
datapoints=size;
}
void Histogram::generate(ColumnVector exclude)
{
Tracer ts("Histogram::generate");
......@@ -81,13 +81,13 @@ namespace MISCMATHS {
histMin=sourceData(i);
}
}
histogram.ReSize(bins);
histogram=0;
datapoints=size;
for(int i=1; i<=size; i++)
{
if(exclude(i)>0){
......@@ -95,7 +95,7 @@ namespace MISCMATHS {
}else datapoints--;
}
}
......@@ -108,19 +108,19 @@ namespace MISCMATHS {
// smooth in i direction
newhist=0;
ColumnVector kernel(3);
ColumnVector kernel(3);
// corresponds to Gaussian with sigma=0.8 voxels
// kernel(1)=0.5;
// kernel(2)=0.2283;
// kernel(2)=0.2283;
// kernel(3)=0.0219;
// corresponds to Gaussian with sigma=0.6 voxels
// kernel(1)=0.6638;
// kernel(2)=0.1655;
// kernel(2)=0.1655;
// kernel(3)=0.0026;
//gauss(0.5,5,1)
kernel(1)=0.7866;
kernel(2)=0.1065;
kernel(2)=0.1065;
kernel(3)=0.0003;
for(int i=1; i<=bins; i++)
......@@ -136,7 +136,7 @@ namespace MISCMATHS {
if(i>2)
{
val+=kernel(3)*(histogram(i-2));
norm+=kernel(3);
norm+=kernel(3);
}
if(i<bins)
{
......@@ -146,7 +146,7 @@ namespace MISCMATHS {
if(i<bins-1)
{
val+=kernel(3)*(histogram(i+2));
norm+=kernel(3);
norm+=kernel(3);
}
val/=norm;
......@@ -188,7 +188,7 @@ namespace MISCMATHS {
void Histogram::match(Histogram &target){
int bin, newbin;
double cdfval, val,dist;
ColumnVector CDF_ref;
......@@ -200,52 +200,52 @@ namespace MISCMATHS {
olddata=sourceData;
histnew=0; dist=0;
if(exclusion.Nrows()==0){cout << " Histogram::match no excl " << endl; exclusion.ReSize(sourceData.Nrows()); exclusion=1;}
float binwidthtarget=(target.getHistMax() - target.getHistMin())/target.getNumBins();
for (int i=1;i<=sourceData.Nrows();i++){
if(exclusion(i)>0){
bin=getBin(sourceData(i));
cdfval=CDF(bin);
newbin=1;
val=sourceData(i)-histMin;
if(bin==target.getNumBins()){
newbin=bin;
}else if( (cdfval- CDF_ref(bin)) > 1e-20){
for (int j=bin;j<target.getNumBins();j++){
if((cdfval - CDF_ref(j)) > 1e-20 && (cdfval - CDF_ref(j+1)) <= 1e-20 ){
newbin=j+1;
dist=(cdfval-CDF_ref(j))/(CDF_ref(j+1)-CDF_ref(j)); // find distance between current bins as proportion of bin width
break;
}
}
}
else{
//search backwards
int j=bin-1;
while (j>0){
// if(sourceData(i)<1.2 && sourceData(i)>1)
if((cdfval - CDF_ref(1)) < -1e-20 && fabs(CDF_ref(bin) - CDF_ref(1)) < 1e-20){newbin=j+1; break; }
else if ((cdfval - CDF_ref(1)) > -1e-20 && fabs(CDF_ref(bin) - CDF_ref(1)) < 1e-20){newbin=j+1; break; }
else if((cdfval - CDF_ref(j)) > 1e-20 && (cdfval - CDF_ref(j+1)) <= 1e-20 ){
newbin=j+1;
dist=(cdfval-CDF_ref(j))/(CDF_ref(j+1)-CDF_ref(j)); // find distance between current bins as proportion of bin width
break;
}
......@@ -253,51 +253,17 @@ namespace MISCMATHS {
}
//search forwards
}
sourceData(i)=target.getHistMin() + (newbin-1)*binwidthtarget+dist*binwidthtarget;
histnew(newbin)++;
if(sourceData(i) < target.getHistMin()) sourceData(i) = target.getHistMin();
if(sourceData(i) > target.getHistMax()) sourceData(i) = target.getHistMax();
}
}
}
}
kernel.cc
+ 25
− 26
View file @ e63ee18d
......@@ -13,7 +13,7 @@ using namespace NEWMAT;
namespace MISCMATHS {
set<kernelstorage*, kernelstorage::comparer> kernel::existingkernels;
set<kernelstorage*, kernelstorage::comparer> kernel::existingkernels;
//////// Support function /////////
......@@ -28,42 +28,42 @@ namespace MISCMATHS {
dn -= n;
if (n>(kernel.Nrows()-1)) return 0.0;
if (n<1) return 0.0;
return kernel(n)*(1.0-dn) + kernel(n+1)*dn;
}
inline bool in_bounds(const ColumnVector& data, int index)
{ return ( (index>=1) && (index<=data.Nrows())); }
inline bool in_bounds(const ColumnVector& data, float index)
{ return ( ((int)floor(index)>=1) && ((int)ceil(index)<=data.Nrows())); }
float sincfn(float x)
{
if (fabs(x)<1e-7) { return 1.0-fabs(x); }
float y=M_PI*x;
return sin(y)/y;
}
float hanning(float x, int w)
{ // w is half-width
if (fabs(x)>w)
if (fabs(x)>w)
return 0.0;
else
return (0.5 + 0.5 *cos(M_PI*x/w));
}
float blackman(float x, int w)
{ // w is half-width
if (fabs(x)>w)
if (fabs(x)>w)
return 0.0;
else
return (0.42 + 0.5 *cos(M_PI*x/w) + 0.08*cos(2.0*M_PI*x/w));
}
float rectangular(float x, int w)
{ // w is half-width
if (fabs(x)>w)
if (fabs(x)>w)
return 0.0;
else
return 1.0;
......@@ -94,8 +94,8 @@ namespace MISCMATHS {
return ker;
}
kernel sinckernel(const string& sincwindowtype, int w, int nstore)
kernel sinckernel(const string& sincwindowtype, int w, int nstore)
{
kernel sinck;
sinck = sinckernel(sincwindowtype,w,w,w,nstore);
......@@ -119,7 +119,7 @@ namespace MISCMATHS {
kx = sinckernel1D(sincwindowtype,wx,nstore);
ky = sinckernel1D(sincwindowtype,wy,nstore);
kz = sinckernel1D(sincwindowtype,wz,nstore);
sinckern.setkernel(kx,ky,kz,hwx,hwy,hwz);
return sinckern;
}
......@@ -140,13 +140,13 @@ namespace MISCMATHS {
return extrapval;
}
// basic trilinear call
float interpolate_1d(const ColumnVector& data, const float index)
{
float interpval;
int low_bound = (int)floor(index);
int high_bound = (int)ceil(index);
int high_bound = (int)ceil(index);
if (in_bounds(data, index))
interpval = data(low_bound) + (index - low_bound)*(data(high_bound) - data(low_bound));
......@@ -156,7 +156,7 @@ namespace MISCMATHS {
return interpval;
}
//////// Spline Support /////////
float hermiteinterpolation_1d(const ColumnVector& data, int p1, int p4, float t)
......@@ -165,7 +165,7 @@ namespace MISCMATHS {
// inputs: points P_1, P_4; tangents R_1, R_4; interpolation index t;
float retval, r1 = 0.0, r4 = 0.0;
if (!in_bounds(data,p1) || !in_bounds(data,p4)) {
cerr << "Hermite Interpolation - ERROR: One or more indicies lie outside the data range. Returning ZERO" << endl;
retval = 0.0;
......@@ -176,17 +176,17 @@ namespace MISCMATHS {
retval = data(p1);
} else if (t == 1.0) {
retval = data(p4);
*/
*/
} else {
r1 = 0.5 * (extrapolate_1d(data, p1) - extrapolate_1d(data, p1 - 1)) + 0.5 * (extrapolate_1d(data, p1 + 1) - extrapolate_1d(data, p1));// tangent @ P_1
r4 = 0.5 * (extrapolate_1d(data, p4) - extrapolate_1d(data, p4 - 1)) + 0.5 * (extrapolate_1d(data, p4 + 1) - extrapolate_1d(data, p4));// tangent @ P_4
float t2 = t*t; float t3 = t2*t;
retval = (2*t3 - 3*t2 + 1)*data(p1) + (-2*t3 + 3*t2)*data(p4) + (t3 - 2*t2 + t)*r1 + (t3 - t2)*r4;
}
// cerr << "p1, p4, t, r1, r4 = " << p1 << ", " << p4 << ", " << t << ", " << r1 << ", " << r4 << endl;
return retval;
}
......@@ -199,14 +199,14 @@ namespace MISCMATHS {
// kernel half-width (i.e. range is +/- w)
int ix0;
ix0 = (int) floor(index);
int wx(widthx);
vector<float> storex(2*wx+1);
for (int d=-wx; d<=wx; d++)
for (int d=-wx; d<=wx; d++)
storex[d+wx] = kernelval((index-ix0+d),wx,userkernel);
float convsum=0.0, interpval=0.0, kersum=0.0;
int xj;
for (int x1=ix0-wx; x1<=ix0+wx; x1++) {
if (in_bounds(data, x1)) {
......@@ -224,7 +224,7 @@ namespace MISCMATHS {
}
return interpval;
}
////// Kernel wrappers //////
float kernelinterpolation_1d(const ColumnVector& data, float index)
......@@ -239,4 +239,3 @@ namespace MISCMATHS {
return kernelinterpolation_1d(data.t(), index, userkernel, 7);
}
}
optimise.cc
+ 21
− 29
View file @ e63ee18d
......@@ -28,7 +28,7 @@ namespace MISCMATHS {
using std::exp;
using std::log;
bool estquadmin(float &xnew, float x1, float xmid, float x2,
bool estquadmin(float &xnew, float x1, float xmid, float x2,
float y1, float ymid, float y2)
{
// Finds the estimated quadratic minimum's position
......@@ -64,7 +64,7 @@ namespace MISCMATHS {
}
return xnew;
}
float nextpt(float x1, float xmid, float x2, float y1, float ymid, float y2)
......@@ -82,10 +82,10 @@ namespace MISCMATHS {
return xnew;
}
void findinitialbound(float &x1, float &xmid, float &x2,
float &y1, float &ymid, float &y2,
void findinitialbound(float &x1, float &xmid, float &x2,
float &y1, float &ymid, float &y2,
float (*func)(const ColumnVector &),
const ColumnVector &unitdir, const ColumnVector &pt)
{
......@@ -109,15 +109,15 @@ namespace MISCMATHS {
y2 = (*func)(x2*unitdir + pt);
while (ymid > y2) { // note: must maintain y1 >= ymid
// cout << " <" << Min(x1,x2) << "," << xmid
// cout << " <" << Min(x1,x2) << "," << xmid
// << "," << Max(x1,x2) << ">" << endl;
maxx2 = xmid + maxextrap*(x2 - xmid);
quadok = estquadmin(newx2,x1,xmid,x2,y1,ymid,y2);
if ((!quadok) || ((newx2 - x1)*dir<0) || ((newx2 - maxx2)*dir>0)) {
newx2 = xmid + extrapolationfactor*(x2-x1);
}
newy2 = (*func)(newx2*unitdir + pt);
if ((newx2 - xmid)*(newx2 - x1)<0) { // newx2 is between x1 and xmid
......@@ -141,19 +141,19 @@ namespace MISCMATHS {
x2 = newx2; y2 = newy2;
}
}
}
if ( (y2<ymid) || (y1<ymid) ) {
cerr << "findinitialbound failed to bracket: current triplet is" << endl;
}
}
float optimise1d(ColumnVector &pt, const ColumnVector dir,
const ColumnVector tol, int &iterations_done,
float optimise1d(ColumnVector &pt, const ColumnVector dir,
const ColumnVector tol, int &iterations_done,
float (*func)(const ColumnVector &), int max_iter,
float &init_value, float boundguess)
float &init_value, float boundguess)
{
// Golden Search Routine
// Must pass in the direction vector in N-space (dir), the initial
......@@ -247,10 +247,10 @@ namespace MISCMATHS {
float optimise(ColumnVector &pt, int numopt, const ColumnVector &tol,
float (*func)(const ColumnVector &), int &iterations_done,
int max_iter, const ColumnVector& boundguess,
float optimise(ColumnVector &pt, int numopt, const ColumnVector &tol,
float (*func)(const ColumnVector &), int &iterations_done,
int max_iter, const ColumnVector& boundguess,
const string type)
{
// Note that numopt can be less than pt.Nrows() - e.g. 6 dof optimisation
......@@ -291,12 +291,12 @@ namespace MISCMATHS {
if (avtol < 1.0) break;
// if continuing then change the directions if using Powell's method
if (type=="powell")
if (type=="powell")
{
// find direction of maximal change
int bestm=1;
for (int m=1; m<=numopt; m++) {
if (deltaf(m)<deltaf(bestm)) bestm=m;
for (int m=1; m<=numopt; m++) {
if (deltaf(m)<deltaf(bestm)) bestm=m;
}
fend=fval;
fextrap=(*func)(initpt + 2*(pt-initpt));
......@@ -317,7 +317,7 @@ namespace MISCMATHS {
}
}
}
}
// cout << endl << "Major iterations = " << it << endl;
iterations_done = littot;
......@@ -326,11 +326,3 @@ namespace MISCMATHS {
}
rungekutta.cc
+ 25
− 25
View file @ e63ee18d
......@@ -14,34 +14,34 @@ using namespace NEWMAT;
namespace MISCMATHS {
void rk(ColumnVector& ret, const ColumnVector& y, const ColumnVector& dy, float x, float h, const Derivative& deriv,const ColumnVector& paramvalues)
{
Tracer tr("rk");
{
Tracer tr("rk");
float hh=h*0.5;
float xh=x+hh;
//first step
ColumnVector yt=y+hh*dy;
//second step
ColumnVector dyt = deriv.evaluate(xh,yt,paramvalues);
yt=y+hh*dyt;
//third step
ColumnVector dym = deriv.evaluate(xh,yt,paramvalues);
yt=y+h*dym;
dym=dym+dyt;
//fourth step
dyt = deriv.evaluate(x+h,yt,paramvalues);
//addup
ret = y+h*(dy+dyt+2*dym)/6;
}
void rkqc(ColumnVector& y, float& x, float& hnext, ColumnVector& dy, float htry, float eps, const Derivative& deriv,const ColumnVector& paramvalues)
{
Tracer tr("rkqc");
Tracer tr("rkqc");
float xsav = x;
ColumnVector dysav = dy;
......@@ -53,12 +53,12 @@ void rkqc(ColumnVector& y, float& x, float& hnext, ColumnVector& dy, float htry,
while(true)
{
// take 2 1/2 step sizes
// first 1/2 step
float hh=h*0.5;
rk(ytemp,ysav,dysav,xsav,hh,deriv,paramvalues);
// second 1/2 step
x=xsav+hh;
dy = deriv.evaluate(x,ytemp,paramvalues);
......@@ -68,32 +68,32 @@ void rkqc(ColumnVector& y, float& x, float& hnext, ColumnVector& dy, float htry,
// take large step size
rk(ytemp,ysav,dysav,xsav,h,deriv,paramvalues);
// eval accuracy
float errmax = 0.0;
for(int i=1; i<=y.Nrows(); i++)
{
//errmax=max(abs((y-ytemp)./y));
float tmp = fabs((y(i)-ytemp(i))/y(i));
if(tmp > errmax) errmax = tmp;
}
errmax=errmax/eps;
if(errmax <=1.0)
if(errmax <=1.0)
{
// step OK, compute step size for next step
hdid=h;
if(errmax>6e-4)
hnext=h*std::exp(-0.2*std::log(errmax));
else
hnext=4*h;
break;
}
else
else
{
// step too large,
h=h*std::exp(-0.25*std::log(errmax));
......@@ -105,7 +105,7 @@ void rkqc(ColumnVector& y, float& x, float& hnext, ColumnVector& dy, float htry,
void runge_kutta(Matrix& yp, ColumnVector& xp, ColumnVector& hp, const ColumnVector& ystart, float x1, float x2, float eps, float hmin, const Derivative& deriv,const ColumnVector& paramvalues)
{
Tracer tr("runge_kutta");
Tracer tr("runge_kutta");
int MAXSTEP=1000;
......@@ -124,20 +124,20 @@ void runge_kutta(Matrix& yp, ColumnVector& xp, ColumnVector& hp, const ColumnVec
yp = 0;
int kout=1;
ColumnVector dy;
for(int k=1; k <= MAXSTEP; k++)
{
{
dy = deriv.evaluate(x,y,paramvalues);
// store results:
xp(kout)=x;
yp.Row(kout)=y;
hp(kout)=h;
kout=kout+1;
// stop overshoot of step past x2:
if((x+h-x2)*(x+h-x1)>0) h=x2-x;
......@@ -150,7 +150,7 @@ void runge_kutta(Matrix& yp, ColumnVector& xp, ColumnVector& hp, const ColumnVec
yp.Row(kout)=y;
hp(kout)=h;
//kout=kout+1;
xp = xp.Rows(1,kout);
yp = yp.Rows(1,kout);
......@@ -160,8 +160,8 @@ void runge_kutta(Matrix& yp, ColumnVector& xp, ColumnVector& hp, const ColumnVec
{
if(hnext<=hmin) cerr << "step size too small" << endl;
h=hnext;
}
}
}
cerr << "too many steps" << endl;
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment