The R Project SVN R

Rev

Rev 12976 | Rev 13942 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 12976 Rev 13214
Line 27... Line 27...
27 
 
 27 
 
28 
/* "GetRowNames" and "GetColNames" are utility routines which */
 28 
/* "GetRowNames" and "GetColNames" are utility routines which */
29 
/* locate and return the row names and column names from the */
 29 
/* locate and return the row names and column names from the */
30 
/* dimnames attribute of a matrix.  They are useful because */
 30 
/* dimnames attribute of a matrix.  They are useful because */
31 
/* old versions of R used pair-based lists for dimnames */
 31 
/* old versions of R used pair-based lists for dimnames */
32 
/* whereas recent versions use vector bassed lists */
 32 
/* whereas recent versions use vector based lists */
33 
 
 33 
 
34 
/* FIXME : This is nonsense.  When the "dimnames" attribute is */
 34 
/* FIXME : This is nonsense.  When the "dimnames" attribute is */
35 
/* grabbed off an array it is always adjusted to be a vector. */
 35 
/* grabbed off an array it is always adjusted to be a vector. */
36 
 
 36 
 
37 
SEXP GetRowNames(SEXP dimnames)
 37 
SEXP GetRowNames(SEXP dimnames)
Line 699... Line 699...
699 
    errorcall(call, "argument is not a matrix");
 699 
    errorcall(call, "argument is not a matrix");
700 
    return call;/* never used; just for -Wall */
 700 
    return call;/* never used; just for -Wall */
701 
}
 701 
}
702 
 
 702 
 
703 
 
 703 
 
- 
 
 704 
#ifdef OLD
704 
/* swap works by finding for a index i, the position */
 705 
/* swap works by finding for a index i, the position */
705 
/* in the array with dimensions dims1 in terms of */
 706 
/* in the array with dimensions dims1 in terms of */
706 
/* (i, j, k, l, m, ...); i.e. component-wise position, */
 707 
/* (i, j, k, l, m, ...); i.e. component-wise position, */
707 
/* then permute these to the order of the array with */
 708 
/* then permute these to the order of the array with */
708 
/* dimensions dims2 and work backwards to an integer */
 709 
/* dimensions dims2 and work backwards to an integer */
Line 727... Line 728...
727 
	t1 *= INTEGER(dims2)[i];
 728 
	t1 *= INTEGER(dims2)[i];
728 
	t1 += INTEGER(ind2)[i];
 729 
	t1 += INTEGER(ind2)[i];
729 
    }
 730 
    }
730 
    return t1;
 731 
    return t1;
731 
}
 732 
}
- 
 
 733 
#endif
732 
 
 734 
 
- 
 
 735 
/*
- 
 
 736 
 New version of aperm, using strides for speed.
- 
 
 737 
 Jonathan Rougier <J.C.Rougier@durham.ac.uk>
- 
 
 738 
 
- 
 
 739 
 v1.0 30.01.01
- 
 
 740 
 
- 
 
 741 
 M.Maechler : expanded	all ../include/Rdefines.h macros
- 
 
 742 
 */
- 
 
 743 
 
- 
 
 744 
/* this increments iip and sets j using strides */
- 
 
 745 
 
- 
 
 746 
#define CLICKJ						\
- 
 
 747 
    for (itmp=0; itmp<n; itmp++)			\
- 
 
 748 
	if (iip[itmp] == INTEGER(dimsr)[itmp]-1)	\
- 
 
 749 
	    iip[itmp] = 0;				\
- 
 
 750 
	else {						\
- 
 
 751 
	    iip[itmp]++;				\
- 
 
 752 
	    break;					\
- 
 
 753 
	}						\
- 
 
 754 
    for (j=0, itmp=0; itmp<n; itmp++)			\
- 
 
 755 
	j += iip[itmp] * stride[itmp];
- 
 
 756
- 
 
 757 
/* aperm (a, perm, resize = TRUE) */
733 
SEXP do_aperm(SEXP call, SEXP op, SEXP args, SEXP rho)
 758 
