/*

 * int_lin_eq.c

 *
 *  Created on: 22/04/2017

 *      Author: Pedro

 */

#ifndef __OPENCL_VERSION__

#include <stddef.h>
#include <stdio.h>
#include <math.h>

#include "int_lin_eq.h"

#include "../bitmaps.h"
#include "../config.h"
#include "../variables.h"

#endif

#include "../kernels/cl_aux_functions.h"
#if CL_D_TYPE == CL_BITMAP
#include "../kernels/cl_bitmaps.h"
#elif CL_D_TYPE == CL_INTERVAL
#include  "../kernels/cl_intervals.h"
#endif
#include "../kernels/cl_constraints.h"
#include "../kernels/cl_variables.h"
#include "../kernels/cl_ttl.h"

#ifndef __OPENCL_VERSION__

/*

 * Creates a new constraint of the int_lin_eq type and return the constraint ID

 * K · Y = c
 * K - constant unsigned integers for this constraint
 * Y_ids - IDs of the variables constrained by this constraint
 * n - number of constants (or variables)
 * c - result of the equation
 */

unsigned int c_int_lin_eq(int *K, unsigned int *Y_ids, unsigned int n, int c) {

	unsigned int i;

	// set to include in kernel compilation

	USE_CS[INT_LIN_EQ] = 1;
	USE_NON_CS_REIFI[INT_LIN_EQ] = 1;

	REV = 1;

	unsigned int *c_vs = malloc(n * sizeof(unsigned int));

	for (i = 0; i < n; i++) {
		c_vs[i] = Y_ids[i];
	}

	// creates a new generic constraint
	unsigned int c_id = c_new(c_vs, n, K, n, -1);

	// pointers to this type of constraint functions

	CS[c_id].kind = INT_LIN_EQ;
	CS[c_id].check_sol_f = &int_lin_eq_check;

	CS[c_id].constant_val = c;

	free(c_vs);

	return c_id;
}

/*

 * Creates a new reified constraint of the int_lin_eq type and return the constraint ID

 * K · Y = c
 * K - constant unsigned integers for this constraint
 * Y_ids - IDs of the variables constrained by this constraint
 * n - number of constants (or variables)
 * c - result of the equation
 * reif_v_id - ID of the reification variable
 */

unsigned int c_int_lin_eq_reif(int *K, unsigned int *Y_ids, unsigned int n, int c, int reif_v_id) {

	unsigned int i;

	if (VS[reif_v_id].max > 1) {
		v_del_gt(&VS[reif_v_id], 1);

		if (VS[reif_v_id].n_vals == 0) {

			printf("\nConstraint INT_LIN_EQ_REIF makes model inconsistent at creation. No solution found.\n");

#if defined(WIN32) || defined(_WIN32) || defined(__WIN32) && !defined(__CYGWIN__)
			printf("\nPress any key to exit\n");
			int a = getchar();
#endif

			exit(0);

		}
	}

	// set to include in kernel compilation

	USE_CS[INT_LIN_EQ] = 1;
	USE_CS_REIFI[INT_LIN_EQ] = 1;

	REV = 1;

	unsigned int *c_vs = malloc(n * sizeof(unsigned int));

	for (i = 0; i < n; i++) {
		c_vs[i] = Y_ids[i];
	}

	// creates a new generic constraint
	unsigned int c_id = c_new(c_vs, n, K, n, reif_v_id);

	// pointers to this type of constraint functions

	CS[c_id].kind = INT_LIN_EQ;
	CS[c_id].check_sol_f = &int_lin_eq_check;

	CS[c_id].constant_val = c;

	free(c_vs);

	return c_id;
}

/*

 * Return true if the int_lin_eq constraint is respected or false if not

 * K · Y = c
 * c - constraint to check if is respected
 * explored - if the CSP was already explored, which mean that all the variables must already be singletons
 * */

