0

我正在寻找使用Armadillo (C++) 的专家。我是一名想成为天文学家的人，我正在尝试在望远镜的特定滤光片中计算 SSP（简单恒星群，一组化学同质和同时代的恒星）的 ab 视星等 m_ab 。

compute_ab执行此计算的函数的输入是SSP 的波长wavelength和相应的光谱能量分布sed（基本上它是 SSP 在一定波长范围内每单位波长的光度）；输入fresp和fwaves是望远镜的吞吐量（基本上是某个波段的过滤器）和吞吐量分布的相应波长范围。它们是std::vector<long double>一维的。我需要输出的是一个数字，可能是双精度或长双精度。

我肯定需要从这些信息中计算 m_ab 是插值，因为 SED 和吞吐量可以在非常不同的波长上；它是卷积和积分。从物理上讲，我在这个函数中所做的段落是正确的，所以我请求一些帮助来使用犰狳，我做得对吗？如何将输出设置为双精度？此外，当我运行我的代码时，我现在遇到了这个错误：

    error: trapz(): length of X must equal the number of rows in Y when dim=0
terminate called after throwing an instance of 'std::logic_error'
  what():  trapz(): length of X must equal the number of rows in Y when dim=0

这是功能：

mat compute_ab(vector <long double> waves,vector <long double> sed, vector <long double> fwaves,vector <long double> fresp){

  colvec filter_interp;

  colvec csed = conv_to< colvec >::from(sed);
  colvec cwaves = conv_to< colvec >::from(waves);
  colvec cfwaves = conv_to< colvec >::from(fwaves);
  colvec cfresp = conv_to< colvec >::from(fresp);

  arma::interp1(cfwaves,cfresp,cwaves,filter_interp);

  colvec filterSpec = arma::conv(filter_interp,csed);

  colvec i1 = arma::conv(filterSpec, cwaves);
  colvec i2 = arma::conv(filter_interp, 1./cwaves); 

  mat I1=arma::trapz(i1,cwaves);
  mat I2=arma::trapz(i2,cwaves);

  mat fnu=I1/I2/speedcunitas;
  mat mAB= -2.5 * log10(fnu) -48.6;
  return mAB;

}

感谢大家的帮助！

4

1 回答 1

1

此代码存在一些问题，并且可能有一些操作没有按照您的意愿进行。

  1. 您正在处理输入向量，这可能非常昂贵且完全没有必要。改变
mat compute_ab(vector <long double> waves,vector <long double> sed, vector <long double> fwaves,vector <long double> fresp){

到

mat compute_ab(const vector <long double>& waves,const vector <long double>& sed, const vector <long double>& fwaves,const vector <long double>& fresp){
  1. colvec并且vec是一回事。双打的列向量。vec如果你愿意，你可以使用 just 。您正在将 转换std::vector<double>为arma::colvecusing conv_to<colvec>::from。这很好，但你又在复制元素。由于您只想对这些std::vector<double>输入执行线性代数运算并且它们仅在函数内部使用，因此您可以做得更好。

例如，考虑下面的代码

std::vector<double> v;
v.push_back(2.4);
v.push_back(-1.4);
v.push_back(10);

arma::vec c = arma::conv_to<arma::colvec>::from(v);
arma::colvec n;

std::cout << &c[0] << std::endl;
std::cout << &v[0] << std::endl;

s将cout打印原始向量和转换后的 vec 中第一个元素的内存地址，这确实是不同的。更好的方法是使用一个特殊的构造函数从arma::vec内存地址初始化向量并告诉它不要复制元素。也就是说，将arma::vec c行更改为

arma::vec c = arma::vec(v.data(), v.size(), false);

的data成员给出std::vector了一个指向其元素的指针，我们可以将它arma::vec与元素的数量一起传递给构造函数。现在couts 将打印相同的地址&c[0]并&v[0]确认元素没有被复制。

  1. 查看arma::conv调用的输出。例如，传递两个每个包含 3 个元素的输入向量将产生一个包含 5 个元素的输出向量。看文档arma::conv。特别是，还有第三个参数称为“形状”，它可能对您有用。
  1. arma::trapz( X, Y )文件说

计算 Y 相对于 X 中间距的梯形积分。X 必须是向量；它的长度必须等于 Y 中的行数。

您传递的输入的尺寸arma::trapz不匹配，您会收到上述错误。原因是 的输出arma::conv可能不是您所期望的，您可能希望将“相同”作为arma::conv.

TLDR；

我的建议是你打印每一行的结果，compute_ab看看它是否是你所期望的。如果不是，请查看该特定犰狳函数的文档。这样做，您将很快对使用犰狳库更有信心。通常犰狳中的一个函数有几个重载。

于 2020-10-12T15:29:10.043 回答