test.cpp

#include "eistype.h"
#include "log.h"
#include <cstdio>
#include <iostream>
#include <eisgenerator/model.h>
#include <eisgenerator/normalize.h>
#include <eisgenerator/basicmath.h>
#include <kissinference/kissinference.h>
#include <stddef.h>
#include <filesystem>

#define STD_MODEL "r-cr"

static std::filesystem::path filepath;

static void filterData(std::vector<eis::DataPoint>& data, size_t outputSize)
{
	data = eis::reduceRegion(data);

	if(data.size() < outputSize/8)
	{
		data = std::vector<eis::DataPoint>();
		return;
	}
	data = eis::rescale(data, outputSize/2);
}

static void printCompare(float* re, float* im, const std::vector<eis::DataPoint>& referance)
{
	for(size_t i = 0; i < referance.size(); ++i)
		std::cout<<i<<":\t"<<re[i]<<'+'<<im[i]<<"i\t=\t"<<referance[i].im.real()<<'+'<<referance[i].im.imag()<<"i\n";
}

static bool compare(float* re, float* im, const std::vector<eis::DataPoint>& referance)
{
	for(size_t i = 0; i <referance.size(); ++i)
	{
		if(!eis::fvalueEq(re[i], referance[i].im.real()))
		{
			std::cout<<re[i]<<" at "<<i<<" is not close to "<<referance[i].im.real()<<'\n';
			printCompare(re, im, referance);
			return false;
		}
		if(!eis::fvalueEq(im[i], referance[i].im.imag()))
		{
			std::cout<<im[i]<<" at "<<i<<" is not close to "<<referance[i].im.imag()<<'\n';
			printCompare(re, im, referance);
			return false;
		}
	}
	return true;
}

static void getArrays(const std::vector<eis::DataPoint>& data, float **re, float **im, float **omega)
{
	(*re) = new float[data.size()];
	(*im) = new float[data.size()];
	(*omega) = new float[data.size()];

	for(size_t i = 0; i < data.size(); ++i)
	{
		(*re)[i] = data[i].im.real();
		(*im)[i] = data[i].im.imag();
		(*omega)[i] = data[i].omega;
	}
}

static bool testFilter()
{
	eis::EisSpectra spectra = eis::EisSpectra::loadFromDisk(filepath);
	if(spectra.data.empty())
	{
		std::cout<<"Could not load "<<filepath<<'\n';
		return false;
	}

	float *re, *im, *omega;
	getArrays(spectra.data, &re, &im, &omega);

	std::vector<eis::DataPoint> dataRef = spectra.data;
	filterData(dataRef, 50);

	float *re_filtered, *im_filtered;
	bool ret = kiss_filter_spectra(re, im, omega, spectra.data.size(), &re_filtered, &im_filtered, 25);
	if(!ret)
	{
		std::cout<<"kiss_filter_spectra returned false\n";
		return false;
	}

	return compare(re_filtered, im_filtered, dataRef);
}