bool int_lin_eq_check(constr *c, bool explored) {

	unsigned int n = c->n_c_consts;		// number of constants and variables constrained by this constraint

	int *K = c->c_consts;	// constants constrained by this constraint
	var **Y = c->c_vs;	// variables constrained by this constraint

	int equat_result = 0;
	unsigned int i;

	if (!explored) {
		for (i = 0; i < c->n_c_vs; i++) {
			if (c->c_vs[i]->n_vals > 1) {
				return false;

			}

		}
	}

	if (c->reified && VS[c->reif_v_id].n_vals > 1) {
		if (explored) {
			fprintf(stderr, "\nError: Reification variable of constraint INT_LIN_EQ_REIF (%d) has 2 values.\n", c->c_id);

			return false;
		}
	}

	for (i = 0; i < n; i++) {
		equat_result += K[i] * Y[i]->min;
	}

	if (((!c->reified || (c->reified && VS[c->reif_v_id].min == 1)) && equat_result != c->constant_val)
			|| (c->reified && VS[c->reif_v_id].min == 0 && equat_result == c->constant_val)) {

		if (explored) {

			if (c->reified) {
				fprintf(stderr, "\nError: Constraint INT_LIN_EQ_REIF (%d) not respected:\n", c->c_id);
				fprintf(stderr, "Reif ID=%u -> minimum=%u, maximum=%u, number of values=%u\n\n", c->reif_v_id, b_get_min_val(&VS[c->reif_v_id].domain_b),
						b_get_max_val(&VS[c->reif_v_id].domain_b), b_cnt_vals(&VS[c->reif_v_id].domain_b));

			} else {
				fprintf(stderr, "\nError: Constraint INT_LIN_EQ (%d) not respected:\n", c->c_id);
			}

			for (i = 0; i < c->n_c_vs; i++) {
				fprintf(stderr, "Variable ID=%u -> minimum=%u, maximum=%u, number of values=%u\n\n", c->c_vs[i]->v_id, b_get_min_val(&c->c_vs[i]->domain_b),
						b_get_max_val(&c->c_vs[i]->domain_b), b_cnt_vals(&c->c_vs[i]->domain_b));
			}
		}
		return false;
	}

	return true;
}

#endif

#if CS_INT_LIN_EQ == 1

/*

 * Propagate the domain of the variable with the ID prop_v_id through all the other variables on the same c_numb ID int_lin_eq constraint

 * K · Y = c
 * vs_per_c_idx - vector with all constrained variables ID per constraint, per constraint ID order

 * c_consts - constant values used by this constraint

 * vs_prop_ - all CSP variables with current step values

 * current_cs - constraint that should be propagated for the variable with prop_v_id ID
 * vs_id_to_prop_ - circular vector with the ids of the variables to propagate

 * prop_ok - will be set to 1 or 0 if the constraint is respected or not
 * terms_mem - auxiliary buffer

 */

