/*
 * International Color Consortium color transform expanded support
 *
 * Author:  Graeme W. Gill
 * Date:    2/7/00
 * Version: 1.00
 *
 * Copyright 2000 - 2007 Graeme W. Gill
 * All rights reserved.
 * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
 * see the License.txt file for licencing details.
 *
 * Based on the xlut.c code.
 */

/*
 * TTBD:
 *
 *		Need to use this for B2A tables rather than inverting
 *		A2B curves. Need to add grid sizing to cover just gamut range
 *      (including level axis gamut, but watch out for devices that
 *      have values below the black point or above the white point),
 *      3x3 matrix optimization, and white point to grid node mapping for B2A.
 *		Currently code assumes output is always PCS ?? - would need to fix for opposite.
 *
 *      Currently the Lab A2B output tables are adjusted for ab symetry
 *		to make the B2A white point land on a grid point, given that
 *		the icc code forces a symetric ab range. This only works
 *      because the B2A is using the same per channel curves.
 *		Done properly, it should be possible to know where grid
 *		points land within the range, and to modify the B2A input curves
 *		to make the white point land on a grid point.
 *		(For B2A the pseudo least squares adjustment needs to be turned
 *		off for that grid point too.)
 *
 *		Note that one quandry is that the curve fitting doesn't
 *		fit well when the input data has an offset and/or plateaus.
 */

/*
 * This module provides curve and matrix fitting functionality used to
 * create per channel input and/or output curves for clut profiles,
 * and optionally creates the rspl to fit the data as well.
 * 
 * The approach used is to initialy creating input and output shaper
 * curves that minimize the overall delta E of the test point set,
 * when a linear matrix is substituted for a clut is used.
 * (This is the same as pre-V0.70 approach. )
 *
 * The residual error between the test point set and the shaper/matrix/shaper
 * model is then computed, and mapped against each input channel, and
 * a positioning mapping curve created that aims to map rspl grid
 * locations more densly where residual errors are high, and more
 * sparsely where they are low.
 * The input shaper and positioning curves are then combined together,
 * so that the positioning curve determines which input values map
 * to rspl grid points, and the shaper curve determines the mapping
 * between the grid points. The width of the grid cells is computed
 * in the shaper mapped input space, and then fed to the rspl fitting
 * code, so that its smoothness evaluation function can take account
 * of the non-uniform spacing of the grid points.
 */

#include <sys/types.h>
#include <string.h>
#include <ctype.h>
#ifdef __sun
#include <unistd.h>
#endif
#include "copyright.h"
#include "aconfig.h"
#include "numlib.h"
#include "icc.h"
#include "rspl.h"
#include "plot.h"
#include "xicc.h"
#include "xfit.h"
#include "sort.h"

#undef USE_XYZ_Y2LCURVE     /* [Und] Use underlying L* curve for XYZ encoding */
							/* This seems to work badly, even with high smoothness. Why ? */
							/* It does speed up 1D lut creation though. */

#undef DEBUG 			/* Debug information */
#undef DEBUG_PROGRESS 	/* Show powell progress */
#undef DEBUG_PLOT 		/* Plot in & out curves */
#undef SPECIAL_FORCE	/* Check rspl nodes against linear XYZ model */
#undef SPECIAL_FORCE_GAMMA		/* Force correct gamma shaper curves */
#undef SPECIAL_TEST_GAMMA		/* Use gamma on reference model */
#undef SPECIAL_TEST_LAB	

#define RSPLFLAGS (0 /* | RSPL_2PASSSMTH */ /* | RSPL_EXTRAFIT2 */)

#undef EXTEND_GRID			/* [Undef] Use extended RSPL grid around interpolation */
#define EXTEND_GRID_BYN 2	/* Rows to extend grid. */

#undef NODDV				/* Use slow non d/dv powell else use conjgrad */
#define CURVEPOW 1.0    	/* Power to raise deltaE squared to in setting in/out curves */
							/* This provides a means of punishing high maximum errors. */

#define POWTOL1 1e-3		/* Shaper Powell optimiser tollerance for first passes */
#define MAXITS1 1000		/* Shaper number of itterations for first passes */
#define POWTOL 1e-5			/* Shaper Powell optimiser tollerance in delta E squared ^ CURVEPOW */
#define MAXITS 4000			/* Shaper number of itterations before giving up */
#define PDDEL  1e-6			/* Fake partial derivative del */

/* Weights for shaper in/out curve parameters, to minimise unconstrained "wiggles" */
#define SHAPE_WEIGHT	1.0		/* Overal shaper weight contribution - err on side of smoothness */
#define SHAPE_HW01		0.002	/* 0 & 1 harmonic weights */
#define SHAPE_HBREAK    4		/* Harmonic that has HWBR */
#define SHAPE_HWBR  	20.0	/* Base weight of harmonics HBREAK up */
#define SHAPE_HWINC  	60.0	/* Increase in weight for each harmonic above HWBR */

/* Weights for the positioning curve parameters */
#define PSHAPE_MINE 0.02		/* Minum background residual error level */
#define PSHAPE_DIST 1.0			/* Agressivness of grid distribution */


/* - - - - - - - - - - - - - - - - - */

/* Extra non-linearity used as base for XYZ output curves. */
/* This makes the XYZ grid values more perceptual, and asks less */
/* of the automatically created output curve shape. */
/* (We assume XYZ is in 0..1 scale */

/* Transfer function with offset and scale + Y2L curve */
static double icxSTransFuncY2L(
double *v,			/* Pointer to first parameter */
int    luord,		/* Number of parameters */
double vv,			/* Source of value */
double min,			/* Scale values */
double max
) {
	max -= min;

	vv = (vv - min)/max;

#ifdef USE_XYZ_Y2LCURVE
	if (vv > 0.008856451586)
		vv = 1.16 * pow(vv,1.0/3.0) - 0.16;
	else
		vv = 9.032962896 * vv;
#endif

	vv = icxTransFunc(v, luord, vv);

	vv = (vv * max) + min;
	return vv;
}

/* Inverse Transfer function with offset and scale + Y2L */
static double icxInvSTransFuncY2L(
double *v,			/* Pointer to first parameter */
int    luord,		/* Number of parameters */
double vv,			/* Source of value */
double min,			/* Scale values */
double max
) {
	max -= min;

	vv = (vv - min)/max;
	vv = icxInvTransFunc(v, luord, vv);

#ifdef USE_XYZ_Y2LCURVE
	if (vv > 0.08)
		vv = pow((vv + 0.16)/1.16, 3.0);
	else
		vv = vv/9.032962896;
#endif

	vv = (vv * max) + min;

	return vv;
}

/* Transfer function with offset and scale, and */
/* partial derivative with respect to the */
/* parameters and the input value. */
static double icxdpdiSTransFuncY2L(
double *v,			/* Pointer to first parameter */
double *dv,			/* Return derivative wrt each parameter */
double *pdin,		/* Return derivative wrt source value */
int    luord,		/* Number of parameters */
double vv,			/* Source of value */
double min,			/* Scale values */
double max
) {
	int i;
	double idv = 1.0;
	max -= min;

#ifdef USE_XYZ_Y2LCURVE
	if (vv > 0.008856451586) {
		vv = 1.16 * pow(vv,1.0/3.0) - 0.16;
		idv = 1.16 / (3.0 * pow(vv, 2.0/3.0));
	} else {
		vv = 9.032962896 * vv;
		idv = 9.032962896; 
	}
#endif

	vv = (vv - min)/max;

	vv = icxdpdiTransFunc(v, dv, pdin, luord, vv);

	*pdin *= idv;		/* Account for input multiplier */

	vv = (vv * max) + min;

	for (i = 0; i < luord; i++) {
		dv[i] *= max;
	}
	return vv;
}

/* - - - - - - - - - - - - - - - - - */

#ifdef DEBUG
static void dump_xfit(xfit *p);
#endif

#ifdef DEBUG_PLOT	/* Not currently used in runtime code*/

/* Lookup a value though an input position curve */
static double xfit_poscurve(xfit *p, double in, int chan) {
	double rv = in;
	if (p->tcomb & oc_p)
		rv = icxSTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], rv,  
		                       p->in_min[chan], p->in_max[chan]);
	return rv;
}

/* Lookup di values though the input position curves */
static void xfit_poscurves(xfit *p, double *out, double *in) {
	int e;
	
	for (e = 0; e < p->di; e++) {
		double val = in[e];
		if (p->tcomb & oc_i)
			val = icxSTransFunc(p->v + p->pos_offs[e], p->iluord[e], val,  
			                       p->in_min[e], p->in_max[e]);
		out[e] = val;
	}
}

/* Inverse Lookup a value though an input position curve */
static double xfit_invposcurve(xfit *p, double in, int chan) {
	double rv = in;
	if (p->tcomb & oc_i)
		rv = icxInvSTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], rv,  
		                       p->in_min[chan], p->in_max[chan]);
	return rv;
}

/* Inverse Lookup di values though input position curves */
static void xfit_invposcurves(xfit *p, double *out, double *in) {
	int e;
	
	for (e = 0; e < p->di; e++) {
		double val = in[e];
		if (p->tcomb & oc_i)
			val = icxInvSTransFunc(p->v + p->pos_offs[e], p->iluord[e], val,  
			                       p->in_min[e], p->in_max[e]);
		out[e] = val;
	}
}
#endif /* DEBUG_PLOT */

/* - - - - - - - - - - - - - - - - - */
/* Lookup a value though input shape curve */
static double xfit_shpcurve(xfit *p, double in, int chan) {
	double rv = in;

	if (p->tcomb & oc_i) {
#ifdef SPECIAL_FORCE_GAMMA	
		double gam;
		if (chan == 0)
			gam = 1.9;
		else if (chan == 1)
			gam = 2.0;
		else if (chan == 2)
			gam = 2.1;
		else
			gam = 1.0;
		rv = pow(in, gam);
#else
		rv = icxSTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], rv,  
			                       p->in_min[chan], p->in_max[chan]);
#endif
	}
	return rv;
}

#ifdef NEVER	/* Not currently used */
/* Lookup a value though shape curves */
static void xfit_shpcurves(xfit *p, double *out, double *in) {
	int e;
	
	for (e = 0; e < p->di; e++) {
		double val = in[e];
		if (p->tcomb & oc_i)
			val = icxSTransFunc(p->v + p->shp_offs[e], p->iluord[e], val,  
			                       p->in_min[e], p->in_max[e]);
		out[e] = val;
	}
}
#endif /* NEVER */

/* Inverse Lookup a value though a shape curve */
static double xfit_invshpcurve(xfit *p, double in, int chan) {
	double rv = in;

	if (p->tcomb & oc_i) {
#ifdef SPECIAL_FORCE_GAMMA	
		double gam;
		if (chan == 0)
			gam = 1.9;
		else if (chan == 1)
			gam = 2.0;
		else if (chan == 2)
			gam = 2.1;
		else
			gam = 1.0;
		rv = pow(rv, 1.0/gam);
#else
		rv = icxInvSTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], rv,  
		                       p->in_min[chan], p->in_max[chan]);
#endif
	}
	return rv;
}

#ifdef NEVER	/* Not currently used */
/* Inverse Lookup a value though shape curves */
static void xfit_invshpcurves(xfit *p, double *out, double *in) {
	int e;
	
	for (e = 0; e < p->di; e++) {
		double val = in[e];
		if (p->tcomb & oc_i)
			val = icxInvSTransFunc(p->v + p->shp_offs[e], p->iluord[e], val,  
			                       p->in_min[e], p->in_max[e]);
		out[e] = val;
	}
}
#endif /* NEVER */

/* - - - - - - - - - - - - - - - - - */

/* Lookup values through the shaper/matrix/shaper model */
static void xfit_shmatsh(xfit *p, double *out, double *in) {
	double tin[MXDI];
	int e, f;
	
	for (e = 0; e < p->di; e++)
		tin[e] = icxSTransFunc(p->v + p->shp_offs[e], p->iluord[e], in[e],  
			                       p->in_min[e], p->in_max[e]);

	icxCubeInterp(p->v + p->mat_off, p->fdi, p->di, out, tin);

	if (p->flags & XFIT_OUT_LAB) {
		for (f = 0; f < p->fdi; f++)
			out[f] = icxSTransFunc(p->v + p->out_offs[f], p->oluord[f], out[f],
				                       p->out_min[f], p->out_max[f]);
	} else {
		for (f = 0; f < p->fdi; f++)
			out[f] = icxSTransFuncY2L(p->v + p->out_offs[f], p->oluord[f], out[f],
				                       p->out_min[f], p->out_max[f]);
	}
}

/* - - - - - - - - - - - - - - - - - - - - */
/* Combined input positioning & shaper transfer curve functions */

//int db = 0;

/* Lookup a value though the input positioning and shaper curves */
static double xfit_inpscurve(xfit *p, double in, int chan) {
	double rv;
	/* Just shaper curve */
	if ((p->tcomb & oc_ip) == oc_i) {
		rv = icxSTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], in,  
			                       p->in_min[chan], p->in_max[chan]);

	/* shaper and positioning */
	} else if ((p->tcomb & oc_ip) == oc_ip) {
		double nin, npind;		/* normalized in, normalized position in' */
		int six;				/* Span index */
		double npind0, npind1;	/* normalized position in' span values */
		double npin0, npin1;	/* normalized position in span values */
		double nsind0, nsind1;	/* normalized shaper in' span values */
		double nsind;			/* normalised shaper in' value */

		/* Normalize */
		nin = (in - p->in_min[chan])/(p->in_max[chan] - p->in_min[chan]);

//if (db) printf("\n~1 inpscurve: cha %d, input value %f, norm %f\n",chan,in,nin);

		/* Locate the span the input point will be in after positioning lookup */
		npind = icxTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], nin); 

		/* Quantize position space value to grid */
		six = (int)floor(npind * (p->gres[chan]-1.0));
		if (six > (p->gres[chan]-2))
			six = (p->gres[chan]-2);

		/* Compute span position position in' values */
		npind0 = six / (p->gres[chan]-1.0);
		npind1 = (six + 1.0) / (p->gres[chan]-1.0);