SEXP do_aperm(SEXP call, SEXP op, SEXP args, SEXP rho)
734 
{
 759 
{
735 
    SEXP a, perm, r, dimsa, dimsr, ind1, ind2;
 760 
    SEXP a, perm, resize, r, dimsa, dimsr, dna;
736 
    int i, j, len;
 761 
    int i, j, n, len, itmp;
- 
 
 762 
    int *pp, *iip, *stride;
- 
 
 763 
    char *vmax;
737 
 
 764 
 
738 
    checkArity(op, args);
 765 
    checkArity(op, args);
739 
 
 766 
 
740 
    a = CAR(args);
 767 
    a = CAR(args);
- 
 
 768 
    if (!isArray(a))
- 
 
 769 
	errorcall(call,"invalid first argument, must be an array");
- 
 
 770 
 
741 
    PROTECT(dimsa = getAttrib(a, R_DimSymbol));
 771 
    PROTECT(dimsa = getAttrib(a, R_DimSymbol));
742 
    if (dimsa == R_NilValue)
 772 
    n = LENGTH(dimsa);
- 
 
 773 
 
743 
	error("aperm: invalid first argument, must be an array");
 774 
    /* check the permutation */
744 
 
 775 
 
745 
    PROTECT(perm = coerceVector(CADR(args), INTSXP));
 776 
    PROTECT(perm = coerceVector(CADR(args), INTSXP));
- 
 
 777 
#ifdef OLD
746 
    if (!isVector(perm) || (length(perm) != length(dimsa)))
 778 
    if (!isVector(perm) || (length(perm) != length(dimsa)))
- 
 
 779 
	errorcall(call,
747 
	error("aperm: invalid second argument, must be a vector of the appropriate length");
 780 
	 "invalid second argument, must be a vector of the appropriate length");
- 
 
 781 
#else
- 
 
 782 
    vmax = vmaxget();
- 
 
 783 
    pp = (int *) R_alloc(n, sizeof(int));
- 
 
 784 
    if (length(perm) == 0) {
- 
 
 785 
	for (i=0; i<n; i++)
- 
 
 786 
	    pp[i] = n-1-i;
- 
 
 787 
    } else if (length(perm) == n) {
- 
 
 788 
	for (i=0; i<n; i++)
- 
 
 789 
	    pp[i] = INTEGER(perm)[i] - 1; /* no offset! */
- 
 
 790 
    } else
- 
 
 791 
	errorcall(call, "`perm' is of wrong length");
- 
 
 792 
 
- 
 
 793 
    iip = (int *) R_alloc(n, sizeof(int));
- 
 
 794 
    for (i=0; i<n; iip[i++] = 0);
- 
 
 795 
    for (i=0; i<n; i++)
- 
 
 796 
	if (pp[i] >= 0 && pp[i] < n)
- 
 
 797 
	    iip[pp[i]]++;
- 
 
 798 
	else
- 
 
 799 
	    errorcall(call, "value of out range in `perm'");
- 
 
 800 
    for (i=0; i<n; i++)
- 
 
 801 
	if (iip[i]==0)
- 
 
 802 
	    errorcall(call, "invalid permutation (`perm')");
- 
 
 803 
#endif
748 
 
 804 
 
- 
 
 805 
    /* create the stride object and permute */
- 
 
 806 
 
749 
    len = length(a);
 807 
    stride = (int *) R_alloc(n, sizeof(int));
- 
 
 808 
 
- 
 
 809 
    for (iip[0] = 1, i = 1; i<n; i++)
- 
 
 810 
	iip[i] = iip[i-1] * INTEGER(dimsa)[i-1];
750 
 
 811 
 
751 
    PROTECT(dimsr = allocVector(INTSXP, length(dimsa)));
 - 
 
752 
    for (i = 0; i < length(dimsa); i++)
 812 
    for (i=0; i<n; i++)
753 
	INTEGER(dimsr)[i] = INTEGER(dimsa)[(INTEGER(perm)[i] - 1)];
 813 
	stride[i] = iip[pp[i]];
754 
 
 814 
 
- 
 
 815 
    /* also need to have the dimensions of r */
- 
 
 816 
 
- 
 
 817 
    PROTECT(dimsr = allocVector(INTSXP,n));
- 
 
 818 
    for (i=0; i<n; i++)
- 
 
 819 
	INTEGER(dimsr)[i] = INTEGER(dimsa)[pp[i]];
- 
 
 820 
 
- 
 
 821 
    /* and away we go! iip will hold the incrementer */
- 
 
 822 
 
- 
 
 823 
    len = LENGTH(a);
- 
 
 824 
    len = length(a);
755 
    PROTECT(r = allocVector(TYPEOF(a), len));
 825 
    PROTECT(r = allocVector(TYPEOF(a), len));
756 
    PROTECT(ind1 = allocVector(INTSXP, LENGTH(dimsa)));
 - 
 
- 
 
 826 
 
757 
    PROTECT(ind2 = allocVector(INTSXP, LENGTH(dimsa)));
 827 
    for (i=0; i<n; iip[i++] = 0);
758 
 
 828 
 
759 
    switch (TYPEOF(a)) {
 829 
    switch (TYPEOF(a)) {
- 
 
 830 
 
760 
    case INTSXP:
 831 
    case INTSXP:
761 
    case LGLSXP:
 832 
    case LGLSXP:
762 
	for (i = 0; i < len; i++) {
 833 
	for (j=0, i=0; i<len; i++) {
763 
	    j = swap(i, dimsa, dimsr, perm, ind1, ind2);
 - 
 
764 
	    INTEGER(r)[j] = INTEGER(a)[i];
 834 
	    INTEGER(r)[i] = INTEGER(a)[j];
- 
 
 835 
	    CLICKJ;
765 
	}
 836 
	}
766 
	break;
 837 
	break;
- 
 
 838 
 
767 
    case REALSXP:
 839 
    case REALSXP:
768 
	for (i = 0; i < len; i++) {
 840 
	for (j=0, i=0; i<len; i++) {
769 
	    j = swap(i, dimsa, dimsr, perm, ind1, ind2);
 - 
 
770 
	    REAL(r)[j] = REAL(a)[i];
 841 
	    REAL(r)[i] = REAL(a)[j];
- 
 
 842 
	    CLICKJ;
771 
	}
 843 
	}
772 
	break;
 844 
	break;
- 
 
 845 
 
773 
    case CPLXSXP:
 846 
    case CPLXSXP:
774 
	for (i = 0; i < len; i++) {
 847 
	for (j=0, i=0; i<len; i++) {
775 
	    j = swap(i, dimsa, dimsr, perm, ind1, ind2);
 848 
	    COMPLEX(r)[i].r = COMPLEX(a)[j].r;
776 
	    COMPLEX(r)[j] = COMPLEX(a)[i];
 849 
	    COMPLEX(r)[i].i = COMPLEX(a)[j].i;
- 
 
 850 
	    CLICKJ;
777 
	}
 851 
	}
778 
	break;
 852 
	break;
- 
 
 853 
 
779 
    case STRSXP:
 854 
    case STRSXP:
780 
	for (i = 0; i < len; i++) {
 855 
	for (j=0, i=0; i<len; i++) {
781 
	    j = swap(i, dimsa, dimsr, perm, ind1, ind2);
 - 
 
782 
	    SET_STRING_ELT(r, j, STRING_ELT(a, i));
 856 
	    SET_STRING_ELT(r, i, STRING_ELT(a, j));
- 
 
 857 
	    CLICKJ;
783 
	}
 858 
	}
784 
	break;
 859 
	break;
- 
 
 860 
 
785 
    case VECSXP:
 861 
    case VECSXP:
786 
	for (i = 0; i < len; i++) {
 862 
	for (j=0, i=0; i<len; i++) {
787 
	    j = swap(i, dimsa, dimsr, perm, ind1, ind2);
 - 
 
788 
	    SET_VECTOR_ELT(r, j, VECTOR_ELT(a, i));
 863 
	    SET_VECTOR_ELT(r, i, VECTOR_ELT(a, j));
- 
 
 864 
	    CLICKJ;
789 
	}
 865 
	}
790 
	break;
 866 
	break;
- 
 
 867 
 
791 
    default:
 868 
    default:
792 
	errorcall(call, R_MSG_IA);
 869 
	errorcall(call, "unsupported type of array");
793 
    }
 870 
    }
794 
 
 871
- 
 
 872 
    /* handle the resize */
- 
 
 873 
    PROTECT(resize = coerceVector(CADDR(args), INTSXP));
795 
    if (INTEGER(CAR(CDDR(args)))[0])
 874 
    if (LOGICAL(resize)[0])
796 
	setAttrib(r, R_DimSymbol, dimsr);
 875 
	setAttrib(r, R_DimSymbol, dimsr);
797 
    else
 876 
    else
798 
	setAttrib(r, R_DimSymbol, dimsa);
 877 
	setAttrib(r, R_DimSymbol, dimsa);
- 
 
 878 
 
- 
 
 879 
    /* and handle the dimnames */
- 
 
 880 
 
- 
 
 881 
    PROTECT(dna = getAttrib(a, R_DimNamesSymbol));
- 
 
 882 
 
- 
 
 883 
    if (LOGICAL(resize)[0] && dna != R_NilValue) {
- 
 
 884 
 
- 
 
 885 
	SEXP dnna, dnr, dnnr;
- 
 
 886 
 
- 
 
 887 
	PROTECT(dnna = getAttrib(dna, R_NamesSymbol));
- 
 
 888 
	PROTECT(dnnr = allocVector(STRSXP,n));
- 
 
 889 
	PROTECT(dnr  = allocVector(VECSXP,n));
- 
 
 890 
 
- 
 
 891 
	for (i=0; i<n; i++) {
- 
 
 892 
	    SET_VECTOR_ELT(dnr, i, VECTOR_ELT(dna, pp[i]));
- 
 
 893 
	    if (dnna != R_NilValue)
- 
 
 894 
		SET_STRING_ELT(dnnr, i, STRING_ELT(dnna, pp[i]));
- 
 
 895 
	}
- 
 
 896 
 
- 
 
 897 
	if (dnna != R_NilValue)
- 
 
 898 
	    setAttrib(dnr, R_NamesSymbol, dnnr);
- 
 
 899 
	setAttrib(r, R_DimNamesSymbol, dnr);
- 
 
 900 
	UNPROTECT(3); /* dnna, dnr, dnnr */
- 
 
 901 
    }
- 
 
 902 
    /* free temporary memory */
799 
    UNPROTECT(6);
 903 
    vmaxset(vmax);
- 
 
 904
- 
 
 905 
    UNPROTECT(6); /* dimsa, perm, r, dimsr, resize, dna */
800 
    return r;
 906 
    return r;
801 
}
 907 
}