text:
/*****************************************************************\
           __
          / /
		 / /                     __  __
		/ /______    _______    / / / / ________   __       __
	   / ______  \  /_____  \  / / / / / _____  | / /      / /
	  / /      | / _______| / / / / / / /____/ / / /      / /
	 / /      / / / _____  / / / / / / _______/ / /      / /
	/ /      / / / /____/ / / / / / / |______  / |______/ /
   /_/      /_/ |________/ / / / /  \_______/  \_______  /
                          /_/ /_/                     / /
			                                         / /
		       High Level Game Framework            /_/

  ---------------------------------------------------------------

  Copyright (c) 2007-2011 - Rodrigo Braz Monteiro.
  This file is subject to the terms of halley_license.txt.

\*****************************************************************/

#pragma once

#include <algorithm>
#include <limits>
#include <list>
#include <cmath>

#include <memory>
#include <array>
#include <functional>
#include <halley/data_structures/vector.h>
#include <gsl/gsl_assert>

#ifdef _MSC_VER
#include <xmmintrin.h>
#endif

namespace Halley {
	// General aliases
	using Byte = unsigned char;
	using Bytes = Vector<Byte>;

	template <typename T>
	constexpr inline T clamp(T value, T minValue, T maxValue)
	{
		return std::min(std::max(minValue, value), maxValue);
	}

	template <typename T>
	constexpr inline T clamp2(T value, T minValue, T maxValue)
	{
		return std::max(minValue, std::min(value, maxValue));
	}

	template <typename T>
	constexpr inline T maxAbs(T a, T b)
	{
		return abs(a) > abs(b) ? a : b;
	}

	template <typename T>
	constexpr inline T minAbs(T a, T b)
	{
		return abs(a) < abs(b) ? a : b;
	}

	template <typename T>
	constexpr inline bool rangeIntersection(T s1, T e1, T s2, T e2)
	{
		return (s1 < e2) && (s2 < e1);
	}

	constexpr double pi()
	{
		return 3.1415926535897932384626433832795;
	}

	// Endianness conversion
	//const bool IS_BIG_ENDIAN = *((short*)"AB") == 0x4243;
	constexpr bool IS_BIG_ENDIAN = false;

	template <typename T>
	constexpr inline T FixEndian(T t)
	{
		if (IS_BIG_ENDIAN) {
			T temp = t;
			unsigned char *v = reinterpret_cast<unsigned char*>(&temp);
			for (size_t i=0; i<sizeof(T)/2; i++) {
				std::swap(v[i], v[sizeof(T)-1-i]);
			}
			return temp;
		} else {
			return t;
		}
	}

	template<>
	constexpr inline char FixEndian<char>(char t) {
		return t;
	}

	template<>
	constexpr inline unsigned char FixEndian<unsigned char>(unsigned char t) {
		return t;
	}

	// True modulo definition
	template <typename T> constexpr inline T modulo (T a, T b)
	{
		static_assert(std::is_signed<T>::value, "Must be signed to use modulo operation");
		T res = a % b;
		if (res < 0) res = b+res;
		return res;
	}

	template <typename T> constexpr inline T floatModulo (T a, T b)
	{
		if (b == 0) return a;
		return a - b * std::floor(a/b);
	}
	template <> constexpr inline float modulo (float a, float b) { return floatModulo(a, b); }
	template <> constexpr inline double modulo (double a, double b) { return floatModulo(a, b); }

	// Float floor division (e.g. -0.5 / 2 = -1.0)
	template <typename T>
	constexpr inline T floorDivFloat (T a, T b)
	{
		return floor(a / b);
	}

	// Int floor division (e.g. -2 / 3 = -1)
	template <typename T>
	constexpr inline T floorDivInt (T a, T b)
	{
		T result = a / b;
		if (a % b < 0) {
			--result;
		}
		return result;
	}

	// Generic floor divs
	constexpr inline float floorDiv (float a, float b) { return floorDivFloat(a, b); }
	constexpr inline double floorDiv (double a, double b) { return floorDivFloat(a, b); }
	constexpr inline long floorDiv (long a, long b) { return floorDivInt(a, b); }
	constexpr inline int floorDiv (int a, int b) { return floorDivInt(a, b); }
	constexpr inline short floorDiv (short a, short b) { return floorDivInt(a, b); }
	constexpr inline char floorDiv (char a, char b) { return floorDivInt(a, b); }

	// Interpolation
	template <typename T>
	constexpr inline T interpolate(T a, T b, float factor) {
		return static_cast<T>(a * (1 - factor) + b * factor);
	}

	template <typename T>
	constexpr inline T lerp(T a, T b, float factor) {
		return static_cast<T>(a * (1 - factor) + b * factor);
	}

	// Smoothing
	template <typename T>
	constexpr inline T smoothCos(T a) {
		return T((1-cos(a * 3.1415926535897932384626433832795))*0.5);
	}

	// ASR (attack-sustain-release) envelope
	// Returns 0 at start of attack and end of release, 1 during sustain
	// Attack and release are linear interpolations
	template <typename T>
	constexpr inline T asr(T x, T a, T s, T r) {
		if (x < a) return x/a;
		T as = a+s;
		if (x < as) return 1;
		return 1 - ((x-as)/r);
	}

	// Next power of 2
	template<typename T>
	constexpr static T nextPowerOf2(T val)
	{
		--val;
		val = (val >> 1) | val;
		val = (val >> 2) | val;
		val = (val >> 4) | val;
		val = (val >> 8) | val;
		val = (val >> 16) | val;
		return val+1;
	}

	constexpr inline int fastLog2Floor (uint32_t value) {
		// From https://stackoverflow.com/questions/11376288/fast-computing-of-log2-for-64-bit-integers
		constexpr int tab32[32] = { 0, 9, 1, 10, 13, 21, 2, 29, 11, 14, 16, 18, 22, 25, 3, 30, 8, 12, 20, 28, 15, 17, 24, 7, 19, 27, 23, 6, 26, 5, 4, 31 };

		Expects(value > 0);
		value |= value >> 1;
		value |= value >> 2;
		value |= value >> 4;
		value |= value >> 8;
		value |= value >> 16;
		return tab32[uint32_t(value * 0x07C4ACDD) >> 27];
	}

	constexpr inline int fastLog2Floor (uint64_t value)
	{
		// From https://stackoverflow.com/questions/11376288/fast-computing-of-log2-for-64-bit-integers
		constexpr int tab64[64] = { 63,  0, 58,  1, 59, 47, 53,  2, 60, 39, 48, 27, 54, 33, 42,  3, 61, 51, 37, 40, 49, 18, 28, 20,  55, 30, 34, 11, 43, 14, 22,  4, 62, 57, 46, 52, 38, 26, 32, 41, 50, 36, 17, 19, 29, 10, 13, 21, 56, 45, 25, 31, 35, 16,  9, 12, 44, 24, 15,  8, 23,  7,  6,  5};

		Expects(value > 0);
		value |= value >> 1;
		value |= value >> 2;
		value |= value >> 4;
		value |= value >> 8;
		value |= value >> 16;
		value |= value >> 32;
		return tab64[uint64_t((value - (value >> 1)) * 0x07EDD5E59A4E28C2) >> 58];
	}

	constexpr inline int fastLog2Ceil (uint32_t value)
	{
		return fastLog2Floor(value - 1) + 1;
	}

	constexpr inline int fastLog2Ceil (uint64_t value)
	{
		return fastLog2Floor(value - 1) + 1;
	}

	// Advance a to b by up to inc
	template<typename T>
	constexpr static T advance(T a, T b, T inc)
	{
		if (a < b) return std::min(a+inc, b);
		else return std::max(a-inc, b);
	}

	// Align address
	template <typename T>
	constexpr T alignUp(T val, T align)
	{
		return val + (align - (val % align)) % align;
	}

	template <typename T>
	constexpr T alignDown(T val, T align)
	{
		return (val / align) * align;
	}

	// Prefetch data from memory
	static inline void prefetchL1(void* p) {
#ifdef _MSC_VER
		_mm_prefetch(static_cast<const char*>(p), _MM_HINT_T0);
#else
		__builtin_prefetch(p);
#endif
	}
	static inline void prefetchL2(void* p) {
#ifdef _MSC_VER
		_mm_prefetch(static_cast<const char*>(p), _MM_HINT_T1);
#else
		__builtin_prefetch(p);
#endif
	}
}

max_stars_repo_name: Healthire/halley

max_issues_repo_name: Healthire/halley

max_forks_repo_name: Healthire/halley


text:
/* vector/gsl_vector_complex_double.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
 *
 * This program 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 of the License, or (at
 * your option) any later version.
 *
 * This program 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 this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_VECTOR_COMPLEX_DOUBLE_H__
#define __GSL_VECTOR_COMPLEX_DOUBLE_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_double.h>
#include <gsl/gsl_vector_complex.h>
#include <gsl/gsl_block_complex_double.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t size;
  size_t stride;
  double *data;
  gsl_block_complex *block;
  int owner;
} gsl_vector_complex;

typedef struct
{
  gsl_vector_complex vector;
} _gsl_vector_complex_view;

typedef _gsl_vector_complex_view gsl_vector_complex_view;

typedef struct
{
  gsl_vector_complex vector;
} _gsl_vector_complex_const_view;

typedef const _gsl_vector_complex_const_view gsl_vector_complex_const_view;

/* Allocation */

GSL_FUN gsl_vector_complex *gsl_vector_complex_alloc (const size_t n);
GSL_FUN gsl_vector_complex *gsl_vector_complex_calloc (const size_t n);

GSL_FUN gsl_vector_complex *
gsl_vector_complex_alloc_from_block (gsl_block_complex * b,
                                           const size_t offset,
                                           const size_t n,
                                           const size_t stride);

GSL_FUN gsl_vector_complex *
gsl_vector_complex_alloc_from_vector (gsl_vector_complex * v,
                                             const size_t offset,
                                             const size_t n,
                                             const size_t stride);

GSL_FUN void gsl_vector_complex_free (gsl_vector_complex * v);

/* Views */

GSL_FUN _gsl_vector_complex_view
gsl_vector_complex_view_array (double *base,
                                     size_t n);

GSL_FUN _gsl_vector_complex_view
gsl_vector_complex_view_array_with_stride (double *base,
                                                 size_t stride,
                                                 size_t n);

GSL_FUN _gsl_vector_complex_const_view
gsl_vector_complex_const_view_array (const double *base,
                                           size_t n);

GSL_FUN _gsl_vector_complex_const_view
gsl_vector_complex_const_view_array_with_stride (const double *base,
                                                       size_t stride,
                                                       size_t n);

GSL_FUN _gsl_vector_complex_view
gsl_vector_complex_subvector (gsl_vector_complex *base,
                                         size_t i,
                                         size_t n);

GSL_FUN _gsl_vector_complex_view
gsl_vector_complex_subvector_with_stride (gsl_vector_complex *v,
                                                size_t i,
                                                size_t stride,
                                                size_t n);

GSL_FUN _gsl_vector_complex_const_view
gsl_vector_complex_const_subvector (const gsl_vector_complex *base,
                                               size_t i,
                                               size_t n);

GSL_FUN _gsl_vector_complex_const_view
gsl_vector_complex_const_subvector_with_stride (const gsl_vector_complex *v,
                                                      size_t i,
                                                      size_t stride,
                                                      size_t n);

GSL_FUN _gsl_vector_view
gsl_vector_complex_real (gsl_vector_complex *v);

GSL_FUN _gsl_vector_view
gsl_vector_complex_imag (gsl_vector_complex *v);

GSL_FUN _gsl_vector_const_view
gsl_vector_complex_const_real (const gsl_vector_complex *v);

GSL_FUN _gsl_vector_const_view
gsl_vector_complex_const_imag (const gsl_vector_complex *v);

/* Operations */

GSL_FUN void gsl_vector_complex_set_zero (gsl_vector_complex * v);
GSL_FUN void gsl_vector_complex_set_all (gsl_vector_complex * v,
                                       gsl_complex z);
GSL_FUN int gsl_vector_complex_set_basis (gsl_vector_complex * v, size_t i);

GSL_FUN int gsl_vector_complex_fread (FILE * stream,
                                    gsl_vector_complex * v);
GSL_FUN int gsl_vector_complex_fwrite (FILE * stream,
                                     const gsl_vector_complex * v);
GSL_FUN int gsl_vector_complex_fscanf (FILE * stream,
                                     gsl_vector_complex * v);
GSL_FUN int gsl_vector_complex_fprintf (FILE * stream,
                                      const gsl_vector_complex * v,
                                      const char *format);

GSL_FUN int gsl_vector_complex_memcpy (gsl_vector_complex * dest, const gsl_vector_complex * src);

GSL_FUN int gsl_vector_complex_reverse (gsl_vector_complex * v);

GSL_FUN int gsl_vector_complex_swap (gsl_vector_complex * v, gsl_vector_complex * w);
GSL_FUN int gsl_vector_complex_swap_elements (gsl_vector_complex * v, const size_t i, const size_t j);

GSL_FUN int gsl_vector_complex_isnull (const gsl_vector_complex * v);
GSL_FUN int gsl_vector_complex_ispos (const gsl_vector_complex * v);
GSL_FUN int gsl_vector_complex_isneg (const gsl_vector_complex * v);
GSL_FUN int gsl_vector_complex_isnonneg (const gsl_vector_complex * v);

GSL_FUN int gsl_vector_complex_add (gsl_vector_complex * a, const gsl_vector_complex * b);
GSL_FUN int gsl_vector_complex_sub (gsl_vector_complex * a, const gsl_vector_complex * b);
GSL_FUN int gsl_vector_complex_mul (gsl_vector_complex * a, const gsl_vector_complex * b);
GSL_FUN int gsl_vector_complex_div (gsl_vector_complex * a, const gsl_vector_complex * b);
GSL_FUN int gsl_vector_complex_scale (gsl_vector_complex * a, const gsl_complex x);
GSL_FUN int gsl_vector_complex_add_constant (gsl_vector_complex * a, const gsl_complex x);

GSL_FUN INLINE_DECL gsl_complex gsl_vector_complex_get (const gsl_vector_complex * v, const size_t i);
GSL_FUN INLINE_DECL void gsl_vector_complex_set (gsl_vector_complex * v, const size_t i, gsl_complex z);
GSL_FUN INLINE_DECL gsl_complex *gsl_vector_complex_ptr (gsl_vector_complex * v, const size_t i);
GSL_FUN INLINE_DECL const gsl_complex *gsl_vector_complex_const_ptr (const gsl_vector_complex * v, const size_t i);