//if (db) printf("~1 npind %f, six %d, npind0 %f, npind1 %f\n",npind,six,npind0,npind1);

		/* Compute span in values */
		npin0 = icxInvTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], npind0); 
		npin1 = icxInvTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], npind1); 
	
//if (db) printf("~1 npin0 %f, npin1 %f\n",npin0,npin1);

		/* Compute shaper space values of in' and spane */
#ifdef NEVER
		nsind = icxTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], nin); 
		nsind0 = icxTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], npin0); 
		nsind1 = icxTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], npin1); 
#else
		nsind = xfit_shpcurve(p, nin, chan); 
		nsind0 = xfit_shpcurve(p, npin0, chan); 
		nsind1 = xfit_shpcurve(p, npin1, chan); 
#endif

//if (db) printf("~1 nsind %f, nsind0 %f, nsind1 %f\n",nsind,nsind0,nsind1);

		/* Offset and scale shaper in' value to match position span */
		rv = (nsind - nsind0)/(nsind1 - nsind0) * (npind1 - npind0) + npind0;

//if (db) printf("~1 scale offset ind %f\n",rv);

		/* de-normalize */
		rv = rv * (p->in_max[chan] - p->in_min[chan]) + p->in_min[chan];
//if (db) printf("~1 returning %d\n",rv);
	
	/* Just positioning curve */
	} else if ((p->tcomb & oc_ip) == oc_p) {
		rv = icxSTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], in,  
				                       p->in_min[chan], p->in_max[chan]);
	} else {
		rv = in;
	}
	return rv;
}

/* Lookup di values though the input positioning and shaper curves */
static void xfit_inpscurves(xfit *p, double *out, double *in) {
	int e;

	for (e = 0; e < p->di; e++)
		out[e] = xfit_inpscurve(p, in[e], e);
}

/* Inverse Lookup a value though the input positioning and shaper curves */
static double xfit_invinpscurve(xfit *p, double in, int chan) {
	double rv;
	/* Just shaper curve */
	if ((p->tcomb & oc_ip) == oc_i) {
		rv = icxInvSTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], in,  
			                       p->in_min[chan], p->in_max[chan]);

	/* shaper and positioning */
	} else if ((p->tcomb & oc_ip) == oc_ip) {
		double nind, nin;		/* normalized in', normalized in */
		int six;				/* Span index */
		double npind0, npind1;	/* normalized position in' span values */
		double npin0, npin1;	/* normalized position in span values */
		double nsind0, nsind1;	/* normalized shaper in' span values */
		double nsind;			/* normalized shaper in' value */

		/* Normalize */
		nind = (in - p->in_min[chan])/(p->in_max[chan] - p->in_min[chan]);

//if (db) printf("\n~1 invinpscurve: cha %d, input value %f, norm %f\n",chan,in,nind);

		/* Quantize to grid */
		six = (int)floor(nind * (p->gres[chan]-1.0));
		if (six > (p->gres[chan]-2))
			six = (p->gres[chan]-2);

		/* Compute span in' values */
		npind0 = six / (p->gres[chan]-1.0);
		npind1 = (six + 1.0) / (p->gres[chan]-1.0);

//if (db) printf("~1 six %d, npind0 %f, npind1 %f\n",six,npind0,npind1);

		/* Lookup span in values through position curve */
		npin0 = icxInvTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], npind0); 
		npin1 = icxInvTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], npind1); 
	
//if (db) printf("~1 npin0 %f, npin1 %f\n",npin0,npin1);

		/* Compute span shaper in' values */
#ifdef NEVER
		nsind0 = icxTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], npin0); 
		nsind1 = icxTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], npin1); 
#else
		nsind0 = xfit_shpcurve(p, npin0, chan); 
		nsind1 = xfit_shpcurve(p, npin1, chan); 
#endif

		/* Offset and scale position in' value to match shaper span */
		nsind = (nind - npind0)/(npind1 - npind0) * (nsind1 - nsind0) + nsind0;

//if (db) printf("~1 nsind %f, nsind0 %f, nsind1 %f\n",nsind,nsind0,nsind1);

		/* Invert through shaper curve */
#ifdef NEVER
		nin = icxInvTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], nsind); 
#else
		nin = xfit_invshpcurve(p, nsind, chan); 
#endif

		/* de-normalize */
		rv = nin * (p->in_max[chan] - p->in_min[chan]) + p->in_min[chan];

//if (db) printf("\n~1 nin = %f, returning %f\n",nin,rv);
	
	/* Just positioning curve */
	} else if ((p->tcomb & oc_ip) == oc_p) {
		rv = icxInvSTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], in,  
			                       p->in_min[chan], p->in_max[chan]);
	} else {
		rv = in;
	}
	return rv;
}

#ifdef NEVER
/* Check that inverse is working */
static double _xfit_inpscurve(xfit *p, double in, int chan) {
	double inv, rv, iinv;

	inv = in;
	rv = xfit_inpscurve(p, in, chan); 
	iinv = xfit_invinpscurve(p, rv, chan);

	if (fabs(in - iinv) > 1e-5)
		warning("xfit_inpscurve check, got %f, should be %f\n",iinv,in);
	
	return rv;
}
#define xfit_inpscurve _xfit_inpscurve
#endif /* NEVER */

/* Inverse Lookup di values though the input positioning and shaper curves */
static void xfit_invinpscurves(xfit *p, double *out, double *in) {
	int e;

	for (e = 0; e < p->di; e++)
		out[e] = xfit_invinpscurve(p, in[e], e);
}

/* - - - - - - - - - - - - - - - - - - - - */
/* Combined output transfer curve functions */

/* Lookup a value though an output curve */
static double xfit_outcurve(xfit *p, double in, int chan) {
	double rv;
	if (p->tcomb & oc_o) {
		if (p->flags & XFIT_OUT_LAB)
			rv = icxSTransFunc(p->v + p->out_offs[chan], p->oluord[chan], in,  
			                       p->out_min[chan], p->out_max[chan]);
		else
			rv = icxSTransFuncY2L(p->v + p->out_offs[chan], p->oluord[chan], in,  
			                       p->out_min[chan], p->out_max[chan]);
	} else
		rv = in;
	return rv;
}

/* Lookup fdi values though the output curves */
static void xfit_outcurves(xfit *p, double *out, double *in) {
	int f;
	
	if (p->flags & XFIT_OUT_LAB) {
		for (f = 0; f < p->fdi; f++) {
			double val = in[f];
			if (p->tcomb & oc_o)
				val = icxSTransFunc(p->v + p->out_offs[f], p->oluord[f], val,  
				                       p->out_min[f], p->out_max[f]);
			out[f] = val;
		}
	} else {
		for (f = 0; f < p->fdi; f++) {
			double val = in[f];
			if (p->tcomb & oc_o)
				val = icxSTransFuncY2L(p->v + p->out_offs[f], p->oluord[f], val,  
				                       p->out_min[f], p->out_max[f]);
			out[f] = val;
		}
	}
}

/* Inverse Lookup a value though an output curve */
static double xfit_invoutcurve(xfit *p, double in, int chan) {
	double rv;
	if (p->tcomb & oc_o) {
		if (p->flags & XFIT_OUT_LAB)
			rv = icxInvSTransFunc(p->v + p->out_offs[chan], p->oluord[chan], in,  
			                       p->out_min[chan], p->out_max[chan]);
		else
			rv = icxInvSTransFuncY2L(p->v + p->out_offs[chan], p->oluord[chan], in,  
			                       p->out_min[chan], p->out_max[chan]);
	} else
		rv = in;
	return rv;
}

/* Inverse Lookup fdi values though output curves */
static void xfit_invoutcurves(xfit *p, double *out, double *in) {
	int f;
	
	if (p->flags & XFIT_OUT_LAB) {
		for (f = 0; f < p->fdi; f++) {
			double val = in[f];
			if (p->tcomb & oc_o)
				val = icxInvSTransFunc(p->v + p->out_offs[f], p->oluord[f], val,  
				                       p->out_min[f], p->out_max[f]);
			out[f] = val;
		}
	} else {
		for (f = 0; f < p->fdi; f++) {
			double val = in[f];
			if (p->tcomb & oc_o)
				val = icxInvSTransFuncY2L(p->v + p->out_offs[f], p->oluord[f], val,  
				                       p->out_min[f], p->out_max[f]);
			out[f] = val;
		}
	}
}

/* - - - - - - - - - */

/* Convert an output value from absolute */
/* to relative using the current white point. */
static void xfit_abs_to_rel(xfit *p, double *out, double *in) {
	if (p->flags & XFIT_OUT_WP_REL) {
		if (p->flags & XFIT_OUT_LAB) {
			icmLab2XYZ(&icmD50, out, in);
			icmMulBy3x3(out, p->fromAbs, out);
			icmXYZ2Lab(&icmD50, out, out);
		} else {
			icmMulBy3x3(out, p->fromAbs, in);
		}
	} else {
		out[0] = in[0];
		out[1] = in[1];
		out[2] = in[2];
	}
}

/* Convert an XYZ output value from absolute */
/* to cLut relative using the current white point. */
static void xfit_XYZ_abs_to_rel(xfit *p, double *out, double *in) {
	if (p->flags & XFIT_OUT_WP_REL) {
		if (p->flags & XFIT_OUT_LAB) {
			icmMulBy3x3(out, p->fromAbs, in);
			icmXYZ2Lab(&icmD50, out, out);
		} else {
			icmMulBy3x3(out, p->fromAbs, in);
		}
	} else {
		if (p->flags & XFIT_OUT_LAB) {
			icmXYZ2Lab(&icmD50, out, in);
		} else {
			out[0] = in[0];
			out[1] = in[1];
			out[2] = in[2];
		}
	}
}

/* - - - - - - - - - */

/* return a weighting for the magnitude of the in and out */
/* shaping parameters squared. This is to reduce unconstrained "wiggles" */
static double shapmag(
xfit  *p			/* Base of optimisation structure */
) {
	double tt, w;
	double *b;			/* Base of parameters for this section */
	int di =  p->di;
	int fdi = p->fdi;
	int e, f, k;
	double iparam = 0.0;
	double oparam = 0.0;
	double dd;

	if (p->opt_msk & oc_i) {
		dd = SHAPE_WEIGHT/(double)(di);
		b = p->v + p->shp_off;
		for (e = 0; e < di; e++) {
			for (k = 0; k < p->iluord[e]; k++) {
				if (k <= 1) {
					w = SHAPE_HW01;
				} else if (k <= SHAPE_HBREAK) {
					double bl = (k - 1.0)/(SHAPE_HBREAK - 1.0);
					w = (1.0 - bl) * SHAPE_HW01 + bl * SHAPE_HWBR;
					w *= p->shp_smooth[e];
				} else {
					w = SHAPE_HWBR + (k-SHAPE_HBREAK) * SHAPE_HWINC;
					w *= p->shp_smooth[e];
				}
				tt = *b++;
				tt *= tt;		/* Squared */
				iparam += w * tt;
			}
		}
		iparam *= dd;
	}

	if (p->opt_msk & oc_o) {
		dd = SHAPE_WEIGHT/(double)(fdi);
		b = p->v + p->out_off;
		for (f = 0; f < fdi; f++) {
			for (k = 0; k < p->oluord[f]; k++) {
				if (k <= 1) {
					w = SHAPE_HW01;
				} else if (k <= SHAPE_HBREAK) {
					double bl = (k - 1.0)/(SHAPE_HBREAK - 1.0);
					w = (1.0 - bl) * SHAPE_HW01 + bl * SHAPE_HWBR;
					w *= p->out_smooth[f];
				} else {
					w = SHAPE_HWBR + (k-SHAPE_HBREAK) * SHAPE_HWINC;
					w *= p->out_smooth[f];
				}
				tt = *b++;
				tt *= tt;		/* Squared */
				oparam += w * tt;
			}
		}
		oparam *= dd;
	}
	return iparam + oparam;
}

/* return a weighting for the magnitude of the in and out */
/* shaping parameters. This is to reduce unconstrained "wiggles" */
/* Also sum the partial derivative for the parameters involved */
static double dshapmag(
xfit  *p,			/* Base of optimisation structure */
double *dav			/* Sum del's */
) {
	double tt, w;
	double *b, *c;			/* Base of parameters for this section */
	int di =  p->di;
	int fdi = p->fdi;
	int e, f, k;
	double iparam = 0.0;
	double oparam = 0.0;
	double dd;

	if (p->opt_msk & oc_i) {
		dd = SHAPE_WEIGHT/(double)(di);
		b = p->v + p->shp_off;
		c = dav + p->shp_off;
		for (e = 0; e < di; e++) {
			for (k = 0; k < p->iluord[e]; k++) {
				if (k <= 1) {
					w = SHAPE_HW01;
				} else if (k <= SHAPE_HBREAK) {	
					double bl = (k - 1.0)/(SHAPE_HBREAK - 1.0);
					w = (1.0 - bl) * SHAPE_HW01 + bl * SHAPE_HWBR;
					w *= p->shp_smooth[e];
				} else {
					w = SHAPE_HWBR + (k-SHAPE_HBREAK) * SHAPE_HWINC;
					w *= p->shp_smooth[e];
				}
				tt = *b++;
				*c++ += 2.0 * dd * w * tt;
				tt *= tt;			/* Squared */
				iparam += w * tt;
				
			}
		}
		iparam *= dd;
	}

	if (p->opt_msk & oc_o) {
		dd = SHAPE_WEIGHT/(double)(fdi);
		b = p->v + p->out_off;
		c = dav + p->out_off;
		for (f = 0; f < fdi; f++) {
			for (k = 0; k < p->oluord[f]; k++) {
				if (k <= 1) {
					w = SHAPE_HW01;
				} else if (k <= SHAPE_HBREAK) {
					double bl = (k - 1.0)/(SHAPE_HBREAK - 1.0);
					w = (1.0 - bl) * SHAPE_HW01 + bl * SHAPE_HWBR;
					w *= p->out_smooth[f];
				} else {
					w = SHAPE_HWBR + (k-SHAPE_HBREAK) * SHAPE_HWINC;
					w *= p->out_smooth[f];
				}
				tt = *b++;
				*c++ += 2.0 * dd * w * tt;
				tt *= tt;			/* Squared */
				oparam += w * tt;
			}
		}
		oparam *= dd;
	}
	return iparam + oparam;
}

