The R Project SVN R

Rev

Rev 37468 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 37468 Rev 37513
Line 76... Line 76...
76 
	REAL(x)[i] = p[i] * (OS->parscale[i]);
 76 
	REAL(x)[i] = p[i] * (OS->parscale[i]);
77 
    }
 77 
    }
78 
    SETCADR(OS->R_fcall, x);
 78 
    SETCADR(OS->R_fcall, x);
79 
    PROTECT_WITH_INDEX(s = eval(OS->R_fcall, OS->R_env), &ipx);
 79 
    PROTECT_WITH_INDEX(s = eval(OS->R_fcall, OS->R_env), &ipx);
80 
    REPROTECT(s = coerceVector(s, REALSXP), ipx);
 80 
    REPROTECT(s = coerceVector(s, REALSXP), ipx);
81 
    if (LENGTH(s) != 1)
81 
    if (LENGTH(s) != 1)
82 
	error(_("objective function in optim evaluates to length %d not 1"),
82 
	error(_("objective function in optim evaluates to length %d not 1"),
83 
	      LENGTH(s));
83 
	      LENGTH(s));
84 
    val = REAL(s)[0]/(OS->fnscale);
 84 
    val = REAL(s)[0]/(OS->fnscale);
85 
    UNPROTECT(2);
 85 
    UNPROTECT(2);
86 
    return val;
 86 
    return val;
87 
}
 87 
}
88 
 
 88 
 
Line 170... Line 170...
170 
	UNPROTECT(1); /* x */
 170 
	UNPROTECT(1); /* x */
171 
    }
 171 
    }
172 
}
 172 
}
173 
 
 173 
 
174 
static void genptry(int n, double *p, double *ptry, double scale, void *ex)
 174 
static void genptry(int n, double *p, double *ptry, double scale, void *ex)
175 
{
175 
{
176 
    SEXP s, x;
 176 
    SEXP s, x;
177 
    int i;
 177 
    int i;
178 
    OptStruct OS = (OptStruct) ex;
 178 
    OptStruct OS = (OptStruct) ex;
179 
    PROTECT_INDEX ipx;
 179 
    PROTECT_INDEX ipx;
180 
 
 180 
 
181 
    if (!isNull(OS->R_gcall)) {
181 
    if (!isNull(OS->R_gcall)) {
182 
	/* user defined generation of candidate point */
 182 
	/* user defined generation of candidate point */
183 
      	PROTECT(x = allocVector(REALSXP, n));
 183 
      	PROTECT(x = allocVector(REALSXP, n));
184 
	for (i = 0; i < n; i++) {
 184 
	for (i = 0; i < n; i++) {
185 
	    if (!R_FINITE(p[i]))
 185 
	    if (!R_FINITE(p[i]))
186 
		error(_("non-finite value supplied by optim"));
 186 
		error(_("non-finite value supplied by optim"));
Line 193... Line 193...
193 
	    error(_("candidate point in optim evaluated to length %d not %d"),
 193 
	    error(_("candidate point in optim evaluated to length %d not %d"),
194 
		  LENGTH(s), n);
 194 
		  LENGTH(s), n);
195 
	for (i = 0; i < n; i++)
 195 
	for (i = 0; i < n; i++)
196 
	    ptry[i] = REAL(s)[i] / (OS->parscale[i]);
 196 
	    ptry[i] = REAL(s)[i] / (OS->parscale[i]);
197 
	UNPROTECT(2);
 197 
	UNPROTECT(2);
198 
    }
198 
    }
199 
    else {  /* default Gaussian Markov kernel */
 199 
    else {  /* default Gaussian Markov kernel */
200 
        for (i = 0; i < n; i++)
 200 
        for (i = 0; i < n; i++)
201 
            ptry[i] = p[i] + scale * norm_rand();  /* new candidate point */
 201 
            ptry[i] = p[i] + scale * norm_rand();  /* new candidate point */
202 
    }
 202 
    }
203 
}
 203 
}
Line 265... Line 265...
265 
	for (i = 0; i < npar; i++)
 265 
	for (i = 0; i < npar; i++)
266 
	    REAL(par)[i] = opar[i] * (OS->parscale[i]);
 266 
	    REAL(par)[i] = opar[i] * (OS->parscale[i]);
267 
	grcount = NA_INTEGER;
 267 
	grcount = NA_INTEGER;
268 
 
 268 
 
269 
    }
 269 
    }
