text:
#pragma once

#include <utility>
#include <array>
#include <vector>
#include <unordered_map>
#include <DirectXMathExtend.h>
#include <gsl.h>
#include <VertexTraits.h>
#include <minmax>

namespace Geometrics
{
	template <typename IndexType>
	struct Triangle
	{
		static constexpr size_t VertexCount = 3;
		typedef std::array<IndexType, VertexCount> container_type;
		container_type v;

		typedef IndexType value_type;
		typedef size_t size_type;
		typedef ptrdiff_t difference_type;
		typedef IndexType *pointer;
		typedef const IndexType *const_pointer;
		typedef IndexType& reference;
		typedef const IndexType& const_reference;

		typedef typename container_type::iterator iterator;
		typedef typename container_type::const_iterator const_iterator;

		typedef typename container_type::reverse_iterator reverse_iterator;
		typedef typename container_type::const_reverse_iterator const_reverse_iterator;

		Triangle() {}
		Triangle(const IndexType &v0, const IndexType &v1, const IndexType &v2)
		{
			v[0] = v0; v[1] = v1; v[2] = v2;
		}
		Triangle(const IndexType *v_)
		{
			v[0] = v_[0]; v[1] = v_[1]; v[2] = v_[2];
		}
		template <class S> explicit Triangle(const S &x)
		{
			v[0] = x[0];  v[1] = x[1];  v[2] = x[2];
		}
		IndexType &operator[] (int i) { return v[i]; }
		const IndexType &operator[] (int i) const { return v[i]; }
		operator const IndexType * () const { return &(v[0]); }
		operator IndexType * () { return &(v[0]); }

		void rewind() {
			std::swap(v[0], v[2]);
		}

		int indexof(IndexType v_) const
		{
			return (v[0] == v_) ? 0 :
				(v[1] == v_) ? 1 :
				(v[2] == v_) ? 2 : -1;
		}

		// contrainer access
		auto begin() { return v.begin(); }
		auto end() { return v.end(); }
		auto begin() const { return v.begin(); }
		auto end() const { return v.end(); }

		auto rbegin() { return v.rbegin(); }
		auto rend() { return v.rend(); }
		auto rbegin() const { return v.rbegin(); }
		auto rend() const { return v.rend(); }

		size_t size() const { return VertexCount; }
	};

	// generate per vertex normal for given triangle mesh
	template <typename VertexType, typename IndexType>
	bool generate_normal(gsl::span<VertexType> vertices, gsl::span<const Triangle<IndexType>> facets)
	{
		using namespace DirectX::VertexTraits;
		using namespace DirectX;

		if (!has_normal<VertexType>::value)
			return false;

		//static_assert(has_normal<VertexType>::value, "The vertex type dose not contains normal field");

		std::vector<XMVECTOR, XMAllocator> normals(vertices.size());

		// set zero
		std::memset(normals.data(), 0, normals.size() * sizeof(XMVECTOR));

		XMVECTOR v0, v1, v2, n;
		for (const auto& face : facets)
		{
			v0 = get_position(vertices[face[0]]);
			v1 = get_position(vertices[face[1]]);
			v2 = get_position(vertices[face[2]]);
			v1 -= v0;
			v2 -= v0;

			v1 = XMVector3Normalize(v1);
			v2 = XMVector3Normalize(v2);
			n = XMVector3Cross(v1, v2); // weighted normal

			normals[face[0]] += n;
			normals[face[1]] += n;
			normals[face[2]] += n;
		}

		for (size_t i = 0; i < vertices.size(); i++)
		{
			n = XMVector3Normalize(normals[i]);
			set_normal(vertices[i], n);
		}

		return true;
	}

	// generate per vertex tangent for given triangle mesh
	template <typename VertexType, typename IndexType>
	bool generate_tangent(gsl::span<VertexType> vertices, gsl::span<const Triangle<IndexType>> facets)
	{
		using namespace DirectX::VertexTraits;
		using namespace DirectX;

		if (!has_tangent<VertexType>::value)
			return false;
		//static_assert(has_tangent<VertexType>::value, "The vertex type dose not contains tangent field");

		std::vector<XMVECTOR, XMAllocator> tan1(vertices.size() * 2);
		XMVECTOR* tan2 = &tan1[vertices.size()];
		std::memset(tan1.data(), 0, tan1.size() * sizeof(XMVECTOR));

		for (const auto& face : facets)
		{
			{
				XMVECTOR v0, v1, v2, w0, w1, w2;
				v0 = get_position(vertices[face[0]]);
				w0 = get_uv(vertices[face[0]]);

				v1 = get_position(vertices[face[1]]);
				w1 = get_uv(vertices[face[1]]);

				v2 = get_position(vertices[face[2]]);
				w2 = get_uv(vertices[face[2]]);
			}
			XMFLOAT4A v_1, v_2, w_1, w_2;

			XMStoreA(v_1, v1 - v0);
			XMStoreA(v_2, v2 - v0);
			XMStoreA(w_1, w1 - w0);
			XMStoreA(w_2, w2 - w0);

			float x1 = v_1.x; float x2 = v_2.x;
			float y1 = v_1.y; float y2 = v_2.y;
			float z1 = v_1.z; float z2 = v_2.z;

			float s1 = w_1.x; float s2 = w_2.x;
			float t1 = w_1.y; float t2 = w_2.y;

			float r = 1.0F / (s1 * t2 - s2 * t1);
			XMVECTOR sdir = XMVectorSet(
				(t2 * x1 - t1 * x2) * r,
				(t2 * y1 - t1 * y2) * r,
				(t2 * z1 - t1 * z2) * r, .0f);
			XMVECTOR tdir = XMVectorSet(
				(s1 * x2 - s2 * x1) * r,
				(s1 * y2 - s2 * y1) * r,
				(s1 * z2 - s2 * z1) * r, .0f);

			tan1[face[0]] += sdir;
			tan1[face[1]] += sdir;
			tan1[face[2]] += sdir;

			tan2[face[0]] += tdir;
			tan2[face[1]] += tdir;
			tan2[face[2]] += tdir;
		}

		for (size_t i = 0; i < vertices.size(); i++)
		{
			XMVECTOR n = get_normal(vertices[i], n);
			XMVECTOR t = tan1[i];

			// Gram-Schmidt orthogonalize
			XMVECTOR nt = XMVector3Normalize(t - n * XMVector3Dot(n, t));

			XMVECTOR w = XMVectorLess(XMVector3Dot(XMVector3Cross(n, t), tan2[i]), XMVectorZero());
			w = XMVectorSelect(g_XMNegativeOne.v, g_XMOne.v, w);

			nt = _DXMEXT XMVectorSelect<0, 0, 0, 1>(nt, w);
			set_tangent(vertices[i], nt);
		}

		return true;
	}

	template <typename _VertexType, typename _IndexType = uint16_t, typename _FaceType = Triangle<_IndexType>>
	class PolygonSoup
	{
	public:
		typedef _VertexType VertexType;
		typedef _IndexType	IndexType;
		typedef _FaceType	FaceType;

		static constexpr size_t PolygonVertex = FaceType::VertexCount;

		PolygonSoup() {}

		bool empty() const {
			return vertices.empty() || indices.empty();
		}

		void clear()
		{
			vertices.clear();
			indices.clear();
		}

		inline const FaceType& facet(int idx) const
		{
			return reinterpret_cast<const FaceType&>(indices[idx * FaceType::VertexCount]);
		}

		inline FaceType& facet(int idx)
		{
			return reinterpret_cast<FaceType&>(indices[idx * FaceType::VertexCount]);
		}

		void add_facet(const FaceType& new_facet)
		{
			indices.insert(indices.end(), new_facet.begin(), new_facet.end());
		}

		template <class... _TIndecies>
		std::enable_if_t<sizeof...(_TIndecies) == FaceType::VertexCount> add_facet(_TIndecies... _indices)
		{
			indices.insert(indices.end(), { static_cast<IndexType>(_indices)... });
		}

		gsl::span<FaceType> facets() {
			return gsl::span<FaceType>(
				reinterpret_cast<FaceType*>(indices.data()),
				indices.size() / FaceType::VertexCount);
		}

		gsl::span<const FaceType> facets() const {
			return gsl::span<const FaceType>(
				reinterpret_cast<const FaceType*>(indices.data()),
				indices.size() / FaceType::VertexCount);
		}

		const VertexType& vertex(int facet, int vidx) const
		{
			return this->vertices[facet * FaceType::VertexCount + vidx];
		}

		VertexType& vertex(int facet, int vidx)
		{
			return this->vertices[facet * FaceType::VertexCount + vidx];
		}

		// flip all the polygons' winding and vertices' normal
		void flip()
		{
			using namespace DirectX::VertexTraits;

			if (has_normal<VertexType>::value)
			{
				for (auto& v : vertices)
					set_normal(v, -get_normal(v));
			}

			static constexpr int vc = FaceType::VertexCount;
			for (int i = 0; i < indices.size() / vc; i++)
				std::reverse(indices.begin() + i * vc, indices.begin() + (i + 1)* vc);
		}

		void generate_normal()
		{
			static_assert(PolygonVertex == 3, "This mesh is not trigle mesh, please trianglize it first");
			generate_normal<VertexType,IndexType>(this->vertices, this->facets());
		}

		void generate_tangent()
		{
			static_assert(PolygonVertex == 3, "This mesh is not trigle mesh, please trianglize it first");
			generate_tangent<VertexType, IndexType>(this->vertices, this->facets());
		}

		// applies an uniform transform to all vertices in the mesh
		void XM_CALLCONV transform(DirectX::FXMMATRIX M)
		{
			using namespace DirectX::VertexTraits;
			for (auto& v : vertices)
			{
				XMVECTOR p = get_position(v);
				p = _DXMEXT XMVector3TransformCoord(p, M);
				set_position(v, p);

				if (has_normal<VertexType>::value)
				{
					p = get_normal(v);
					p = _DXMEXT XMVector3TransformNormal(p, M);
					set_normal(v, p);

					if (has_tangent<VertexType>::value)
					{
						p = get_tangent(v);
						p = _DXMEXT XMVector3TransformNormal(p, M);
						set_tangent(v, p);
					}
				}
			}
		}

		std::vector<VertexType> vertices;
		std::vector<IndexType>	indices;
	};

	template <class _VertexType, class _IndexType>
	int XM_CALLCONV intersect(const PolygonSoup<_VertexType, _IndexType, Triangle<_IndexType>> &Mesh, DirectX::FXMVECTOR Origin, DirectX::FXMVECTOR Direction, std::vector<float>* distances = nullptr)
	{
		using namespace DirectX;
		using namespace DirectX::VertexTraits;
		int count = 0;
		XMVECTOR vDir = XMVector3Normalize(Direction);
		for (const auto& tri : Mesh.facets())
		{
			float distance;

			XMVECTOR v0 = get_position(&Mesh.vertices[tri[0]]);
			XMVECTOR v1 = get_position(&Mesh.vertices[tri[1]]);
			XMVECTOR v2 = get_position(&Mesh.vertices[tri[2]]);

			bool hr = DirectX::TriangleTests::Intersects(Origin, vDir, v0, v1, v2, distance);
			if (hr) {
				++count;
				if (distances) {
					distances->push_back(distances);
				}
			}
		}
		return count;
	}

	struct MeshRayIntersectionInfo
	{
		using Vector3 = DirectX::Vector3;

		Vector3 position;
		float	distance;
		Vector3 barycentric;
		int		facet; // facet id

		inline bool operator < (const MeshRayIntersectionInfo& rhs)
		{
			return this->distance < rhs.distance;
		}
	};

	/// <summary>
	/// Basic triangle mesh, each index represent an edge, which is the edge oppsite to the vertex in it's owner triangle
	/// </summary>
	template <typename _VertexType, typename _IndexType = uint16_t>
	class TriangleMesh : public PolygonSoup<_VertexType, _IndexType, Triangle<_IndexType>>
	{
	public:
		static const size_t VertexCount = FaceType::VertexCount;

	public:
		// edge's reverse edge
		// stores the adjacent edges of a edge in a triangle
		std::vector<IndexType> revedges;

		// vertex's first adjacant edge
		// std::vector<IndexType> vedges;

		union EdgeType
		{
			_IndexType v[2];
			struct {
				_IndexType v0, v1;
			};
		};

