docs: def decompose(
  			cls,
  			heights,
  			t            = None,
  			t0           = None,
  			interval     = None,
  			constituents = constituent.noaa,
  			initial      = None,
  			n_period     = 2,
  			callback     = None,
  			full_output  = False
  		):
  		"""
  		Return an instance of Tide which has been fitted to a series of tidal observations.
  		Arguments:
  		It is not necessary to provide t0 or interval if t is provided.
  		heights -- ndarray of tidal observation heights
  		t -- ndarray of tidal observation times
  		t0 -- datetime representing the time at which heights[0] was recorded
  		interval -- hourly interval between readings
  		constituents -- list of constituents to use in the fit (default: constituent.noaa)
  		initial -- optional Tide instance to use as first guess for least squares solver
  		n_period -- only include constituents which complete at least this many periods (default: 2)
  		callback -- optional function to be called at each iteration of the solver
  		full_output -- whether to return the output of scipy's leastsq solver (default: False)
  		"""
  		if t is not None:
  			if isinstance(t[0], datetime):
  				hours = Tide._hours(t[0], t)
  				t0 = t[0]
  			elif t0 is not None:
  				hours = t
  			else:
  				raise ValueError("t can be an array of datetimes, or an array "
  				                 "of hours since t0 in which case t0 must be "
  				                 "specified.")
  		elif None not in [t0, interval]:
  			hours = np.arange(len(heights)) * interval
  		else:
  			raise ValueError("Must provide t(datetimes), or t(hours) and "
  			                 "t0(datetime), or interval(hours) and t0(datetime) "
  			                 "so that each height can be identified with an "
  			                 "instant in time.")
  
  		#Remove duplicate constituents (those which travel at exactly the same
  		#speed, irrespective of phase)
  		constituents = list(OrderedDict.fromkeys(constituents))
  
  		#No need for least squares to find the mean water level constituent z0,
  		#work relative to mean
  		constituents = [c for c in constituents if not c == constituent._Z0]
  		z0 = np.mean(heights)
  		heights = heights - z0
  
  		#Only analyse frequencies which complete at least n_period cycles over
  		#the data period.
  		constituents = [
  			c for c in constituents
  			if 360.0 * n_period < hours[-1] * c.speed(astro(t0))
  		]
  		n = len(constituents)
  
  		sort = np.argsort(hours)
  		hours = hours[sort]
  		heights = heights[sort]
  
  		#We partition our time/height data into intervals over which we consider
  		#the values of u and f to assume a constant value (that is, their true
  		#value at the midpoint of the interval).  Constituent
  		#speeds change much more slowly than the node factors, so we will
  		#consider these constant and equal to their speed at t0, regardless of
  		#the length of the time series.
  
  		partition = 240.0
  
  		t     = Tide._partition(hours, partition)
  		times = Tide._times(t0, [(i + 0.5)*partition for i in range(len(t))])
  
  		speed, u, f, V0 = Tide._prepare(constituents, t0, times, radians = True)
  
  		#Residual to be minimised by variation of parameters (amplitudes, phases)
  		def residual(hp):
  			H, p = hp[:n, np.newaxis], hp[n:, np.newaxis]
  			s = np.concatenate([
  				Tide._tidal_series(t_i, H, p, speed, u_i, f_i, V0)
  				for t_i, u_i, f_i in izip(t, u, f)
  			])
  			res = heights - s
  			if callback:
  				callback(res)
  			return res
  
  		#Analytic Jacobian of the residual - this makes solving significantly
  		#faster than just using gradient approximation, especially with many
  		#measurements / constituents.
  		def D_residual(hp):
  			H, p = hp[:n, np.newaxis], hp[n:, np.newaxis]
  			ds_dH = np.concatenate([
  				f_i*np.cos(speed*t_i+u_i+V0-p)
  				for t_i, u_i, f_i in izip(t, u, f)],
  				axis = 1)
  
  			ds_dp = np.concatenate([
  				H*f_i*np.sin(speed*t_i+u_i+V0-p)
  				for t_i, u_i, f_i in izip(t, u, f)],
  				axis = 1)
  
  			return np.append(-ds_dH, -ds_dp, axis=0)
  
  		#Initial guess for solver, haven't done any analysis on this since the
  		#solver seems to converge well regardless of the initial guess We do
  		#however scale the initial amplitude guess with some measure of the
  		#variation
  		amplitudes = np.ones(n) * (np.sqrt(np.dot(heights, heights)) / len(heights))
  		phases     = np.ones(n)
  
  		if initial:
  			for (c0, amplitude, phase) in initial.model:
  				for i, c in enumerate(constituents):
  					if c0 == c:
  						amplitudes[i] = amplitude
  						phases[i] = d2r*phase
  
  		initial = np.append(amplitudes, phases)
  
  		lsq = leastsq(residual, initial, Dfun=D_residual, col_deriv=True, ftol=1e-7)
  
  		model = np.zeros(1+n, dtype=cls.dtype)
  		model[0] = (constituent._Z0, z0, 0)
  		model[1:]['constituent'] = constituents[:]
  		model[1:]['amplitude'] = lsq[0][:n]
  		model[1:]['phase'] = lsq[0][n:]
  
  		if full_output:
  			return cls(model = model, radians = True), lsq
  		return cls(model = model, radians = True)