#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#		"End of archive 3 (of 4)."
# Contents:  error.c error.h install macro.h misc.c printimage.c
# Wrapped by cyee@dec14.cs.monash.edu.au on Wed Dec  9 15:24:12 1992
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f error.c -a "${1}" != "-c" ; then 
  echo shar: Will not over-write existing file \"error.c\"
else
echo shar: Extracting \"error.c\" \(2088 characters\)
sed "s/^X//" >error.c <<'END_OF_error.c'
X
X/*
X *  Copyright  (c)  L.Allison, C.S.Wallace, C.N.Yee.
X *     Dept. Computer Science, Monash University, Australia 3168,    May 1991.
X *     xxx@bruce.cs.monash.edu.au    where xxx = {lloyd, csw, cyee}
X *  1. These routines may only be used for non-commercial, non-classified,
X *     research purposes.
X *  2. No liability of any kind is accepted for their use.
X *  3. As a courtesy, please inform us of any interesting uses of the routines,
X *     of any bugs that you find in them and of how to correct
X *     such bugs, if known.
X *  4. Do not delete this copyright notice.
X *  5. Please include the following references in any publication based wholly
X *     or in part on the use of these routines or on the ideas contained in them:
X *   . Allison L. & C.N.Yee.
X *        Minimum Message Length Encoding and the Comparison of Macro-Molecules.
X *        Bull. Math. Bio. 52(3) 431-453 1990
X *   . Allison L., C.S.Wallace & C.N.Yee. Finite-State Models in the Alignment
X *        of Macromolecules. J. Mol. Evol. (1992) 35:77-89
X *  6. You may distribute the routines to others if and only if
X *     they agree to accept and maintain these conditions and no fee is charged.
X *  7. There is no rule 7.
X */
X/*
X *	"Error.c"
X *	The setting up of default Error handling routines table
X */
X#include <stdio.h>
X#include "error.h"
X
X/*
X *	returns the error string that corresponds the ErrCode
X */
Xchar *ErrMesgStr(ErrCode)
X{
X	switch (ErrCode) {
X
X	case ERRLineTooLong : return "String Too Long" ;
X	case ERRBadChar : return "Bad Character Read" ;
X	default :
X		printf("ERROR TRAP: illegal Error Code %d\n", ErrCode) ;
X		exit(99) ;
X	}
X	return NULL ;
X}
X
Xint ErrMesgExit(ErrCode)
X	int	ErrCode ;
X{
X	char *errmesg ;
X
X	errmesg = ErrMesgStr(ErrCode) ;
X	printf("%s\n", errmesg) ;
X		
X	exit(ErrCode) ;
X}
X
Xint ErrMesgRet(ErrCode)
X	int	ErrCode ;
X{
X	char *errmesg ;
X
X	errmesg = ErrMesgStr(ErrCode) ;
X	printf("%s\n", errmesg) ;
X}
X
XErrFileOpen(filename)
X	char filename[] ;
X{
X	printf("Error opening : %s\n",filename) ;
X	exit(0) ;
X}
X
Xint (*ErrorHandling[NError])() = {NULL, ErrMesgExit, ErrMesgRet, ErrFileOpen};
END_OF_error.c
if test 2088 -ne `wc -c <error.c`; then
    echo shar: \"error.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f error.h -a "${1}" != "-c" ; then 
  echo shar: Will not over-write existing file \"error.h\"
else
echo shar: Extracting \"error.h\" \(1419 characters\)
sed "s/^X//" >error.h <<'END_OF_error.h'
X
X/*
X *  Copyright  (c)  L.Allison, C.S.Wallace, C.N.Yee.
X *     Dept. Computer Science, Monash University, Australia 3168,    May 1991.
X *     xxx@bruce.cs.monash.edu.au    where xxx = {lloyd, csw, cyee}
X *  1. These routines may only be used for non-commercial, non-classified,
X *     research purposes.
X *  2. No liability of any kind is accepted for their use.
X *  3. As a courtesy, please inform us of any interesting uses of the routines,
X *     of any bugs that you find in them and of how to correct
X *     such bugs, if known.
X *  4. Do not delete this copyright notice.
X *  5. Please include the following references in any publication based wholly
X *     or in part on the use of these routines or on the ideas contained in them:
X *   . Allison L. & C.N.Yee.
X *        Minimum Message Length Encoding and the Comparison of Macro-Molecules.
X *        Bull. Math. Bio. 52(3) 431-453 1990
X *   . Allison L., C.S.Wallace & C.N.Yee. Finite-State Models in the Alignment
X *        of Macromolecules. J. Mol. Evol. (1992) 35:77-89
X *  6. You may distribute the routines to others if and only if
X *     they agree to accept and maintain these conditions and no fee is charged.
X *  7. There is no rule 7.
X */
X/*
X *	"Error.h"
X *	Error code definitions
X */
X#define NError		4
X/*			0 is to signify no error */
X#define ERRLineTooLong	1
X#define ERRBadChar	2
X#define ERRFileOpening	3
X
Xextern int (*ErrorHandling[NError])() ;
END_OF_error.h
if test 1419 -ne `wc -c <error.h`; then
    echo shar: \"error.h\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f install -a "${1}" != "-c" ; then 
  echo shar: Will not over-write existing file \"install\"