CUDA_FUNC void int_lin_eq_prop( CL_INTS_MEM int *vs_per_c_idx, CL_INTS_MEM int *c_consts, CL_MEMORY VARS_PROP *vs_prop_, CL_CS_MEM cl_constr *current_cs,
		CL_MEMORY unsigned short *vs_id_to_prop_, bool *prop_ok, __global int *terms_mem CS_IGNORE_FUNC TTL_CTR) {

	int terms = current_cs->n_c_consts;	// number of constants and variables constrained by this constraint

	CL_INTS_MEM int *K = c_consts;	// constants constrained by this constraint

	int equat_result = current_cs->constant_val;

	__global int *mins = terms_mem;
	__global int *maxs = &terms_mem[current_cs->n_c_consts];

	int y_id;

	int vl, vh;
	int min, max;

	int ymin, ymax;

	bool changed = 0;
	int c;
	int i;

	min = 0;
	max = 0;
	for (i = 0; i < terms; i++) {
		CHECK_TTL(ttl_ctr, 69)
		y_id = vs_per_c_idx[i];

		if (K[i] != 0) {

			if (K[i] > 0) {
				vl = mins[i] = V_MIN(vs_prop_[y_id]);
				vh = maxs[i] = V_MAX(vs_prop_[y_id]);
			} else {
				vl = maxs[i] = V_MAX(vs_prop_[y_id]);
				vh = mins[i] = V_MIN(vs_prop_[y_id]);
			}

			min += K[i] * vl;
			max += K[i] * vh;
		}

	}

	if (min > equat_result || max < equat_result) {
		*prop_ok = 0;
		return;
	}

	if (min == max) {

#if CL_CS_IGNORE
		cs_ignore[current_cs->c_id] = 1;
#endif
		return;
	}

	if (min == equat_result) {
		for (i = 0; i < terms; i++) {
			CHECK_TTL(ttl_ctr, 70)
			y_id = vs_per_c_idx[i];

			c = K[i];

			if (V_N_VALS(vs_prop_[y_id]) > 1 && c != 0) {

				if (c > 0) {

					cl_v_del_all_except_val_no_tests_m(&vs_prop_[y_id], mins[i] TTL_CTR_V);

					v_add_to_prop(vs_id_to_prop_, vs_prop_, y_id);

				} else {
					cl_v_del_all_except_val_no_tests_m(&vs_prop_[y_id], maxs[i] TTL_CTR_V);

					v_add_to_prop(vs_id_to_prop_, vs_prop_, y_id);
				}
			}
		}

#if CL_CS_IGNORE
		cs_ignore[current_cs->c_id] = 1;
#endif

		return;
	}

	if (max == equat_result) {

		for (i = 0; i < terms; i++) {

			CHECK_TTL(ttl_ctr, 71)
			y_id = vs_per_c_idx[i];

			c = K[i];

			if (V_N_VALS(vs_prop_[y_id]) > 1 && c != 0) {

				if (c > 0) {

					cl_v_del_all_except_val_no_tests_m(&vs_prop_[y_id], maxs[i] TTL_CTR_V);

					v_add_to_prop(vs_id_to_prop_, vs_prop_, y_id);

				} else {
					cl_v_del_all_except_val_no_tests_m(&vs_prop_[y_id], mins[i] TTL_CTR_V);

					v_add_to_prop(vs_id_to_prop_, vs_prop_, y_id);
				}
			}
		}

#if CL_CS_IGNORE
		cs_ignore[current_cs->c_id] = 1;
#endif
		return;
	}

#if CL_USE_BOOLEAN_VS
	if (current_cs->boolean) {

		for (i = 0; i < terms; i++) {
			CHECK_TTL(ttl_ctr, 72)
			y_id = vs_per_c_idx[i];
			c = K[i];

			if (V_N_VALS(vs_prop_[y_id]) > 1 && c != 0) {

				if ((c > 0 && min + c > equat_result) || (c < 0 && max + c < equat_result)) {

					cl_v_bool_del_val_m(&vs_prop_[y_id], 1 TTL_CTR_V);
					v_add_to_prop(vs_id_to_prop_, vs_prop_, y_id);

				} else if ((c > 0 && max - c < equat_result) || (c < 0 && min - c > equat_result)) {

					cl_v_bool_del_val_m(&vs_prop_[y_id], 0 TTL_CTR_V);
					v_add_to_prop(vs_id_to_prop_, vs_prop_, y_id);
				}
			}
		}
		return;
	}
#endif

	if (max > equat_result) {

		for (i = 0; i < terms; i++) {

			CHECK_TTL(ttl_ctr, 72)
			y_id = vs_per_c_idx[i];

			c = K[i];

			if (V_N_VALS(vs_prop_[y_id]) > 1 && c != 0) {

				ymin = mins[i];
				ymax = maxs[i];

				if (c > 0) {

					if ((ymax - ymin) * c > equat_result - min) {
						ymax = (equat_result - min) / c + ymin;

						cl_v_del_gt_m(&changed, &vs_prop_[y_id], ymax TTL_CTR_V);

					}

				} else {
					if ((ymin - ymax) * c > equat_result - min) {
						ymin = convert_int (ceil((equat_result - min * 1.0) / c + ymax));

						cl_v_del_lt_m(&changed, &vs_prop_[y_id], ymin TTL_CTR_V);

					}
				}

				if (changed) {
					if (V_IS_EMPTY(vs_prop_[y_id])) {
						*prop_ok = 0;
						return;
					}
					v_add_to_prop(vs_id_to_prop_, vs_prop_, y_id);
				}

			}
		}
	}

	if (min < equat_result) {

		for (i = 0; i < terms; i++) {

			CHECK_TTL(ttl_ctr, 73)
			y_id = vs_per_c_idx[i];

			c = K[i];

			if (V_N_VALS(vs_prop_[y_id]) > 1 && c != 0) {

				ymin = mins[i];
				ymax = maxs[i];

				if (c > 0) {

					if ((ymax - ymin) * c > max - equat_result) {
						ymin = convert_int (ceil((equat_result - max * 1.0) / c + ymax));

						cl_v_del_lt_m(&changed, &vs_prop_[y_id], ymin TTL_CTR_V);

					}

				} else {
					if ((ymax - ymin) * c < equat_result - max) {
						ymax = (equat_result - max) / c + ymin;

						cl_v_del_gt_m(&changed, &vs_prop_[y_id], ymax TTL_CTR_V);

					}
				}

				if (changed) {
					if (V_IS_EMPTY(vs_prop_[y_id])) {
						*prop_ok = 0;
						return;
					}
					v_add_to_prop(vs_id_to_prop_, vs_prop_, y_id);
				}

			}
		}
	}
}

