Recipes_Dist.cpp
上传用户:szb0815
上传日期:2007-06-13
资源大小:338k
文件大小:7k
- /*
- Software and source code Copyright (C) 1998-2000 Stanford University
- Written by Michael Eisen (eisen@genome.stanford.edu)
- 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 (eisen@genome.stanford.edu).
- 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