C+
C NAME:
C	ice_analyze
C PURPOSE:
C	Finds the optimum parameters for compressing the specified integer array with href=ice_pack=
C CATEGORY:
C	Bits and bytes
C CALLING SEQUENCE:
	function ice_analyze(kmax, kshift, klenbit, ksign, kperm, nOrig, Orig)
C INPUTS:
C	nOrig		integer		# elements in Orig
C	Orig(nOrig)	integer		array to be analyzed
C OUTPUTS:
C	ice_analyze	integer		# bits needed to hold compressed data for the optimum parameters
C
C	kmax		integer		cutoff for applying compression
C	kshift		integer		minimum # bits copied for each difference
C
C	klenbit		integer		# bits encoded  when the full value instead
C					of a difference is used
C	ksign		integer		-1: all values in Orig are <  zero
C					+1: all values in Orig are >= zero
C					 0: both positive and negative values are present
C	kperm(0:nbit-1)	integer		permutation table
C					nbit = 32 (# bits in integer*4)
C					Only elements 0...kmax of the returned array should be used.
C CALLS:
C	ArrI4GetMinMax, ArrI4Zero, IndexI4, ArrI4Reverse, iArrI4Total, ArrI4Copy
C SEE ALSO:
C	ice_pack, ice_unpack, ice_write, ice_read
C PROCEDURE:
C	kmax, kshift, klenbit, ksign, kperm are ready to be passed to ice_pack for compression
C	of Orig. The compression factor is ibit/(nOrig*nbit). (use nbit=16 instead of nbit=32
C	if the original data were actually integer*2 stored in integer*4 form).
C	See href=ice_pack= for more information.
C MODIFICATION HISTORY:
C	AUG-2000, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C-
	parameter	(nbit  = 32)
	parameter	(nbit1 = nbit-1)

	    integer	kmax
	    integer	kshift
	    integer	klenbit
	    integer	ksign
	    integer	kperm(0:nbit1)			! Permutation table
	    integer	nOrig
	    integer	Orig (0:nOrig-1)

	integer		nrank(0:nbit1)
	integer		nhist(0:nbit1)
	integer		nhist_shift(0:nbit1)

	!-------
	! Check properties of original data values (we check the
	! differences below)
	!
	! Find the optimum ksign. If the data are all >= 0 or < 0 then the we know that
	! the sign bit is always set or never set, so we don't need to store it.

	call ArrI4GetMinMax(nOrig,Orig,nmin,nmax)

	if (nmin .ge. 0) then				! >= 0 (sign bit never set)
	    ksign =  1
	else if (nmax .lt. 0) then			! <  0 (sign bit always set)
	    ksign = -1
	else
	    ksign =  0
	end if

	!-------
	! ice_pack always stores absolute values, and adds a sign bit if needed.
	! Find out how many bits are needed to store the absolute values for a full nbit value

	nmax = max(abs(nmin),abs(nmax))

	if (nmax .eq. 0) then				! Array of zeros
	    klenbit = -1
	else
	    klenbit = nbit1				! Find highest 1-bit in absolute value
	    do while (.not. btest(nmax, klenbit))
		klenbit = klenbit-1
	    end do
	end if

	!-------
	! When storing full values (rather than differences) ice_pack needs to use klenbit bits.
	! This includes a sign bit only if ksign = 0 (see above)
	! For an array of zeros klenbit = 0 (this results in kshift = 0, kmax = 1 and a maximum
	! compression factor of nbit, since only terminator 1-bits are stored for each element).

	klenbit = klenbit+1+1-abs(ksign)


	!-------
	! Check the properties of differences of neighbouring elements.
	! First build a histogram of differences as a function of how many bits are needed
	! to store them.

	iprev = 0					! The '0th element' of Orig

	call ArrI4Zero(nbit, nhist)

	do i=0,nOrig-1					! Loop over all elements of input array

	    idel = Orig(i)-iprev			! Difference with previous element

	    if (idel .eq. 0) then			! Zero difference

		ipad = 0

	    else

		idel = abs(idel)			! Absolute difference

		itop = nbit1-1				! Find highest 1-bit in absolute difference
		do while (.not. btest(idel, itop))
		    itop = itop-1
		end do

		ipad = itop+1				

	    end if

	    !-------
	    ! ipad is the number of bits needed to store the difference: 0<= ipad <= nbit-1

	    nhist(ipad) = nhist(ipad)+1			! Accumulate histogram

	    iprev = Orig(i)

	end do

	!-------
	! Find the optimum permutation table.
	! kshift is the # of low order bits of the absolute difference that are stored unmodified.


	ktop = min(klenbit+1,nbit1)

	nbits_opt  = -1

	do kshift=0,ktop-1

	    !-------
	    ! nhist(i) (i=0,kmax) contains the number of differences that need i bits to store them.
	    !
	    ! If the first kshift bits are removed, then i-kshift bits remain. For i=0,kshift zero bits
	    ! remain. Accumulate a second histogram by summing hist(i) (i=0,kshift) into hist_shift(0),
	    ! and shifting hist(i) (i=kshift+1,nbit1) by kshift places. The permutation table is based
	    ! on this modified histogram (for kshift=0 both histograms are identical).

	    nhist_shift(0) = iArrI4Total(kshift+1, nhist, i)			! Sum 0..kshift
	    call ArrI4Copy(nbit1-kshift, nhist(kshift+1), nhist_shift(1))	! Shift kshift+1..nbit1 to  1..nbit1-kshift

	    !-------
	    ! Find the optimum permutation table (only use the first nbit-kshift elements).

	    call IndexI4(0,nbit1-kshift,0,nbit1-kshift,nhist_shift,kperm)
	    call ArrI4Reverse(nbit-kshift,kperm,kperm)

	    do i=0,nbit1-kshift					! Invert permutation table
		nrank(kperm(i)) = i
	    end do

	    !-------
	    ! Find the optimum kmax. For each kmax calculate how many bits the
	    ! compressed data will have. Pick the kmax giving the highest compression.
	    ! For a given nshift we only need consider kmax >= kshift+1.

	    do kmax=kshift+1,ktop
		ib = (nrank(0)+1+1*min(1,kshift))*nhist(0)	! 0-bit

		do i=1,kshift
		    ib = ib+(nrank(0)+1+1+1+kshift)*nhist(i)
		end do

		do i=kshift+1,kmax-1
		    ib = ib+(nrank(i-kshift)+1+i)*nhist(i)
		end do

		nbits = ib+(kmax-kshift+1+klenbit)*iArrI4Total(nbit-kmax, nhist(kmax), i)

		!print *, kshift, kmax, nbits, nOrig*16, 1.*nbits/(nOrig*16)

		if (nbits_opt .eq. -1 .or. nbits .lt. nbits_opt) then
		    kshift_opt = kshift
		    kmax_opt   = kmax
		    nbits_opt  = nbits
		end if

	    end do

	end do

	ice_analyze = nbits_opt

	kshift = kshift_opt
	kmax   = kmax_opt

	!print *, '>>', kshift, kmax, nbits_opt, nOrig*16, 1.*nbits_opt/(nOrig*16)

	!-------
	! All absolute differences requiring kmax bits
	! or more are lumped together and will be encoded as full values instead.

	nhist(kmax) = iArrI4Total(nbit-kmax, nhist(kmax), i)

	!-------
	! Update permutation table for optimum kmax. Only nhist(i), i=0,kmax is used.

	nhist_shift(0) = iArrI4Total(kshift+1, nhist, i)		! Sum 0..kshift
	call ArrI4Copy(kmax-kshift, nhist(kshift+1), nhist_shift(1))	! Shift kshift+1..kmax to 1..kmax-kshift

	call IndexI4(0,kmax-kshift, 0, kmax-kshift, nhist_shift, kperm)
	call ArrI4Reverse(kmax-kshift+1, kperm, kperm)			! Only 0..kmax-kshift should be used.

	return
	end