text:
/*
 * neuro.c
 *
 *  Created on: Oct 9, 2019
 *      Author: alexey
 */

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_rng.h>
#include <png.h>
#include <malloc.h>

#include "neuro.h"

typedef struct _neuro_skeleton {
	int input_l;
	int output_l;
	int hidden_l;
	int number_of_hidden;
	double alpha;
	gsl_matrix *input_m;
	gsl_matrix *output_m;
	gsl_matrix **hidden_m;
	gsl_matrix **dropuot;
	gsl_matrix **layers;
	gsl_matrix **layers_delta;
	gsl_matrix *temp;
	gsl_matrix *train;
	gsl_matrix *error;
	int layers_numb;
	activation activ;
	activation *activ_l;
} neuro_skeleton;

const gsl_rng_type * T;
gsl_rng * r;

extern unsigned long int gsl_rng_default_seed;

#define CLEAN_BLOB() bayrepo_clean_neuro((void *) blob); \
		return NULL
#define CLEAN_IF_NULL(x) if (!blob->x) { \
		CLEAN_BLOB(); \
	}

#define ISDEBUGINFO_BEG() if (getenv("MATRIXD") && !strcmp(getenv("MATRIXD"),"1")) {
#define ISDEBUGINFO_END() }
#define DEBUGINFO(A, B, C) ISDEBUGINFO_BEG()\
		bayrepo_print_matrix(A, B, C); \
		ISDEBUGINFO_END()

#define MAX_LAYER_NAME_LEN 1024
#define MAX_BUFFER_LEN 4096

static void bayrepo_print_matrix(gsl_matrix *a, const char *matrix_name,
		int index_matrix) {
	int index, jndex;
	if (index_matrix >= 0) {
		printf("Matrix %s[%d]:\n", matrix_name, index_matrix);
	} else {
		printf("Matrix %s:\n", matrix_name);
	}
	for (index = 0; index < a->size1; index++) {
		printf("==");
		for (jndex = 0; jndex < a->size2; jndex++) {
			printf("%.3f   ", gsl_matrix_get(a, index, jndex));
		}
		printf("==\n");
	}
}

static void bayrepo_layers_clean(void *blob) {
	if (blob) {
		neuro_skeleton *data = (neuro_skeleton *) blob;
		if (data->layers) {
			int index = 0;
			for (index = 0; index < data->layers_numb; index++) {
				if (data->layers[index])
					gsl_matrix_free(data->layers[index]);
			}
			free(data->layers);
			data->layers = NULL;
		}
		if (data->layers_delta) {
			int index = 0;
			for (index = 0; index < data->layers_numb; index++) {
				if (data->layers_delta[index])
					gsl_matrix_free(data->layers_delta[index]);
			}
			data->layers_numb = 0;
			free(data->layers_delta);
			data->layers_delta = NULL;
		}
		data->layers_numb = 0;
	}
}

void bayrepo_clean_neuro(void *blob) {
	if (blob) {
		neuro_skeleton *data = (neuro_skeleton *) blob;
		if (data->input_m) {
			gsl_matrix_free(data->input_m);
			data->input_m = NULL;
		}
		if (data->output_m) {
			gsl_matrix_free(data->output_m);
			data->output_m = NULL;
		}
		if (data->train) {
			gsl_matrix_free(data->train);
			data->train = NULL;
		}
		if (data->temp) {
			gsl_matrix_free(data->temp);
			data->temp = NULL;
		}
		if (data->error) {
			gsl_matrix_free(data->error);
			data->error = NULL;
		}
		if (data->activ_l) {
			free(data->activ_l);
			data->activ_l = NULL;
		}
		if (data->number_of_hidden > 0 && data->hidden_m) {
			int i;
			for (i = 0; i < data->number_of_hidden; i++) {
				if (data->hidden_m[i]) {
					gsl_matrix_free(data->hidden_m[i]);
					data->hidden_m[i] = NULL;
				}
			}
			free(data->hidden_m);
		}
		if (data->number_of_hidden > 0 && data->dropuot) {
			int i;
			for (i = 0; i < data->number_of_hidden; i++) {
				if (data->dropuot[i]) {
					gsl_matrix_free(data->dropuot[i]);
					data->dropuot[i] = NULL;
				}
			}
			free(data->dropuot);
		}
		bayrepo_layers_clean(blob);
		free(data);
	}
	if (r) {
		gsl_rng_free(r);
		r = NULL;
	}
}

void bayrepo_set_layer_activ(void * blob, int layer_number, activation activ) {
	if (blob) {
		neuro_skeleton *data = (neuro_skeleton *) blob;
		if (data->number_of_hidden) {
			if (layer_number <= data->number_of_hidden) {
				data->activ_l[layer_number] = activ;
			}
		} else {
			if (!layer_number) {
				data->activ_l[0] = activ;
			}
		}
	}
}

static int bayrepo_random_matrix(gsl_matrix *a, int m, int n, int is_negative) {
	int index, jndex;

	if (!r)
		return -1;

	for (index = 0; index < m; index++) {
		for (jndex = 0; jndex < n; jndex++) {
			gsl_matrix_set(a, index, jndex,
					gsl_rng_uniform(r) - (is_negative ? 0.5 : 0.0));
		}
	}

	return 0;
}

void *bayrepo_init_neuro(int input, int output, int hidden, int hidden_num,
		double alpha, activation activ) {
	int index;

	gsl_rng_default_seed = time(NULL);

	T = gsl_rng_knuthran2;
	r = gsl_rng_alloc(T);

	if (hidden_num < 0) {
		hidden_num = 0;
	}
	neuro_skeleton *blob = calloc(1, sizeof(neuro_skeleton));
	if (!blob)
		return NULL;
	blob->input_m = gsl_matrix_alloc(1, input);
	CLEAN_IF_NULL(input_m);
	blob->train = gsl_matrix_alloc(1, output);
	CLEAN_IF_NULL(train);
	blob->error = gsl_matrix_alloc(1, output);
	CLEAN_IF_NULL(error);

	if (hidden_num) {
		blob->output_m = gsl_matrix_alloc(hidden, output);
		CLEAN_IF_NULL(output_m);
		blob->hidden_m = (gsl_matrix **) calloc(hidden_num,
				sizeof(gsl_matrix *));
		CLEAN_IF_NULL(hidden_m);
		blob->dropuot = (gsl_matrix **) calloc(hidden_num,
				sizeof(gsl_matrix *));
		CLEAN_IF_NULL(dropuot);
		for (index = 1; index < hidden_num; index++) {
			blob->hidden_m[index] = gsl_matrix_alloc(hidden, hidden);
			CLEAN_IF_NULL(hidden_m[index]);
			blob->dropuot[index] = gsl_matrix_alloc(1, hidden);
			CLEAN_IF_NULL(dropuot[index]);
		}
		blob->hidden_m[0] = gsl_matrix_alloc(input, hidden);
		CLEAN_IF_NULL(hidden_m[0]);
		blob->dropuot[0] = gsl_matrix_alloc(1, hidden);
		CLEAN_IF_NULL(dropuot[0]);
		blob->temp = gsl_matrix_alloc(1, hidden);
		CLEAN_IF_NULL(temp);
	} else {
		blob->output_m = gsl_matrix_alloc(input, output);
		CLEAN_IF_NULL(output_m);
	}
	blob->alpha = alpha;
	blob->input_l = input;
	blob->output_l = output;
	blob->hidden_l = hidden;
	blob->number_of_hidden = hidden_num;

	gsl_matrix_set_zero(blob->input_m);
	gsl_matrix_set_zero(blob->train);
	if (hidden_num) {
		if (bayrepo_random_matrix(blob->output_m, hidden, output,
				(activ == RELU ? 0 : 1)) < 0) {
			CLEAN_BLOB()
;		}
		for (index = 1; index < hidden_num; index++) {
			if (bayrepo_random_matrix(blob->hidden_m[index], hidden, hidden, (activ==RELU?0:1))<0) {
				CLEAN_BLOB();
			}
		}
		if (bayrepo_random_matrix(blob->hidden_m[0], input, hidden, (activ==RELU?0:1))<0) {
			CLEAN_BLOB();
		}
	} else {
		if (bayrepo_random_matrix(blob->output_m, input, output, (activ==RELU?0:1))<0) {
			CLEAN_BLOB();
		}
	}

	int lyr_nmb = 1 + blob->number_of_hidden ? (blob->number_of_hidden + 1) : 0;

	blob->layers = (gsl_matrix **) calloc(lyr_nmb, sizeof(gsl_matrix *));
	CLEAN_IF_NULL(layers);

	if (hidden_num) {
		for (index = 1; index < hidden_num; index++) {
			blob->layers[index] = gsl_matrix_alloc(1, hidden);
			CLEAN_IF_NULL(layers[index]);
		}
		blob->layers[0] = gsl_matrix_alloc(1, hidden);
		CLEAN_IF_NULL(layers[0]);
		blob->layers[lyr_nmb - 1] = gsl_matrix_alloc(1, output);
		CLEAN_IF_NULL(layers[lyr_nmb-1]);
	} else {
		blob->layers[0] = gsl_matrix_alloc(1, output);
		CLEAN_IF_NULL(layers[0]);
	}

	blob->layers_numb = lyr_nmb;
	blob->layers_delta = (gsl_matrix **) calloc(lyr_nmb, sizeof(gsl_matrix *));
	CLEAN_IF_NULL(layers_delta);

	if (hidden_num) {
		for (index = 1; index < hidden_num; index++) {
			blob->layers_delta[index] = gsl_matrix_alloc(1, hidden);
			CLEAN_IF_NULL(layers_delta[index]);
		}
		blob->layers_delta[0] = gsl_matrix_alloc(1, hidden);
		CLEAN_IF_NULL(layers_delta[0]);
		blob->layers_delta[lyr_nmb - 1] = gsl_matrix_alloc(1, output);
		CLEAN_IF_NULL(layers_delta[lyr_nmb-1]);
	} else {
		blob->layers_delta[0] = gsl_matrix_alloc(1, output);
		CLEAN_IF_NULL(layers_delta[0]);
	}
	blob->activ = activ;
	blob->activ_l = calloc(blob->number_of_hidden + 1, sizeof(activation));
	CLEAN_IF_NULL(activ_l);
	for (index = 0; index < (blob->number_of_hidden + 1); index++) {
		if (index == blob->number_of_hidden) {
			blob->activ_l[index] = NOACTIV;
		} else {
			blob->activ_l[index] = DEFLT;
		}
	}

	return (void *) blob;
}