#if CS_R_INT_LIN_EQ == 1

/*

 * Validate int_lin_eq constraint to be normally propagated, when reified

 * K · Y = c
 * vs_per_c_idx - vector with all constrained variables ID per constraint, per constraint ID order

 * c_consts - constant values used by this constraint

 * vs_prop_ - all CSP variables with current step values
 * current_cs - constraint that should be propagated for the variable with prop_v_id ID
 * vs_id_to_prop_ - circular vector with the ids of the variables to propagate
 */

CUDA_FUNC void int_lin_eq_reif( CL_INTS_MEM int *vs_per_c_idx, CL_INTS_MEM int *c_consts, CL_MEMORY VARS_PROP *vs_prop_, CL_CS_MEM cl_constr *current_cs,
		CL_MEMORY unsigned short *vs_id_to_prop_ CS_IGNORE_FUNC TTL_CTR) {

	int terms = current_cs->n_c_consts;	// number of constants and variables constrained by this constraint

	CL_INTS_MEM int *K = c_consts;	// constants constrained by this constraint

	int equat_result = current_cs->constant_val;
	int y_id;
	int vl, vh;
	int min, max;
	int i;

	min = 0;
	max = 0;

	for (i = 0; i < terms; i++) {

		CHECK_TTL(ttl_ctr, 74)
		y_id = vs_per_c_idx[i];
		if (K[i] > 0) {
			vl = V_MIN(vs_prop_[y_id]);
			vh = V_MAX(vs_prop_[y_id]);
		} else {
			vl = V_MAX(vs_prop_[y_id]);
			vh = V_MIN(vs_prop_[y_id]);
		}

		min += K[i] * vl;
		max += K[i] * vh;
	}

	if (min > equat_result || max < equat_result) {
		cl_v_bool_del_val_m(&vs_prop_[current_cs->reif_var_id], 1 TTL_CTR_V);

		v_add_to_prop(vs_id_to_prop_, vs_prop_, convert_int (current_cs->reif_var_id));

#if CL_CS_IGNORE

		cs_ignore[current_cs->c_id] = 1;

#endif
		return;
	}

	// constraint already fixed
	if (min == max && min == equat_result) {
		cl_v_bool_del_val_m(&vs_prop_[current_cs->reif_var_id], 0 TTL_CTR_V);

		v_add_to_prop(vs_id_to_prop_, vs_prop_, convert_int (current_cs->reif_var_id));

#if CL_CS_IGNORE

		cs_ignore[current_cs->c_id] = 1;

#endif
		return;
	}
}

/*

 * Propagate the domain of the variable with the ID prop_v_id through all the other variables on the same c_numb ID int_lin_eq opposite constraint

 * K · Y != c
 * vs_per_c_idx - vector with all constrained variables ID per constraint, per constraint ID order

 * c_consts - constant values used by this constraint

 * vs_prop_ - all CSP variables with current step values

 * current_cs - constraint that should be propagated for the variable with prop_v_id ID
 * vs_id_to_prop_ - circular vector with the ids of the variables to propagate

 * prop_ok - will be set to 1 or 0 if the constraint is respected or not
 * terms_mem - auxiliary buffer

 */