/* Scale the shaper derivatives */
static void dshapscale(
xfit  *p,			/* Base of optimisation structure */
double *dav,		/* del's */
double scale		/* Scale factor */
) {
	double tt, w;
	double *b, *c;			/* Base of parameters for this section */
	int di =  p->di;
	int fdi = p->fdi;
	int e, f, k;

	if (p->opt_msk & oc_i) {
		c = dav + p->shp_off;
		for (e = 0; e < di; e++) {
			for (k = 0; k < p->iluord[e]; k++) {
				*c++ *= scale;
			}
		}
	}

	if (p->opt_msk & oc_o) {
		c = dav + p->out_off;
		for (f = 0; f < fdi; f++) {
			for (k = 0; k < p->oluord[f]; k++) {
				*c++ *= scale;
			}
		}
	}
}

/* Progress function */
static void xfitprog(void *pdata, int perc) {
	xfit *p = (xfit *)pdata;

	if (p->verb) {
		printf("%c% 3d%%",cr_char,perc); 
		if (perc == 100)
			printf("\n");
		fflush(stdout);
	}
}


int xfitfunc_trace = 1;

/* Shaper+Matrix optimisation function handed to powell() */
/* We simply minimize the total delta E squared, consistent with smoothness */
static double xfitfunc(void *edata, double *v) {
	xfit *p = (xfit *)edata;
	double tw = 0.0;				/* Total weight */
	double ev = 0.0, rv, smv;
	double tin[MXDI], out[MXDO];
	int di = p->di;
	int fdi = p->fdi;
	int i, e, f;

	/* Copy the parameters being optimised into xfit structure */

	/* Special case - a single shaper curve. The first sm_iluord params */
	/* are the common curve parameters, and the remainder are the matrix onwards */
	if (p->opt_ssch) {

		for (e = 0; e < di; e++) {	/* Duplicate and extend to per channel curve params */
			for (i = 0; i < p->sm_iluord; i++)
				p->v[p->shp_offs[e] + i] = v[i];
			for (; i < p->iluord[e]; i++)
				p->v[p->shp_offs[e] + i] = 0.0;
		}
		for (i = p->sm_iluord; i < p->opt_cnt; i++)
			p->v[p->mat_off + i - p->sm_iluord] = v[i];
	} else {
		for (i = 0; i < p->opt_cnt; i++) { 
//printf("~1 param %d = %f\n",i,v[i]);
			p->v[p->opt_off + i] = v[i];
		}
	}

	/* For all our data points */
	for (i = 0; i < p->nodp; i++) {
		double del;

		/* Apply input shaper channel curves */
		for (e = 0; e < di; e++)
			tin[e] = icxSTransFunc(p->v + p->shp_offs[e], p->iluord[e], p->rpoints[i].p[e],  
			                       p->in_min[e], p->in_max[e]);

		/* Apply matrix cube interpolation */
		icxCubeInterp(p->v + p->mat_off, fdi, di, out, tin);

		/* Apply output channel curves */
		for (f = 0; f < fdi; f++) {
			if (p->flags & XFIT_OUT_LAB) {
				out[f] = icxSTransFunc(p->v + p->out_offs[f], p->oluord[f], out[f],  
				                       p->out_min[f], p->out_max[f]);
			} else {
				out[f] = icxSTransFuncY2L(p->v + p->out_offs[f], p->oluord[f], out[f],  
				                       p->out_min[f], p->out_max[f]);
			}
		}

		/* Evaluate the error squared */
		if (p->flags & XFIT_FM_INPUT) {
			double pp[MXDI];
			for (e = 0; e < di; e++)
				pp[e] = p->rpoints[i].p[e];
			for (f = 0; f < fdi; f++) {
				double t1 = p->rpoints[i].v[f] - out[f];	/* Error in output */
	
				/* Create input point offset by equivalent delta to output point */
				/* error, in proportion to the partial derivatives for that output. */
				for (e = 0; e < di; e++)
					pp[e] += t1 * p->piv[i].ide[f][e];
			}
			del = p->to_de2(p->cntx2, pp, p->rpoints[i].p);
		} else {
			del = p->to_de2(p->cntx2, out, p->rpoints[i].v);
		}
		if (CURVEPOW > 1.0)
			del = pow(del, CURVEPOW);
		tw += p->rpoints[i].w;
		ev += p->rpoints[i].w * del;
	}

	/* Normalise error to be an average delta E squared */
	ev /= tw;

	/* Sum with shaper parameters squared, to */
	/* minimise unsconstrained "wiggles" */
	smv = shapmag(p);
	if (CURVEPOW > 1.0)
		smv = pow(smv, CURVEPOW);
	rv = ev + smv;

#ifdef DEBUG_PROGRESS
if (xfitfunc_trace)
fprintf(stdout,"~1(sm %f, ev %f)xfitfunc returning %f\n",smv,ev,rv);
#endif

	return rv;
}

/* Shaper+Matrix optimisation function with partial derivatives, */
/* handed to conjgrad() */
static double dxfitfunc(void *edata, double *dv, double *v) {
	xfit *p = (xfit *)edata;
	double tw = 0.0;				/* Total weight */
	double ev = 0.0, rv, smv;
	double tin[MXDI], out[MXDO];

	double dav[MXPARMS];				/* Overall del due to del param vals */
	double sdav[MXPARMS];				/* Overall del due to del smooth param vals */

	double dtin_iv[MXDI * MXLUORD];		/* Del in itrans out due to del itrans param vals */
	double dmato_mv[1 << MXDI];			/* Del in mat out due to del in matrix param vals */
	double dmato_tin[MXDO * MXDI];		/* Del in mat out due to del in matrix input values */
	double dout_ov[MXDO * MXLUORD];		/* Del in otrans out due to del in otrans param values */
	double dout_mato[MXDO];				/* Del in otrans out due to del in otrans input values */

	double dout_de[2][MXDIDO];			/* Del in DE due to two output values */

	int di = p->di;
	int fdi = p->fdi;
	int i, jj, k, e, ee, f, ff;

	/* Copy the parameters being optimised into xfit structure */

	/* Special case - a single shaper curve. The first sm_iluord params */
	/* are the common curve parameters, and the remainder are the matrix onwards */
	if (p->opt_ssch) {
		for (e = 0; e < di; e++) {	/* Duplicate and extend to per channel curve params */
			for (i = 0; i < p->sm_iluord; i++)
				p->v[p->shp_offs[e] + i] = v[i];
			for (; i < p->iluord[e]; i++)
				p->v[p->shp_offs[e] + i] = 0.0;
		}
		for (i = p->sm_iluord; i < p->opt_cnt; i++) 
			p->v[p->mat_off + i - p->sm_iluord] = v[i];

	} else {
		for (i = 0; i < p->opt_cnt; i++) { 
			p->v[p->opt_off + i] = v[i];
		}
	}

	/* Zero the accumulated partial derivatives */
	/* We compute deriv for all parameters (not just current optimised) */
	for (i = 0; i < p->tot_cnt; i++)
		dav[i] = 0.0;

	/* For all our data points */
	for (i = 0; i < p->nodp; i++) {
		double del;

		/* Apply input channel curves */
		for (e = 0; e < di; e++)
			tin[e] = icxdpSTransFunc(p->v + p->shp_offs[e], &dtin_iv[p->shp_offs[e] - p->shp_off],
		                         p->iluord[e], p->rpoints[i].p[e], p->in_min[e], p->in_max[e]);

		/* Apply matrix cube interpolation */
		icxdpdiCubeInterp(p->v + p->mat_off, dmato_mv, dmato_tin, fdi, di, out, tin);

		/* Apply output channel curves */
		for (f = 0; f < fdi; f++) {
			if (p->flags & XFIT_OUT_LAB)
				out[f] = icxdpdiSTransFunc(p->v + p->out_offs[f],
				                           &dout_ov[p->out_offs[f] - p->out_off], &dout_mato[f],
				                           p->oluord[f], out[f], p->out_min[f], p->out_max[f]);
			else
				out[f] = icxdpdiSTransFuncY2L(p->v + p->out_offs[f],
				                           &dout_ov[p->out_offs[f] - p->out_off], &dout_mato[f],
				                           p->oluord[f], out[f], p->out_min[f], p->out_max[f]);
		}

		/* Convert to Delta E and compute pde's into dout_de squared */
		if (p->flags & XFIT_FM_INPUT) {
			double tdout_de[2][MXDIDO];
			double pp[MXDI];
			for (e = 0; e < di; e++)
				pp[e] = p->rpoints[i].p[e];
			for (f = 0; f < fdi; f++) {
				double t1 = p->rpoints[i].v[f] - out[f];	/* Error in output */
	
				/* Create input point offset by equivalent delta to output point */
				/* error, in proportion to the partial derivatives for that output. */
				for (e = 0; e < di; e++)
					pp[e] += t1 * p->piv[i].ide[f][e];
			}
			del = p->to_dde2(p->cntx2, tdout_de, pp, p->rpoints[i].p);
			if (CURVEPOW > 1.0) {
				double dadj;
				dadj = CURVEPOW * pow(del, CURVEPOW - 1.0); /* Adjust derivative accordingly */
				del = pow(del, CURVEPOW);
				for (e = 0; e < di; e++)
					tdout_de[0][e] *= dadj;
			}

			/* Compute partial derivative */
			for (e = 0; e < di; e++) {
				dout_de[0][e] = 0.0;
				for (f = 0; f < fdi; f++) {
					dout_de[0][e] += tdout_de[0][e] * p->piv[i].ide[f][e];
				}
			}
		} else {
			del = p->to_dde2(p->cntx2, dout_de, out, p->rpoints[i].v);
			if (CURVEPOW > 1.0) {
				double dadj;
				dadj = CURVEPOW * pow(del, CURVEPOW - 1.0); /* Adjust derivative accordingly */
				del = pow(del, CURVEPOW);
				for (f = 0; f < fdi; f++)
					dout_de[0][f] *= dadj;
			}
		}

		/* Accumulate total weighted delta E squared */
		tw += p->rpoints[i].w;
		ev += p->rpoints[i].w * del;

		/* Compute and accumulate partial difference values for each parameter value */
		if (p->opt_msk & oc_i) {
			/* Input transfer parameters */
			for (ee = 0; ee < di; ee++) {				/* Parameter input chanel */
				for (k = 0; k < p->iluord[ee]; k++) {	/* Param within channel */
					double vv = 0.0;
					jj = p->shp_offs[ee] - p->shp_off + k;	/* Overall input trans param */

//					for (ff = 0; ff < 3; ff++) {		/* Lab channels */
					for (ff = 0; ff < fdi; ff++) {		/* Output channels */
						vv += dout_de[0][ff] * dout_mato[ff]
						    * dmato_tin[ff * di + ee] * dtin_iv[jj];
					}
					dav[p->shp_off + jj] += p->rpoints[i].w * vv;
				}
			}
		}

		if (p->opt_msk & oc_m) {
			/* Matrix parameters */
			for (ff = 0; ff < fdi; ff++) {				/* Parameter output chanel */
				for (ee = 0; ee < (1 << di); ee++) {	/* Matrix input combination chanel */
					double vv = 0.0;
					jj = ff * (1 << di) + ee;			/* Matrix Parameter index */

					vv += dout_de[0][ff] * dout_mato[ff] * dmato_mv[ee];
					dav[p->mat_off + jj] += p->rpoints[i].w * vv;
				}
			}
		}

		if (p->opt_msk & oc_o) {
			/* Output transfer parameters */
			for (ff = 0; ff < fdi; ff++) {				/* Parameter output chanel */
				for (k = 0; k < p->oluord[ff]; k++) {	/* Param within channel */
					double vv = 0.0;
					jj = p->out_offs[ff] - p->out_off + k;	/* Overall output trans param */

					vv += dout_de[0][ff] * dout_ov[jj];
					dav[p->out_off + jj] += p->rpoints[i].w * vv;
				}
			}
		}
	}

	/* Normalise error to be an average delta E squared */
	ev /= tw;
	for (i = 0; i < p->tot_cnt; i++) {
		dav[i] /= tw;
		sdav[i] = 0.0;
	}

	/* Sum with shaper parameters squared, to */
	/* minimise unsconstrained "wiggles" */
	/* Compute partial derivative wrt those parameters too */
	smv = dshapmag(p, sdav);
	if (CURVEPOW > 1.0) {
		double dadj;
		dadj = CURVEPOW * pow(smv, CURVEPOW - 1.0); /* Adjust derivative accordingly */
		smv = pow(smv, CURVEPOW);
		dshapscale(p, sdav, dadj);		/* Scale the partial derivatives */
	}
	rv = ev + smv;

	/* Sum the del for parameters being optimised and copy to return array */
	if (p->opt_ssch) {
		for (i = 0; i < p->sm_iluord; i++)
			dv[i] = 0.0;
		for (e = 0; e < di; e++) {	/* Combine per channel curve de's */
			for (i = 0; i < p->sm_iluord; i++)
				dv[i] += dav[p->shp_offs[e] + i] + sdav[p->shp_offs[e] + i];
		}
		for (i = p->sm_iluord; i < p->opt_cnt; i++)	/* matrix and rest de's */ 
			dv[i] = dav[p->mat_off + i - p->sm_iluord] + sdav[p->mat_off + i - p->sm_iluord];

	} else {
		for (i = 0; i < p->opt_cnt; i++)
			dv[i] = dav[p->opt_off + i] + sdav[p->opt_off + i];
	}

#ifdef DEBUG_PROGRESS
fprintf(stdout,"~1(sm %f, ev %f)dxfitfunc returning %f\n",smv,ev,rv);
#endif

	return rv;
}

#ifdef NEVER
/* Check partial derivative function within xfitfunc() [Intensive check] */