		inline EdgeType edge(int facet, int edge) const
		{
			EdgeType e;
			auto& tri = this->facet(facet);
			switch (edge)
			{
			case 0:
				e.v0 = tri[1];
				e.v1 = tri[2];
				break;
			case 1:
				e.v0 = tri[2];
				e.v1 = tri[0];
				break;
			case 2:
				e.v0 = tri[0];
				e.v1 = tri[1];
				break;
			default:
				e.v0 = -1;
				e.v1 = -1;
				break;
			}
			return e;
		}

		inline EdgeType edge(int eid) const
		{
			return edge(eid / VertexCount, eid % VertexCount);
		}

		inline int adjacentFacet(int facet, int edge) const
		{
			return revedges[facet * VertexCount + edge] / VertexCount;
		}

		inline int adjacentFacet(int eid) const
		{
			return revedges[eid] / VertexCount;
		}

		// build the adjacent map so we can access all the 1 rings
		void build()
		{
			// intialize all adjacant edges to -1
			revedges.resize(this->indices.size(), -1);
			_IndexType vsize = this->vertices.size();
			int esize = this->indices.size();

			using hash_type = uint64_t;
			// make sure int for hash is enough
			assert(vsize*vsize < std::numeric_limits<hash_type>::max());
			std::unordered_map<hash_type, _IndexType> edges(esize * 2);
			// max items inside this table should be less than (esize / 2)

			_IndexType fid = 0;
			for (const Triangle<_IndexType>& tri : this->facets())
			{
				for (_IndexType i = 0; i < VertexCount; i++)
				{
					_IndexType eid = i + fid * 3;
					auto e = edge(fid, i);

					hash_type ehash = e.v0 * vsize + e.v1;
					hash_type revehash = e.v0 + e.v1 * vsize;

					auto revItr = edges.find(revehash);
					if (revItr == edges.end())
					{
						edges[ehash] = fid + 0;
					}
					else // find reversed edge, remove from edges map
					{
						auto revEid = revItr->second;
						revedges[eid] = revEid;
						revedges[revEid] = eid;
						edges.erase(revItr);
					}
				}

				++fid;
			}
		}

		int XM_CALLCONV intersect(DirectX::FXMVECTOR Origin, DirectX::FXMVECTOR Direction, std::vector<MeshRayIntersectionInfo>* output) const
		{
			//assert(revedges.size() == indices.size());
			using namespace DirectX;
			using namespace DirectX::VertexTraits;
			size_t count = 0;
			XMVECTOR vDir = XMVector3Normalize(Direction);
			auto fid = 0;
			for (const auto& tri : this->facets())
			{
				float distance;

				XMVECTOR v0 = get_position(&this->vertices[tri[0]]);
				XMVECTOR v1 = get_position(&this->vertices[tri[1]]);
				XMVECTOR v2 = get_position(&this->vertices[tri[2]]);

				bool hr = DirectX::TriangleTests::Intersects(Origin, vDir, v0, v1, v2, distance);
				if (hr) {
					++count;
					if (output) {
						output->emplace_back();
						auto& info = output->back();
						XMVECTOR pos = distance * vDir + Origin;
						XMVECTOR bc = DirectX::TriangleTests::BarycentricCoordinate(pos, v0, v1, v2);
						info.facet = fid;
						info.position = pos;
						info.distance = distance;
						info.barycentric = bc;
					}
				}
				++fid;
			}

			if (output)
				std::sort(output->begin(), output->end());
			return count;
		}
	};

	namespace Internal
	{
#ifndef min
		template <typename T>
		inline T clamp(T value, T minV, T maxV)
		{
			return std::max(std::min(value, maxV), minV);
		}
#else
		template <typename T>
		inline T clamp(T value, T minV, T maxV)
		{
			return max(min(value, maxV), minV);
		}
#endif
	}

	/// <summary>
	/// Compute the closest projection point from P0 to triangle(V0,V1,V2).
	/// </summary>
	/// <param name="P0">The p0.</param>
	/// <param name="V0">The v0.</param>
	/// <param name="V1">The v1.</param>
	/// <param name="V2">The v2.</param>
	/// <returns></returns>
	inline DirectX::XMVECTOR XM_CALLCONV Projection(DirectX::FXMVECTOR P0, DirectX::FXMVECTOR V0, DirectX::FXMVECTOR V1, DirectX::GXMVECTOR V2)
	{
		using namespace DirectX;
		using namespace Geometrics::Internal;
		XMVECTOR edge0 = V1 - V0;
		XMVECTOR edge1 = V2 - V0;
		XMVECTOR p0 = V0 - P0;

		//XMVectorBaryCentric();

		float a = XMVectorGetX(_DXMEXT XMVector3LengthSq(edge0));
		float b = XMVectorGetX(_DXMEXT XMVector3Dot(edge0, edge1));
		float c = XMVectorGetX(_DXMEXT XMVector3LengthSq(edge1));
		float d = XMVectorGetX(_DXMEXT XMVector3Dot(edge0, p0));
		float e = XMVectorGetX(_DXMEXT XMVector3Dot(edge1, p0));

		float det = a*c - b*b;
		float s = b*e - c*d;
		float t = b*d - a*e;

		if (s + t < det)
		{
			if (s < 0.f)
			{
				if (t < 0.f)
				{
					if (d < 0.f)
					{
						s = clamp(-d / a, 0.f, 1.f);
						t = 0.f;
					}
					else
					{
						s = 0.f;
						t = clamp(-e / c, 0.f, 1.f);
					}
				}
				else
				{
					s = 0.f;
					t = clamp(-e / c, 0.f, 1.f);
				}
			}
			else if (t < 0.f)
			{
				s = clamp(-d / a, 0.f, 1.f);
				t = 0.f;
			}
			else
			{
				float invDet = 1.f / det;
				s *= invDet;
				t *= invDet;
			}
		}
		else
		{
			if (s < 0.f)
			{
				float tmp0 = b + d;
				float tmp1 = c + e;
				if (tmp1 > tmp0)
				{
					float numer = tmp1 - tmp0;
					float denom = a - 2 * b + c;
					s = clamp(numer / denom, 0.f, 1.f);
					t = 1 - s;
				}
				else
				{
					t = clamp(-e / c, 0.f, 1.f);
					s = 0.f;
				}
			}
			else if (t < 0.f)
			{
				if (a + d > b + e)
				{
					float numer = c + e - b - d;
					float denom = a - 2 * b + c;
					s = clamp(numer / denom, 0.f, 1.f);
					t = 1 - s;
				}
				else
				{
					s = clamp(-e / c, 0.f, 1.f);
					t = 0.f;
				}
			}
			else
			{
				float numer = c + e - b - d;
				float denom = a - 2 * b + c;
				s = clamp(numer / denom, 0.f, 1.f);
				t = 1.f - s;
			}
		}

		return V0 + s * edge0 + t * edge1;
	}

	inline float XM_CALLCONV Distance(DirectX::FXMVECTOR P0, DirectX::FXMVECTOR V0, DirectX::FXMVECTOR V1, DirectX::GXMVECTOR V2)
	{
		using namespace DirectX;
		XMVECTOR vProj = Projection(P0, V0, V1, V2);
		vProj -= P0;
		return XMVectorGetX(_DXMEXT XMVector3Length(vProj));
	}

	template <typename _VertexType, typename _IndexType>
	inline float XM_CALLCONV Distance(const TriangleMesh<_VertexType, _IndexType> &Mesh, DirectX::FXMVECTOR Point)
	{
		using namespace DirectX;
		float minDis = numeric_limits<float>::max();
		for (const auto& tri : Mesh.facets())
		{
			float dis;

			XMVECTOR v0 = XMLoadFloat3(&Mesh.vertices[tri[0]].position);
			XMVECTOR v1 = XMLoadFloat3(&Mesh.vertices[tri[1]].position);
			XMVECTOR v2 = XMLoadFloat3(&Mesh.vertices[tri[2]].position);

			dis = Distance(Point, v0, v1, v2);
			if (dis < minDis) minDis = dis;
			//minDis = std::min(dis,minDis);
		}

		//XMVECTOR vDis = XMVectorReplicate(numeric_limits<float>::max());
		//for (const auto& vertex : Mesh.vertices)
		//{
		//	XMVECTOR vPos = XMLoadFloat3(&vertex.position);
		//	vPos -= Point;
		//	XMVECTOR dis = XMVector3Length(vPos);
		//	vDis = XMVectorMin(vDis,dis);
		//}

		//float minDisV = XMVectorGetX(vDis);
		//assert(minDisV >= minDis);
		return minDis;
		//return XMVectorGetX(minDis);
	}

	template <typename _VertexType, typename _IndexType>
	bool XM_CALLCONV Inside(const TriangleMesh<_VertexType, _IndexType> &Mesh, DirectX::FXMVECTOR Point)
	{
		XMFLOAT3A Direction;
		Direction.x = (float)std::rand() / (RAND_MAX + 1);
		Direction.y = (float)std::rand() / (RAND_MAX + 1);
		Direction.z = (float)std::rand() / (RAND_MAX + 1);
		XMVECTOR vDir = XMLoadFloat3A(&Direction);
		auto count = intersect(Mesh, Point, vDir, nullptr);
		return count & 1; //count % 2
	}

}

max_stars_repo_name: ArcEarth/PPARM

max_issues_repo_name: ArcEarth/PPARM

max_forks_repo_name: ArcEarth/PPARM


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>

//THESE ARE GSL FUNCTIONS
//YOU DO NOT NEED TO INCLUDE ALL THESE HEADER FILES IN YOUR CODE
//JUST THE ONES YOU ACTUALLY NEED;
//IN THIS APPLICATION, WE ONLY NEED gsl/gsl_matrix.h
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_sort_double.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_errno.h>

void printmatrix(const char* filename,gsl_matrix* m);
gsl_matrix* transposematrix(gsl_matrix* m);
void matrixproduct(gsl_matrix* m1,gsl_matrix* m2,gsl_matrix* m);
void getInverse(gsl_matrix* K, gsl_matrix* inverse);
gsl_matrix* MakeSubmatrix(gsl_matrix* M,
			  int* IndRow,int lenIndRow,
			  int* IndColumn,int lenIndColumn);
double logdet(gsl_matrix* K);
gsl_matrix* diagMatrix(int n);
gsl_matrix* getcholesky(gsl_matrix* m);
gsl_matrix* getsimusample(int size, int p, gsl_rng* r, gsl_matrix* mat_chol);
gsl_vector*  getInverseLogit(gsl_vector*  x);
gsl_vector* getInverseLogit2(gsl_vector* x);
gsl_vector* getPi(gsl_vector* x, gsl_vector* beta);
gsl_vector* getPi2(gsl_vector* x, gsl_vector* beta);
double getLoglik(gsl_vector* y, gsl_vector* x, gsl_vector* beta);
double getLstar(gsl_vector* y, gsl_vector* x, gsl_vector* beta);
void getGradient(gsl_vector* y, gsl_vector* x, gsl_vector* beta, gsl_matrix* Gradient);
void getHessian(gsl_vector* y, gsl_vector* x, gsl_vector* beta, gsl_matrix* Hessian);
double getMaxBetaDiff(gsl_vector* betaCurrent,  gsl_vector* betaNew);
gsl_vector* getcoefNR(int response, int explanatory, gsl_matrix* data, double error);
gsl_vector* getMHinteration(gsl_rng* r, gsl_vector* y, gsl_vector* x, gsl_vector* betaCurrent, gsl_matrix* NegInvHessian);
double getLaplace(int response, int explanatory, gsl_matrix* data, gsl_vector* betaMode);
double getMonteCarlo(int NumIterations, gsl_rng* r, int response, int explanatory, gsl_matrix* data);
gsl_vector* getcoefMH(int NumIterations, gsl_rng* r, int response, int explanatory, gsl_matrix* data, gsl_vector* betaMode);

ext: h


text:
#pragma once
/*
 * (C) Copyright 2020-2021 UCAR
 *
 * This software is licensed under the terms of the Apache Licence Version 2.0
 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
 */
/*! \addtogroup ioda_cxx_variable
 *
 * @{
 * \file Has_Variables.h
 * \brief Interfaces for ioda::Has_Variables and related classes.
 */

#include <cstring>
#include <gsl/gsl-lite.hpp>
#include <iostream>
#include <map>
#include <memory>
#include <string>
#include <utility>
#include <vector>

