[Rcpp-devel] Rccp code with vector and matrix inputs and matrix output

Petre Caraiani petre.caraiani at gmail.com
Mon Mar 31 09:51:25 CEST 2014


Hello everybody,
Thanks for your suggestions. I used the cppFunction as follows:
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]

 NumericMatrix fun (NumericVector a, NumericVector b, NumericMatrix Am ,
int T  ){
int nrow = Am.nrow();
int ncol  = Am.ncol();

for (int ii = 0; ii <= nrow-1; ii++)
      for (int jj = 0; jj<= ncol-1; jj++)

         Am(ii,jj) = exp((jj+1)) ;
return Am;
}

A<-matrix(data=NA,nrow=10,ncol=13)
a <- rep(1:10)
b<-  rep(1:10)
T<-15

fun(a,b,A,T)



Which works fine, however when I write:
#include <Rcpp.h>
using namespace Rcpp;

// For more on using Rcpp click the Help button on the editor toolbar

// [[Rcpp::export]]

 NumericMatrix fun (NumericVector a, NumericVector b, NumericMatrix Am ,
int T  ){
int nrow = Am.nrow();
int ncol  = Am.ncol();

for (int ii = 0; ii <= nrow-1; ii++)
      for (int jj = 0; jj<= ncol-1; jj++)

         Am(ii,jj) = exp((jj+1)/(T-1)) ;
return Am;
}



and use  fun(a,b,A,T), I get a matrix of 1s. Is there something wrong with
the division in this Cpp code?

I get the same thing when I write: Am(ii,jj) = exp((jj+1)/(T-1)*log(a[ii]))
; but, again, the results are normal if I just write: Am(ii,jj) = log(a[ii])

Thanks!






On Fri, Mar 28, 2014 at 4:30 PM, Hadley Wickham <h.wickham at gmail.com> wrote:

> I'd recommend:
>
> a) using cppFunction instead of the older cxxfunction
>
> b) creating an example that anyone can copy and paste out of their
> email client and into R.
>
> You are more likely to get helpful responses if you reduce the burden
> on the potential helpers as much as possible
>
> Hadley
>
> On Fri, Mar 28, 2014 at 6:51 AM, Petre Caraiani
> <petre.caraiani at gmail.com> wrote:
> > Thank you for your quick reply!
> > I corrected the issue mentioned by you, but I get the same error.
> > It must be a beginner's issue.
> >
> > The data is quite big, but this kind of output could be produced with any
> > random data, I guess.
> >
> > Namely, using the Rcpp, I get Am as:
> >        [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
> [,13]
> >   [1,] 0.1825   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> > NA
> >   [2,] 0.2050   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> > NA
> >   [3,] 0.2083   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> > NA
> >   [4,] 0.2100   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> > NA
> >   [5,] 0.2133   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> > NA
> >   [6,] 0.2100   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> > NA
> >   [7,] 0.2025   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> > NA
> >   [8,] 0.2020   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> > NA
> >   [9,] 0.2025   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> > NA
> >  [10,] 0.1867   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> > NA
> >  [11,] 0.1867   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> > NA
> >  [12,] 0.1867   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> > NA
> >  [13,] 0.2000   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> > NA
> >  [14,]     NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> > NA
> >  and only NAs everywhere else
> >
> > while I get this using the pure R code:
> >   X0   X0.1   X0.2   X0.3   X0.4   X0.5   X0.6   X0.7   X0.8   X0.9
>  X0.10
> > X0.11  X0.12
> > 1  0.1811 0.1797 0.1784 0.1770 0.1757 0.1743 0.1730 0.1717 0.1704 0.1691
> > 0.1678 0.1665 0.1653
> > 2  0.2021 0.1993 0.1964 0.1937 0.1909 0.1883 0.1856 0.1830 0.1804 0.1779
> > 0.1753 0.1729 0.1704
> > 3  0.2044 0.2007 0.1969 0.1933 0.1897 0.1862 0.1828 0.1794 0.1760 0.1728
> > 0.1696 0.1664 0.1634
> > 4  0.2056 0.2013 0.1970 0.1929 0.1888 0.1849 0.1810 0.1772 0.1735 0.1698
> > 0.1662 0.1627 0.1593
> > 5  0.2085 0.2039 0.1993 0.1949 0.1905 0.1863 0.1821 0.1780 0.1741 0.1702
> > 0.1664 0.1626 0.1590
> > 6  0.2052 0.2006 0.1960 0.1915 0.1872 0.1829 0.1788 0.1747 0.1707 0.1669
> > 0.1631 0.1594 0.1557
> > 7  0.2004 0.1982 0.1961 0.1941 0.1920 0.1900 0.1880 0.1860 0.1840 0.1820
> > 0.1801 0.1782 0.1763
> > 8  0.1993 0.1967 0.1941 0.1916 0.1891 0.1866 0.1841 0.1817 0.1793 0.1769
> > 0.1746 0.1723 0.1701
> > 9  0.2005 0.1985 0.1966 0.1947 0.1928 0.1909 0.1890 0.1871 0.1853 0.1835
> > 0.1817 0.1799 0.1781
> > 10 0.1854 0.1840 0.1827 0.1814 0.1801 0.1788 0.1775 0.1763 0.1750 0.1737
> > 0.1725 0.1713 0.1700
> > 11 0.1855 0.1843 0.1831 0.1819 0.1807 0.1795 0.1784 0.1772 0.1761 0.1749
> > 0.1738 0.1727 0.1715
> > 12 0.1849 0.1831 0.1813 0.1795 0.1778 0.1760 0.1743 0.1726 0.1709 0.1692
> > 0.1676 0.1659 0.1643
> > 13 0.1964 0.1928 0.1893 0.1859 0.1826 0.1793 0.1760 0.1728 0.1697 0.1666
> > 0.1636 0.1607 0.1578
> > 14 0.1733 0.1716 0.1700 0.1683 0.1667 0.1651 0.1635 0.1619 0.1603 0.1587
> > 0.1572 0.1557 0.1542
> > 15 0.1814 0.1795 0.1777 0.1759 0.1741 0.1723 0.1705 0.1687 0.1670 0.1653
> > 0.1636 0.1619 0.1602
> > 16 0.1809 0.1784 0.1761 0.1737 0.1714 0.1691 0.1669 0.1646 0.1625 0.1603
> > 0.1582 0.1560 0.1540
> > 17 0.1828 0.1796 0.1765 0.1734 0.1704 0.1675 0.1646 0.1617 0.1589 0.1562
> > 0.1534 0.1508 0.1482
> > 18 0.1768 0.1744 0.1720 0.1696 0.1673 0.1651 0.1628 0.1606 0.1584 0.1563
> > 0.1541 0.1520 0.1500
> > 19 0.1537 0.1534 0.1531 0.1528 0.1525 0.1522 0.1519 0.1516 0.1513 0.1510
> > 0.1507 0.1504 0.1501
> > 20 0.1412 0.1425 0.1437 0.1450 0.1463 0.1476 0.1489 0.1502 0.1515 0.1529
> > 0.1542 0.1556 0.1570
> > 21 0.1460 0.1469 0.1479 0.1489 0.1498 0.1508 0.1518 0.1528 0.1538 0.1548
> > 0.1558 0.1569 0.1579
> > 22 0.1312 0.1323 0.1335 0.1347 0.1359 0.1372 0.1384 0.1396 0.1409 0.1421
> > 0.1434 0.1447 0.1460
> > 23 0.1555 0.1544 0.1532 0.1521 0.1509 0.1498 0.1487 0.1476 0.1465 0.1454
> > 0.1443 0.1432 0.1421
> > 24 0.1562 0.1556 0.1551 0.1546 0.1541 0.1535 0.1530 0.1525 0.1520 0.1515
> > 0.1509 0.1504 0.1499
> > 25 0.1561 0.1556 0.1551 0.1545 0.1540 0.1534 0.1529 0.1523 0.1518 0.1513
> > 0.1507 0.1502 0.1497
> > 26 0.1332 0.1338 0.1345 0.1351 0.1358 0.1365 0.1371 0.1378 0.1385 0.1392
> > 0.1399 0.1406 0.1413
> > 27 0.1548 0.1546 0.1544 0.1542 0.1540 0.1538 0.1536 0.1534 0.1532 0.1530
> > 0.1528 0.1526 0.1524
> > 28 0.1547 0.1543 0.1540 0.1537 0.1534 0.1530 0.1527 0.1524 0.1521 0.1517
> > 0.1514 0.1511 0.1508
> > 29 0.1552 0.1554 0.1556 0.1558 0.1560 0.1563 0.1565 0.1567 0.1569 0.1571
> > 0.1573 0.1575 0.1577
> > 30 0.1571 0.1575 0.1579 0.1583 0.1587 0.1591 0.1594 0.1598 0.1602 0.1606
> > 0.1610 0.1614 0.1618
> > 31 0.1571 0.1575 0.1580 0.1584 0.1588 0.1593 0.1597 0.1601 0.1605 0.1610
> > 0.1614 0.1618 0.1623
> > 32 0.1574 0.1581 0.1588 0.1594 0.1601 0.1608 0.1615 0.1622 0.1629 0.1636
> > 0.1644 0.1651 0.1658
> > 33 0.1580 0.1594 0.1607 0.1621 0.1635 0.1648 0.1662 0.1677 0.1691 0.1705
> > 0.1720 0.1734 0.1749
> >
> >
> > On Fri, Mar 28, 2014 at 1:41 PM, Petre Caraiani <
> petre.caraiani at gmail.com>
> > wrote:
> >>
> >> Thank you for your quick reply!
> >> I corrected the issue mentioned by you, but I get the same error.
> >> It must be a beginner's issue.
> >>
> >> The data is quite big, but this kind of output could be produced with
> any
> >> random data, I guess.
> >>
> >> Namely, using the Rcpp, I get Am as:
> >>        [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
> >> [,13]
> >>   [1,] 0.1825   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> >> NA
> >>   [2,] 0.2050   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> >> NA
> >>   [3,] 0.2083   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> >> NA
> >>   [4,] 0.2100   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> >> NA
> >>   [5,] 0.2133   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> >> NA
> >>   [6,] 0.2100   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> >> NA
> >>   [7,] 0.2025   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> >> NA
> >>   [8,] 0.2020   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> >> NA
> >>   [9,] 0.2025   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> >> NA
> >>  [10,] 0.1867   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> >> NA
> >>  [11,] 0.1867   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> >> NA
> >>  [12,] 0.1867   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> >> NA
> >>  [13,] 0.2000   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> >> NA
> >>  [14,]     NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
> >> NA
> >>  and only NAs everywhere else
> >>
> >> while I get this using the pure R code:
> >>   X0   X0.1   X0.2   X0.3   X0.4   X0.5   X0.6   X0.7   X0.8   X0.9
>  X0.10
> >> X0.11  X0.12
> >> 1  0.1811 0.1797 0.1784 0.1770 0.1757 0.1743 0.1730 0.1717 0.1704 0.1691
> >> 0.1678 0.1665 0.1653
> >> 2  0.2021 0.1993 0.1964 0.1937 0.1909 0.1883 0.1856 0.1830 0.1804 0.1779
> >> 0.1753 0.1729 0.1704
> >> 3  0.2044 0.2007 0.1969 0.1933 0.1897 0.1862 0.1828 0.1794 0.1760 0.1728
> >> 0.1696 0.1664 0.1634
> >> 4  0.2056 0.2013 0.1970 0.1929 0.1888 0.1849 0.1810 0.1772 0.1735 0.1698
> >> 0.1662 0.1627 0.1593
> >> 5  0.2085 0.2039 0.1993 0.1949 0.1905 0.1863 0.1821 0.1780 0.1741 0.1702
> >> 0.1664 0.1626 0.1590
> >> 6  0.2052 0.2006 0.1960 0.1915 0.1872 0.1829 0.1788 0.1747 0.1707 0.1669
> >> 0.1631 0.1594 0.1557
> >> 7  0.2004 0.1982 0.1961 0.1941 0.1920 0.1900 0.1880 0.1860 0.1840 0.1820
> >> 0.1801 0.1782 0.1763
> >> 8  0.1993 0.1967 0.1941 0.1916 0.1891 0.1866 0.1841 0.1817 0.1793 0.1769
> >> 0.1746 0.1723 0.1701
> >> 9  0.2005 0.1985 0.1966 0.1947 0.1928 0.1909 0.1890 0.1871 0.1853 0.1835
> >> 0.1817 0.1799 0.1781
> >> 10 0.1854 0.1840 0.1827 0.1814 0.1801 0.1788 0.1775 0.1763 0.1750 0.1737
> >> 0.1725 0.1713 0.1700
> >> 11 0.1855 0.1843 0.1831 0.1819 0.1807 0.1795 0.1784 0.1772 0.1761 0.1749
> >> 0.1738 0.1727 0.1715
> >> 12 0.1849 0.1831 0.1813 0.1795 0.1778 0.1760 0.1743 0.1726 0.1709 0.1692
> >> 0.1676 0.1659 0.1643
> >> 13 0.1964 0.1928 0.1893 0.1859 0.1826 0.1793 0.1760 0.1728 0.1697 0.1666
> >> 0.1636 0.1607 0.1578
> >> 14 0.1733 0.1716 0.1700 0.1683 0.1667 0.1651 0.1635 0.1619 0.1603 0.1587
> >> 0.1572 0.1557 0.1542
> >> 15 0.1814 0.1795 0.1777 0.1759 0.1741 0.1723 0.1705 0.1687 0.1670 0.1653
> >> 0.1636 0.1619 0.1602
> >> 16 0.1809 0.1784 0.1761 0.1737 0.1714 0.1691 0.1669 0.1646 0.1625 0.1603
> >> 0.1582 0.1560 0.1540
> >> 17 0.1828 0.1796 0.1765 0.1734 0.1704 0.1675 0.1646 0.1617 0.1589 0.1562
> >> 0.1534 0.1508 0.1482
> >> 18 0.1768 0.1744 0.1720 0.1696 0.1673 0.1651 0.1628 0.1606 0.1584 0.1563
> >> 0.1541 0.1520 0.1500
> >> 19 0.1537 0.1534 0.1531 0.1528 0.1525 0.1522 0.1519 0.1516 0.1513 0.1510
> >> 0.1507 0.1504 0.1501
> >> 20 0.1412 0.1425 0.1437 0.1450 0.1463 0.1476 0.1489 0.1502 0.1515 0.1529
> >> 0.1542 0.1556 0.1570
> >> 21 0.1460 0.1469 0.1479 0.1489 0.1498 0.1508 0.1518 0.1528 0.1538 0.1548
> >> 0.1558 0.1569 0.1579
> >> 22 0.1312 0.1323 0.1335 0.1347 0.1359 0.1372 0.1384 0.1396 0.1409 0.1421
> >> 0.1434 0.1447 0.1460
> >> 23 0.1555 0.1544 0.1532 0.1521 0.1509 0.1498 0.1487 0.1476 0.1465 0.1454
> >> 0.1443 0.1432 0.1421
> >> 24 0.1562 0.1556 0.1551 0.1546 0.1541 0.1535 0.1530 0.1525 0.1520 0.1515
> >> 0.1509 0.1504 0.1499
> >> 25 0.1561 0.1556 0.1551 0.1545 0.1540 0.1534 0.1529 0.1523 0.1518 0.1513
> >> 0.1507 0.1502 0.1497
> >> 26 0.1332 0.1338 0.1345 0.1351 0.1358 0.1365 0.1371 0.1378 0.1385 0.1392
> >> 0.1399 0.1406 0.1413
> >> 27 0.1548 0.1546 0.1544 0.1542 0.1540 0.1538 0.1536 0.1534 0.1532 0.1530
> >> 0.1528 0.1526 0.1524
> >> 28 0.1547 0.1543 0.1540 0.1537 0.1534 0.1530 0.1527 0.1524 0.1521 0.1517
> >> 0.1514 0.1511 0.1508
> >> 29 0.1552 0.1554 0.1556 0.1558 0.1560 0.1563 0.1565 0.1567 0.1569 0.1571
> >> 0.1573 0.1575 0.1577
> >> 30 0.1571 0.1575 0.1579 0.1583 0.1587 0.1591 0.1594 0.1598 0.1602 0.1606
> >> 0.1610 0.1614 0.1618
> >> 31 0.1571 0.1575 0.1580 0.1584 0.1588 0.1593 0.1597 0.1601 0.1605 0.1610
> >> 0.1614 0.1618 0.1623
> >> 32 0.1574 0.1581 0.1588 0.1594 0.1601 0.1608 0.1615 0.1622 0.1629 0.1636
> >> 0.1644 0.1651 0.1658
> >> 33 0.1580 0.1594 0.1607 0.1621 0.1635 0.1648 0.1662 0.1677 0.1691 0.1705
> >> 0.1720 0.1734 0.1749
> >>
> >>
> >> On Fri, Mar 28, 2014 at 1:25 PM, Dirk Eddelbuettel <edd at debian.org>
> wrote:
> >>>
> >>>
> >>> On 28 March 2014 at 10:30, Petre Caraiani wrote:
> >>> | Hello,
> >>> | The following code works well in R:
> >>> | attach(dataqtr)
> >>> |
> >>> | dataqtr <- data.table(dataqtr)
> >>> | setkeyv(dataqtr,c("gvkey","qtr"))
> >>> |
> >>> | vec_growth <- data.frame(0,0,0,0,0,0,0,0,0,0,0,0,0)
> >>> | vec_eps    <- data.frame(0,0,0,0,0,0,0,0,0,0,0,0,0)
> >>> |
> >>> | T <- 15
> >>> |
> >>> | for (i in 1:nrow(dataqtr)) {
> >>> |   vec_growth[i,] <- ( dataqtr[i,LTG] * exp(1:(T-2)/(T-1)*log( dataqtr
> >>> | [i,meanLTG] / dataqtr[i,LTG] )))
> >>> |
> >>> | }
> >>> |
> >>> |
> >>> | However I am not able to reproduce it using the following Rccp code:
> >>> |
> >>> | a<- dataqtr[,LTG]
> >>> | b<- dataqtr[,meanLTG]
> >>> |
> >>> | src <-'
> >>> |  Rcpp::NumericVector a(aa);
> >>> |  Rcpp::NumericVector b(bb);
> >>> |  Rcpp::NumericMatrix Am (A);
> >>> |  int n = a.size();
> >>> |  int m = b.size();
> >>> |  int T=15;
> >>> |
> >>> |  int nrows = Am.nrow();
> >>> |  int ncol  = Am.ncol();
> >>> |
> >>> |      for (int ii = 1; ii < nrows; ii++) {
> >>> |         for (int jj = 1; jj<ncol; jj++)  {
> >>> |         Am[ii,jj] =  a[jj] * exp(jj/(T-1)*log( b[jj] / a[jj] ));
> >>> |         }}
> >>> | return Am;
> >>> | '
> >>> | fun <- cxxfunction(signature(aa="numeric", bb="numeric",A="numeric"),
> >>> body =
> >>> | src, plugin="Rcpp")
> >>> |
> >>> | A<-matrix(data=NA,nrow=100,ncol=13)
> >>> | fun(a,b,A)
> >>> |
> >>> | I don't understand the error.
> >>>
> >>> What is the error you are getting and do not understand? Can you share
> >>> it?
> >>>
> >>> Is it __build-time__ ?  Is it __run-time__ ?
> >>>
> >>> At a first glance both your for loops are wrong as indices in C and C++
> >>> have
> >>> to go from 0 to n-1, not 1 to n.
> >>>
> >>> Dirk
> >>>
> >>> --
> >>> Dirk Eddelbuettel | edd at debian.org | http://dirk.eddelbuettel.com
> >>
> >>
> >
> >
> > _______________________________________________
> > Rcpp-devel mailing list
> > Rcpp-devel at lists.r-forge.r-project.org
> > https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
>
>
>
> --
> http://had.co.nz/
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20140331/903ecd67/attachment-0001.html>


More information about the Rcpp-devel mailing list