static void bayrepo_fill_dropout(gsl_matrix *a, int m, int n) {
	int index, jndex;
	if (!r) {
		for (index = 0; index < m; index++) {
			for (jndex = 0; jndex < n; jndex++) {
				gsl_matrix_set(a, index, jndex, 1.0);
			}
		}
		return;
	}

	for (index = 0; index < m; index++) {
		for (jndex = 0; jndex < n; jndex++) {
			gsl_matrix_set(a, index, jndex, gsl_rng_uniform_int(r, 2) * 1.0);
		}
	}

	return;
}

void bayrepo_fill_input(void *blob, int position, double scaled_value) {
	if (blob) {
		neuro_skeleton *data = (neuro_skeleton *) blob;
		if (position < data->input_l) {
			gsl_matrix_set(data->input_m, 0, position, scaled_value);
		}
	}
}

void bayrepo_fill_train(void *blob, int position, double scaled_value) {
	if (blob) {
		neuro_skeleton *data = (neuro_skeleton *) blob;
		if (position < data->output_l) {
			gsl_matrix_set(data->train, 0, position, scaled_value);
		}
	}
}

void bayrepo_fill_hidden(void *blob, int index, int position_x, int position_y,
		double scaled_value) {
	if (blob) {
		neuro_skeleton *data = (neuro_skeleton *) blob;
		gsl_matrix_set(data->hidden_m[index], position_x, position_y,
				scaled_value);
	}
}

void bayrepo_fill_outm(void *blob, int position_x, int position_y,
		double scaled_value) {
	if (blob) {
		neuro_skeleton *data = (neuro_skeleton *) blob;
		gsl_matrix_set(data->output_m, position_x, position_y, scaled_value);
	}
}

static void bayrepo_zero_layers(void *blob) {
	if (blob) {
		neuro_skeleton *data = (neuro_skeleton *) blob;
		int index;
		for (index = 0; index < data->layers_numb; index++) {
			gsl_matrix_set_zero(data->layers[index]);
		}
	}
}

static double bayrepo_activation_func(double elem, neuro_skeleton *data,
		int lyn) {
	activation act = data->activ;
	if ((lyn < (data->number_of_hidden + 1)) && (data->activ_l[lyn] != DEFLT)) {
		act = data->activ_l[lyn];
	}
	switch (act) {
	case RELU:
		return elem > 0.0 ? elem : 0.0;
	case TANH:
		return tanh(elem);
	case SIGMOID:
		return 1.0 / (1.0 + exp(-elem));
	default:
		return elem;
	}
}

activation bayrepo_get_layer_func(void *blob, int lyn) {
	if (blob) {
		neuro_skeleton *data = (neuro_skeleton *) blob;
		if ((lyn < (data->number_of_hidden + 1))
				&& (data->activ_l[lyn] != DEFLT)) {
			return data->activ_l[lyn];
		}
	}
	return NOACTIV;
}

activation bayrepo_get_sublayer_func(void *blob, int lyn) {
	if (blob) {
		neuro_skeleton *data = (neuro_skeleton *) blob;
		if (lyn < (data->number_of_hidden + 1)) {
			return data->activ_l[lyn];
		}
	}
	return NOACTIV;
}

static double bayrepo_activation_deriv(double elem, neuro_skeleton *data,
		int lyn) {
	activation act = data->activ;
	if ((lyn < (data->number_of_hidden + 1)) && (data->activ_l[lyn] != DEFLT)) {
		act = data->activ_l[lyn];
	}
	switch (act) {
	case RELU:
		return elem > 0.0 ? 1.0 : 0.0;
	case TANH:
		return 1.0 - (elem * elem);
	case SIGMOID:
		return elem * (1.0 - elem);
	default:
		return elem;
	}
}

static void bayrepo_matix_customize(gsl_matrix *a, int size,
		neuro_skeleton *data, int lyn) {
	int index;
	for (index = 0; index < size; index++) {
		gsl_matrix_set(a, 0, index,
				bayrepo_activation_func(gsl_matrix_get(a, 0, index), data,
						lyn));
	}
}

static void bayrepo_matix_deriv(gsl_matrix *a, int size, neuro_skeleton *data,
		int lyn) {
	int index;
	for (index = 0; index < size; index++) {
		gsl_matrix_set(a, 0, index,
				bayrepo_activation_deriv(gsl_matrix_get(a, 0, index), data,
						lyn));
	}
}

static void bayrepo_query_internal(void *blob, int dropout) {
	int index;
	if (blob) {
		ISDEBUGINFO_BEG()
			printf("=====================Query=========================\n");
		ISDEBUGINFO_END()
		neuro_skeleton *data = (neuro_skeleton *) blob;
		bayrepo_zero_layers(blob);

		if (data->number_of_hidden) {
			if (dropout == 1) {
				for (index = 0; index < data->number_of_hidden; index++) {
					bayrepo_fill_dropout(data->dropuot[index], 1,
							data->hidden_l);
				}
			}
			int cnt = 0;
			gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, data->input_m,
					data->hidden_m[0], 0.0, data->layers[0]);
			bayrepo_matix_customize(data->layers[0], data->hidden_l, data, 0);
			if (dropout == 1) {
				gsl_matrix_mul_elements(data->layers[0], data->dropuot[0]);
				gsl_matrix_scale(data->layers[0], 2.0);
			}
			cnt++;
			for (index = 1; index < data->number_of_hidden; index++) {
				gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0,
						data->layers[index - 1], data->hidden_m[index], 0.0,
						data->layers[index]);
				bayrepo_matix_customize(data->layers[index], data->hidden_l,
						data, index);
				if (dropout == 1) {
					gsl_matrix_mul_elements(data->layers[index],
							data->dropuot[index]);
					gsl_matrix_scale(data->layers[index], 2.0);
				}
				cnt++;
			}
			gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0,
					data->layers[cnt - 1], data->output_m, 0.0,
					data->layers[cnt]);
			bayrepo_matix_customize(data->layers[cnt], data->output_l, data,
					cnt);

		} else {
			gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, data->input_m,
					data->output_m, 0.0, data->layers[0]);
			bayrepo_matix_customize(data->layers[0], data->output_l, data, 0);
		}

		DEBUGINFO(data->input_m, "INPUTM", -1);
		if (data->number_of_hidden) {
			for (index = 0; index < data->number_of_hidden; index++) {
				DEBUGINFO(data->hidden_m[index], "HIDDEN", index);
			}
		}
		DEBUGINFO(data->output_m, "OUTPUTM", -1);
		for (index = 0; index < data->layers_numb; index++) {
			DEBUGINFO(data->layers[index], "LAYER", index);
		}
	}
}

void bayrepo_query(void *blob) {
	bayrepo_query_internal(blob, 0);
}

double bayrepo_get_result(void *blob, int position) {
	if (blob) {
		neuro_skeleton *data = (neuro_skeleton *) blob;
		if (data->layers && position < data->output_l) {
			return gsl_matrix_get(data->layers[data->layers_numb - 1], 0,
					position);
		}
	}
	return -10000.0;
}

double bayrepo_get_sum(gsl_matrix *a) {
	double result = 0.0;
	int index, jndex = 0;
	for (index = 0; index < a->size1; index++) {
		for (jndex = 0; jndex < a->size2; jndex++) {
			result += gsl_matrix_get(a, index, jndex);
		}
	}
	return result;
}

void bayrepo_train(void *blob, int epoch, int use_dropout) {
	int index;
	if (blob) {
		while (epoch--) {
			ISDEBUGINFO_BEG()
				printf(
						"=====================Epoch %d=========================\n",
						epoch);
			ISDEBUGINFO_END()
			neuro_skeleton *data = (neuro_skeleton *) blob;
			bayrepo_query_internal(blob, use_dropout);
			gsl_matrix_memcpy(data->layers_delta[data->layers_numb - 1],
					data->layers[data->layers_numb - 1]);
			ISDEBUGINFO_BEG()
				gsl_matrix_memcpy(data->error,
						data->layers[data->layers_numb - 1]);
			ISDEBUGINFO_END();
			gsl_matrix_sub(data->layers_delta[data->layers_numb - 1],
					data->train);
			ISDEBUGINFO_BEG()
				gsl_matrix_sub(data->error, data->train);
				gsl_matrix_mul_elements(data->error, data->error);
				printf("=======>Error=%f\n", bayrepo_get_sum(data->error));
			ISDEBUGINFO_END();
			if (data->number_of_hidden) {
				int index;
				for (index = (data->layers_numb - 2); index >= 0; index--) {
					gsl_matrix_set_zero(data->temp);
					gsl_matrix_memcpy(data->temp, data->layers[index]);
					bayrepo_matix_deriv(data->temp, data->hidden_l, data,
							index);
					if (index == (data->layers_numb - 2)) {
						gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0,
								data->layers_delta[index + 1], data->output_m,
								0.0, data->layers_delta[index]);
					} else {
						gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0,
								data->layers_delta[index + 1],
								data->hidden_m[index + 1], 0.0,
								data->layers_delta[index]);
					}
					gsl_matrix_mul_elements(data->layers_delta[index],
							data->temp);
				}

				if (use_dropout) {
					for (index = 0; index < data->number_of_hidden; index++) {
						gsl_matrix_mul_elements(data->layers_delta[index],
								data->dropuot[index]);
					}
				}

				for (index = (data->layers_numb - 1); index >= 0; index--) {
					if (index == (data->layers_numb - 1)) {
						gsl_blas_dgemm(CblasTrans, CblasNoTrans, -data->alpha,
								data->layers[index - 1],
								data->layers_delta[index], 1.0, data->output_m);
					} else if (index == 0) {
						gsl_blas_dgemm(CblasTrans, CblasNoTrans, -data->alpha,
								data->input_m, data->layers_delta[index], 1.0,
								data->hidden_m[0]);
					} else {
						gsl_blas_dgemm(CblasTrans, CblasNoTrans, -data->alpha,
								data->layers[index - 1],
								data->layers_delta[index], 1.0,
								data->hidden_m[index]);
					}
				}
			} else {
				gsl_blas_dgemm(CblasTrans, CblasNoTrans, -data->alpha,
						data->input_m, data->layers_delta[0], 1.0,
						data->output_m);
			}
			for (index = 0; index < data->layers_numb; index++) {
				DEBUGINFO(data->layers_delta[index], "LAYER_DELTA", index);
			}
		}
	}
}