#include "ioda/Attributes/Attribute_Creator.h"
#include "ioda/Exception.h"
#include "ioda/Layout.h"
#include "ioda/Misc/Eigen_Compat.h"
#include "ioda/Misc/MergeMethods.h"
#include "ioda/Types/Type.h"
#include "ioda/Variables/FillPolicy.h"
#include "ioda/Variables/Variable.h"
#include "ioda/defs.h"

namespace ioda {
class Has_Variables;
class ObsGroup;
struct Named_Variable;
namespace detail {
class Has_Variables_Backend;
class Has_Variables_Base;
class DataLayoutPolicy;
class Group_Base;

}  // namespace detail

/// \brief A few chunking strategies for Variables
namespace chunking {
/// Convenience function for setting default chunking parameters.
inline bool Chunking_Max(const std::vector<Dimensions_t>& in, std::vector<Dimensions_t>& out) {
  out = in;
  return true;
}
}  // namespace chunking

/// \brief Used to specify Variable creation-time properties.
/// \ingroup ioda_cxx_variable
struct IODA_DL VariableCreationParameters {
private:
  std::vector<std::pair<unsigned int, Variable> > dimsToAttach_;
  std::string dimScaleName_;

public:
  /// @name Fill Value
  /// @{

  detail::FillValueData_t fillValue_;

  template <class DataType>
  VariableCreationParameters& setFillValue(DataType fill) {
    detail::assignFillValue<DataType>(fillValue_, fill);
    return *this;
  }
  inline VariableCreationParameters& unsetFillValue() {
    fillValue_.set_ = false;
    return *this;
  }

  /// @}
  /// @name Chunking and compression
  /// @{

  /// \brief Do we chunk this variable? Required for extendible / compressible Variables.
  /// \details Requires a chunking strategy.
  bool chunk = false;

  /// \brief Manually specify the chunks. Never directly use. Use getChunks(...) instead.
  std::vector<Dimensions_t> chunks;
  /// Set variable chunking strategy. Used only if chunk == true and chunks.size() == 0.
  std::function<bool(const std::vector<Dimensions_t>&, std::vector<Dimensions_t>&)>
    fChunkingStrategy = chunking::Chunking_Max;
  /// Figure out the chunking size
  /// \param cur_dims are the current dimensions
  std::vector<Dimensions_t> getChunks(const std::vector<Dimensions_t>& cur_dims) const {
    if (chunks.size()) return chunks;
    std::vector<Dimensions_t> res;
    if (fChunkingStrategy(cur_dims, res)) return res;
    throw Exception("Cannot figure out an appropriate chunking size.", ioda_Here());
  }

  bool gzip_                        = false;
  bool szip_                        = false;
  int gzip_level_                   = 6;  // 1 (fastest) - 9 (most compression)
  unsigned int szip_PixelsPerBlock_ = 16;
  unsigned int szip_options_        = 4;  // Defined as H5_SZIP_EC_OPTION_MASK in hdf5.h;

  void noCompress();
  void compressWithGZIP(int level = 6);
  void compressWithSZIP(unsigned PixelsPerBlock = 16, unsigned options = 4);

  /// @}
  /// @name General Functions
  /// @{

  /// Set any initial attributes here
  Attribute_Creator_Store atts;

  VariableCreationParameters();
  VariableCreationParameters(const VariableCreationParameters&);
  VariableCreationParameters& operator=(const VariableCreationParameters&);

  template<class DataType>
  static VariableCreationParameters defaulted() {
    VariableCreationParameters ret;
    ret.chunk = true;
    ret.compressWithGZIP();
    FillValuePolicies::applyFillValuePolicy<DataType>(FillValuePolicy::NETCDF4, ret.fillValue_);
    return ret;
  }
  template <class DataType>
  static VariableCreationParameters defaults() {
    return defaulted<DataType>();
  }

  /// Finalize routine to make sure struct members are intact (e.g. for fill values)
  detail::FillValueData_t::FillValueUnion_t finalize() const { return fillValue_.finalize(); }

  detail::python_bindings::VariableCreationFillValues<VariableCreationParameters> _py_setFillValue;

private:
  friend class detail::Has_Variables_Base;
  /// Apply the properties to a Variable (second pass; after Variable is created).
  Variable applyImmediatelyAfterVariableCreation(Variable h) const;

  /// @}
};

typedef std::vector<Variable> NewVariables_Scales_t;
/// \brief Used to specify a new variable with the collective createWithScales function.
struct IODA_DL NewVariable_Base : std::enable_shared_from_this<NewVariable_Base> {
  /// Name of the variable.
  std::string name_;
  /// Type of the new dimension. Int, char, etc. Used if a type is not passed directly.
  std::type_index dataType_;
  /// Type of the new dimension. Used if a type is passed directly.
  Type dataTypeKnown_;
  /// Dimension scales
  NewVariables_Scales_t scales_;
  /// Var creation params
  VariableCreationParameters vcp_;

  virtual ~NewVariable_Base() {}

  NewVariable_Base(const std::string& name, const Type& dataType,
                   const NewVariables_Scales_t& scales,
                   const VariableCreationParameters& params)
      : name_(name), dataType_(typeid(void)), dataTypeKnown_(dataType),
        scales_(scales),
        vcp_(params) {}

  NewVariable_Base(const std::string& name, const std::type_index& dataType,
                   const NewVariables_Scales_t& scales,
                   const VariableCreationParameters& params)
      : name_(name),
        dataType_(dataType),
        scales_(scales),
        vcp_(params) {}
};
typedef std::vector<std::shared_ptr<NewVariable_Base>> NewVariables_t;

template <class DataType>
inline std::shared_ptr<NewVariable_Base> NewVariable(
  const std::string& name, const NewVariables_Scales_t& scales,
  const VariableCreationParameters& params = VariableCreationParameters::defaulted<DataType>()) {
  return std::make_shared<NewVariable_Base>(name, typeid(DataType), scales, params);
}

inline std::shared_ptr<NewVariable_Base> NewVariable(const std::string& name,
                                                     const Type& DataType,
                                                     const NewVariables_Scales_t& scales,
                                                     const VariableCreationParameters& params
                                                     = VariableCreationParameters()) {
  return std::make_shared<NewVariable_Base>(name, DataType, scales, params);
}

namespace detail {

/// \ingroup ioda_cxx_variable
class IODA_DL Has_Variables_Base {
  // friend class Group_Base;
private:
  /// Using an opaque object to implement the backend.
  std::shared_ptr<Has_Variables_Backend> backend_;
  /// Set by ObsGroup.
  std::shared_ptr<const detail::DataLayoutPolicy> layout_;
  std::vector<ComplementaryVariableCreationParameters> complementaryVariables_;
  /// \brief FillValuePolicy helper
  /// \details Hides the template function calls, so that the headers are smaller.
  static void _py_fvp_helper(BasicTypes dataType, FillValuePolicy& fvp,
                             VariableCreationParameters& params);

  ComplementaryVariableCreationParameters createDerivedVariableParameters(
      const std::string &inputName, const std::string &outputName, size_t position);
  std::vector<std::vector<std::string>> loadComponentVariableData(
      const ComplementaryVariableCreationParameters& derivedVariableParams);

protected:
  Has_Variables_Base(std::shared_ptr<Has_Variables_Backend>,
                     std::shared_ptr<const DataLayoutPolicy> = nullptr);

public:
  virtual ~Has_Variables_Base();

  /// Set the mapping policy to determine the Layout of Variables stored under this Group.
  /// Usually only set by ObsGroup when we create / open.
  virtual void setLayout(std::shared_ptr<const detail::DataLayoutPolicy>);

  /// Query the backend and get the type provider.
  virtual Type_Provider* getTypeProvider() const;

  /// \brief Get the fill value policy used for Variables within this Group
  /// \details The backend has to be consulted for this operation. Storage of this policy is
  /// backend-dependent.
  virtual FillValuePolicy getFillValuePolicy() const;

  /// @name General Functions
  /// @{
  ///

  /// \brief Does a Variable with the specified name exist?
  /// \param name is the name of the Variable that we are looking for.
  /// \returns true if it exists.
  /// \returns false otherwise.
  virtual bool exists(const std::string& name) const;
  /// \brief Delete an Attribute with the specified name.
  /// \param attname is the name of the Variable that we are deleting.
  /// \throws ioda::Exception if no such attribute exists.
  virtual void remove(const std::string& name);
  /// \brief Open a Variable by name
  /// \param name is the name of the Variable to be opened.
  /// \returns An instance of a Variable that can be queried (with getDimensions()) and read.
  virtual Variable open(const std::string& name) const;
  /// \brief Open a Variable by name
  /// \param name is the name of the Variable to be opened.
  /// \returns An instance of a Variable that can be queried (with getDimensions()) and read.
  inline Variable operator[](const std::string& name) const { return open(name); }

  /// List all Variables under this group (one-level search).
  /// \see Group_Base::listObjects if you need to enumerate both Groups and Variables, or
  ///   if you need recursion.
  virtual std::vector<std::string> list() const;
  /// Convenience function to list all Variables under this group (one-level search).
  /// \see Group_Base::listObjects if you need to enumerate both Groups and Variables, or
  ///   if you need recursion.
  inline std::vector<std::string> operator()() const { return list(); }

  /// \brief Combines all complementary variables as specified in the mapping file, opens them,
  /// and optionally removes the originals from the ObsGroup.
  ///
  /// \p removeOriginals determines if the original complementary variables should be removed from
  /// the ObsGroup. Later functionality will ensure that the original complementary variables can
  /// be recreated on writing back to the original file.
  void stitchComplementaryVariables(bool removeOriginals = true);

  /// \brief Converts unit to SI for all eligible variables. If conversion function not defined,
  /// stores unit as attribute.
  ///
  /// Makes the conversion if the variable's unit is defined in the mapping file and the unit conversion
  /// is defined in UnitConversions.h.
  void convertVariableUnits(std::ostream &out = std::cerr);

  /// \brief Create a Variable without setting its data.
  /// \param attrname is the name of the Variable.
  /// \param dimensions is a vector representing the size of the metadata.
  ///   Each element of the vector is a dimension with a certain size.
  /// \param in_memory_datatype is the runtime description of the Attribute's data type.
  /// \returns A Variable that can be written to.
  virtual Variable create(const std::string& name, const Type& in_memory_dataType,
                          const std::vector<Dimensions_t>& dimensions     = {1},
                          const std::vector<Dimensions_t>& max_dimensions = {},
                          const VariableCreationParameters& params = VariableCreationParameters());

  /// Python compatability function
  /// \note Multiple ways to specify dimensions to match possible
  ///   Python function signatures.
  Variable _create_py(const std::string& name, BasicTypes dataType,
                             const std::vector<Dimensions_t>& cur_dimensions = {1},
                             const std::vector<Dimensions_t>& max_dimensions = {},
                             const std::vector<Variable>& dimension_scales   = {},
                             const VariableCreationParameters& params
                             = VariableCreationParameters());

  inline Variable create(const std::string& name, const Type& in_memory_dataType,
                         const ioda::Dimensions& dims,
                         const VariableCreationParameters& params = VariableCreationParameters()) {
    return create(name, in_memory_dataType, dims.dimsCur, dims.dimsMax, params);
  }

  /// \brief Create a Variable without setting its data.
  /// \tparam DataType is the type of the data. I.e. float, int32_t, uint16_t, std::string, etc.
  /// \param name is the name of the Variable.
  /// \param dimensions is a vector representing the size of the metadata. Each element of the
  ///   vector is a dimension with a certain size.
  /// \returns A Variable that can be written to.
  template <class DataType>
  Variable create(const std::string& name, const std::vector<Dimensions_t>& dimensions = {1},
                  const std::vector<Dimensions_t>& max_dimensions = {},
                  const VariableCreationParameters& params        = VariableCreationParameters::defaulted<DataType>()) {
    try {
      VariableCreationParameters params2 = params;
      FillValuePolicies::applyFillValuePolicy<DataType>(getFillValuePolicy(), params2.fillValue_);
      Type in_memory_dataType = Types::GetType<DataType>(getTypeProvider());
      auto var                = create(name, in_memory_dataType, dimensions,
        max_dimensions, params2);
      return var;
    } catch (...) {
      std::throw_with_nested(Exception(ioda_Here()));
    }
  }

