I'm trying to calculate residuals after fitting linear regression, and I've got some code in C using the gsl which should do it. Everything works fine if I use static arrays (below, defining X[15], y[5] etc.). Trouble is, I won't know the number of individuals or covariates until runtime, so I'm stuck using dynamic arrays. But passing them to C causes a seg fault. I'm very sure I'm doing something (lots of things) stupid. If someone could point out what, I'd be very grateful. Thanks very much.

Andrew

Here's a toy D example:

import std.stdio;

extern(C) {
  void regress(int nInd, int nCov, ref double[] x, ref double[] y, ref double[] rOut);
}

void main(){
  int nInd = 5;
  int nCov = 3;
  double[] x = new double[nCov * nInd];
  x = [1, 4, 3, 1, 4, 3, 1, 4, 2, 1, 6, 7, 1, 3, 2];
  double[] y = new double[nInd];
  y = [5, 3, 4, 1, 5];
  double[] residuals = new double[nInd];
  regress(5, 3, x, y, residuals);

  writeln(residuals);
}

and the C code it calls:

#include <gsl/gsl_multifit.h>

void regress(int nInd, int nCov, double *x, double *y, double *rOut){
  int i, j;
  gsl_matrix *xMat, *cov;
  gsl_vector *yVec, *c, *r;
  double chisq;

  xMat = gsl_matrix_alloc(nInd, nCov);
  yVec = gsl_vector_alloc(nInd);
  r = gsl_vector_alloc(nInd);
  c = gsl_vector_alloc(nCov);
  cov = gsl_matrix_alloc(nCov, nCov);
  for(i = 0; i < nInd; i++)
    {
      gsl_vector_set(yVec, i, *(y+i));
      for(j = 0; j < nCov; j++)
        gsl_matrix_set(xMat, i, j, *(x + i * nCov + j));
    }

  gsl_multifit_linear_workspace *work =  gsl_multifit_linear_alloc(nInd, nCov);
  gsl_multifit_linear(xMat, yVec, c, cov, &chisq, work);
  gsl_multifit_linear_residuals(xMat, yVec, c, r);
  gsl_multifit_linear_free(work);


  for(i = 0; i < nInd; i++)
    rOut[i] = gsl_vector_get(r, i);
}

On 05/15/2014 01:55 PM, Andrew Brown wrote:
> extern(C) {
>    void regress(int nInd, int nCov, ref double[] x, ref double[] y, ref
> double[] rOut);
> }

I don't think that should even be allowed. C functions should not know or be compatible with 'ref' D parameters. Define the arguments as simple 'double *' or 'const double *'.

That makes sense because your C function is defined that way anyway.

> void main(){
>    int nInd = 5;
>    int nCov = 3;
>    double[] x = new double[nCov * nInd];
// ...
>    regress(5, 3, x, y, residuals);

You want to pass the address of the first array member: x.ptr (&(x[0]) would work as well). (Same for y.)

Otherwise, what ends up happening is that the address of the x and y slices are passed and C has no idea of what that is.

Ali

That worked a treat! Thank you very much!

On Thursday, 15 May 2014 at 21:11:54 UTC, Ali Çehreli wrote:
> On 05/15/2014 01:55 PM, Andrew Brown wrote:
> > extern(C) {
> >    void regress(int nInd, int nCov, ref double[] x, ref
> double[] y, ref
> > double[] rOut);
> > }
>
> I don't think that should even be allowed. C functions should not know or be compatible with 'ref' D parameters. Define the arguments as simple 'double *' or 'const double *'.
>
> That makes sense because your C function is defined that way anyway.
>
> > void main(){
> >    int nInd = 5;
> >    int nCov = 3;
> >    double[] x = new double[nCov * nInd];
> // ...
> >    regress(5, 3, x, y, residuals);
>
> You want to pass the address of the first array member: x.ptr (&(x[0]) would work as well). (Same for y.)
>
> Otherwise, what ends up happening is that the address of the x and y slices are passed and C has no idea of what that is.
>
> Ali

Ali Çehreli:

> I don't think that should even be allowed. C functions should not know or be compatible with 'ref' D parameters. Define the arguments as simple 'double *' or 'const double *'.

Time ago I suggested something like that, that is not meaningful to accept extern(C) function signatures with ref, [], fixed-sized arrays, lazy arguments, associative arrays, pure, and so on, because the C side doesn't know about those things or can't enforce them (like purity). Walter didn't agree with me, but I didn't understand his answer (or I don't remember it now).

Bye,
bearophile

For example, windows headers do use C++ &-references in function signatures and msdn provides code examples using that convention, the equivalent in D is ref.

On Friday, 16 May 2014 at 11:42:35 UTC, Kagamin wrote:
> For example, windows headers do use C++ &-references in function signatures and msdn provides code examples using that convention, the equivalent in D is ref.

But that's extern(C++), not extern(C)...

On Friday, 16 May 2014 at 14:52:17 UTC, Marc Schütz wrote:
> On Friday, 16 May 2014 at 11:42:35 UTC, Kagamin wrote:
>> For example, windows headers do use C++ &-references in function signatures and msdn provides code examples using that convention, the equivalent in D is ref.
>
> But that's extern(C++), not extern(C)...

I guess my confusion came about because in the page about interfacing with C, there's a static array example where parameters are given in terms D understands:

extern (C)
{
  void foo(ref int[3] a); // D prototype
}

I guess D has no problem translating that into a simple pointer that C can deal with. I assumed the same would be true of dynamic arrays, but maybe the leap is too far?

>
> I guess my confusion came about because in the page about interfacing with C, there's a static array example where parameters are given in terms D understands:
>
> extern (C)
> {
>   void foo(ref int[3] a); // D prototype
> }
>
> I guess D has no problem translating that into a simple pointer that C can deal with. I assumed the same would be true of dynamic arrays, but maybe the leap is too far?

And I've finally got round to seeing the table above, it clearly says that C array is equivalent to D pointer. Sorry for being a time waster.

On 05/16/2014 08:15 AM, Andrew Brown wrote:

> I guess my confusion came about because in the page about interfacing
> with C, there's a static array example where parameters are given in
> terms D understands:
>
> extern (C)
> {
>    void foo(ref int[3] a); // D prototype
> }
>
> I guess D has no problem translating that into a simple pointer that C
> can deal with. I assumed the same would be true of dynamic arrays, but
> maybe the leap is too far?

There is a major difference. A static array is direct equivalent of C arrays when it comes to how they are stored in memory. Static arrays are simply consecutive elements. (Of course, static arrays are superior to C arrays in many other aspects. :))

One difference between C arrays is the fact that static arrays are by-value when passed even to functions. (No more "decaying to pointer to first element" confusion.)

Since we know that references are implemented as pointers, 'ref int[3] a' is passed as a.ptr. Since the memory layout is the same as a C array, it works perfectly.

On the other hand, dynamic arrays (aka slices) are the equivalent of the following struct:

struct Slice(T)
{
    size_t length;
    T * ptr;    // points to an array of elements
}

Ali

On 05/16/2014 08:24 AM, Ali Çehreli wrote:

> On 05/16/2014 08:15 AM, Andrew Brown wrote:

>  >    void foo(ref int[3] a); // D prototype

> Since we know that references are implemented as pointers, 'ref int[3]
> a' is passed as a.ptr.

Sorry, that's confusing. Yes, it ends up being equal to a.ptr but conceptually, the compiler does not pass .ptr directly; it passes the address of the entire array. Since the address of the entire static array is the same as the address of its first element it equals a.ptr and works perfectly in the C land.

Ali