#include <stdlib.h>
#include "../pnheaders/matrix.h"
#include "../pnheaders/pns.h"
#include "spnbox.h"
#include "extendedmatrix.h"

static void ExtractConstraint(matrix* L, matrix* l, int* B, int* b, int Constraint);
static void RestoreConstraint(matrix* L, matrix* l, int* B, int b, int Constraint);

reduce_r reduce(matrix* L0, int* B0)
{
	reduce_r result;
	chkcons_r ChkResult;
	matrix l, L;
	int i, j, b, Rows;
	int *KeepRow, *B;
	
	/*Initialize the result to a null error return value.*/
	memset(&result, 0, sizeof(reduce_r));
	
	/*L and B are both required.*/
	if (! (L0 && B0))
	{
		merror(0, "reduce: One of L or B has not been given");
		return result;
	}
	
	/*Make a copy of the constraints that we can edit.*/
	AllocateMatrixType(2, &L, NumberOfRows(*L0), NumberOfColumns(*L0));
	CopyMatrix(L0, &L);
	AllocateMatrixType(2, &l, 1, NumberOfColumns(*L0));
	B = tcalloc(NumberOfRows(*L0), sizeof(int));
	memcpy(B, B0, NumberOfRows(*L0) * sizeof(int));
	
	/*Allocate space for a flag array to be used to indicate if a row should be
	removed.*/
	KeepRow = tcalloc(NumberOfRows(*L0), sizeof(int));
	
	/*The original matlab implementation removed constraint rows found to be
	redundant as it went along. In order to optimize for speed we will simply
	zero such constraint rows here. chkcons uses ipsolve to test for redundancy,
	and ipsolve will automatically ignore null rows.*/
	
	/*Iterate through the constraint rows. Test each one for redundancy using
	chkcons, and if it is redundant mark it for removal and zero it.*/
	Rows = 0;
	for (i = 0; i < NumberOfRows(*L0); i++)
	{
		if (is_verbose() >= VRB_REDUCE) printf("REDUCE: Constraint %3d\n", i);
		
		/*chkcons tests a given set of constraints against a single constraint.
		Copy the current constraint row into the row variables and then zero it
		in the constraint variables. If it turns out not to be redundant it will be
		restored.*/
		ExtractConstraint(&L, &l, B, &b, i);
		
		/*Check the constraint*/
		ChkResult = chkcons(&L, B, &l, b);
		
		/*chkcons will return a 0 in resCount if the constraint is redundant.*/
		if (ChkResult.resCount)
		{
			/*The constraint was not redundant. Restore it, mark it for keeping, and
			increment the count of rows kept.*/
			RestoreConstraint(&L, &l, B, b, i);
			KeepRow[i] = 1;
			Rows++;
		}
		DeallocateChkcons(&ChkResult);
		/*Otherwise the row is redundant. Leave it nulled out.*/
	}	
	
	/*Allocate space for the result and copy it.*/
	AllocateMatrixType(2, &result.Lf, Rows, NumberOfColumns(L));
	result.Bf = tcalloc(Rows, sizeof(int));
	result.indf = tcalloc(Rows, sizeof(int));
	j = 0;
	for (i = 0; i < NumberOfRows(L); i++)
	{
		if (KeepRow[i])
		{
			CopyBlock(1, NumberOfColumns(L), &L, i, 0, &result.Lf, j, 0);
			result.Bf[j] = B[i];
			/*Note the original indices of the rows we kept.*/
			result.indf[j] = i;
			j++;			
		}
	}
	/*If we removed any rows, set how to 1.*/
	if (Rows != NumberOfRows(L))
	{
		result.how = 1;
	}
	/*Free memory.*/
	DeallocateMatrix(&L);
	DeallocateMatrix(&l);
	free(B);
	free(KeepRow);
	
	return result;
}
/*******************************************************************************
ExtractConstraint nulls out l and swaps its single row with the given constraint
row from L. It swaps out the proper element of B and b as well. This is in
preparation for redundancy tests with chkcons.*/
void ExtractConstraint(matrix* L, matrix* l, int* B, int* b, int Constraint)
{
	int i;
	MakeZeroRow(l, 0);
	SwapRows(L, l, Constraint, 0);
	*b = B[Constraint];
	B[Constraint] = 0;
}

/*******************************************************************************
RestoreConstraint reverses the operations carried out in ExtractConstraint. Used
when the given constraint is determined to be non-redundant.*/
void RestoreConstraint(matrix* L, matrix* l, int* B, int b, int Constraint)
{
	int i;
	SwapRows(L, l, Constraint, 0);
	B[Constraint] = b;
}