270 
    else if (strcmp(tn, "SANN") == 0) {
270 
    else if (strcmp(tn, "SANN") == 0) {
271 
        tmax = asInteger(getListElement(options, "tmax"));
 271 
        tmax = asInteger(getListElement(options, "tmax"));
272 
        temp = asReal(getListElement(options, "temp"));
 272 
        temp = asReal(getListElement(options, "temp"));
273 
        if (tmax == NA_INTEGER) error(_("'tmax' is not an integer"));
 273 
        if (tmax == NA_INTEGER) error(_("'tmax' is not an integer"));
274 
        if (!isNull(gr)) {
 274 
        if (!isNull(gr)) {
275 
            if (!isFunction(gr)) error(_("'gr' is not a function"));
 275 
            if (!isFunction(gr)) error(_("'gr' is not a function"));
Line 651... Line 651...
651 
	   double alpha, double bet, double gamm, int trace,
 651 
	   double alpha, double bet, double gamm, int trace,
652 
	   int *fncount, int maxit)
 652 
	   int *fncount, int maxit)
653 
{
 653 
{
654 
    char action[50];
 654 
    char action[50];
655 
    int C;
 655 
    int C;
656 
    Rboolean calcvert, shrinkfail = FALSE;
 656 
    Rboolean calcvert;
657 
    double convtol, f;
 657 
    double convtol, f;
658 
    int funcount=0, H, i, j, L=0;
 658 
    int funcount=0, H, i, j, L=0;
659 
    int n1=0;
 659 
    int n1=0;
660 
    double oldsize;
 660 
    double oldsize;
661 
    double **P;
 661 
    double **P;
Line 710... Line 710...
710 
	    }
 710 
	    }
711 
	    size += trystep;
 711 
	    size += trystep;
712 
	}
 712 
	}
713 
	oldsize = size;
 713 
	oldsize = size;
714 
	calcvert = TRUE;
 714 
	calcvert = TRUE;
715 
	shrinkfail = FALSE;
 - 
 