else
echo shar: Extracting \"install\" \(171 characters\)
sed "s/^X//" >install <<'END_OF_install'
Xcc -o rtheory 1State.c 1StateR0.c 3State.c 3StateR0.c 5State.c IO.c driver.c error.c misc.c sumlg.c -lm
Xcc -o scale scale.c misc.o -lm 
Xcc -o pimg printimage.c misc.o -lm
END_OF_install
if test 171 -ne `wc -c <install`; then
    echo shar: \"install\" unpacked with wrong size!
fi
chmod +x install
# end of overwriting check
fi
if test -f macro.h -a "${1}" != "-c" ; then 
  echo shar: Will not over-write existing file \"macro.h\"
else
echo shar: Extracting \"macro.h\" \(2339 characters\)
sed "s/^X//" >macro.h <<'END_OF_macro.h'
X
X/*
X *  Copyright  (c)  L.Allison, C.S.Wallace, C.N.Yee.
X *     Dept. Computer Science, Monash University, Australia 3168,    May 1991.
X *     xxx@bruce.cs.monash.edu.au    where xxx = {lloyd, csw, cyee}
X *  1. These routines may only be used for non-commercial, non-classified,
X *     research purposes.
X *  2. No liability of any kind is accepted for their use.
X *  3. As a courtesy, please inform us of any interesting uses of the routines,
X *     of any bugs that you find in them and of how to correct
X *     such bugs, if known.
X *  4. Do not delete this copyright notice.
X *  5. Please include the following references in any publication based wholly
X *     or in part on the use of these routines or on the ideas contained in them:
X *   . Allison L. & C.N.Yee.
X *        Minimum Message Length Encoding and the Comparison of Macro-Molecules.
X *        Bull. Math. Bio. 52(3) 431-453 1990
X *   . Allison L., C.S.Wallace & C.N.Yee. Finite-State Models in the Alignment
X *        of Macromolecules. J. Mol. Evol. (1992) 35:77-89
X *  6. You may distribute the routines to others if and only if
X *     they agree to accept and maintain these conditions and no fee is charged.
X *  7. There is no rule 7.
X */
X#include <stdio.h>
X
X/*
X *	general macro definitions
X */
X#ifdef DEBUG
X#define	debug(X)	X
X#else
X#define debug(X)
X#endif
X
X#define MAXINT	999999
X#define infinity 1e20
X
X#define LINELENGTH		75
X#define SCREENHEIGHT 	40
X#define SCREENWIDTH 	80
X
X#define bool	int
X#define	FALSE	0
X#define	TRUE	1
X
Xchar	*malloc() ;
XFILE	*popen(), *fopen() ;
Xchar	*strcat(), *strcpy() ;
Xint	strncmp() ;
X
X#define min(x,y)	((x)<(y) ? (x):(y))
X#define max(x,y)	((x)>(y) ? (x):(y))
X
X/*  make sure that temp is declared with the appropriate type when
X    this macro is used 
X */
X#define swap(x,y)	{temp = x ; x = y ; y = temp ;}  
X
X#define odd(n)	((n) % 2 ? TRUE : FALSE)
X#define even(n)	((n) % 2 ? FALSE : TRUE)
X
X#define abs(n) ((n) > 0 ? (n) : -(n))
X
Xint	min3() ;
Xfloat fmin2(), fmin3(), fmax3(), fmin5(), flabs() ;
Xlong lmin3() ;
Xdouble log2(), exp2() ;
Xdouble lloyd(), riss(), logstar(), yee() ;
Xfloat  sumlg(), sumlg3(), sumlg5() ;
Xfloat  paramcost2(), paramcost3(), paramcost4(), paramcost6(), U() ;
Xint	randomnum() ;
Xfloat NullMesgLength(), DiffCode() ;
X
Xtypedef	struct {
X	int	x, y ; 	/* used as integer */
X	} coord ;
X
Xextern double ln2 ;
X
Xfloat RCount() ;
END_OF_macro.h
if test 2339 -ne `wc -c <macro.h`; then
    echo shar: \"macro.h\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f misc.c -a "${1}" != "-c" ; then 
  echo shar: Will not over-write existing file \"misc.c\"