static double _xfitfunc(void *edata, double *v) {
	xfit *p = (xfit *)edata;
	int i;
	double dv[MXPARMS];
	double rv, drv;
	double trv;
	int verb;
	
	rv = xfitfunc(edata, v);
	verb = p->verb;
	p->verb = 0;
	drv = dxfitfunc(edata, dv, v);
	p->verb = verb;

	if (fabs(rv - drv) > 1e-6)
		printf("######## RV MISMATCH is %f should be %f ########\n",rv,drv);

	/* Check each parameter delta */
	xfitfunc_trace = 0;
	for (i = 0; i < p->opt_cnt; i++) {
		double del;

		v[i] += 1e-7;
		trv = xfitfunc(edata, v);
		v[i] -= 1e-7;
		
		/* Check that del is correct */
		del = (trv - rv)/1e-7;
		if (fabs(dv[i] - del) > 0.04) {
//printf("~1 del = %f from (trv %f - rv %f)/0.1\n",del,trv,rv);
			printf("######## EXCESSIVE at v[%d] is %f should be %f ########\n",i,dv[i],del);
		}
	}
	xfitfunc_trace = 1;
	return rv;
}

#define xfitfunc _xfitfunc
#endif	/* NEVER */

#ifdef NEVER
/* Check partial derivative function within dxfitfunc() [Less intensive check] */

static double _dxfitfunc(void *edata, double *dv, double *v) {
	xfit *p = (xfit *)edata;
	int i;
	double rv, drv;
	double trv;
	int verb;
	int exec = 0;
	
	rv = xfitfunc(edata, v);
	verb = p->verb;
	p->verb = 0;
	drv = dxfitfunc(edata, dv, v);
	p->verb = verb;

	if (fabs(rv - drv) > 1e-6)
		printf("######## RV MISMATCH is %f should be %f ########\n",rv,drv);

	/* Check each parameter delta */
	xfitfunc_trace = 0;
	for (i = 0; i < p->opt_cnt; i++) {
		double del;

		v[i] += 1e-7;
		trv = xfitfunc(edata, v);
		v[i] -= 1e-7;
		
		/* Check that del is correct */
		del = (trv - rv)/1e-7;
		if (fabs(dv[i] - del) > 0.04) {
//printf("~1 del = %f from (trv %f - rv %f)/0.1\n",del,trv,rv);
			printf("######## EXCESSIVE at v[%d] is %f should be %f ########\n",i,dv[i],del);
			exec = 1;
		}
	}
#ifdef NEVER
	if (exec) {
		printf("~1 parameters are:\n");
		for (i = 0; i < p->opt_cnt; i++)
			printf("p->wv[%d] = %f;\n",i,v[i]);
		exit(1);
	}
#endif
	xfitfunc_trace = 1;
	return rv;
}

#define dxfitfunc _dxfitfunc
#endif	/* NEVER */

/* - - - - - - - - - */

/* Output curve symetry optimisation function handed to powell() */
/* Just the order 0 value will be adjusted */
static double symoptfunc(void *edata, double *v) {
	xfit *p = (xfit *)edata;
	double out[1], in[1] = { 0.0 };
	int ch = p->opt_ch;		/* Output channel being adjusted for symetry */
	double rv;

	/* Copy the parameter being tested back into xfit */
	p->v[p->out_offs[ch]] = v[0];
	if (p->flags & XFIT_OUT_LAB)
		*out = icxSTransFunc(p->v + p->out_offs[ch], p->oluord[ch], *in,
							   p->out_min[ch], p->out_max[ch]);
	else
		*out = icxSTransFuncY2L(p->v + p->out_offs[ch], p->oluord[ch], *in,
							   p->out_min[ch], p->out_max[ch]);

	rv = out[0] * out[0];

#ifdef DEBUG_PROGRESS
printf("~1symoptfunc returning %f\n",rv);
#endif
	return rv;
}

/* - - - - - - - - - */

/* Set up for an optimisation run: */
/* Figure out parameters being optimised, */
/* copy them to start values, */
/* init and scale the search radius */
static void setup_xfit(
xfit *p,
double *wv,		/* Return parameters to hand to optimiser */
double *sa,		/* Return search radius to hand to optimiser */
double transrad,/* Nominal transfer curve radius, 0.0 - 3.0 */
double matrad	/* Nominal matrix radius, 0.0 - 1.0 */
) {
	int i;
	p->opt_off = -1;
	p->opt_cnt = 0;
	
	if (p->opt_msk & oc_i) {

		if (p->opt_ssch) {	/* Special case - should only be used first, */
							/* Fitting a sigle common input shaper curve. */
			if (p->opt_off < 0)
				p->opt_off = p->mat_off - p->sm_iluord;	/* Shouldn't be used... */
			p->opt_cnt += p->sm_iluord;
		
			for (i = 0; i < p->sm_iluord; i++) {
				*wv++ = 0.0;
				*sa++ = transrad;
			}

		} else {	/* Initial or continuing fitting of all the curves */
			if (p->opt_off < 0)
				p->opt_off = p->shp_off;
			p->opt_cnt += p->shp_cnt;
		
			for (i = 0; i < p->shp_cnt; i++) {
				*wv++ = p->v[p->shp_off + i];
				*sa++ = transrad;
			}
		}
	}
	if (p->opt_msk & oc_m) {
		if (p->opt_off < 0)
			p->opt_off = p->mat_off;
		p->opt_cnt += p->mat_cnt;

		for (i = 0; i < p->mat_cnt; i++) {
			*wv++ = p->v[p->mat_off + i];
			*sa++ = matrad;
		}
	}
	if (p->opt_msk & oc_o) {
		if (p->opt_off < 0)
			p->opt_off = p->out_off;
		p->opt_cnt += p->out_cnt;

		for (i = 0; i < p->out_cnt; i++) {
			*wv++ = p->v[p->out_off + i];
			*sa++ = transrad;
		}
	}
	if (p->opt_cnt > MXPARMS)
		error("setup_xfit: asert, %d exceeded MXPARMS %d",p->opt_cnt,MXPARMS);

#ifdef DEBUG
	printf("setup_xfit() got opt_msk 0x%x, opt_off %d, opt_cnt = %d\n",
	                                  p->opt_msk,p->opt_off,p->opt_cnt);
#endif /* DEBUG */
}

#ifdef DEBUG
/* Diagnostic */
static void dump_xfit(
xfit *p
) {
	int i, e, f;
	double *b;			/* Base of parameters for this section */
	int di, fdi;
	di   = p->di;
	fdi  = p->fdi;

	/* Input positioning curve */
	b = p->v + p->pos_off;
	for (e = 0; e < di; b += p->iluord[e], e++) {
		printf("pos %d = ",e);
		for (i = 0; i < p->iluord[e]; i++)
			printf("%f ",b[i]);
		printf("\n");
	}

	/* Input shaper curve */
	b = p->v + p->shp_off;
	for (e = 0; e < di; b += p->iluord[e], e++) {
		printf("shp %d = ",e);
		for (i = 0; i < p->iluord[e]; i++)
			printf("%f ",b[i]);
		printf("\n");
	}

	/* Matrix */
	b = p->v + p->mat_off;
	for (e = 0; e < (1 << di); e++) {
		printf("mx %d = ",e);
		for (f = 0; f < fdi; f++)
			printf("%f ",*b++);
		printf("\n");
	}

	/* Output curve */
	b = p->v + p->out_off;
	for (f = 0; f < fdi; b += p->oluord[f], f++) {
		printf("out %d = ",f);
		for (i = 0; i < p->oluord[f]; i++)
			printf("%f ",b[i]);
		printf("\n");
	}
}
#endif /* DEBUG */

/* - - - - - - - - - */

/* Setup the pseudo inverse information for each test point, */
/* using the current model. */
static void setup_piv(xfit *p) {
	int di =  p->di;
	int fdi = p->fdi;
	int i, e, f;

	/* Estimate in -> out partial derivatives */
	for (i = 0; i < p->nodp; i++) {
		double pd[MXDO][MXDI];
		double pp[MXDI];
		double vv[MXDIDO];

		/* Estimate in -> out partial derivatives */
		for (e = 0; e < di; e++)
			pp[e] = p->ipoints[i].p[e];

		/* Apply input shaper channel curves */
		for (e = 0; e < di; e++)
			vv[e] = icxSTransFunc(p->v + p->shp_offs[e], p->iluord[e], pp[e],  
			                       p->in_min[e], p->in_max[e]);

		/* Apply matrix cube interpolation */
		icxCubeInterp(p->v + p->mat_off, fdi, di, vv, vv);

		/* Apply output channel curves */
		for (f = 0; f < fdi; f++) {
			if (p->flags & XFIT_OUT_LAB)
				vv[f] = icxSTransFunc(p->v + p->out_offs[f], p->oluord[f], vv[f],  
				                       p->out_min[f], p->out_max[f]);
			else
				vv[f] = icxSTransFuncY2L(p->v + p->out_offs[f], p->oluord[f], vv[f],  
				                       p->out_min[f], p->out_max[f]);
		}
		
		for (e = 0; e < di; e++) {
			double tt[MXDIDO];
	
			pp[e] += 1e-4;

			/* Apply input shaper channel curves */
			for (e = 0; e < di; e++)
				tt[e] = icxSTransFunc(p->v + p->shp_offs[e], p->iluord[e], pp[e],  
				                       p->in_min[e], p->in_max[e]);
	
			/* Apply matrix cube interpolation */
			icxCubeInterp(p->v + p->mat_off, fdi, di, tt, tt);
	
			/* Apply output channel curves */
			for (f = 0; f < fdi; f++) {
				if (p->flags & XFIT_OUT_LAB)
					tt[f] = icxSTransFunc(p->v + p->out_offs[f], p->oluord[f], tt[f],  
					                       p->out_min[f], p->out_max[f]);
				else
					tt[f] = icxSTransFuncY2L(p->v + p->out_offs[f], p->oluord[f], tt[f],  
					                       p->out_min[f], p->out_max[f]);
			}
	
			for (f = 0; f < p->fdi; f++)
				pd[f][e] = (tt[f] - vv[f])/1e-4;
		
			pp[e] -= 1e-4;
		}

		/* Compute a psudo inverse matrix to map rout delta E to */
		/* in delta E in proportion to the pd magnitude. */
		for (f = 0; f < fdi; f++) {
			double ss = 0.0;

			for (e = 0; e < di; e++)		/* Sum of pd's ^4 */
				ss += pd[f][e] * pd[f][e] * pd[f][e] * pd[f][e];
			ss = sqrt(ss);
			if (ss > 1e-8) {
				for (e = 0; e < di; e++)
					p->piv[i].ide[f][e] = pd[f][e]/ss;
			} else {		/* Hmm. */
				for (e = 0; e < di; e++)
					p->piv[i].ide[f][e] = 0.0;
			}
		}
	}
}

/* - - - - - - - - - */

/* Function to pass to rspl to re-set output values, */
/* to account for skeleton model offset. */
static void
skm_rspl_out(
	void *pp,			/* relativectx structure */
	double *out,		/* output value */
	double *in			/* input value */
) {
	xfit *p    = (xfit *)pp;
	int f, fdi = p->fdi;
	double inval[MXDI];
	double skval[MXDO];

	/* Look up the skeleton value for this grid point */
	xfit_invinpscurves(p, inval, in);		/* Back to input values */
	p->skm->lookup(p->skm, skval, inval);	/* Skm */
	xfit_abs_to_rel(p, skval, skval);
	xfit_invoutcurves(p, skval, skval);

	for (f = 0; f < fdi; f++)
		out[f] += skval[f]; 			/* Add it back */
}

/* Weak function rspl callback (not used) */
void skm_weak(void *cbntx, double *out, double *in) {
	xfit *p    = (xfit *)cbntx;

#ifndef NEVER
	int f, fdi = p->fdi;

	for (f = 0; f < fdi; f++)
		out[f] = 0.0;			/* Deviation from skeleton should tend to zero */

#else		/* Skeleton as weak atractor */
	int f, fdi = p->fdi;
	double inval[MXDI];

	/* Look up the skeleton value for this grid point */
	xfit_invinpscurves(p, inval, in);		/* Back to input values */
	p->skm->lookup(p->skm, out, inval);	/* Skm */
	xfit_abs_to_rel(p, out, out);
	xfit_invoutcurves(p, out, out);

#endif
}

/* - - - - - - - - - */

/* Function to pass to rspl to re-set output values, */
/* to make them relative to the white and black points */
static void
conv_rspl_out(
	void *pp,			/* relativectx structure */
	double *out,		/* output value */
	double *in			/* input value */
) {
	xfit *p    = (xfit *)pp;
	double tt[3];

	/* Convert the clut values to output values */
	xfit_outcurves(p, tt, out);

	if (p->flags & XFIT_OUT_LAB) {
		icmLab2XYZ(&icmD50, tt, tt);
		icmMulBy3x3(out, p->cmat, tt);
		icmXYZ2Lab(&icmD50, out, out);

	} else {	/* We are all in XYZ */
		icmMulBy3x3(out, p->cmat, tt);
	}

	/* And then convert them back to clut values */
	xfit_invoutcurves(p, out, out);
}

/* Function to pass to rspl to re-set output values, */
/* to clip any with Y over 1.0 to D50 */
static void
clip_rspl_out(
	void *pp,			/* relativectx structure */
	double *out,		/* output value */
	double *in			/* input value */
) {
	xfit *p    = (xfit *)pp;
	double tt[3];

	/* Convert the clut values to output values */
	xfit_outcurves(p, tt, out);

	if (p->flags & XFIT_OUT_LAB) {
		if (tt[0] > 100.0)
			icmCpy3(out, p->cmat[0]);
	} else {
		if (tt[1] > 1.0)
			icmCpy3(out, p->cmat[0]);
	}
}

//#ifdef SPECIAL_TEST
/* - - - - - - - - - */