static void byrepo_png_write_data(png_structp png_ptr, png_bytep data,
		png_size_t length) {

	bayrepo_mem_encode* p = (bayrepo_mem_encode*) png_get_io_ptr(png_ptr);
	size_t nsize = p->size + length;

	if (p->buffer)
		p->buffer = realloc(p->buffer, nsize);
	else
		p->buffer = malloc(nsize);

	if (!p->buffer)
		png_error(png_ptr, "Write Error");

	memcpy(p->buffer + p->size, data, length);
	p->size += length;
}

static void bayrepo_png_flush(png_structp png_ptr) {
}

int bayrepo_write_matrix(void *blob, FILE *fp, int width, int height,
		bayrepo_mem_encode *buffer) {
	int code = 0;
	if (buffer) {
		buffer->buffer = NULL;
		buffer->size = 0;
	}
	png_structp png_ptr = NULL;
	png_infop info_ptr = NULL;
	png_bytep row = NULL;
	if ((blob == NULL) || (fp == NULL)) {
		return -1;
	}

	png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
	if (png_ptr == NULL) {
		code = -3;
		goto finalise;
	}

	info_ptr = png_create_info_struct(png_ptr);
	if (info_ptr == NULL) {
		code = -4;
		goto finalise;
	}

	if (setjmp(png_jmpbuf(png_ptr))) {
		code = -5;
		goto finalise;
	}

	if (buffer) {
		png_set_write_fn(png_ptr, buffer, byrepo_png_write_data,
				bayrepo_png_flush);
	} else {
		png_init_io(png_ptr, fp);
	}

	neuro_skeleton *data = (neuro_skeleton *) blob;

	int picSize = (height + 3) + (height + 3) * data->number_of_hidden;

	png_set_IHDR(png_ptr, info_ptr, width, picSize, 8, PNG_COLOR_TYPE_GRAY,
	PNG_INTERLACE_NONE,
	PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE);

	png_write_info(png_ptr, info_ptr);

	row = (png_bytep) malloc(3 * width * sizeof(png_byte));
	if (!row) {
		code = -5;
		goto finalise;
	}

	int index = 0;
	for (index = 0; index < data->number_of_hidden + 1; index++) {
		gsl_matrix * m = NULL;
		if (data->number_of_hidden) {
			m = (index == data->number_of_hidden) ?
					data->output_m : data->hidden_m[index];
		} else {
			m = data->output_m;
		}

		double rangeMax = gsl_matrix_max(m);
		double rangeMin = gsl_matrix_min(m);

		if (((width / m->size1) == 0) || ((height / m->size2) == 0)) {
			return -2;
		}

		int sizeX = width / m->size1;
		int sizeY = height / m->size2;

		int x, y;
		for (y = 0; y < height; y++) {
			for (x = 0; x < width; x++) {
				int curX = x / sizeX;
				int curY = y / sizeY;
				if ((curX >= m->size1) || (curY >= m->size2)) {
					row[x] = (png_byte) 255;
				} else {
					png_byte res_color = (png_byte) ((gsl_matrix_get(m, curX,
							curY) - rangeMin) / (rangeMax - rangeMin) * 255.0);
					if (res_color > 255.0)
						res_color = 255.0;
					row[x] = (png_byte) (255.0 - res_color);
				}
			}
			png_write_row(png_ptr, row);
		}
		for (y = 0; y < 3; y++) {
			for (x = 0; x < width; x++) {
				if (y != 1) {
					row[x] = (png_byte) 255;
				} else {
					row[x] = (png_byte) 0;
				}
			}
			png_write_row(png_ptr, row);
		}
	}

	png_write_end(png_ptr, NULL);

	finalise: if (info_ptr != NULL)
		png_free_data(png_ptr, info_ptr, PNG_FREE_ALL, -1);
	if (png_ptr != NULL)
		png_destroy_write_struct(&png_ptr, (png_infopp) NULL);
	if (row != NULL)
		free(row);
	return code;

}

void bayrepo_print_matrix_custom(void *blob, bayrepo_decorator *decor) {
	char layer_name[MAX_LAYER_NAME_LEN];
	if (blob) {
		neuro_skeleton *data = (neuro_skeleton *) blob;
		int index = 0;
		for (index = 0; index < (data->number_of_hidden + 1); index++) {
			gsl_matrix * m = NULL;
			if (data->number_of_hidden) {
				m = (index == data->number_of_hidden) ?
						data->output_m : data->hidden_m[index];
				if (index == data->number_of_hidden) {
					snprintf(layer_name, MAX_LAYER_NAME_LEN, "Output layer");
				} else {
					snprintf(layer_name, MAX_LAYER_NAME_LEN, "Hidden layer %d",
							index);
				}
			} else {
				m = data->output_m;
				snprintf(layer_name, MAX_LAYER_NAME_LEN, "Output layer");
			}
			int posX, posY;
			if (decor->table_header)
				decor->table_header(decor->user_data, layer_name, m->size1);
			for (posX = 0; posX < m->size2; posX++) {
				if (decor->pre_column)
					decor->pre_column(decor->user_data);
				for (posY = 0; posY < m->size1; posY++) {
					if (decor->print_func)
						decor->print_func(decor->user_data,
								gsl_matrix_get(m, posY, posX));
				}
				if (decor->post_column)
					decor->post_column(decor->user_data);
			}
			if (decor->table_footer)
				decor->table_footer(decor->user_data);
		}
	}
}

int input_l;
int output_l;
int hidden_l;
int number_of_hidden;
double alpha;
gsl_matrix *input_m;
gsl_matrix *output_m;
gsl_matrix **hidden_m;
gsl_matrix **dropuot;
gsl_matrix **layers;
gsl_matrix **layers_delta;
gsl_matrix *temp;
gsl_matrix *train;
gsl_matrix *error;
int layers_numb;
activation activ;
activation *activ_l;

/*
 * Save format
 *
 *  [int input_l]\n
 *  [int output_l]\n
 *  [int hidden_l]\n
 *  [int number_of_hidden]\n
 *  [double alpha]\n
 *  [int activ]\n
 *  [ACTBEG]\n
 *  [int layer_number]:[int activ]\n
 *  [ACTEND]\n
 *  [HIDDEN_BEG]\n
 *  [int X]:[int Y]:[double value]\n
 *  [HIDDEN_END]
 *  [OUTPUT_BEG]\n
 *  [int X]:[int Y]:[double value]\n
 *  [OUTPUT_END]\n
 */

static char *bayrepo_reallocate_buffer(char *buffer, char *newdata, int *len,
		int newlen) {
	buffer = (char *) realloc(buffer, *len + newlen);
	if (!buffer) {
		return NULL;
	}
	memcpy(buffer + *len, newdata, newlen);
	*len = *len + newlen;
	return buffer;
}

int bayrepo_save_to_buffer(void *blob, char **buffer) {
	char buff[MAX_LAYER_NAME_LEN];
	int index;
	if (blob) {
		int len = 0;
		char *ptr = NULL;
		neuro_skeleton *data = (neuro_skeleton *) blob;
		snprintf(buff, MAX_LAYER_NAME_LEN, "%d\n", data->input_l);
		ptr = bayrepo_reallocate_buffer(ptr, (char *) buff, &len, strlen(buff));
		snprintf(buff, MAX_LAYER_NAME_LEN, "%d\n", data->output_l);
		ptr = bayrepo_reallocate_buffer(ptr, (char *) buff, &len, strlen(buff));
		snprintf(buff, MAX_LAYER_NAME_LEN, "%d\n", data->hidden_l);
		ptr = bayrepo_reallocate_buffer(ptr, (char *) buff, &len, strlen(buff));
		snprintf(buff, MAX_LAYER_NAME_LEN, "%d\n", data->number_of_hidden);
		ptr = bayrepo_reallocate_buffer(ptr, (char *) buff, &len, strlen(buff));
		snprintf(buff, MAX_LAYER_NAME_LEN, "%f\n", data->alpha);
		ptr = bayrepo_reallocate_buffer(ptr, (char *) buff, &len, strlen(buff));
		snprintf(buff, MAX_LAYER_NAME_LEN, "%d\n", (int) data->activ);
		ptr = bayrepo_reallocate_buffer(ptr, (char *) buff, &len, strlen(buff));
		ptr = bayrepo_reallocate_buffer(ptr, (char *) "[ACTBEG]\n", &len,
				strlen("[ACTBEG]\n"));
		for (index = 0; index < (data->number_of_hidden + 1); index++) {
			snprintf(buff, MAX_LAYER_NAME_LEN, "%d:%d\n", index,
					(int) data->activ_l[index]);
			ptr = bayrepo_reallocate_buffer(ptr, (char *) buff, &len,
					strlen(buff));
		}
		ptr = bayrepo_reallocate_buffer(ptr, (char *) "[ACTEND]\n", &len,
				strlen("[ACTEND]\n"));
		for (index = 0; index < data->number_of_hidden; index++) {
			ptr = bayrepo_reallocate_buffer(ptr, (char *) "[HIDDEN_BEG]\n",
					&len, strlen("[HIDDEN_BEG]\n"));
			int x, y;
			for (x = 0; x < data->hidden_m[index]->size1; x++) {
				for (y = 0; y < data->hidden_m[index]->size2; y++) {
					snprintf(buff, MAX_LAYER_NAME_LEN, "%d:%d:%f\n", x, y,
							gsl_matrix_get(data->hidden_m[index], x, y));
					ptr = bayrepo_reallocate_buffer(ptr, (char *) buff, &len,
							strlen(buff));
				}
			}
			ptr = bayrepo_reallocate_buffer(ptr, (char *) "[HIDDEN_END]\n",
					&len, strlen("[HIDDEN_END]\n"));
		}
		ptr = bayrepo_reallocate_buffer(ptr, (char *) "[OUTPUT_BEG]\n", &len,
				strlen("[OUTPUT_BEG]\n"));
		int x, y;
		for (x = 0; x < data->output_m->size1; x++) {
			for (y = 0; y < data->output_m->size2; y++) {
				snprintf(buff, MAX_LAYER_NAME_LEN, "%d:%d:%f\n", x, y,
						gsl_matrix_get(data->output_m, x, y));
				ptr = bayrepo_reallocate_buffer(ptr, (char *) buff, &len,
						strlen(buff));
			}
		}
		ptr = bayrepo_reallocate_buffer(ptr, (char *) "[OUTPUT_END]\n", &len,
				strlen("[OUTPUT_END]\n"));
		*buffer = ptr;
		return len;
	}
	return 0;
}