716 
	do {
 715 
	do {
717 
	    if (calcvert) {
 716 
	    if (calcvert) {
718 
		for (j = 0; j < n1; j++) {
 717 
		for (j = 0; j < n1; j++) {
719 
		    if (j + 1 != L) {
 718 
		    if (j + 1 != L) {
720 
			for (i = 0; i < n; i++)
 719 
			for (i = 0; i < n; i++)
Line 744... Line 743...
744 
			VH = f;
 743 
			VH = f;
745 
		    }
 744 
		    }
746 
		}
 745 
		}
747 
	    }
 746 
	    }
748 
 
 747 
 
749 
	    if (VH > VL + convtol && VL > abstol) {
 748 
	    if (VH <= VL + convtol || VL <= abstol) break;
750 
		sprintf(tstr, "%5d", funcount);
 - 
 
751 
		if (trace) Rprintf("%s%s %f %f\n", action, tstr, VH, VL);
 - 
 
752 
 
 749 
 
- 
 
 750 
	    sprintf(tstr, "%5d", funcount);
- 
 
 751 
	    if (trace) Rprintf("%s%s %f %f\n", action, tstr, VH, VL);
- 
 
 752 
 
- 
 
 753 
	    for (i = 0; i < n; i++) {
- 
 
 754 
		temp = -P[i][H - 1];
- 
 
 755 
		for (j = 0; j < n1; j++)
- 
 
 756 
		    temp += P[i][j];
- 
 
 757 
		P[i][C - 1] = temp / n;
- 
 
 758 
	    }
- 
 
 759 
	    for (i = 0; i < n; i++)
- 
 
 760 
		Bvec[i] = (1.0 + alpha) * P[i][C - 1] - alpha * P[i][H - 1];
- 
 
 761 
	    f = fminfn(n, Bvec, ex);
- 
 
 762 
	    if (!R_FINITE(f)) f = big;
- 
 
 763 
	    funcount++;
- 
 
 764 
	    strcpy(action, "REFLECTION     ");
- 
 
 765 
	    VR = f;
- 
 
 766 
	    if (VR < VL) {
- 
 
 767 
		P[n1 - 1][C - 1] = f;
753 
		for (i = 0; i < n; i++) {
 768 
		for (i = 0; i < n; i++) {
754 
		    temp = -P[i][H - 1];
 769 
		    f = gamm * Bvec[i] + (1 - gamm) * P[i][C - 1];
755 
		    for (j = 0; j < n1; j++)
 770 
		    P[i][C - 1] = Bvec[i];
756 
			temp += P[i][j];
 771 
		    Bvec[i] = f;
757 
		    P[i][C - 1] = temp / n;
 - 
 
758 
		}
 772 
		}
759 
		for (i = 0; i < n; i++)
 - 
 
760 
		    Bvec[i] = (1.0 + alpha) * P[i][C - 1] - alpha * P[i][H - 1];
 - 
 
761 
		f = fminfn(n, Bvec, ex);
 773 
		f = fminfn(n, Bvec, ex);
762 
		if (!R_FINITE(f)) f = big;
 774 
		if (!R_FINITE(f)) f = big;
763 
		funcount++;
 775 
		funcount++;
764 
		strcpy(action, "REFLECTION     ");
 - 
 
765 
		VR = f;
 - 
 
766 
		if (VR < VL) {
 776 
		if (f < VR) {
767 
		    P[n1 - 1][C - 1] = f;
 - 
 
768 
		    for (i = 0; i < n; i++) {
 777 
		    for (i = 0; i < n; i++)
769 
			f = gamm * Bvec[i] + (1 - gamm) * P[i][C - 1];
 - 
 
770 
			P[i][C - 1] = Bvec[i];
 778 
			P[i][H - 1] = Bvec[i];
771 
			Bvec[i] = f;
 - 
 
772 
		    }
 - 
 
773 
		    f = fminfn(n, Bvec, ex);
 - 
 
774 
		    if (!R_FINITE(f)) f = big;
 - 
 
775 
		    funcount++;
 - 
 
776 
		    if (f < VR) {
 - 
 
777 
			for (i = 0; i < n; i++)
 - 
 
778 
			    P[i][H - 1] = Bvec[i];
 - 
 
779 
			P[n1 - 1][H - 1] = f;
 779 
		    P[n1 - 1][H - 1] = f;
780 
			strcpy(action, "EXTENSION      ");
 780 
		    strcpy(action, "EXTENSION      ");
781 
		    } else {
 - 
 
782 
			for (i = 0; i < n; i++)
 - 
 
783 
			    P[i][H - 1] = P[i][C - 1];
 - 
 
784 
			P[n1 - 1][H - 1] = VR;
 - 
 
785 
		    }
 - 
 
786 
		} else {
 781 
		} else {
787 
		    strcpy(action, "HI-REDUCTION   ");
 - 
 
788 
		    if (VR < VH) {
 - 
 
789 
			for (i = 0; i < n; i++)
 - 
 
790 
			    P[i][H - 1] = Bvec[i];
 - 
 
791 
			P[n1 - 1][H - 1] = VR;
 - 
 
792 
			strcpy(action, "LO-REDUCTION   ");
 - 
 
793 
		    }
 - 
 
794 
 
 - 
 
795 
		    for (i = 0; i < n; i++)
 782 
		    for (i = 0; i < n; i++)
796 
			Bvec[i] = (1 - bet) * P[i][H - 1] + bet * P[i][C - 1];
 783 
			P[i][H - 1] = P[i][C - 1];
797 
		    f = fminfn(n, Bvec, ex);
 784 
		    P[n1 - 1][H - 1] = VR;
- 
 
 785 
		}
- 
 
 786 
	    } else {
- 
 
 787 
		strcpy(action, "HI-REDUCTION   ");
- 
 
 788 
		if (VR < VH) {
798 
		    if (!R_FINITE(f)) f = big;
 789 
		    for (i = 0; i < n; i++)
- 
 
 790 
			P[i][H - 1] = Bvec[i];
799 
		    funcount++;
 791 
		    P[n1 - 1][H - 1] = VR;
- 
 
 792 
		    strcpy(action, "LO-REDUCTION   ");
- 
 
 793 
		}
800 
 
 794 
 
- 
 
 795 
		for (i = 0; i < n; i++)
- 
 
 796 
		    Bvec[i] = (1 - bet) * P[i][H - 1] + bet * P[i][C - 1];
- 
 
 797 
		f = fminfn(n, Bvec, ex);
- 
 
 798 
		if (!R_FINITE(f)) f = big;
- 
 
 799 
		funcount++;
- 
 
 800 
 
801 
		    if (f < P[n1 - 1][H - 1]) {
 801 
		if (f < P[n1 - 1][H - 1]) {
802 
			for (i = 0; i < n; i++)
 802 
		    for (i = 0; i < n; i++)
803 
			    P[i][H - 1] = Bvec[i];
 803 
			P[i][H - 1] = Bvec[i];
804 
			P[n1 - 1][H - 1] = f;
 804 
		    P[n1 - 1][H - 1] = f;
805 
		    } else {
 805 
		} else {
806 
			if (VR >= VH) {
 806 
		    if (VR >= VH) {
807 
			    strcpy(action, "SHRINK         ");
 807 
			strcpy(action, "SHRINK         ");
808 
			    calcvert = TRUE;
 808 
			calcvert = TRUE;
809 
			    size = 0.0;
 809 
			size = 0.0;
810 
			    for (j = 0; j < n1; j++) {
 810 
			for (j = 0; j < n1; j++) {
811 
				if (j + 1 != L) {
 811 
			    if (j + 1 != L) {
812 
				    for (i = 0; i < n; i++) {
 812 
				for (i = 0; i < n; i++) {
813 
					P[i][j] = bet * (P[i][j] - P[i][L - 1])
 813 
				    P[i][j] = bet * (P[i][j] - P[i][L - 1])
814 
					    + P[i][L - 1];
 814 
					+ P[i][L - 1];
815 
					size += fabs(P[i][j] - P[i][L - 1]);
 815 
				    size += fabs(P[i][j] - P[i][L - 1]);
816 
				    }
 - 
 
817 
				}
 816 
				}
818 
			    }
 817 
			    }
- 
 
 818 
			}
819 
			    if (size < oldsize) {
 819 
			if (size < oldsize) {
820 
				shrinkfail = FALSE;
 - 
 
821 
				oldsize = size;
 820 
			    oldsize = size;
822 
			    } else {
 821 
			} else {
823 
				if (trace)
 822 
			    if (trace)
824 
				    Rprintf("Polytope size measure not decreased in shrink\n");
 823 
				Rprintf("Polytope size measure not decreased in shrink\n");
825 
				shrinkfail = TRUE;
 824 
			    *fail = 10;
826 
			    }
 825 
			    break;
827 
			}
 826 
			}
828 
		    }
 827 
		    }
829 
		}
 828 
		}
830 
	    }
 829 
	    }
831 
 
 830 
 
832 
	} while (!(VH <= VL + convtol || VL <= abstol ||
 - 
 
833 
		   shrinkfail || funcount > maxit));
 831 
	} while (funcount <= maxit);
834 
 
 832 
 
835 
    }
 833 
    }
836 
 
 834 
 
837 
    if (trace) {
 835 
    if (trace) {
838 
	Rprintf("Exiting from Nelder Mead minimizer\n");
 836 
	Rprintf("Exiting from Nelder Mead minimizer\n");
839 
	Rprintf("    %d function evaluations used\n", funcount);
 837 
	Rprintf("    %d function evaluations used\n", funcount);
840 
    }
 838 
    }
841 
    *Fmin = P[n1 - 1][L - 1];
 839 
    *Fmin = P[n1 - 1][L - 1];
842 
    for (i = 0; i < n; i++) X[i] = P[i][L - 1];
 840 
    for (i = 0; i < n; i++) X[i] = P[i][L - 1];
843 
    if (shrinkfail) *fail = 10;
 - 
 
844 
    if (funcount > maxit) *fail = 1;
 841 
    if (funcount > maxit) *fail = 1;
845 
    *fncount = funcount;
 842 
    *fncount = funcount;
846 
}
 843 
}
847 
 
 844 
 
848 
void cgmin(int n, double *Bvec, double *X, double *Fmin,
 845 
void cgmin(int n, double *Bvec, double *X, double *Fmin,