// Supplementary to "Software Design ...", John A Robinson, Newnes, 2004,
//
// Program to visualize operation of simplex.
// Uses CLIP - Class library for image processing
//
// To generate pictures where successive moves of simplex are overlaid
// in increasing luminance, make changes as per *** comments. Such
// pictures were used in Software Design book.
//
#include "simplex.h"
#include "picture.h"
#include "estring.h"
picture_of_int *pin, *pout;
double picvalat(const double *dims) {
	double x = dims[0];
	double y = dims[1];
	double retval = (double) (*pin)[y][x];
	return retval;
	}
void triangle(int *x, int *y, int *starts, int *ends) {
	// Sort
	int i;
	for (i = 1; i < 3; i++) {
		if ((y[i] < y[0])||((y[i] == y[0])&&(x[i] < x[0]))) {
			int temp = x[0]; x[0] = x[i]; x[i] = temp;
			temp = y[0]; y[0] = y[i]; y[i] = temp;
			}
		}
	if ((y[1] > y[2])||((y[1] == y[2])&&(x[1] > x[2]))) {
		int temp = x[2]; x[2] = x[1]; x[1] = temp;
		temp = y[2]; y[2] = y[1]; y[1] = temp;
		}
	// Now x[0],y[0] contains lowest vals and x[2],y[2] highest vals
	// Trace down both sides
	int linenum = y[0];
	float xposstart = (float) x[0];
	float xposend = (float) x[0];
	if (y[1] == y[0])
		xposend = x[1];
	float xinc01, xinc02, xinc12;
	if (y[1] == y[0])
		xinc01 = x[1]-x[0];
	else
		xinc01 = (float)(x[1]-x[0])/(y[1]-y[0]);
	if (y[2] == y[0])
		xinc02 = x[2]-x[0];
	else
		xinc02 = (float)(x[2]-x[0])/(y[2]-y[0]);
	if (y[2] == y[1])
		xinc12 = x[2]-x[1];
	else
		xinc12 = (float)(x[2]-x[1])/(y[2]-y[1]);
	for (i = 0; i < y[2]-y[0]+1; i++) {
		if (xposend > xposstart) {
			starts[i] = (int) xposstart;
			ends[i] = (int) xposend;
			}
		else {
			starts[i] = (int) xposend;
			ends[i] = (int) xposstart;
			}
		if (i + y[0] == y[1])
			xposend = x[1];
		else if (i + y[0] > y[1])
			xposend += xinc12;
		else
			xposend += xinc01;
		if (i + y[0] == y[2])
			xposstart = x[2];
		else
			xposstart += xinc02;
		}
	}
void showsimplex(const simplex& s, const int paintval) {
	// *** Comment out next line if make a cumulative series of pictures
	*pout = *pin;
	int x[3], y[3];
	int startsarray[(*pin).nrows()];
	int endsarray[(*pin).nrows()];
	double *v;
	v = s.getv(0);
	x[0] = (int) v[0]; y[0] = (int) v[1];
	v = s.getv(1);
	x[1] = (int) v[0]; y[1] = (int) v[1];
	v = s.getv(2);
	x[2] = (int) v[0]; y[2] = (int) v[1];
	triangle(x,y,startsarray,endsarray);
	for (int i = 0; i < y[2]-y[0]+1; i++) {
		for (int j = startsarray[i]; j <= endsarray[i]; j++)
			(*pout)[i+y[0]][j] = paintval;
		}
	}

	
int main(int argc, char *argv[]) {
	if (argc < 2) {
		cerr << "Usage: vizsimplex function_as_picture xstart ystart\n";
		return -1;
		}
	picture_of_int in(argv[1]);
	in *= 1000;	// To get good variation when interpolated.
	pin = &in;
	// int paintval = in.min();
	// *** Uncomment previous line and comment next if making cumulative
	int paintval = in.max() + 20;
	simplex s(2, picvalat);
	s.bound(0,0,in.ncols()-1);
	s.bound(1,0,in.nrows()-1);
	if (argc > 2) {
		double startcoords[2];
		startcoords[0] = atof(argv[2]);
		startcoords[1] = atof(argv[3]);
		if (argc > 4) {
			double size = atof(argv[4]);
			s.initequ(startcoords,size);
			}
		else
			s.initat(startcoords);
		}
	else
		s.initat(); // Initialize at random position
	s.print();
	picture_of_int out(in);
	pout = &out;
	showsimplex(s,paintval);
	out.show("Simplex picture");
	out.inspect();
	int movetype;
	estring e("outpic000");
	while(movetype = s.move(4)) {
		// *** Uncomment next line if making a cumulative series
		// paintval += in.max()/40;
		showsimplex(s,paintval);
		out.reshow();
		// cout << "After move type " << movetype << endl;
		s.print();
		out.inspect();
		out.write(e.estr2str());
		++e;
		}
	cout << "Found minimum at\n";
	double *verts = s.getmin();
	cout << "f(" << verts[0] << "," << verts[1];
	cout << ") = " << s.getminval() << endl;
	// out.inspect();
	}