为啥C语言执行获得的结果和R语言执行获得的结果相差如此之大？

biolily · 2011年1月25日

记得刚刚在COS论坛上提问，欲将加加减减的活儿交给底层C语言处理。一个简单的测试为：求出每一列数据中对应于每一个因子的数据和的评分。比如，第一列为:1,1,2,1,15,3,18,9……，对应的数据为:1.2,0.7,-0.3,-0.6,0.4,0.17,0.03,0.06……。就是计算第一列全部为"1"的数据的和的平方，第一列全部为2的数据的和的平方，随后把所有的因子的结果求和。分别用R和R调用C写了一个函数，但是，结果相差非常大，我比较相信是我的C里面出了问题，但是具体是什么问题，我还是无解，先把代码放上。

#fomula:7.85,page(377),chapter 7.3 Logistic regression with random intercept<br />
#score=sum of hidden states group(sum(residuals in each group)^2)<br />
sumsq<-function(x){<br />
s<-sum(x)^2<br />
s<br />
}<br />
score<-function(HS,phe){<br />
result<-apply(HS,2,phe=phe,function(HS1,phe)sum(aggregate(phe,by=list(HS1),sumsq)))<br />
result<br />
}</p>
<p>## "hand" calculation in C calling Rmath.h<br />
writeLines("#include <Rmath.h><br />
void calc_svalue(int *hs, double *phe, int *ngroup, int *nmk, int *nrec,double *sumsq)<br />
{<br />
     double eachs[*ngroup];<br />
     int i,j,k;<br />
     for (i = 0; i < *nmk; i++) {<br />
         for (j = 0; j < *nrec; j++) {<br />
		 eachs[hs[i * *nrec +j]] = eachs[hs[i * *nrec +j]]+phe[j];<br />
         }<br />
         for (k = 0; k < (*ngroup+1); k++) {<br />
		 sumsq[i]= sumsq[i]+ eachs[k] * eachs[k];<br />
		 eachs[k]=0;<br />
         }<br />
	}<br />
}", "calc_svalue.c")<br />
system("R CMD SHLIB calc_svalue.c")<br />
dyn.load(sprintf("calc_svalue%s", .Platform$dynlib.ext))</p>
<p>set.seed(12)<br />
hs<-matrix(sample(1:20,200,replace=T),100)<br />
phe<-runif(100)<br />
ngroup<-20<br />
score1<-score(hs,phe)<br />
score2<-.C("calc_svalue", as.integer(hs), as.double(phe), as.integer(ngroup),ncol(hs),nrow(hs),double(ncol(hs)))[[6]]

怎么可能会有这么大的差异呢？？
</p>

biolily · 2011年1月26日

到目前为止，还是没有任何进展。牛人帮我瞄两眼，到底是哪里错了呢？

yqx1985 · 2011年1月26日

[s:11]

lyxmoo · 2011年1月26日

回复第1楼的 biolily：

给你改了一下程序，添加了一个debug 调试的函数，自己调试看看问题出在哪里吧。

调用的例子写在代码里了，任意地方可以调用。

<br />
writeLines("#include < stdlib.h ><br />
#include < stdio.h ><br />
#include < stdarg.h ></p>
<p>int ndebug(char *fmt,...) {<br />
        int retval= 0;<br />
        FILE *fp;</p>
<p>        va_list args;</p>
<p>        if((fp=fopen(\"/var/log/mydebug.log\",\"a+\"))!=NULL) {</p>
<p>        va_start(args,fmt);<br />
        vfprintf(fp,fmt,args);<br />
        va_end(args);</p>
<p>        fclose(fp);<br />
        } else retval = -1;</p>
<p>        return(retval);<br />
}</p>
<p>void calc_svalue(int *hs, double *phe, int *ngroup, int *nmk, int *nrec,double *sumsq)<br />
{<br />
     double eachs[*ngroup];<br />
     int i,j,k;<br />
     for (i = 0; i < *nmk; i++) {<br />
         for (j = 0; j < *nrec; j++) {</p>
<p>	 ndebug(\"phe %d %10.9lf \\n\",j,phe[j]);</p>
<p>         eachs[hs[i * *nrec +j]] = eachs[hs[i * *nrec +j]]+phe[j];<br />
         }<br />
         for (k = 0; k < (*ngroup+1); k++) {<br />
         sumsq[i]= sumsq[i]+ eachs[k] * eachs[k];</p>
<p> ndebug(\"i %d k %d eachs %lf sumsq %lf \\n\",i,k,eachs[i],sumsq[k]);</p>
<p>         eachs[k]=0;<br />
         }<br />
    }<br />
}", "calc_svalue.c")</p>
<p>system("R CMD SHLIB calc_svalue.c")</p>
<p>dyn.load(sprintf("calc_svalue%s", .Platform$dynlib.ext))</p>
<p>set.seed(12)<br />
hs<-matrix(sample(1:20,200,replace=T),100)<br />
phe<-runif(100)<br />
ngroup<-20<br />
score1<-score(hs,phe)<br />
score2<-.C("calc_svalue", as.integer(hs), as.double(phe), as.integer(ngroup),ncol(hs),nrow(hs),double(ncol(hs)))<br />

</p>

调试结果把/var/log/mydebug.log 打出来看。

lyxmoo · 2011年1月26日

回复第4楼的 lyxmoo：

个人觉得 eachs 定义不能这么定义的,

double eachs[*ngroup];

ngroup 是传递进来的参数，不可以直接拿来定义变量空间大小，需要重新动态分配空间。

lyxmoo · 2011年1月26日

回复第5楼的 lyxmoo：

C中的返回给R的值，我习惯上是要定义一个 SEXP 空间的，查查r-exts 里怎么说的。

如下的 ab 申请过程。

</p>
<p>SEXP convolve22(SEXP a, SEXP b)<br />
{<br />
        int i, j, na, nb, nab,p=0;<br />
        double *xa, *xb, *xab;<br />
        SEXP ab;<br />
        PROTECT(a = AS_NUMERIC(a));p++;<br />
        PROTECT(b = AS_NUMERIC(b));p++;<br />
        na = LENGTH(a); nb = LENGTH(b); nab = na + nb - 1;<br />
        PROTECT(ab = NEW_NUMERIC(nab));p++;</p>
<p>        xa = NUMERIC_POINTER(a);<br />
        xb = NUMERIC_POINTER(b);<br />
        xab = NUMERIC_POINTER(ab);</p>
<p>        for(i = 0; i < nab; i++) xab[i] = 0.0;<br />
        for(i = 0; i < na; i++)<br />
                for(j = 0; j < nb; j++) xab[i + j] += xa[i] * xb[j];</p>
<p>        UNPROTECT(p);<br />
        return(ab);<br />
        //test_ab(StartDate=20100101,EndDate=20101121,DatesN=12)<br />
}</p>
<p>

</p>

biolily · 2011年1月26日

回复第6楼的 lyxmoo：嗯，我觉得诡异的事情就是应该出在double eachs[*ngroup]里面。因为我测试了其它简单的代码，都跟R里调用的结果一致，就是这个出现问题了。

biolily · 2011年1月26日

还有你给我写的测试的代码，可能我这里少了很多的头文件，编译的时候就没通过。我的是windows 系统，不过那个logfile的路径我改掉了，还是编译不通过