#ifdef HAVE_INLINE

INLINE_FUN
gsl_complex
gsl_vector_complex_get (const gsl_vector_complex * v,
                              const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      gsl_complex zero = {{0, 0}};
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, zero);
    }
#endif
  return *GSL_COMPLEX_AT (v, i);
}

INLINE_FUN
void
gsl_vector_complex_set (gsl_vector_complex * v,
                              const size_t i, gsl_complex z)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  *GSL_COMPLEX_AT (v, i) = z;
}

INLINE_FUN
gsl_complex *
gsl_vector_complex_ptr (gsl_vector_complex * v,
                              const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return GSL_COMPLEX_AT (v, i);
}

INLINE_FUN
const gsl_complex *
gsl_vector_complex_const_ptr (const gsl_vector_complex * v,
                                    const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return GSL_COMPLEX_AT (v, i);
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_COMPLEX_DOUBLE_H__ */

ext: h

max_stars_repo_name: berkus/music-cs

max_issues_repo_name: berkus/music-cs

max_forks_repo_name: berkus/music-cs


text:
#pragma once
#define VULKAN_H_ 1
#include <vulkan/vk_platform.h>
#include <vulkan/vulkan_core.h>

#if defined(VK_USE_PLATFORM_WIN32_KHR)
// This is to avoid including Windows.h
#ifndef DECLARE_HANDLE
using HANDLE = void*;

#define DECLARE_HANDLE(name) struct name##__; typedef struct name##__ *name
DECLARE_HANDLE(HINSTANCE);
DECLARE_HANDLE(HMONITOR);
DECLARE_HANDLE(HWND);
#undef DECLARE_HANDLE

struct _SECURITY_ATTRIBUTES;
using SECURITY_ATTRIBUTES = _SECURITY_ATTRIBUTES;

using LPCWSTR = const wchar_t*;
using PCWSTR = const wchar_t*;
using DWORD = unsigned long; // NOLINT(google-runtime-int)
#endif

// ReSharper disable once CppUnusedIncludeDirective
#include <vulkan/vulkan_win32.h>
#endif

#if defined(VK_USE_PLATFORM_XCB_KHR)
typedef struct xcb_connection_t xcb_connection_t;
typedef uint32_t xcb_window_t;
typedef uint32_t xcb_visualid_t;

// ReSharper disable once CppUnusedIncludeDirective
#include <vulkan/vulkan_xcb.h>
#endif

#if defined(VK_USE_PLATFORM_XLIB_KHR)
typedef struct _XDisplay Display;
typedef unsigned long XID;
typedef XID Window;
typedef unsigned long VisualID;

// ReSharper disable once CppUnusedIncludeDirective
#include <vulkan/vulkan_xlib.h>
#endif

#if defined(VK_USE_PLATFORM_XLIB_XRANDR_EXT)
typedef struct _XDisplay Display;
typedef unsigned long XID;
typedef XID Window;
typedef unsigned long VisualID;
typedef XID RROutput;

// ReSharper disable once CppUnusedIncludeDirective
#include <vulkan/vulkan_xlib_xrandr.h>
#endif

#include "Config.h"

#include <vulkan/vk_icd.h>

#include <gsl/gsl>

#if !defined(_MSC_VER)
#define __debugbreak() __asm__("int3")

inline void strcpy_s(char* destination, const char* source)
{
	strcpy(destination, source);
}
#endif
#define TODO_ERROR() if (1) { __debugbreak(); abort(); } else (void)0
#define FATAL_ERROR() if (1) { __debugbreak(); abort(); } else (void)0

#if defined(_MSC_VER)
#define CP_DLL_EXPORT __declspec(dllexport)
#else
#define CP_DLL_EXPORT
#endif

#pragma warning(push)
#pragma warning(disable: 26490 26474 26408 26409)

struct DeviceMemory;

class Buffer;
class BufferView;
class CommandBuffer;
class CommandPool;
class CompiledShaderModule;
class ComputePipeline;
class DescriptorPool;
class DescriptorSet;
class DescriptorSetLayout;
class Device;
class Event;
class Fence;
class Framebuffer;
class GraphicsPipeline;
class Image;
class ImageView;
class Instance;
class PhysicalDevice;
class Pipeline;
class PipelineCache;
class PipelineLayout;
class QueryPool;
class Queue;
class RenderPass;
class Sampler;
class Semaphore;
class ShaderModule;

#if defined(VK_NV_ray_tracing)
class RayTracingPipeline;
#endif

template<typename LocalType>
struct VulkanTypeHelper;

#define VK_TYPE_HELPER(VulkanType, LocalType, NonDispatchable) template<> struct VulkanTypeHelper<LocalType> { static constexpr auto IsNonDispatchable = NonDispatchable; using Type = VulkanType; }

VK_TYPE_HELPER(VkInstance, Instance, false);

VK_TYPE_HELPER(VkPhysicalDevice, PhysicalDevice, false);

VK_TYPE_HELPER(VkDevice, Device, false);

VK_TYPE_HELPER(VkQueue, Queue, false);

VK_TYPE_HELPER(VkSemaphore, Semaphore, true);

VK_TYPE_HELPER(VkCommandBuffer, CommandBuffer, false);

VK_TYPE_HELPER(VkFence, Fence, true);

VK_TYPE_HELPER(VkDeviceMemory, DeviceMemory, true);

VK_TYPE_HELPER(VkBuffer, Buffer, true);

VK_TYPE_HELPER(VkImage, Image, true);

VK_TYPE_HELPER(VkEvent, Event, true);

VK_TYPE_HELPER(VkQueryPool, QueryPool, true);

VK_TYPE_HELPER(VkBufferView, BufferView, true);

VK_TYPE_HELPER(VkImageView, ImageView, true);

VK_TYPE_HELPER(VkShaderModule, ShaderModule, true);

VK_TYPE_HELPER(VkPipelineCache, PipelineCache, true);

VK_TYPE_HELPER(VkPipelineLayout, PipelineLayout, true);

VK_TYPE_HELPER(VkRenderPass, RenderPass, true);

VK_TYPE_HELPER(VkPipeline, Pipeline, true);

VK_TYPE_HELPER(VkPipeline, ComputePipeline, true);

VK_TYPE_HELPER(VkPipeline, GraphicsPipeline, true);

#if defined(VK_NV_ray_tracing)
VK_TYPE_HELPER(VkPipeline, RayTracingPipeline, true);
#endif

VK_TYPE_HELPER(VkDescriptorSetLayout, DescriptorSetLayout, true);

VK_TYPE_HELPER(VkSampler, Sampler, true);

VK_TYPE_HELPER(VkDescriptorPool, DescriptorPool, true);

VK_TYPE_HELPER(VkDescriptorSet, DescriptorSet, true);

VK_TYPE_HELPER(VkFramebuffer, Framebuffer, true);

VK_TYPE_HELPER(VkCommandPool, CommandPool, true);

#if defined(VK_VERSION_1_1)
class DescriptorUpdateTemplate;
class SamplerYcbcrConversion;

VK_TYPE_HELPER(VkSamplerYcbcrConversion, SamplerYcbcrConversion, true);

VK_TYPE_HELPER(VkDescriptorUpdateTemplate, DescriptorUpdateTemplate, true);
#endif

#if defined(VK_KHR_surface)
VK_TYPE_HELPER(VkSurfaceKHR, VkIcdSurfaceBase, true);

#if defined(VK_USE_PLATFORM_WIN32_KHR)
VK_TYPE_HELPER(VkSurfaceKHR, VkIcdSurfaceWin32, true);
#endif

#if defined(VK_USE_PLATFORM_XCB_KHR)
VK_TYPE_HELPER(VkSurfaceKHR, VkIcdSurfaceXcb, true);
#endif

#if defined(VK_USE_PLATFORM_XLIB_KHR)
VK_TYPE_HELPER(VkSurfaceKHR, VkIcdSurfaceXlib, true);
#endif
#endif

#if defined(VK_KHR_swapchain)
class Swapchain;
VK_TYPE_HELPER(VkSwapchainKHR, Swapchain, true);
#endif

#if defined(VK_KHR_display)
class VDisplay;
class DisplayMode;

VK_TYPE_HELPER(VkDisplayKHR, VDisplay, true);
VK_TYPE_HELPER(VkDisplayModeKHR, DisplayMode, true);
#endif

#if defined(VK_EXT_debug_report)
class DebugReportCallback;
VK_TYPE_HELPER(VkDebugReportCallbackEXT, DebugReportCallback, true);
#endif

#if defined(VK_NVX_device_generated_commands)
class IndirectCommandsLayout;
class ObjectTable;

VK_TYPE_HELPER(VkObjectTableNVX, ObjectTable, true);

VK_TYPE_HELPER(VkIndirectCommandsLayoutNVX, IndirectCommandsLayout, true);
#endif

#if defined(VK_EXT_debug_utils)
class DebugUtilsMessenger;
VK_TYPE_HELPER(VkDebugUtilsMessengerEXT, DebugUtilsMessenger, true);
#endif

#if defined(VK_EXT_validation_cache)
class ValidationCache;
VK_TYPE_HELPER(VkValidationCacheEXT, ValidationCache, true);
#endif

#if defined(VK_NV_ray_tracing)
class AccelerationStructure;
VK_TYPE_HELPER(VkAccelerationStructureNV, AccelerationStructure, true);
#endif

#if defined(VK_INTEL_performance_query)
class PerformanceConfiguration;
VK_TYPE_HELPER(VkPerformanceConfigurationINTEL, PerformanceConfiguration, true);
#endif

#undef VK_TYPE_HELPER

template<typename T, class... Types>
T* Allocate(const VkAllocationCallbacks* pAllocator, VkSystemAllocationScope allocationScope, Types&& ... args)
{
	constexpr auto size = VulkanTypeHelper<T>::IsNonDispatchable ? sizeof(T) : sizeof(T) + 16;
	const auto data = pAllocator
		                  ? static_cast<uint8_t*>(pAllocator->pfnAllocation(pAllocator->pUserData, size, 16, allocationScope))
		                  : static_cast<uint8_t*>(malloc(size));

	if (!data)
	{
		return nullptr;
	}

	if (VulkanTypeHelper<T>::IsNonDispatchable)
	{
		return new(data) T(std::forward<Types>(args)...);
	}

	*static_cast<uintptr_t*>(static_cast<void*>(data)) = ICD_LOADER_MAGIC;
	return new(data + 16) T(std::forward<Types>(args)...);
}

template<typename T>
void Free(T* value, const VkAllocationCallbacks* pAllocator) noexcept
{
	value->OnDelete(pAllocator);
	value->~T();
	const auto data = static_cast<uint8_t*>(static_cast<void*>(value)) - (VulkanTypeHelper<T>::IsNonDispatchable ? 0 : 16);

	if (pAllocator)
	{
		pAllocator->pfnFree(pAllocator->pUserData, data);
	}
	else
	{
		free(data);
	}
}

#if defined(VK_USE_PLATFORM_WIN32_KHR)
template<>
inline void Free(VkIcdSurfaceWin32* value, const VkAllocationCallbacks* pAllocator) noexcept
{
	if (pAllocator)
	{
		pAllocator->pfnFree(pAllocator->pUserData, value);
	}
	else
	{
		free(value);
	}
}
#endif

#if defined(VK_USE_PLATFORM_XCB_KHR)
template<>
inline void Free(VkIcdSurfaceXcb* value, const VkAllocationCallbacks* pAllocator) noexcept
{
	if (pAllocator)
	{
		pAllocator->pfnFree(pAllocator->pUserData, value);
	}
	else
	{
		free(value);
	}
}
#endif

#if defined(VK_USE_PLATFORM_XLIB_KHR)
template<>
inline void Free(VkIcdSurfaceXlib* value, const VkAllocationCallbacks* pAllocator) noexcept
{
	if (pAllocator)
	{
		pAllocator->pfnFree(pAllocator->pUserData, value);
	}
	else
	{
		free(value);
	}
}
#endif

template<bool IsNonDispatchable>
void WrapVulkan(void* local, uint64_t* vulkan) noexcept;

template<bool IsNonDispatchable>
void* UnwrapVulkan(uint64_t vulkanValue) noexcept;

template<>
inline void WrapVulkan<false>(void* local, uint64_t* vulkan) noexcept
{
	assert(local);
	assert(vulkan);
	const auto data = static_cast<uint8_t*>(local);
	*vulkan = reinterpret_cast<uint64_t>(data - 16);
}

template<>
inline void* UnwrapVulkan<false>(uint64_t vulkanValue) noexcept
{
	assert(vulkanValue);
	return reinterpret_cast<void*>(vulkanValue + 16);
}

template<>
inline void WrapVulkan<true>(void* local, uint64_t* vulkan) noexcept
{
	assert(local);
	assert(vulkan);
	*vulkan = reinterpret_cast<uint64_t>(local);
}

template<>
inline void* UnwrapVulkan<true>(uint64_t vulkanValue) noexcept
{
	assert(vulkanValue);
	return reinterpret_cast<void*>(vulkanValue);
}

template<typename LocalType>
void WrapVulkan(LocalType* local, typename VulkanTypeHelper<LocalType>::Type* vulkan) noexcept
{
	WrapVulkan<VulkanTypeHelper<LocalType>::IsNonDispatchable>(local, reinterpret_cast<uint64_t*>(vulkan));
}

template<typename LocalType>
LocalType* UnwrapVulkan(typename VulkanTypeHelper<LocalType>::Type vulkanValue) noexcept
{
	return reinterpret_cast<LocalType*>(UnwrapVulkan<VulkanTypeHelper<LocalType>::IsNonDispatchable>(reinterpret_cast<uint64_t>(vulkanValue)));
}
#pragma warning(pop)

constexpr auto PI = 3.14159265358979323846;

max_stars_repo_name: MatthewSmit/CPVulkan

max_issues_repo_name: MatthewSmit/CPVulkan

max_forks_repo_name: MatthewSmit/CPVulkan


text:
#include <float.h>

#include <gsl/gsl_blas.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_spmatrix.h>
#include <gsl/gsl_vector.h>
#include <osqp/osqp.h>

#include "qdm.h"

#define QDM_THETA_MIN_THRESHOLD 0.0001

int
qdm_theta_optimize(
    gsl_vector *result,
    gsl_vector *emperical_quantiles,
    gsl_matrix *ix
)
{
  int status = 0;

  gsl_matrix *p = gsl_matrix_alloc(ix->size2, ix->size2);
  gsl_matrix *a = gsl_matrix_alloc(ix->size2, ix->size2);
  gsl_vector *q = gsl_vector_alloc(ix->size2);
  gsl_vector *l = gsl_vector_alloc(ix->size2);
  gsl_vector *u = gsl_vector_alloc(ix->size2);

  // FIXME check c_malloc...
  OSQPData *data         = c_malloc(sizeof(OSQPData));
  OSQPSettings *settings = c_malloc(sizeof(OSQPSettings));
  OSQPWorkspace *work    = NULL;

  /* Ensure all data fields are initialized before use. */
  *data = (OSQPData) {
    .n = 0,
    .m = 0,
    .P = NULL,
    .A = NULL,
    .q = NULL,
    .l = NULL,
    .u = NULL,
  };

  status = qdm_matrix_tmm(ix, p);
  if (status != 0) {
    goto cleanup;
  }

  qdm_matrix_select_upper_triangle(p);

  gsl_matrix_set_identity(a);

  gsl_vector_set_all(l, QDM_THETA_MIN_THRESHOLD);
  gsl_vector_set_all(u, INFINITY);

  status = gsl_blas_dgemv(CblasTrans, -1.0, ix, emperical_quantiles, 0, q);
  if (status != 0) {
    goto cleanup;
  }

  /* Solve... */
  data->n = ix->size2;
  data->m = ix->size2;

  status = qdm_matrix_to_csc_matrix(&data->P, p);
  if (status != 0) {
    goto cleanup;
  }

  status = qdm_matrix_to_csc_matrix(&data->A, a);
  if (status != 0) {
    goto cleanup;
  }

  data->q = q->data;
  data->l = l->data;
  data->u = u->data;

  osqp_set_default_settings(settings);
  settings->verbose = 0;

  status = osqp_setup(&work, data, settings);
  if (status != 0) {
    goto cleanup;
  }

  status = osqp_solve(work);
  if (status != 0) {
    goto cleanup;
  }

  gsl_vector_view beta0 = gsl_vector_view_array(work->solution->x, ix->size2);
  status = gsl_vector_memcpy(result, &beta0.vector);
  if (status != 0) {
    goto cleanup;
  }

cleanup:
  if (data != NULL) {
    if (data->A != NULL) {
      free(data->A->p);
      free(data->A->i);
      free(data->A->x);
    }

    c_free(data->A);

    if (data->P != NULL) {
      free(data->P->p);
      free(data->P->i);
      free(data->P->x);
    }

    c_free(data->P);

    c_free(data);
  }

  c_free(settings);
  osqp_cleanup(work);

  gsl_vector_free(u);
  gsl_vector_free(l);
  gsl_vector_free(q);
  gsl_matrix_free(a);
  gsl_matrix_free(p);

  return status;
}

void
qdm_theta_matrix_constrain(
    gsl_matrix *theta,
    double min
)
{
  /* First row is constrained to be greater than zero. */
  for (size_t j = 0; j < theta->size2; j++) {
    if (gsl_matrix_get(theta, 0, j) < 0) {
      gsl_matrix_set(theta, 0, j, min);
    }
  }

  /* Remaining cells are set to zero if the sum of the column is negative.
   *
   * NOTE: The first column is not modified.
   */
  for (size_t i = 1; i < theta->size1; i++) {
    for (size_t j = 1; j < theta->size2; j++) {
      gsl_vector_view column = gsl_matrix_column(theta, j);
      double sum = qdm_vector_sum(&column.vector);

      if (sum < 0) {
        gsl_matrix_set(theta, i, j, 0);
      }
    }
  }
}

ext: c

max_stars_repo_name: calebcase/qdm

max_issues_repo_name: calebcase/qdm

max_forks_repo_name: calebcase/qdm


text:
#pragma once
#include <gsl/assert>

constexpr char deathstring[] = "Expected Death";
constexpr char failed_set_terminate_deathstring[] = ".*";

// This prevents a failed call to set_terminate from failing the test suite.
constexpr const char* GetExpectedDeathString(std::terminate_handler handle)
{
    return handle ? deathstring : failed_set_terminate_deathstring;
}

max_stars_repo_name: Batfing/GSL

max_issues_repo_name: Batfing/GSL

max_forks_repo_name: Batfing/GSL


text:
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include "fermi.h"

int main(int argc, char *argv[]) {
    int num_bands = 2;
    double t = 1.0;
    double E0 = 6.0;
    double deltaE = 14.0;
    void Efn(double k[3], gsl_vector *energies) {
        double cx = cos(2.0 * M_PI * k[0]);
        double cy = cos(2.0 * M_PI * k[1]);
        double cz = cos(2.0 * M_PI * k[2]);
        double tk = -2.0 * t * (cx + cy + cz);
        int i;
        double E0_band;
        for (i = 0; i < num_bands; i++) {
            E0_band = E0 + ((double)i) * deltaE;
            gsl_vector_set(energies, i, E0_band + tk);
        }
    }
    int na = 8;
    int nb = 8;
    int nc = 8;
    int G_order[3] = {0, 1, 2};
    int G_neg[3] = {1, 1, 1};
    bool use_cache = true;

    EnergyCache *Ecache = init_EnergyCache(na, nb, nc, num_bands, G_order, G_neg, Efn, use_cache);

    double num_electrons = 1.0;
    double E_Fermi = 0.0;
    int err = FindFermi(num_electrons, Ecache, &E_Fermi);
    printf("Got E_Fermi = %f\n", E_Fermi);
    if (err != CTETRA_BISECT_OK) {
        printf("Error = %d returned from FindFermi.\n", err);
        return err;
    }

    double expected = 13.0;
    if (E_Fermi < 12.0 || E_Fermi > 14.0) {
        printf("Incorrect E_Fermi; got %f, expected %f\n", E_Fermi, expected);
        return 1;
    }

    num_electrons = 0.5;
    err = FindFermi(num_electrons, Ecache, &E_Fermi);
    printf("Got E_Fermi = %f\n", E_Fermi);
    if (err != CTETRA_BISECT_OK) {
        printf("Error = %d returned from FindFermi.\n", err);
        return err;
    }

    double tol = 1e-8;
    expected = 6.0;
    if (fabs(E_Fermi - expected) > tol) {
        printf("Incorrect E_Fermi; got %f, expected %f\n", E_Fermi, expected);
        return 1;
    }

    printf("Fermi test passed.\n");
    return 0;
}

max_stars_repo_name: tflovorn/ctetra

max_issues_repo_name: tflovorn/ctetra

max_forks_repo_name: tflovorn/ctetra


text:
/**
 * @file common.h
 * @brief Public header file containing the common objects, API used by different applications
 */

#ifndef COMMON_H
#define COMMON_H

#include <petsc.h>

/**
 * @brief The communicator context
 */
struct _p_COMM {
  MPI_Comm     type; /**< MPI communicator SELF or WORLD */
  PetscMPIInt  rank; /**< Process rank */
  PetscMPIInt  size; /**< Communicator size */
  PetscInt     refct; /**< Reference count to know how many objects are sharing the communicator */
};

typedef struct _p_COMM *COMM;

/**
 * @brief Creates the communicator object COMM
 * @param [in] MPI_Comm mpicomm - The MPI communicator
 * @param [out] COMM* outcomm - The COMM object
 */
extern PetscErrorCode COMMCreate(MPI_Comm,COMM*);
/**
 * @brief Destroys the communicator object COMM created with COMMCreate
 * @param [in] COMM* outcomm - The COMM object
 */
extern PetscErrorCode COMMDestroy(COMM*);
/**
 * @brief NOT IMPLIMENTED
 */
extern PetscErrorCode SetMatrixValues(Mat,PetscInt,PetscInt[],PetscInt,PetscInt[],PetscScalar[]);
#endif

max_stars_repo_name: pflow-team/PFLOW

max_issues_repo_name: pflow-team/PFLOW

max_forks_repo_name: pflow-team/PFLOW


text:
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_math.h>

#define pi M_PI

//#############################################################################
//#############################################################################
//Ray path solver

//Data structure for pointers to julia cavity data and functions
typedef struct {
    //pointer to Boundary object
    void *bnd;
    //pointer to radius function
    double (*rfunc_p)(void *bnd,double theta);
    //pointer to (mutating) radius and normal vector angle function
    void (*rsys_p)(void *bnd, double theta, double results[]);
    //pointer to RefractiveIndex object
    void *idx;
    //pointer to refractive index value function
    double (*nfunc_p)(void *idx, double r, double theta);
    //pointer to (mutating) refractive index value and derivative function
    void (*nderiv_p)(void *idx, double r, double theta, double results[]);
} cavinfo;

//#############################################################################
//Functions for ODE

//ODE derivatives function
int odefunc(double t, const double y[], double f[], void *params){
    //Get Julia pointers
    cavinfo *cip = (cavinfo *)params;

    //Calculate n and its derivatives and store in results array
    double results[3];
    (*(*cip).nderiv_p)((*cip).idx,y[0],y[1],results);
    //results[0] = n, results[1] = dn/dr, results[2] = dn/dtheta

    //Debugging purposes
    //printf("n = %.5f, dr_n = %.5f, dtheta_n = %.5f\n",n,dr_n,dtheta_n);

    //Calculate and store derivative of ODE coordinate vector
    //Equations are:
    //r' = p_r/n
    //theta' = p_theta/(r^2*n)
    //p_r' = p_theta^2/(r^3*n) + dn/dr
    //p_theta' = dn/dtheta
    f[0] = y[2]/results[0];
    f[1] = y[3]/(y[0]*y[0]*results[0]);
    f[2] = y[3]*y[3]/(y[0]*y[0]*y[0]*results[0]) + results[1];
    f[3] = results[2];

    return GSL_SUCCESS;
}

//ODE Jacobian function
int odejac(double t, const double y[], double *dfdy, double dfdt[], void *params){
    //RK8PD does not require the Jacobian
    return GSL_SUCCESS;
}

//#############################################################################
//Solver functions

//Ray reflection
//Re-compute a ray coordinate vector to simulate a bounce, using the ODE coordinate vectors S0 and S immediately before and after the bounce, where S0 = (r0,theta0,pr0,ptheta0), and S = (r,theta,pr,ptheta).
double raybounce(cavinfo *cip, double S0[], double S[]){

    //Binary search for intersection of trajectory with cavity boundary
    //Assume that light trajectory between (r0,theta0) and (r,theta) is a straight line
    //(fastest way to interpolate, after all this is already within 1 stepsize of time)
    const double x0 = S0[0]*cos(S0[1]), y0 = S0[0]*sin(S0[1]);
    const double x1 = S[0]*cos(S[1]), y1 = S[0]*sin(S[1]);
    double uA = 0.0, uB = 1.0; //bounds of the binary search
    double xC,yC,rC,thetaC,RC,uC;
    do{
        uC = 0.5*(uA+uB); //get middle point
        xC = (1-uC)*x0 + uC*x1; yC = (1-uC)*y0 + uC*y1;
        rC = hypot(xC,yC); thetaC = atan2(yC,xC);
        RC = (*(*cip).rfunc_p)((*cip).bnd,thetaC);
        //Change boundary
        if(rC > RC) uB = uC;
        else uA = uC;
    } while(fabs(rC-RC) > 1e-12);

    //Get angle of incidence and store in results array
    double results[2];
    (*(*cip).rsys_p)((*cip).bnd,thetaC,results);
    //results[0] = rC, results[1] = alpha
    //Use linear interpolation to find ray angle phi at intersection
    const double phi0 = S0[1] + atan2(S0[3],S0[0]*S0[2]);
    const double phi1 = S[1] + atan2(S[3],S[0]*S[2]);
    const double phiC = (1-uC)*phi0 + uC*phi1;
    //Angle of incidence
    const double chi = phiC - results[1];

    //Store interpolated and reflected ODE coordinate vector
    const double phi = pi - chi + results[1]; //ray angle after reflection
    const double n = (*(*cip).nfunc_p)((*cip).idx,rC,thetaC);
    const double pr = n*cos(phi-thetaC);
    const double ptheta = n*rC*sin(phi-thetaC);
    S[0] = rC; S[1] = thetaC; S[2] = pr; S[3] = ptheta;

    //Report bounce information
    return chi;
}

//Ray evolution
void rayevolve(
    //Storage arrays
    double raypath_r[], double raypath_theta[], long bounceindices[],
    double bouncepts_chi[], long lengths[],
    //Initial conditions
    double r0, double theta0, double pr0, double ptheta0,
    //Simulation parameters
    double tmax, long bouncemax, double reltol, double abstol,
    //Cavity parameters
    void *bnd, double (*rfunc_p)(void *bnd,double theta),
    void (*rsys_p)(void *bnd, double theta,double results[]),
    void *idx, double (*nfunc_p)(void *idx, double r, double theta),
    void (*nderiv_p)(void *idx, double r, double theta,double results[])){

        //Condense input cavity data into cavinfo struct
        cavinfo ci = {bnd,rfunc_p,rsys_p,idx,nfunc_p,nderiv_p};

        //Initialize solver
        gsl_odeiv2_system sys = {odefunc,odejac,4,&ci};
        const gsl_odeiv2_step_type *steptype = gsl_odeiv2_step_rk8pd;
        gsl_odeiv2_step *step = gsl_odeiv2_step_alloc(steptype,4);
        gsl_odeiv2_control *control = gsl_odeiv2_control_y_new(abstol,reltol);
        gsl_odeiv2_evolve *evolve = gsl_odeiv2_evolve_alloc(4);

        //Initialize parameters
        double t = 0.0, dt = 0.0001;
        double y0[4], y[4] = {r0,theta0,pr0,ptheta0};

        //Prepare results record
        long bouncenum = 0, stepnum = 1; //indicates postion to record next
        const long prealloc = 250*ceil(tmax); //length of preallocated raypath array
        raypath_r[0] = y[0]; raypath_theta[0] = y[1];

        //Solver loop
        while(t < tmax && bouncenum < bouncemax && stepnum < prealloc){
            //Record initial position
            y0[0] = y[0], y0[1] = y[1], y0[2] = y[2], y0[3] = y[3];

            //Run Solver
            int status =
                gsl_odeiv2_evolve_apply(evolve,control,step,&sys,&t,tmax,&dt,y);
            if(status != GSL_SUCCESS) break;

            //Check Hamiltonian once in a while for sanity check
            if(stepnum%1000 == 0){
                double H = y[2]*y[2]+y[3]*y[3]/(y[0]*y[0]) - gsl_pow_2((*ci.nfunc_p)(ci.idx,y[0],y[1]));
                if(H > 1e-9) printf("Warning: Error in Hamiltonian is %.f\n",H);
            }

            //Check difference in ray and boundary radial positions
            double dr = y[0] - (*ci.rfunc_p)(ci.bnd,y[1]);
            if(dr > 0){
                //Boundary crossing!
                //get chi and corrected y
                bouncepts_chi[bouncenum] = raybounce(&ci,y0,y);
                //store (Julia's 1-based) index for thetaC values recorded in
                //raypath_theta array
                bounceindices[bouncenum] = stepnum+1;
                bouncenum += 1;
            }

            //Record position
            raypath_r[stepnum] = y[0]; raypath_theta[stepnum] = y[1];
            stepnum += 1;

            //Display progress for debugging
            //printf("t = %.5f, y = [%.5f,%.5f,%.5f,%.5f]\n",t,y[0],y[1],y[2],y[3]);

        }
        //Store lengths of arrays
        lengths[0] = stepnum; lengths[1] = bouncenum;

        //Free memory
        gsl_odeiv2_step_free(step);
        gsl_odeiv2_control_free(control);
        gsl_odeiv2_evolve_free(evolve);

}


text:
#include "../precision.h"

#undef SKIP

#ifdef GSL_FLOAT
#ifdef MATLAB_FLOAT_H
#define SKIP
#else
#define MATLAB_FLOAT_H
#endif
#endif

#ifdef GSL_DOUBLE
#ifdef MATLAB_H
#define SKIP
#else
#define MATLAB_H
#endif
#endif

#ifdef GSL_LONG_DOUBLE
#ifdef MATLAB_LONG_DOUBLE_H
#define SKIP
#else
#define MATLAB_LONG_DOUBLE_H
#endif
#endif

#ifndef SKIP

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_vector.h>
#include <limits>
#include <string>

namespace MATLAB_NAMESPACE {

	// Functions
	VECTOR_T* abs(const VECTOR_T* v);
	MATRIX_T* abs(const MATRIX_T* m);
	int all(const VECTOR_T* v);
	VECTOR_T* all(const MATRIX_T* m, int dim = 1);
	int any(const VECTOR_T* v);
	VECTOR_T* any(const MATRIX_T* m, int dim = 1);
	std::string dec2bin(int n);
	std::string dec2bin(int n, int len);
	MATRIX_T* diag(const VECTOR_T* v, int k = 0);
	VECTOR_T* diag(const MATRIX_T* m, int k = 0);
	MATRIX_T* eye(int size);
	MATRIX_T* eye(int size1, int size2);
	VECTOR_T* find(const VECTOR_T* v, int n = std::numeric_limits<int>::max(), const std::string& direction = "first");
	VECTOR_T* find(const MATRIX_T* m, int n = std::numeric_limits<int>::max(), const std::string& direction = "first");
	MATRIX_T* find_ij(const MATRIX_T* m, int n = std::numeric_limits<int>::max(), const std::string& direction = "first");
	VECTOR_T* hist(const VECTOR_T* v, int n = 10);
	VECTOR_T* hist(const VECTOR_T* v, const VECTOR_T* centers);
	MATRIX_T* inv(const MATRIX_T* m);
	int length(const VECTOR_T* v);
	int length(const MATRIX_T* m);
	FP_T max(FP_T x, FP_T y);
	FP_T max(const VECTOR_T* v);
	VECTOR_T* max(const MATRIX_T* m, int dim = 1);
	FP_T mean(const VECTOR_T* v, const std::string& opt = "a");
	VECTOR_T* mean(const MATRIX_T* m, int dim = 1, const std::string& opt = "a");
	FP_T min(FP_T x, FP_T y);
	FP_T min(const VECTOR_T* v);
	VECTOR_T* min(const MATRIX_T* m, int dim = 1);
	int nnz(const VECTOR_T* v);
	int nnz(const MATRIX_T* m);
	VECTOR_T* nonzeros(const MATRIX_T* m);
	FP_T norm(const VECTOR_T* v, int p);
	VECTOR_T* normpdf(const VECTOR_T* v, FP_T mean, FP_T stdev);
	MATRIX_T* ones(int size);
	MATRIX_T* ones(int size1, int size2);
	VECTOR_T* ones_vector(int size);
	FP_T prod(const VECTOR_T* v);
	VECTOR_T* prod(const MATRIX_T* m, int dim = 1);
	MATRIX_T* rand(int size);
	MATRIX_T* rand(int size1, int size2);
	VECTOR_T* rand_vector(int size);
	gsl_permutation* randperm(int size);
	VECTOR_T* reverse(const VECTOR_T* v);
	VECTOR_T* setxor(const VECTOR_T* v1, const VECTOR_T* v2);
	VECTOR_T* sort(const VECTOR_T* v, const std::string& mode = "ascend", VECTOR_T** ind = NULL);
	MATRIX_T* sort(const MATRIX_T* m, int dim = 1, const std::string& mode = "ascend", MATRIX_T** ind = NULL);
	VECTOR_T* sortrows(const VECTOR_T* v, VECTOR_T** ind = NULL);
	MATRIX_T* sortrows(const MATRIX_T* m, VECTOR_T** ind = NULL);
	FP_T std(const VECTOR_T* v, int opt = 0);
	VECTOR_T* std(const MATRIX_T* m, int opt = 0, int dim = 1);
	FP_T sum(const VECTOR_T* v);
	VECTOR_T* sum(const MATRIX_T* m, int dim = 1);
	MATRIX_T* toeplitz(const VECTOR_T* column, const VECTOR_T* row = NULL);
	MATRIX_T* tril(const MATRIX_T* m, int k = 0);
	MATRIX_T* triu(const MATRIX_T* m, int k = 0);
	VECTOR_T* unique(const VECTOR_T* v, const std::string& first_or_last = "last", VECTOR_T** i = NULL, VECTOR_T** j = NULL);
	VECTOR_T* unique(const MATRIX_T* m, const std::string& first_or_last = "last", VECTOR_T** i = NULL, VECTOR_T** j = NULL);
	MATRIX_T* unique_rows(const MATRIX_T* m, const std::string& first_or_last = "last", VECTOR_T** i = NULL, VECTOR_T** j = NULL);
	MATRIX_T* zeros(int size);
	MATRIX_T* zeros(int size1, int size2);
	VECTOR_T* zeros_vector(int size);

	// Operators
	VECTOR_T* concatenate(const VECTOR_T* v, FP_T x);
	VECTOR_T* concatenate(FP_T x, const VECTOR_T* v);
	VECTOR_T* concatenate(const VECTOR_T* v1, const VECTOR_T* v2);
	MATRIX_T* concatenate_columns(const VECTOR_T* v1, const VECTOR_T* v2);
	MATRIX_T* concatenate_columns(const MATRIX_T* m, const VECTOR_T* v);
	MATRIX_T* concatenate_columns(const VECTOR_T* v, const MATRIX_T* m);
	MATRIX_T* concatenate_columns(const MATRIX_T* m1, const MATRIX_T* m2);
	MATRIX_T* concatenate_rows(const VECTOR_T* v1, const VECTOR_T* v2);
	MATRIX_T* concatenate_rows(const MATRIX_T* m, const VECTOR_T* v);
	MATRIX_T* concatenate_rows(const VECTOR_T* v, const MATRIX_T* m);
	MATRIX_T* concatenate_rows(const MATRIX_T* m1, const MATRIX_T* m2);
	VECTOR_T* copy(const VECTOR_T* v);
	MATRIX_T* copy(const MATRIX_T* m);
	MATRIX_T* div_left(const MATRIX_T* m1, const MATRIX_T* m2);
	MATRIX_T* div_right(const MATRIX_T* m1, const MATRIX_T* m2);
	VECTOR_T* logical_and(const VECTOR_T* v1, const VECTOR_T* v2);
	MATRIX_T* logical_and(const MATRIX_T* m1, const MATRIX_T* m2);
	VECTOR_T* logical_not(const VECTOR_T* v);
	MATRIX_T* logical_not(const MATRIX_T* m);
	VECTOR_T* logical_or(const VECTOR_T* v1, const VECTOR_T* v2);
	MATRIX_T* logical_or(const MATRIX_T* m1, const MATRIX_T* m2);
	MATRIX_T* mul(const MATRIX_T* m1, const MATRIX_T* m2);
	MATRIX_T* pow(const MATRIX_T* m, int power);
	VECTOR_T* pow_elements(const VECTOR_T* v, FP_T power);
	VECTOR_T* pow_elements(const VECTOR_T* v, const VECTOR_T* powers);
	MATRIX_T* pow_elements(const MATRIX_T* m, FP_T power);
	MATRIX_T* pow_elements(const MATRIX_T* m, const MATRIX_T* powers);
	VECTOR_T* sequence(int start, int end);
	VECTOR_T* sequence(int start, int step, int end);

	// Floating-point comparison
	extern FP_T epsilon;
	int fp_compare(FP_T x, FP_T y);
	bool fp_zero(FP_T x);
	bool fp_nonzero(FP_T x);
	bool fp_equal(FP_T x, FP_T y);
	bool fp_not_equal(FP_T x, FP_T y);
	bool fp_less(FP_T x, FP_T y);
	bool fp_less_or_equal(FP_T x, FP_T y);
	bool fp_greater(FP_T x, FP_T y);
	bool fp_greater_or_equal(FP_T x, FP_T y);

	// Vector/matrix comparison
	typedef bool (*comparator)(FP_T, FP_T);
	int compare_vectors(const VECTOR_T* v1, const VECTOR_T* v2);
	bool vector_less(VECTOR_T* v1, VECTOR_T* v2);
	int compare_matrices(const MATRIX_T* m1, const MATRIX_T* m2);
	bool matrix_less(MATRIX_T* m1, MATRIX_T* m2);
	VECTOR_T* compare_elements(const VECTOR_T* v, comparator compare, FP_T x);
	VECTOR_T* compare_elements(const VECTOR_T* v1, comparator compare, const VECTOR_T* v2);
	MATRIX_T* compare_elements(const MATRIX_T* m, comparator compare, FP_T x);
	MATRIX_T* compare_elements(const MATRIX_T* m1, comparator compare, const MATRIX_T* m2);

	// Vector-by-vector indexing
	VECTOR_T* ordinal_index(const VECTOR_T* v, const VECTOR_T* indices);
	void ordinal_index_assign(VECTOR_T* v, const VECTOR_T* indices, FP_T value);
	void ordinal_index_assign(VECTOR_T* v, const VECTOR_T* indices, const VECTOR_T* values);
	VECTOR_T* logical_index(const VECTOR_T* v, const VECTOR_T* logical_v);
	void logical_index_assign(VECTOR_T* v, const VECTOR_T* logical_v, FP_T value);
	void logical_index_assign(VECTOR_T* v, const VECTOR_T* logical_v, const VECTOR_T* values);

	// Matrix-by-integer indexing
	FP_T ordinal_index(const MATRIX_T* m, int index);
	void ordinal_index_assign(MATRIX_T* m, int index, FP_T value);

	// Matrix-by-vector indexing
	VECTOR_T* ordinal_index(const MATRIX_T* m, const VECTOR_T* indices);
	void ordinal_index_assign(MATRIX_T* m, const VECTOR_T* indices, FP_T value);
	void ordinal_index_assign(MATRIX_T* m, const VECTOR_T* indices, const VECTOR_T* values);
	VECTOR_T* logical_index(const MATRIX_T* m, const VECTOR_T* logical_v);
	void logical_index_assign(MATRIX_T* m, const VECTOR_T* logical_v, FP_T value);
	void logical_index_assign(MATRIX_T* m, const VECTOR_T* logical_v, const VECTOR_T* values);

	// Matrix-by-two-vectors indexing (non-mixed)
	MATRIX_T* ordinal_index(const MATRIX_T* m, const VECTOR_T* rows, const VECTOR_T* columns);
	void ordinal_index_assign(MATRIX_T* m, const VECTOR_T* rows, const VECTOR_T* columns, FP_T value);
	void ordinal_index_assign(MATRIX_T* m, const VECTOR_T* rows, const VECTOR_T* columns, const MATRIX_T* values);
	MATRIX_T* logical_index(const MATRIX_T* m, const VECTOR_T* logical_rows, const VECTOR_T* logical_columns);
	void logical_index_assign(MATRIX_T* m, const VECTOR_T* logical_rows, const VECTOR_T* logical_columns, FP_T value);
	void logical_index_assign(MATRIX_T* m, const VECTOR_T* logical_rows, const VECTOR_T* logical_columns, const MATRIX_T* values);

	// Matrix-by-two-vectors indexing (mixed)
	MATRIX_T* ord_log_index(const MATRIX_T* m, const VECTOR_T* rows, const VECTOR_T* logical_columns);
	void ord_log_index_assign(MATRIX_T* m, const VECTOR_T* rows, const VECTOR_T* logical_columns, FP_T value);
	void ord_log_index_assign(MATRIX_T* m, const VECTOR_T* rows, const VECTOR_T* logical_columns, const MATRIX_T* values);
	MATRIX_T* log_ord_index(const MATRIX_T* m, const VECTOR_T* logical_rows, const VECTOR_T* columns);
	void log_ord_index_assign(MATRIX_T* m, const VECTOR_T* logical_rows, const VECTOR_T* columns, FP_T value);
	void log_ord_index_assign(MATRIX_T* m, const VECTOR_T* logical_rows, const VECTOR_T* columns, const MATRIX_T* values);

	// Matrix-by-matrix indexing
	MATRIX_T* ordinal_index(const MATRIX_T* m, const MATRIX_T* indices);
	void ordinal_index_assign(MATRIX_T* m, const MATRIX_T* indices, FP_T value);
	void ordinal_index_assign(MATRIX_T* m, const MATRIX_T* indices, const MATRIX_T* values);
	VECTOR_T* logical_index(const MATRIX_T* m, const MATRIX_T* logical_m);
	void logical_index_assign(MATRIX_T* m, const MATRIX_T* logical_m, FP_T value);
	void logical_index_assign(MATRIX_T* m, const MATRIX_T* logical_m, const VECTOR_T* values);

	// Vector/matrix conversion
	void to_array(const VECTOR_T* v, FP_T* array);
	bool to_bool(const VECTOR_T* v);
	bool to_bool(const MATRIX_T* m);
	VECTOR_T* to_vector(const gsl_vector* v_d);
	gsl_vector* to_vector_double(const VECTOR_T* v);
	VECTOR_T* to_vector(const MATRIX_T* m);
	MATRIX_T* to_column_matrix(const VECTOR_T* v);
	MATRIX_T* to_row_matrix(const VECTOR_T* v);
	MATRIX_T* to_matrix(const gsl_matrix* m_d);
	gsl_matrix* to_matrix_double(const MATRIX_T* m);
	VECTOR_T* to_vector(const gsl_permutation* p);
	gsl_permutation* to_permutation(const VECTOR_T* v);

	// Utility
	gsl_rng* get_rng();
	void seed_rng(const gsl_rng* rng, unsigned long seed);
	VECTOR_T* permute(const gsl_permutation* p, const VECTOR_T* v);
	MATRIX_T* permute_columns(const gsl_permutation* p, const MATRIX_T* m);
	MATRIX_T* permute_rows(const gsl_permutation* p, const MATRIX_T* m);
}

#endif

max_stars_repo_name: devuci/bct-cpp

max_issues_repo_name: devuci/bct-cpp

max_forks_repo_name: devuci/bct-cpp


text:
// The MIT License (MIT)
//
// Copyright (c) 2018 Mateusz Pusz
//
// 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

#include <units/bits/external/hacks.h>
#include <units/bits/math_concepts.h>
#include <units/bits/ratio_maths.h>
#include <units/bits/pow.h>
#include <units/bits/root.h>
#include <cstdint>
#include <numeric>
#include <type_traits>
#include <tuple>
#include <gsl/gsl-lite.hpp>

namespace units {

struct ratio;
constexpr ratio inverse(const ratio& r);

/**
 * @brief Provides compile-time rational arithmetic support.
 *
 * This class is really similar to @c std::ratio but gets an additional `Exp`
 * template parameter that defines the exponent of the ratio. Another important
 * difference is the fact that the objects of that class are used as class NTTPs
 * rather then a type template parameter kind.
 */
struct ratio {
  std::intmax_t num;
  std::intmax_t den;
  std::intmax_t exp;

  constexpr explicit ratio(std::intmax_t n, std::intmax_t d = 1, std::intmax_t e = 0): num(n), den(d), exp(e)
  {
    gsl_Expects(den != 0);
    detail::normalize(num, den, exp);
  }

  [[nodiscard]] friend constexpr bool operator==(const ratio&, const ratio&) = default;

  [[nodiscard]] friend constexpr ratio operator*(const ratio& lhs, const ratio& rhs)
  {
    const std::intmax_t gcd1 = std::gcd(lhs.num, rhs.den);
    const std::intmax_t gcd2 = std::gcd(rhs.num, lhs.den);
    return ratio(detail::safe_multiply(lhs.num / gcd1, rhs.num / gcd2),
                 detail::safe_multiply(lhs.den / gcd2, rhs.den / gcd1),
                 lhs.exp + rhs.exp);
  }

  [[nodiscard]] friend constexpr ratio operator/(const ratio& lhs, const ratio& rhs)
  {
    return lhs * inverse(rhs);
  }
};

[[nodiscard]] constexpr ratio inverse(const ratio& r)
{
  return ratio(r.den, r.num, -r.exp);
}

[[nodiscard]] constexpr bool is_integral(const ratio& r)
{
  if(r.exp < 0) {
    return false;
  } else {
    return detail::gcdpow(r.num, r.exp, r.den) == r.den;
  }
}

namespace detail {

[[nodiscard]] constexpr auto make_exp_align(const ratio& r, std::intmax_t alignment)
{
  gsl_Expects(alignment > 0);
  const std::intmax_t rem = r.exp % alignment;

  if (rem == 0) {  // already aligned
    return std::array{r.num, r.den, r.exp};
  }

  if (r.exp > 0) {  // remainder is positive
    return std::array{r.num * ipow10(rem), r.den, r.exp - rem};
  }

  // remainder is negative
  return std::array{r.num, r.den * ipow10(-rem), r.exp - rem};
}

template<std::intmax_t N>
  requires gt_zero<N>
[[nodiscard]] constexpr ratio root(const ratio& r)
{
  if constexpr (N == 1) {
    return r;
  } else {
    if (r.num == 0) {
      return ratio(0);
    }

    const auto aligned = make_exp_align(r, N);
    return ratio(iroot<N>(aligned[0]), iroot<N>(aligned[1]), aligned[2] / N);
  }
}

}  // namespace detail

template<std::intmax_t Num, std::intmax_t Den = 1>
  requires detail::non_zero<Den>
[[nodiscard]] constexpr ratio pow(const ratio& r)
{
  if constexpr (Num == 0) {
    return ratio(1);
  } else if constexpr (Num == Den) {
    return r;
  } else {
    // simplify factors first and compute power for positive exponent
    constexpr std::intmax_t gcd = std::gcd(Num, Den);
    constexpr std::intmax_t num = detail::abs(Num / gcd);
    constexpr std::intmax_t den = detail::abs(Den / gcd);

    // integer root loses precision so do pow first
    const ratio result = detail::root<den>(detail::pow_impl<num>(r));

    if constexpr (Num * Den < 0) {  // account for negative exponent
      return inverse(result);
    } else {
      return result;
    }
  }
}

[[nodiscard]] constexpr ratio sqrt(const ratio& r) { return pow<1, 2>(r); }

[[nodiscard]] constexpr ratio cbrt(const ratio& r) { return pow<1, 3>(r); }

// common_ratio
[[nodiscard]] constexpr ratio common_ratio(const ratio& r1, const ratio& r2)
{
  const auto res = detail::gcd_frac(r1.num, r1.den, r1.exp, r2.num, r2.den, r2.exp);
  return ratio(res[0], res[1], res[2]);
}

}  // namespace units

max_stars_repo_name: HazardyKnusperkeks/units

max_issues_repo_name: HazardyKnusperkeks/units

max_forks_repo_name: HazardyKnusperkeks/units


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "fitsio.h"
#include <gsl/gsl_fit.h>
#include <gsl/gsl_multifit.h>

//The Function Declarations
int read_image(char[],float[],long[]);
int read_table(char[]);
int write_image(char[],float[],long[]);
int image_size(char[], long[]);
int fit(double[],double[],double[],float[],int);
int calc_pr(float[],float[],float,int,int,long[]);
void printerror(int);

max_stars_repo_name: boada/ICD

max_issues_repo_name: boada/ICD

max_forks_repo_name: boada/ICD


text:
//Gets deltas (1st order) and delta-deltas (2nd order) differences of X.
//Only the deltas and delta-deltas are output in Y.
//Thus, Y must be pre-allocated to have twice the size of X.

//I implement this just like FIR for speed and shorter code,
//except non-causal and mid-sample of B is 0,
//so I don't explicitly make B (e.g., B[n] just equals sc*n).

//Note that this may treat edge samples differently than other code.
//But this could be changed with a few lines of additional code here,
//without changing the super-efficient FIR implementation.
//That is, since out-of-range samps were assumed to be 0, nothing was
//added to Y for them, so can just add something later.

#include <stdio.h>
#include <cblas.h>

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

int get_delta_deltas_s (float *Y, const float *X, const int iscolmajor, const int R, const int C, const int dim, const int N);
int get_delta_deltas_d (double *Y, const double *X, const int iscolmajor, const int R, const int C, const int dim, const int N);

int get_delta_deltas_s (float *Y, const float *X, const int iscolmajor, const int R, const int C, const int dim, const int N)
{
    const float z = 0.0f;
    int r, c, n;
    float sc = 1.0f;

    //Checks
    if (R<1) { fprintf(stderr,"error in get_delta_deltas_s: R (nrows Y) must be positive\n"); return 1; }
    if (C<1) { fprintf(stderr,"error in get_delta_deltas_s: C (ncols Y) must be positive\n"); return 1; }
    if (N<1) { fprintf(stderr,"error in get_delta_deltas_s: N (delta winlength) must be positive\n"); return 1; }

    //Get sc (normalizer)
    for (n=2; n<=N; n++) { sc += n*n; }
    sc = 0.5f/sc;

    //Initialize Y
    cblas_scopy(2*R*C,&z,0,&Y[0],1);

    if (dim==0)
    {
        if (iscolmajor)
        {
            for (c=0; c<C; c++)
            {
                for (n=1; n<=N; n++)
                {
                    cblas_saxpy(n,-sc*n,&X[c*R],0,&Y[c*R],1);            //beg edge samps
                    cblas_saxpy(R-n,-sc*n,&X[c*R],1,&Y[c*R+n],1);        //past samps
                    cblas_saxpy(R-n,sc*n,&X[c*R+n],1,&Y[c*R],1);         //future samps
                    cblas_saxpy(n,sc*n,&X[c*R+R-1],0,&Y[c*R+R-n],1);     //end edge samps
                }
                for (n=1; n<=N; n++)
                {
                    cblas_saxpy(n,-sc*n,&Y[c*R],0,&Y[R*C+c*R],1);        //beg edge samps
                    cblas_saxpy(R-n,-sc*n,&Y[c*R],1,&Y[R*C+c*R+n],1);    //past samps
                    cblas_saxpy(R-n,sc*n,&Y[c*R+n],1,&Y[R*C+c*R],1);     //future samps
                    cblas_saxpy(n,sc*n,&Y[c*R+R-1],0,&Y[R*C+c*R+R-n],1); //end edge samps
                }
            }
        }
        else
        {
            for (c=0; c<C; c++)
            {
                for (n=1; n<=N; n++)
                {
                    cblas_saxpy(n,-sc*n,&X[c],0,&Y[c],2*C);                      //beg edge samps
                    cblas_saxpy(R-n,-sc*n,&X[c],C,&Y[c+n*2*C],2*C);              //past samps
                    cblas_saxpy(R-n,sc*n,&X[c+n*C],C,&Y[c],2*C);                 //future samps
                    cblas_saxpy(n,sc*n,&X[c+C*(R-1)],0,&Y[c+2*C*(R-n)],2*C);     //end edge samps
                }
                for (n=1; n<=N; n++)
                {
                    cblas_saxpy(n,-sc*n,&Y[c],0,&Y[C+c],2*C);                    //beg edge samps
                    cblas_saxpy(R-n,-sc*n,&Y[c],2*C,&Y[C+c+n*2*C],2*C);          //past samps
                    cblas_saxpy(R-n,sc*n,&Y[c+n*2*C],2*C,&Y[C+c],2*C);           //future samps
                    cblas_saxpy(n,sc*n,&Y[c+2*C*(R-1)],0,&Y[C+c+2*C*(R-n)],2*C); //end edge samps
                }
            }
        }
    }
    else if (dim==1)
    {
        if (iscolmajor)
        {
            for (r=0; r<R; r++)
            {
                for (n=1; n<=N; n++)
                {
                    cblas_saxpy(n,-sc*n,&X[r],0,&Y[r],2*R);                      //beg edge samps
                    cblas_saxpy(C-n,-sc*n,&X[r],R,&Y[r+n*2*R],2*R);              //past samps
                    cblas_saxpy(C-n,sc*n,&X[r+n*R],R,&Y[r],2*R);                 //future samps
                    cblas_saxpy(n,sc*n,&X[r+R*(C-1)],0,&Y[r+2*R*(C-n)],2*R);     //end edge samps
                }
                for (n=1; n<=N; n++)
                {
                    cblas_saxpy(n,-sc*n,&Y[r],0,&Y[R+r],2*R);                    //beg edge samps
                    cblas_saxpy(C-n,-sc*n,&Y[r],2*R,&Y[R+r+n*2*R],2*R);          //past samps
                    cblas_saxpy(C-n,sc*n,&Y[r+n*2*R],2*R,&Y[R+r],2*R);           //future samps
                    cblas_saxpy(n,sc*n,&Y[r+2*R*(C-1)],0,&Y[R+r+2*R*(C-n)],2*R); //end edge samps
                }
            }
        }
        else
        {
            for (r=0; r<R; r++)
            {
                for (n=1; n<=N; n++)
                {
                    cblas_saxpy(n,-sc*n,&X[r*C],0,&Y[r*C],1);            //beg edge samps
                    cblas_saxpy(C-n,-sc*n,&X[r*C],1,&Y[r*C+n],1);        //past samps
                    cblas_saxpy(C-n,sc*n,&X[r*C+n],1,&Y[r*C],1);         //future samps
                    cblas_saxpy(n,sc*n,&X[r*C+C-1],0,&Y[r*C+C-n],1);     //end edge samps
                }
                for (n=1; n<=N; n++)
                {
                    cblas_saxpy(n,-sc*n,&Y[r*C],0,&Y[R*C+r*C],1);        //beg edge samps
                    cblas_saxpy(C-n,-sc*n,&Y[r*C],1,&Y[R*C+r*C+n],1);    //past samps
                    cblas_saxpy(C-n,sc*n,&Y[r*C+n],1,&Y[R*C+r*C],1);     //future samps
                    cblas_saxpy(n,sc*n,&Y[r*C+C-1],0,&Y[R*C+r*C+C-n],1); //end edge samps
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in get_delta_deltas_s: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

int get_delta_deltas_d (double *Y, const double *X, const int iscolmajor, const int R, const int C, const int dim, const int N)
{
    const double z = 0.0;
    int r, c, n;
    double sc = 1.0;

    //Checks
    if (R<1) { fprintf(stderr,"error in get_delta_deltas_d: R (nrows Y) must be positive\n"); return 1; }
    if (C<1) { fprintf(stderr,"error in get_delta_deltas_d: C (ncols Y) must be positive\n"); return 1; }
    if (N<1) { fprintf(stderr,"error in get_delta_deltas_d: N (delta winlength) must be positive\n"); return 1; }

    //Get sc (normalizer)
    for (n=2; n<=N; n++) { sc += n*n; }
    sc = 0.5/sc;

    //Initialize Y
    cblas_dcopy(2*R*C,&z,0,&Y[0],1);

    if (dim==0)
    {
        if (iscolmajor)
        {
            for (c=0; c<C; c++)
            {
                for (n=1; n<=N; n++)
                {
                    cblas_daxpy(n,-sc*n,&X[c*R],0,&Y[c*R],1);            //beg edge samps
                    cblas_daxpy(R-n,-sc*n,&X[c*R],1,&Y[c*R+n],1);        //past samps
                    cblas_daxpy(R-n,sc*n,&X[c*R+n],1,&Y[c*R],1);         //future samps
                    cblas_daxpy(n,sc*n,&X[c*R+R-1],0,&Y[c*R+R-n],1);     //end edge samps
                }
                for (n=1; n<=N; n++)
                {
                    cblas_daxpy(n,-sc*n,&Y[c*R],0,&Y[R*C+c*R],1);        //beg edge samps
                    cblas_daxpy(R-n,-sc*n,&Y[c*R],1,&Y[R*C+c*R+n],1);    //past samps
                    cblas_daxpy(R-n,sc*n,&Y[c*R+n],1,&Y[R*C+c*R],1);     //future samps
                    cblas_daxpy(n,sc*n,&Y[c*R+R-1],0,&Y[R*C+c*R+R-n],1); //end edge samps
                }
            }
        }
        else
        {
            for (c=0; c<C; c++)
            {
                for (n=1; n<=N; n++)
                {
                    cblas_daxpy(n,-sc*n,&X[c],0,&Y[c],2*C);                      //beg edge samps
                    cblas_daxpy(R-n,-sc*n,&X[c],C,&Y[c+n*2*C],2*C);              //past samps
                    cblas_daxpy(R-n,sc*n,&X[c+n*C],C,&Y[c],2*C);                 //future samps
                    cblas_daxpy(n,sc*n,&X[c+C*(R-1)],0,&Y[c+2*C*(R-n)],2*C);     //end edge samps
                }
                for (n=1; n<=N; n++)
                {
                    cblas_daxpy(n,-sc*n,&Y[c],0,&Y[C+c],2*C);                    //beg edge samps
                    cblas_daxpy(R-n,-sc*n,&Y[c],2*C,&Y[C+c+n*2*C],2*C);          //past samps
                    cblas_daxpy(R-n,sc*n,&Y[c+n*2*C],2*C,&Y[C+c],2*C);           //future samps
                    cblas_daxpy(n,sc*n,&Y[c+2*C*(R-1)],0,&Y[C+c+2*C*(R-n)],2*C); //end edge samps
                }
            }
        }
    }
    else if (dim==1)
    {
        if (iscolmajor)
        {
            for (r=0; r<R; r++)
            {
                for (n=1; n<=N; n++)
                {
                    cblas_daxpy(n,-sc*n,&X[r],0,&Y[r],2*R);                      //beg edge samps
                    cblas_daxpy(C-n,-sc*n,&X[r],R,&Y[r+n*2*R],2*R);              //past samps
                    cblas_daxpy(C-n,sc*n,&X[r+n*R],R,&Y[r],2*R);                 //future samps
                    cblas_daxpy(n,sc*n,&X[r+R*(C-1)],0,&Y[r+2*R*(C-n)],2*R);     //end edge samps
                }
                for (n=1; n<=N; n++)
                {
                    cblas_daxpy(n,-sc*n,&Y[r],0,&Y[R+r],2*R);                    //beg edge samps
                    cblas_daxpy(C-n,-sc*n,&Y[r],2*R,&Y[R+r+n*2*R],2*R);          //past samps
                    cblas_daxpy(C-n,sc*n,&Y[r+n*2*R],2*R,&Y[R+r],2*R);           //future samps
                    cblas_daxpy(n,sc*n,&Y[r+2*R*(C-1)],0,&Y[R+r+2*R*(C-n)],2*R); //end edge samps
                }
            }
        }
        else
        {
            for (r=0; r<R; r++)
            {
                for (n=1; n<=N; n++)
                {
                    cblas_daxpy(n,-sc*n,&X[r*C],0,&Y[r*C],1);            //beg edge samps
                    cblas_daxpy(C-n,-sc*n,&X[r*C],1,&Y[r*C+n],1);        //past samps
                    cblas_daxpy(C-n,sc*n,&X[r*C+n],1,&Y[r*C],1);         //future samps
                    cblas_daxpy(n,sc*n,&X[r*C+C-1],0,&Y[r*C+C-n],1);     //end edge samps
                }
                for (n=1; n<=N; n++)
                {
                    cblas_daxpy(n,-sc*n,&Y[r*C],0,&Y[R*C+r*C],1);        //beg edge samps
                    cblas_daxpy(C-n,-sc*n,&Y[r*C],1,&Y[R*C+r*C+n],1);    //past samps
                    cblas_daxpy(C-n,sc*n,&Y[r*C+n],1,&Y[R*C+r*C],1);     //future samps
                    cblas_daxpy(n,sc*n,&Y[r*C+C-1],0,&Y[R*C+r*C+C-n],1); //end edge samps
                }
            }
        }
    }
    else
    {
        fprintf(stderr,"error in get_delta_deltas_d: dim must be 0 or 1.\n"); return 1;
    }

    return 0;
}

#ifdef __cplusplus
}
}
#endif


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

#include <gsl/gsl_errno.h>

#include "ccl.h"

// Debug mode policy: whether to print error messages as they are raised.
// Defualt is ON.
static CCLDebugModePolicy _ccl_debug_mode_policy = CCL_DEBUG_MODE_ON;

// Set debug mode policy
void ccl_set_debug_policy(CCLDebugModePolicy debug_policy) {
    _ccl_debug_mode_policy = debug_policy;
}

// Convenience function to handle warnings
void ccl_raise_warning(int err, const char* msg, ...) {
  char message[256];

  va_list va;
  va_start(va, msg);
  vsnprintf(message, 250, msg, va);
  va_end(va);

  // For now just print warning to stderr if debug is enabled.
  // TODO: Implement some kind of error stack that can be passed on to, e.g.,
  // the python binding.
  if (_ccl_debug_mode_policy == CCL_DEBUG_MODE_ON) {
    fprintf(stderr, "WARNING %d: %s\n", err, message);
  }
}

// Convenience function to handle warnings
void ccl_raise_gsl_warning(int gslstatus, const char* msg, ...) {
  char message[256];

  va_list va;
  va_start(va, msg);
  vsnprintf(message, 250, msg, va);
  va_end(va);

  ccl_raise_warning(gslstatus, "%s: GSL ERROR: %s", message, gsl_strerror(gslstatus));
  return;
}

max_stars_repo_name: Jappenn/CCL

max_issues_repo_name: Jappenn/CCL

max_forks_repo_name: Jappenn/CCL


text:
#include "asf.h"
#include "asf_meta.h"
#include "dateUtil.h"

#include <assert.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>

static FILE *fopen_workreport_ext(const char *fileName, report_level_t level)
{
  char *path = getPath(fileName);
  char *basename = get_basename(fileName);
  char *txtFile = MALLOC(sizeof(char) * (strlen(path)+strlen(fileName)+20));

  char *p;
  if (strncmp_case(basename, "LED-", 4)==0)
    p = basename + 4;
  else
    p = basename;

  if (strlen(path) > 0)
    sprintf(txtFile, "%s/%s", path, p);
  else
    strcpy(txtFile, p);

  char *workreport_filename = appendExt(txtFile, ".txt");

  // first attempt: basename.txt
  FILE *fp;
  if (!fileExists(workreport_filename)) {
    // second attempt: path/'workreport'
    FREE(workreport_filename);
    workreport_filename = MALLOC(sizeof(char) * (strlen(path) + 20));

    if (strlen(path) > 0)
      sprintf(workreport_filename, "%s/workreport", path);
    else
      strcpy(workreport_filename, "workreport");

    if (!fileExists(workreport_filename)) {

      // third attempt: path/'summary.txt'
      if (strlen(path) > 0)
        sprintf(workreport_filename, "%s/summary.txt", path);
      else
        strcpy(workreport_filename, "summary.txt");

      if (!fileExists(workreport_filename)) {
        // failed!
        fp = NULL;
      }
      else {
        // success with 'summary.txt'
        asfReport(level, "workreport file found as: summary.txt\n");
        fp = FOPEN(workreport_filename, "r");
      }
    }
    else {
      // success with 'workreport'
      asfReport(level, "workreport file found as: workreport\n");
      fp = FOPEN(workreport_filename, "r");
    }
  }
  else {
    // success with 'basename.txt'
    char *basename = get_basename(workreport_filename);
    asfReport(level, "workreport file found as: %s\n", basename);
    fp = FOPEN(workreport_filename, "r");
    FREE(basename);
  }

  FREE(workreport_filename);
  FREE(txtFile);
  FREE(path);
  FREE(basename);

  return fp;
}

FILE *fopen_workreport(const char *fileName)
{
  return fopen_workreport_ext(fileName, REPORT_LEVEL_STATUS);
}

// Get the delta image time for ALOS data out of the summary file
int get_alos_delta_time (const char *fileName, double *delta)
{
  struct dataset_sum_rec dssr;
  hms_time dssr_time, summary_time, start_time, end_time;
  ymd_date dssr_date, summary_date, start_date, end_date;
  char line[512], dateStr[30], *str;

  get_dssr(fileName, &dssr);
  date_dssr2date(dssr.inp_sctim, &dssr_date, &dssr_time);

  FILE *fp = fopen_workreport(fileName);
  if (!fp) {
    // no workreport file...
    *delta = 0;
    return FALSE;
  }

  while (fgets(line, 512, fp)) {
    if (strstr(line, "Img_SceneCenterDateTime")) {
      str = strchr(line, '"');
      sprintf(dateStr, "%s", str+1);
      dateStr[strlen(dateStr)-2] = '\0';
      date_alos2date(dateStr, &summary_date, &summary_time);
      // bumped up the tolerance to 2 seconds... can't see how this would
      // introduce any false positives, and definitely reduces false negatives
      if (date_difference(&dssr_date, &dssr_time,
          &summary_date, &summary_time) > 2.0)
      {
        asfPrintWarning("Summary file does not correspond to leader file.\n"
                        "DSSR: %s\nSummary: %s\n",
                        dssr.inp_sctim, dateStr);
        *delta = 0;
        FCLOSE(fp);
        return FALSE;
      }
    }
    else if (strstr(line, "Img_SceneStartDateTime")) {
      str = strchr(line, '"');
      sprintf(dateStr, "%s", str+1);
      dateStr[strlen(dateStr)-2] = '\0';
      date_alos2date(dateStr, &start_date, &start_time);
    }
    else if (strstr(line, "Img_SceneEndDateTime")) {
      str = strchr(line, '"');
      sprintf(dateStr, "%s", str+1);
      dateStr[strlen(dateStr)-2] = '\0';
      date_alos2date(dateStr, &end_date, &end_time);
    }
  }

  *delta = date_difference(&start_date, &start_time, &end_date, &end_time);
  FCLOSE(fp);
  return TRUE;
}

// Get the SAR processor version from the workreport
int get_alos_processor_version (const char *fileName, double *version)
{
  char line[512], versionStr[10], *str;

  FILE *fp = fopen_workreport_ext(fileName, REPORT_LEVEL_NONE);
  if (!fp) {
    // no workreport file...
    *version = 0.0;
    return FALSE;
  }

  while (fgets(line, 512, fp)) {
    if (strstr(line, "Ver_PSR_CorPrcSigmaSAR")) {
      str = strchr(line, '"');
      sprintf(versionStr, "%s", str+1);
      versionStr[strlen(versionStr)-2] = '\0';
      *version = atof(versionStr);
    }
  }

  FCLOSE(fp);
  return TRUE;
}

// ------------------------------------------------------------------------
// helper code for refine_slc_geolocation_from_workreport()

struct refine_shift_params {
    meta_parameters *meta;
    double center_lat, center_lon;
    double ul_lat, ul_lon;
    double ur_lat, ur_lon;
    double ll_lat, ll_lon;
    double lr_lat, lr_lon;
};

static double lldist(double lat1, double lat2, double lon1, double lon2)
{
    double dlat = lat1 - lat2;
    double dlon;

    // some kludgery to handle crossing the meridian
    // get lat1,lon1 to be on the same side as lat2,lon2
    if (fabs(lon1-lon2) > 300) {
        if (lon2 < 0 && lon1 > 0) lon1 -= 360;
        if (lon2 > 0 && lon1 < 0) lon1 += 360;
    }

    dlon = lon1 - lon2;

    // Scale longitude difference to take into accound the fact
    // that longitude lines are a lot closer at the pole.
    dlon *= cos (lat2 * PI / 180.0);

    return dlat * dlat + dlon * dlon;
}

static double err_at_pixel(struct refine_shift_params *p, double off_t,
                           double off_x, double line, double samp,
                           double real_lat, double real_lon)
{
    const double huge_err = 999999.9;

    meta_parameters *meta = p->meta;
    double lat, lon;
    int bad=0;

    bad = meta_get_latLon(meta, line, samp, 0.0, &lat, &lon);
    if (bad) return huge_err;

    return lldist(lat, real_lat, lon, real_lon);
}

static int
getObjective(const gsl_vector *x, void *params, gsl_vector *f)
{
    double dt = gsl_vector_get(x,0);
    double ds = gsl_vector_get(x,1);

    if (!meta_is_valid_double(ds) || !meta_is_valid_double(dt)) {
        // This does happen sometimes, when we've already found the root
	return GSL_FAILURE;
    }

    struct refine_shift_params *p = (struct refine_shift_params *)params;
    meta_parameters *meta = p->meta;

    int nl = meta->general->line_count;
    int ns = meta->general->sample_count;

    double saved_timeOffset = meta->sar->time_shift;
    double saved_slant = meta->sar->slant_shift;
    meta->sar->time_shift += dt;
    meta->sar->slant_shift += ds;

    double err = 0.0;

    // add up the errors for the 5 known points
    err += err_at_pixel(p, dt, ds, nl/2, ns/2, p->center_lat, p->center_lon);

    if (meta->general->orbit_direction=='A') {
      err += err_at_pixel(p, dt, ds, 0,    0,    p->ll_lat, p->ll_lon);
      err += err_at_pixel(p, dt, ds, nl-1, 0,    p->ul_lat, p->ul_lon);
      err += err_at_pixel(p, dt, ds, 0,    ns-1, p->lr_lat, p->lr_lon);
      err += err_at_pixel(p, dt, ds, nl-1, ns-1, p->ur_lat, p->ur_lon);
    }
    else {
      err += err_at_pixel(p, dt, ds, 0,    0,    p->ul_lat, p->ur_lon);
      err += err_at_pixel(p, dt, ds, nl-1, 0,    p->ll_lat, p->lr_lon);
      err += err_at_pixel(p, dt, ds, 0,    ns-1, p->ur_lat, p->ul_lon);
      err += err_at_pixel(p, dt, ds, nl-1, ns-1, p->lr_lat, p->ll_lon);
    }

    meta->sar->time_shift = saved_timeOffset;
    meta->sar->slant_shift = saved_slant;

    gsl_vector_set(f,0,err);
    gsl_vector_set(f,1,err);
    return GSL_SUCCESS;
}

static void coarse_search(double t_extent_min, double t_extent_max,
                          double s_extent_min, double s_extent_max,
                          double *t_min, double *s_min,
                          struct refine_shift_params *params)
{
    double the_min = 9999999;
    double min_t=99, min_s=99;
    int i,j,k=6;
    double t_extent = t_extent_max - t_extent_min;
    double s_extent = s_extent_max - s_extent_min;
    gsl_vector *v = gsl_vector_alloc(2);
    gsl_vector *u = gsl_vector_alloc(2);
    //printf("           ");
    //for (j = 0; j <= k; ++j) {
    //    double s = s_extent_min + ((double)j)/k*s_extent;
    //    printf("%9.3f ", s);
    //}
    //printf("\n           ");
    //for (j = 0; j <= k; ++j)
    //    printf("--------- ");
    //printf("\n");
    for (i = 0; i <= k; ++i) {
        double t = t_extent_min + ((double)i)/k*t_extent;
        //printf("%9.3f | ", t);

        for (j = 0; j <= k; ++j) {
            double s = s_extent_min + ((double)j)/k*s_extent;

            gsl_vector_set(v, 0, t);
            gsl_vector_set(v, 1, s);
            getObjective(v,(void*)params, u);
            double n = gsl_vector_get(u,0);
            //printf("%9.3f ", n);
            if (n<the_min) {
                the_min=n;
                min_t=gsl_vector_get(v,0);
                min_s=gsl_vector_get(v,1);
            }
        }
        //printf("\n");
    }

    *t_min = min_t;
    *s_min = min_s;

    gsl_vector_free(v);
    gsl_vector_free(u);
}

/* this one can be used to do a time-only search (no slant adjustment)
static void coarse_search_t(double t_extent_min, double t_extent_max,
                            double *t_min, struct refine_shift_params *params)
{
    double the_min = 9999999;
    double min_t=99;
    int i,k=10;
    double t_extent = t_extent_max - t_extent_min;
    gsl_vector *v = gsl_vector_alloc(2);
    gsl_vector *u = gsl_vector_alloc(2);
    for (i = 0; i <= k; ++i) {
        double t = t_extent_min + ((double)i)/k*t_extent;

        gsl_vector_set(v, 0, t);
        gsl_vector_set(v, 1, 0);
        getObjective(v,(void*)params, u);
        double n = gsl_vector_get(u,0);
        if (n<the_min) {
            the_min=n;
            min_t=gsl_vector_get(v,0);
        }
    }

    *t_min = min_t;

    gsl_vector_free(v);
    gsl_vector_free(u);
}
*/

static void generate_start(struct refine_shift_params *params,
                           double *start_t, double *start_s)
{
    int i;

    double extent_t_min = -20;
    double extent_t_max = 20;

    double extent_s_min = -5000;
    double extent_s_max = 5000;

    double t_range = extent_t_max - extent_t_min;
    double s_range = extent_s_max - extent_s_min;

    for (i = 0; i < 12; ++i)
    {
        coarse_search(extent_t_min, extent_t_max, extent_s_min, extent_s_max,
                      start_t, start_s, params);
        //coarse_search_t(extent_t_min, extent_t_max, start_t, params);

        t_range /= 3;
        s_range /= 3;

        extent_t_min = *start_t - t_range/2;
        extent_t_max = *start_t + t_range/2;

        extent_s_min = *start_s - s_range/2;
        extent_s_max = *start_s + s_range/2;

        //printf("refining search to region: (%9.3f,%9.3f)\n"
        //       "                           (%9.3f,%9.3f)\n",
        //       extent_t_min, extent_t_max,
        //       extent_s_min, extent_s_max);
    }
}

/* small debug func
static void print_state(int iter, gsl_multiroot_fsolver *s)
{
    printf("iter = %3d   x = (%.3f %.3f)   f(x) = %.8f\n", iter,
           gsl_vector_get(s->x, 0), gsl_vector_get(s->x, 1),
           gsl_vector_get(s->f, 0));
}
*/

// This won't change the passed-in metadata -- afterwards, if you want
// to apply the found corrections, do this:
//   meta->sar->slant_shift += slant_shift_adjustment;
//   meta->sar->time_shift += time_shift_adjustment;
// usually you will only want to do this if the function returns TRUE
// (on success), though if it returns FALSE (failed) these two parameters
// will be 0
int refine_slc_geolocation_from_workreport(const char *metaName,
                                           const char *basename,
                                           meta_parameters *meta,
                                           double *time_shift_adjustment,
                                           double *slant_shift_adjustment)
{
  char *fileName = MALLOC(sizeof(char)*(10+strlen(basename)+strlen(metaName)));
  char *dirname = get_dirname(metaName);
  if (strlen(dirname)>0)
    sprintf(fileName, "%s%s", dirname, basename);
  else
    strcpy(fileName, basename);
  FILE *fp = fopen_workreport(fileName);
  FREE(dirname);
  if (!fp) {
    // failed to open workreport file
    asfPrintWarning(
"Attempted to refine the geolocation with the workreport file.  However,\n"
"this file could not be found.  To get the most accurate geolocations with\n"
"Palsar SLC data, include the workreport file in the same directory as the\n"
"data file.  First <basename>.txt is checked, then 'workreport'.\n\n"
"Proceeding... however, extremely poor geolocations may result.\n");
    *slant_shift_adjustment = 0.0;
    *time_shift_adjustment = 0.0;
    FREE(fileName);
    return FALSE;
  }

  // Read the workreport file, find the 4 corners and scene center
  // lat and long info.  We also check to make sure this is the right
  // workreport file...

  struct dataset_sum_rec dssr;
  hms_time dssr_time, summary_time;
  ymd_date dssr_date, summary_date;
  char line[512], tmp[512], *p;
  get_dssr(metaName, &dssr);
  date_dssr2date(dssr.inp_sctim, &dssr_date, &dssr_time);

  double center_lat=-999, center_lon=-999;
  double ul_lat=-999, ul_lon=-999;
  double ur_lat=-999, ur_lon=-999;
  double ll_lat=-999, ll_lon=-999;
  double lr_lat=-999, lr_lon=-999;

  while (fgets(line, 512, fp)) {
    if (strstr(line, "Img_SceneCenterDateTime")) {
      p = strchr(line, '"');
      sprintf(tmp, "%s", p+1);
      tmp[strlen(tmp)-2] = '\0';
      date_alos2date(tmp, &summary_date, &summary_time);
      if (date_difference(&dssr_date, &dssr_time,
          &summary_date, &summary_time) > 10.0)
      {
        asfPrintWarning(
"Summary file does not correspond to leader file.\n"
"  DSSR: %s\n  Summary: %s\n"
"This is most likely because you have a mismatched workreport file.  Since\n"
"by default these files are all named 'workreport' this file may have been\n"
"overwritten with another scene's workreport file.  You may wish to rename\n"
"your workreport files using the <basename>.txt convention, for example:\n"
"  %s.txt\n"
"Proceeding... however, extremely poor geolocations may result.\n",
                        dssr.inp_sctim, tmp, fileName);
        *slant_shift_adjustment = 0;
        *time_shift_adjustment = 0;
        FCLOSE(fp);
        FREE(fileName);
        return 0;
      }
    }
    else if (strstr(line, "ImageSceneCenterLatitude")) {
      p = strchr(line, '"');
      if (p) center_lat = atof(p+1);
    }
    else if (strstr(line, "ImageSceneCenterLongitude")) {
      p = strchr(line, '"');
      if (p) center_lon = atof(p+1);
    }
    else if (strstr(line, "ImageSceneLeftTopLatitude")) {
      p = strchr(line, '"');
      if (p) ul_lat = atof(p+1);
    }
    else if (strstr(line, "ImageSceneLeftTopLongitude")) {
      p = strchr(line, '"');
      if (p) ul_lon = atof(p+1);
    }
    else if (strstr(line, "ImageSceneRightTopLatitude")) {
      p = strchr(line, '"');
      if (p) ur_lat = atof(p+1);
    }
    else if (strstr(line, "ImageSceneRightTopLongitude")) {
      p = strchr(line, '"');
      if (p) ur_lon = atof(p+1);
    }
    else if (strstr(line, "ImageSceneLeftBottomLatitude")) {
      p = strchr(line, '"');
      if (p) ll_lat = atof(p+1);
    }
    else if (strstr(line, "ImageSceneLeftBottomLongitude")) {
      p = strchr(line, '"');
      if (p) ll_lon = atof(p+1);
    }
    else if (strstr(line, "ImageSceneRightBottomLatitude")) {
      p = strchr(line, '"');
      if (p) lr_lat = atof(p+1);
    }
    else if (strstr(line, "ImageSceneRightBottomLongitude")) {
      p = strchr(line, '"');
      if (p) lr_lon = atof(p+1);
    }
  }

  FREE(fileName);
  fclose(fp);

  // Did we get all that we needed?
  int ok = TRUE;
  if (center_lat == -999) {
    asfPrintStatus("Missing: ImageSceneCenterLatitude\n");
    ok = FALSE;
  }
  if (center_lon == -999) {
    asfPrintStatus("Missing: ImageSceneCenterLongitude\n");
    ok = FALSE;
  }
  if (ul_lat == -999) {
    asfPrintStatus("Missing: ImageSceneTopLeftLatitude\n");
    ok = FALSE;
  }
  if (ul_lon == -999) {
    asfPrintStatus("Missing: ImageSceneTopLeftLongitude\n");
    ok = FALSE;
  }
  if (ur_lat == -999) {
    asfPrintStatus("Missing: ImageSceneTopRightLatitude\n");
    ok = FALSE;
  }
  if (ur_lon == -999) {
    asfPrintStatus("Missing: ImageSceneTopRightLongitude\n");
    ok = FALSE;
  }
  if (ll_lat == -999) {
    asfPrintStatus("Missing: ImageSceneBottomLeftLatitude\n");
    ok = FALSE;
  }
  if (ll_lon == -999) {
    asfPrintStatus("Missing: ImageSceneBottomLeftLongitude\n");
    ok = FALSE;
  }
  if (lr_lat == -999) {
    asfPrintStatus("Missing: ImageSceneBottomRightLatitude\n");
    ok = FALSE;
  }
  if (lr_lon == -999) {
    asfPrintStatus("Missing: ImageSceneBottomRightLongitude\n");
    ok = FALSE;
  }
  if (!ok) {
    asfPrintWarning("Not all necessary information is available in the "
                    "workreport file.\nProceeding... however, extremely "
                    "poor geolocations may result.\n");
    *slant_shift_adjustment = 0.0;
    *time_shift_adjustment = 0.0;
    return FALSE;
  }

  // all required info is present-- search for min slant and time shifts

  struct refine_shift_params params;
  params.meta = meta;
  params.center_lat = center_lat;
  params.center_lon = center_lon;
  params.ul_lat = ul_lat;
  params.ul_lon = ul_lon;
  params.ur_lat = ur_lat;
  params.ur_lon = ur_lon;
  params.ll_lat = ll_lat;
  params.ll_lon = ll_lon;
  params.lr_lat = lr_lat;
  params.lr_lon = lr_lon;

  //int status;
  //int iter = 0, max_iter = 1000;
  //const gsl_multiroot_fsolver_type *T;
  //gsl_multiroot_fsolver *s;
  //gsl_error_handler_t *prev;
  //const size_t n = 2;

  //gsl_multiroot_function F = {&getObjective, n, &params};
  //gsl_vector *x = gsl_vector_alloc(n);

  double start_t, start_s;
  generate_start(&params, &start_t, &start_s);
  //printf("Starting point at (%f, %f)\n", start_t, start_s);

/*
  gsl_vector_set (x, 0, start_t);
  gsl_vector_set (x, 1, start_s);

  T = gsl_multiroot_fsolver_hybrid;
  s = gsl_multiroot_fsolver_alloc(T, n);
  gsl_multiroot_fsolver_set(s, &F, x);

  prev = gsl_set_error_handler_off();

  do {
    ++iter;
    status = gsl_multiroot_fsolver_iterate(s);

    //print_state(iter, s);

    // abort if stuck
    if (status) break;

    status = gsl_multiroot_test_residual (s->f, 1e-8);
  } while (status == GSL_CONTINUE && iter < max_iter);

  *time_shift_adjustment = gsl_vector_get(s->x, 0);
  *slant_shift_adjustment = gsl_vector_get(s->x, 1);

  gsl_vector *retrofit = gsl_vector_alloc(n);
  gsl_vector_set(retrofit, 0, *time_shift_adjustment);
  gsl_vector_set(retrofit, 1, *slant_shift_adjustment);
  gsl_vector *output = gsl_vector_alloc(n);
  getObjective(retrofit, (void*)&params, output);
  double val= gsl_vector_get(output,0);
  printf("GSL Result: %f at (%f,%f)\n",
         val, *time_shift_adjustment, *slant_shift_adjustment);
  gsl_vector_free(retrofit);
  gsl_vector_free(output);

  gsl_multiroot_fsolver_free(s);
  gsl_vector_free(x);
  gsl_set_error_handler(prev);
*/
  *time_shift_adjustment = start_t;
  *slant_shift_adjustment = start_s;

  // the rest of this is a bunch of debug code
/*
  int nl = meta->general->line_count;
  int ns = meta->general->sample_count;
  double lat,lon;

  if (meta->general->orbit_direction=='A') {

    printf("BEFORE-->\n");
    meta_get_latLon(meta, nl/2, ns/2, 0, &lat, &lon);
    printf("  CN: mgll: %f,%f - wr:%f,%f - err:%f\n",
           lat,lon,center_lat,center_lon,
           lldist(lat,center_lat,lon,center_lon));
    meta_get_latLon(meta, nl-1, 0, 0, &lat, &lon);
    printf("  UL: mgll: %f,%f - wr:%f,%f - err:%f\n",
           lat,lon,ul_lat,ul_lon,lldist(lat,ul_lat,lon,ul_lon));
    meta_get_latLon(meta, nl-1, ns-1, 0, &lat, &lon);
    printf("  UR: mgll: %f,%f - wr:%f,%f - err:%f\n",
           lat,lon,ur_lat,ur_lon,lldist(lat,ur_lat,lon,ur_lon));
    meta_get_latLon(meta, 0, 0, 0, &lat, &lon);
    printf("  LL: mgll: %f,%f - wr:%f,%f - err:%f\n",
           lat,lon,ll_lat,ll_lon,lldist(lat,ll_lat,lon,ll_lon));
    meta_get_latLon(meta, 0, ns-1, 0, &lat, &lon);
    printf("  LR: mgll: %f,%f - wr:%f,%f - err:%f\n",
           lat,lon,lr_lat,lr_lon,lldist(lat,lr_lat,lon,lr_lon));

    double saved_timeOffset = meta->sar->time_shift;
    double saved_slant = meta->sar->slant_shift;
    meta->sar->time_shift += *time_shift_adjustment;
    meta->sar->slant_shift += *slant_shift_adjustment;

    printf("AFTER-->\n");
    meta_get_latLon(meta, nl/2, ns/2, 0, &lat, &lon);
    printf("  CN: mgll: %f,%f - wr:%f,%f - err:%f\n",
           lat,lon,center_lat,center_lon,
           lldist(lat,center_lat,lon,center_lon));
    meta_get_latLon(meta, nl-1, 0, 0, &lat, &lon);
    printf("  UL: mgll: %f,%f - wr:%f,%f - err:%f\n",
           lat,lon,ul_lat,ul_lon,lldist(lat,ul_lat,lon,ul_lon));
    meta_get_latLon(meta, nl-1, ns-1, 0, &lat, &lon);
    printf("  UR: mgll: %f,%f - wr:%f,%f - err:%f\n",
           lat,lon,ur_lat,ur_lon,lldist(lat,ur_lat,lon,ur_lon));
    meta_get_latLon(meta, 0, 0, 0, &lat, &lon);
    printf("  LL: mgll: %f,%f - wr:%f,%f - err:%f\n",
           lat,lon,ll_lat,ll_lon,lldist(lat,ll_lat,lon,ll_lon));
    meta_get_latLon(meta, 0, ns-1, 0, &lat, &lon);
    printf("  LR: mgll: %f,%f - wr:%f,%f - err:%f\n",
           lat,lon,lr_lat,lr_lon,lldist(lat,lr_lat,lon,lr_lon));

    meta->sar->time_shift = saved_timeOffset;
    meta->sar->slant_shift = saved_slant;
  }
  else {
    printf("BEFORE-->\n");
    meta_get_latLon(meta, nl/2, ns/2, 0, &lat, &lon);
    printf("  CN: mgll: %f,%f - wr:%f,%f - err:%f\n",
           lat,lon,center_lat,center_lon,
           lldist(lat,center_lat,lon,center_lon));
    meta_get_latLon(meta, 0, ns-1, 0, &lat, &lon);
    printf("  UL: mgll: %f,%f - wr:%f,%f - err:%f\n",
           lat,lon,ul_lat,ul_lon,lldist(lat,ul_lat,lon,ul_lon));
    meta_get_latLon(meta, 0, 0, 0, &lat, &lon);
    printf("  UR: mgll: %f,%f - wr:%f,%f - err:%f\n",
           lat,lon,ur_lat,ur_lon,lldist(lat,ur_lat,lon,ur_lon));
    meta_get_latLon(meta, nl-1, ns-1, 0, &lat, &lon);
    printf("  LL: mgll: %f,%f - wr:%f,%f - err:%f\n",
           lat,lon,ll_lat,ll_lon,lldist(lat,ll_lat,lon,ll_lon));
    meta_get_latLon(meta, nl-1, 0, 0, &lat, &lon);
    printf("  LR: mgll: %f,%f - wr:%f,%f - err:%f\n",
           lat,lon,lr_lat,lr_lon,lldist(lat,lr_lat,lon,lr_lon));

    double saved_timeOffset = meta->sar->time_shift;
    double saved_slant = meta->sar->slant_shift;
    meta->sar->time_shift += *time_shift_adjustment;
    meta->sar->slant_shift += *slant_shift_adjustment;

    printf("AFTER-->\n");
    meta_get_latLon(meta, nl/2, ns/2, 0, &lat, &lon);
    printf("  CN: mgll: %f,%f - wr:%f,%f - err:%f\n",
           lat,lon,center_lat,center_lon,
           lldist(lat,center_lat,lon,center_lon));
    meta_get_latLon(meta, 0, ns-1, 0, &lat, &lon);
    printf("  UL: mgll: %f,%f - wr:%f,%f - err:%f\n",
           lat,lon,ul_lat,ul_lon,lldist(lat,ul_lat,lon,ul_lon));
    meta_get_latLon(meta, 0, 0, 0, &lat, &lon);
    printf("  UR: mgll: %f,%f - wr:%f,%f - err:%f\n",
           lat,lon,ur_lat,ur_lon,lldist(lat,ur_lat,lon,ur_lon));
    meta_get_latLon(meta, nl-1, ns-1, 0, &lat, &lon);
    printf("  LL: mgll: %f,%f - wr:%f,%f - err:%f\n",
           lat,lon,ll_lat,ll_lon,lldist(lat,ll_lat,lon,ll_lon));
    meta_get_latLon(meta, nl-1, 0, 0, &lat, &lon);
    printf("  LR: mgll: %f,%f - wr:%f,%f - err:%f\n",
           lat,lon,lr_lat,lr_lon,lldist(lat,lr_lat,lon,lr_lon));

    meta->sar->time_shift = saved_timeOffset;
    meta->sar->slant_shift = saved_slant;
  }
*/
  return TRUE;
}

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady


text:
// Copyright (C) 2021 Intel Corporation
// SPDX-License-Identifier: Apache-2.0

#ifndef _HEBench_ClearText_EltAdd_H_7e5fa8c2415240ea93eff148ed73539b
#define _HEBench_ClearText_EltAdd_H_7e5fa8c2415240ea93eff148ed73539b

#include <gsl/gsl>

#include "hebench/api_bridge/cpp/hebench.hpp"

#include "clear_benchmark.h"

template <class T>
class EltAdd_Benchmark : public ClearTextBenchmark
{
private:
    HEBERROR_DECLARE_CLASS_NAME(EltAdd_Benchmark)

public:
    EltAdd_Benchmark(hebench::cpp::BaseEngine &engine,
                     const hebench::APIBridge::BenchmarkDescriptor &bench_desc,
                     const hebench::APIBridge::WorkloadParams &bench_params);
    ~EltAdd_Benchmark() override;

    hebench::APIBridge::Handle encode(const hebench::APIBridge::DataPackCollection *p_parameters) override;
    void decode(hebench::APIBridge::Handle encoded_data, hebench::APIBridge::DataPackCollection *p_native) override;

    hebench::APIBridge::Handle load(const hebench::APIBridge::Handle *p_local_data, std::uint64_t count) override;
    void store(hebench::APIBridge::Handle remote_data,
               hebench::APIBridge::Handle *p_local_data, std::uint64_t count) override;

    hebench::APIBridge::Handle operate(hebench::APIBridge::Handle h_remote_packed,
                                       const hebench::APIBridge::ParameterIndexer *p_param_indexers) override;

protected:
    std::uint64_t m_vector_size;

private:
    static void eltAdd(gsl::span<T> &result,
                       const gsl::span<const T> &M0, const gsl::span<const T> &M1,
                       std::size_t element_count);
};

#include "inl/bench_eltadd.inl"

#endif // defined _HEBench_ClearText_EltAdd_H_7e5fa8c2415240ea93eff148ed73539b

max_stars_repo_name: hebench/backend-cpu-cleartext

max_issues_repo_name: hebench/backend-cpu-cleartext

max_forks_repo_name: hebench/backend-cpu-cleartext


text:
#include <gsl/gsl_statistics_double.h>

#include "qdm.h"

gsl_matrix *
qdm_bias_values(
    qdm_evaluation *h,
    qdm_evaluation *f
)
{
  size_t n = f->mcmc->r.tau->size2;
  size_t s = f->mcmc->r.s;

  double tau_low = h->parameters.tau_low;
  double tau_high = h->parameters.tau_high;

  double lower_bound = h->lower_bound;
  double upper_bound = h->upper_bound;

  gsl_matrix *r = gsl_matrix_calloc(n, s);

  fprintf(stderr, "r[%zu, %zu]\n",
      r->size1,
      r->size2
  );

  fprintf(stderr, "r.tau[%zu, %zu, %zu]\n",
      f->mcmc->r.tau->size1,
      f->mcmc->r.tau->size2,
      f->mcmc->r.tau->size3
  );

  fprintf(stderr, "r.theta[%zu, %zu, %zu]\n",
      f->mcmc->r.theta->size1,
      f->mcmc->r.theta->size2,
      f->mcmc->r.theta->size3
  );

  fprintf(stderr, "r.xi[%zu, %zu, %zu]\n",
      f->mcmc->r.xi->size1,
      f->mcmc->r.xi->size2,
      f->mcmc->r.xi->size3
  );

  for (size_t j = 0; j < n; j++) {
    for (size_t k = 0; k < s; k++) {
      double tau = qdm_ijk_get(f->mcmc->r.tau, 0, j, k);
      double value = 0;

      gsl_matrix_view theta = qdm_ijk_get_ij(h->mcmc->r.theta, k);
      gsl_vector_view mmm = gsl_matrix_row(&theta.matrix, 0);

      if (tau_low <= tau && tau <= tau_high) {
        value = qdm_tau_ispline_mmm(
            h->t,
            tau,
            &mmm.vector
        );
      } else if (tau < tau_low) {
        double is = qdm_tau_ispline_mmm(
            h->t,
            tau_low,
            &mmm.vector
        );

        double ms = qdm_tau_ispline_mmm(
            h->t,
            tau_low,
            &mmm.vector
        );

        double xi_low = qdm_ijk_get(h->mcmc->r.xi, 0, 0, k);

        value = is - tau_low / ms * (pow(tau / tau_low, -xi_low) - 1);
        value = fmax(value, lower_bound);
      } else { // tau_high < tau
        double is = qdm_tau_ispline_mmm(
            h->t,
            tau_high,
            &mmm.vector
        );

        double ms = qdm_tau_ispline_mmm(
            h->t,
            tau_high,
            &mmm.vector
        );

        double xi_high = qdm_ijk_get(h->mcmc->r.xi, 0, 1, k);

        value = is - tau_high / ms * (pow((1 - tau) / (1 - tau_high), -xi_high) - 1);
        value = fmin(value, upper_bound);
      }

      gsl_matrix_set(r, j, k, value);
    }
  }

  return r;
}

gsl_matrix *
qdm_bias_correct(
    const gsl_vector *years,
    const double month,
    const gsl_vector *days,

    const gsl_vector *y,

    qdm_evaluation *fp,
    qdm_evaluation *ho,
    qdm_evaluation *hc
)
{
  gsl_matrix *hov = qdm_bias_values(ho, fp);
  gsl_matrix *hcv = qdm_bias_values(hc, fp);
  gsl_vector *row = gsl_vector_calloc(hov->size2);

  gsl_matrix *r = gsl_matrix_calloc(hov->size1, 5);

  fprintf(stderr, "years[%zu]\n", years->size);
  fprintf(stderr, "month %f\n", month);
  fprintf(stderr, "days[%zu]\n", days->size);
  fprintf(stderr, "y[%zu]\n", y->size);
  fprintf(stderr, "hov[%zu, %zu]\n", hov->size1, hov->size2);
  fprintf(stderr, "hcv[%zu, %zu]\n", hcv->size1, hcv->size2);
  fprintf(stderr, "row[%zu]\n", row->size);

  for (size_t i = 0; i < hov->size1; i++) {
    for (size_t j = 0; j < hov->size2; j++) {
      gsl_vector_set(row, j, gsl_matrix_get(hov, i, j) + (gsl_vector_get(y, i) - gsl_matrix_get(hcv, i, j)));

      double mean = gsl_stats_mean(row->data, row->stride, row->size);
      double sd = gsl_stats_sd(row->data, row->stride, row->size);

      gsl_matrix_set(r, i, 0, gsl_vector_get(years, i));
      gsl_matrix_set(r, i, 1, month);
      gsl_matrix_set(r, i, 2, gsl_vector_get(days, i));
      gsl_matrix_set(r, i, 3, mean);
      gsl_matrix_set(r, i, 4, sd);
    }
  }

  gsl_vector_free(row);
  gsl_matrix_free(hcv);
  gsl_matrix_free(hov);

  return r;
}


text:
#ifndef BANDED_FLOAT_IMAGE_H
#define BANDED_FLOAT_IMAGE_H

#include <stdio.h>
#include <sys/types.h>

#include <glib.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_histogram.h>

#include "float_image.h"

typedef struct
{
        int nbands;
        FloatImage **images;
} BandedFloatImage;

BandedFloatImage *
banded_float_image_new(int nbands, size_t size_x, size_t size_y);

BandedFloatImage *
banded_float_image_new_with_value(int nBands, ssize_t size_x, ssize_t size_y,
				  float value);

int
banded_float_image_store (BandedFloatImage *self, const char *file,
			  float_image_byte_order_t byte_order);

void
banded_float_image_free(BandedFloatImage *self);

float
banded_float_image_get_pixel(BandedFloatImage *self, int nband,
                             ssize_t x, ssize_t y);

void
banded_float_image_set_pixel(BandedFloatImage *self, int nband,
                             ssize_t x, ssize_t y, float value);

FloatImage *
banded_float_image_get_band(BandedFloatImage *self, int nband);

BandedFloatImage *
banded_float_image_new_from_model_scaled (BandedFloatImage *model,
                                          ssize_t scale_factor);

void
banded_float_image_export_as_jpeg(BandedFloatImage *self,
                                  const char *output_name);

#endif