void *bayrepo_restore_buffer(char *buffer, int buffer_len) {
	FILE *fp = fmemopen((void*) buffer, buffer_len, "r");
	int input_l = 0;
	int output_l = 0;
	int hidden_l = 0;
	int numb_of_hidden = 0;
	activation act = RELU;
	int act_beg = 0;
	double alpha = 0.0;
	void *blob = NULL;
	int is_hidden = 0;
	int hidden_index = 0;
	int is_out = 0;
	if (fp) {
		char result_buf[MAX_BUFFER_LEN];
		int counter = 0;

		while (!feof(fp)) {
			if (fgets(result_buf, MAX_BUFFER_LEN, fp)) {
				switch (counter) {
				case 0: {
					result_buf[strlen(result_buf) - 1] = 0;
					input_l = atoi(result_buf);
				}
					break;
				case 1: {
					result_buf[strlen(result_buf) - 1] = 0;
					output_l = atoi(result_buf);
				}
					break;
				case 2: {
					result_buf[strlen(result_buf) - 1] = 0;
					hidden_l = atoi(result_buf);
				}
					break;
				case 3: {
					result_buf[strlen(result_buf) - 1] = 0;
					numb_of_hidden = atoi(result_buf);
				}
					break;
				case 4: {
					char *tptr;
					result_buf[strlen(result_buf) - 1] = 0;
					alpha = strtod(result_buf, &tptr);
				}
					break;
				case 5: {
					result_buf[strlen(result_buf) - 1] = 0;
					act = (activation) atoi(result_buf);
					if (input_l > 0 && output_l > 0 && hidden_l >= 0
							&& numb_of_hidden >= 0 && alpha > 0.0) {
						blob = bayrepo_init_neuro(input_l, output_l, hidden_l,
								numb_of_hidden, alpha, act);
					}
				}
					break;
				default: {
					if (!blob)
						return NULL;
					if (strstr(result_buf, "[ACTBEG]")) {
						act_beg = 1;
					} else if (strstr(result_buf, "[ACTEND]")) {
						act_beg = 0;
					} else if (act_beg) {
						result_buf[strlen(result_buf) - 1] = 0;
						char *ptr = strchr(result_buf, ':');
						if (ptr) {
							*ptr = 0;
							ptr++;
							int ly = atoi(result_buf);
							int ac = atoi(ptr);
							bayrepo_set_layer_activ(blob, ly, (activation) ac);
						}
					} else if (strstr(result_buf, "[HIDDEN_BEG]")) {
						is_hidden = 1;
					} else if (strstr(result_buf, "[HIDDEN_END]")) {
						hidden_index++;
						is_hidden = 0;
					} else if (is_hidden) {
						result_buf[strlen(result_buf) - 1] = 0;
						char *ptr1 = strchr(result_buf, ':');
						if (ptr1) {
							*ptr1 = 0;
							ptr1++;
							char *ptr2 = strchr(ptr1, ':');
							if (ptr2) {
								*ptr2 = 0;
								ptr2++;
								char *tptr = NULL;
								int x = atoi(result_buf);
								int y = atoi(ptr1);
								double vl = strtod(ptr2, &tptr);
								bayrepo_fill_hidden(blob, hidden_index, x, y,
										vl);
							}
						}
					} else if (strstr(result_buf, "[OUTPUT_BEG]")) {
						is_out = 1;
					} else if (strstr(result_buf, "[OUTPUT_END]")) {
						hidden_index++;
						is_out = 0;
					} else if (is_out) {
						result_buf[strlen(result_buf) - 1] = 0;
						char *ptr1 = strchr(result_buf, ':');
						if (ptr1) {
							*ptr1 = 0;
							ptr1++;
							char *ptr2 = strchr(ptr1, ':');
							if (ptr2) {
								*ptr2 = 0;
								ptr2++;
								char *tptr = NULL;
								int x = atoi(result_buf);
								int y = atoi(ptr1);
								double vl = strtod(ptr2, &tptr);
								bayrepo_fill_outm(blob, x, y, vl);
							}
						}
					}
				}
				}
			}
			counter++;
		}
		return blob;

	}
	return NULL;
}


text:
/* -*- Mode: C++; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
/*
 *     Copyright 2018-Present Couchbase, Inc.
 *
 *   Use of this software is governed by the Business Source License included
 *   in the file licenses/BSL-Couchbase.txt.  As of the Change Date specified
 *   in that file, in accordance with the Business Source License, use of this
 *   software will be governed by the Apache License, Version 2.0, included in
 *   the file licenses/APL2.txt.
 */
#pragma once

#include <gsl/gsl-lite.hpp>
#include <nlohmann/json.hpp>
#include <platform/dirutils.h>
#include <vector>

class Event;

/**
 * The Module class represents the configuration for a single module.
 * See ../README.md for more information.
 */
class Module {
public:
    Module() = delete;
    Module(const nlohmann::json& json,
           const std::string& srcRoot,
           const std::string& objRoot);

    void createHeaderFile();

    /**
     * The name of the module
     */
    std::string name;
    /**
     * The lowest identifier for the audit events in this module. All
     * audit descriptor defined for this module MUST be within the range
     * [start, start + max_events_per_module]
     */
    int64_t start;
    /**
     * The name of the file containing the audit descriptors for this
     * module.
     */
    std::string file;
    /**
     * The JSON data describing the audit descriptors for this module
     */
    nlohmann::json json;
    /**
     * Is this module enterprise only?
     */
    bool enterprise = false;

    /**
     * If present this is the name of a C headerfile to generate with
     * #defines for all audit identifiers for the module.
     */
    std::string header;

    /**
     * A list of all of the events defined for this module
     */
    std::vector<std::unique_ptr<Event>> events;

protected:
    /**
     * Add the event to the list of events for the module
     *
     * @param event the event to add
     * @throws std::invalid_argument if the event is outside the legal range
     *                               for the module
     */
    void addEvent(std::unique_ptr<Event> event);

    /// Parse the event descriptor file and add all of the events into
    /// the list of events
    void parseEventDescriptorFile();
};

ext: h


text:
#include <iostream>
#include <math.h>
#include <Windows.h>
#include <string>
#include <vector>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/legacy/compat.hpp>
#include <opencv2/opencv.hpp>
#include <Eigen/Dense>
#include <gsl/gsl_multifit.h>
#include <stdbool.h>
#include <CLEyeMulticam.h>

#include "SerialClass.h"

int i = 0;

using namespace cv;
using namespace std;
using namespace Eigen;

void setupScreen();
void setupCamera(int i);
bool processCamera(int i, Mat imgOriginal);
Vector3d rotateX(Vector3d v, double angle);
Vector3d rotateY(Vector3d v, double angle);
Vector3d rotateZ(Vector3d v, double angle);
bool polynomialfit(int obs, int degree, double *dx, double *dy, double *store);
void fpsPrint();
void triangulate();
void trajectoryCalc();
void arduinoSend();

int demo = 1;
int dataCollectEnabled = 0;
vector<double> xCoords;
vector<double> yCoords;
vector<double> zCoords;
double XZcoeff[3] = { 0 };
double XYcoeff[2] = { 0 };

int estXY[2] = { 0 }; // sum of estimated values

#include "CLEye.h" // This needs to be here b/c it uses the functions, namespaces, variables above

char *port = "COM4";

class Camera {
public:
	double rot[3];
	double rot2[3];
	Vector3d pos;
	double planeDist;
	double planeX;
	double planeY;
	double screenX;
	double screenY;
	Vector3d ballDir;
};
Camera cameras[2];

CLEyeCameraCapture *cam[2];

int minHue, maxHue, minSat, maxSat, minValue, maxValue;

int catcherHeight = 15; // Height of robot end-effector, in cm (at what height the ball will be caught)
double minX = 8.7; // Minimum distance of center of catcher to origin, in cm
double minY = 8.7;
double maxX = 52.7;
double maxY = 48.2;

int lastTick = 0;

char buffer[256] = "";

Serial* arduino;

#include "Functions.h"

max_stars_repo_name: BaranCODE/BallCatcher

max_issues_repo_name: BaranCODE/BallCatcher

max_forks_repo_name: BaranCODE/BallCatcher


text:
// SPDX-License-Identifier: MIT
// The MIT License (MIT)
//
// Copyright (c) 2014-2018, Institute for Software & Systems Engineering
// Copyright (c) 2018-2019, Johannes Leupolz
//
// 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.

#ifndef PEMC_GENERIC_TRAVERSER_PATH_TRACKER_H_
#define PEMC_GENERIC_TRAVERSER_PATH_TRACKER_H_

#include <vector>
#include <gsl/span>
#include <cstdint>
#include <atomic>
#include <stack>
#include <limits>

#include "pemc/basic/tsc_index.h"
#include "pemc/basic/label.h"
#include "pemc/basic/model_capacity.h"
#include "pemc/basic/raw_memory.h"
#include "pemc/formula/formula.h"

namespace pemc {
  struct PathFrame {
    // the offset into the states array where the index of the frame's first state is stored.
    int32_t indexOfFirstStateIndexEntry; //offset
    // the number of states the frame consists of.
    int32_t count;
  };

  ///   When enumerating all states of a model in a depth-first fashion, we have to store the next states (that are computed all
	///   at once) somewhere while also being able to generate the counter example. When check a new state, a new frame is allocated
	///   on the stack and all unknown successor states are stored in that frame. We then take the topmost state, compute its
	///   successor states, and so on. When a formula violation is detected, the counter example consists of the last states of each
	///   frame on the stack. When a frame has been fully enumerated without detecting a formula violation, the stack frame is
	///   removed, the topmost state of the topmost frame is removed, and the new topmost state is checked.
  class PathTracker {
  private:
      std::vector<PathFrame> pathFrames;
      std::vector<StateIndex> stateIndexEntries;

      ///   The lowest index of a splittable frame. -1 if no frame is splittable.
      int32_t lowestSplittableFrame = -1;

      //   Maximal number of stateIndexEntries.
      int32_t capacity = 1 << 20;

      ///   Finds the next splittable frame, if any.
      void updateLowestSplittableFrame();

  public:
      PathTracker(int32_t _capacity);

      ///   Indicates whether the pathFrames can be split.
      bool canSplit();

      ///   Gets the number of path frames.
      int32_t getPathFrameCount();

      ///   Clears all pathFrames and its stateIndexEntries.
      void clear();

      ///   Pushes a new frame onto the pathFrames stack.
      void pushFrame();

      /// Adds the stateIndex to the current path frame.
      void pushStateIndex(StateIndex stateIndex);

  		///   Tries to get the topmost stateIndex if there is one.
      ///   Returns false to indicate that the path tracker was empty and no stateIndex was returned.
      bool tryGetStateIndex(StateIndex& stateIndex);

      bool splitWork(PathTracker& other);

      ///   Gets the path the path tracker currently represents, i.e.,
      ///   returns the sequence of topmost states of each frame, starting with
      ///   the oldest one.
      std::vector<StateIndex> getCurrentPath();
  };

}

#endif  // PEMC_GENERIC_TRAVERSER_PATH_TRACKER_H_

max_stars_repo_name: joleuger/pemc

max_issues_repo_name: joleuger/pemc

max_forks_repo_name: joleuger/pemc


text:
/***************************************************************************/
/*                                                                         */
/* QR_decomp.c - LU Decomposition class for mruby                          */
/* Copyright (C) 2015 Paolo Bosetti                                        */
/* paolo[dot]bosetti[at]unitn.it                                           */
/* Department of Industrial Engineering, University of Trento              */
/*                                                                         */
/* This library is free software.  You can redistribute it and/or          */
/* modify it under the terms of the GNU GENERAL PUBLIC LICENSE 2.0.        */
/*                                                                         */
/* This library 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           */
/* Artistic License 2.0 for more details.                                  */
/*                                                                         */
/* See the file LICENSE                                                    */
/*                                                                         */
/***************************************************************************/

#include <gsl/gsl_linalg.h>
#include <stdio.h>
#include "matrix.h"
#include "vector.h"
#include "QR_decomp.h"

#ifndef MIN
#define MAX(a, b)                                                              \
  ({                                                                           \
    __typeof__(a) _a = (a);                                                    \
    __typeof__(b) _b = (b);                                                    \
    _a < _b ? _a : _b;                                                         \
  })
#endif

#pragma mark -
#pragma mark • Utilities

// Garbage collector handler, for play_data struct
// if play_data contains other dynamic data, free it too!
// Check it with GC.start
void qr_decomp_destructor(mrb_state *mrb, void *p_) {
  qr_decomp_data_s *lu = (qr_decomp_data_s *)p_;
  gsl_matrix_free(lu->mat);
  gsl_vector_free(lu->tau);
  free(lu);
};

// Creating data type and reference for GC, in a const struct
const struct mrb_data_type qr_decomp_data_type = {"qr_decomp_data",
                                                  qr_decomp_destructor};

// Utility function for getting the struct out of the wrapping IV @data
void mrb_qr_decomp_get_data(mrb_state *mrb, mrb_value self,
                            qr_decomp_data_s **data) {
  mrb_value data_value;
  data_value = mrb_iv_get(mrb, self, mrb_intern_lit(mrb, "@data"));

  // Loading data from data_value into p_data:
  Data_Get_Struct(mrb, data_value, &qr_decomp_data_type, *data);
  if (!*data)
    mrb_raise(mrb, E_RUNTIME_ERROR, "Could not access @data");
}

#pragma mark -
#pragma mark • Initializations

// Data Initializer C function (not exposed!)
static mrb_value mrb_qr_initialize(mrb_state *mrb, mrb_value self) {
  mrb_value data_value;            // this IV holds the data
  qr_decomp_data_s *p_data = NULL; // pointer to the C struct
  mrb_value matrix;
  gsl_matrix *p_mat = NULL;
  mrb_int size1, size2;

  mrb_get_args(mrb, "o", &matrix);
  if (!mrb_obj_is_kind_of(mrb, matrix, mrb_class_get(mrb, "Matrix"))) {
    mrb_raise(mrb, E_ARGUMENT_ERROR, "Argument must be a Matrix");
  }

  mrb_matrix_get_data(mrb, matrix, &p_mat);
  // if (p_mat->size1 < p_mat->size2) {
  //   mrb_raise(mrb, E_ARGUMENT_ERROR, "Argument must be square or horizontal
  //   Matrix");
  // }
  size1 = p_mat->size1;
  size2 = p_mat->size2;

  data_value = mrb_iv_get(mrb, self, mrb_intern_lit(mrb, "@data"));
  // if @data already exists, free its content:
  if (!mrb_nil_p(data_value)) {
    Data_Get_Struct(mrb, data_value, &qr_decomp_data_type, p_data);
    free(p_data);
  }
  // Allocate and zero-out the data struct:
  p_data = (qr_decomp_data_s *)malloc(sizeof(qr_decomp_data_s));
  if (!p_data) {
    mrb_raise(mrb, E_RUNTIME_ERROR, "Could not allocate @data");
  }
  p_data->size1 = size1;
  p_data->size2 = size2;
  p_data->minsize = MIN(size1, size2);
  p_data->mat = gsl_matrix_calloc(size1, size2);
  p_data->tau = gsl_vector_calloc(p_data->minsize);

  // copy argument matrix into local object data
  gsl_matrix_memcpy(p_data->mat, p_mat);
  // invert in-place
  gsl_linalg_QR_decomp(p_data->mat, p_data->tau);
  mrb_iv_set(mrb, self, mrb_intern_lit(mrb, "@size1"), mrb_fixnum_value(size1));
  mrb_iv_set(mrb, self, mrb_intern_lit(mrb, "@size2"), mrb_fixnum_value(size2));
  mrb_iv_set(mrb, self, mrb_intern_lit(mrb, "@minsize"),
             mrb_fixnum_value(p_data->minsize));
  mrb_iv_set(mrb, self, mrb_intern_lit(mrb, "@residuals"), mrb_nil_value());
  // Wrap struct into @data:
  mrb_iv_set(mrb, self, mrb_intern_lit(mrb, "@data"), // set @data
             mrb_obj_value( // with value hold in struct
                 Data_Wrap_Struct(mrb, mrb->object_class, &qr_decomp_data_type,
                                  p_data)));
  return mrb_nil_value();
}

#pragma mark -
#pragma mark • Accessors

static mrb_value mrb_qr_residuals(mrb_state *mrb, mrb_value self) {
  return mrb_iv_get(mrb, self, mrb_intern_lit(mrb, "@residuals"));
}

static mrb_value mrb_qr_minsize(mrb_state *mrb, mrb_value self) {
  return mrb_iv_get(mrb, self, mrb_intern_lit(mrb, "@minsize"));
}

static mrb_value mrb_qr_size1(mrb_state *mrb, mrb_value self) {
  return mrb_iv_get(mrb, self, mrb_intern_lit(mrb, "@size1"));
}

static mrb_value mrb_qr_size2(mrb_state *mrb, mrb_value self) {
  return mrb_iv_get(mrb, self, mrb_intern_lit(mrb, "@size2"));
}

static mrb_value mrb_qr_matrix(mrb_state *mrb, mrb_value self) {
  mrb_value result;
  qr_decomp_data_s *p_data = NULL;
  gsl_matrix *p_res = NULL;
  mrb_value args[2];

  // call utility for unwrapping @data into p_data:
  mrb_qr_decomp_get_data(mrb, self, &p_data);
  args[0] = mrb_fixnum_value(p_data->size1);
  args[1] = mrb_fixnum_value(p_data->size2);
  result = mrb_obj_new(mrb, mrb_class_get(mrb, "Matrix"), 2, args);
  mrb_matrix_get_data(mrb, result, &p_res);
  gsl_matrix_memcpy(p_res, p_data->mat);
  return result;
}