else
echo shar: Extracting \"misc.c\" \(10885 characters\)
sed "s/^X//" >misc.c <<'END_OF_misc.c'
X
X/*
X *  Copyright  (c)  L.Allison, C.S.Wallace, C.N.Yee.
X *     Dept. Computer Science, Monash University, Australia 3168,    May 1991.
X *     xxx@bruce.cs.monash.edu.au    where xxx = {lloyd, csw, cyee}
X *  1. These routines may only be used for non-commercial, non-classified,
X *     research purposes.
X *  2. No liability of any kind is accepted for their use.
X *  3. As a courtesy, please inform us of any interesting uses of the routines,
X *     of any bugs that you find in them and of how to correct
X *     such bugs, if known.
X *  4. Do not delete this copyright notice.
X *  5. Please include the following references in any publication based wholly
X *     or in part on the use of these routines or on the ideas contained in them:
X *   . Allison L. & C.N.Yee.
X *        Minimum Message Length Encoding and the Comparison of Macro-Molecules.
X *        Bull. Math. Bio. 52(3) 431-453 1990
X *   . Allison L., C.S.Wallace & C.N.Yee. Finite-State Models in the Alignment
X *        of Macromolecules. J. Mol. Evol. (1992) 35:77-89
X *  6. You may distribute the routines to others if and only if
X *     they agree to accept and maintain these conditions and no fee is charged.
X *  7. There is no rule 7.
X */
X/*
X *	"misc.c"  miscellaneous routines
X *
X *	SellersProb(PMatches,PChanges,PInsDel)
X *		float	*PMatches,*PChanges,*PInsDel ;
X *		returns a set of probability that is equivalent to Seller's metric 
X *
X *	changetoprobability3(n1,n2,n3,p1,p2,p3)
X *		int	n1,n2,n3 ; 	float	*p1,*p2,*p3 ;
X *		change a set of counts to a set of probability
X *
X *	changetoprobability4(n1,n2,n3,n4,p1,p2,p3,p4)
X *		int	n1,n2,n3,n4 ; 	float	*p1,*p2,*p3,*p4 ;
X *		change a set of counts to a set of probability
X *
X *	long lmin3(x, y, z)
X *		long x, y, z ;
X *
X *	int min3(x, y, z)
X *		int x, y, z ;
X *
X *	float fmin3(x, y, z), flabs(x, y)
X *		float x, y, z ;
X *
X *	double log2(x)
X *		double x ;
X *
X *	double logstar(x)
X *		int x ;
X *
X *	double riss(x)
X *		int x ;
X *
X *	double yee(x)
X *		int x ;
X *
X *	double lloyd(x)
X *		int x ;
X *		Hybrid of smoothen logstar and yee
X *
X *	double geometric(P, x)
X *		float	P ;
X *		int x ;
X *
X *	int randomnum(n)
X *	int n ;
X *		this routine returns a random number of range 0..n-1
X *		randomseed must be called before calling this routine
X *
X *	randomseed(seed)
X *	int seed ; 
X *		provide the random number with a starting seed
X */	
X#include <math.h>
X#include "Strings.h"
X#include "macro.h"
X
X/*
Xdouble	ln2, LogRissConstant, RissConstant, YeeConstant, LogYeeConstant ;
X*/
Xlong time() ;
Xdouble	ln2 = 0.69314718055994530941 ;
Xdouble	logDelta ;
Xstatic int tablesize = 0 ;
X/*
Xstatic double	risstable[MAXSEQLENGTH] ;
Xstatic double	yeetable[MAXSEQLENGTH] ;
X*/
Xstatic double	lloydtable[MAXSEQLENGTH] ;
X
Xlong lmin3(x,y,z)
X	long x,y,z ;
X{
X	long t ;
X	t = (x < y) ? x : y ;
X	return (t < z)? t : z ;
X}
X
Xint min7(a,b,c,d,e,f,g)
X	int a,b,c,d,e,f,g ;
X{
X	int	t,u,v,x,y ;
X
X	t = min(a,b) ; u = min(c,d) ; v = min(e,f) ;
X
X	x = min(t,u) ; y = min(v,g) ;
X	return min(x,y) ;
X}
X
Xint min3(x, y, z)
X	int	x, y, z ;
X{
X	int	t ;
X
X	t = (x < y) ? x : y ;
X	return (t < z)? t : z ;
X}
X
Xfloat flabs(a)
X	float a ;
X{
X	if (a < 0) return -a ;
X	else return a ;
X}
X
Xfloat fmin2(x, y)
X	float	x, y ;
X{
X	return x < y ? x : y;
X}
X
Xfloat fmin3(x, y, z)
X	float	x, y, z ;
X{
X	float	t ;
X
X	t = (x < y) ? x : y ;
X	return (t < z)? t : z ;
X}
X
Xfloat fmin5(a, b, c, d, e)
X	float	a, b, c, d, e ;
X{ return fmin3(min(a,b),min(c,d),e) ; }
X
Xfloat fmax3(x, y, z)
X	float	x, y, z ;
X{
X	float	t ;
X
X	t = (x > y) ? x : y ;
X	return (t > z)? t : z ;
X}
X
X/*
X *	generates a random number from 0..n-1
X */
Xint randomnum(n)
X	int n ;
X{
X	int t ;
X
X	t = rand()/1000 ; /* get rid of bottom 3 digits */
X	return t % n ;
X}
X
Xrandomseed(seed)
X	int	seed ;
X{ srand(seed) ; }
X
X/*
X *	use system clock to initialize the random number generator
X */
Xgeneraterandomseed(seed)
X	int *seed ;
X{
X	long clock ;
X	clock = time((long *) 0) ;
X	*seed = (int) clock ;
X	randomseed(*seed) ;
X}
X
Xdouble	exp2(x)
X	double	x ;
X{ return exp(x*ln2) ; }
X
Xdouble log2(x)
X	double x ;
X{ return (double) (log(x) / ln2) ; }
X
Xdouble geometric(P, x) 
X	float	P ;
X	int	x ;
X{
X	return (x-1)*log2((double)P) + log2((double)1-P) ;
X}
X
Xdouble logstar(x)
X	int	x ;
X{
X	double	y, u ;
X
X	u = (double) x ;
X	y = 0 ;
X	while (u > 1.0) {
X		u = log2(u) ;
X		y += u ;
X	}
X	return y ;
X}
X
X/*
Xstatic double ris(x)
X	int	x ;
X{ 
X	return logstar(x) + LogRissConstant ;
X}
X
Xdouble riss(x)
X	int	x ;
X{
X	if ( x < tablesize) return risstable[x] ;
X	else return ris(x) ;
X}
X*/
X
X
X/* log(x) + log^2(x+1) + log^3(x+3) + log^4(x+15) + ... */
Xstatic double yeelog(x) 
X	int x ; /* x>=1 */
X{
X	int 	i, j ;
X	double	z, y ;
X	static double	power[] = {1, 2, 4, 16, 65536 /*64*1024*/} ;
X
X	y = 0;
X	for (i = 0 ; i < 5; i++) {
X		z = (double) (x + power[i] - 1) ;
X	    for (j = 0 ; j <= i; j++) z = log2(z);
X	    y += z ;
X	}
X	return y ;
X}
X
X/*
Xdouble yees(x)
X	int x ;
X{
X	return yeelog(x) + LogYeeConstant ;
X}
X
Xdouble yee(x)
X	int	x ;
X{
X	if ( x < tablesize) return yeetable[x] ;
X	else return yees(x) ;
X}
X*/
X
Xdouble lloyds(x)
X	int	x ;
X{
X   double	ans ;
X
X	if (x <= 8)  {
X		if (x == 2 || x == 4)
X			ans = 0.45*logstar(x-1)+0.55*logstar(x+1) ; /* concave */
X	 	else ans = logstar(x);
X	 	ans += (yeelog(9)-logstar(9)) ; /* smooth join */
X	} else /* x>8 */
X         ans = yeelog(x);
X
X    return ans + logDelta ;
X}
X
Xdouble lloyd(x)
X	int	x ;
X{
X	if ( x < tablesize) return lloydtable[x] ;
X	else return lloyds(x) ;
X}
X
Xloginit(size)
X	int	size ;
X{
X	int	i ;
X
X	if (size > MAXSEQLENGTH) {
X		printf("loginit: too long\n") ;
X		exit(99) ;
X	}
X
X/*
X	RissConstant = 2.865064 ;
X	LogRissConstant = log2(RissConstant) ;
X	YeeConstant = 2.37564 ;
X	LogYeeConstant = log2(YeeConstant) ;
X*/
X
X   	logDelta = log2((double)2.279);
X	tablesize = size ;
X
X/*
X	risstable[0] = 0 ;
X	for (i=1 ; i <= size; i++)
X		risstable[i] = ris(i) ;
X
X	yeetable[0] = 0 ;
X	for (i=1 ; i <= size; i++)
X		yeetable[i] = yees(i) ;
X*/
X
X	lloydtable[0] = 0 ;
X	for (i=1 ; i <= size; i++)
X		lloydtable[i] = lloyds(i) ;
X}
X
X/*--------------------------------------------------------------------------*/
X
X/*
X *	this routine return a set of probability that for R-model
X *	will give the same alignment as Seller's metric
X */
XSellersProb(PMatches,PChanges,PInsDel)
X	float *PMatches,*PChanges,*PInsDel ;
X{
X	*PMatches = 0.5743492 ;
X	*PChanges = 0.1542580 ;
X	*PInsDel = 1 - *PMatches - *PChanges ;
X}
X
X#define SENSITIVITY 0.5
X
Xchangetoprobability3(n1,n2,n3,p1,p2,p3)
X	int	n1,n2,n3 ;
X	float	*p1,*p2,*p3 ;
X{
X	float sum ;
X	sum = n1 + n2 + n3 + 3*SENSITIVITY ;
X	*p1 = (n1+SENSITIVITY)/sum ;
X	*p2 = (n2+SENSITIVITY)/sum ;
X	*p3 = (n3+SENSITIVITY)/sum ;
X}
X
Xfchangetoprobability2(n1,n2,p1,p2)
X	float	n1,n2,*p1,*p2 ;
X{ 
X	*p1 = (n1+SENSITIVITY)/(n1+n2+2*SENSITIVITY) ;
X	*p2 = (n2+SENSITIVITY)/(n1+n2+2*SENSITIVITY) ;
X}
X
Xfchangetoprobability3(n1,n2,n3,p1,p2,p3)
X	float	n1,n2,n3,*p1,*p2,*p3 ;
X{
X	float sum ;
X	sum = n1 + n2 + n3 + 3*SENSITIVITY ;
X	*p1 = (n1+SENSITIVITY)/sum ;
X	*p2 = (n2+SENSITIVITY)/sum ;
X	*p3 = (n3+SENSITIVITY)/sum ;
X}
X
Xfchangetoprobability4(n1,n2,n3,n4,p1,p2,p3,p4)
X	float	n1,n2,n3,n4 ;
X	float	*p1,*p2,*p3,*p4 ;
X{
X	float sum ;
X	sum = n1+n2+n3+n4+4*SENSITIVITY ;
X	*p1 = (n1+SENSITIVITY)/sum ; *p2 = (n2+SENSITIVITY)/sum ;
X	*p3 = (n3+SENSITIVITY)/sum ; *p4 = (n4+SENSITIVITY)/sum ;
X}
X
Xfchangetoprobability6(n1,n2,n3,n4,n5,n6,p1,p2,p3,p4,p5,p6)
X	float	n1,n2,n3,n4,n5,n6 ;
X	float	*p1,*p2,*p3,*p4,*p5,*p6 ;
X{
X	float sum ;
X	sum = n1+n2+n3+n4+n5+n6+6*SENSITIVITY ;
X	*p1 = (n1+SENSITIVITY)/sum ; *p2 = (n2+SENSITIVITY)/sum ;
X	*p3 = (n3+SENSITIVITY)/sum ; *p4 = (n4+SENSITIVITY)/sum ;
X	*p5 = (n5+SENSITIVITY)/sum ; *p6 = (n6+SENSITIVITY)/sum ;
X}
X
X
X/*---------------------------------------------------------------------*/
Xfloat lgstar(x)
X	float x ;
X{
X	if (x > 1) {
X		x = (float)log2((double)x) ;
X		return x + lgstar(x) ;
X	} else return 0 ;
X}
X
X/*
X *	increase in mes' len' due to using an est' of prob's rather than true
X *	prob's which are of course unknown
X */
Xstatic float errorcost(K)
X	int K ;
X{ return (float)((K-1.0)/2.0 * log2(exp((double)1))) ;}
X
X
X/*
X *	estimate cost of transmitting probability of params in a header.
X *	The parameter has 6 values, 5 degrees of freedom.
X */
Xfloat paramcost6(fm/*atch*/, fc/*hange*/, fh/*horiz*/, fv/*ert*/,fa,fb)
X	float	fm, fc, fh, fv, fa, fb ;
X{
X	float V, deltaVsq, fsum, pm, pc, ph, pv, pa, pb ;
X	int K ;
X
X	K = 6 ;
X	fsum = fm+fc+fh+fv+fa+fb ;
X	fsum += 0.05 ; /* in case fsum is 0 */
X	fchangetoprobability6(fm,fc,fh,fv,fa,fb,&pm,&pc,&ph,&pv,&pa,&pb) ;
X	deltaVsq =  (pm*pc*ph*pv*pa*pb) / exp((double)(K-1)*log((double)fsum/12));
X	V = 1.0/(5*4*3*2*1) ; /*=1/(K-1)!*/
X
X	return  (float)(log((double)1 + V*V/deltaVsq)/2.0/log((double)2) + errorcost(K)) ;
X}
X
X/*
X *	estimate cost of transmitting probability of params in a header.
X *	The parameter has 4 values, 3 degrees of freedom.
X */
Xfloat paramcost4(/*frequencies*/ fm/*atch*/, fc/*hange*/, fh/*horiz*/, fv/*ert*/)
X	float	fm, fc, fh, fv ;
X{
X	float V, deltaVsq, fsum, pm, pc, ph, pv ;
X	int	K ;
X
X	K = 4 ;
X	fsum = fm+fc+fh+fv ;
X	fsum += 0.05 ; /* in case fsum is 0 */
X	fchangetoprobability4(fm,fc,fh,fv,&pm,&pc,&ph,&pv) ;
X	deltaVsq =  (pm*pc*ph*pv) / exp((double)(K-1)*log((double)fsum/12));
X	V = 1.0/(3.0*2.0) ; /*=1/(4-1)!*/
X	return  (float)(log((double)1 + V*V/deltaVsq)/2.0/log((double)2)+errorcost(K)) ;
X}
X
Xfloat paramcost2(fm, fc)
X	float fm, fc ;
X{
X	float V, deltaVsq, fsum, pm, pc ;
X	int K ;
X
X	K = 2 ;
X	fsum = fm+fc;
X	fsum += 0.05 ; /* in case fsum == 0 */
X	fchangetoprobability2(fm,fc,&pm,&pc) ;
X	deltaVsq = (float)(pm*pc/exp((double)(K-1)*log((double)fsum/12)));
X	V = 1.0/1.0 ; /* = 1/(K-1)! */
X	return (float)(log((double)1 + V*V/deltaVsq)/2.0/log((double)2) + errorcost(K)) ;
X}
X
X
Xfloat paramcost3(fm, fc, fid)
X	float fm, fc, fid ;
X{
X	float V, deltaVsq, fsum, pm, pc, pid ;
X	int K ;
X
X	K = 3 ;
X	fsum = fm+fc+fid;
X	fsum += 0.05 ; /* in case fsum == 0 */
X	fchangetoprobability3(fm,fc,fid,&pm,&pc,&pid) ;
X	deltaVsq = (float)(pm*pc*pid/exp((double)(K-1)*log((double)fsum/12)));
X	V = 1.0/2.0 ; /* = 1/(K-1)! */
X	return (float)(log((double)1 + V*V/deltaVsq)/2.0/log((double)2)+errorcost(K)) ;
X}
X
X/* 
X *	Code for integers ~ universal prior
X */
Xfloat U( L )
X	float L ; /* L>=1 */
X{ return (float) (lgstar(L) + log2((double)2.865)) ;}
X
X
X/*
X *	Binomial Code for len1 - len2
X *
X *       1          1                 n=0
X *     1   1        2                 n=1
X *   1   2   1      4                 n=2
X * 1   3   3   1  = 8 = 2**3          n=3
X */
Xfloat DiffCode(len1, len2 )
X	int	len1, len2 ;
X{
X	int i, n, r, sum, diff ;
X	float	codeword ;
X
X	sum = len1+len2; diff = len1-len2 /*can be + or -*/;
X	/*NB. sum is even <=> diff is even*/
X	n = sum; r = (sum+diff) / 2 /* >= 0 and <= n */;
X	codeword = n;  /*p(r) = C(n,r)/2**n*/
X	for (i = r+1; i <= n; i++ ) codeword = codeword - log2((double)i);
X	for (i = 1 ; i <= n-r; i++) codeword = codeword + log2((double)i);
X	return codeword ;
X}
X
X/*
X *	Compute the null theory message length for two strings of length m and n
X */
Xfloat NullMesgLength(m,n)
X{
X	return (float)U((float)m+n)+DiffCode(m,n)+(m+n)*2;
X}
END_OF_misc.c
if test 10885 -ne `wc -c <misc.c`; then
    echo shar: \"misc.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f printimage.c -a "${1}" != "-c" ; then 
  echo shar: Will not over-write existing file \"printimage.c\"