/* Execute the linear XYZ device model */
static void domodel(double *out, double *in) {
	double tmp[3];
	int i, j;
	double col[3][3];	/* sRGB additive colorant values in XYZ :- [out][in]  */

	col[0][0] = 0.412424;	/* X from R */
	col[0][1] = 0.357579;	/* X from G */
	col[0][2] = 0.180464;	/* X from B */
	col[1][0] = 0.212656;	/* Y from R */
	col[1][1] = 0.715158;	/* Y from G */
	col[1][2] = 0.0721856;	/* Y from B */
	col[2][0] = 0.0193324;	/* Z from R */
	col[2][1] = 0.119193;	/* Z from G */
	col[2][2] = 0.950444;	/* Z from B */

#ifdef SPECIAL_TEST_GAMMA	
	tmp[0] = pow(in[0], 1.9);
	tmp[1] = pow(in[1], 2.0);
	tmp[2] = pow(in[2], 2.1);
#else
	tmp[0] = in[0];
	tmp[1] = in[1];
	tmp[2] = in[2];
#endif

	for (j = 0; j < 3; j++) {
		out[j] = 0.0;
		for (i = 0; i < 3; i++)
			out[j] += col[j][i] * tmp[i];
	}
}

#ifdef SPECIAL_FORCE
/* Function to pass to rspl to set nodes against */
/* synthetic model. */
static void
set_rspl_out1(
	void *pp,			/* relativectx structure */
	double *out,		/* output value */
	double *in			/* input value */
) {
	xfit *p  = (xfit *)pp;
	double tt[3], tout[3];

	/* Convert the input' values to input values */
	xfit_invinpscurves(p, tt, in);

	/* Synthetic linear rgb->XYZ model */
	domodel(tout, tt);

	/* Apply abs->rel white point adjustment */
	icmMulBy3x3(tout, p->cmat, tout);

#ifdef SPECIAL_TEST_LAB
	icmXYZ2Lab(&icmD50, tout, tout);
#endif
	/* And then convert them back to clut values */
	xfit_invoutcurves(p, tout, tout);

#ifdef DEBUG
printf("~1 changing %f %f %f -> %f %f %f\n", out[0], out[1], out[2], tout[0], tout[1], tout[2]);
#endif

	out[0] = tout[0];
	out[1] = tout[1];
	out[2] = tout[2];
}

#endif /* SPECIAL_FORCE */

/* - - - - - - - - - */
/* Do the fitting. */
/* return nz on error */
/* 1 = malloc or other error */
static int xfit_fit(
	struct _xfit *p, 
	int flags,				/* Flag values */
	int di,					/* Input dimensions */
	int fdi,				/* Output dimensions */
	int rsplflags,			/* clut rspl creation flags */
	double *wp,				/* if flags & XFIT_OUT_WP_REL or XFIT_OUT_WP_REL_US, */
							/* Initial white point, returns final wp */
	double *dw,				/* Device white value to adjust to be D50 */
	double wpscale,			/* If >= 0.0 scale final wp */  
	double *dgw,			/* Device space gamut boundary white for XFIT_OUT_WP_REL_US */
							/* (ie. RGB 1,1,1 CMYK 0,0,0,0, etc) */
	cow *ipoints,			/* Array of data points to fit - referece taken */
	int nodp,				/* Number of data points */
	icxMatrixModel *skm,	/* Optional skeleton model (used for input profiles) */
	double in_min[MXDI],	/* Input value scaling/domain minimum */
	double in_max[MXDI],	/* Input value scaling/domain maximum */
	int gres[MXDI],			/* clut resolutions being optimised for/returned */
	double out_min[MXDO],	/* Output value scaling/range minimum */
	double out_max[MXDO],	/* Output value scaling/range maximum */
//	co *bpo,				/* If != NULL, black point override in same spaces as ipoints */
	double smooth,			/* clut rspl smoothing factor */
	double oavgdev[MXDO],	/* Average output value deviation */
	double demph,			/* dark emphasis factor for cLUT grid res. */
	int iord[],				/* Order of input pos/shaper curve for each dimension */
	int sord[],				/* Order of input sub-grid shaper curve (not used) */
	int oord[],				/* Order of output shaper curve for each dimension */
	double shp_smooth[MXDI],/* Smoothing factors for each curve, nom = 1.0 */
	double out_smooth[MXDO],
	optcomb tcomb,			/* Flag - target elements to fit. */
	void *cntx2,			/* Context of callbacks */
							/* Callback to convert two fit values delta E squared */
	double (*to_de2)(void *cntx, double *in1, double *in2),
							/* Same as above, with partial derivatives */
	double (*to_dde2)(void *cntx, double dout[2][MXDIDO], double *in1, double *in2)
) {
	int i, e, f;
	double *b;				/* Base of parameters for this section */
	int poff;

	double powtol = POWTOL1;		/* powell/conjgrad initial tollerance */
	int maxits = MAXITS1;			/* powell/conjgrad initial maximum itterations */

	if (tcomb & oc_io)		/* If we're doing anything, we need the matrix */
		tcomb |= oc_m;

	p->flags  = flags;
	if (flags & XFIT_VERB)
		p->verb = 1;
	else
		p->verb = 0;
	p->di      = di;
	p->fdi     = fdi;
	p->wp      = wp;		/* Take reference, so modified wp can be returned */
	p->dw      = dw;
	p->nodp    = nodp;
	p->skm     = skm;		/* This isn't current used by profin, because it doesn't help.. */
	p->ipoints = ipoints;
	p->tcomb   = tcomb;
	p->cntx2   = cntx2;
	for (e = 0; e < di; e++)
		p->gres[e] = gres[e];
	p->to_de2  = to_de2;
	p->to_dde2 = to_dde2;

#ifdef DEBUG
	printf("xfit_fit called with flags = 0x%x, di = %d, fdi = %d, nodp = %d, tcomb = 0x%x\n",flags,di,fdi,nodp,tcomb);
#endif

//printf("~1 out min = %f %f %f max = %f %f %f\n", out_min[0], out_min[1], out_min[2], out_max[0], out_max[1], out_max[2]);

	/* Sanity protect shaper orders */
	/* and save scaling and smoothness factors. */
	p->sm_iluord = MXLUORD+1;
	for (e = 0; e < di; e++) {
		if (iord[e] > MXLUORD)
			p->iluord[e] = MXLUORD;
		else
			p->iluord[e] = iord[e];
		if (p->iluord[e] < p->sm_iluord)
			p->sm_iluord = p->iluord[e];
		p->in_min[e] = in_min[e];
		p->in_max[e] = in_max[e];
		p->shp_smooth[e] = shp_smooth[e];
	}
	for (f = 0; f < fdi; f++) {
		if (oord[f] > MXLUORD)
			p->oluord[f] = MXLUORD;
		else
			p->oluord[f] = oord[f];
		p->out_min[f] = out_min[f];
		p->out_max[f] = out_max[f];
		p->out_smooth[f] = out_smooth[f];
	}


	/* Compute parameter offset and count information */
	p->shp_off = 0;
	for (poff = p->shp_off, p->shp_cnt = 0, e = 0; e < di; e++) {
		p->shp_offs[e] = poff;
		p->shp_cnt += p->iluord[e];
		poff       += p->iluord[e];
	}

	p->mat_off = p->shp_off + p->shp_cnt;
	for (poff = p->mat_off, p->mat_cnt = 0, f = 0; f < fdi; f++) {
		p->mat_offs[f] = poff;
		p->mat_cnt += (1 << di);
		poff       += (1 << di);
	}

	p->out_off = p->mat_off + p->mat_cnt;
	for (poff = p->out_off, p->out_cnt = 0, f = 0; f < fdi; f++) {
		p->out_offs[f] = poff;
		p->out_cnt += p->oluord[f];
		poff       += p->oluord[f];
	}

	p->pos_off = p->out_off + p->out_cnt;
	for (poff = p->pos_off, p->pos_cnt = 0, e = 0; e < di; e++) {
		p->pos_offs[e] = poff;
		p->pos_cnt += p->iluord[e];
		poff       += p->iluord[e];
	}

	p->tot_cnt = p->shp_cnt + p->mat_cnt + p->out_cnt + p->pos_cnt;
	if (p->tot_cnt > MXPARMS)
		error("xfit_fit: assert tot_cnt exceeds MXPARMS");

	/* Allocate space for parameter values */
	if (p->v != NULL) {
		free(p->v);
		p->v = NULL;
	}
	if (p->wv != NULL) {
		free(p->wv);
		p->wv = NULL;
	}
	if (p->sa != NULL) {
		free(p->sa);
		p->sa = NULL;
	}
	if ((p->v = (double *)calloc(p->tot_cnt, sizeof(double))) == NULL)
		return 1;
	if ((p->wv = (double *)calloc(p->tot_cnt, sizeof(double))) == NULL)
		return 1;
	if ((p->sa = (double *)calloc(p->tot_cnt, sizeof(double))) == NULL)
		return 1;

	/* Setup initial white point abs->rel conversions */
	if ((p->flags & XFIT_OUT_WP_REL) != 0) { 
		icmXYZNumber _wp;

		icmAry2XYZ(_wp, p->wp);

		/* Absolute->Aprox. Relative Adaptation matrix and */
		/* Aproximate relative to absolute conversion matrix */
		if (p->picc != NULL) {
			p->picc->chromAdaptMatrix(p->picc, ICM_CAM_NONE, p->toAbs, p->fromAbs, icmD50, _wp);
		} else {
			icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, _wp, p->fromAbs);
			icmChromAdaptMatrix(ICM_CAM_BRADFORD, _wp, icmD50, p->toAbs);
		}
	
		if (p->verb) {
			double lab[3];
			icmXYZ2Lab(&icmD50, lab, p->wp);
			printf("Initial White Point XYZ %f %f %f, Lab %f %f %f\n",
			    p->wp[0], p->wp[1], p->wp[2], lab[0], lab[1], lab[2]);
		}
	} else {
		icmSetUnity3x3(p->fromAbs);
		icmSetUnity3x3(p->toAbs);
		icmCpy3(p->wp, icmD50_ary3);		/* Default to forced absolute */
	}

	/* Setup input position/shape curves to be linear initially */
	b = p->v + p->shp_off;
	for (e = 0; e < di; b += p->iluord[e], e++) {
		for (i = 0; i < p->iluord[e]; i++) {
			b[i] = 0.0;
		}
	}

	/* Setup matrix to be pure colorant' values initially */
	b = p->v + p->mat_off;
	for (e = 0; e < (1 << di); e++) {	/* For each colorant combination */
		int j, k, bk = 0;
		double bdif = 1e6;
		double ov[MXDO];
	
		/* Search the patch list to find the one closest to this input combination */
		for (k = 0; k < p->nodp; k++) {
			double dif = 0.0;
			for (j = 0; j < di; j++) {
				double tt;
				if (e & (1 << j))
					tt = p->in_max[j] - p->ipoints[k].p[j];
				else
					tt = p->in_min[j] - p->ipoints[k].p[j];
				dif += tt * tt;
			}
			if (dif < bdif) {		/* best so far */
				bdif = dif;
				bk = k;
				if (dif < 0.001)
					break;			/* Don't bother looking further */
			}
		}
		xfit_abs_to_rel(p, ov, p->ipoints[bk].v);

		for (f = 0; f < fdi; f++)
			b[f * (1 << di) + e] = ov[f];
	}

	/* Setup output curves to be linear initially */
	b = p->v + p->out_off;
	for (f = 0; f < fdi; b += p->oluord[f], f++) {
		for (i = 0; i < p->oluord[f]; i++) {
			b[i] = 0.0;
		}
	}

	/* Setup positioning curves to be linear initially */
	b = p->v + p->pos_off;
	for (e = 0; e < di; b += p->iluord[e], e++) {
		for (i = 0; i < p->iluord[e]; i++) {
			b[i] = 0.0;
		}
	}

	/* Create copy of input points with output converted to white relative */
	if (p->rpoints == NULL) {
		if ((p->rpoints = (cow *)malloc(p->nodp * sizeof(cow))) == NULL)
			return 1;
	}
	for (i = 0; i < p->nodp; i++) {
		p->rpoints[i].w = p->ipoints[i].w;
		for (e = 0; e < di; e++)
			p->rpoints[i].p[e] = p->ipoints[i].p[e];
		for (f = 0; f < fdi; f++)
			p->rpoints[i].v[f] = p->ipoints[i].v[f];

		/* out -> rout */
		xfit_abs_to_rel(p, p->rpoints[i].v, p->rpoints[i].v);
	}
  
	/* Allocate array of pseudo-inverse matricies */
	if ((p->flags & XFIT_FM_INPUT) != 0 && p->piv == NULL) {
		if ((p->piv = (xfit_piv *)malloc(p->nodp * sizeof(xfit_piv))) == NULL)
			return 1;
	}

	/* Allocate array of span DE's for current opt channel */
	{
		int lres = 0;
		for (e = 0; e < di; e++) {
			if (p->gres[e] > lres)
				lres = p->gres[e];
		}
		if ((p->uerrv = (double *)malloc(lres * sizeof(double))) == NULL)
			return 1;
	}

	/* Do the fitting one part at a time, then together */
	/* Shaper curves are created if position or shaper curves are requested */

	/* Fit just the matrix */
	if ((p->tcomb & oc_ipo) != 0
	 && (p->tcomb & oc_m) == oc_m) {	/* Only bother with matrix if in and/or out */
		double rerr;

		if (p->verb)
			printf("About to optimise temporary matrix\n");

		/* Setup pseudo-inverse if we need it */
		if (p->flags & XFIT_FM_INPUT)
			setup_piv(p);

#ifdef DEBUG
printf("\nBefore matrix opt:\n");
dump_xfit(p);
#endif
		/* Optimise matrix on its own */
		p->opt_ssch = 0;
		p->opt_ch = -1;
		p->opt_msk = oc_m;
		setup_xfit(p, p->wv, p->sa, 0.0, 0.5); 

#ifdef NODDV
		if (powell(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
		                                 xfitfunc, (void *)p, xfitprog, (void *)p) != 0)
			warning("xfit_fit: Powell failed to converge, residual error = %f",rerr);
#else
		if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
		                       xfitfunc, dxfitfunc, (void *)p, xfitprog, (void *)p) != 0)
			warning("xfit_fit: Conjgrad failed to converge, residual error = %f", rerr);
#endif
		for (i = 0; i < p->opt_cnt; i++)		/* Copy optimised values back */
			p->v[p->opt_off + i] = p->wv[i];

#ifdef DEBUG
printf("\nAfter matrix opt:\n");
dump_xfit(p);

