使用ArmadilloC ++计算视在大小

问题描述

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

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

我肯定需要从这些信息中计算出m_ab是一个插值，因为SED和吞吐量可能处于非常不同的波长。这是一个卷积和整合。从物理上讲，我在此功能中所做的段落是正确的，所以我要求使用Armadillo的帮助，我做得正确吗？如何将输出设置为双精度？而且，当我运行我的代码时，我现在收到此错误：

    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;

}

谢谢大家的帮助！

解决方法

此代码存在一些问题，也许还有一些操作没有按照您的要求进行。

---

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。您正在使用std::vector<double>将arma::colvec转换为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;

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

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

std::vector的{​​{3}}给出了指向其元素的指针，我们可以沿着元素数量将其传递给arma::vec构造函数。 现在，cout将为&c[0]和&v[0]打印相同的地址，以确认没有复制元素。

---

1. 看看您的arma::conv呼叫的输出。例如，传递两个分别包含3个元素的输入向量将导致包含5个元素的输出向量。查看data member的文档。特别是，还有第三个参数“ shape”可能对您有用。

---

1. arma::trapz( X,Y )文档说

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

您传递给arma::trapz的输入的尺寸不匹配，您会遇到上述错误。原因是arma::conv的输出可能不是您期望的，并且您可能希望传递与“ arma::conv的第三个参数相同的值。

---

TLDR;

我的建议是，您在compute_ab中打印每一行的结果，以查看其是否符合您的期望。如果不是，请查阅该特定犰狳功能的文档。这样做，您将很快使用armadillo库更加自信。犰狳中的函数通常有多个重载。