else
echo shar: Extracting \"printimage.c\" \(6982 characters\)
sed "s/^X//" >printimage.c <<'END_OF_printimage.c'
X
X/*
X *  Copyright  (c)  L.Allison, C.S.Wallace, C.N.Yee.
X *     Dept. Computer Science, Monash University, Australia 3168,    May 1991.
X *     xxx@bruce.cs.monash.edu.au    where xxx = {lloyd, csw, cyee}
X *  1. These routines may only be used for non-commercial, non-classified,
X *     research purposes.
X *  2. No liability of any kind is accepted for their use.
X *  3. As a courtesy, please inform us of any interesting uses of the routines,
X *     of any bugs that you find in them and of how to correct
X *     such bugs, if known.
X *  4. Do not delete this copyright notice.
X *  5. Please include the following references in any publication based wholly
X *     or in part on the use of these routines or on the ideas contained in them:
X *   . Allison L. & C.N.Yee.
X *        Minimum Message Length Encoding and the Comparison of Macro-Molecules.
X *        Bull. Math. Bio. 52(3) 431-453 1990
X *   . Allison L., C.S.Wallace & C.N.Yee. Finite-State Models in the Alignment
X *        of Macromolecules. J. Mol. Evol. (1992) 35:77-89
X *  6. You may distribute the routines to others if and only if
X *     they agree to accept and maintain these conditions and no fee is charged.
X *  7. There is no rule 7.
X */
X/*
X *	"Printimage.c"
X *	transform a greyscale image into poscript for lazer writer
X *
X *	input :
X *		height width   of image 
X *		Astring
X *		Bstring
X *		[height-1,0], ... , [height-1,width-1]
X *				.
X *				.
X *				.
X *		[0,0],        ... , [0, width-1]
X *		Caption - optional,  if present will be printed at bottom of page
X */
X#include <stdio.h>
X#include "macro.h"
X
X#define ImageSize 450	/* ~6.5 inches,  72 pixels/inch */
X#define LegendSpace 216 /* 3 inches,     ""          */
X#define Xmargin 80
X#define Ymargin 80
X#define LegendDot 25
X
Xchar *FontType = "Courier-Bold" ;
X
X#define MAXSEQLENGTH 2000
X
Xint width, height ;
Xint usersize ;
Xfloat PixelsPerDot ;
Xint x, y, Ymax, fontsize ;
Xchar A[MAXSEQLENGTH], B[MAXSEQLENGTH] ;
Xchar Label[100] ;
X
Xint i, j ;
X
Xreadln()
X{
X	while (getchar() != '\n') ;
X}
X
Xreadstring(str)
X	char *str ;
X{
X	int i ;
X
X	while ((*str = getchar()) != '\n' && *str != EOF) str++ ;
X	*str = '\0' ;
X}
X
X
X/*
X *	determines the pointsize to use for labelling the density plot 
X */
Xint FontSize(size)
X	float size ;
X{ return (int) size+1 ; }
X
Xprinthexdigit(n)
X	int n ;
X{
X	if (n <= 9) putchar((char)(n+'0')) ;
X	else putchar((char)(n-10+'a')) ;
X}
X
X/*
X *	change a decimal number to hex ... works only for two digits number!!
X */
Xprinthex(value)
X	int value ;
X{
X	printhexdigit(value/16) ;
X	printhexdigit(value%16) ;
X}
X
Xsetfont(fonttype, fontsize)
X	char *fonttype ;
X	int fontsize ;
X{
X	printf("/%s findfont\n",fonttype) ;
X	printf("%d scalefont setfont\n", fontsize) ;
X}
X
X/*
X *	left justify a string at x,y
X */
XLeftJustify(s, x, y) 
X	char *s ;
X	int x, y ;
X{
X	printf("%d ",x) ;
X	printf("(%s) stringwidth\n",s);
X	printf("pop sub\n") ;
X	printf("%d moveto\n",y) ; /* y coordinate */
X	printf("(%s) show\n", s);
X}
X
X/*
X *	generates postscript command to center a string at x, y in a space of width
X */
Xcenterstring(s, x, y, width)
X	char *s ;
X	int x, y, width ;
X{
X	printf("%d\n",width) ;
X	printf("(%s) stringwidth\n",s) ;
X	printf("pop sub 2 div %d add\n",x) ;
X	printf("%d moveto\n", y) ;
X	printf("(%s) show\n",s) ;
X}
X
XCaption()
X{
X	char string[100] ;
X
X	readstring(string) ;
X	setfont("Times-Roman",14) ;
X	centerstring(string,Xmargin, 60, ImageSize) ;
X	readstring(string) ;
X	centerstring(string,Xmargin, 40, ImageSize) ;
X}
X
Xdolabel()
X{
X	setfont("Helvetica-Bold",14) ;
X	readstring(Label) ;
X	centerstring(Label,Xmargin, Ymargin+LegendSpace-40, ImageSize) ;
X}
X
XStringAt(s, x, y)
X	char *s ;
X	int x, y ;
X{
X	printf("%d %d moveto\n",x, y) ;
X	printf("(%s) show\n", s) ;
X}
X
XDoLegend(greyvalue, x, y, label)
X{
X	printf("%d %d translate\n",x, y) ;
X	printf("%d %d scale\n",LegendDot, LegendDot) ;
X	printf("1 1 8 [1 0 0 1 0 1]\n") ;
X	printf("{currentfile dot readhexstring pop} image\n") ;
X	printhex(greyvalue) ;
X	putchar('\n') ;
X	printf("grestore\n") ;
X	setfont("Helvetica", 12) ;
X	StringAt(label, x+LegendDot+15, y-LegendDot/2) ;
X}
X
XLegend()
X{
X	int x, y, greyvalue ;
X	int label[100] ;
X
X	printf("/dot 1 string def\n") ;
X	x = Xmargin + 20 ;
X	y = Ymargin + LegendSpace - 50 ;
X	while (scanf("%d",&greyvalue) == 1) {
X		readstring(label) ;
X		DoLegend(greyvalue,x, y, label) ;
X		y = y - LegendDot - 5 ;
X		if (y < Ymargin) break ;
X	}
X	readln() ;
X}
X
Xskipblanks()
X{
X	int c ;
X	while ((c = getchar()) == ' ') ;
X	ungetc(c,stdin) ;
X}
X
Xgetfield(s)
X	char *s ;
X{
X	skipblanks() ;
X	while ((*s = getchar()) != ' ' && *s != '\n') s++ ;
X	*s = '\0' ;
X}
X/*
X *	int the form of horizontal keys
X */
X#define keyheight 14
X#define labelwidth 	50
X#define labelgap	15
XLegend2()
X{
X	char label[100], keystr[10];
X	int nkeys, i, keylength, keyval ;
X
X	readstring(label) ;
X	setfont("Helvetica-Oblique",11) ;
X	x = Xmargin + labelwidth ;
X	y = Ymargin + LegendSpace - 50 ;
X	LeftJustify("Key", x, y) ;
X	LeftJustify(label, x, y - 15) ;
X
X	keylength = ImageSize -  labelwidth -labelgap ;
X
X	printf("/keystr 256 string def\n") ;
X	printf("%d %d translate\n", x+labelgap, y-4) ;
X	printf("%d %d scale\n", keylength, keyheight) ;
X	printf("256 1 8 ") ;
X	printf("[256 0 0 -1 0 1]\n") ;
X	printf("{currentfile keystr readhexstring pop}\n") ;
X	printf("image\n") ;
X	for ( i = 0 ; i < 256; i++) {
X		printhex(i) ;
X	}
X	putchar('\n') ;
X
X	printf("\ngrestore\n") ;
X	y = y - 15 ;
X	x = x + labelgap ;
X	setfont("Helvetica",9) ;
X	scanf("%d",&nkeys) ;
X	for ( i = 0; i <= nkeys; i++) {
X		getfield(keystr) ;
X		StringAt(keystr,x,y) ;
X		x = x+(int)((float)keylength/nkeys+0.5) ;
X	}
X	readln() ;
X}
X
Xmain()
X{
X	int greyvalue ;
X
X	scanf("%d %d",&height, &width) ;
X	usersize = width > height ? width : height ;
X	readln() ;
X	readstring(A) ;
X	readstring(B) ;
X	PixelsPerDot = (float)ImageSize/usersize ;
X	printf("initgraphics\n") ;
X	printf("/picstr %d string def\n", width) ;
X	printf("%d %d translate\n", Xmargin, (int)(Ymargin+LegendSpace-(usersize-height)*PixelsPerDot)) ;
X	printf("%d %d scale\n", ImageSize, ImageSize) ;
X	printf("%d %d 8 ",width, height) ;
X	printf("[%d 0 0 %d 0 %d]\n",usersize,-usersize,usersize) ;
X	printf("{currentfile picstr readhexstring pop}\n") ;
X	printf("image\n") ;
X	for ( i = 0 ; i < height; i++) {
X		for ( j = 0; j < width; j++) {
X			scanf("%d", &greyvalue) ;
X			printhex(greyvalue) ;
X		}
X		putchar('\n') ;
X		readln() ;
X	}
X	printf("\ngrestore\n") ;
X
X	/*
X	 *	Now print the strings alongside the density plot
X	 */
X	fontsize = FontSize(PixelsPerDot) ;
X	setfont(FontType, fontsize) ;
X
X	y = (int)(Ymargin + LegendSpace + PixelsPerDot*height+fontsize/2.0) ;
X	for ( i = 0; i < width-1; i++) {
X		x = (int)(Xmargin + (i+1)*PixelsPerDot-fontsize/4.0+0.5);
X		printf("%d %d moveto",x,y) ;
X		printf(" (%c) show\n",B[i]) ;
X	}
X	Ymax = (int)(Ymargin + LegendSpace + PixelsPerDot*height) ;
X
X	x = (int)(Xmargin - fontsize);
X	for ( i = 0 ; i < height-1; i++) {
X		y = (int)(Ymargin + LegendSpace + PixelsPerDot*(height-i-1) - fontsize/3.0+0.5) ; 
X		printf("%d %d moveto",x,y) ;
X		printf(" (%c) show\n",A[i]) ;
X	}
X	Legend2() ;
X	Caption() ;
X	printf("showpage\n") ;
X}
END_OF_printimage.c
if test 6982 -ne `wc -c <printimage.c`; then
    echo shar: \"printimage.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
echo shar: End of archive 3 \(of 4\).
cp /dev/null ark3isdone
MISSING=""
for I in 1 2 3 4 ; do
    if test ! -f ark${I}isdone ; then
	MISSING="${MISSING} ${I}"
    fi
done
if test "${MISSING}" = "" ; then
    echo You have unpacked all 4 archives.
    rm -f ark[1-9]isdone
else
    echo You still need to unpack the following archives:
    echo "        " ${MISSING}
fi
##  End of shell archive.
exit 0