source: tags/initial/GDE/SUPPORT/lsadt.c

Last change on this file was 2, checked in by oldcode, 24 years ago

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