  template <class DataType>
  Variable create(const std::string& name, const ioda::Dimensions& dims,
                  const VariableCreationParameters& params
                  = VariableCreationParameters::defaulted<DataType>()) {
    try {
      VariableCreationParameters params2 = params;
      FillValuePolicies::applyFillValuePolicy<DataType>(getFillValuePolicy(), params2.fillValue_);
      return create<DataType>(name, dims.dimsCur, dims.dimsMax, params2);
    } catch (...) {
      std::throw_with_nested(Exception(ioda_Here()));
    }
  }

  /// \brief Convenience function to create a Variable from certain dimension scales.
  /// \tparam DataType is the type of the data. I.e. int, int32_t, uint16_t, std::string, etc.
  /// \param name is the name of the Variable.
  /// \param dimensions is a vector representing the size of the metadata. Each element of the
  ///   vector is a dimension with a certain size.
  /// \returns A Variable that can be written to.
  template <class DataType>
  Variable createWithScales(const std::string& name,
                            const std::vector<Variable>& dimension_scales,
                            const VariableCreationParameters& params
                            = VariableCreationParameters::defaulted<DataType>()) {
    try {
      Type in_memory_dataType = Types::GetType<DataType>(getTypeProvider());

      NewVariables_t newvars{NewVariable(name, in_memory_dataType, dimension_scales, params)};
      createWithScales(newvars);
      return open(name);
    } catch (...) {
      std::throw_with_nested(Exception(ioda_Here()));
    }
  }

  /// @brief Collective function optimized to mass-construct variables and attach scales.
  /// @param newvars is a vector of the new variables to be created.
  /// @see NewVariable for the signature of the objects to add.
  void createWithScales(const NewVariables_t& newvars);

  /// @}
  /// @name Collective functions
  /// @brief These functions apply the an operation to a *set* of variables in situations where
  ///   such an operation would produce better performance results than a loop of serial
  ///   function calls.
  /// @{

  /// @brief Attach dimension scales to many Dimension Numbers in a set of Variables.
  /// @param DimensionNumber
  /// @param mapping is the scale mappings for each variable. The first part of the pair refers
  ///   to the variable that you are attaching scales to. The second part is a sequence of
  ///   scales that are attached along each dimension (indexed by the vector).
  /// @details
  /// For some backends, particularly HDF5, attaching a dimension scale to a variable is a slow
  /// procedure when you have many variables. This function batches low-level calls and avoids
  /// loops.
  virtual void attachDimensionScales(
    const std::vector<std::pair<Variable, std::vector<Variable>>>& mapping);

  /// @}
};

class IODA_DL Has_Variables_Backend : public detail::Has_Variables_Base {
protected:
  Has_Variables_Backend();

public:
  virtual ~Has_Variables_Backend();
  FillValuePolicy getFillValuePolicy() const override;
  void attachDimensionScales(
    const std::vector<std::pair<Variable, std::vector<Variable>>>& mapping) override;
};
}  // namespace detail

/// \brief This class exists inside of ioda::Group and provides the interface to manipulating
///   Variables.
/// \ingroup ioda_cxx_variable
///
/// \note It should only be constructed inside of a Group. It has no meaning elsewhere.
/// \see ioda::Variable for the class that represents individual variables.
/// \throws ioda::Exception on all exceptions.
class IODA_DL Has_Variables : public detail::Has_Variables_Base {
public:
  virtual ~Has_Variables();
  Has_Variables();
  Has_Variables(std::shared_ptr<detail::Has_Variables_Backend>,
                std::shared_ptr<const detail::DataLayoutPolicy> = nullptr);
};
}  // namespace ioda

/// @}

max_stars_repo_name: NOAA-EMC/ioda

max_issues_repo_name: NOAA-EMC/ioda

max_forks_repo_name: NOAA-EMC/ioda


text:
/*
 * BRAINS
 * (B)LR (R)everberation-mapping (A)nalysis (I)n AGNs with (N)ested (S)ampling
 * Yan-Rong Li, liyanrong@ihep.ac.cn
 * Thu, Aug 4, 2016
 */

/*! \file dnest_sa1d.c
 *  \brief run dnest sampling for sa and 2d rm analysis.
 */

#ifdef SpecAstro

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <stddef.h>
#include <math.h>
#include <float.h>
#include <string.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_interp.h>
#include <mpi.h>

#include "brains.h"

DNestFptrSet *fptrset_sa2d;

int dnest_sa2d(int argc, char **argv)
{
  int i;
  double logz_sa2d;

  if(parset.flag_sa_par_mutual == 0) /* SA and RM have the same BLR */
  {
    set_blr_model2d();
    num_params_sa_blr_model = 0;
  }
  else if(parset.flag_sa_par_mutual == 1) /* SA and RM have different BLRs but share the same mbh and inc. */
  {
    set_blr_model2d();
    set_sa_blr_model();
  }
  else
  {
    set_blr_model2d();
    set_sa_blr_model();
  }

  /* parameter order: RM BLR, SA BLR, RM continuum */
  /* RM */
  num_params_blr = num_params_blr_model + num_params_nlr
                 + num_params_res + num_params_linecenter + 1; /* include and line sys err */
  num_params_rm = parset.n_con_recon + num_params_blr + num_params_var;

  /* SA */
  num_params_sa_blr = num_params_sa_blr_model + num_params_sa_extpar;
  num_params_sa = num_params_sa_blr;

  /* total */
  num_params_blr_tot = num_params_blr + num_params_sa_blr;
  num_params = num_params_sa + num_params_rm;
  idx_resp = num_params_blr_tot + num_params_drw + num_params_trend;
  idx_difftrend = idx_resp + num_params_resp;
  idx_linecenter = num_params_blr_model + num_params_nlr + num_params_res;

  rnd_frac = fmax(0.2, 1.0*(num_params_blr_tot+num_params_var)/num_params);

  par_fix = (int *) malloc(num_params * sizeof(int));
  par_fix_val = (double *) malloc(num_params * sizeof(double));
  par_range_model = malloc( num_params * sizeof(double *));
  par_prior_gaussian = malloc( num_params * sizeof(double *));
  for(i=0; i<num_params; i++)
  {
    par_range_model[i] = malloc(2*sizeof(double));
    par_prior_gaussian[i] = malloc(2*sizeof(double));
  }
  par_prior_model = malloc( num_params * sizeof(int));

  fptrset_sa2d = dnest_malloc_fptrset();
  /* setup functions used for dnest*/
  fptrset_sa2d->from_prior = from_prior_sa2d;
  fptrset_sa2d->print_particle = print_particle_sa2d;
  fptrset_sa2d->restart_action = restart_action_sa;
  fptrset_sa2d->accept_action = accept_action_sa2d;
  fptrset_sa2d->kill_action = kill_action_sa2d;
  fptrset_sa2d->perturb = perturb_sa2d;
  fptrset_sa2d->read_particle = read_particle_sa2d;

  fptrset_sa2d->log_likelihoods_cal_initial = log_likelihoods_cal_initial_sa2d;
  fptrset_sa2d->log_likelihoods_cal_restart = log_likelihoods_cal_restart_sa2d;
  fptrset_sa2d->log_likelihoods_cal = log_likelihoods_cal_sa2d;

  set_par_range_sa2d();
  set_par_fix_blrmodel();
  set_par_fix_sa_blrmodel();

  /* the rest parameters */
  for(i=num_params_blr_model; i<num_params_blr; i++)
  {
    par_fix[i] = 0;
    par_fix_val[i] = -DBL_MAX;
  }
  for(i=num_params_blr + num_params_sa_blr_model; i<num_params; i++)
  {
    par_fix[i] = 0;
    par_fix_val[i] = -DBL_MAX;
  }

  /* fix non-linear response */
  if(parset.flag_nonlinear !=1)
  {
    par_fix[idx_resp + 1] = 1;
    par_fix_val[idx_resp + 1] = 0.0;
  }

  /* fix systematic error of line */
  if(parset.flag_line_sys_err != 1)
  {
    par_fix[num_params_blr-1] = 1;
    par_fix_val[num_params_blr-1] = log(1.0);
  }

  /* fix mbh and inc of SA BLR if SA and RM have different BLR */
  if(parset.flag_sa_par_mutual == 1)
  {
    /* get the index for mbh and inclination parameters */
    set_idx_par_mutual();
    /* mbh */
    par_fix[num_params_blr + idx_sa_par_mutual[0]] = 1;
    par_fix_val[num_params_blr + idx_sa_par_mutual[0]] = 0.0;

    /* inc */
    par_fix[num_params_blr + idx_sa_par_mutual[1]] = 1;
    par_fix_val[num_params_blr + idx_sa_par_mutual[1]] = 0.0;
  }
  else if(parset.flag_sa_par_mutual == 2)
  {
    /* get the index for mbh and inclination parameters */
    set_idx_par_mutual();
    /* inc and opening angle */
    par_fix[num_params_blr + idx_sa_par_mutual[1]] = 1;
    par_fix_val[num_params_blr + idx_sa_par_mutual[1]] = 0.0;
    par_fix[num_params_blr + idx_sa_par_mutual[2]] = 1;
    par_fix_val[num_params_blr + idx_sa_par_mutual[2]] = 0.0;
    /* dynamical parameters from mass on */
    for(i=idx_sa_par_mutual[0]; i<num_params_sa_blr_model; i++)
    {
      par_fix[num_params_blr + i] = 1;
      par_fix_val[num_params_blr + i] = 0.0;
    }
  }

  /* fix FA */
  par_fix[num_params_blr + num_params_sa_blr_model+2] = 1;
  par_fix_val[num_params_blr + num_params_sa_blr_model+2] = log(1.0);

  /* fix systematic error of continuum */
  if(parset.flag_con_sys_err != 1)
  {
    par_fix[num_params_blr_tot] = 1;
    par_fix_val[num_params_blr_tot] = log(1.0);
  }

  /* fix continuum variation parameter sigma and tau if flag_fixvar is true */
  if(parset.flag_fixvar == 1)
  {
    par_fix[num_params_blr_tot + 1] = 1;
    par_fix_val[num_params_blr_tot + 1] = var_param[1];
    par_fix[num_params_blr_tot + 2] = 1;
    par_fix_val[num_params_blr_tot + 2] = var_param[2];
  }

  print_par_names_sa2d();

  force_update = parset.flag_force_update;
  if(parset.flag_para_name != 1)
    logz_sa2d = dnest(argc, argv, fptrset_sa2d, num_params, NULL, NULL, NULL, "data/", dnest_options_file, NULL, NULL);

  dnest_free_fptrset(fptrset_sa2d);

  return 0;
}

