recipes_dist.cpp

著名的基因芯片表达数据聚类软件。esein编写

/*

Software and source code Copyright (C) 1998-2000 Stanford University

Written by Michael Eisen ([e-mail])

This software is copyright under the following conditions:

Permission to use, copy, and modify this software and its documentation
is hereby granted to all academic and not-for-profit institutions
without fee, provided that the above copyright notice and this permission
notice appear in all copies of the software and related documentation.
Permission to distribute the software or modified or extended versions
thereof on a not-for-profit basis is explicitly granted, under the above
conditions. However, the right to use this software in conjunction with
for profit activities, and the right to distribute the software or modified or
extended versions thereof for profit are *NOT* granted except by prior
arrangement and written consent of the copyright holders.

Use of this source code constitutes an agreement not to criticize, in any
way, the code-writing style of the author, including any statements regarding
the extent of documentation and comments present.

The software is provided "AS-IS" and without warranty of ank kind, express,
implied or otherwise, including without limitation, any warranty of
merchantability or fitness for a particular purpose.

In no event shall Stanford University or the authors be liable for any special,
incudental, indirect or consequential damages of any kind, or any damages
whatsoever resulting from loss of use, data or profits, whether or not
advised of the possibility of damage, and on any theory of liability,
arising out of or in connection with the use or performance of this software.

This code was written using Borland C++ Builder 4 (Inprise Inc., www.inprise.com)
and may be subject to certain additional restrictions as a result.

*/

//---------------------------------------------------------------------------
#pragma hdrstop

#include "Recipes_Dist.h"

float *circx,*circy,*circw;
int circn;

/*
Becuase of restrictions in the Numerical Recipes License, we are unable
to distribute NR source code used in this program. For information on how to
obtain this code, see www.nr.com, or contact the author ([e-mail]).
You should obtain and insert the standard NR c code for the following routines
here:

medfit
rofunc
select
frprmn
linmin
brent
mnbrak
f1dim
kstwo
probks
sort
pearsn
betai
gammln
betacf

I am in the process of writing these routines out of the program. Some are, in
fact, already unnecessary.

*/


/*

The following are extensions/modifications of NR code. I think I am allowed
to release these

*/

float fcircfit(float p[])
{
        int i;
        float val = 0;
        float val1,val2;

        for (i=1;i<=circn;i++)
        {
                val1 = sqrt(pow(circx[i]-p[1],2)+pow(circy[i]-p[2],2));
                val2 = val1 - p[3];
                val += circw[i] * pow(val2,2);
        }

        return val;
}

void fcircderiv(float p[], float xi[])
{
        int i;
        float val1,val2;

        xi[1] = 0;
        xi[2] = 0;
        xi[3] = 0;

        for (i=1;i<=circn;i++)
        {
                val1 = sqrt(pow(circx[i]-p[1],2)+pow(circy[i]-p[2],2));
                val2 = val1 - p[3];

                if (val1 > 0)
                {

                xi[1] += circw[i] * 2*p[1]*val2*(circx[i] - p[1])/val1;
                xi[2] += circw[i] * 2*p[2]*val2*(circy[i] - p[2])/val1;
                xi[3] += circw[i] * 2 * val2;
                }
        }
}


#define NRANSI
#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;
#define M 7
#define NSTACK 50