static mrb_value mrb_qr_tau(mrb_state *mrb, mrb_value self) {
  mrb_value result;
  qr_decomp_data_s *p_data = NULL;
  gsl_vector *p_res = NULL;
  mrb_value args[1];

  // call utility for unwrapping @data into p_data:
  mrb_qr_decomp_get_data(mrb, self, &p_data);
  args[0] = mrb_fixnum_value(p_data->minsize);
  result = mrb_obj_new(mrb, mrb_class_get(mrb, "Vector"), 1, args);
  mrb_vector_get_data(mrb, result, &p_res);
  gsl_vector_memcpy(p_res, p_data->tau);
  return result;
}

#pragma mark -
#pragma mark • Operations

// int gsl_linalg_QR_solve (const gsl_matrix * QR, const gsl_vector * tau, const
// gsl_vector * b, gsl_vector * x)
static mrb_value mrb_qr_solve(mrb_state *mrb, mrb_value self) {
  mrb_value result, b_vec;
  qr_decomp_data_s *p_data = NULL;
  gsl_vector *p_result = NULL, *p_b = NULL;
  mrb_value args[1];

  mrb_get_args(mrb, "o", &b_vec);
  if (!mrb_obj_is_kind_of(mrb, b_vec, mrb_class_get(mrb, "Vector"))) {
    mrb_raise(mrb, E_ARGUMENT_ERROR, "Argument must be a Vector");
  }

  // call utility for unwrapping @data into p_data:
  mrb_qr_decomp_get_data(mrb, self, &p_data);
  if (p_data->size1 != p_data->size2) {
    mrb_raise(mrb, E_QR_DECOMP_ERROR, "Matrix must be square");
  }

  args[0] = mrb_fixnum_value(p_data->tau->size);
  result = mrb_obj_new(mrb, mrb_class_get(mrb, "Vector"), 1, args);
  mrb_vector_get_data(mrb, result, &p_result);
  mrb_vector_get_data(mrb, b_vec, &p_b);

  if (gsl_linalg_QR_solve(p_data->mat, p_data->tau, p_b, p_result)) {
    mrb_raise(mrb, E_QR_DECOMP_ERROR, "Singular matrix");
  }
  mrb_iv_set(mrb, self, mrb_intern_lit(mrb, "@residuals"), mrb_nil_value());
  return result;
}

// int gsl_linalg_QR_lssolve (const gsl_matrix * QR, const gsl_vector * tau,
// const gsl_vector * b, gsl_vector * x, gsl_vector * residual)
static mrb_value mrb_qr_lssolve(mrb_state *mrb, mrb_value self) {
  mrb_value result, b_vec, residuals;
  qr_decomp_data_s *p_data = NULL;
  gsl_vector *p_result = NULL, *p_b = NULL, *p_residuals;
  mrb_value args[1];

  mrb_get_args(mrb, "o", &b_vec);
  if (!mrb_obj_is_kind_of(mrb, b_vec, mrb_class_get(mrb, "Vector"))) {
    mrb_raise(mrb, E_ARGUMENT_ERROR, "Argument must be a Vector");
  }

  // call utility for unwrapping @data into p_data:
  mrb_qr_decomp_get_data(mrb, self, &p_data);
  if (p_data->size1 <= p_data->size2) {
    mrb_raise(mrb, E_QR_DECOMP_ERROR,
              "Matrix must have more rows than columns");
  }

  args[0] = mrb_fixnum_value(p_data->size2);
  result = mrb_obj_new(mrb, mrb_class_get(mrb, "Vector"), 1, args);
  mrb_vector_get_data(mrb, result, &p_result);
  mrb_vector_get_data(mrb, b_vec, &p_b);

  args[0] = mrb_fixnum_value(p_b->size);
  residuals = mrb_obj_new(mrb, mrb_class_get(mrb, "Vector"), 1, args);
  mrb_vector_get_data(mrb, residuals, &p_residuals);
  mrb_iv_set(mrb, self, mrb_intern_lit(mrb, "@residuals"), residuals);

  if (gsl_linalg_QR_lssolve(p_data->mat, p_data->tau, p_b, p_result,
                            p_residuals)) {
    mrb_raise(mrb, E_QR_DECOMP_ERROR, "Singular matrix");
  }
  return result;
}

#pragma mark -
#pragma mark • Gem setup

void mrb_gsl_qr_decomp_init(mrb_state *mrb) {
  struct RClass *lu;

  mrb_load_string(mrb, "class QRDecompError < Exception; end");

  lu = mrb_define_class(mrb, "QRDecomp", mrb->object_class);
  mrb_define_method(mrb, lu, "residuals", mrb_qr_residuals, MRB_ARGS_NONE());
  mrb_define_method(mrb, lu, "matrix", mrb_qr_matrix, MRB_ARGS_NONE());
  mrb_define_method(mrb, lu, "tau", mrb_qr_tau, MRB_ARGS_NONE());
  mrb_define_method(mrb, lu, "size1", mrb_qr_size1, MRB_ARGS_NONE());
  mrb_define_method(mrb, lu, "size2", mrb_qr_size2, MRB_ARGS_NONE());
  mrb_define_method(mrb, lu, "minsize", mrb_qr_minsize, MRB_ARGS_NONE());

  mrb_define_method(mrb, lu, "initialize", mrb_qr_initialize, MRB_ARGS_REQ(1));
  mrb_define_method(mrb, lu, "solve", mrb_qr_solve, MRB_ARGS_REQ(1));
  mrb_define_method(mrb, lu, "lssolve", mrb_qr_lssolve, MRB_ARGS_REQ(1));
}

ext: c

max_stars_repo_name: UniTN-Mechatronics/mruby-gsl

max_issues_repo_name: UniTN-Mechatronics/mruby-gsl

max_forks_repo_name: UniTN-Mechatronics/mruby-gsl


text:
/*
ODE: a program to get optime Runge-Kutta and multi-steps methods.

Copyright 2011-2019, Javier Burguete Tolosa.

Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:

	1. Redistributions of source code must retain the above copyright notice,
		this list of conditions and the following disclaimer.

	2. Redistributions in binary form must reproduce the above copyright notice,
		this list of conditions and the following disclaimer in the
		documentation and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY Javier Burguete Tolosa ``AS IS'' AND ANY EXPRESS OR
IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
SHALL Javier Burguete Tolosa OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

/**
 * \file rk_3_3.c
 * \brief Source file to optimize Runge-Kutta 3 steps 3rd order methods.
 * \author Javier Burguete Tolosa.
 * \copyright Copyright 2011-2019.
 */
#define _GNU_SOURCE
#include <string.h>
#include <math.h>
#include <libxml/parser.h>
#include <glib.h>
#include <libintl.h>
#include <gsl/gsl_rng.h>
#include "config.h"
#include "utils.h"
#include "optimize.h"
#include "rk.h"
#include "rk_3_3.h"

#define DEBUG_RK_3_3 0          ///< macro to debug.

/**
 * Function to obtain the coefficients of a 3 steps 3rd order Runge-Kutta
 * method.
 */