void set_par_range_sa2d()
{
  int i, j, i1, i2;

  /* setup parameter range, BLR parameters first */
  for(i=0; i<num_params_blr_model; i++)
  {
    par_range_model[i][0] = blr_range_model[i][0];
    par_range_model[i][1] = blr_range_model[i][1];

    par_prior_model[i] = UNIFORM;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 0.0;
  }

  /* cope with narrow line */
  for(i=num_params_blr_model; i<num_params_blr_model + num_params_nlr; i++)
  {
    par_range_model[i][0] = nlr_range_model[i - num_params_blr_model][0];
    par_range_model[i][1] = nlr_range_model[i - num_params_blr_model][1];

    par_prior_model[i] = nlr_prior_model[i - num_params_blr_model];
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 1.0; /* note that for logarithm prior of flux, this value is not used, so does not matter */
  }
  /* cope with spectral broadening */
  for(i=num_params_blr_model + num_params_nlr; i<num_params_blr_model + num_params_nlr + num_params_res; i++)
  {
    par_range_model[i][0] = -10.0;
    par_range_model[i][1] =  10.0;

    par_prior_model[i] = GAUSSIAN;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 1.0;
  }
  /* cope with line center */
  for(i=num_params_blr-num_params_linecenter-1; i< num_params_blr-1; i++)
  {
    par_range_model[i][0] = -10.0;
    par_range_model[i][1] =  10.0;

    par_prior_model[i] = GAUSSIAN;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 1.0;
  }

  /* the last is systematic error */
  i = num_params_blr-1;
  par_range_model[i][0] = sys_err_line_range[0];
  par_range_model[i][1] = sys_err_line_range[1];

  par_prior_model[i] = UNIFORM;
  par_prior_gaussian[i][0] = 0.0;
  par_prior_gaussian[i][1] = 0.0;

  /* SA BLR parameters */
  for(i=num_params_blr; i<num_params_blr + num_params_sa_blr_model; i++)
  {
    par_range_model[i][0] = sa_blr_range_model[i-num_params_blr][0];
    par_range_model[i][1] = sa_blr_range_model[i-num_params_blr][1];

    par_prior_model[i] = UNIFORM;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 0.0;
  }
  /* SA extra parameters */
  for(i=num_params_blr + num_params_sa_blr_model; i<num_params_blr_tot; i++)
  {
    par_range_model[i][0] = sa_extpar_range[i-(num_params_blr + num_params_sa_blr_model)][0];
    par_range_model[i][1] = sa_extpar_range[i-(num_params_blr + num_params_sa_blr_model)][1];

    par_prior_model[i] = UNIFORM;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 0.0;
  }

  /* variability parameters */
   /* first systematic error */
  i = num_params_blr_tot;
  par_range_model[i][0] = var_range_model[0][0];
  par_range_model[i][1] = var_range_model[0][1];

  par_prior_model[i] = UNIFORM;
  par_prior_gaussian[i][0] = 0.0;
  par_prior_gaussian[i][1] = 0.0;

  for(i=num_params_blr_tot+1; i<num_params_drw + num_params_blr_tot; i++)
  {
    if(var_param_std[i-num_params_blr_tot] > 0.0)
    {
      par_range_model[i][0] = var_param[i-num_params_blr_tot] - 5.0 * var_param_std[i-num_params_blr_tot];
      par_range_model[i][1] = var_param[i-num_params_blr_tot] + 5.0 * var_param_std[i-num_params_blr_tot];

      /* make sure that the range lies within the initial range */
      par_range_model[i][0] = fmax(par_range_model[i][0], var_range_model[i-num_params_blr_tot][0]);
      par_range_model[i][1] = fmin(par_range_model[i][1], var_range_model[i-num_params_blr_tot][1]);

      par_prior_model[i] = GAUSSIAN;
      par_prior_gaussian[i][0] = var_param[i-num_params_blr_tot];
      par_prior_gaussian[i][1] = var_param_std[i-num_params_blr_tot];
    }
    else
    {
      par_range_model[i][0] = var_range_model[i-num_params_blr_tot][0];
      par_range_model[i][1] = var_range_model[i-num_params_blr_tot][1];

      par_prior_model[i] = UNIFORM;
      par_prior_gaussian[i][0] = 0.0;
      par_prior_gaussian[i][1] = 0.0;
    }
  }
  /* long-term trend of continuum */
  for(i=num_params_drw + num_params_blr_tot; i< num_params_drw + num_params_trend + num_params_blr_tot; i++)
  {
    par_range_model[i][0] = var_range_model[3][0];
    par_range_model[i][1] = var_range_model[3][1];

    par_prior_model[i] = GAUSSIAN;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 1.0;
  }

  /* response A and Ag */
  j = 0;
  i1 = idx_resp;
  i2 = idx_resp + num_params_resp;
  for(i=i1; i<i2; i++)
  {
    par_range_model[i][0] = resp_range[j][0];
    par_range_model[i][1] = resp_range[j][1];

    par_prior_model[i] = UNIFORM;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 0.0;

    j++;
  }

  /* different trend in continuum and line */
  j = 0;
  i1 = idx_difftrend;
  i2 = idx_difftrend + num_params_difftrend;
  for(i= i1; i< i2; i++)
  {
    par_range_model[i][0] = var_range_model[4 + j][0];
    par_range_model[i][1] = var_range_model[4 + j][1];

    par_prior_model[i] = UNIFORM;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 0.0;

    j++;
  }

  /* continuum ligth curve values */
  for(i=num_params_blr_tot+num_params_var; i<num_params; i++)
  {
    par_range_model[i][0] = var_range_model[4+num_params_difftrend][0];
    par_range_model[i][1] = var_range_model[4+num_params_difftrend][1];

    par_prior_model[i] = GAUSSIAN;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 1.0;
  }

  return;
}

/*!
 *  print names and prior ranges for parameters
 *
 */
