Move the solver to Colour and update all tests

1 files changed, 0 insertions, 131 deletions

diff --git a/gsoc_common.py b/gsoc_common.py
index 5d4aea0..b73d81c 100644
--- a/gsoc_common.py
+++ b/gsoc_common.py
@@ -24,137 +24,6 @@ def model_sd(ccp, primed=True):
 	grid = wvlp if primed else wvl.range()
 	return SpectralDistribution(model(grid, ccp), wvl.range(), name="Model")
 
-# The goal is to minimize the color difference between a given distrbution
-# and the one computed from the model above.
-# This function also calculates the first derivatives with respect to c's.
-def error_function(c, target, wvl, cmfs, illuminant, illuminant_XYZ):
-	U = c[0] * wvl**2 + c[1] * wvl + c[2]
-	t1 = np.sqrt(1 + U**2)
-	R = 1 / 2 + U / (2 * t1)
-
-	t2 = 1 / (2 * t1) - U**2 / (2 * t1**3)
-	dR_dc = [wvl**2 * t2, wvl * t2, t2]
-
-	E = illuminant.values * R / 100
-	dE_dc = illuminant.values * dR_dc / 100
-
-	XYZ = np.empty(3)
-	dXYZ_dc = np.empty((3, 3))
-
-	dlambda = cmfs.wavelengths[1] - cmfs.wavelengths[0]
-	for i in range(3):
-		XYZ[i] = np.dot(E, cmfs.values[:, i]) * dlambda
-		for j in range(3):
-			dXYZ_dc[i, j] = np.dot(dE_dc[j], cmfs.values[:, i]) * dlambda
-
-	# FIXME: this isn't the full CIE 1976 lightness function
-	f = (XYZ / illuminant_XYZ)**(1/3)
-
-	# FIXME: this can be vectorized
-	df_dc = np.empty((3, 3))
-	for i in range(3):
-		for j in range(3):
-			df_dc[i, j] = 1 / (3 * illuminant_XYZ[i]**(1/3)
-			                   * XYZ[i]**(2/3)) * dXYZ_dc[i, j]
-
-	Lab = np.array([
-		116 * f[1] - 16,
-		500 * (f[0] - f[1]),
-		200 * (f[1] - f[2])
-	])
-
-	dLab_dc = np.array([
-		116 * df_dc[1],
-		500 * (df_dc[0] - df_dc[1]),
-		200 * (df_dc[1] - df_dc[2])
-	])
-
-	error = np.sqrt(np.sum((Lab - target)**2))
-
-	derror_dc = np.zeros(3)
-	for i in range(3):
-		for j in range(3):
-			derror_dc[i] += dLab_dc[j, i] * (Lab[j] - target[j])
-		derror_dc[i] /= error
-
-	#print("c=[%.3g, %.3g, %.3g], XYZ=[%.3g, %.3g, %.3g], Lab=[%.3g, %.3g, %.3g], error=%g" % (*c, *XYZ, *Lab, error))
-
-	return error, derror_dc
-
-# Finds the parameters for Jakob and Hanika's model
-def jakob_hanika(target_XYZ, ill_sd, ill_xy, ccp0=(0, 0, 0), verbose=True, try_hard=True):
-	ill_sd = ill_sd.align(wvl)
-	ill_XYZ = sd_to_XYZ(ill_sd) / 100
-
-	def do_optimize(XYZ, ccp0):
-		Lab = XYZ_to_Lab(XYZ, ill_xy)
-
-		def fun(*args):
-			error, derror_dc = error_function(*args)
-			return error
-
-		def jac(*args):
-			error, derror_dc = error_function(*args)
-			return derror_dc
-
-		cmfs = STANDARD_OBSERVER_CMFS["CIE 1931 2 Degree Standard Observer"].align(wvl)
-
-		opt = optimize.minimize(
-			fun, ccp0, (Lab, wvlp, cmfs, ill_sd, ill_XYZ), jac=jac,
-			method="L-BFGS-B", options={"disp": verbose}
-		)
-		if verbose:
-			print(opt)
-		return opt
-
-	# A special case that's hard to solve numerically
-	if np.allclose(target_XYZ, [0, 0, 0]):
-		return np.array([0, 0, -1e+9]), 0 # FIXME: dtype?
-
-	if verbose:
-		print("Trying the target directly, XYZ=%s" % target_XYZ)
-	opt = do_optimize(target_XYZ, ccp0)
-	if opt.fun < 0.1 or not try_hard:
-		return opt.x, opt.fun
-
-	good_XYZ = (1/3, 1/3, 1/3)
-	good_ccp = (2.1276356, -1.07293026, -0.29583292) # FIXME: valid only for D65
-
-	divisions = 3
-	while divisions < 30:
-		if verbose:
-			print("Trying with %d divisions" % divisions)
-
-		keep_divisions = False
-		ref_XYZ = good_XYZ
-		ref_ccp = good_ccp
-
-		ccp0 = ref_ccp
-		for i in range(1, divisions):
-			intermediate_XYZ = ref_XYZ + (target_XYZ - ref_XYZ) * i / (divisions - 1)
-			if verbose:
-				print("Intermediate step %d/%d, XYZ=%s with ccp0=%s" %
-				      (i + 1, divisions, intermediate_XYZ, ccp0))
-			opt = do_optimize(intermediate_XYZ, ccp0)
-			if opt.fun > 1e-3:
-				if verbose:
-					print("WARNING: intermediate optimization failed")
-				break
-			else:
-				good_XYZ = intermediate_XYZ
-				good_ccp = opt.x
-				keep_divisions = True
-
-			ccp0 = opt.x
-		else:
-			return opt.x, opt.fun
-
-		if not keep_divisions:
-			divisions += 2
-
-	raise Exception("optimization failed for target_XYZ=%s, ccp0=%s" \
-	                % (target_XYZ, ccp0))
-
 # Makes a comparison plot with SDs and swatches
 def plot_comparison(target, matched_sd, label, error, ill_sd, show=True):
 	if type(target) is SpectralDistribution: