Skip to content
GitLab
Explore
Sign in
Primary navigation
Search or go to…
Project
M
miscmaths
Manage
Activity
Members
Labels
Plan
Issues
Issue boards
Milestones
Wiki
Code
Merge requests
Repository
Branches
Commits
Tags
Repository graph
Compare revisions
Snippets
Build
Pipelines
Jobs
Pipeline schedules
Artifacts
Deploy
Releases
Container Registry
Model registry
Operate
Environments
Monitor
Incidents
Analyze
Value stream analytics
Contributor analytics
CI/CD analytics
Repository analytics
Model experiments
Help
Help
Support
GitLab documentation
Compare GitLab plans
Community forum
Contribute to GitLab
Provide feedback
Keyboard shortcuts
?
Snippets
Groups
Projects
Show more breadcrumbs
FSL
miscmaths
Commits
d86209cb
Commit
d86209cb
authored
19 years ago
by
David Flitney
Browse files
Options
Downloads
Patches
Plain Diff
Changes to support MS Visual Studio brokeness
parent
fe86bcb1
No related branches found
No related tags found
No related merge requests found
Changes
1
Hide whitespace changes
Inline
Side-by-side
Showing
1 changed file
miscmaths.cc
+39
-36
39 additions, 36 deletions
miscmaths.cc
with
39 additions
and
36 deletions
miscmaths.cc
+
39
−
36
View file @
d86209cb
...
@@ -8,26 +8,29 @@
...
@@ -8,26 +8,29 @@
// Miscellaneous maths functions
// Miscellaneous maths functions
#define NOMINMAX
#include
<cstdlib>
#include
<cmath>
#include
"miscmaths.h"
#include
"miscmaths.h"
#include
"miscprob.h"
#include
"miscprob.h"
#include
"stdlib.h"
#include
"stdlib.h"
#include
"newmatio.h"
#include
"newmatio.h"
using
namespace
std
;
namespace
MISCMATHS
{
namespace
MISCMATHS
{
// The following lines are ignored by the current SGI compiler
// The following lines are ignored by the current SGI compiler
// (version egcs-2.91.57)
// (version egcs-2.91.57)
// A temporary fix of including the std:: in front of all abs() etc
// A temporary fix of including the std:: in front of all abs() etc
// has been done for now
// has been done for now
/*
using std::abs;
using std::abs;
using std::sqrt;
using std::sqrt;
using std::exp;
using std::exp;
using std::log;
using std::log;
// using std::pow;
// using std::pow;
using std::atan2;
using std::atan2;
*/
string
size
(
const
Matrix
&
mat
)
string
size
(
const
Matrix
&
mat
)
{
{
...
@@ -935,7 +938,7 @@ namespace MISCMATHS {
...
@@ -935,7 +938,7 @@ namespace MISCMATHS {
for
(
int
i
=
1
;
i
<=
3
;
i
++
)
{
params
(
i
+
3
)
=
transl
(
i
);
}
for
(
int
i
=
1
;
i
<=
3
;
i
++
)
{
params
(
i
+
3
)
=
transl
(
i
);
}
ColumnVector
rotparams
(
3
);
ColumnVector
rotparams
(
3
);
(
*
rotmat2params
)(
rotparams
,
rotmat
);
(
*
rotmat2params
)(
rotparams
,
rotmat
);
for
(
int
i
=
1
;
i
<=
3
;
i
++
)
{
params
(
i
)
=
rotparams
(
i
);
}
for
(
int
i
i
=
1
;
i
i
<=
3
;
i
i
++
)
{
params
(
i
i
)
=
rotparams
(
i
i
);
}
return
0
;
return
0
;
}
}
...
@@ -1067,7 +1070,7 @@ float median(const ColumnVector& x)
...
@@ -1067,7 +1070,7 @@ float median(const ColumnVector& x)
void
cart2sph
(
const
ColumnVector
&
dir
,
float
&
th
,
float
&
ph
)
void
cart2sph
(
const
ColumnVector
&
dir
,
float
&
th
,
float
&
ph
)
{
{
float
mag
=
sqrt
(
dir
(
1
)
*
dir
(
1
)
+
dir
(
2
)
*
dir
(
2
)
+
dir
(
3
)
*
dir
(
3
));
float
mag
=
std
::
sqrt
(
dir
(
1
)
*
dir
(
1
)
+
dir
(
2
)
*
dir
(
2
)
+
dir
(
3
)
*
dir
(
3
));
if
(
mag
==
0
){
if
(
mag
==
0
){
ph
=
M_PI
/
2
;
ph
=
M_PI
/
2
;
th
=
M_PI
/
2
;
th
=
M_PI
/
2
;
...
@@ -1076,13 +1079,13 @@ void cart2sph(const ColumnVector& dir, float& th, float& ph)
...
@@ -1076,13 +1079,13 @@ void cart2sph(const ColumnVector& dir, float& th, float& ph)
if
(
dir
(
1
)
==
0
&&
dir
(
2
)
>=
0
)
ph
=
M_PI
/
2
;
if
(
dir
(
1
)
==
0
&&
dir
(
2
)
>=
0
)
ph
=
M_PI
/
2
;
else
if
(
dir
(
1
)
==
0
&&
dir
(
2
)
<
0
)
ph
=-
M_PI
/
2
;
else
if
(
dir
(
1
)
==
0
&&
dir
(
2
)
<
0
)
ph
=-
M_PI
/
2
;
else
if
(
dir
(
1
)
>
0
)
ph
=
atan
(
dir
(
2
)
/
dir
(
1
));
else
if
(
dir
(
1
)
>
0
)
ph
=
std
::
atan
(
dir
(
2
)
/
dir
(
1
));
else
if
(
dir
(
2
)
>
0
)
ph
=
atan
(
dir
(
2
)
/
dir
(
1
))
+
M_PI
;
else
if
(
dir
(
2
)
>
0
)
ph
=
std
::
atan
(
dir
(
2
)
/
dir
(
1
))
+
M_PI
;
else
ph
=
atan
(
dir
(
2
)
/
dir
(
1
))
-
M_PI
;
else
ph
=
std
::
atan
(
dir
(
2
)
/
dir
(
1
))
-
M_PI
;
if
(
dir
(
3
)
==
0
)
th
=
M_PI
/
2
;
if
(
dir
(
3
)
==
0
)
th
=
M_PI
/
2
;
else
if
(
dir
(
3
)
>
0
)
th
=
atan
(
sqrt
(
dir
(
1
)
*
dir
(
1
)
+
dir
(
2
)
*
dir
(
2
))
/
dir
(
3
));
else
if
(
dir
(
3
)
>
0
)
th
=
std
::
atan
(
std
::
sqrt
(
dir
(
1
)
*
dir
(
1
)
+
dir
(
2
)
*
dir
(
2
))
/
dir
(
3
));
else
th
=
atan
(
sqrt
(
dir
(
1
)
*
dir
(
1
)
+
dir
(
2
)
*
dir
(
2
))
/
dir
(
3
))
+
M_PI
;
else
th
=
std
::
atan
(
std
::
sqrt
(
dir
(
1
)
*
dir
(
1
)
+
dir
(
2
)
*
dir
(
2
))
/
dir
(
3
))
+
M_PI
;
}
}
}
}
...
@@ -1099,7 +1102,7 @@ void cart2sph(const Matrix& dir,ColumnVector& th,ColumnVector& ph)
...
@@ -1099,7 +1102,7 @@ void cart2sph(const Matrix& dir,ColumnVector& th,ColumnVector& ph)
}
}
for
(
int
i
=
1
;
i
<=
dir
.
Ncols
();
i
++
)
{
for
(
int
i
=
1
;
i
<=
dir
.
Ncols
();
i
++
)
{
float
mag
=
sqrt
(
dir
(
1
,
i
)
*
dir
(
1
,
i
)
+
dir
(
2
,
i
)
*
dir
(
2
,
i
)
+
dir
(
3
,
i
)
*
dir
(
3
,
i
));
float
mag
=
std
::
sqrt
(
dir
(
1
,
i
)
*
dir
(
1
,
i
)
+
dir
(
2
,
i
)
*
dir
(
2
,
i
)
+
dir
(
3
,
i
)
*
dir
(
3
,
i
));
if
(
mag
==
0
){
if
(
mag
==
0
){
ph
(
i
)
=
M_PI
/
2
;
ph
(
i
)
=
M_PI
/
2
;
th
(
i
)
=
M_PI
/
2
;
th
(
i
)
=
M_PI
/
2
;
...
@@ -1107,13 +1110,13 @@ void cart2sph(const Matrix& dir,ColumnVector& th,ColumnVector& ph)
...
@@ -1107,13 +1110,13 @@ void cart2sph(const Matrix& dir,ColumnVector& th,ColumnVector& ph)
else
{
else
{
if
(
dir
(
1
,
i
)
==
0
&&
dir
(
2
,
i
)
>=
0
)
ph
(
i
)
=
M_PI
/
2
;
if
(
dir
(
1
,
i
)
==
0
&&
dir
(
2
,
i
)
>=
0
)
ph
(
i
)
=
M_PI
/
2
;
else
if
(
dir
(
1
,
i
)
==
0
&&
dir
(
2
,
i
)
<
0
)
ph
(
i
)
=-
M_PI
/
2
;
else
if
(
dir
(
1
,
i
)
==
0
&&
dir
(
2
,
i
)
<
0
)
ph
(
i
)
=-
M_PI
/
2
;
else
if
(
dir
(
1
,
i
)
>
0
)
ph
(
i
)
=
atan
(
dir
(
2
,
i
)
/
dir
(
1
,
i
));
else
if
(
dir
(
1
,
i
)
>
0
)
ph
(
i
)
=
std
::
atan
(
dir
(
2
,
i
)
/
dir
(
1
,
i
));
else
if
(
dir
(
2
,
i
)
>
0
)
ph
(
i
)
=
atan
(
dir
(
2
,
i
)
/
dir
(
1
,
i
))
+
M_PI
;
else
if
(
dir
(
2
,
i
)
>
0
)
ph
(
i
)
=
std
::
atan
(
dir
(
2
,
i
)
/
dir
(
1
,
i
))
+
M_PI
;
else
ph
(
i
)
=
atan
(
dir
(
2
,
i
)
/
dir
(
1
,
i
))
-
M_PI
;
else
ph
(
i
)
=
std
::
atan
(
dir
(
2
,
i
)
/
dir
(
1
,
i
))
-
M_PI
;
if
(
dir
(
3
,
i
)
==
0
)
th
(
i
)
=
M_PI
/
2
;
if
(
dir
(
3
,
i
)
==
0
)
th
(
i
)
=
M_PI
/
2
;
else
if
(
dir
(
3
,
i
)
>
0
)
th
(
i
)
=
atan
(
sqrt
(
dir
(
1
,
i
)
*
dir
(
1
,
i
)
+
dir
(
2
,
i
)
*
dir
(
2
,
i
))
/
dir
(
3
,
i
));
else
if
(
dir
(
3
,
i
)
>
0
)
th
(
i
)
=
std
::
atan
(
std
::
sqrt
(
dir
(
1
,
i
)
*
dir
(
1
,
i
)
+
dir
(
2
,
i
)
*
dir
(
2
,
i
))
/
dir
(
3
,
i
));
else
th
(
i
)
=
atan
(
sqrt
(
dir
(
1
,
i
)
*
dir
(
1
,
i
)
+
dir
(
2
,
i
)
*
dir
(
2
,
i
))
/
dir
(
3
,
i
))
+
M_PI
;
else
th
(
i
)
=
std
::
atan
(
std
::
sqrt
(
dir
(
1
,
i
)
*
dir
(
1
,
i
)
+
dir
(
2
,
i
)
*
dir
(
2
,
i
))
/
dir
(
3
,
i
))
+
M_PI
;
}
}
}
}
...
@@ -1145,7 +1148,7 @@ ReturnMatrix repmat(const Matrix &mat, const int rows, const int cols)
...
@@ -1145,7 +1148,7 @@ ReturnMatrix repmat(const Matrix &mat, const int rows, const int cols)
Matrix
res
=
mat
;
Matrix
res
=
mat
;
for
(
int
ctr
=
1
;
ctr
<
cols
;
ctr
++
){
res
|=
mat
;}
for
(
int
ctr
=
1
;
ctr
<
cols
;
ctr
++
){
res
|=
mat
;}
Matrix
tmpres
=
res
;
Matrix
tmpres
=
res
;
for
(
int
ctr
=
1
;
ctr
<
rows
;
ctr
++
){
res
&=
tmpres
;}
for
(
int
ctr
1
=
1
;
ctr
1
<
rows
;
ctr
1
++
){
res
&=
tmpres
;}
res
.
Release
();
res
.
Release
();
return
res
;
return
res
;
}
}
...
@@ -1395,8 +1398,8 @@ ReturnMatrix stdev(const Matrix& mat, const int dim)
...
@@ -1395,8 +1398,8 @@ ReturnMatrix stdev(const Matrix& mat, const int dim)
ReturnMatrix
gt
(
const
Matrix
&
mat1
,
const
Matrix
&
mat2
)
ReturnMatrix
gt
(
const
Matrix
&
mat1
,
const
Matrix
&
mat2
)
{
{
int
ctrcol
=
std
::
min
(
mat1
.
Ncols
(),
mat2
.
Ncols
());
int
ctrcol
=
X
min
(
mat1
.
Ncols
(),
mat2
.
Ncols
());
int
ctrrow
=
std
::
min
(
mat1
.
Nrows
(),
mat2
.
Nrows
());
int
ctrrow
=
X
min
(
mat1
.
Nrows
(),
mat2
.
Nrows
());
Matrix
res
(
ctrrow
,
ctrcol
);
Matrix
res
(
ctrrow
,
ctrcol
);
res
=
0.0
;
res
=
0.0
;
...
@@ -1415,8 +1418,8 @@ ReturnMatrix gt(const Matrix& mat1,const Matrix& mat2)
...
@@ -1415,8 +1418,8 @@ ReturnMatrix gt(const Matrix& mat1,const Matrix& mat2)
ReturnMatrix
lt
(
const
Matrix
&
mat1
,
const
Matrix
&
mat2
)
ReturnMatrix
lt
(
const
Matrix
&
mat1
,
const
Matrix
&
mat2
)
{
{
int
ctrcol
=
std
::
min
(
mat1
.
Ncols
(),
mat2
.
Ncols
());
int
ctrcol
=
X
min
(
mat1
.
Ncols
(),
mat2
.
Ncols
());
int
ctrrow
=
std
::
min
(
mat1
.
Nrows
(),
mat2
.
Nrows
());
int
ctrrow
=
X
min
(
mat1
.
Nrows
(),
mat2
.
Nrows
());
Matrix
res
(
ctrrow
,
ctrcol
);
Matrix
res
(
ctrrow
,
ctrcol
);
res
=
0.0
;
res
=
0.0
;
...
@@ -1435,8 +1438,8 @@ ReturnMatrix lt(const Matrix& mat1,const Matrix& mat2)
...
@@ -1435,8 +1438,8 @@ ReturnMatrix lt(const Matrix& mat1,const Matrix& mat2)
ReturnMatrix
geqt
(
const
Matrix
&
mat1
,
const
Matrix
&
mat2
)
ReturnMatrix
geqt
(
const
Matrix
&
mat1
,
const
Matrix
&
mat2
)
{
{
int
ctrcol
=
std
::
min
(
mat1
.
Ncols
(),
mat2
.
Ncols
());
int
ctrcol
=
X
min
(
mat1
.
Ncols
(),
mat2
.
Ncols
());
int
ctrrow
=
std
::
min
(
mat1
.
Nrows
(),
mat2
.
Nrows
());
int
ctrrow
=
X
min
(
mat1
.
Nrows
(),
mat2
.
Nrows
());
Matrix
res
(
ctrrow
,
ctrcol
);
Matrix
res
(
ctrrow
,
ctrcol
);
res
=
0.0
;
res
=
0.0
;
...
@@ -1455,8 +1458,8 @@ ReturnMatrix geqt(const Matrix& mat1,const Matrix& mat2)
...
@@ -1455,8 +1458,8 @@ ReturnMatrix geqt(const Matrix& mat1,const Matrix& mat2)
ReturnMatrix
leqt
(
const
Matrix
&
mat1
,
const
Matrix
&
mat2
)
ReturnMatrix
leqt
(
const
Matrix
&
mat1
,
const
Matrix
&
mat2
)
{
{
int
ctrcol
=
std
::
min
(
mat1
.
Ncols
(),
mat2
.
Ncols
());
int
ctrcol
=
X
min
(
mat1
.
Ncols
(),
mat2
.
Ncols
());
int
ctrrow
=
std
::
min
(
mat1
.
Nrows
(),
mat2
.
Nrows
());
int
ctrrow
=
X
min
(
mat1
.
Nrows
(),
mat2
.
Nrows
());
Matrix
res
(
ctrrow
,
ctrcol
);
Matrix
res
(
ctrrow
,
ctrcol
);
res
=
0.0
;
res
=
0.0
;
...
@@ -1475,8 +1478,8 @@ ReturnMatrix leqt(const Matrix& mat1,const Matrix& mat2)
...
@@ -1475,8 +1478,8 @@ ReturnMatrix leqt(const Matrix& mat1,const Matrix& mat2)
ReturnMatrix
eq
(
const
Matrix
&
mat1
,
const
Matrix
&
mat2
)
ReturnMatrix
eq
(
const
Matrix
&
mat1
,
const
Matrix
&
mat2
)
{
{
int
ctrcol
=
std
::
min
(
mat1
.
Ncols
(),
mat2
.
Ncols
());
int
ctrcol
=
X
min
(
mat1
.
Ncols
(),
mat2
.
Ncols
());
int
ctrrow
=
std
::
min
(
mat1
.
Nrows
(),
mat2
.
Nrows
());
int
ctrrow
=
X
min
(
mat1
.
Nrows
(),
mat2
.
Nrows
());
Matrix
res
(
ctrrow
,
ctrcol
);
Matrix
res
(
ctrrow
,
ctrcol
);
res
=
0.0
;
res
=
0.0
;
...
@@ -1495,8 +1498,8 @@ ReturnMatrix eq(const Matrix& mat1,const Matrix& mat2)
...
@@ -1495,8 +1498,8 @@ ReturnMatrix eq(const Matrix& mat1,const Matrix& mat2)
ReturnMatrix
neq
(
const
Matrix
&
mat1
,
const
Matrix
&
mat2
)
ReturnMatrix
neq
(
const
Matrix
&
mat1
,
const
Matrix
&
mat2
)
{
{
int
ctrcol
=
std
::
min
(
mat1
.
Ncols
(),
mat2
.
Ncols
());
int
ctrcol
=
X
min
(
mat1
.
Ncols
(),
mat2
.
Ncols
());
int
ctrrow
=
std
::
min
(
mat1
.
Nrows
(),
mat2
.
Nrows
());
int
ctrrow
=
X
min
(
mat1
.
Nrows
(),
mat2
.
Nrows
());
Matrix
res
(
ctrrow
,
ctrcol
);
Matrix
res
(
ctrrow
,
ctrcol
);
res
=
0.0
;
res
=
0.0
;
...
@@ -1675,7 +1678,7 @@ void element_mod_n(Matrix& Mat,double n)
...
@@ -1675,7 +1678,7 @@ void element_mod_n(Matrix& Mat,double n)
int
nextpow2
(
int
n
)
int
nextpow2
(
int
n
)
{
{
return
(
int
)
pow
(
2
,
ceil
(
log
(
float
(
n
))
/
log
(
float
(
2
))));
return
(
int
)
std
::
pow
(
2
,
ceil
(
std
::
log
(
float
(
n
))
/
std
::
log
(
float
(
2
))));
}
}
void
xcorr
(
const
Matrix
&
p_ts
,
Matrix
&
ret
,
int
lag
,
int
p_zeropad
)
void
xcorr
(
const
Matrix
&
p_ts
,
Matrix
&
ret
,
int
lag
,
int
p_zeropad
)
...
@@ -1718,10 +1721,10 @@ void xcorr(const Matrix& p_ts, Matrix& ret, int lag, int p_zeropad)
...
@@ -1718,10 +1721,10 @@ void xcorr(const Matrix& p_ts, Matrix& ret, int lag, int p_zeropad)
float
varx
=
var
(
x
.
Rows
(
1
,
sizeTS
)).
AsScalar
();
float
varx
=
var
(
x
.
Rows
(
1
,
sizeTS
)).
AsScalar
();
ret
.
Column
(
i
)
=
realifft
.
Rows
(
1
,
lag
);
ret
.
Column
(
i
)
=
realifft
.
Rows
(
1
,
lag
);
for
(
int
j
=
1
;
j
<=
lag
-
1
;
j
++
)
for
(
int
j
j
=
1
;
j
j
<=
lag
-
1
;
j
j
++
)
{
{
// Correction to make autocorr unbiased and normalised
// Correction to make autocorr unbiased and normalised
ret
(
j
,
i
)
=
ret
(
j
,
i
)
/
((
sizeTS
-
j
)
*
varx
);
ret
(
j
j
,
i
)
=
ret
(
j
j
,
i
)
/
((
sizeTS
-
j
j
)
*
varx
);
}
}
}
}
}
}
...
@@ -1749,15 +1752,15 @@ void detrend(Matrix& p_ts, int p_level)
...
@@ -1749,15 +1752,15 @@ void detrend(Matrix& p_ts, int p_level)
for
(
int
t
=
1
;
t
<=
sizeTS
;
t
++
)
for
(
int
t
=
1
;
t
<=
sizeTS
;
t
++
)
{
{
for
(
int
l
=
0
;
l
<=
p_level
;
l
++
)
for
(
int
l
=
0
;
l
<=
p_level
;
l
++
)
a
(
t
,
l
+
1
)
=
pow
((
float
)
t
/
sizeTS
,
l
);
a
(
t
,
l
+
1
)
=
std
::
pow
((
float
)
t
/
sizeTS
,
l
);
}
}
// Form residual forming matrix R:
// Form residual forming matrix R:
Matrix
R
=
Identity
(
sizeTS
)
-
a
*
pinv
(
a
);
Matrix
R
=
Identity
(
sizeTS
)
-
a
*
pinv
(
a
);
for
(
int
t
=
1
;
t
<=
sizeTS
;
t
++
)
for
(
int
t
t
=
1
;
t
t
<=
sizeTS
;
t
t
++
)
{
{
p_ts
.
Column
(
t
)
=
R
*
p_ts
.
Column
(
t
);
p_ts
.
Column
(
t
t
)
=
R
*
p_ts
.
Column
(
t
t
);
}
}
}
}
...
@@ -1993,7 +1996,7 @@ float csevl(const float x, const ColumnVector& cs, const int n)
...
@@ -1993,7 +1996,7 @@ float csevl(const float x, const ColumnVector& cs, const int n)
else
else
{
{
const
float
aux
=
csevl
(
8
/
(
Sqr
(
y
))
-
1
,
apsics
,
ntapsi
);
const
float
aux
=
csevl
(
8
/
(
Sqr
(
y
))
-
1
,
apsics
,
ntapsi
);
psi
=
log
(
x
)
-
0.5
/
x
+
aux
;
psi
=
std
::
log
(
x
)
-
0.5
/
x
+
aux
;
}
}
return
psi
;
return
psi
;
...
...
This diff is collapsed.
Click to expand it.
Preview
0%
Loading
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Save comment
Cancel
Please
register
or
sign in
to comment