#endif
	}

	/* Optimise input and matrix together */
	if ((p->tcomb & oc_im) == oc_im) {
		double rerr;
		int sm_iluord = p->sm_iluord;

		if (p->verb)
			printf("About to optimise a common ord 0 input curve and matrix\n");

		/* Setup pseudo-inverse if we need it */
		if (p->flags & XFIT_FM_INPUT)
			setup_piv(p);

		p->opt_ssch = 1;
		p->sm_iluord = 1;		/* Do a single order for first up */
		p->opt_ch = -1;
		p->opt_msk = oc_im;
		setup_xfit(p, p->wv, p->sa, 0.5, 0.3); 

#ifdef NODDV
		if (powell(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
		                                xfitfunc, (void *)p, xfitprog, (void *)p) != 0) {
#ifdef DEBUG
			warning("xfit_fit: Powell failed to converge, residual error = %f",rerr);
#endif
		}
#else
		if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
		                     xfitfunc, dxfitfunc, (void *)p, xfitprog, (void *)p) != 0) {
#ifdef DEBUG
			warning("xfit_fit: Conjgrad failed to converge, residual error = %f",rerr);
#endif
		}
#endif	/* !NODDV */
		for (e = 0; e < di; e++) {				/* Copy optimised values back */
			for (i = 0; i < p->sm_iluord; i++)
				p->v[p->shp_offs[e] + i] = p->wv[i];
			for (; i < p->iluord[e]; i++)
				p->v[p->shp_offs[e] + i] = 0.0;
		}
		for (i = p->sm_iluord; i < p->opt_cnt; i++) 
			p->v[p->mat_off + i - p->sm_iluord] = p->wv[i];
#ifdef DEBUG
printf("\nAfter single input and matrix opt:\n");
dump_xfit(p);
#endif

		/* - - - - - - - - - - - */
		if (p->verb)
			printf("About to optimise a common input curve and matrix\n");

		/* Setup pseudo-inverse if we need it */
		if (p->flags & XFIT_FM_INPUT)
			setup_piv(p);

		p->opt_ssch = 1;
		p->sm_iluord = sm_iluord;		/* restore this */
		p->opt_ch = -1;
		p->opt_msk = oc_im;
		setup_xfit(p, p->wv, p->sa, 0.5, 0.3); 

#ifdef NODDV
		if (powell(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
		                                xfitfunc, (void *)p, xfitprog, (void *)p) != 0) {
#ifdef DEBUG
			warning("xfit_fit: Powell failed to converge, residual error = %f",rerr);
#endif
		}
#else
		if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
		                     xfitfunc, dxfitfunc, (void *)p, xfitprog, (void *)p) != 0) {
#ifdef DEBUG
			warning("xfit_fit: Conjgrad failed to converge, residual error = %f",rerr);
#endif
		}
#endif	/* !NODDV */
		for (e = 0; e < di; e++) {				/* Copy optimised values back */
			for (i = 0; i < p->sm_iluord; i++)
				p->v[p->shp_offs[e] + i] = p->wv[i];
			for (; i < p->iluord[e]; i++)
				p->v[p->shp_offs[e] + i] = 0.0;
		}
		for (i = p->sm_iluord; i < p->opt_cnt; i++) 
			p->v[p->mat_off + i - p->sm_iluord] = p->wv[i];
#ifdef DEBUG
printf("\nAfter single input and matrix opt:\n");
dump_xfit(p);
#endif

		/* - - - - - - - - - - - */
		if (p->verb)
			printf("About to optimise input curves and matrix\n");

		if ((p->tcomb & oc_mo) != oc_mo) {	/* If this will be last fit */
			powtol = POWTOL;
			maxits = MAXITS;
		}

		/* Setup pseudo-inverse if we need it */
		if (p->flags & XFIT_FM_INPUT)
			setup_piv(p);

		p->opt_ssch = 0;
		p->opt_ch = -1;
		p->opt_msk = oc_im;
		setup_xfit(p, p->wv, p->sa, 0.5, 0.3); 
		/* Suppress the warnings the first time through - it's better to cut off the */
		/* itterations and move on to the output curve, and worry about it not */
		/* converging the second time through. */
#ifdef NODDV
		if (powell(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
		                                xfitfunc, (void *)p, xfitprog, (void *)p) != 0) {
#ifdef DEBUG
			warning("xfit_fit: Powell failed to converge, residual error = %f",rerr);
#endif
		}
#else
		if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
		                     xfitfunc, dxfitfunc, (void *)p, xfitprog, (void *)p) != 0) {
#ifdef DEBUG
			warning("xfit_fit: Conjgrad failed to converge, residual error = %f",rerr);
#endif
		}
#endif	/* !NODDV */
		for (i = 0; i < p->opt_cnt; i++)		/* Copy optimised values back */
			p->v[p->opt_off + i] = p->wv[i];
#ifdef DEBUG
printf("\nAfter input and matrix opt:\n");
dump_xfit(p);
#endif
	}

	/* Optimise the matrix and output curves together */
	if ((p->tcomb & oc_mo) == oc_mo) {
		double rerr;

		if (p->verb)
			printf("About to optimise output curves and matrix\n");

		if ((p->tcomb & oc_im) != oc_im) {	/* If this will be last fit */
			powtol = POWTOL;
			maxits = MAXITS;
		}

		/* Setup pseudo-inverse if we need it */
		if (p->flags & XFIT_FM_INPUT)
			setup_piv(p);

		p->opt_ssch = 0;
		p->opt_ch = -1;
		p->opt_msk = oc_mo;
		setup_xfit(p, p->wv, p->sa, 0.3, 0.3); 
#ifdef NODDV
		if (powell(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
		                                    xfitfunc, (void *)p, xfitprog, (void *)p) != 0)
			warning("xfit_fit: Powell failed to converge, residual error = %f",rerr);
#else
		if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits, xfitfunc,
		                                   dxfitfunc, (void *)p, xfitprog, (void *)p) != 0)
			warning("xfit_fit: Conjgrad failed to converge, residual error = %f",rerr);
#endif
		for (i = 0; i < p->opt_cnt; i++)		/* Copy optimised values back */
			p->v[p->opt_off + i] = p->wv[i];
#ifdef DEBUG
printf("\nAfter output opt:\n");
dump_xfit(p);
#endif

		/* Optimise input and matrix together again, after altering matrix */
		if ((p->tcomb & oc_im) == oc_im) {

			if (p->verb)
				printf("About to optimise input curves and matrix again\n");


#ifndef NODDV
			if ((p->tcomb & oc_imo) != oc_imo)	/* If this will be last fit */
#endif
			{
				powtol = POWTOL;
				maxits = MAXITS;
			}

			/* Setup pseudo-inverse if we need it */
			if (p->flags & XFIT_FM_INPUT)
				setup_piv(p);

			p->opt_ssch = 0;
			p->opt_ch = -1;
			p->opt_msk = oc_im;
			setup_xfit(p, p->wv, p->sa, 0.2, 0.2); 
#ifdef NODDV
			if (powell(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
			                               xfitfunc, (void *)p, xfitprog, (void *)p) != 0)
				warning("xfit_fit: Powell failed to converge, residual error = %f",rerr);
#else
			if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
			                           xfitfunc, dxfitfunc, (void *)p, xfitprog, (void *)p) != 0)
				warning("xfit_fit: Conjgrad failed to converge, residual error = %f",rerr);
#endif
			for (i = 0; i < p->opt_cnt; i++)		/* Copy optimised values back */
				p->v[p->opt_off + i] = p->wv[i];
		}
#ifdef DEBUG
printf("\nAfter 2nd input and matrix opt:\n");
dump_xfit(p);
#endif

#ifndef NODDV
		/* Optimise all together */
		/* (This is very slow using powell) */
		if ((p->tcomb & oc_imo) == oc_imo) {

			if (p->verb)
				printf("About to optimise input, matrix and output together\n");

			/* Setup pseudo-inverse if we need it */
			if (p->flags & XFIT_FM_INPUT)
				setup_piv(p);

			p->opt_ssch = 0;
			p->opt_ch = -1;
			p->opt_msk = oc_imo;
			setup_xfit(p, p->wv, p->sa, 0.1, 0.1); 
			if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, POWTOL, MAXITS,
			                    xfitfunc, dxfitfunc, (void *)p, xfitprog, (void *)p) != 0)
				warning("xfit_fit: Conjgrad failed to converge, residual error = %f",rerr);
			for (i = 0; i < p->opt_cnt; i++)		/* Copy optimised values back */
				p->v[p->opt_off + i] = p->wv[i];

			/* Setup final pseudo-inverse from shaper/matrix/out */
			if (p->flags & XFIT_FM_INPUT)
				setup_piv(p);
		}

#ifdef DEBUG
printf("\nAfter all together opt:\n");
dump_xfit(p);
#endif

#endif /* !NODDV */

		/* Adjust output curve white point. */
		/* This is for the benefit of the B2A table */
		if (p->flags & XFIT_OUT_ZERO) {

			if (p->verb)
				printf("About to adjust a and b output curves for white point\n");

			for (f = 1; f < 3 && f < p->fdi; f++) {
				p->opt_ch = f;
				p->wv[0] = p->v[p->out_offs[f]];	/* Current parameter value */
				p->sa[0] = 0.1;						/* Search radius */
				if (powell(&rerr, 1, p->wv, p->sa, 0.0000001, 1000,
				                                   symoptfunc, (void *)p, NULL, NULL) != 0)
					error("xfit_fit: Powell failed to converge, residual error = %f",rerr);
				p->v[p->out_offs[f]] = p->wv[0];	/* Copy results back */
			}
		}
	}

	/* In case we don't generate position curves, */
	/* copy the input curves to the position, so that */
	/* ipos is computed correctly */
	{
		double *bb;
		
		b = p->v + p->shp_off;
		bb = p->v + p->pos_off;
		for (e = 0; e < di; b += p->iluord[e], bb += p->iluord[e], e++) {
			for (i = 0; i < p->iluord[e]; i++) {
				bb[i] = b[i];
			}
		}
	}


	/* If we want position curves, generate them */
	/* (This could possibly be improved by using some sort */
	/* of optimization drivel approach rather than the predictive */
	/* method used here.) */
	if (p->tcomb & oc_p) {
		int ee;

		if (p->verb)
			printf("About to create grid position input curves\n");

		/* Allocate in->rout duplicate point set */
		if (p->rpoints == NULL) {
			if ((p->rpoints = (cow *)malloc(p->nodp * sizeof(cow))) == NULL)
				return 1;
		}

		/* Create a set of 1D rspl setup points that contains */
		/* the residual error */
		for (i = 0; i < p->nodp; i++) {
			double tv[MXDO];		/* Target output value */
			double mv[MXDO];		/* Model output value */
			double ev;

			xfit_abs_to_rel(p, tv, p->ipoints[i].v);
			xfit_shmatsh(p, mv, p->ipoints[i].p);

			/* Evaluate the error squared */
			if (p->flags & XFIT_FM_INPUT) {
				double pp[MXDI];
				for (e = 0; e < di; e++)
					pp[e] = p->ipoints[i].p[e];
				for (f = 0; f < fdi; f++) {
					double t1 = tv[f] - mv[f];	/* Error in output */
		
					/* Create input point offset by equivalent delta to output point */
					/* error, in proportion to the partial derivatives for that output. */
					for (e = 0; e < di; e++)
						pp[e] += t1 * p->piv[i].ide[f][e];
				}
				ev = p->to_de2(p->cntx2, pp, p->ipoints[i].p);
			} else {
				ev = p->to_de2(p->cntx2, mv, tv);
			}
			p->rpoints[i].v[0] = ev;
			p->rpoints[i].w    = p->ipoints[i].w;
		}

		/* Do each input axis in turn */
		for (ee = 0; ee < p->di; ee++) {
			rspl *resid;
			double imin[1],imax[1],omin[1],omax[1]; 
			int resres[1] = { 1024 };
#define NPGP 100
			mcv *posc;
			mcvco pgp[NPGP];
			double vo, vs;
			double *pms;
			
			/* Create a rspl that plots the residual error */
			/* vs the axis value */
			for (i = 0; i < p->nodp; i++)
				p->rpoints[i].p[0] = p->ipoints[i].p[ee];

			imin[0] = in_min[ee];
			imax[0] = in_max[ee];
			omin[0] = 0.0;
			omax[0] = 0.0;
			
			if ((resid = new_rspl(RSPL_NOFLAGS, 1, 1)) == NULL)
				return 1;

			resid->fit_rspl_w(resid, RSPLFLAGS, p->rpoints, p->nodp, imin, imax, resres,
				omin, omax, 2.0, NULL, NULL);

#ifdef DEBUG_PLOT
			{
#define	XRES 100
				double xx[XRES];
				double y1[XRES];

				printf("Input residual error channel %d\n",ee);
				for (i = 0; i < XRES; i++) {
					co pp;
					double x;
					x = i/(double)(XRES-1);
					xx[i] = x = x * (imax[0] - imin[0]) + imin[0];
					pp.p[0] = xx[i];
					resid->interp(resid, &pp);
					y1[i] = pp.v[0];
					if (y1[i] < 0.0)
						y1[i] = 0.0;
					y1[i] = pow(y1[i], 0.5);	/* Convert from error^2 to error */
				}
				do_plot(xx,y1,NULL,NULL,XRES);
			}
#endif /* DEBUG_PLOT */

			/* Create a set of guide points that contain */
			/* the accumulated residual error vs. the axis value */
			for (i = 0; i < NPGP; i++) {
				co pp;
				double vv;

				pp.p[0] = i/(NPGP-1.0);
				resid->interp(resid, &pp);

				pgp[i].p = (pp.p[0] - in_min[ee])/(in_max[ee] - in_min[ee]);
				pgp[i].w = 1.0; 

				vv = pp.v[0];
				if (vv < 0.0)
					vv = 0.0;
				vv = pow(vv, 0.5);		/* Convert from error^2 to error */
				vv += PSHAPE_MINE;		/* In case error is near zero */
				vv = pow(vv, PSHAPE_DIST);		/* Agressivness of grid distribution */

				if (i == 0)
					pgp[i].v = vv;
				else
					pgp[i].v = pgp[i-1].v + vv;

			}
			resid->del(resid);

			/* Normalize the output range */
			vo = pgp[0].v;
			vs = pgp[NPGP-1].v - vo;
			
			for (i = 0; i < NPGP; i++) {
				pgp[i].v = (pgp[i].v - vo)/vs;
				
				/* Apply any dark emphasis */
				if (demph > 1.0) {
					pgp[i].v = icx_powlike(pgp[i].v, 1.0/demph);
				}
			}
			/* Fit the non-monotonic parameters to the guide points */
			if ((posc = new_mcv_noos()) == NULL)
				return 1;

			posc->fit(posc, 0, p->iluord[ee], pgp, NPGP, 0.1);	// ~~99

#ifdef DEBUG_PLOT
			{
#define	XRES 100
				double xx[XRES];
				double y1[XRES];

				printf("Position curve %d\n",ee);
				for (i = 0; i < XRES; i++) {
					xx[i] = i/(double)(XRES-1);
					y1[i] = posc->interp(posc, xx[i]);
				}
				do_plot(xx,y1,NULL,NULL,XRES);
			}
#endif /* DEBUG_PLOT */

			/* Transfer parameters to xfit pos (skip offset and scale) */
			posc->get_params(posc, &pms);

			for (i = 0; i < p->iluord[ee]; i++) {
				p->v[p->pos_offs[ee] + i] = pms[i+2];
			}
//p->v[p->in_offs[ee]] = -1.5;

			free(pms);
			posc->del(posc);
		}
	}

