sunwen/Aurora
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add column vector support to min, max, sum, prod.

Browse Source
This commit is contained in:
Krad
2023-04-26 14:58:18 +08:00
parent 6fd22ffc89
commit b4423b756e
3 changed files with 90 additions and 30 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/Function1D.h
View File

@@ -83,6 +83,8 @@ namespace Aurora {
* @return 查询结果
*/
Matrix polyval(const Matrix& aP, const Matrix& aX);
void nantoval(Matrix& aMatrix,double val);
};

91
src/Function2D.cpp
View File

@@ -117,18 +117,26 @@ Matrix Aurora::interpn(const Matrix& aX, const Matrix& aY, const Matrix& aV, con
}
Matrix Aurora::std(const Matrix &aMatrix) {
auto std = Aurora::malloc(aMatrix.getDimSize(1) * aMatrix.getDimSize(2));
if (aMatrix.getDimSize(2) > 1 || aMatrix.isComplex()) {
std::cerr
<< (aMatrix.getDimSize(2) > 1 ? "std() not support 3D data!" : "std() not support complex value type!")
<< std::endl;
return Matrix();
}
Matrix src = aMatrix.isComplex() ? Aurora::abs(aMatrix) : aMatrix;
int calc_size = src.getDimSize(0) == 1 ? src.getDimSize(1) : src.getDimSize(0);
auto std = Aurora::malloc(aMatrix.getDimSize(1));
for (int i = 0; i < src.getDimSize(1); ++i) {
double *p = src.getData() + i * src.getDimSize(0);
double mean = cblas_dasum(src.getDimSize(0), p, 1) / src.getDimSize(0);
vdSubI(src.getDimSize(0), p, 1, &mean, 0, p, 1);
vdSqr(src.getDimSize(0), p, p);
std[i] = cblas_dasum(src.getDimSize(0), p, 1) / (src.getDimSize(0) - 1);
double *p = src.getData() + i * calc_size;
double mean = cblas_dasum(calc_size, p, 1) / calc_size;
vdSubI(calc_size, p, 1, &mean, 0, p, 1);
vdSqr(calc_size, p, p);
std[i] = cblas_dasum(calc_size, p, 1) / (calc_size - 1);
}
vdSqrt(src.getDimSize(1), std, std);
return Matrix::New(std, 1, aMatrix.getDimSize(1), aMatrix.getDimSize(2));
return Matrix::New(std,1,aMatrix.getDimSize(1), aMatrix.getDimSize(2));
}
Matrix Aurora::min(const Matrix &aMatrix, FunctionDirection direction, long& rowIdx, long& colIdx) {
@@ -138,6 +146,10 @@ Matrix Aurora::min(const Matrix &aMatrix, FunctionDirection direction, long& row
<< std::endl;
return Matrix();
}
//针对向量行等于列
if (direction == Column && aMatrix.getDimSize(0)==1){
direction = All;
}
switch (direction)
{
case All: {
@@ -248,6 +260,10 @@ Matrix Aurora::max(const Matrix &aMatrix, FunctionDirection direction, long& row
<< std::endl;
return Matrix();
}
//针对向量行等于列
if (direction == Column && aMatrix.getDimSize(0)==1){
direction = All;
}
switch (direction)
{
case All:
@@ -293,6 +309,10 @@ Matrix Aurora::sum(const Matrix &aMatrix, FunctionDirection direction) {
<< std::endl;
return Matrix();
}
//针对向量行等于列
if (direction == Column && aMatrix.getDimSize(0)==1){
direction = Row;
}
if (aMatrix.isComplex()){
switch (direction)
{
@@ -360,6 +380,10 @@ Matrix Aurora::mean(const Matrix &aMatrix, FunctionDirection direction, bool aIn
<< std::endl;
return Matrix();
}
//针对向量行等于列
if (direction == Column && aMatrix.getDimSize(0)==1){
direction = All;
}
if (aIncludeNan){
switch (direction)
{
@@ -465,6 +489,7 @@ Matrix Aurora::sort(const Matrix &aMatrix) {
}
return sort(std::forward<Matrix &&>(aMatrix.deepCopy()));
}
Matrix Aurora::sort(Matrix &&aMatrix) {
if (aMatrix.getDimSize(2)>1 || aMatrix.isComplex()) {
std::cerr
@@ -532,17 +557,20 @@ Matrix Aurora::median(const Matrix &aMatrix) {
<< std::endl;
return Matrix();
}
Matrix sorted = sort(aMatrix);
Eigen::Map<Eigen::MatrixXd> srcM(sorted.getData(),sorted.getDimSize(0),sorted.getDimSize(1));
bool horVector = aMatrix.getDimSize(0)==1;
Matrix sorted = horVector?sortrows(aMatrix):sort(aMatrix);
int rows = horVector?sorted.getDimSize(1):sorted.getDimSize(0);
int cols = horVector?sorted.getDimSize(0):sorted.getDimSize(1);
Eigen::Map<Eigen::MatrixXd> srcM(sorted.getData(),rows,cols);
bool flag = aMatrix.getDimSize(0) % 2 == 1;
double* ret = malloc(aMatrix.getDimSize(1));
Eigen::Map<Eigen::VectorXd> retV(ret,aMatrix.getDimSize(1));
double* ret = malloc(cols);
Eigen::Map<Eigen::VectorXd> retV(ret,cols);
if (flag) {
retV = srcM.row(aMatrix.getDimSize(0)/2);
return Matrix::New(ret,1,aMatrix.getDimSize(1));
retV = srcM.row(rows/2);
return Matrix::New(ret,1,cols);
} else {
retV = (srcM.row(aMatrix.getDimSize(0)/2-1).array()+srcM.row(aMatrix.getDimSize(0)/2).array())/2;
return Matrix::New(ret,1,aMatrix.getDimSize(1));
retV = (srcM.row(rows/2-1).array()+srcM.row(rows/2).array())/2;
return Matrix::New(ret,1,cols);
}
}
@@ -667,43 +695,46 @@ Matrix Aurora::hilbert(const Matrix &aMatrix) {
vdMulI(aMatrix.getDimSize(0), p + 1, 2, h, 1, p + 1, 2);
}
auto result = ifft( x);
delete[] h;
free(h);
return result;
}
Matrix Aurora::prod(const Matrix &aMatrix) {
if (aMatrix.getDimSize(2) > 1 ) {
std::cerr<< "prod() not support 3D data!"
<< std::endl;
std::cerr<< "prod() not support 3D data!"<< std::endl;
return Matrix();
}
int row = aMatrix.getDimSize(0)==1?aMatrix.getDimSize(1):aMatrix.getDimSize(0);
int col = aMatrix.getDimSize(0)==1?1:aMatrix.getDimSize(1);
if (aMatrix.isComplex()){
Eigen::Map<Eigen::MatrixXcd> srcM((std::complex<double>*)aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
auto ret = malloc(aMatrix.getDimSize(1),true);
Eigen::Map<Eigen::VectorXcd> retV((std::complex<double>*)ret,aMatrix.getDimSize(1));
Eigen::Map<Eigen::MatrixXcd> srcM((std::complex<double>*)aMatrix.getData(),row,col);
auto ret = malloc(col,true);
Eigen::Map<Eigen::VectorXcd> retV((std::complex<double>*)ret,col);
retV = srcM.colwise().prod();
return Matrix::New(ret,1,aMatrix.getDimSize(1),1,Complex);
return Matrix::New(ret,1,col,1,Complex);
}
else{
Eigen::Map<Eigen::MatrixXd> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
auto ret = malloc(aMatrix.getDimSize(1));
Eigen::Map<Eigen::VectorXd> retV(ret,aMatrix.getDimSize(1));
Eigen::Map<Eigen::MatrixXd> srcM(aMatrix.getData(),row,col);
auto ret = malloc(col);
Eigen::Map<Eigen::VectorXd> retV(ret,col);
retV = srcM.colwise().prod();
return Matrix::New(ret,1,aMatrix.getDimSize(1));
return Matrix::New(ret,1,col);
}
}
Matrix Aurora::dot(const Matrix &aMatrix,const Matrix& aOther,FunctionDirection direction ) {
if ( direction == All){
std::cerr<< "dot() not support 3D data!"
<< std::endl;
std::cerr<< "dot() not support 3D data!"<< std::endl;
return Matrix();
}
if (!aMatrix.compareShape(aOther)){
std::cerr<< "dot() matrix must be same shape!"
<< std::endl;
std::cerr<< "dot() matrix must be same shape!"<< std::endl;
return Matrix();
}
//针对向量行等于列
if (direction == Column && aMatrix.getDimSize(0)==1){
direction = Row;
}
if (aMatrix.isComplex()){
return sum(conj(aMatrix)*aOther,direction);
}