#include "nutpier_density_kernel_v1.h"
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#include <ctype.h>

typedef struct { size_t n; double *x; double *y; } bound_t;
static int fail(char *out,size_t cap,const char *s){if(cap){size_t n=strlen(s);if(n>=cap)n=cap-1;memcpy(out,s,n);out[n]=0;}return 2;}
static void ws(const char **p,const char *e){while(*p<e && isspace((unsigned char)**p))++*p;}
static int string(const char **p,const char *e,const char **s,size_t *n){ws(p,e);if(*p>=e||*(*p)++!='"')return 0;*s=*p;while(*p<e&&**p!='"'){if(**p=='\\')return 0;++*p;}if(*p>=e)return 0;*n=(size_t)(*p-*s);++*p;return 1;}
static int number(const char **p,const char *e,double *v){char *q;ws(p,e);*v=strtod(*p,&q);if(q==*p||q>e||!isfinite(*v))return 0;*p=q;return 1;}
static int array(const char **p,const char *e,double **out,size_t *n){size_t cap=16,k=0;double *a;ws(p,e);if(*p>=e||*(*p)++!='[')return 0;a=(double*)malloc(cap*sizeof(double));if(!a)return 0;ws(p,e);if(*p<e&&**p==']'){++*p;*out=a;*n=0;return 1;}for(;;){double v;if(!number(p,e,&v)){free(a);return 0;}if(k==cap){cap*=2;double *z=(double*)realloc(a,cap*sizeof(double));if(!z){free(a);return 0;}a=z;}a[k++]=v;ws(p,e);if(*p>=e){free(a);return 0;}if(**p==']'){++*p;break;}if(*(*p)++!=','){free(a);return 0;}}*out=a;*n=k;return 1;}
static int skip_value(const char **p,const char *e){double v,*a=NULL;size_t n;if(number(p,e,&v))return 1;if(array(p,e,&a,&n)){free(a);return 1;}return 0;}
static int parse(const char *json,size_t len,bound_t *b){const char *p=json,*e=json+len;int seen=0;size_t nx=0,ny=0;ws(&p,e);if(p>=e||*p++!='{')return 0;for(;;){const char *s;size_t k;double v;ws(&p,e);if(p<e&&*p=='}'){++p;break;}if(!string(&p,e,&s,&k)){return 0;}ws(&p,e);if(p>=e||*p++!=':')return 0;if(k==1&&s[0]=='N'){if((seen&1)||!number(&p,e,&v)||v<1||v!=floor(v))return 0;b->n=(size_t)v;seen|=1;}else if(k==1&&s[0]=='x'){if((seen&2)||!array(&p,e,&b->x,&nx))return 0;seen|=2;}else if(k==1&&s[0]=='y'){if((seen&4)||!array(&p,e,&b->y,&ny))return 0;seen|=4;}else if(!skip_value(&p,e))return 0;ws(&p,e);if(p>=e)return 0;if(*p=='}'){++p;break;}if(*p++!=',')return 0;}ws(&p,e);return p==e&&seen==7&&nx==b->n&&ny==b->n;}
NUTPIER_DENSITY_KERNEL_EXPORT uint32_t nutpier_density_kernel_abi_version(void){return 1;}
NUTPIER_DENSITY_KERNEL_EXPORT int32_t nutpier_density_kernel_bind(const char *json,size_t len,size_t ndim,const char *layout,size_t layout_len,void **out,char *err,size_t cap){bound_t *b=(bound_t*)calloc(1,sizeof(*b));const char expected[]="alpha\nbeta\nsigma";*out=NULL;if(!b)return fail(err,cap,"allocation failed");if(!parse(json,len,b)){free(b->x);free(b->y);free(b);return fail(err,cap,"expected JSON N and numeric arrays x,y with matching lengths");}if(ndim!=3||layout_len!=sizeof(expected)-1||memcmp(layout,expected,sizeof(expected)-1)){free(b->x);free(b->y);free(b);return fail(err,cap,"expected unconstrained layout alpha\\nbeta\\nsigma");}*out=b;return 0;}
NUTPIER_DENSITY_KERNEL_EXPORT void nutpier_density_kernel_destroy(void *v){bound_t*b=(bound_t*)v;if(b){free(b->x);free(b->y);free(b);}}
NUTPIER_DENSITY_KERNEL_EXPORT int32_t nutpier_density_kernel_workspace(void*b,void**out,char*err,size_t cap){(void)b;(void)err;(void)cap;*out=NULL;return 0;}
NUTPIER_DENSITY_KERNEL_EXPORT void nutpier_density_kernel_workspace_destroy(void*b,void*w){(void)b;(void)w;}
NUTPIER_DENSITY_KERNEL_EXPORT int32_t nutpier_density_kernel_evaluate(void *v,void*w,const double*q,size_t ndim,double*lp,double*g,char*err,size_t cap){const bound_t*b=(const bound_t*)v;double a,be,r,s,s2,ss=0,sr=0,srx=0;(void)w;if(ndim!=3)return fail(err,cap,"evaluation dimension mismatch");a=q[0];be=q[1];r=q[2];s=exp(r);s2=s*s;if(!isfinite(s)||s==0)return 1;for(size_t i=0;i<b->n;i++){double z=b->y[i]-a-be*b->x[i];ss+=z*z;sr+=z;srx+=z*b->x[i];}/* Stan propto=TRUE,jacobian=TRUE: prior kernels; likelihood -N log(sigma) and quadratic; plus lower-bound Jacobian r. */*lp=-0.5*(a*a/25.0+be*be/4.0+s2+ss/s2)-(double)b->n*r+r;g[0]=-a/25.0+sr/s2;g[1]=-be/4.0+srx/s2;g[2]=-s2-(double)b->n+ss/s2+1.0;if(!isfinite(*lp)||!isfinite(g[0])||!isfinite(g[1])||!isfinite(g[2]))return 1;return 0;}