#ifdef DEBUG
printf("Final parameters:\n");
dump_xfit(p);
#endif

#ifdef DEBUG_PLOT
	{
#define	XRES 100
		double xx[XRES];
		double y1[XRES];

		for (e = 0; e < p->di; e++) {
			printf("Input position curve channel %d\n",e);
			for (i = 0; i < XRES; i++) {
				double x;
				x = i/(double)(XRES-1);
				xx[i] = x = x * (p->in_max[e] - p->in_min[e]) + p->in_min[e];
				y1[i] = xfit_poscurve(p, x, e);
			}
			do_plot(xx,y1,NULL,NULL,XRES);
		}

		for (e = 0; e < p->di; e++) {
			printf("Input shape curve channel %d\n",e);
			for (i = 0; i < XRES; i++) {
				double x;
				x = i/(double)(XRES-1);
				xx[i] = x = x * (p->in_max[e] - p->in_min[e]) + p->in_min[e];
				y1[i] = xfit_shpcurve(p, x, e);
			}
			do_plot(xx,y1,NULL,NULL,XRES);
		}
	
		for (e = 0; e < p->di; e++) {
			printf("Combined input curve channel %d\n",e);
			for (i = 0; i < XRES; i++) {
				double x;
				x = i/(double)(XRES-1);
				xx[i] = x = x * (p->in_max[e] - p->in_min[e]) + p->in_min[e];
				y1[i] = p->incurve(p, x, e);
			}
			do_plot(xx,y1,NULL,NULL,XRES);
		}
	
		for (f = 0; f < p->fdi; f++) {
			printf("Output curve channel %d\n",f);
			for (i = 0; i < XRES; i++) {
				double x;
				x = i/(double)(XRES-1);
				xx[i] = x = x * (p->out_max[f] - p->out_min[f]) + p->out_min[f];
				y1[i] = p->outcurve(p, x, f);
			}
			do_plot(xx,y1,NULL,NULL,XRES);
		}
	}
#endif /* DEBUG_PLOT */

	/* Create final clut rspl using the established pos/shape/output curves */
	/* and white point */ 
	if (flags & XFIT_MAKE_CLUT) {
		double *ipos[MXDI];

		/* Create an in' -> rout' scattered test point set */
		if (p->rpoints == NULL) {
			if ((p->rpoints = (cow *)malloc(p->nodp * sizeof(cow))) == NULL)
				return 1;
		}
		for (i = 0; i < p->nodp; i++) {
			p->rpoints[i].w = p->ipoints[i].w;
			xfit_inpscurves(p, p->rpoints[i].p, p->ipoints[i].p);
			for (f = 0; f < fdi; f++)
				p->rpoints[i].v[f] = p->ipoints[i].v[f];
			xfit_abs_to_rel(p, p->rpoints[i].v, p->rpoints[i].v);
			xfit_invoutcurves(p, p->rpoints[i].v, p->rpoints[i].v);

			if (p->skm) {
				/* Look up the skeleton value */
				double skval[MXDO];
				p->skm->lookup(p->skm, skval, p->ipoints[i].p);
				xfit_abs_to_rel(p, skval, skval);
				xfit_invoutcurves(p, skval, skval);
//printf("~1 point %d at %f %f %f, targ %f %f %f skm %f %f %f\n", i,p->ipoints[i].p[0],p->ipoints[i].p[1],p->ipoints[i].p[2], p->rpoints[i].v[0],p->rpoints[i].v[1],p->rpoints[i].v[2], skval[0], skval[1], skval[2]);
				/* Subtract it from value at this point, */
				/* so rspl will fit difference to skeleton model */
				for (f = 0; f < fdi; f++)
					p->rpoints[i].v[f] -= skval[f]; 
			}
//printf("~1 point %d, w %f, %f %f %f %f -> %f %f %f\n",i,p->rpoints[i].w,p->rpoints[i].p[0], p->rpoints[i].p[1], p->rpoints[i].p[2], p->rpoints[i].p[3],p->rpoints[i].v[0], p->rpoints[i].v[1], p->rpoints[i].v[2]);
		}

		/* Create ipos[] arrays, that hold the shaper space */
		/* grid position due to the positioning curves. */
		/* This tells the rspl scattered data interpolator */
		/* the grid spacing that smoothness should be */
		/* measured against. */
#ifdef DEBUG
printf("~1 about to setup ipos\n");
#endif
		for (e = 0; e < p->di; e++) {
//printf("~1 e = %d\n",e);
			if ((ipos[e] = (double *)malloc((p->gres[e]) * sizeof(double))) == NULL)
				return 1;
//printf("~1 about to do %d spans\n",p->gres[e]);
			for (i = 0; i < p->gres[e]; i++) {
				double cv;
				cv = (double)i/p->gres[e];

//printf("~1 i = %d, pos space = %f\n",i,cv);
				/* Inverse lookup grid position through positioning curve */
				/* to give device space type value */
				cv = icxInvTransFunc(p->v + p->pos_offs[e], p->iluord[e], cv); 

//printf("~1 dev space = %f\n",cv);
				/* Forward lookup device type value through the shaper curve */
				/* to give value in shape linearized space. */
				cv = icxTransFunc(p->v + p->shp_offs[e], p->iluord[e], cv); 
//printf("~1 shape space = %f\n",cv);
				ipos[e][i] = cv;
#ifdef DEBUG
printf("~1 ipos[%d][%d] = %f\n",e,i,cv);
#endif
			}
		}

		if (p->clut != NULL)
			p->clut->del(p->clut);
		if ((p->clut = new_rspl(RSPL_NOFLAGS, di, fdi)) == NULL)
			return 1;

		if (p->verb)
			printf("Create final clut from scattered data\n");

#ifdef EXTEND_GRID
#define XN EXTEND_GRID_BYN
		/* Try increasing the grid by one row all around */
		{
#pragma message("!!!!!!!!!!!! Experimental rspl fitting resolution !!!!!!!!!")
			double xin_min[MXDI];
			double xin_max[MXDI];
			int xgres[MXDI];
			double del;
			double *xipos[MXDI];

			for (e = 0; e < p->di; e++) {
				del = (in_max[e] - in_min[e])/(gres[e]-1.0);		/* Extension */
				xin_min[e] = in_min[e] - XN * del;
				xin_max[e] = in_max[e] + XN * del;
				xgres[e] = gres[e] + 2 * XN;
//printf("~1 xgres %d, gres %d\n",xgres[e], gres[e]);

				if ((xipos[e] = (double *)malloc((xgres[e]) * sizeof(double))) == NULL)
					return 1;
				for (i = 0; i < xgres[e]; i++) {
					if (i < XN) {					/* Extrapolate bottom */
						xipos[e][i] = ipos[e][0] - (XN - i) * (ipos[e][1] - ipos[e][0]);
//printf("~1 xipos[%d] %f from ipos[%d] %f and ipos[%d] %f\n",i,xipos[e][i],0,ipos[e][0],1,ipos[e][1]);
					} else if (i >= (xgres[e]-XN)) {	/* Extrapolate top */
						xipos[e][i] = ipos[e][gres[e]-1] + (i - xgres[e] + XN + 1) * (ipos[e][gres[e]-1] - ipos[e][gres[e]-2]);
//printf("~1 xipos[%d] %f from ipos[%d] %f and ipos[%d] %f\n",i,xipos[e][i],gres[e]-1,ipos[e][gres[e]-1],gres[e]-2,ipos[e][gres[e]-2]);
					} else {
						xipos[e][i] = ipos[e][i-XN];
//printf("~1 xipos[%d] %f from ipos[%d] %f\n",i,xipos[e][i],i-XN,ipos[e][i-XN]);
					}
				}
			}

			p->clut->fit_rspl_w(p->clut, rsplflags, p->rpoints, p->nodp, xin_min, xin_max, xgres,
				out_min, out_max, smooth, oavgdev, xipos);

			for (e = 0; e < p->di; e++) {
				free(xipos[e]);
			}
		}
#undef XN
#else
//		if (p->skm) {
//			/* This doesn't seem to work as well as some explicit neutral axis points.. */
//			p->clut->fit_rspl_w_df(p->clut, rsplflags, p->rpoints, p->nodp, in_min, in_max, gres,
//				out_min, out_max, smooth, oavgdev, ipos, 1.0, (void *)p, skm_weak);
//		} else
		/* Normal multi-d scattered point fitting */
		p->clut->fit_rspl_w(p->clut, rsplflags, p->rpoints, p->nodp, in_min, in_max, gres,
				out_min, out_max, smooth, oavgdev, ipos);
#endif
		if (p->verb)
			printf("\n");

		for (e = 0; e < p->di; e++)
			free(ipos[e]);

		/* If we used a skeleton model, add it into the resulting rspl values */
		if (p->skm) {

			/* Undo the input point change to allow diagnostic code to work */
			for (i = 0; i < p->nodp; i++) {
				/* Look up the skeleton value */
				double skval[MXDO];
				p->skm->lookup(p->skm, skval, p->ipoints[i].p);
				xfit_abs_to_rel(p, skval, skval);
				xfit_invoutcurves(p, skval, skval);

				/* Subtract it from value at this point, */
				/* so rspl will fit difference to skeleton model */
				for (f = 0; f < fdi; f++)
					p->rpoints[i].v[f] += skval[f]; 
			}

			/* Undo the skm from the resultant rspl */
			p->clut->re_set_rspl(
				p->clut,		/* this */
				0,				/* Combination of flags */
				(void *)p,		/* Opaque function context */
				skm_rspl_out	/* Function to set from */
			);

		}

		/* The overall device to absolute conversion is now what we want */
		/* (as dictated by the points, weighting and best fit), */
		/* but we need to adjust the device to relative conversion */
		/* to make device white map exactly to D50, without touching */
		/* the overall absolute behaviour. */
		if (p->flags & XFIT_OUT_WP_REL) {
			co wcc;						/* device white + aprox rel. white */
			icmXYZNumber _wp;			/* Uncorrected dw maps to _wp */ 

			if (p->verb)
				printf("Doing White point fine tune:\n");
				
			/* See what the relative and absolute white point has turned out to be, */
			/* by looking up the device white in the current conversion */
			xfit_inpscurves(p, wcc.p, dw);
			p->clut->interp(p->clut, &wcc);
			xfit_outcurves(p, wcc.v, wcc.v);
			if (p->flags & XFIT_OUT_LAB)
				icmLab2XYZ(&icmD50, wcc.v, wcc.v);

			if (p->verb) {
				double labwp[3];
				icmXYZ2Lab(&icmD50, labwp, wcc.v);
				printf("Before fine tune, rel WP = XYZ %f %f %f, Lab %f %f %f\n",
				wcc.v[0], wcc.v[1],wcc.v[2], labwp[0], labwp[1], labwp[2]);
			}

			/* Matrix needed to correct approx rel wp to target D50 */
			icmAry2XYZ(_wp, wcc.v);		/* Aprox relative target white point */
			if (p->picc != NULL)	/* Correction */
				p->picc->chromAdaptMatrix(p->picc, ICM_CAM_NONE, NULL, p->cmat, icmD50, _wp);
			else
				icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, _wp, p->cmat);

			/* Compute the actual white point, and return it to caller */
			icmMulBy3x3(wp, p->toAbs, wcc.v);
			icmAry2Ary(p->wp, wp);

			/* Apply correction to fine tune rspl data. */
			/* NOTE: this doesn't always give us a perfect D50 white for */
			/* Lab PCS input profiles because the dev white may land */
			/* within a cell, and the clipping of Lab PCS values in the grid */
			/* may introduce errors in the interpolated value. */
			p->clut->re_set_rspl(
				p->clut,		/* this */
				0,				/* Combination of flags */
				(void *)p,		/* Opaque function context */
				conv_rspl_out	/* Function to set from */
			);

			/* Fix absolute conversions to leave absolute response unchanged. */
			icmAry2XYZ(_wp, wp);		/* Actual white point */
			if (p->picc != NULL) {
				p->picc->chromAdaptMatrix(p->picc, ICM_CAM_NONE, p->toAbs, p->fromAbs, icmD50, _wp);
			} else {
				icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, _wp, p->fromAbs);
				icmChromAdaptMatrix(ICM_CAM_BRADFORD, _wp, icmD50, p->toAbs);
			}

			if (p->verb) {
				double labwp[3];

				/* Lookup white again */
				xfit_inpscurves(p, wcc.p, dw);
				p->clut->interp(p->clut, &wcc);
				xfit_outcurves(p, wcc.v, wcc.v);
				if (p->flags & XFIT_OUT_LAB)
					icmLab2XYZ(&icmD50, wcc.v, wcc.v);
				icmXYZ2Lab(&icmD50, labwp, wcc.v);
				printf("After fine tune, rel WP = XYZ %f %f %f, Lab %f %f %f\n",
				wcc.v[0], wcc.v[1], wcc.v[2], labwp[0], labwp[1], labwp[2]);
				printf("                 abs WP = XYZ %s, Lab %s\n", icmPdv(3, wp), icmPLab(wp));
			}
		}

		/* Create default wpscale */
		if (wpscale < 0.0) {
			wpscale = 1.0;
		} else {
			if (p->verb) {
				printf("White manual point scale %f\n", wpscale);
			}
		}

		/* If we are going to auto scale the WP to avoid clipping */
		/* cLUT values above the WP: */
		if ((p->flags & XFIT_OUT_WP_REL_US) == XFIT_OUT_WP_REL_US) {
			co wcc;
			double bw[3];
			icmXYZNumber _wp;
			double uswpscale = 1.0;
			double mxd, mxY;
			double ndw[3];

			/* See what device space gamut boundary white (ie. 1,1,1) maps to */
			xfit_inpscurves(p, wcc.p, dgw);
			p->clut->interp(p->clut, &wcc);
			xfit_outcurves(p, wcc.v, wcc.v);
			if (p->flags & XFIT_OUT_LAB)
				icmLab2XYZ(&icmD50, wcc.v, wcc.v);
			icmMulBy3x3(wcc.v, p->toAbs, wcc.v);	/* Convert to absolute */

			mxY = wcc.v[1];
			icmCpy3(bw, wcc.v);
//printf("~1 1,1,1 Y = %f\n",wcc.v[1]);

			/* See what the device white point value scaled to 1 produces */
			mxd = -1.0;
			for (e = 0; e < p->di; e++) {
				if (dw[e] > mxd)
					mxd = dw[e];
			}
			for (e = 0; e < p->di; e++)
				ndw[e] = dw[e]/mxd;

			xfit_inpscurves(p, wcc.p, ndw);
			p->clut->interp(p->clut, &wcc);
			xfit_outcurves(p, wcc.v, wcc.v);
			if (p->flags & XFIT_OUT_LAB)
				icmLab2XYZ(&icmD50, wcc.v, wcc.v);
			icmMulBy3x3(wcc.v, p->toAbs, wcc.v);	/* Convert to absolute */

//printf("~1 ndw = %f %f %f Y = %f\n",ndw[0],ndw[1],ndw[2],wcc.v[1]);
			if (wcc.v[1] > mxY) {
				mxY = wcc.v[1];
				icmCpy3(bw, wcc.v);
			}

			/* Compute WP scale factor needed to fit mxY */
			if (mxY > wp[1]) {
				uswpscale = mxY/wp[1];
				wpscale *= uswpscale;
				if (p->verb) {
					printf("Dev boundary white XYZ %s, scale WP by %f, total WP scale %f\n",
					icmPdv(3, bw), uswpscale, wpscale);
				}
			}
		}

		/* If the scaled WP would have Y > 1.0, clip it to 1.0 */
		if (p->flags & XFIT_CLIP_WP) {

			if ((wp[1] * wpscale) > 1.0) {
				wpscale = 1.0/wp[1];		/* Make wp Y = 1.0 */		
				if (p->verb) {
					printf("WP Y would ve > 1.0. scale by %f to clip it\n",wpscale);
				}
			}
		}

		/* Apply our total wp scale factor */
		if (wpscale != 1.0) {
			icmXYZNumber _wp;

			/* Create inverse scaling matrix for relative rspl data */
			icmSetUnity3x3(p->cmat);
			icmScale3x3(p->cmat, p->cmat, 1.0/wpscale);

			/* Inverse scale the rspl */
			p->clut->re_set_rspl(
				p->clut,		/* this */
				0,				/* Combination of flags */
				(void *)p,		/* Opaque function context */
				conv_rspl_out	/* Function to set from */
			);

			/* Scale the WP */
			icmScale3(wp, wp, wpscale);

			/* return scaled white point to caller */
			icmAry2Ary(p->wp, wp);

			/* Fix absolute conversions to leave absolute response unchanged. */
			icmAry2XYZ(_wp, wp);		/* Actual white point */
			icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, _wp, p->fromAbs);
			icmChromAdaptMatrix(ICM_CAM_BRADFORD, _wp, icmD50, p->toAbs);
		}

		/* Clip any values in the grid over D50 L to D50 */
		if ((p->flags & XFIT_OUT_WP_REL_C) == XFIT_OUT_WP_REL_C) {

			/* Compute the rspl D50 value to avoid calc in clip_rspl_out() */
			if (p->flags & XFIT_OUT_LAB) {
				p->cmat[0][0] = 100.0;
				p->cmat[0][1] = 0.0;
				p->cmat[0][2] = 0.0;
			} else {
				icmXYZ2Ary(p->cmat[0], icmD50);
			}
			xfit_invoutcurves(p, p->cmat[0], p->cmat[0]);

			if (p->verb)
				printf("Clipping any cLUT grid points with Y > 1 to D50\n");

			p->clut->re_set_rspl(
				p->clut,		/* this */
				0,				/* Combination of flags */
				(void *)p,		/* Opaque function context */
				clip_rspl_out	/* Function to set from */
			);
		}

		/* Force black point to given value */