static bool testGrad()
{
	bool ret = true;
	eis::Model model(STD_MODEL, 100, false);
	eis::Range omegaRange = eis::Range(1, 1e6, 50, true);
	std::vector<eis::DataPoint> data = model.executeSweep(omegaRange);
	float *re, *im, *omega;
	getArrays(data, &re, &im, &omega);

	std::vector<std::complex<float>> gradsRef(data.size());
	for(size_t i = 0; i < data.size(); ++i)
		gradsRef[i] = eis::absGrad(data, i);

	std::vector<std::complex<float>> grads(data.size());
	for(size_t i = 0; i < data.size(); ++i)
	{
		float *grad = kiss_absgrad(re, im, omega, data.size(), i);
		grads[i] = std::complex<float>(grad[0], grad[1]);
		free(grad);
	}

	for(size_t i = 0; i < data.size(); ++i)
	{
		if(!eis::fvalueEq(gradsRef[i].real(), grads[i].real()) || !eis::fvalueEq(gradsRef[i].imag(), grads[i].imag()))
		{
			std::cerr<<gradsRef[i]<<" != "<<grads[i]<<" at "<<i<<'\n';
			ret = false;
		}
	}

	std::vector<float> realData(data.size());
	for(size_t i = 0; i < data.size(); ++i)
		realData[i] = std::abs(data[i].im);

	std::vector<float> realAbsRef(data.size());
	for(size_t i = 0; i < data.size(); ++i)
		realAbsRef[i] = eis::grad(realData, omegaRange.getRangeVector(), i);

	std::vector<float> realAbs(data.size());
	for(size_t i = 0; i < data.size(); ++i)
		realAbs[i] = kiss_grad(realData.data(), omegaRange.getRangeVector().data(), realAbs.size(), i);

	for(size_t i = 0; i < realAbs.size(); ++i)
	{
		if(!eis::fvalueEq(realAbs[i], realAbsRef[i]))
		{
			std::cerr<<realAbs[i]<<" != "<<realAbsRef[i]<<" at "<<i<<'\n';
			ret = false;
		}
	}

	delete [] re;
	delete [] im;
	delete [] omega;
	return ret;
}

static bool testNorm()
{
	bool ret = true;
	eis::Model model(STD_MODEL, 100, false);
	eis::Range omegaRange = eis::Range(1, 1e6, 50, true);
	std::vector<eis::DataPoint> data = model.executeSweep(omegaRange);
	float *re, *im, *omega;
	getArrays(data, &re, &im, &omega);

	eis::normalize(data);
	kiss_normalize_spectra(re, im, data.size());

	for(size_t i = 0; i < data.size(); ++i)
	{
		if(!eis::fvalueEq(data[i].im.real(), re[i]) || !eis::fvalueEq(data[i].im.imag(), im[i]))
		{
			std::cerr<<i<<":\t"<<re[i]<<'+'<<im[i]<<"i\t!=\t"<<data[i].im.real()<<'+'<<data[i].im.imag()<<"i\n";
			ret = false;
		}
	}

	delete [] re;
	delete [] im;
	delete [] omega;
	return ret;
}

static bool testResample()
{
	bool ret = true;
	eis::EisSpectra spectra = eis::EisSpectra::loadFromDisk(filepath);
	if(spectra.data.empty())
	{
		std::cout<<"Could not load "<<filepath<<'\n';
		return false;
	}
	std::vector<eis::DataPoint> data = spectra.data;
	float *re, *im, *omega;
	size_t size = data.size();
	getArrays(data, &re, &im, &omega);

	data = eis::rescale(data, 10);
	assert(data.size() == 10);

	float *re_resampled, *im_resampled;
	kiss_resample_spectra(re, im, size, &re_resampled, &im_resampled, 10);

	for(size_t i = 0; i < data.size(); ++i)
	{
		if(!eis::fvalueEq(data[i].im.real(), re_resampled[i]) || !eis::fvalueEq(data[i].im.imag(), im_resampled[i]))
		{
			std::cerr<<i<<":\t"<<re_resampled[i]<<'+'<<im_resampled[i]<<"i\t!=\t"<<data[i].im.real()<<'+'<<data[i].im.imag()<<"i\n";
			ret = false;
		}
	}

	free(re_resampled);
	free(im_resampled);

	delete [] re;
	delete [] im;
	delete [] omega;
	return ret;
}