CUDA_FUNC void int_lin_eq_prop_opposite( CL_INTS_MEM int *vs_per_c_idx, CL_INTS_MEM int *c_consts, CL_MEMORY VARS_PROP *vs_prop_,
		CL_CS_MEM cl_constr *current_cs, CL_MEMORY unsigned short *vs_id_to_prop_, bool *prop_ok, __global int *terms_mem CS_IGNORE_FUNC TTL_CTR) {

	int terms = current_cs->n_c_consts;	// number of constants and variables constrained by this constraint

	CL_INTS_MEM int *K = c_consts;	// constants constrained by this constraint

	int equat_result = current_cs->constant_val;
	int y_id;

	__global int *mins = terms_mem;

	int vl, vh;
	int min, max;
	bool changed = 0;
	int not_singl = 0;
	int not_singl_idx = 0;
	int not_singl_id = -1;
	int val_to_rem;
	int sum = 0;
	int i;

	min = 0;
	max = 0;

	for (i = 0; i < terms; i++) {

		CHECK_TTL(ttl_ctr, 69)
		y_id = vs_per_c_idx[i];
		if (K[i] > 0) {
			vl = mins[i] = V_MIN(vs_prop_[y_id]);
			vh = V_MAX(vs_prop_[y_id]);
		} else {
			vl = V_MAX(vs_prop_[y_id]);
			vh = mins[i] = V_MIN(vs_prop_[y_id]);
		}

		if (K[i] != 0) {
			if (vl != vh) {
				not_singl_idx = i;
				not_singl_id = y_id;
				not_singl++;

			} else {
				sum += K[i] * mins[i];
			}
		}

		min += K[i] * vl;
		max += K[i] * vh;
	}
	sum -= equat_result;

	if (min == max && min == equat_result) {
		*prop_ok = 0;
		return;
	}

	if (min > equat_result || max < equat_result) {

#if CL_CS_IGNORE
		cs_ignore[current_cs->c_id] = 1;
#endif
		return;
	}

	// if all but one variable are already singleton, remove the only value from the one that is not singleton that would lead to equality
	if (not_singl == 1) {
		val_to_rem = (-1) * (sum / K[not_singl_idx]);

		if (sum % K[not_singl_idx] == 0 && val_to_rem >= 0) {

			cl_v_del_val_m(&changed, &vs_prop_[not_singl_id], val_to_rem TTL_CTR_V);

			if (changed) {
				v_add_to_prop(vs_id_to_prop_, vs_prop_, not_singl_id);
			}
		}

#if CL_CS_IGNORE
		cs_ignore[current_cs->c_id] = 1;
#endif
	}
}

#endif

/*
 * Decides the propagator to call for this constraint
 * vs_per_c_idx - vector with all constrained variables ID per constraint, per constraint ID order
 * c_consts - constant values used by this constraint
 * vs_prop_ - all CSP variables with current step values
 * current_cs - constraint that should be propagated for the variable with prop_v_id ID
 * vs_id_to_prop_ - circular vector with the ids of the variables to propagate
 * prop_ok - will be set to 1 or 0 if the constraint is respected or not
 * terms_mem - auxiliary buffer
 */
CUDA_FUNC void int_lin_eq_propagate( CL_INTS_MEM int *vs_per_c_idx, CL_INTS_MEM int *c_consts, CL_MEMORY VARS_PROP *vs_prop_,
		CL_CS_MEM cl_constr *current_cs, CL_MEMORY unsigned short *vs_id_to_prop_, bool *prop_ok, __global int *terms_mem
		PROPAGATED_FUNC CS_IGNORE_FUNC TTL_CTR) {

#if CS_R_INT_LIN_EQ == 0
	int_lin_eq_prop(vs_per_c_idx, c_consts, vs_prop_, current_cs, vs_id_to_prop_, prop_ok, terms_mem CS_IGNORE_CALL TTL_CTR_V);

#if CL_STATS == 1
	*propagated = true;
#endif

#elif CS_R_INT_LIN_EQ == 1

	if (current_cs->reified == 1) {
		if (V_N_VALS(vs_prop_[current_cs->reif_var_id]) > 1) {

			int_lin_eq_reif(vs_per_c_idx, c_consts, vs_prop_, current_cs, vs_id_to_prop_ CS_IGNORE_CALL TTL_CTR_V);

		} else {
			if (V_MIN(vs_prop_[current_cs->reif_var_id]) == 1) {

				int_lin_eq_prop(vs_per_c_idx, c_consts, vs_prop_, current_cs, vs_id_to_prop_, prop_ok, terms_mem CS_IGNORE_CALL TTL_CTR_V);

			} else {

				int_lin_eq_prop_opposite(vs_per_c_idx, c_consts, vs_prop_, current_cs, vs_id_to_prop_, prop_ok, terms_mem CS_IGNORE_CALL TTL_CTR_V);

			}
#if CL_STATS == 1
			*propagated = true;
#endif
		}
	} else {

		int_lin_eq_prop(vs_per_c_idx, c_consts, vs_prop_, current_cs, vs_id_to_prop_, prop_ok, terms_mem CS_IGNORE_CALL TTL_CTR_V);

#if CL_STATS == 1
		*propagated = true;
#endif
	}
#endif
}

#endif