void print_par_names_sa2d()
{
  if(thistask != roottask)
    return;

  int i, j;
  FILE *fp;
  char fname[BRAINS_MAX_STR_LENGTH], str_fmt[BRAINS_MAX_STR_LENGTH];

  sprintf(fname, "%s/%s", parset.file_dir, "data/para_names_sa2d.txt");
  fp = fopen(fname, "w");
  if(fp == NULL)
  {
    fprintf(stderr, "# Error: Cannot open file %s.\n", fname);
    exit(0);
  }

  strcpy(str_fmt, "%4d %-15s %10.6f %10.6f %4d %4d %15.6e\n");

  printf("# Print parameter name in %s\n", fname);

  fprintf(fp, "#*************************************************\n");
  fprint_version(fp);
  fprintf(fp, "#*************************************************\n");

  fprintf(fp, "%4s %-15s %10s %10s %4s %4s %15s\n", "#", "Par", "Min", "Max", "Prior", "Fix", "Val");

  i=-1;
  for(j=0; j<num_params_blr_model; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "BLR model", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  }

  for(j=0; j<num_params_nlr; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "narrow line", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  }
  for(j=0; j<num_params_res; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "line broaden", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  }
  for(j=0; j<num_params_linecenter; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "line center", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  }

  i++;
  fprintf(fp, str_fmt, i, "sys_err_line", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);

  for(j=0; j<num_params_sa_blr_model; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "SA BLR model", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  }

  for(j=0; j<num_params_sa_extpar; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "SA Extra Par", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  }

  i++;
  fprintf(fp, str_fmt, i, "sys_err_con", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  i++;
  fprintf(fp, str_fmt, i, "sigmad", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  i++;
  fprintf(fp, str_fmt, i, "taud", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);

  for(j=0; j<num_params_trend; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "trend", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  }

  i++;
  fprintf(fp, str_fmt, i, "A", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);

  i++;
  fprintf(fp, str_fmt, i, "Ag", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);

  for(j=0; j<num_params_difftrend; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "diff trend", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  }

  for(j=0; j<parset.n_con_recon; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "time series", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  }

  fclose(fp);
}

/*!
 * this function generates a sample from prior.
 */
void from_prior_sa2d(void *model)
{
  int i;
  double *pm = (double *)model;

  for(i=0; i<num_params; i++)
  {
    if(par_prior_model[i] == GAUSSIAN)
    {
      pm[i] = dnest_randn()*par_prior_gaussian[i][1] + par_prior_gaussian[i][0];
      dnest_wrap(&pm[i], par_range_model[i][0], par_range_model[i][1]);
    }
    else
    {
      pm[i] = par_range_model[i][0] + dnest_rand() * (par_range_model[i][1] - par_range_model[i][0]);
    }
  }

  /* cope with fixed parameters */
  for(i=0; i<num_params; i++)
  {
    if(par_fix[i] == 1)
      pm[i] = par_fix_val[i];
  }

  which_parameter_update = -1;

  return;
}

/*!
 * this function calculates likelihood at initial step.
 */
double log_likelihoods_cal_initial_sa2d(const void *model)
{
  double logL;
  logL = prob_initial_sa2d(model);
  return logL;
}

/*!
 * this function calculates likelihood at initial step.
 */
double log_likelihoods_cal_restart_sa2d(const void *model)
{
  double logL;
  logL = prob_initial_sa2d(model);
  return logL;
}

/*!
 * this function calculates likelihood.
 */
double log_likelihoods_cal_sa2d(const void *model)
{
  double logL;
  logL = prob_sa2d(model);
  return logL;
}

/*!
 * this function prints out parameters.
 */
void print_particle_sa2d(FILE *fp, const void *model)
{
  int i;
  double *pm = (double *)model;

  for(i=0; i<num_params; i++)
  {
    fprintf(fp, "%e ", pm[i] );
  }
  fprintf(fp, "\n");
  return;
}

/*!
 * This function read the particle from the file.
 */
void read_particle_sa2d(FILE *fp, void *model)
{
  int j;
  double *psample = (double *)model;

  for(j=0; j < num_params; j++)
  {
    if(fscanf(fp, "%lf", psample+j) < 1)
    {
      printf("%f\n", *psample);
      fprintf(stderr, "#Error: Cannot read file sample_sa2d.txt.\n");
      exit(0);
    }
  }
  return;
}

/*!
 * this function perturbs parameters.
 */
double perturb_sa2d(void *model)
{
  double *pm = (double *)model;
  double logH = 0.0, limit1, limit2, width, rnd;
  int which, which_level, size_levels, count_saves;

  /*
   * fixed parameters need not to update
   * perturb important parameters more frequently
   */
  do
  {
    rnd = dnest_rand();
    if(rnd < rnd_frac)
      which = dnest_rand_int(num_params_blr_tot + num_params_var);
    else
      which = dnest_rand_int(parset.n_con_recon) + num_params_blr_tot + num_params_var;

  }while(par_fix[which] == 1);

  which_parameter_update = which;

  /* level-dependent width */
  count_saves = dnest_get_count_saves();
  which_level_update = dnest_get_which_level_update();
  size_levels = dnest_get_size_levels();
  which_level = which_level_update > (size_levels-10)?(size_levels-10):which_level_update;

  if( which_level > 0)
  {
    limit1 = limits[(which_level-1) * num_params *2 + which *2];
    limit2 = limits[(which_level-1) * num_params *2 + which *2 + 1];
    width = (limit2 - limit1);
  }
  else
  {
    limit1 = par_range_model[which][0];
    limit2 = par_range_model[which][1];
    width = (par_range_model[which][1] - par_range_model[which][0]);
  }
  width /= (2.35);

  if(par_prior_model[which] == GAUSSIAN)
  {
    logH -= (-0.5*pow((pm[which] - par_prior_gaussian[which][0])/par_prior_gaussian[which][1], 2.0) );
    pm[which] += dnest_randh() * width;
    dnest_wrap(&pm[which], par_range_model[which][0], par_range_model[which][1]);
    //dnest_wrap(&pm[which], limit1, limit2);
    logH += (-0.5*pow((pm[which] - par_prior_gaussian[which][0])/par_prior_gaussian[which][1], 2.0) );
  }
  else
  {
    pm[which] += dnest_randh() * width;
    dnest_wrap(&(pm[which]), par_range_model[which][0], par_range_model[which][1]);
    //dnest_wrap(&pm[which], limit1, limit2);
  }

  return logH;
}

/*!
 * this function calculates likelihood.
 */
double log_likelihoods_cal_sa2d_exam(const void *model)
{
  return 0.0;
}

void accept_action_sa2d()
{
  int param;
  double *ptemp;

  // the parameter previously updated
  param = which_parameter_update;

  /* continuum parameter is updated */
  if( param >= num_params_blr_tot )
  {
    /*
     *note that (response) Fline is also changed as long as Fcon is changed.
     *num_params_blr-the parameter is the systematic error of continuum.
     *the change of this parameter also changes continuum reconstruction.
     */

    ptemp = Fcon_rm_particles[which_particle_update];
    Fcon_rm_particles[which_particle_update] = Fcon_rm_particles_perturb[which_particle_update];
    Fcon_rm_particles_perturb[which_particle_update] = ptemp;
  }
  else if((
           param < num_params_blr_model ||
          (param >= num_params_blr && param < num_params_blr_tot)
          ) && force_update != 1 )
  {
    /* RM and SA BLR parameter is updated
     * Note a) that the (num_par_blr-1)-th parameter is systematic error of line.
     * when this parameter is updated, Trans2D, Fline, phase_sa, and Fline_sa are unchanged.
     *      b) Fline is always changed, except param = num_params_blr-1.
     */

    ptemp = TransTau_particles[which_particle_update];
    TransTau_particles[which_particle_update] = TransTau_particles_perturb[which_particle_update];
    TransTau_particles_perturb[which_particle_update] = ptemp;

    ptemp = Trans2D_at_veldata_particles[which_particle_update];
    Trans2D_at_veldata_particles[which_particle_update] = Trans2D_at_veldata_particles_perturb[which_particle_update];
    Trans2D_at_veldata_particles_perturb[which_particle_update] = ptemp;

    prob_sa_particles[which_particle_update] = prob_sa_particles_perturb[which_particle_update];
  }

  //Fline is always updated except param = num_params_blr - 1
  if( param != num_params_blr-1 && force_update != 1)
  {
    ptemp = Fline_at_data_particles[which_particle_update];
    Fline_at_data_particles[which_particle_update] = Fline_at_data_particles_perturb[which_particle_update];
    Fline_at_data_particles_perturb[which_particle_update] = ptemp;
  }
  return;
}

/*
 * action when particle i is killed in cdnest sampling.
 * particle i_copy's properties is copyed to particle i.
 */
void kill_action_sa2d(int i, int i_copy)
{
  memcpy(Fcon_rm_particles[i], Fcon_rm_particles[i_copy], parset.n_con_recon * sizeof(double));
  memcpy(Fline_at_data_particles[i], Fline_at_data_particles[i_copy], n_line_data * n_vel_data_ext * sizeof(double));
  memcpy(TransTau_particles[i], TransTau_particles[i_copy], parset.n_tau*sizeof(double));
  memcpy(Trans2D_at_veldata_particles[i], Trans2D_at_veldata_particles[i_copy], parset.n_tau * n_vel_data_ext * sizeof(double));

  prob_sa_particles[i] = prob_sa_particles[i_copy];

  return;
}

#endif

ext: c

max_stars_repo_name: LiyrAstroph/BRAINS

max_issues_repo_name: LiyrAstroph/BRAINS

max_forks_repo_name: LiyrAstroph/BRAINS


text:
/* -*- c -*- */
/*
 * Copyright (C) 2017 Free Software Foundation, Inc.
 *
 * This file is part of GNU Radio
 *
 * GNU Radio is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3, or (at your option)
 * any later version.
 *
 * GNU Radio is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with GNU Radio; see the file COPYING.  If not, write to
 * the Free Software Foundation, Inc., 51 Franklin Street,
 * Boston, MA 02110-1301, USA.
 */

/*
 * generate MMSE interpolating differentiator FIR table values
 */

#include <math.h>
#include <assert.h>
#include <gsl/gsl_integration.h>

#define	MU	0.5			/* the MU for which we're computing coeffs */

#define	Ts	(1.0)			/* sampling period 		*/
#define	B	(1.0/(4*Ts))		/* one-sided signal bandwidth 	*/
//#define	B	(1.0/(8./3*Ts))		/* one-sided signal bandwidth 	*/

static unsigned global_n;
static double  *global_h;
double		global_mu = MU;
double		global_B  = B;

gsl_integration_workspace *global_gsl_int_workspace = NULL;

/*
 * This function computes the difference squared between the ideal
 * interpolating differentiator frequency response at frequency OMEGA and the
 * approximation defined by the FIR coefficients in global_h[]
 *
 * See eqn (9-7), "Digital Communication Receivers", Meyr, Moeneclaey
 * and Fechtel, Wiley, 1998. for the jist of what going on here.
 *
 * See eqns (7.99) and (7.100), "Discrete Time Signal Processing",
 * Oppenheim, Alan V., Schafer, Ronald W., Prentice Hall 1989. for the
 * (non-bandlimited) ideal differentiator.
 *
 * Andy Walls derived the expression for the h[n] for the bandlimited
 * differentiator by performing an Inverse DTFT on the bandlimited ideal
 * response. (Integration by parts.)
 */

static double
integrand (double omega, void *params)
{
  double real_ideal;
  double real_approx;
  double real_diff;
  double imag_ideal;
  double imag_approx;
  double imag_diff;

  int	 i, n;
  int	 I1;
  double *h;

  /*
   * Ideal differentiator frequency response,
   * with freq modulation to effect a time shift left (advance),
   * and band limited to (-fc, fc).
   * N.B. fc is global_B, mu is global_mu, Fs is 1/Ts, 2*pi*f is omega.
   * fc and Fs are in Hz, and omega is in radians/second.
   *
   * H(f) = j*2*pi*f*Ts * exp(j*2*pi*f*Ts*mu)  for |f| <= fc <= Fs/2
   */
  real_ideal = -omega * Ts * sin (omega * Ts * global_mu);
  imag_ideal =  omega * Ts * cos (omega * Ts * global_mu);

  n = global_n;
  h = global_h;
  I1 = -(n / 2);

  real_approx = 0;
  imag_approx = 0;

  for (i = 0; i < n; i++){
    real_approx += h[i] * cos (-omega * Ts * (i + I1));
    imag_approx += h[i] * sin (-omega * Ts * (i + I1));
  }

  real_diff = real_ideal - real_approx;
  imag_diff = imag_ideal - imag_approx;

  return real_diff * real_diff + imag_diff * imag_diff;
}

/*
 * Integrate the difference squared over all frequencies of interest.
 */
double
c_fcn (double *x, int n)
{
  gsl_function F;
  double result, error;

  F.function = integrand;
  F.params = NULL;

  assert ((n & 1) == 0);	/* assert n is even */
  global_n = n;
  global_h = x;
  // global_B is fc, the cutoff frequency in Hz.  It defaults to Fs/4.0.
  gsl_integration_qag(&F, -2 * M_PI * global_B, 2 * M_PI * global_B,
                      0.0, 1e-12, 1000, GSL_INTEG_GAUSS61,
                      global_gsl_int_workspace, &result, &error);
  return result;
}

/* this is the interface expected by the calling fortran code */

double
objective (double x[], int *ndim)
{
  return c_fcn (x, *ndim);
}

/*
 * starting guess for optimization
 */
void initpt (double x[], int ndim)
{
  int i;
  double w_c;
  double start;
  double t;
  double arg;
  double denom;

  /* Bandlimited, time shifted, differentiator filter (infinitely long)
   *
   * Andy Walls derived the expression for the h[n] for the bandlimited
   * differentiator by performing an Inverse DTFT on the bandlimited ideal
   * response. (Integration by parts.)
   *
   * Ideal differentiator frequency response,
   * with freq modulation to effect a time shift left (advance),
   * and band limited to (-fc, fc).
   * N.B. fc is global_B, mu is global_mu, Fs is 1/Ts.
   * fc and Fs are in Hz, and omega is in radians/second.
   *
   *   H(f) = j*2*pi*f*Ts * exp(j*2*pi*f*Ts*mu)  for |f| <= fc <= Fs/2
   *   H(f) = 0                                  for |f| >  fc <= Fs/2
   *
   * Substituting normalized radian frequency, w, for 2*pi*f*Ts, we get:
   *
   *   H(w) = j*w * exp(j*w*mu)  for |w| <= w_c = 2*pi*fc*Ts <= pi
   *   H(w) = 0                  for |w| >  w_c = 2*pi*fc*Ts <= pi
   *
   * and the Inverse Discrete Time Fourier Transform is:
   *
   *   h[n] = (1/(2*pi)) * integral(-pi,pi)[ H(w)*exp(j*w*n)*dw ]
   *
   * Performing integration by parts of the above expression, using these hints:
   *   integral [u*dv] = u*v - integral [v*du]
   *   u = j*w
   *   dv = exp(j*w*(n+mu))*dw
   *
   * We eventually get an infinitely long, but decaying impulse response:
   *
   * h[n] = (w_c/pi)/(n+mu) * [cos(pi*(w_c/pi)*(n+mu)) - sinc((w_c/pi)*(n+mu))]
   *
   *      = [ w_c*(n+mu)*cos(w_c*(n+mu)) - sin(w_c*(n+mu)) ] / (pi*(n+mu)^2)
   *
   * w_c is the normalized radian cutoff frequency:
   * 2*pi*fc*Ts, with fc <= Fs/2 and Fs = 1/Ts, so w_c <= pi
   *
   * For the case of w_c = pi and mu = 0, the cos() term alternates
   * between +/- 1 and the sinc() term goes to 0 except at n == 0, so
   * we get the following:
   *
   *    h[n] = (-1)^n / n for n != 0
   *         = 0          for n == 0
   * or
   *
   *    h[n] = ..., 1/5, -1/4, 1/3, -1/2, 1, 0, -1, 1/2, -1/3, 1/4, -1/5, ...
   *
   * Compare with:
   * Eqns (7.99) and (7.100), "Discrete Time Signal Processing",
   * Oppenheim, Alan V., Schafer, Ronald W., Prentice Hall 1989. for the
   * non-bandlimited (w_c = pi) ideal differentiator.
   */

  /* Truncated, bandlimited, time shifted, differentiator filter */
  w_c = 2.0 * M_PI * global_B * Ts;
  start = (double)(-ndim/2) + global_mu;

  for (i = 0; i < ndim; i++){
    t = (double)(i) + start;
    arg = w_c * t;
    denom = M_PI*t*t; /* always >= 0.0 */
    if (denom < 1e-9)
       x[i] = 0.0;
    else
       x[i] = (arg * cos(arg) - sin(arg))/denom;
  }
}


text:
/* SO3.c */

// SO(3) Lie group implemented with GNU GSL
// cf. wikipedia, "Rotation_group_SO(3)#Infinitesimal_rotations"

// EY note : 20160519 Bizarrely, matrix multiplication isn't implemented in a straightforward manner with the gsl_matrix class; see blas (which I don't know what it is yet) and dgemm for the libraries, and SO3blas.c in this repo for that implementation

#include <stdio.h>
#include <gsl/gsl_matrix.h>

// Function prototype
// gsl_matrix commute(gsl_matrix, gsl_matrix);

int main(void)
{
  gsl_matrix *L_x = gsl_matrix_calloc(3,3);
  gsl_matrix *L_y = gsl_matrix_calloc(3,3);
  gsl_matrix *L_z = gsl_matrix_calloc(3,3);

  gsl_matrix_set(L_x,1,2,-1);
  gsl_matrix_set(L_x,2,1,1);

  gsl_matrix_set(L_y,0,2,1);
  gsl_matrix_set(L_y,2,0,-1);

  gsl_matrix_set(L_z,0,1,-1);
  gsl_matrix_set(L_z,1,0,1);

  for (int i=0;i<3;i++)
    for (int j=0;j<3;j++)
      printf("m(%d,%d) = %g\n", i, j,
	     gsl_matrix_get(L_x,i,j));

  for (int i=0;i<3;i++)
    for (int j=0;j<3;j++)
      printf("m(%d,%d) = %g\n", i, j,
	     gsl_matrix_get(L_y,i,j));

 for (int i=0;i<3;i++)
   for (int j=0;j<3;j++)
      printf("m(%d,%d) = %g\n", i, j,
	     gsl_matrix_get(L_z,i,j));

  gsl_matrix_free(L_x);
  gsl_matrix_free(L_y);
  gsl_matrix_free(L_z);

  return 0;
}

//gsl_matrix commute(gsl_matrix x, gsl_matrix y)
//{

//}

max_stars_repo_name: ernestyalumni/CompPhys

max_issues_repo_name: ernestyalumni/CompPhys

max_forks_repo_name: ernestyalumni/CompPhys


text:
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_errno.h>

/* Setup */
gsl_rng *mgsl_rng_setup(int type)
{
  const gsl_rng_type *types[] = {
    gsl_rng_borosh13, gsl_rng_coveyou, gsl_rng_cmrg, gsl_rng_fishman18, gsl_rng_fishman20,
    gsl_rng_fishman2x, gsl_rng_gfsr4, gsl_rng_knuthran, gsl_rng_knuthran2, gsl_rng_knuthran2002,
    gsl_rng_lecuyer21, gsl_rng_minstd, gsl_rng_mrg, gsl_rng_mt19937, gsl_rng_mt19937_1999,
    gsl_rng_mt19937_1998, gsl_rng_r250, gsl_rng_ran0, gsl_rng_ran1, gsl_rng_ran2, gsl_rng_ran3,
    gsl_rng_rand, gsl_rng_rand48, gsl_rng_random128_bsd, gsl_rng_random128_glibc2,
    gsl_rng_random128_libc5, gsl_rng_random256_bsd, gsl_rng_random256_glibc2,
    gsl_rng_random256_libc5, gsl_rng_random32_bsd, gsl_rng_random32_glibc2,
    gsl_rng_random32_libc5, gsl_rng_random64_bsd, gsl_rng_random64_glibc2,
    gsl_rng_random64_libc5, gsl_rng_random8_bsd, gsl_rng_random8_glibc2, gsl_rng_random8_libc5,
    gsl_rng_random_bsd, gsl_rng_random_glibc2, gsl_rng_random_libc5, gsl_rng_randu,
    gsl_rng_ranf, gsl_rng_ranlux, gsl_rng_ranlux389, gsl_rng_ranlxd1, gsl_rng_ranlxd2,
    gsl_rng_ranlxs0, gsl_rng_ranlxs1, gsl_rng_ranlxs2, gsl_rng_ranmar, gsl_rng_slatec,
    gsl_rng_taus, gsl_rng_taus2, gsl_rng_taus113, gsl_rng_transputer, gsl_rng_tt800,
    gsl_rng_uni, gsl_rng_uni32, gsl_rng_vax, gsl_rng_waterman14, gsl_rng_zuf, gsl_rng_default
  };
  const gsl_rng_type *T;
  gsl_rng *r;
  gsl_rng_env_setup();
  T = types[type];
  r = gsl_rng_alloc(T);
  return r;
}

/* RNG IO*/
int mgsl_rng_fwrite(const char *filename, const gsl_rng *r)
{
  FILE *fp;
  if((fp = fopen(filename, "w")) == NULL) return GSL_EFAILED;
  if(gsl_rng_fwrite(fp, r) != GSL_SUCCESS) return GSL_EFAILED;
  fclose(fp);
  return GSL_SUCCESS;
}

int mgsl_rng_fread(const char *filename, gsl_rng *r)
{
  FILE *fp;
  if((fp = fopen(filename, "r")) == NULL) return GSL_EFAILED;
  if(gsl_rng_fread(fp, r) != GSL_SUCCESS) return GSL_EFAILED;
  fclose(fp);
  return GSL_SUCCESS;
}

max_stars_repo_name: frithnanth/raku-Math-Libgsl-Random

max_issues_repo_name: frithnanth/raku-Math-Libgsl-Random

max_forks_repo_name: frithnanth/raku-Math-Libgsl-Random


text:
#include <stdio.h>
#include <time.h>
#include <cblas.h>// libatlas-dev (3.10.2-7)

#include "include/Matrix.h"
#include "include/HelpFunctions.h"

// file output
char* RESULT_FOLDER_PATH = "results/";
struct stat ST = {0};

int main(int argc, const char **argv) {

    // VARIABLES

    clock_t start, end; // measure time
    bool calc_correct[3] = {false, false, false}; // correct calculations? (cache, avx, blas)
    double overall_times[4] = {0,0,0,0}; // save different overall times in array

    // arguments from command line, a.out = argv[0]
    unsigned int N = (unsigned int) atoi(argv[1]); // dimension (N x N)
    unsigned int BS = (unsigned int) atoi(argv[2]); // block size
    int REPETITIONS = atoi(argv[3]);

    printf("cmd args: %d %d %d\n",N,BS,REPETITIONS);

    // INIT.

    float* m1_f = createRandomizedMatrix_f(N);
    float* m2_f = createRandomizedMatrix_f(N);
    float* stdAlgorithm_F = calloc(N * N, sizeof (float));

    // TEST FOR DIMENSIONS

    if (N % AVX_VECTOR_SIZE != 0 || N < AVX_VECTOR_SIZE) {
        printf("> ERROR! dimension N (%d) has to be: N mod %d == 0\n",N,AVX_VECTOR_SIZE);
        return EXIT_FAILURE;
    }
    else if (BS % 2 != 0 || BS < 2 || BS > N) {
        printf("> ERROR! block size BS (%d) has to be: BS mod 2 == 0\n && <= N",BS);
        return EXIT_FAILURE;
    }
    else if (REPETITIONS <= 0) {
        printf("> ERROR! REPETITIONS (%d) has to be > 0 \n",REPETITIONS);
        return EXIT_FAILURE;
    }
    printHeader(N, BS, REPETITIONS);

    // === CALCULATIONS ==================================================================

    // STANDARD ALGORITHM

    printf("executing standard algorithm multiplications …\n");

    for (int i = 0; i < REPETITIONS; ++i) {
        float* result_f = calloc(N * N, sizeof (float));

        // one multiplication
        start = clock();
        standardMatrixMul_f(m1_f, m2_f, result_f, N);
        end = clock();

        // copy values from other matrix (later: comparison)
        if (i == 0) {
            for (unsigned int j = 0; j < N; ++j) {
                for (unsigned int k = 0; k < N; ++k) {
                    stdAlgorithm_F[(N * j) + k] = result_f[(N * j) + k];
                }
            }
        }

        overall_times[0] += ((double) (end - start)) / CLOCKS_PER_SEC;

        free(result_f);
    }

    // CACHE OPTIMIZED

    printf("executing cache optimized multiplications …\n");

    for (int i = 0; i < REPETITIONS; ++i) {
        float* result_cache_f = calloc(N * N, sizeof (float));

        start = clock();
        optimizedMatrixMul_f(m1_f, m2_f, result_cache_f, N, BS);
        end = clock();

        overall_times[1] += ((double) (end - start)) / CLOCKS_PER_SEC;
        calc_correct[0] = compareResultMatrices_f(stdAlgorithm_F, result_cache_f, N);

        free(result_cache_f);
    }

    // PARALLEL (AVX)

    printf("executing parallel multiplications …\n");

    for (int i = 0; i < REPETITIONS; ++i) {
        float* result_para_f = calloc(N * N, sizeof (float));

        start = clock();
        parallelMatrixMul_AVX_f(m1_f, m2_f, result_para_f, N);
        end = clock();

        overall_times[2] += ((double) (end - start)) / CLOCKS_PER_SEC;
        calc_correct[1] = compareResultMatrices_f(stdAlgorithm_F, result_para_f, N);

        free(result_para_f);
    }

    // BLAS (ATLAS)

    printf("executing multiplications using BLAS …\n");
    float ALPHA = 1;
    float BETA = 0;

    // matrix dimensions
    int LDA=N;
    int LDB=N;
    int LDC=N;

    // Single Precision, General Matrix Multiplication
    for (int i = 0; i < REPETITIONS; ++i) {
        float* result_cblas_f = calloc(N * N, sizeof (float));

        start = clock();
        cblas_sgemm(CblasRowMajor,CblasNoTrans,CblasNoTrans,N,N,N,ALPHA,
                    m1_f,LDA,
                    m2_f,LDB,BETA,
                    result_cblas_f,LDC);
        end = clock();

        overall_times[3] += ((double) (end - start)) / CLOCKS_PER_SEC;
        calc_correct[2] = compareResultMatrices_f(stdAlgorithm_F, result_cblas_f, N);

        free(result_cblas_f);
    }

    // ===== cleanup allocated memory ===================================

    free(m1_f);
    free(m2_f);
    free(stdAlgorithm_F);

    // === FILE OUTPUT ===================================================================

    createResultFolder(RESULT_FOLDER_PATH, ST);
    char* filename =
            createResultFile(RESULT_FOLDER_PATH, N, BS, REPETITIONS, overall_times, calc_correct);
    printf("\nwriting results to %s …\n", filename);
    printf("done.\n");

    free(filename);

    return 0;
}

max_stars_repo_name: matthiastz/FPS-Matrixmultiplikation

max_issues_repo_name: matthiastz/FPS-Matrixmultiplikation

max_forks_repo_name: matthiastz/FPS-Matrixmultiplikation


text:
/**
 *
 * @file core_zgetrf_incpiv.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.8.0
 * @author Hatem Ltaief
 * @author Mathieu Faverge
 * @author Jakub Kurzak
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup dplasma_cores_complex64
 *
 *  CORE_zgetrf_incpiv computes an LU factorization of a general M-by-N tile A
 *  using partial pivoting with row interchanges.
 *
 *  The factorization has the form
 *
 *    A = P * L * U
 *
 *  where P is a permutation matrix, L is lower triangular with unit
 *  diagonal elements (lower trapezoidal if m > n), and U is upper
 *  triangular (upper trapezoidal if m < n).
 *
 *  This is the right-looking Level 2.5 BLAS version of the algorithm.
 *
 *******************************************************************************
 *
 * @param[in] M
 *          The number of rows of the tile A.  M >= 0.
 *
 * @param[in] N
 *         The number of columns of the tile A.  N >= 0.
 *
 * @param[in] IB
 *         The inner-blocking size.  IB >= 0.
 *
 * @param[in,out] A
 *         On entry, the M-by-N tile to be factored.
 *         On exit, the factors L and U from the factorization
 *         A = P*L*U; the unit diagonal elements of L are not stored.
 *
 * @param[in] LDA
 *         The leading dimension of the array A.  LDA >= max(1,M).
 *
 * @param[out] IPIV
 *         The pivot indices; for 1 <= i <= min(M,N), row i of the
 *         tile was interchanged with row IPIV(i).
 *
 * @param[out] INFO
 *         See returned value.
 *
 *******************************************************************************
 *
 * @return
 *         \retval PLASMA_SUCCESS successful exit
 *         \retval <0 if INFO = -k, the k-th argument had an illegal value
 *         \retval >0 if INFO = k, U(k,k) is exactly zero. The factorization
 *              has been completed, but the factor U is exactly
 *              singular, and division by zero will occur if it is used
 *              to solve a system of equations.
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zgetrf_incpiv = PCORE_zgetrf_incpiv
#define CORE_zgetrf_incpiv PCORE_zgetrf_incpiv
#define CORE_zgessm PCORE_zgessm
int  CORE_zgessm(int M, int N, int K, int IB,
                 int *IPIV,
                 const PLASMA_Complex64_t *L, int LDL,
                 PLASMA_Complex64_t *A, int LDA);
#endif
int CORE_zgetrf_incpiv(int M, int N, int IB,
                       PLASMA_Complex64_t *A, int LDA,
                       int *IPIV, int *INFO)
{
    int i, j, k, sb;
    int iinfo;

    /* Check input arguments */
    *INFO = 0;
    if (M < 0) {
        coreblas_error(1, "Illegal value of M");
        return -1;
    }
    if (N < 0) {
        coreblas_error(2, "Illegal value of N");
        return -2;
    }
    if (IB < 0) {
        coreblas_error(3, "Illegal value of IB");
        return -3;
    }
    if ((LDA < max(1,M)) && (M > 0)) {
        coreblas_error(5, "Illegal value of LDA");
        return -5;
    }

    /* Quick return */
    if ((M == 0) || (N == 0) || (IB == 0))
        return PLASMA_SUCCESS;

    k = min(M, N);

    for(i =0 ; i < k; i += IB) {
        sb = min(IB, k-i);
        /*
         * Factor diagonal and subdiagonal blocks and test for exact singularity.
         */
        iinfo = LAPACKE_zgetf2_work(LAPACK_COL_MAJOR, M-i, sb, &A[LDA*i+i], LDA, &IPIV[i]);
        /*
         * Adjust INFO and the pivot indices.
         */
        if((*INFO == 0) && (iinfo > 0))
            *INFO = iinfo + i;

        if (i+sb < N) {
            CORE_zgessm(
                M-i, N-(i+sb), sb, sb,
                &IPIV[i],
                &A[LDA*i+i], LDA,
                &A[LDA*(i+sb)+i], LDA);
        }

        for(j = i; j < i+sb; j++) {
            IPIV[j] = i + IPIV[j];
        }
    }
    return PLASMA_SUCCESS;
}

max_stars_repo_name: therault/dplasma

max_issues_repo_name: therault/dplasma

max_forks_repo_name: therault/dplasma


text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#pragma once

#ifndef style_ac50d5be_0f99_48dc_908f_ab8a309b2a27_h
#define style_ac50d5be_0f99_48dc_908f_ab8a309b2a27_h

#include <gslib/string.h>
#include <gslib/std.h>
#include <ariel/type.h>
#include <ariel/painter.h>

__ariel_begin__

enum style_sheet_type
{
    sst_unknown     = -1,
    sst_color       = 0,
    sst_integer,
    sst_boolean,
    sst_float,
    sst_string,
    sst_capacity,
};

extern const string& get_style_sheet_type_name(style_sheet_type sst);
typedef std::pair<style_sheet_type, string> style_sheet_def;
typedef unordered_map<string, int> style_sheet_info_map;

struct accel_key;

class __gs_novtable style_sheet abstract
{
public:
    static const int npos = -1;

public:
    style_sheet(const style_sheet_def* ssp, int len);
    virtual ~style_sheet() {}
    virtual bool get_value(const string& name, string& value) = 0;
    virtual void set_value(const string& name, const string& value) = 0;
    virtual int get_content_size() const { return _ss_length; }
    virtual style_sheet_type get_content_type(int index) const;
    virtual const string& get_content_name(int index) const;
    virtual void flush_style() = 0;

protected:
    const style_sheet_def*  _ss_pairs;
    int                     _ss_length;
    style_sheet_info_map    _ss_info;

public:
    void initialize_style_sheet(const style_sheet_def* ssp, int len);
    int get_style_sheet_index(const string& name) const;

protected:
    static bool from_color(string& str, const color& cr);
    static bool to_color(color& cr, const string& str);
    static bool from_integer(string& str, int i);
    static bool to_integer(int& i, const string& str);
    static bool from_boolean(string& str, bool b);
    static bool to_boolean(bool& b, const string& str);
    static bool from_float(string& str, float f);
    static bool to_float(float& f, const string& str);
    static bool from_accel_key(string& str, const accel_key& k);
    static bool to_accel_key(accel_key& k, const string& str);
    static void setup_brush_by_color(painter_brush& brush, const color& cr);
    static void setup_pen_by_color(painter_pen& pen, const color& cr);
    static void setup_font(font& ft, const string& name, int size);
};

__ariel_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
//Gets cepstral coefficients (CCs). This is the same as dct + lifter.
//This computes "the" DCT (discrete cosine transformation) along one dimension of matrix X.
//This uses fftw3 to compute the DCT-II.
//Then it applies a lifter (cepstral domain window) to the same dimension of X.
//The lifter has a single parameter, Q, which is default 22.0 in HTK, Kaldi and Ellis.
//Only the first K CCs are output.

#include <stdio.h>
#include <float.h>
#include <math.h>
#include <cblas.h>
#include <fftw3.h>

#ifndef M_PIf
   #define M_PIf 3.141592653589793238462643383279502884f
#endif

#ifndef M_PI
   #define M_PI 3.141592653589793238462643383279502884
#endif

#ifdef __cplusplus
namespace ov {
extern "C" {
#endif

int get_ccs_s (float *Y, const float *X, const int iscolmajor, const int R, const int C, const int dim, const int ndct, const float Q, const int K);
int get_ccs_d (double *Y, const double *X, const int iscolmajor, const int R, const int C, const int dim, const int ndct, const double Q, const int K);

int get_ccs_s (float *Y, const float *X, const int iscolmajor, const int R, const int C, const int dim, const int ndct, const float Q, const int K)
{
    const float z = 0.0f;
    const float Q2 = 0.5f*Q, PQ = M_PIf/Q;
    int r, c, k;
    float *X1, *Y1; //*lift;
    fftwf_plan plan;
    //struct timespec tic, toc;

    //Checks
    if (R<1) { fprintf(stderr,"error in get_ccs_s: R (nrows X) must be positive\n"); return 1; }
    if (C<1) { fprintf(stderr,"error in get_ccs_s: C (ncols X) must be positive\n"); return 1; }
    if (ndct<R && dim==0) { fprintf(stderr,"error in get_ccs_s: ndct must be >= R for dim==0\n"); return 1; }
    if (ndct<C && dim==1) { fprintf(stderr,"error in get_ccs_s: ndct must be >= C for dim==1\n"); return 1; }
    if (Q<0.0f) { fprintf(stderr,"error in get_ccs_s: Q must be positive (or 0 to skip lifter)\n"); return 1; }
    if (K<1) { fprintf(stderr,"error in get_ccs_s: K must be positive\n"); return 1; }
    if (K>ndct) { fprintf(stderr,"error in get_ccs_s: K must be <= ndct>\n"); return 1; }

    //Initialize fftwf
    X1 = fftwf_alloc_real((size_t)ndct);
    Y1 = fftwf_alloc_real((size_t)ndct);
    plan = fftwf_plan_r2r_1d(ndct,X1,Y1,FFTW_REDFT10,FFTW_ESTIMATE);
    if (!plan) { fprintf(stderr,"error in get_ccs_s: problem creating fftw plan\n"); return 1; }
    cblas_scopy(ndct,&z,0,&X1[0],1); //zero-pad

    //clock_gettime(CLOCK_REALTIME,&tic);

    //Initialize lifter
    //if (!(lift=(float *)malloc((size_t)K*sizeof(float)))) { fprintf(stderr,"error in get_ccs_s: problem with malloc for lifter. "); perror("malloc"); return 1; }
    //for (k=0; k<K; k++) { lift[k] = fmaf(Q2,sinf(k*PQ),1.0f); }

    if (dim==0)
    {
        if (iscolmajor)
        {
            for (c=0; c<C; c++)
            {
                cblas_scopy(R,&X[c*R],1,&X1[0],1);
                fftwf_execute(plan);
                //cblas_ssbmv(CblasColMajor,CblasUpper,K,0,1.0f,&lift[0],1,&Y1[0],1,0.0f,&Y[c*K],1);
                cblas_scopy(K,&Y1[0],1,&Y[c*K],1);
            }
            if (Q>FLT_EPSILON)
            {
                for (k=0; k<K; k++) { cblas_sscal(C,fmaf(Q2,sinf(PQ*k),1.0f),&Y[k],K); }
            }
        }
        else
        {
            for (c=0; c<C; c++)
            {
                cblas_scopy(R,&X[c],C,&X1[0],1);
                fftwf_execute(plan);
                //cblas_ssbmv(CblasRowMajor,CblasUpper,K,0,1.0f,&lift[0],1,&Y1[0],1,0.0f,&Y[c],C);
                cblas_scopy(K,&Y1[0],1,&Y[c],C);
            }
            if (Q>FLT_EPSILON)
            {
                for (k=0; k<K; k++) { cblas_sscal(C,fmaf(Q2,sinf(PQ*k),1.0f),&Y[k*C],1); }
            }
        }
    }
    else if (dim==1)
    {
        if (iscolmajor)
        {
            for (r=0; r<R; r++)
            {
                cblas_scopy(C,&X[r],R,&X1[0],1);
                fftwf_execute(plan);
                //cblas_ssbmv(CblasColMajor,CblasUpper,K,0,1.0f,&lift[0],1,&Y1[0],1,0.0f,&Y[r],R);
                cblas_scopy(K,&Y1[0],1,&Y[r],R);
            }
            if (Q>FLT_EPSILON)
            {
                for (k=0; k<K; k++) { cblas_sscal(R,fmaf(Q2,sinf(PQ*k),1.0f),&Y[k*R],1); }
            }
        }
        else
        {
            for (r=0; r<R; r++)
            {
                cblas_scopy(C,&X[r*C],1,&X1[0],1);
                fftwf_execute(plan);
                //cblas_ssbmv(CblasRowMajor,CblasUpper,K,0,1.0f,&lift[0],1,&Y1[0],1,0.0f,&Y[r*K],1);
                cblas_scopy(K,&Y1[0],1,&Y[r*K],1);
            }
            if (Q>FLT_EPSILON)
            {
                for (k=0; k<K; k++) { cblas_sscal(R,fmaf(Q2,sinf(PQ*k),1.0f),&Y[k],K); }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in get_ccs_s: dim must be 0 or 1.\n"); return 1;
    }
    //clock_gettime(CLOCK_REALTIME,&toc); fprintf(stderr,"elapsed time = %f ms\n",(toc.tv_sec-tic.tv_sec)*1e3+(toc.tv_nsec-tic.tv_nsec)/1e6);

    fftwf_destroy_plan(plan); fftwf_free(X1); fftwf_free(Y1);
    return 0;
}

int get_ccs_d (double *Y, const double *X, const int iscolmajor, const int R, const int C, const int dim, const int ndct, const double Q, const int K)
{
    const double z = 0.0;
    const double Q2 = 0.5*Q, PQ = M_PI/Q;
    int r, c, k;
    double *X1, *Y1;
    fftw_plan plan;

    //Checks
    if (R<1) { fprintf(stderr,"error in get_ccs_d: R (nrows X) must be positive\n"); return 1; }
    if (C<1) { fprintf(stderr,"error in get_ccs_d: C (ncols X) must be positive\n"); return 1; }
    if (ndct<R && dim==0) { fprintf(stderr,"error in get_ccs_d: ndct must be >= R for dim==0\n"); return 1; }
    if (ndct<C && dim==1) { fprintf(stderr,"error in get_ccs_d: ndct must be >= C for dim==1\n"); return 1; }
    if (Q<0.0) { fprintf(stderr,"error in get_ccs_d: Q must be positive (or 0 to skip lifter)\n"); return 1; }
    if (K<1) { fprintf(stderr,"error in get_ccs_d: K must be positive\n"); return 1; }
    if (K>ndct) { fprintf(stderr,"error in get_ccs_d: K must be <= ndct>\n"); return 1; }

    //Initialize fftw
    X1 = fftw_alloc_real((size_t)ndct);
    Y1 = fftw_alloc_real((size_t)ndct);
    plan = fftw_plan_r2r_1d(ndct,X1,Y1,FFTW_REDFT10,FFTW_ESTIMATE);
    if (!plan) { fprintf(stderr,"error in get_ccs_s: problem creating fftw plan\n"); return 1; }
    cblas_dcopy(ndct,&z,0,&X1[0],1); //zero-pad

    if (dim==0)
    {
        if (iscolmajor)
        {
            for (c=0; c<C; c++)
            {
                cblas_dcopy(R,&X[c*R],1,&X1[0],1);
                fftw_execute(plan);
                cblas_dcopy(K,&Y1[0],1,&Y[c*K],1);
            }
            if (Q>DBL_EPSILON)
            {
                for (k=0; k<K; k++) { cblas_dscal(C,fma(Q2,sin(PQ*k),1.0),&Y[k],K); }
            }
        }
        else
        {
            for (c=0; c<C; c++)
            {
                cblas_dcopy(R,&X[c],C,&X1[0],1);
                fftw_execute(plan);
                cblas_dcopy(K,&Y1[0],1,&Y[c],C);
            }
            if (Q>DBL_EPSILON)
            {
                for (k=0; k<K; k++) { cblas_dscal(C,fma(Q2,sin(PQ*k),1.0),&Y[k*C],1); }
            }
        }
    }
    else if (dim==1)
    {
        if (iscolmajor)
        {
            for (r=0; r<R; r++)
            {
                cblas_dcopy(C,&X[r],R,&X1[0],1);
                fftw_execute(plan);
                cblas_dcopy(K,&Y1[0],1,&Y[r],R);
            }
            if (Q>DBL_EPSILON)
            {
                for (k=0; k<K; k++) { cblas_dscal(R,fma(Q2,sin(PQ*k),1.0),&Y[k*R],1); }
            }
        }
        else
        {
            for (r=0; r<R; r++)
            {
                cblas_dcopy(C,&X[r*C],1,&X1[0],1);
                fftw_execute(plan);
                cblas_dcopy(K,&Y1[0],1,&Y[r*K],1);
            }
            if (Q>DBL_EPSILON)
            {
                for (k=0; k<K; k++) { cblas_dscal(R,fma(Q2,sin(PQ*k),1.0),&Y[k],K); }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in get_ccs_d: dim must be 0 or 1.\n"); return 1;
    }

    fftw_destroy_plan(plan); fftw_free(X1); fftw_free(Y1);
    return 0;
}

#ifdef __cplusplus
}
}
#endif


text:
#include <stdio.h>
#include <gsl/gsl_linalg.h>

double A[] =
{
  -1.0, 0.0, 2.0, 1.0,
  0.0, -1.0, 3.0, 2.0,
  1.0, 2.0, 3.0, 0.0,
  -3.0, 1.0, 4.0, 2.0
};

double y[] = { 0.0, 1.0, 2.0, 1.0 };

int main() {
  gsl_matrix_view m = gsl_matrix_view_array(A, 4, 4);
  gsl_vector_view b = gsl_vector_view_array(y, 4);
  gsl_vector* x = gsl_vector_alloc(4);

  int s;
  gsl_permutation* p = gsl_permutation_alloc(4);
  gsl_linalg_LU_decomp(&m.matrix, p, &s);
  gsl_linalg_LU_solve(&m.matrix, p, &b.vector, x);

  printf("x = \n");
  gsl_vector_fprintf(stdout, x, "%5.1f");

  gsl_permutation_free(p);
  gsl_vector_free(x);
  return 0;
}

max_stars_repo_name: julnamoo/practice-linux

max_issues_repo_name: julnamoo/practice-linux

max_forks_repo_name: julnamoo/practice-linux