;+
; NAME:
;	lsqLinearFit
; PURPOSE:
;	Least-squares fit of function of N independent variables
; CATEGORY:
;	gen/idl/toolbox/math
; CALLING SEQUENCE:
	FUNCTION lsqLinearFit, f,x,y,z,t,	$
		yzero	= yzero 	, $
		unit	= unit		, $
		torigin = torigin	, $
		weights = weights
; INPUTS:
;	f			array[n]; type: float
;					function values
;					f can be a time structure (see keywords unit and torigin)
; OPTIONAL PARAMETERS:
;	x			array[n] or array[m,n]; type: float; default: findgen(n)
;					array[n  ]: n 1-dimensional position vectors
;					array[m,n]; n m-dimensional position vectors
;
;					If x is the only argument specified (i.e. y,z,t
;					are NOT present) then x can be a time structure
;					(see keywords unit and torigin)
;
;	y,z,t		array[n] (same as x); type: float
;					n position vectors [x,y,z,t].
;
;	weights=weights
;				array[n] (same as x); type: float
;					weights for each data points
;					(usually proportional to the inverse of the variance)
;	/yzero		fits to a linear function with no leading constant
;				(i.e. it fits to a linear function that goes through
;				the origin x=y=z=t=0)
;				The intercept Result[0] is explicitly set to zero
;				in this case.
;
;	torigin=torigin
;				array[1]; type: time structure; default: none
;	unit=unit	scalar; type: integer; default: TimeUnit(/day)
;					if f and/or x are time structures then the lsq fit
;					is done using the specified time units relative
;					to time origin 'torigin'.
;					If torigin is not defined then it is defined
;					internally (as the earliest time specified in f
;					and x, truncated to time units of 'unit')
; OUTPUTS:
;	Result		array[m+1]; type: float
;					least-square fit parameters of the linear
;					function fitted to the points.
; OPTIONAL PARAMETERS:
;	torigin=torigin
;				array[1]; type: time structure
;					If f and/or x are time structures and torigin
;					is not defined as input then the internally
;					defined value is returned here.
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; SEE ALSO:
;	normalfit, circlefit
; CALLS:
;	SuperArray, InitVar, IsType, IsTime, TimeLimits, TimeGet
; PROCEDURE:
;	For 1,2,3 and 4-dim fits the solution is evaluated directly
;	from analytic expressions. For higher dimensions a matrix
;	inversion is used.
; MODIFICATION HISTORY:
;	FEB-2003, Paul Hick (UCSD/CASS)
;	SEP-2006, Paul Hick (UCSD/CASS)
;		Added /yzero.
;	DEC-2006, Paul Hick (UCSD/CASS)
;		Added keywords unit and torigin
;	MAY-2007, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Added weights keyword.
;-

InitVar, yzero, /key
InitVar, unit , TimeUnit(/day)

; If weights are specified, normalize them

have_weights = IsType(weights,/defined)
IF have_weights THEN ww = weights/total(weights)

n = n_params()
IF n EQ 0 THEN message, 'syntax: r = linearlsqfit(f [,x,y,z,t, /yzero])'

nf = n_elements(f)
ff = f

ff_IsTime = IsTime(ff)
need_origin = NOT IsTime(torigin)

IF ff_IsTime THEN		$
	IF need_origin THEN	$
		torigin = TimeGet( TimeLimits(ff,/min), unit, /botime )

CASE n OF

1: xx = findgen(nf)			; 1-dim lsq; only fnc values specified; use array index for x

2: BEGIN					; 1-dim lsq

	x0_IsTime = IsTime(x)

	CASE x0_IsTime OF
	0: xx = x
	1: BEGIN
		IF need_origin THEN BEGIN
			xx = TimeGet( TimeLimits(x,/min), unit, /botime )
			IF IsType(torigin,/defined) THEN torigin = TimeLimits( [torigin,xx], /min ) ELSE torigin = xx
		ENDIF
		xx = TimeGet(/subtract, x, torigin, unit)
	END
	ENDCASE

END

