source: trunk/GDE/SUPPORT/lsadt.c

Last change on this file was 19480, checked in by westram, 17 months ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 27.3 KB
Line 
1/*
2  PROGRAM LSADT
3  Fortran->C conversion by Mike Maciukenas, CPGA, Microbiology at University
4  of Illinois.
5  C-----------------------------------------------------------------------
6  C
7  C       LEAST SQUARES ALGORITHM FOR FITTING ADDITIVE TREES TO
8  C       PROXIMITY DATA
9  C
10  C       GEERT DE SOETE  --  VERSION 1.01 - FEB. 1983
11  C                           VERSION 1.02 - JUNE 1983
12  C                           VERSION 1.03 - JULY 1983
13  C
14  C       REFERENCE: DE SOETE, G. A LEAST SQUARES ALGORITHM FOR FITTING
15  C          ADDITIVE TREES TO PROXIMITY DATA. PSYCHOMETRIKA, 1983, 48,
16  C          621-626.
17  C          DE SOETE, G. ADDITIVE TREE REPRESENTATIONS OF INCOMPLETE
18  C          DISSIMILARITY DATA. QUALITY AND QUANTITY, 1984, 18,
19  C          387-393.
20  C       REMARKS
21  C       -------
22  C       2) UNIFORMLY DISTRIBUTED RANDOM NUMBERS ARE GENERATED BY A
23  C          PROCEDURE DUE TO SCHRAGE. CF.
24  C          SCHRAGE, L.  A MORE PORTABLE FORTRAN RANDOM NUMBER GENERATOR.
25  C          ACM TRANS. ON MATH. SOFTW., 1979, 5, 132-138.
26  C       3) SUBROUTINES VA14AD AND VA14AC (translated into minfungra) ARE
27  C          ADAPTED FROM THE HARWELL SUBROUTINE LIBRARY (1979 EDITION).
28  C       4) ALTHOUGH THIS PROGRAM HAS BEEN CAREFULLY TESTED, THE
29  C          AUTHOR DISCLAIMS ANY RESPONSABILITY FOR POSSIBLE
30  C          ERRORS.
31  C
32  C-----------------------------------------------------------------------
33*/
34#include <stdio.h>
35#include <stdlib.h>
36#include <string.h>
37#include <unistd.h>
38#include <math.h>
39
40#ifdef DARWIN
41#include <float.h>
42#define MAXDOUBLE DBL_MAX
43#define MINDOUBLE DBL_MIN
44#else
45#include <values.h>
46#endif
47
48#include <ctype.h>
49#include <fcntl.h>
50#include <errno.h>
51
52
53#define BUFLEN  1024
54#define MAXLEAVES       256
55
56static int m, n, dissim, pr, start, save, seed, nempty;
57static double ps1, ps2, GLOBAL_f, empty, tol;
58static char fname[1000];
59static char *names[MAXLEAVES];
60static double *delta[MAXLEAVES];
61static double **d;
62static double **g;
63static double **dold;
64static FILE *reportf;
65static int report;
66extern int errno;
67double nfac;
68
69extern double strtod();
70
71double dabs(a)
72     double a;
73{
74        return((a<0.0) ? -a : a);
75}
76
77double sqr(a)
78     double a;
79{
80        return(a*a);
81}
82
83double max(a, b)
84     double a;
85     double b;
86{
87        return((a>b)?a:b);
88}
89
90int imin(a, b)
91     int a;
92     int b;
93{
94        return((a<b)?a:b);
95}
96
97double min(a, b)
98     double a;
99     double b;
100{
101        return((a<b)?a:b);
102}
103
104void show_error(message)
105     char *message;
106{
107        printf("\n>>>>ERROR:\n>>>>%s\n", message);
108        exit(0);
109}
110
111
112float runif()
113{
114#define A 16807
115#define P 2147483647
116#define B15 32768
117#define B16 65536
118
119        int xhi, xalo, leftlo, fhi, k;
120        float t1, t2;
121
122        xhi = seed/B16;
123        xalo = (seed - (xhi * B16)) * A;
124        leftlo = xalo / B16;
125        fhi = xhi*A+leftlo;
126        k = fhi/B15;
127        seed = (((xalo - leftlo*B16) - P) + (fhi - k*B15)*B16) + k;
128        if(seed<0) seed += P;
129        t1 = seed;
130        t2 = t1 * 4.656612875e-10;
131        return (t2); /* seed * (1/(2**31-1))*/
132}
133
134float rnorm()
135{
136        static float t1, t2;
137        static float n1, n2;
138        static int second = 0;
139
140        if(second)
141        {
142                second = 0;
143                n1 = sin(t2);
144                n2 = t1*n1;
145                n1 = t1*sin(t2);
146                return(n2);
147        }
148        else
149        {
150                t1 = runif();
151                t1 = sqrt(-2.0*log(t1));
152                n2 = 6.283185308;
153                t2 = runif()*n2;
154                n1 = cos(t2);
155                n2 = t1*n1;
156                n1 = t1*cos(t2);
157                second = 1;
158                return(n2);
159        }
160}
161#define gd(i,j)  ( (j<i)? d[i][j] : d[j][i] )
162
163#if 0
164double gd(i, j)
165     int i, j;
166{
167        if(j<i)
168                return(d[i][j]);
169        else if(j>i)
170                return(d[j][i]);
171        else
172                show_error("gd: i=j -- programmer screwed up!");
173}
174#endif
175
176char *repeatch(string, ch, num)
177     char *string;
178     int ch;
179     int num;
180{
181        for(string[num--] = '\0'; num >= 0; string[num--] = ch);
182        return(string);
183}
184
185int getachar()
186     /* skips comments! */
187{
188        static int oldchar = '\0';
189        int ch;
190        int more=1;
191
192        while(more)
193        {
194                ch = getchar();
195                if(oldchar == '\n' && ch == '#')
196                {
197                        while(ch!='\n'&&ch!=EOF)
198                                ch=getchar();
199                        oldchar = ch;
200                }
201                else if(oldchar == '\n' && isspace(ch))
202                        ;
203                else more=0;
204        }
205        oldchar = ch;
206        return(ch);
207}
208
209int skip_space()
210{
211        int ch;
212
213        while(isspace(ch=getachar()));
214        return(ch);
215}
216
217int getaword(string, len)
218     /* 0 if failed, 1 if data was read, -1 if data read to end of file */
219     char *string;
220     int len;
221{
222        int i;
223        int ch;
224        ch = skip_space();
225        if(ch == EOF)
226                return(0);
227        for(i=0; !isspace(ch) && ch != EOF; i++)
228        {
229                if(i<len-1)
230                {
231                        string[i] = ch;
232                }
233                ch = getachar();
234        }
235        string[i<len?i:len-1] = '\0';
236        if(ch==EOF)
237                return(-1);
238        else
239                return(1);
240}
241
242int readtobar(string, len)
243     /* 0 if failed, 1 if data was read */
244     char *string;
245     int len;
246{
247        int i;
248        int ch;
249
250        ch = skip_space();
251        if(ch==EOF)
252                return(0);
253        for(i=0; ch!=EOF && ch!='|'; i++)
254        {
255                if(ch=='\n'||ch=='\r'||ch=='\t')
256                        i--;
257                else
258                {
259                        if(i<len-1)
260                                string[i]=ch;
261                }
262                ch=getachar();
263        }
264        len = i<len?i:len-1;
265        for(i=len-1;i>=0 && isspace(string[i]); i--);
266        string[i+1] = '\0';
267        if(ch==EOF)
268                return(-1);
269        else
270                return(1);
271}
272
273int readtobarorcolon(string, len)
274     /* 0 if failed, 1 if data was read */
275     char *string;
276     int len;
277{
278        int i;
279        int ch;
280
281        ch = skip_space();
282        if(ch==EOF)
283                return(0);
284        for(i=0; ch!=EOF && ch!='|' && ch!=':'; i++)
285        {
286                if(ch=='\n'||ch=='\r'||ch=='\t')
287                        i--;
288                else
289                {
290                        if(i<len-1)
291                                string[i]=ch;
292                }
293                ch=getachar();
294        }
295        len = i<len?i:len-1;
296        for(i=len-1;i>=0 && isspace(string[i]); i--);
297        string[i+1] = '\0';
298        if(ch==EOF)
299                return(-1);
300        else
301                return(1);
302}
303
304char *getmem(nelem, elsize)
305     unsigned nelem, elsize;
306{
307        char *temp;
308
309        temp = (char *)calloc(nelem+1, elsize);
310        if(temp == NULL) show_error("Couldn't allocate memory.");
311        else return(temp);
312}
313
314void show_help()
315{
316        printf("\nlsadt--options:\n");
317        printf("  -f file  - write report to 'file'\n");
318        printf("  -d       - treat data as dissimilarities (default)\n");
319        printf("  -s       - treat data as similarities\n");
320        printf("  -print 0 - don't print out iteration history (default)\n");
321        printf("  -print n>0 - print out iteration history every n iterations\n");
322        printf("  -print n<0 - print out iteration history every n iterations\n");
323        printf("               (including current distance estimates & gradients)\n");
324        printf("  -init  n - initial parameter estimates (default = 3)\n");
325        printf("         n=1 - uniformly distributed random numbers\n");
326        printf("         n=2 - error-perturbed data\n");
327        printf("         n=3 - original distance data from input matrix\n");
328        printf("  -save    - save final parameter estimates (default is don't save\n");
329        printf("  -seed  n - seed for random number generator (default = 12345)\n");
330        printf("  -ps1   n - convergence criterion for inner iterations (default = 0.0001)\n");
331        printf("  -ps2   n - convergence criterion for major iterations (default = 0.0001)\n");
332        printf("  -empty n - missing data indicator (default = 0.0)\n");
333        printf("  -help    - show this help text\n");
334        exit(0);
335}
336
337
338int get_parms(argc, argv)
339     int argc;
340     char **argv;
341{
342        int i;
343        int cur_arg;
344
345        /* codes for current argument:
346        ** 0 = no current argument
347        ** 1 = pr
348        ** 2 = start
349        ** 3 = seed
350        ** 4 = ps1
351        ** 5 = ps2
352        ** 6 = empty
353        ** 7 = filename */
354
355        dissim = 0;
356        pr = 0;
357        start = 2;
358        save = 0;
359        seed = 12345;
360        ps1 = 0.0001;
361        ps2 = 0.0001;
362        empty = 0.0;
363        n = 0;
364        cur_arg = 0;
365        for(i=1; i<argc; i++)
366                switch(cur_arg) {
367            case 0:
368                if(strcmp(argv[i],"-d")==0)
369                    dissim = 0;
370                else if(strcmp(argv[i],"-s")==0)
371                    dissim = 1;
372                else if(strcmp(argv[i],"-print")==0)
373                    cur_arg = 1;
374                else if(strcmp(argv[i],"-init")==0)
375                    cur_arg = 2;
376                else if(strcmp(argv[i],"-save")==0)
377                    save = 1;
378                else if(strcmp(argv[i],"-seed")==0)
379                    cur_arg = 3;
380                else if(strcmp(argv[i],"-ps1")==0)
381                    cur_arg = 4;
382                else if(strcmp(argv[i],"-ps2")==0)
383                    cur_arg = 5;
384                else if(strcmp(argv[i],"-empty")==0)
385                    cur_arg = 6;
386                else if(strcmp(argv[i],"-f")==0)
387                    cur_arg = 7;
388                else if(strcmp(argv[i],"-help")==0)
389                    show_help();
390                else
391                {
392                    printf("\nlsadt: bad option\n");
393                    show_help();
394                }
395                break;
396            case 1:
397                pr = atoi(argv[i]);
398                cur_arg = 0;
399                break;
400            case 2:
401                start = atoi(argv[i]);
402                cur_arg = 0;
403                break;
404            case 3:
405                seed = atoi(argv[i]);
406                cur_arg = 0;
407                break;
408            case 4:
409                ps1 = atof(argv[i]);
410                cur_arg = 0;
411                break;
412            case 5:
413                ps2 = atof(argv[i]);
414                cur_arg = 0;
415                break;
416            case 6:
417                empty = atof(argv[i]);
418                cur_arg = 0;
419                break;
420            case 7:
421                strcpy(fname, argv[i]);
422                cur_arg = 0;
423                break;
424            default:
425                printf("\nlsadt: bad option\n");
426                show_help();
427                break;
428                }
429
430
431    /* check validity of parameters */
432        if(ps1<0.0) ps1 = 0.0001;
433        if(ps2<0.0) ps2 = 0.0001;
434        if(dissim<0) dissim = 0;
435        if(start < 1 || start > 3) start = 3;
436        if(save != 1) save = 0;
437        if(seed < 0) seed = 12345;
438
439    /*printf("dissim=%d\n", dissim);*/
440    /*printf("pr=%d\n", pr);*/
441    /*printf("start=%d\n", start);*/
442    /*printf("save=%d\n", save);*/
443    /*printf("seed=%d\n", seed);*/
444    /*printf("ps1=%f\n", ps1);*/
445    /*printf("ps2=%f\n", ps2);*/
446    /*printf("empty=%f\n", empty);*/
447}
448
449int get_data()
450{
451        int i, j, more;
452        char buf[BUFLEN];
453        char *ptr;
454        char ch;
455        int result;
456        double temp, nfactor, datmin, datmax;
457
458
459        nempty = n = 0;
460        more = 1;
461        ptr = &ch;
462        while(more)
463        {
464                result=readtobarorcolon(buf, BUFLEN);
465                if(result == 0 || result == -1)
466                        more = 0;
467                else
468                {
469                        n++;
470                        names[n] = getmem(BUFLEN, 1);
471                        result=readtobar(buf, BUFLEN);
472                        if(result != 1)
473                                show_error("get_data: bad name syntax, or missing '|'");
474                        strcpy(names[n], buf);
475                        if(n>1)
476                                delta[n]=(double *)getmem(n, sizeof(double));
477                        else
478                                delta[n]=NULL;
479                        for(j=1; j<n; j++)
480                        {
481                                if(!getaword(buf, BUFLEN))
482                                        show_error("get_data: missing data values.\n");
483                                delta[n][j] = strtod(buf, &ptr);
484                                if(ptr == buf)
485                                        show_error("get_data: invalid data value.\n");
486                                if(delta[n][j] == empty) nempty++;
487                        }
488                }
489        }
490        if(n<4) show_error("Must be at least 4 nodes.");
491        m = n * (n - 1) / 2;
492        /* make all positive */
493        datmin = MAXDOUBLE;
494        for(i=2; i<=n; i++)
495                for(j=1; j<i; j++)
496                        if(delta[i][j] < datmin && delta[i][j] != empty)
497                                datmin = delta[i][j];
498        if(datmin < 0.0)
499                for(i=2; i<=n; i++)
500                        for(j=1; j<i; j++)
501                                if(delta[i][j] != empty)
502                                        delta[i][j] -= datmin;
503        /* convert dissimilarities into similarities if -s option specified */
504        if(dissim)
505        {
506                datmax = MINDOUBLE;
507                for(i=2; i<=n; i++)
508                        for(j=1; j<i; j++)
509                                if(delta[i][j] > datmax && delta[i][j] != empty)
510                                        datmax = delta[i][j];
511                datmax += 1.0;
512                for(i=2; i<=n; i++)
513                        for(j=1; j<i; j++)
514                                if(delta[i][j] != empty)
515                                        delta[i][j] = datmax - delta[i][j];
516        }
517
518
519        /* make sum of squared deviations from delta to average(delta) = 1 */
520        temp = 0.0;
521        for(i=2; i<=n; i++)
522                for(j=1; j<i; j++)
523                        if(delta[i][j] != empty)
524                                temp += delta[i][j];
525        temp = temp/(double)(m-nempty);
526        /* now temp = average of all non-empty cells */
527        nfactor = 0.0;
528        for(i=2; i<=n; i++)
529                for(j=1; j<i; j++)
530                        if(delta[i][j] != empty)
531                                nfactor+=sqr(delta[i][j]-temp);
532        nfactor = 1.0/sqrt(nfactor);
533        nfac = nfactor;
534        /* now nfactor is the normalization factor */
535        if(report)
536                fprintf(reportf, "Normalization factor: %15.10f\n", nfactor);
537        for(i=2; i<=n; i++)
538                for(j=1; j<i; j++)
539                        if(delta[i][j] != empty)
540                                delta[i][j] *= nfactor;
541        /* now all is normalized */
542}
543
544int initd()
545{
546
547        int i, j, k;
548        double t1, t2;
549        int emp, nemp;
550        double dmax;
551        double efac, estd;
552
553
554        efac = 0.333333333333333;
555        if(start == 1)
556        {
557                for(i=2; i<=n; i++)
558                        for(j=1; j<i; j++)
559                                d[i][j] = runif();
560                return 0; 
561        }
562        if(nempty == 0 && start == 2)
563        {
564                estd = sqrt(efac/m);
565                for(i=2; i<=n; i++)
566                        for(j=1; j<i; j++)
567                                d[i][j] = delta[i][j] + rnorm()*estd;
568                return 0;
569        }
570        for(i=2; i<=n; i++)
571                for(j=1; j<i; j++)
572                        d[i][j] = delta[i][j];
573        if(start==3 && nempty==0) return 0;
574        emp=nempty;
575        nemp=0;
576        while(nemp<emp)
577                for(i=2; i<=n; i++)
578                        for(j=1; j<i; j++)
579                                if(d[i][j] == empty)
580                                {
581                                        dmax = MAXDOUBLE;
582                                        for(k=1; k<=n; k++)
583                                        {
584                                                t1 = gd(i,k);
585                                                t2 = gd(j,k);
586                                                if(t1!=empty && t2!=empty)
587                            dmax = min(dmax, max(t1,t2));
588                                        }
589                                        if(dmax!=MAXDOUBLE)
590                                                d[i][j] = dmax;
591                                        else
592                                                nemp++;
593                                }
594        if(nemp!=0)
595        {
596                t1 = 0.0;
597                for(i=2; i<=n; i++)
598                        for(j=1; j<i; j++)
599                                if(d[i][j]!=empty)
600                                        t1+=d[i][j];
601                t1 /= m-nemp;
602                for(i=2; i<=n; i++)
603                        for(j=1; j<i; j++)
604                                if(d[i][j]==empty)
605                                        d[i][j]=t1;
606        }
607        if(start!=3)
608        {
609                estd=sqrt(efac/(m-nempty));
610                for(i=2; i<=n; i++)
611                        for(j=1; j<i; j++)
612                                d[i][j]+= rnorm()*estd;
613        }
614        return 0;
615}
616
617static int nm0, nm1, nm2, nm3;
618static double fitl, fitp, r;
619
620void fungra()
621{
622        register double djilk,dkilj,dlikj;
623        register double fw,gki,gkj,gji;
624        double fac;
625        int i, j, k, l;
626        double wijkl;
627        for(i=2;i<=n;i++)
628                for(j=1;j<i;j++)
629                        g[i][j] = 0.0;
630        fitl = 0.0;
631        fac = 2.0;
632        for(i=2;i<=n;i++)
633                for(j=1;j<i;j++)
634                        if(delta[i][j] != empty)
635                        {
636                                fitl += sqr(d[i][j] - delta[i][j]);
637                                g[i][j] += fac*(d[i][j] - delta[i][j]);
638                        }
639        fitp = 0.0;
640        fac = 2.0*r;
641        for(i=1;i<=nm3;i++){
642                for(j=i+1;j<=nm2;j++){
643                        for(k=j+1;k<=nm1;k++){
644                                gki = gkj = gji = 0.0;
645                                for(l=k+1;l<=nm0;l++){
646                                        djilk = d[j][i]+d[l][k];
647                                        dkilj = d[k][i]+d[l][j];
648                                        dlikj = d[l][i]+d[k][j];
649
650                                        if((djilk>=dlikj)&&
651                                           (dkilj>=dlikj))
652                                        {
653                                                wijkl=djilk-dkilj;
654                                                fitp+=wijkl*wijkl;
655                                                fw = fac*wijkl;
656                                                gji+=fw;
657                                                g[l][k]+=fw;
658                                                gki-=fw;
659                                                g[l][j]-=fw;
660                                        }
661                                        else
662                        if((djilk>=dkilj)&&
663                           (dlikj>=dkilj))
664                        {
665                            wijkl=djilk-dlikj;
666                            fitp+=wijkl*wijkl;
667                            fw = fac*wijkl;
668                            gji+=fw;
669                            g[l][k]+=fw;
670                            gkj-=fw;
671                            g[l][i]-=fw;
672                        }else{
673                            wijkl=dkilj-dlikj;
674                            fitp+=wijkl*wijkl;
675                            fw = fac*wijkl;
676                            gki+=fw;
677                            g[l][j]+=fw;
678                            g[l][i]-=fw;
679                            gkj-=fw;
680                        }
681                                } /* l */
682                                g[k][i] += gki;
683                                g[k][j] += gkj;
684                                g[j][i] += gji;
685                        } /* k */
686                } /* j */
687        } /* i*/
688        GLOBAL_f = fitl+r*fitp;
689}
690
691static double **dr, **dgr, **d1, **gs, **xx, **gg;
692static int iterc, prc;
693void print_iter(maxfnc, f)
694     int maxfnc;
695     double f;
696{
697        int i, j;
698
699        if(pr == 0)
700        {
701                iterc++;
702        }
703        else if(prc < abs(pr))
704        {
705                prc++;
706                iterc++;
707        }
708        else
709        {
710                printf("Iteration %6d", iterc);
711                printf(": function values %6d", maxfnc);
712                printf(" f = %24.16e\n", f);
713                if(pr < 0)
714                {
715                        printf("  d[] looks like this:\n");
716                        for(i=2;i<=n;i++)
717                        {
718                                printf("  ");
719                                for(j=1;j<i;j++)
720                                        printf("%24.16e ", d[i][j]);
721                                printf("\n  ");
722                        }
723                        printf("  g[] looks like this:\n");
724                        for(i=2;i<=n;i++)
725                        {
726                                printf("  ");
727                                for(j=1;j<i;j++)
728                                        printf("%24.16e ", g[i][j]);
729                                printf("\n  ");
730                        }
731                }
732                prc = 1;
733                iterc++;
734        }
735}
736
737void minfungra(dfn, acc, maxfn)
738     double dfn, acc;
739     int maxfn;
740{
741        double beta, c, clt, dginit, dgstep, dtest, finit;
742        double fmin, gamden, gamma, ggstar, gm, gmin, gnew;
743        double gsinit, gsumsq, xbound, xmin, xnew, xold;
744        int itcrs, flag, retry, maxfnk, maxfnc;
745        int i, j;
746
747        flag = 0;
748        retry = -2;
749        iterc = 0;
750        maxfnc = 1;
751        prc = abs(pr);
752        goto L190;
753
754 L10:
755        xnew = dabs(dfn+dfn)/gsumsq;
756        dgstep = -gsumsq;
757        itcrs = m+1;
758 L30:
759        if(maxfnk<maxfnc)
760        {
761                GLOBAL_f = fmin;
762                for(i=2;i<=n;i++)
763                        for(j=1;j<i;j++)
764                        {
765                                d[i][j] = xx[i][j];
766                                g[i][j] = gg[i][j];
767                        }
768        }
769        if(flag > 0)
770        {
771                prc = 0;
772                print_iter(maxfnc, GLOBAL_f);
773                return;
774        }
775        if(itcrs>m)
776        {
777                for(i=2;i<=n;i++)
778                        for(j=1;j<i;j++)
779                                d1[i][j] = -g[i][j];
780    }
781        print_iter(maxfnc, GLOBAL_f);
782        dginit = 0.0;
783        for(i=2;i<=n;i++)
784                for(j=1;j<i;j++)
785                {
786                        gs[i][j] = g[i][j];
787                        dginit += d1[i][j]*g[i][j];
788                }
789        if(dginit >= 0.0)
790        {
791                retry = -retry;
792                if(imin(retry, maxfnk)<2)
793                {
794                        printf("minfungra: \
795                                gradient wrong or acc too small\n");
796                        flag = 2;
797                }
798                else
799                        itcrs = m+1;
800                goto L30;
801        }
802        xmin = 0.0;
803        fmin = GLOBAL_f;
804        finit = GLOBAL_f;
805        gsinit = gsumsq;
806        gmin = dginit;
807        gm = dginit;
808        xbound = -1.0;
809        xnew = xnew * min(1.0, dgstep/dginit);
810        dgstep = dginit;
811 L170:
812        c = xnew-xmin;
813        dtest = 0.0;
814        for(i=1;i<=n;i++)
815                for(j=1;j<i;j++)
816                {
817                        d[i][j] = xx[i][j]+c*d1[i][j];
818                        dtest += dabs(d[i][j]-xx[i][j]);
819                }
820        if(dtest<=0.0)
821        {
822                clt = .7;
823                goto L300;
824        }
825        if(max(xbound, xmin-c) < 0.0) {
826        clt = 0.1;
827        }
828        maxfnc++;
829 L190:
830        fungra();
831
832        if(maxfnc>1)
833        {
834                for(gnew = 0.0,i=1;i<=n;i++)
835                        for(j=1;j<i;j++)
836                                gnew += d1[i][j]*g[i][j];
837    }
838        if(maxfnc<=1 ||
839           (maxfnc>1 && GLOBAL_f<fmin) ||
840           (maxfnc>1 && GLOBAL_f==fmin && dabs(gnew) <= dabs(gmin)))
841        {
842                maxfnk = maxfnc;
843                gsumsq = 0.0;
844                for(i=1;i<=n;i++)
845                        for(j=1;j<i;j++)
846                                gsumsq += sqr(g[i][j]);
847                if(gsumsq <= acc)
848                {
849                        prc = 0;
850                        print_iter(maxfnc, GLOBAL_f);
851                        return;
852                }
853        }
854        if(maxfnc==maxfn)
855        {
856                printf("minfungra: \
857                        maxfn function values have been calculated\n");
858                flag = 1;
859                goto L30;
860        }
861        if(maxfnk < maxfnc) goto L300;
862        for(i=1;i<=n;i++)
863                for(j=1;j<i;j++)
864                {
865                        xx[i][j] = d[i][j];
866                        gg[i][j] = g[i][j];
867                }
868        if(maxfnc <= 1) goto L10;
869        gm = gnew;
870        if(gm<=dginit) goto L310;
871        ggstar = 0.0;
872        for(i=1;i<=n;i++)
873                for(j=1;j<i;j++)
874                        ggstar += gs[i][j]*g[i][j];
875        beta = (gsumsq-ggstar)/(gm-dginit);
876 L300:
877        if(dabs(gm)>clt*dabs(dginit) ||
878           (dabs(gm)<=clt*dabs(dginit) && dabs(gm*beta) >= clt*gsumsq))
879        {
880    L310:
881                clt += 0.3;
882                if(clt>0.8)
883                {
884                        retry = -retry;
885                        if(imin(retry, maxfnk)<2)
886                        {
887                                printf("minfungra: \
888                                        gradient wrong or acc too small\n");
889                                flag = 2;
890                        }
891                        else
892                                itcrs = m+1;
893                        goto L30;
894                }
895                xold = xnew;
896                xnew = .5*(xmin+xold);
897                if(maxfnk >= maxfnc && gmin*gnew > 0.0)
898                {
899                        xnew = 10.0*xold;
900                        if(xbound>=0.0) {
901                xnew = 0.5*(xold+xbound);
902                        }
903                }
904                c = gnew-(3.0*gnew + gmin-4.0*(GLOBAL_f-fmin)/(xold-xmin))*
905            (xold-xnew)/(xold-xmin);
906                if(maxfnk>=maxfnc)
907                {
908                        if(gmin*gnew<=0.0) {
909                xbound = xmin;
910                        }
911                        xmin = xold;
912                        fmin = GLOBAL_f;
913                        gmin = gnew;
914                }
915                else
916                        xbound = xold;
917                if(c*gmin < 0.0)
918                        xnew = (xmin*c-xnew*gmin)/(c-gmin);
919                goto L170;
920        }
921        if(min(GLOBAL_f, fmin)<finit ||
922           (min(GLOBAL_f, fmin)>=finit && gsumsq < gsinit))
923        {
924                if(itcrs<m && dabs(ggstar) < .2*gsumsq)
925                {
926                        gamma = 0.0;
927                        c = 0.0;
928                        for(i=1;i<=n;i++)
929                                for(j=1;j<i;j++)
930                                {
931                                        gamma += gg[i][j] * dgr[i][j];
932                                        c += gg[i][j]*dr[i][j];
933                                }
934                        gamma /= gamden;
935                        if(dabs(beta*gm+gamma*c)<.2*gsumsq)
936                        {
937                                for(i=1;i<=n;i++)
938                                        for(j=1;j<i;j++)
939                                                d1[i][j] = -gg[i][j] +
940                            beta*d1[i][j] +
941                            gamma*dr[i][j];
942                                itcrs++;
943                        }
944                        else
945                        {
946                                itcrs = 2;
947                                gamden = gm-dginit;
948                                for(i=1;i<=n;i++)
949                                        for(j=1;j<i;j++)
950                                        {
951                                            dr[i][j] = d1[i][j];
952                                            dgr[i][j] = gg[i][j]-gs[i][j];
953                                            d1[i][j] = -gg[i][j]+beta*d1[i][j];
954                                        }
955                        }
956                }
957                else
958                {
959                        itcrs = 2;
960                        gamden = gm-dginit;
961                        for(i=1;i<=n;i++)
962                                for(j=1;j<i;j++)
963                                {
964                                        dr[i][j] = d1[i][j];
965                                        dgr[i][j] = gg[i][j]-gs[i][j];
966                                        d1[i][j] = -gg[i][j]+beta*d1[i][j];
967                                }
968                }
969                goto L30;
970        }
971        else
972        {
973                retry = -retry;
974                if(imin(retry, maxfnk)<2)
975                {
976                        printf("minfungra: \
977                                gradient wrong or acc too small\n");
978                        flag = 2;
979                }
980                else
981                        itcrs = m+1;
982                goto L30;
983        }
984}
985
986void get_dist()
987{
988        static double cons = 1.0;
989        static int maxfn = 500;
990        int i, j, iter;
991        double dif;
992
993
994        dr = (double **)getmem(n, sizeof(delta[1]));
995        dgr = (double **)getmem(n, sizeof(delta[1]));
996        d1 = (double **)getmem(n, sizeof(delta[1]));
997        gs = (double **)getmem(n, sizeof(delta[1]));
998        xx = (double **)getmem(n, sizeof(delta[1]));
999        gg = (double **)getmem(n, sizeof(delta[1]));
1000        dr[1] = NULL;
1001        dgr[1] = NULL;
1002        d1[1] = NULL;
1003        gs[1] = NULL;
1004        xx[1] = NULL;
1005        gg[1] = NULL;
1006        for(i=2; i<=n; i++)
1007        {
1008                dr[i] = (double *)getmem(i-1, sizeof(double));
1009                dgr[i] = (double *)getmem(i-1, sizeof(double));
1010                d1[i] = (double *)getmem(i-1, sizeof(double));
1011                gs[i] = (double *)getmem(i-1, sizeof(double));
1012                xx[i] = (double *)getmem(i-1, sizeof(double));
1013                gg[i] = (double *)getmem(i-1, sizeof(double));
1014        }
1015        nm0 = n;
1016        nm1 = n-1;
1017        nm2 = n-2;
1018        nm3 = n-3;
1019        if(start==3)
1020        {
1021                r = 0.001;
1022                fungra();
1023        }
1024        else
1025        {
1026                r = 1.0;
1027                fungra();
1028                r = cons*fitl/fitp;
1029                GLOBAL_f = (1.0+cons)*fitl;
1030        }
1031        iter=1;
1032        do
1033        {
1034                for(i=1;i<=n;i++)
1035                        for(j=1;j<i;j++)
1036                                dold[i][j] = d[i][j];
1037                minfungra(0.5*GLOBAL_f, ps1, maxfn);
1038                dif = 0.0;
1039                for(i=1;i<=n;i++)
1040                        for(j=1;j<i;j++)
1041                                dif +=sqr(d[i][j] - dold[i][j]);
1042                dif = sqrt(dif);
1043                if(dif>ps2)
1044                {
1045                        iter++;
1046                        r*=10.0;
1047                }
1048        } while (dif>ps2);
1049        fungra();
1050        for(i=2; i<=n; i++)
1051        {
1052                free(dr[i]);
1053                free(dgr[i]);
1054                free(d1[i]);
1055                free(gs[i]);
1056                free(xx[i]);
1057                free(gg[i]);
1058        }
1059        free(dr);
1060        free(dgr);
1061        free(d1);
1062        free(gs);
1063        free(xx);
1064        free(gg);
1065}
1066
1067double gttol()
1068{
1069        double result;
1070        int i, j, k, l;
1071
1072        result = 0.0;
1073        nm0 = n;
1074        nm1 = n-1;
1075        nm2 = n-2;
1076        nm3 = n-3;
1077        for(i=1;i<=nm3;i++)
1078            for(j=i+1;j<=nm2;j++)
1079            for(k=j+1;k<=nm1;k++)
1080                for(l=k+1;l<=nm0;l++)
1081                    if((d[j][i]+d[l][k]>=d[l][i]+d[k][j])&&
1082                       (d[k][i]+d[l][j]>=d[l][i]+d[k][j]))
1083                        result=max(result,
1084                                   dabs(d[j][i]+d[l][k]-d[k][i]-d[l][j]));
1085                    else
1086                        if((d[j][i]+d[l][k]>=d[k][i]+d[l][j])&&
1087                           (d[l][i]+d[k][j]>=d[k][i]+d[l][j]))
1088                            result=max(result,
1089                                       dabs(d[j][i]+d[l][k]-d[l][i]-d[k][j]));
1090                        else
1091                            if((d[k][i]+d[l][j]>=d[j][i]+d[l][k])&&
1092                               (d[l][i]+d[k][j]>=d[j][i]+d[l][k]))
1093                                result=max(result,
1094                                           dabs(d[k][i]+d[l][j]-d[l][i]-d[k][j]));
1095        return(result);
1096}
1097
1098void additive_const()
1099{
1100    int i, j, k;
1101
1102    double c=0.0;
1103    for(i=1;i<=n;i++)
1104        for(j=1;j<i;j++)
1105            for(k=1;k<=n;k++)
1106                if(k!=i && k!=j)
1107                    c=max(c,gd(i, j)-
1108                          gd( i, k)-
1109                          gd( j, k));
1110    if(report)
1111        fprintf(reportf, "Additive Constant: %15.10f\n", c);
1112    if(c!=0.0)
1113        for(i=1;i<=n;i++)
1114            for(j=1;j<i;j++)
1115                d[i][j]+=c;
1116}
1117
1118static int *mergei, *mergej;
1119static int ninv;
1120void get_tree()
1121{
1122#define false 0
1123#define true 1
1124        int i, j, k, l, im, jm, km, lm, count;
1125        int maxnode, maxcnt, nnode, nact;
1126        double arci, arcj, arcim, arcjm;
1127        char *act;
1128
1129        maxnode = 2*n-2;
1130        mergei = (int *)getmem(maxnode, sizeof(int));
1131        mergej = (int *)getmem(maxnode, sizeof(int));
1132        act = (char *)getmem(maxnode, sizeof(char));
1133        for(i=1;i<=maxnode;i++)
1134                act[i] = false;
1135        nact=n;
1136        ninv=0;
1137        nnode=n;
1138        while(nact>4)
1139        {
1140            maxcnt=-1;
1141            for(i=1;i<=nnode;i++) if(!act[i])
1142            for(j=1;j<i;j++) if(!act[j])
1143            {
1144                count=0;
1145                arci=0.0;
1146                arcj=0.0;
1147                for(k=2;k<=nnode;k++) if(!act[k] && k!=i && k!=j)
1148                    for(l=1;l<k;l++) if(!act[l] && l!= i && l!= j)
1149                        if(gd(i,j)+gd(k,l)<=gd(i,k)+gd(j,l) &&
1150                           gd(i,j)+gd(k,l)<=gd(i,l)+gd(j,k))
1151                        {
1152                            count++;
1153                            arci+=(gd(i,j)+gd(i,k)-gd(j,k))/2.0;
1154                            arcj+=(gd(i,j)+gd(j,l)-gd(i,l))/2.0;
1155                        }
1156                if(count>maxcnt)
1157                {
1158                    maxcnt = count;
1159                    arcim=max(0.0, arci/count);
1160                    arcjm=max(0.0, arcj/count);
1161                    im=i;
1162                    jm=j;
1163                }
1164            }
1165
1166            nnode++;
1167            if(nnode+2>maxnode)
1168                    show_error("get_tree: number of nodes exceeds 2N-2");
1169            ninv++;
1170            mergei[ninv]=im;
1171            mergej[ninv]=jm;
1172            act[im]=true;
1173            act[jm]=true;
1174            d[nnode]=(double *)getmem(nnode-1, sizeof(double));
1175            d[nnode][im]=arcim;
1176            d[nnode][jm]=arcjm;
1177            for(i=1;i<=nnode-1;i++)
1178                    if(!act[i])
1179                            d[nnode][i] = max(0.0, gd(im,i)-arcim);
1180            nact--;
1181        }
1182
1183        for(i=1;act[i];i++)
1184            if(i>nnode)
1185            show_error("get_tree: can't find last two invisible nodes");
1186        im=i;
1187        for(i=im+1;act[i];i++)
1188            if(i>nnode)
1189            show_error("get_tree: can't find last two invisible nodes");
1190        jm=i;
1191        for(i=jm+1;act[i];i++)
1192            if(i>nnode)
1193            show_error("get_tree: can't find last two invisible nodes");
1194        km=i;
1195        for(i=km+1;act[i];i++)
1196            if(i>nnode)
1197            show_error("get_tree: can't find last two invisible nodes");
1198        lm=i;
1199        if(gd(im,jm)+gd(km,lm)<=gd(im,km)+gd(jm,lm)+tol &&
1200           gd(im,jm)+gd(km,lm)<=gd(im,lm)+gd(jm,km)+tol)
1201        {
1202                i=im;
1203                j=jm;
1204                k=km;
1205                l=lm;
1206        }
1207        else if(gd(im,lm)+gd(jm,km)<=gd(im,km)+gd(jm,lm)+tol &&
1208                gd(im,lm)+gd(jm,km)<=gd(im,jm)+gd(km,lm)+tol)
1209        {
1210                i=im;
1211                j=lm;
1212                k=km;
1213                l=jm;
1214        }
1215        else if(gd(im,km)+gd(jm,lm)<=gd(im,jm)+gd(km,lm)+tol &&
1216                gd(im,km)+gd(jm,lm)<=gd(im,lm)+gd(jm,km)+tol)
1217        {
1218                i=im;
1219                j=km;
1220                k=lm;
1221                l=jm;
1222        }
1223
1224        nnode++;
1225        ninv++;
1226        mergei[ninv]=i;
1227        mergej[ninv]=j;
1228        d[nnode]=(double *)getmem(nnode-1, sizeof(double));
1229        d[nnode][i] = max(0.0, (gd(i,j)+gd(i,k)-gd(j,k))/2.0);
1230        d[nnode][j] = max(0.0, (gd(i,j)+gd(j,l)-gd(i,l))/2.0);
1231        nnode++;
1232        ninv++;
1233        mergei[ninv]=k;
1234        mergej[ninv]=l;
1235        d[nnode]=(double *)getmem(nnode-1, sizeof(double));
1236        d[nnode][k] = max(0.0, (gd(k,l)+gd(i,k)-gd(i,l))/2.0);
1237        d[nnode][l] = max(0.0, (gd(k,l)+gd(j,l)-gd(j,k))/2.0);
1238        d[nnode][nnode-1] = max(0.0, (gd(i,k)+gd(j,l)-gd(i,j)-gd(k,l))/2.0);
1239
1240}
1241
1242void print_node(node, dist, indent)
1243     int node;
1244     double dist;
1245     int indent;
1246{
1247        static char buf[BUFLEN];
1248
1249        if(node<=n)
1250                printf("%s%s:%6.4f",
1251               repeatch(buf, '\t', indent),
1252               names[node],
1253               dist/nfac);
1254        else
1255        {
1256                printf("%s(\n",
1257               repeatch(buf, '\t', indent));
1258                print_node(mergei[node-n], gd(node, mergei[node-n]), indent+1);
1259                printf(",\n");
1260                print_node(mergej[node-n], gd(node, mergej[node-n]), indent+1);
1261                printf("\n%s):%6.4f", repeatch(buf, '\t', indent),
1262               dist/nfac);
1263        }
1264}
1265
1266void show_tree()
1267{
1268        int i, j, current;
1269        int ij[2];
1270
1271        current=0;
1272        for(i=1;current<2;i++)
1273        {
1274                for(j=1;(mergei[j]!=i && mergej[j] != i) && j<=ninv;j++);
1275                if(j>ninv)
1276                        ij[current++]=i;
1277        }
1278        printf("(\n");
1279        print_node(ij[0], gd(ij[0],ij[1])/2.0, 1);
1280        printf(",\n");
1281        print_node(ij[1], gd(ij[0],ij[1])/2.0, 1);
1282        printf("\n);\n");
1283}
1284
1285int main(argc, argv)
1286     int argc;
1287     char **argv;
1288{
1289        int i;
1290
1291        strcpy(fname, "");
1292        get_parms(argc, argv);
1293        if(strcmp(fname, ""))
1294        {
1295                report=1;
1296                reportf = fopen(fname, "w");
1297                if(reportf==NULL)
1298                {
1299                        perror("lsadt");
1300                        return 0;
1301                }
1302
1303        }
1304        else
1305                report=0;
1306        get_data();
1307        d = (double **)getmem(n, sizeof(delta[1]));
1308        g = (double **)getmem(n, sizeof(delta[1]));
1309        dold = (double **)getmem(n, sizeof(delta[1]));
1310        d[1]=NULL;
1311        g[1]=NULL;
1312        dold[1]=NULL;
1313        for(i=2; i<=n; i++)
1314        {
1315                d[i]=(double *)getmem(i-1, sizeof(double));
1316                g[i]=(double *)getmem(i-1, sizeof(double));
1317                dold[i]=(double *)getmem(i-1, sizeof(double));
1318        }
1319        initd();
1320        get_dist();
1321        tol=gttol();
1322        additive_const();
1323        get_tree();
1324        show_tree();
1325        if(report)
1326            fclose(reportf);
1327}
Note: See TracBrowser for help on using the repository browser.