source: branches/items/GDE/SUPPORT/lsadt.c

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