//		if (bpo != NULL) {
//			co tv;
//			int rv;
//
//			xfit_inpscurves(p, tv.p, bpo->p);
//
//			xfit_XYZ_abs_to_rel(p, tv.v, bpo->v);
//			xfit_invoutcurves(p, tv.v, tv.v);
//printf("~1 xfit: fine after curves  black at %f %f %f to %f %f %f\n",
//tv.p[0], tv.p[1], tv.p[2], tv.v[0], tv.v[1], tv.v[2]);
//			rv = p->clut->tune_value(p->clut, &tv);
//			if (rv != 0)
//				warning("Black Point Override failed - clipping");
//		}

#ifdef SPECIAL_FORCE
		/* Replace the rspl nodes with ones directly computed */
		/* from the synthetic linear RGB->XYZ model */
		{
			double twp[3];
			icmXYZNumber _wp;

			/* See what current device white maps to */
			twp[0] = twp[1] = twp[2] = 1.0;
			domodel(twp, twp);

			icmAry2XYZ(_wp, twp);

			/* Matrix needed to correct to D50 */
			icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, _wp, p->cmat);
	
			p->clut->re_set_rspl(
				p->clut,		/* this */
				0,				/* Combination of flags */
				(void *)p,		/* Opaque function context */
				set_rspl_out1	/* Function to set from */
			);
		}
#endif /* SPECIAL_FORCE */

	/* Evaluate the residual error now, with the rspl in place */
#ifdef DEBUG_PLOT
	{
		int ee;
		double maxe = 0.0;
		double avee = 0.0;

		/* Allocate in->rout duplicate point set */
		if (p->rpoints == NULL) {
			if ((p->rpoints = (cow *)malloc(p->nodp * sizeof(cow))) == NULL)
				return 1;
		}

		/* Create a set of 1D rspl setup points that contains */
		/* the residual error */
		for (i = 0; i < p->nodp; i++) {
			double tv[MXDO];		/* Target output value */
			double mv[MXDO];		/* Model output value */
			co pp;
			double ev;

			xfit_abs_to_rel(p, tv, p->ipoints[i].v);

			for (e = 0; e < p->di; e++)
				pp.p[e] = p->ipoints[i].p[e];

			xfit_inpscurves(p, pp.p, pp.p);
			p->clut->interp(p->clut, &pp);
			xfit_outcurves(p, pp.v, pp.v);

			for (f = 0; f < p->fdi; f++)
				mv[f] = pp.v[f];

			/* Evaluate the residual error suqared */
			if (p->flags & XFIT_FM_INPUT) {
				double pp[MXDI];
				for (e = 0; e < di; e++)
					pp[e] = p->ipoints[i].p[e];
				for (f = 0; f < fdi; f++) {
					double t1 = tv[f] - mv[f];	/* Error in output */
		
					/* Create input point offset by equivalent delta to output point */
					/* error, in proportion to the partial derivatives for that output. */
					for (e = 0; e < di; e++)
						pp[e] += t1 * p->piv[i].ide[f][e];
				}
				ev = p->to_de2(p->cntx2, pp, p->ipoints[i].p);
			} else {
				ev = p->to_de2(p->cntx2, mv, tv);
			}
//printf("~1 point %d, loc %f %f %f %f\n",i, p->rpoints[i].p[0], p->rpoints[i].p[1], p->rpoints[i].p[2], p->rpoints[i].p[3]);
//printf("        targ %f %f %f, is %f %f %f\n", tv[0], tv[1], tv[2], mv[0], mv[1], mv[2]);

			p->rpoints[i].v[0] = ev;
			p->rpoints[i].w    = p->ipoints[i].w;

			ev = sqrt(ev);
			if (ev > maxe)
				maxe = ev;
			avee += ev;
		}
		printf("Max resid err = %f, avg err = %f\n",maxe, avee/(double)p->nodp);

		/* Evaluate each input axis in turn */
		for (ee = 0; ee < p->di; ee++) {
			rspl *resid;
			double imin[1],imax[1],omin[1],omax[1]; 
			int resres[1] = { 1024 };
			
			/* Create a rspl that gives the residual error squared */
			/* vs. the axis value */
			for (i = 0; i < p->nodp; i++)
				p->rpoints[i].p[0] = p->ipoints[i].p[ee];

			imin[0] = in_min[ee];
			imax[0] = in_max[ee];
			omin[0] = 0.0;
			omax[0] = 0.0;
			
			if ((resid = new_rspl(RSPL_NOFLAGS, 1, 1)) == NULL)
				return 1;

			resid->fit_rspl_w(resid, RSPLFLAGS, p->rpoints, p->nodp, imin, imax, resres,
				omin, omax, 2.0, NULL, NULL);

			{
#define	XRES 100
				double xx[XRES];
				double y1[XRES];

				printf("Finale residual error vs. input channel %d\n",ee);
				for (i = 0; i < XRES; i++) {
					co pp;
					double x;
					x = i/(double)(XRES-1);
					xx[i] = x = x * (imax[0] - imin[0]) + imin[0];
					pp.p[0] = xx[i];
					resid->interp(resid, &pp);
					y1[i] = sqrt(fabs(pp.v[0]));
				}
				do_plot(xx,y1,NULL,NULL,XRES);
			}
			resid->del(resid);
		}
	}
#endif /* DEBUG_PLOT */

	/* Special test code to figure out what's wrong with position curves */
#ifdef NEVER
	{
		double rgb[3], xyz[3];
		co pp;

extern int rspldb;
db = 1;
rspldb = 1;
		printf("~1 gres = %d %d %d\n", gres[0], gres[1], gres[2]);
		for (i = 0; i < 3; i++) {

			if (i == 0) {
				/* Test point on a grid point */
				printf("\n~1 ##########################################\n");
				printf("~1 testing input at first diagonal grid point\n");
	
				rgb[0] = 1.0/(gres[0]-1.0);
				rgb[1] = 1.0/(gres[0]-1.0);
				rgb[2] = 1.0/(gres[0]-1.0);

				printf("~1 target rgb' = %f %f %f\n", rgb[0], rgb[1], rgb[2]);
				xfit_invinpscurves(p, rgb, rgb);

			} else if (i == 1) {
				/* Test point half way through a grid point */
				printf("\n~1 ##########################################\n");
				printf("\n~1 testing half way through diagonal grid point\n");
	
				rgb[0] = 0.5/(gres[0]-1.0);
				rgb[1] = 0.5/(gres[0]-1.0);
				rgb[2] = 0.5/(gres[0]-1.0);

				printf("~1 target rgb' = %f %f %f\n", rgb[0], rgb[1], rgb[2]);
				xfit_invinpscurves(p, rgb, rgb);

			} else {
				printf("\n~1 ##########################################\n");
				printf("\n~1 testing worst case point\n");
	
				rgb[0] = 0.039915;
				rgb[1] = 0.053148;
				rgb[2] = 0.230610;
			}

			pp.p[0] = rgb[0];
			pp.p[1] = rgb[1];
			pp.p[2] = rgb[2];

			printf("~1 rgb = %f %f %f\n", pp.p[0], pp.p[1], pp.p[2]);

			xfit_inpscurves(p, pp.p, pp.p);
			
			printf("~1 rgb' = %f %f %f\n", pp.p[0], pp.p[1], pp.p[2]);

			p->clut->interp(p->clut, &pp);

			printf("~1 xyz' = %f %f %f\n", pp.v[0], pp.v[1], pp.v[2]);

			xfit_outcurves(p, pp.v, pp.v);

			printf("~1 xyz = %f %f %f\n", pp.v[0], pp.v[1], pp.v[2]);

			/* Synthetic linear rgb->XYZ model */
			domodel(xyz, rgb);

			/* Apply abs->rel white point adjustment */
			icmMulBy3x3(xyz, p->mat, xyz);

			printf("~1 ref = %f %f %f, de = %f\n", xyz[0], xyz[1], xyz[2],icmXYZLabDE(&icmD50,xyz,pp.v));
		}

db = 0;
rspldb = 0;
//		exit(0);
	}
#endif /* NEVER */

	}
	return 0;
}

/* We're done with an xfit */
static void xfit_del(xfit *p) {
	if (p->v != NULL)
		free(p->v);
	if (p->wv != NULL)
		free(p->wv);
	if (p->sa != NULL)
		free(p->sa);
	if (p->rpoints != NULL)
		free(p->rpoints);
	if (p->piv != NULL)
		free(p->piv);
	if (p->uerrv != NULL)
		free(p->uerrv);
	free(p);
}

/* Create a transform fitting object */
/* return NULL on error */
xfit *new_xfit(
icc *picc			/* ICC profile used to set cone space matrix, NULL for Bradford. */
) {
	xfit *p;

	if ((p = (xfit *)calloc(1, sizeof(xfit))) == NULL) {
		return NULL;
	}

	p->picc = picc;

	/* Set method pointers */
	p->fit         = xfit_fit;
	p->incurve     = xfit_inpscurve;
	p->invincurve  = xfit_invinpscurve;
	p->outcurve    = xfit_outcurve;
	p->invoutcurve = xfit_invoutcurve;
	p->del = xfit_del;

	return p;
}