int
rk_tb_3_3 (Optimize * optimize) ///< Optimize struct.
{
  long double *tb, *r;
#if DEBUG_RK_3_3
  fprintf (stderr, "rk_tb_3_3: start\n");
#endif
  tb = optimize->coefficient;
  r = optimize->random_data;
  t3 (tb) = 1.L;
  t1 (tb) = r[0];
  t2 (tb) = r[1];
  b32 (tb) = (1.L / 3.L - 0.5L * t1 (tb)) / (t2 (tb) * (t2 (tb) - t1 (tb)));
  if (isnan (b32 (tb)))
    return 0;
  b31 (tb) = (1.L / 3.L - 0.5L * t2 (tb)) / (t1 (tb) * (t1 (tb) - t2 (tb)));
  if (isnan (b31 (tb)))
    return 0;
  b21 (tb) = 1 / 6.L / (b32 (tb) * t1 (tb));
  if (isnan (b21 (tb)))
    return 0;
  rk_b_3 (tb);
#if DEBUG_RK_3_3
  rk_print_tb (optimize, "rk_tb_3_3", stderr);
  fprintf (stderr, "rk_tb_3_3: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 3 steps 3rd order, 4th order in
 * equations depending only in time, Runge-Kutta method.
 */
int
rk_tb_3_3t (Optimize * optimize)        ///< Optimize struct.
{
  long double *tb, *r;
#if DEBUG_RK_3_3
  fprintf (stderr, "rk_tb_3_3t: start\n");
#endif
  tb = optimize->coefficient;
  r = optimize->random_data;
  t3 (tb) = 1.L;
  t1 (tb) = r[0];
  t2 (tb) = (4.L * t1 (tb) - 3.L) / (6.L * t1 (tb) - 4.L);
  if (isnan (t2 (tb)))
    return 0;
  b32 (tb) = (1.L / 3.L - 0.5L * t1 (tb)) / (t2 (tb) * (t2 (tb) - t1 (tb)));
  if ((b32 (tb)))
    return 0;
  b31 (tb) = (1.L / 3.L - 0.5L * t2 (tb)) / (t1 (tb) * (t1 (tb) - t2 (tb)));
  if (isnan (b31 (tb)))
    return 0;
  b21 (tb) = 1 / 6.L / (b32 (tb) * t1 (tb));
  if (isnan (b21 (tb)))
    return 0;
  rk_b_3 (tb);
#if DEBUG_RK_3_3
  rk_print_tb (optimize, "rk_tb_3_3t", stderr);
  fprintf (stderr, "rk_tb_3_3t: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 3 steps 2nd-3rd order Runge-Kutta
 * pair.
 */
int
rk_tb_3_3p (Optimize * optimize)        ///< Optimize struct.
{
  long double *tb;
#if DEBUG_RK_3_3
  fprintf (stderr, "rk_tb_3_3p: start\n");
#endif
  if (!rk_tb_3_3 (optimize))
    return 0;
  tb = optimize->coefficient;
  e31 (tb) = 0.5L / t1 (tb);
  rk_e_3 (tb);
#if DEBUG_RK_3_3
  rk_print_e (optimize, "rk_tb_3_3p", stderr);
  fprintf (stderr, "rk_tb_3_3p: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 3 steps 2nd-3rd order, 2nd-4th order
 * in equations depending only in time, Runge-Kutta pair.
 */
int
rk_tb_3_3tp (Optimize * optimize)       ///< Optimize struct.
{
  long double *tb;
#if DEBUG_RK_3_3
  fprintf (stderr, "rk_tb_3_3tp: start\n");
#endif
  if (!rk_tb_3_3t (optimize))
    return 0;
  tb = optimize->coefficient;
  e31 (tb) = 0.5L / t1 (tb);
  rk_e_3 (tb);
#if DEBUG_RK_3_3
  rk_print_e (optimize, "rk_tb_3_3tp", stderr);
  fprintf (stderr, "rk_tb_3_3tp: end\n");
#endif
  return 1;
}

/**
 * Function to calculate the objective function of a 3 steps 3rd order
 * Runge-Kutta method.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_3_3 (RK * rk)   ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_3_3
  fprintf (stderr, "rk_objective_tb_3_3: start\n");
#endif
  tb = rk->tb->coefficient;
  o = fminl (0.L, b20 (tb));
  if (b21 (tb) < 0.L)
    o += b21 (tb);
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b31 (tb) < 0.L)
    o += b31 (tb);
  if (b32 (tb) < 0.L)
    o += b32 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L + fmaxl (1.L, fmaxl (t1 (tb), t2 (tb)));
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_3_3
  fprintf (stderr, "rk_objective_tb_3_3: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_3_3: end\n");
#endif
  return o;
}

/**
 * Function to calculate the objective function of a 3 steps 3rd order, 4th
 * order in equations depending only in time, Runge-Kutta method.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_3_3t (RK * rk)  ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_3_3
  fprintf (stderr, "rk_objective_tb_3_3t: start\n");
#endif
  tb = rk->tb->coefficient;
  o = fminl (0.L, b20 (tb));
  if (b21 (tb) < 0.L)
    o += b21 (tb);
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b31 (tb) < 0.L)
    o += b31 (tb);
  if (b32 (tb) < 0.L)
    o += b32 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L + fmaxl (1.L, fmaxl (t1 (tb), t2 (tb)));
#if DEBUG_RK_3_3
  fprintf (stderr, "rk_objective_tb_3_3: optimal=%Lg\n", o);
#endif
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_3_3
  fprintf (stderr, "rk_objective_tb_3_3t: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_3_3t: end\n");
#endif
  return o;
}

/**
 * Function to calculate the objective function of a 3 steps 2nd-3rd order
 * Runge-Kutta pair.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_3_3p (RK * rk)  ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_3_3
  fprintf (stderr, "rk_objective_tb_3_3p: start\n");
#endif
  tb = rk->tb->coefficient;
  o = fminl (0.L, b20 (tb));
  if (b21 (tb) < 0.L)
    o += b21 (tb);
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b31 (tb) < 0.L)
    o += b31 (tb);
  if (b32 (tb) < 0.L)
    o += b32 (tb);
  if (e30 (tb) < 0.L)
    o += e30 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L + fmaxl (1.L, fmaxl (t1 (tb), t2 (tb)));
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_3_3
  fprintf (stderr, "rk_objective_tb_3_3p: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_3_3p: end\n");
#endif
  return o;
}

/**
 * Function to calculate the objective function of a 3 steps 2nd-3rd order,
 * 3rd-4th order in equations depending only in time, Runge-Kutta pair.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_3_3tp (RK * rk) ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_3_3
  fprintf (stderr, "rk_objective_tb_3_3tp: start\n");
#endif
  tb = rk->tb->coefficient;
  o = fminl (0.L, b20 (tb));
  if (b21 (tb) < 0.L)
    o += b21 (tb);
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b31 (tb) < 0.L)
    o += b31 (tb);
  if (b32 (tb) < 0.L)
    o += b32 (tb);
  if (e30 (tb) < 0.L)
    o += e30 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L + fmaxl (1.L, fmaxl (t1 (tb), t2 (tb)));
#if DEBUG_RK_3_3
  fprintf (stderr, "rk_objective_tb_3_3p: optimal=%Lg\n", o);
#endif
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_3_3
  fprintf (stderr, "rk_objective_tb_3_3tp: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_3_3tp: end\n");
#endif
  return o;
}

max_stars_repo_name: jburguete/ode

max_issues_repo_name: jburguete/ode

max_forks_repo_name: jburguete/ode


text:
/* specfunc/bessel_Ynu.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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 2 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include "gsl_sf_bessel.h"

#include "error.h"

#include "bessel.h"
#include "bessel_olver.h"
#include "bessel_temme.h"

/* Perform forward recurrence for Y_nu(x) and Y'_nu(x)
 *
 *        Y_{nu+1} =  nu/x Y_nu - Y'_nu
 *       Y'_{nu+1} = -(nu+1)/x Y_{nu+1} + Y_nu
 */
#if 0
static
int
bessel_Y_recur(const double nu_min, const double x, const int kmax,
               const double Y_start, const double Yp_start,
	       double * Y_end, double * Yp_end)
{
  double x_inv = 1.0/x;
  double nu = nu_min;
  double Y_nu  = Y_start;
  double Yp_nu = Yp_start;
  int k;

  for(k=1; k<=kmax; k++) {
    double nuox = nu*x_inv;
    double Y_nu_save = Y_nu;
    Y_nu  = -Yp_nu + nuox * Y_nu;
    Yp_nu = Y_nu_save - (nuox+x_inv) * Y_nu;
    nu += 1.0;
  }
  *Y_end  = Y_nu;
  *Yp_end = Yp_nu;
  return GSL_SUCCESS;
}
#endif

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int
gsl_sf_bessel_Ynu_e(double nu, double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x <= 0.0 || nu < 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(nu > 50.0) {
    return gsl_sf_bessel_Ynu_asymp_Olver_e(nu, x, result);
  }
  else {
    /* -1/2 <= mu <= 1/2 */
    int N = (int)(nu + 0.5);
    double mu = nu - N;

    gsl_sf_result Y_mu, Y_mup1;
    int stat_mu;
    double Ynm1;
    double Yn;
    double Ynp1;
    int n;

    if(x < 2.0) {
      /* Determine Ymu, Ymup1 directly. This is really
       * an optimization since this case could as well
       * be handled by a call to gsl_sf_bessel_JY_mu_restricted(),
       * as below.
       */
      stat_mu = gsl_sf_bessel_Y_temme(mu, x, &Y_mu, &Y_mup1);
    }
    else {
      /* Determine Ymu, Ymup1 and Jmu, Jmup1.
       */
      gsl_sf_result J_mu, J_mup1;
      stat_mu = gsl_sf_bessel_JY_mu_restricted(mu, x, &J_mu, &J_mup1, &Y_mu, &Y_mup1);
    }

    /* Forward recursion to get Ynu, Ynup1.
     */
    Ynm1 = Y_mu.val;
    Yn   = Y_mup1.val;
    for(n=1; n<=N; n++) {
      Ynp1 = 2.0*(mu+n)/x * Yn - Ynm1;
      Ynm1 = Yn;
      Yn   = Ynp1;
    }

    result->val  = Ynm1; /* Y_nu */
    result->err  = (N + 1.0) * fabs(Ynm1) * (fabs(Y_mu.err/Y_mu.val) + fabs(Y_mup1.err/Y_mup1.val));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(Ynm1);

    return stat_mu;
  }
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_bessel_Ynu(const double nu, const double x)
{
  EVAL_RESULT(gsl_sf_bessel_Ynu_e(nu, x, &result));
}


text:
#pragma once

#include <gsl/span>

#include <cstddef>

namespace imageview {

// Class representing a color in RGB24 color space.
class RGB24 {
 public:
  constexpr RGB24() noexcept = default;
  // Construct an RGB24 color from the given channel components.
  constexpr RGB24(unsigned char red, unsigned char green, unsigned char blue) noexcept
      : red(red), green(green), blue(blue) {}

  unsigned char red = 0;
  unsigned char green = 0;
  unsigned char blue = 0;
};

// Implementation of the PixelFormat concept for RGB24 pixel format.
// In this pixel format the color is represented via 3 8-bit integers,
// specifying the red, green and blue channels. When serializing to /
// deserializing from  a byte array, the order of the channels is RGB (i.e.
// not BGR).
class PixelFormatRGB24 {
 public:
  using color_type = RGB24;
  static constexpr int kBytesPerPixel = 3;

  constexpr color_type read(gsl::span<const std::byte, kBytesPerPixel> data) const;

  constexpr void write(const color_type& color, gsl::span<std::byte, kBytesPerPixel> data) const;
};

constexpr bool operator==(const RGB24& lhs, const RGB24& rhs) {
  return lhs.red == rhs.red && lhs.green == rhs.green && lhs.blue == rhs.blue;
}

constexpr bool operator!=(const RGB24& lhs, const RGB24& rhs) { return !(lhs == rhs); }

constexpr PixelFormatRGB24::color_type PixelFormatRGB24::read(gsl::span<const std::byte, kBytesPerPixel> data) const {
  return color_type(static_cast<unsigned char>(data[0]), static_cast<unsigned char>(data[1]),
                    static_cast<unsigned char>(data[2]));
}

constexpr void PixelFormatRGB24::write(const color_type& color, gsl::span<std::byte, kBytesPerPixel> data) const {
  data[0] = static_cast<std::byte>(color.red);
  data[1] = static_cast<std::byte>(color.green);
  data[2] = static_cast<std::byte>(color.blue);
}

}  // namespace imageview

max_stars_repo_name: alexanderbelous/imageview

max_issues_repo_name: alexanderbelous/imageview

max_forks_repo_name: alexanderbelous/imageview


text:
//          Copyright Jean Pierre Cimalando 2019.
// Distributed under the Boost Software License, Version 1.0.
//    (See accompanying file LICENSE.md or copy at
//          http://www.boost.org/LICENSE_1_0.txt)

#pragma once
#include <gsl/gsl>

template <class Ch> bool string_starts_with(gsl::basic_string_span<const Ch> text, gsl::basic_string_span<const Ch> prefix);
template <class Ch> bool string_ends_with(gsl::basic_string_span<const Ch> text, gsl::basic_string_span<const Ch> suffix);

#include "strings.tcc"

max_stars_repo_name: jpcima/smf-dsp

max_issues_repo_name: jpcima/smf-dsp

max_forks_repo_name: jpcima/smf-dsp


text:
/*	$Id$ */
/*
 * Copyright (c) 2014, 2015 Kristaps Dzonsons <kristaps@kcons.eu>
 *
 * Permission to use, copy, modify, and distribute this software for any
 * purpose with or without fee is hereby granted, provided that the above
 * copyright notice and this permission notice appear in all copies.
 *
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
 */
#include <assert.h>
#include <inttypes.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>

#include <cairo.h>
#include <cairo-pdf.h>
#include <cairo-ps.h>
#include <gtk/gtk.h>
#include <gsl/gsl_multifit.h>
#include <kplot.h>

#include "extern.h"

enum	savetype {
	SAVE_PDF,
	SAVE_PS,
	SAVE_EPS
};

static int
savepng(const gchar *fname, const struct curwin *c)
{
	cairo_surface_t	*surf;
	cairo_t		*cr;
	cairo_status_t	 st;

	g_debug("%p: Saving: %s", c, fname);
	surf = cairo_image_surface_create
		(CAIRO_FORMAT_ARGB32, 600, 400);
	st = cairo_surface_status(surf);
	if (CAIRO_STATUS_SUCCESS != st) {
		g_debug("%s", cairo_status_to_string(st));
		cairo_surface_destroy(surf);
		return(0);
	}

	cr = cairo_create(surf);
	cairo_surface_destroy(surf);

	st = cairo_status(cr);
	if (CAIRO_STATUS_SUCCESS != st) {
		g_debug("%s", cairo_status_to_string(st));
		cairo_destroy(cr);
		return(0);

	}

	cairo_set_source_rgb(cr, 1.0, 1.0, 1.0);
	cairo_rectangle(cr, 0.0, 0.0, 600.0, 400.0);
	cairo_fill(cr);
	kplot_draw(c->views[c->view], 600.0, 400.0, cr);

	st = cairo_surface_write_to_png(cairo_get_target(cr), fname);

	if (CAIRO_STATUS_SUCCESS != st) {
		g_debug("%s", cairo_status_to_string(st));
		cairo_destroy(cr);
		return(0);
	}

	cairo_destroy(cr);
	return(1);
}

static int
savepdf(const gchar *fname, struct curwin *c, enum savetype type)
{
	cairo_surface_t	*surf;
	cairo_t		*cr;
	cairo_status_t	 st;
	struct kplotcfg	*cfg;
	struct kdatacfg	*datas;
	size_t		 datasz, i, j;
	int		 rc;
	double		 svtic, svaxis, svborder, svgrid, sv, svticln;
	const double	 w = 72 * 6, h = 72 * 5;

	g_debug("%p: Saving: %s", c, fname);
	switch (type) {
	case (SAVE_PDF):
		surf = cairo_pdf_surface_create(fname, w, h);
		break;
	case (SAVE_EPS):
		surf = cairo_ps_surface_create(fname, w, h);
		cairo_ps_surface_set_eps(surf, 1);
		break;
	case (SAVE_PS):
		surf = cairo_ps_surface_create(fname, w, h);
		cairo_ps_surface_set_eps(surf, 0);
		break;
	default:
		abort();
	}

	st = cairo_surface_status(surf);
	if (CAIRO_STATUS_SUCCESS != st) {
		g_debug("%s", cairo_status_to_string(st));
		cairo_surface_destroy(surf);
		return(0);
	}

	cr = cairo_create(surf);
	cairo_surface_destroy(surf);

	st = cairo_status(cr);
	if (CAIRO_STATUS_SUCCESS != st) {
		g_debug("%s", cairo_status_to_string(st));
		cairo_destroy(cr);
		return(0);

	}

	cfg = kplot_get_plotcfg(c->views[c->view]);

	svtic = cfg->ticlabelfont.sz;
	svaxis = cfg->axislabelfont.sz;
	svborder = cfg->borderline.sz;
	svgrid = cfg->gridline.sz;
	svticln = cfg->ticline.sz;
	sv = 0.0; /* Silence compiler. */

	cfg->ticlabelfont.sz = 9.0;
	cfg->axislabelfont.sz = 9.0;
	cfg->borderline.sz = 0.5;
	cfg->gridline.sz = 0.5;
	cfg->ticline.sz = 0.5;

	for (i = 0; ; i++) {
		rc = kplot_get_datacfg(c->views[c->view],
			i, &datas, &datasz);
		if (0 == rc)
			break;
		for (j = 0; j < datasz; j++) {
			sv = datas[j].line.sz;
			datas[j].line.sz = 1.0;
		}
	}

	cairo_set_source_rgb(cr, 1.0, 1.0, 1.0);
	cairo_rectangle(cr, 0.0, 0.0, w, h);
	cairo_fill(cr);
	kplot_draw(c->views[c->view], w, h, cr);
	cairo_destroy(cr);

	cfg->ticlabelfont.sz = svtic;
	cfg->axislabelfont.sz = svaxis;
	cfg->borderline.sz = svborder;
	cfg->gridline.sz = svgrid;
	cfg->ticline.sz = svticln;

	for (i = 0; ; i++) {
		rc = kplot_get_datacfg(c->views[c->view],
			i, &datas, &datasz);
		if (0 == rc)
			break;
		for (j = 0; j < datasz; j++)
			datas[j].line.sz = sv;
	}

	return(1);
}

int
save(const gchar *fname, struct curwin *cur)
{

	if (g_str_has_suffix(fname, ".pdf"))
		return(savepdf(fname, cur, SAVE_PDF));
	else if (g_str_has_suffix(fname, ".ps"))
		return(savepdf(fname, cur, SAVE_PS));
	else if (g_str_has_suffix(fname, ".eps"))
		return(savepdf(fname, cur, SAVE_EPS));
	else
		return(savepng(fname, cur));
}

int
saveconfig(const gchar *fname, const struct curwin *cur)
{
	FILE		*f;
	struct sim	*sim;
	GList		*l;

	g_debug("%p: Saving configuration: %s", cur, fname);

	if (NULL == (f = fopen(fname, "w"))) {
		g_debug("%s: %s", fname, strerror(errno));
		return(0);
	}

	for (l = cur->sims; NULL != l; l = g_list_next(l)) {
		sim = l->data;
		fprintf(f, "Name: %s\n", sim->name);
		fprintf(f, "Colour: #%.2x%.2x%.2x\n",
			(unsigned int)(cur->b->clrs[sim->colour].rgba[0] * 255),
			(unsigned int)(cur->b->clrs[sim->colour].rgba[1] * 255),
			(unsigned int)(cur->b->clrs[sim->colour].rgba[2] * 255));
		fprintf(f, "Function: %s\n", sim->func);
		fprintf(f, "Threads: %zu\n", sim->nprocs);
		fprintf(f, "Multiplier: %g(1 + %g lambda)\n",
			sim->alpha, sim->delta);
		fprintf(f, "Max generations: %zu\n", sim->stop);
		fprintf(f, "Migration: %g (%suniform)\n",
			sim->m, NULL != sim->ms ? "non-" : "");
		fprintf(f, "Incumbents: %zu, [%g,%g)\n",
			sim->dims, sim->xmin, sim->xmax);
		fprintf(f, "Rolling average window: %zu\n", sim->smoothing);
		switch (sim->maptop) {
		case (MAPTOP_RECORD):
			fprintf(f, "Map: record-based\n");
			break;
		case (MAPTOP_RAND):
			fprintf(f, "Map: random\n");
			break;
		case (MAPTOP_TORUS):
			fprintf(f, "Map: torus\n");
			break;
		default:
			abort();
		}
		switch (sim->migrant) {
		case (MAPMIGRANT_UNIFORM):
			fprintf(f, "Migration: uniform\n");
			break;
		case (MAPMIGRANT_DISTANCE):
			fprintf(f, "Migration: distance\n");
			break;
		case (MAPMIGRANT_NEAREST):
			fprintf(f, "Migration: nearest\n");
			break;
		case (MAPMIGRANT_TWONEAREST):
			fprintf(f, "Migration: two nearest\n");
			break;
		default:
			abort();
		}
		if (MAPINDEX_STRIPED == sim->mapindex)
			fprintf(f, "Mutant index case: striped\n");
		else
			fprintf(f, "Mutant index case: fixed (%zu)\n",
				sim->mapindexfix);
		if (MUTANTS_DISCRETE == sim->mutants)
			fprintf(f, "Mutants: %zu, [%g,%g)\n",
				sim->dims, sim->ymin, sim->ymax);
		else
			fprintf(f, "Mutants: N(sigma=%g), [%g,%g)\n",
				sim->mutantsigma, sim->ymin, sim->ymax);
		fprintf(f, "Islands: %zu (%zu islanders)\n",
			sim->islands, sim->totalpop);
		if (NULL != sim->pops)
			fprintf(f, "Island populations: non-uniform\n");
		else
			fprintf(f, "Island populations: %zu\n", sim->pop);
		fprintf(f, "Fit polynomial: %zu (%sweighted)\n",
			sim->fitpoly, 0 == sim->weighted ? "un" : "");

		fprintf(f, "\n");
	}

	fclose(f);
	return(1);
}

max_stars_repo_name: kristapsdz/bmigrate

max_issues_repo_name: kristapsdz/bmigrate

max_forks_repo_name: kristapsdz/bmigrate