/* Frozen scalar feature implementation from commit 4945313.
 * Kept independent of patch_compute.h for numerical regression checks. */
/* Channel-major, PyTorch cross-correlation convention and exact GELU. */
static void reference_lp_conv(const float *x,float *y,int length,int ni,int no,int kernel,const float *w,const float *bias) {
    for(int o=0;o<no;o++) for(int t=0;t<length;t++) {
        float s=bias[o];
        for(int i=0;i<ni;i++) for(int k=0;k<kernel;k++) {
            int at=t+k-kernel/2;
            if(at>=0 && at<length) s+=x[i*length+at]*w[(o*ni+i)*kernel+k];
        }
        y[o*length+t]=.5f*s*(1.f+erff(s*.7071067811865475f));
    }
}
static void reference_lp_embed(const unsigned char *bytes,int n,int length,const float *w,float *out) {
    for(int i=0;i<length;i++) for(int c=0;c<LP_E;c++) out[c*length+i]=w[(i<n ? (int)bytes[i]+1:0)*LP_E+c];
}
static float reference_lp_dot(const float *x,const float *w,int n,float b) { for(int i=0;i<n;i++) b+=x[i]*w[i];return b; }
static void reference_lf_features(const PatchModel *m,const char *request,const unsigned char *bytes,int n,float *h,float *joint,float *scores,float *logs) {
    size_t size=strlen(request);
    if(!size || size>96 || n<1 || n>256) die("expected 1..96 command bytes and 1..256 binary bytes");
    float ce[LP_E*96],c1[LP_C*96],c2[LP_C*96],pooled[LP_POOL],c[LP_H];
    float be[LP_E*256],b1[LP_H*256],b2[LP_H*256];
    reference_lp_embed((const unsigned char*)request,(int)size,96,m->ce,ce);
    reference_lp_conv(ce,c1,96,LP_E,LP_C,5,m->c1w,m->c1b);reference_lp_conv(c1,c2,96,LP_C,LP_C,5,m->c2w,m->c2b);
    for(int k=0;k<LP_C;k++) for(int bin=0;bin<6;bin++) {
        float v=-INFINITY; for(int t=bin*16;t<(bin+1)*16;t++) if(c2[k*96+t]>v) v=c2[k*96+t]; pooled[k*6+bin]=v;
    }
    for(int k=0;k<LP_H;k++) c[k]=tanhf(reference_lp_dot(pooled,m->cfw+k*LP_POOL,LP_POOL,m->cfb[k]));
    reference_lp_embed(bytes,n,256,m->be,be);reference_lp_conv(be,b1,256,LP_E,LP_H,9,m->b1w,m->b1b);reference_lp_conv(b1,b2,256,LP_H,LP_H,9,m->b2w,m->b2b);
    for(int t=0;t<256;t++) {
        for(int k=0;k<LP_H;k++) {joint[k]=b2[k*256+t];joint[k+LP_H]=c[k];}
        for(int k=0;k<LP_H;k++) h[t*LP_H+k]=tanhf(reference_lp_dot(joint,m->fw+k*LP_J,LP_J,m->fb[k]));
        scores[t]=t<n ? reference_lp_dot(h+t*LP_H,m->pw,LP_H,m->pb[0]):-1e9f;
    }
    for(int k=0;k<LP_H;k++) { float v=-INFINITY;for(int t=0;t<256;t++) if(b2[k*256+t]>v) v=b2[k*256+t]; joint[k]=v; joint[k+LP_H]=c[k]; }
    scores[256]=reference_lp_dot(joint,m->nw,LP_J,m->nb[0]);
    float max=scores[0],sum=0;
    for(int t=0;t<257;t++) {if(!isfinite(scores[t])) die("nonfinite inference");if(scores[t]>max)max=scores[t];}
    for(int t=0;t<257;t++){logs[t]=scores[t]-max;scores[t]=expf(logs[t]);sum+=scores[t];}
    for(int t=0;t<257;t++){scores[t]/=sum;logs[t]-=logf(sum);}
}