static bool testReduce()
{
	bool ret = true;
	eis::EisSpectra spectra = eis::EisSpectra::loadFromDisk(filepath);
	if(spectra.data.empty())
	{
		std::cout<<"Could not load "<<filepath<<'\n';
		return false;
	}
	std::vector<eis::DataPoint> data = spectra.data;
	float *re, *im, *omega;
	size_t size = data.size();
	getArrays(data, &re, &im, &omega);

	data = eis::reduceRegion(data, 0.01, false);

	float *re_reduced, *im_reduced;
	size_t len;
	kiss_reduce_spectra(re, im, omega, size, 0.01, false, &re_reduced, &im_reduced, &len);

	if(len != data.size())
	{
		std::cerr<<"Wrong output size "<<len<<" vs "<<data.size()<<'\n';
		ret = false;
	}

	if(ret)
	{
		for(size_t i = 0; i < data.size(); ++i)
		{
			if(!eis::fvalueEq(data[i].im.real(), re_reduced[i]) || !eis::fvalueEq(data[i].im.imag(), im_reduced[i]))
			{
				std::cerr<<i<<":\t"<<re_reduced[i]<<'+'<<im_reduced[i]<<"i\t!=\t"<<data[i].im.real()<<'+'<<data[i].im.imag()<<"i\n";
				ret = false;
			}
		}
	}

	free(re_reduced);
	free(im_reduced);

	delete [] re;
	delete [] im;
	delete [] omega;
	return ret;
}

static bool testMedian()
{
	eis::EisSpectra spectra = eis::EisSpectra::loadFromDisk(filepath);
	if(spectra.data.empty())
	{
		std::cout<<"Could not load "<<filepath<<'\n';
		return false;
	}
	std::vector<eis::DataPoint> data = spectra.data;
	std::vector<fvalue> grads;
	grads.reserve(data.size());
	for(size_t i = 0; i < data.size(); ++i)
		grads.push_back(std::abs(eis::absGrad(data, i)));

	float median = kiss_median(grads.data(), grads.size());
	float medianRef = eis::median(grads);

	bool ret = true;
	if(!eis::fvalueEq(median, medianRef))
	{
		ret = false;
		std::cout<<"Grads\n";
		for(float val : grads)
			std::cout<<val<<'\n';
		std::cout<<median<<" is not "<<medianRef<<'\n';
	}

	return ret;
}

bool testTranslate()
{
	const char *model = "r-rc(r-w)";
	char *cdc = kiss_eis_to_cdc(model);
	char *relaxis = kiss_eis_to_relaxis(model);
	std::cout<<"Eis: "<<model<<'\n';
	std::cout<<"CDC: "<<cdc<<'\n';
	std::cout<<"RelaxIS: "<<relaxis<<'\n';

	free(cdc);
	free(relaxis);

	model = "rc-rc-rc";
	cdc = kiss_eis_to_cdc(model);
	relaxis = kiss_eis_to_relaxis(model);
	std::cout<<"Eis: "<<model<<'\n';
	std::cout<<"CDC: "<<cdc<<'\n';
	std::cout<<"RelaxIS: "<<relaxis<<'\n';

	free(cdc);
	free(relaxis);

	return true;
}

int main(int argc, char** argv)
{
	if(argc < 2)
	{
		std::cout<<"A spectra file must be provided\n";
		return 1;
	}

	filepath = argv[1];
	eis::Log::level = eis::Log::ERROR;
	bool same = std::is_same<float, fvalue>::value;
	if(!same)
	{
		std::cerr<<"float and fvalue must be the same type for these tests\n";
		return 1;
	}

	if(!testMedian())
	{
		std::cerr<<"kiss_median not working\n";
		return 2;
	}

	if(!testGrad())
	{
		std::cerr<<"kiss_absgrad or kiss_grad not working\n";
		return 3;
	}

	if(!testNorm())
	{
		std::cerr<<"kiss_normalize_spectra not working\n";
		return 4;
	}

	if(!testResample())
	{
		std::cerr<<"kiss_resample_spectra not working\n";
		return 5;
	}

	if(!testReduce())
	{
		std::cerr<<"kiss_reduce_spectra not working\n";
		return 6;
	}

	if(!testFilter())
	{
		std::cerr<<"kiss_filter_spectra not working\n";
		return 7;
	}

	if(!testTranslate())
	{
		std::cerr<<"kiss_eis_to_cdc or kiss_eis_to_relaxis not working\n";
		return 8;
	}
}