3: BEGIN					; 2-dim lsq

	xx = SuperArray(x,n-1,/lead)

	xx[0,*] = x
	xx[1,*] = y

END

4: BEGIN					; 3-dim lsq

	xx = SuperArray(x,n-1,/lead)

	xx[0,*] = x
	xx[1,*] = y
	xx[2,*] = z

END

5: BEGIN					; 4-dim lsq

	xx = SuperArray(x,n-1,/lead)

	xx[0,*] = x
	xx[1,*] = y
	xx[2,*] = z
	xx[3,*] = t

END

ELSE: ; Nothing to do		; 5-dim lsq and higher

ENDCASE

IF ff_IsTime THEN ff = TimeGet(/subtract, ff, torigin, unit)

nx = n_elements(xx)
nn = nx/nf
IF nn*nf NE nx THEN message, 'incompatible f and x arrays'

ff = reform(ff,nf,/overwrite)
xx = reform(xx,nn,nx/nn,/overwrite)

CASE yzero OF

0: BEGIN

	CASE nn OF

	1: BEGIN				; 1-dim lsq

		tx = mean(xx)
		tf = mean(ff)

		dx = xx-tx
		df = ff-tf

		xscale = max(abs(dx))
		fscale = max(abs(df))

		dx = dx/xscale
		df = df/fscale

		CASE have_weights OF
		0: BEGIN
			sxx = mean(dx*dx)
			sxf = mean(dx*df)
		END
		1: BEGIN
			sxx = total(ww*dx*dx)
			sxf = total(ww*dx*df)
		END
		ENDCASE

		aa = sxf/sxx*(fscale/xscale)

		rr = [tf-aa*tx,aa]

	END

	2: BEGIN				; 2-dim lsq

		tx = mean(xx[0,*])
		ty = mean(xx[1,*])
		tf = mean(ff)

		dx = xx[0,*]-tx
		dy = xx[1,*]-ty
		df = ff 	-tf

		xscale = max(abs(dx))
		yscale = max(abs(dy))
		fscale = max(abs(df))

		dx = dx/xscale
		dy = dy/yscale
		df = df/fscale

		CASE have_weights OF
		0: BEGIN
			sxx = mean(dx*dx)
			sxy = mean(dx*dy)
			syy = mean(dy*dy)
			sxf = mean(dx*df)
			syf = mean(dy*df)
		END
		1: BEGIN
			sxx = total(ww*dx*dx)
			sxy = total(ww*dx*dy)
			syy = total(ww*dy*dy)
			sxf = total(ww*dx*df)
			syf = total(ww*dy*df)
		END
		ENDCASE

		aa = sxf*syy-syf*sxy
		bb = syf*sxx-sxf*sxy

		qq = sxx*syy-sxy*sxy

		aa = aa/qq*(fscale/xscale)
		bb = bb/qq*(fscale/yscale)

		rr = [tf-(aa*tx+bb*ty), aa, bb]

	END

	3: BEGIN			; 3-dim lsq

		tx = mean(xx[0,*])
		ty = mean(xx[1,*])
		tz = mean(xx[2,*])
		tf = mean(ff)

		dx = xx[0,*]-tx
		dy = xx[1,*]-ty
		dz = xx[2,*]-tz
		df = ff 	-tf

		xscale = max(abs(dx))
		yscale = max(abs(dy))
		zscale = max(abs(dz))
		fscale = max(abs(df))

		dx = dx/xscale
		dy = dy/yscale
		dz = dz/zscale
		df = df/fscale

		CASE have_weights OF
		0: BEGIN
			sxx = mean(dx*dx)
			sxy = mean(dx*dy)
			sxz = mean(dx*dz)
			syy = mean(dy*dy)
			syz = mean(dy*dz)
			szz = mean(dz*dz)
			sxf = mean(dx*df)
			syf = mean(dy*df)
			szf = mean(dz*df)
		END
		1: BEGIN
			sxx = total(ww*dx*dx)
			sxy = total(ww*dx*dy)
			sxz = total(ww*dx*dz)
			syy = total(ww*dy*dy)
			syz = total(ww*dy*dz)
			szz = total(ww*dz*dz)
			sxf = total(ww*dx*df)
			syf = total(ww*dy*df)
			szf = total(ww*dz*df)
		END
		ENDCASE

		aa = (syy*szz-syz*syz)*sxf+(sxz*syz-sxy*szz)*syf+(sxy*syz-sxz*syy)*szf
		bb = (syz*sxz-sxy*szz)*sxf+(szz*sxx-sxz*sxz)*syf+(sxy*sxz-syz*sxx)*szf
		cc = (syz*sxy-sxz*syy)*sxf+(sxz*sxy-syz*sxx)*syf+(sxx*syy-sxy*sxy)*szf

		qq = (syy*szz-syz*syz)*sxx+(sxz*syz-sxy*szz)*sxy+(sxy*syz-sxz*syy)*sxz

		aa = aa/qq*(fscale/xscale)
		bb = bb/qq*(fscale/yscale)
		cc = cc/qq*(fscale/zscale)

		rr = [tf-(aa*tx+bb*ty+cc*tz), aa, bb, cc]

	END

	4: BEGIN				; 4-dim lsq

		tx = mean(xx[0,*])
		ty = mean(xx[1,*])
		tz = mean(xx[2,*])
		tt = mean(xx[3,*])
		tf = mean(ff)

		dx = xx[0,*]-tx
		dy = xx[1,*]-ty
		dz = xx[2,*]-tz
		dt = xx[3,*]-tt
		df = ff 	-tf

		xscale = max(abs(dx))
		yscale = max(abs(dy))
		zscale = max(abs(dz))
		tscale = max(abs(dt))
		fscale = max(abs(df))

		dx = dx/xscale
		dy = dy/yscale
		dz = dz/zscale
		dt = dt/tscale
		df = df/fscale

		CASE have_weights OF
		0: BEGIN
			sxx = mean(dx*dx)
			sxy = mean(dx*dy)
			sxz = mean(dx*dz)
			sxt = mean(dx*dt)
			syy = mean(dy*dy)
			syz = mean(dy*dz)
			syt = mean(dy*dt)
			szz = mean(dz*dz)
			szt = mean(dz*dt)
			stt = mean(dt*dt)
			sxf = mean(dx*df)
			syf = mean(dy*df)
			szf = mean(dz*df)
			stf = mean(dt*df)
		END
		1: BEGIN
			sxx = total(ww*dx*dx)
			sxy = total(ww*dx*dy)
			sxz = total(ww*dx*dz)
			sxt = total(ww*dx*dt)
			syy = total(ww*dy*dy)
			syz = total(ww*dy*dz)
			syt = total(ww*dy*dt)
			szz = total(ww*dz*dz)
			szt = total(ww*dz*dt)
			stt = total(ww*dt*dt)
			sxf = total(ww*dx*df)
			syf = total(ww*dy*df)
			szf = total(ww*dz*df)
			stf = total(ww*dt*df)
		END
		ENDCASE

		aa = (  syy*szz*stt-szt*szt*syy-syt*syt*szz-syz*syz*stt+2*syz*szt*syt	 )*sxf+ $
			 ( (szt*szt-szz*stt)*sxy+(syz*stt-szt*syt)*sxz+(szz*syt-szt*syz)*sxt )*syf+ $
			 ( (stt*syz-szt*syt)*syz+(syt*syt-syy*stt)*sxz+(szt*syy-syz*syt)*sxt )*szf+ $
			 ( (syt*szz-syz*szt)*sxy+(syy*szt-syt*syz)*sxz+(syz*syz-szz*syy)*sxt )*stf

		bb = ( (sxz*stt-sxt*szt)*syz+(szz*sxt-sxz*szt)*syt+(szt*szt-stt*szz)*sxy )*sxf+ $
			 (  sxx*szz*stt-sxt*sxt*szz-sxz*sxz*stt-szt*szt*sxx+2*szt*sxt*sxz	 )*syf+ $
			 ( (sxt*sxt-stt*sxx)*syz+(sxx*szt-sxt*sxz)*syt+(stt*sxz-sxt*szt)*sxy )*szf+ $
			 ( (sxx*szt-sxz*sxt)*syz+(sxz*sxz-sxx*szz)*syt+(sxt*szz-szt*sxz)*sxy )*stf

		cc = ( (syy*sxt-syt*sxy)*szt+(syt*syt-syy*stt)*sxz+(sxy*stt-sxt*syt)*syz )*sxf+ $
			 ( (syt*sxx-sxy*sxt)*szt+(stt*sxy-syt*sxt)*sxz+(sxt*sxt-sxx*stt)*syz )*syf+ $
			 (  sxx*syy*stt-sxy*sxy*stt-syt*syt*sxx-sxt*sxt*syy+2*sxt*sxy*syt	 )*szf+ $
			 ( (sxy*sxy-sxx*syy)*szt+(syy*sxt-sxy*syt)*sxz+(sxx*syt-sxy*sxt)*syz )*stf

		dd = ( (syz*syz-syy*szz)*sxt+(szz*sxy-syz*sxz)*syt+(syy*sxz-syz*sxy)*szt )*sxf+ $
			 ( (szz*sxy-sxz*syz)*sxt+(sxz*sxz-sxx*szz)*syt+(syz*sxx-sxy*sxz)*szt )*syf+ $
			 ( (sxz*syy-syz*sxy)*sxt+(sxx*syz-sxz*sxy)*syt+(sxy*sxy-sxx*syy)*szt )*szf+ $
			 (  sxx*syy*szz-syz*syz*sxx-sxz*sxz*syy-sxy*sxy*szz+2*sxy*syz*sxz	 )*stf

		qq = (  syy*szz*stt-szt*szt*syy-syt*syt*szz-syz*syz*stt+2*syz*szt*syt	 )*sxx+ $
			 ( (szt*szt-szz*stt)*sxy+(syz*stt-szt*syt)*sxz+(szz*syt-szt*syz)*sxt )*sxy+ $
			 ( (stt*syz-szt*syt)*syz+(syt*syt-syy*stt)*sxz+(szt*syy-syz*syt)*sxt )*sxz+ $
			 ( (syt*szz-syz*szt)*sxy+(syy*szt-syt*syz)*sxz+(syz*syz-szz*syy)*sxt )*sxt

		aa = aa/qq*(fscale/xscale)
		bb = bb/qq*(fscale/yscale)
		cc = cc/qq*(fscale/zscale)
		dd = dd/qq*(fscale/tscale)

		rr = [tf-(aa*tx+bb*ty+cc*tz+dd*tt), aa, bb, cc, dd]

	END

	ELSE: BEGIN 			; 5-dim lsq and higher

		one = replicate(1,nf)

		tx = total(xx,2)/nf
		tf = mean(ff)

		dx = xx-tx#one
		df = ff-tf

		xscale = replicate(dx[0],nn)
		FOR i=0,nn-1 DO xscale[i] = max(abs(dx[i,*]))
		fscale = max(abs(df))

		dx = dx/(xscale#one)
		df = df/fscale

		CASE have_weights OF
		0: sxx = (dx#transpose(dx))/nf
		1: BEGIN
			ww = replicate(1,nn)#ww
			sxx = ((ww*dx)#transpose(dx))
		END
		ENDCASE

		sxx = LA_INVERT(sxx, /double, status=status) 
		IF status NE 0 THEN message, 'singular matrix?'

		CASE have_weights OF
		0: aa = sxx#((dx#df)/nf)
		1: aa = sxx#((ww*dx)#df)
		ENDCASE

		aa = aa*(fscale/xscale)

		rr = [tf-total(tx*aa),aa]

	END

	ENDCASE

END

1: BEGIN

	CASE nn OF

	1: BEGIN				; 1-dim lsq

		CASE have_weights OF
		0: BEGIN
			sxx = mean(xx*xx)
			sxf = mean(xx*ff)
		END
		1: BEGIN
			sxx = total(ww*xx*xx)
			sxf = total(ww*xx*ff)
		END
		ENDCASE

		rr = [0,sxf/sxx]

	END

	2: BEGIN				; 2-dim lsq

		dx = xx[0,*]
		dy = xx[1,*]
		df = ff

		CASE have_weights OF
		0: BEGIN
			sxx = mean(dx*dx)
			sxy = mean(dx*dy)
			syy = mean(dy*dy)
			sxf = mean(dx*df)
			syf = mean(dy*df)
		END
		1: BEGIN
			sxx = total(ww*dx*dx)
			sxy = total(ww*dx*dy)
			syy = total(ww*dy*dy)
			sxf = total(ww*dx*df)
			syf = total(ww*dy*df)
		END
		ENDCASE

		rr = [0,[sxf*syy-syf*sxy, syf*sxx-sxf*sxy]/(sxx*syy-sxy*sxy)]

	END

	3: BEGIN				; 3-dim lsq

		dx = xx[0,*]
		dy = xx[1,*]
		dz = xx[2,*]
		df = ff

		CASE have_weights OF
		0: BEGIN
			sxx = mean(dx*dx)
			sxy = mean(dx*dy)
			sxz = mean(dx*dz)
			syy = mean(dy*dy)
			syz = mean(dy*dz)
			szz = mean(dz*dz)
			sxf = mean(dx*df)
			syf = mean(dy*df)
			szf = mean(dz*df)
		END
		1: BEGIN
			sxx = total(ww*dx*dx)
			sxy = total(ww*dx*dy)
			sxz = total(ww*dx*dz)
			syy = total(ww*dy*dy)
			syz = total(ww*dy*dz)
			szz = total(ww*dz*dz)
			sxf = total(ww*dx*df)
			syf = total(ww*dy*df)
			szf = total(ww*dz*df)
		END
		ENDCASE

		aa = (syy*szz-syz*syz)*sxf+(sxz*syz-sxy*szz)*syf+(sxy*syz-sxz*syy)*szf
		bb = (syz*sxz-sxy*szz)*sxf+(szz*sxx-sxz*sxz)*syf+(sxy*sxz-syz*sxx)*szf
		cc = (syz*sxy-sxz*syy)*sxf+(sxz*sxy-syz*sxx)*syf+(sxx*syy-sxy*sxy)*szf

		qq = (syy*szz-syz*syz)*sxx+(sxz*syz-sxy*szz)*sxy+(sxy*syz-sxz*syy)*sxz

		rr = [0,[aa, bb, cc]/qq]

	END

	4: BEGIN				; 4-dim lsq

		dx = xx[0,*]
		dy = xx[1,*]
		dz = xx[2,*]
		dt = xx[3,*]
		df = ff

		CASE have_weights OF
		0: BEGIN
			sxx = mean(dx*dx)
			sxy = mean(dx*dy)
			sxz = mean(dx*dz)
			sxt = mean(dx*dt)
			syy = mean(dy*dy)
			syz = mean(dy*dz)
			syt = mean(dy*dt)
			szz = mean(dz*dz)
			szt = mean(dz*dt)
			stt = mean(dt*dt)
			sxf = mean(dx*df)
			syf = mean(dy*df)
			szf = mean(dz*df)
			stf = mean(dt*df)
		END
		1: BEGIN
			sxx = total(ww*dx*dx)
			sxy = total(ww*dx*dy)
			sxz = total(ww*dx*dz)
			sxt = total(ww*dx*dt)
			syy = total(ww*dy*dy)
			syz = total(ww*dy*dz)
			syt = total(ww*dy*dt)
			szz = total(ww*dz*dz)
			szt = total(ww*dz*dt)
			stt = total(ww*dt*dt)
			sxf = total(ww*dx*df)
			syf = total(ww*dy*df)
			szf = total(ww*dz*df)
			stf = total(ww*dt*df)
		END
		ENDCASE

		aa = (  syy*szz*stt-szt*szt*syy-syt*syt*szz-syz*syz*stt+2*syz*szt*syt	 )*sxf+ $
			 ( (szt*szt-szz*stt)*sxy+(syz*stt-szt*syt)*sxz+(szz*syt-szt*syz)*sxt )*syf+ $
			 ( (stt*syz-szt*syt)*syz+(syt*syt-syy*stt)*sxz+(szt*syy-syz*syt)*sxt )*szf+ $
			 ( (syt*szz-syz*szt)*sxy+(syy*szt-syt*syz)*sxz+(syz*syz-szz*syy)*sxt )*stf

		bb = ( (sxz*stt-sxt*szt)*syz+(szz*sxt-sxz*szt)*syt+(szt*szt-stt*szz)*sxy )*sxf+ $
			 (  sxx*szz*stt-sxt*sxt*szz-sxz*sxz*stt-szt*szt*sxx+2*szt*sxt*sxz	 )*syf+ $
			 ( (sxt*sxt-stt*sxx)*syz+(sxx*szt-sxt*sxz)*syt+(stt*sxz-sxt*szt)*sxy )*szf+ $
			 ( (sxx*szt-sxz*sxt)*syz+(sxz*sxz-sxx*szz)*syt+(sxt*szz-szt*sxz)*sxy )*stf

		cc = ( (syy*sxt-syt*sxy)*szt+(syt*syt-syy*stt)*sxz+(sxy*stt-sxt*syt)*syz )*sxf+ $
			 ( (syt*sxx-sxy*sxt)*szt+(stt*sxy-syt*sxt)*sxz+(sxt*sxt-sxx*stt)*syz )*syf+ $
			 (  sxx*syy*stt-sxy*sxy*stt-syt*syt*sxx-sxt*sxt*syy+2*sxt*sxy*syt	 )*szf+ $
			 ( (sxy*sxy-sxx*syy)*szt+(syy*sxt-sxy*syt)*sxz+(sxx*syt-sxy*sxt)*syz )*stf

		dd = ( (syz*syz-syy*szz)*sxt+(szz*sxy-syz*sxz)*syt+(syy*sxz-syz*sxy)*szt )*sxf+ $
			 ( (szz*sxy-sxz*syz)*sxt+(sxz*sxz-sxx*szz)*syt+(syz*sxx-sxy*sxz)*szt )*syf+ $
			 ( (sxz*syy-syz*sxy)*sxt+(sxx*syz-sxz*sxy)*syt+(sxy*sxy-sxx*syy)*szt )*szf+ $
			 (  sxx*syy*szz-syz*syz*sxx-sxz*sxz*syy-sxy*sxy*szz+2*sxy*syz*sxz	 )*stf

		qq = (  syy*szz*stt-szt*szt*syy-syt*syt*szz-syz*syz*stt+2*syz*szt*syt	 )*sxx+ $
			 ( (szt*szt-szz*stt)*sxy+(syz*stt-szt*syt)*sxz+(szz*syt-szt*syz)*sxt )*sxy+ $
			 ( (stt*syz-szt*syt)*syz+(syt*syt-syy*stt)*sxz+(szt*syy-syz*syt)*sxt )*sxz+ $
			 ( (syt*szz-syz*szt)*sxy+(syy*szt-syt*syz)*sxz+(syz*syz-szz*syy)*sxt )*sxt

		rr = [0,[aa, bb, cc, dd]/qq]

	END

	ELSE: BEGIN 				; 5-dim lsq and higher

		CASE have_weights OF
		0: sxx = (xx#transpose(xx))/nf
		1: BEGIN
			ww = replicate(1,nn)#ww
			sxx = ((ww*xx)#transpose(xx))
		END
		ENDCASE

		sxx = LA_INVERT(sxx, /double, status=status) 
		IF status NE 0 THEN message, 'singular matrix?'

		CASE have_weights OF
		0: rr = [0,sxx#((xx#ff)/nf)]
		1: rr = [0,sxx#((ww*xx)#ff)]
		ENDCASE

	END

	ENDCASE

END

ENDCASE

RETURN, rr  &  END