void sort2fi(unsigned long n, float arr[], int brr[])
{
	unsigned long i,ir=n,j,k,l=1;
	int *istack,jstack=0;
	float a,temp;
        int b;

	istack=ivector(1,NSTACK);
	for (;;) {
		if (ir-l < M) {
			for (j=l+1;j<=ir;j++) {
				a=arr[j];
				b=brr[j];
				for (i=j-1;i>=l;i--) {
					if (arr[i] <= a) break;
					arr[i+1]=arr[i];
					brr[i+1]=brr[i];
				}
				arr[i+1]=a;
				brr[i+1]=b;
			}
			if (!jstack) {
				free_ivector(istack,1,NSTACK);
				return;
			}
			ir=istack[jstack];
			l=istack[jstack-1];
			jstack -= 2;
		} else {
			k=(l+ir) >> 1;
			SWAP(arr[k],arr[l+1])
			SWAP(brr[k],brr[l+1])
			if (arr[l] > arr[ir]) {
				SWAP(arr[l],arr[ir])
				SWAP(brr[l],brr[ir])
			}
			if (arr[l+1] > arr[ir]) {
				SWAP(arr[l+1],arr[ir])
				SWAP(brr[l+1],brr[ir])
			}
			if (arr[l] > arr[l+1]) {
				SWAP(arr[l],arr[l+1])
				SWAP(brr[l],brr[l+1])
			}
			i=l+1;
			j=ir;
			a=arr[l+1];
			b=brr[l+1];
			for (;;) {
				do i++; while (arr[i] < a);
				do j--; while (arr[j] > a);
				if (j < i) break;
				SWAP(arr[i],arr[j])
				SWAP(brr[i],brr[j])
			}
			arr[l+1]=arr[j];
			arr[j]=a;
			brr[l+1]=brr[j];
			brr[j]=b;
			jstack += 2;
			if (jstack > NSTACK) nrerror("NSTACK too small in sort2.");
			if (ir-i+1 >= j-l) {
				istack[jstack]=ir;
				istack[jstack-1]=i;
				ir=j-1;
			} else {
				istack[jstack]=j-1;
				istack[jstack-1]=l;
				l=i;
			}
		}
	}
}
#undef M
#undef NSTACK
#undef SWAP
#undef NRANSI

#define NRANSI
#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;
#define M 7
#define NSTACK 50

void sort3fff(unsigned long n, float arr[], float brr[], float crr[])
{
	unsigned long i,ir=n,j,k,l=1;
	int *istack,jstack=0;
	float a,temp;
        float b,c;

	istack=ivector(1,NSTACK);
	for (;;) {
		if (ir-l < M) {
			for (j=l+1;j<=ir;j++) {
				a=arr[j];
				b=brr[j];
                                c=crr[j];
				for (i=j-1;i>=l;i--) {
					if (arr[i] <= a) break;
					arr[i+1]=arr[i];
					brr[i+1]=brr[i];
                                        crr[i+1]=crr[i];
				}
				arr[i+1]=a;
				brr[i+1]=b;
                                crr[i+1]=c;
			}
			if (!jstack) {
				free_ivector(istack,1,NSTACK);
				return;
			}
			ir=istack[jstack];
			l=istack[jstack-1];
			jstack -= 2;
		} else {
			k=(l+ir) >> 1;
			SWAP(arr[k],arr[l+1])
			SWAP(brr[k],brr[l+1])
                        SWAP(crr[k],crr[l+1])
			if (arr[l] > arr[ir]) {
				SWAP(arr[l],arr[ir])
				SWAP(brr[l],brr[ir])
                                SWAP(crr[l],crr[ir])
			}
			if (arr[l+1] > arr[ir]) {
				SWAP(arr[l+1],arr[ir])
				SWAP(brr[l+1],brr[ir])
				SWAP(crr[l+1],crr[ir])
			}
			if (arr[l] > arr[l+1]) {
				SWAP(arr[l],arr[l+1])
				SWAP(brr[l],brr[l+1])
				SWAP(crr[l],crr[l+1])
			}
			i=l+1;
			j=ir;
			a=arr[l+1];
			b=brr[l+1];
                        c=crr[l+1];
			for (;;) {
				do i++; while (arr[i] < a);
				do j--; while (arr[j] > a);
				if (j < i) break;
				SWAP(arr[i],arr[j])
				SWAP(brr[i],brr[j])
				SWAP(crr[i],crr[j])
			}
			arr[l+1]=arr[j];
			arr[j]=a;
			brr[l+1]=brr[j];
			brr[j]=b;
			crr[l+1]=crr[j];
			crr[j]=c;
			jstack += 2;
			if (jstack > NSTACK) nrerror("NSTACK too small in sort2.");
			if (ir-i+1 >= j-l) {
				istack[jstack]=ir;
				istack[jstack-1]=i;
				ir=j-1;
			} else {
				istack[jstack]=j-1;
				istack[jstack-1]=l;
				l=i;
			}
		}
	}
}
#undef M
#undef NSTACK
#undef SWAP
#undef NRANSI