/
vfromsp.c
192 lines (182 loc) · 7.96 KB
/
vfromsp.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
/* Copyright 1990-2011, Jsoftware Inc. All rights reserved. */
/* License in license.txt. */
/* */
/* Verbs: { on sparse arguments */
#include "j.h"
static A jtfromis1(J jt,A ind,A w,A z,I wf){A a,a1,j1,p,q,x,x1,y,y1;C*xu,*xuu,*xv;I an,*av,c,d,
h,i,*iv,j,*jv,k,m,n,*pv,*qu,*qv,r,s,*u,*v,xk,*yu,*yv;P*wp,*zp;
zp=PAV(z); wp=PAV(w); r=AR(ind); n=AN(ind);
a=SPA(wp,a); an=AN(a); av=AV(a); DO(an, if(wf==av[i]){h=i; break;});
y=SPA(wp,i); RZ(q=eps(fromr(sc(h),y),ravel(ind))); RZ(y=repeat(q,y)); RZ(x=repeat(q,SPA(wp,x)));
GA(a1,INT,r+an-1,1,0); v=AV(a1); SPB(zp,a,a1);
k=av[h]; u=av; DO(h, *v++=*u++;); DO(r, *v++=k++;); u++; DO(an-1-h, *v++=*u+++r-1;);
if(!r)
if(AR(z)){GA(q,INT,an-1,1,0); v=AV(q); DO(an, if(i!=h)*v++=i;); SPB(zp,i,fromr(q,y)); SPB(zp,x,x); R z;}
else R reshape(mtv,AN(x)?x:SPA(zp,e));
if(h){q=grade1(fromr(sc(h),y)); RZ(y=ifrom(q,y)); RZ(x=ifrom(q,x));}
RZ(q=odom(2L,r,AS(ind))); iv=AV(q);
m=*AS(y); s=0; j=-1; u=h+AV(y); v=u+an;
GA(p,INT,m,1,0); pv=AV(p); memset(pv,CFF,SZI*m);
GA(q,INT,m,1,0); qu=AV(q);
GA(q,INT,m,1,0); qv=AV(q);
DO(m-1, if(*u!=*v){pv[s]=*u; qu[s]=1+j; qv[s++]=i-j; j=i;} u=v; v+=an;);
if(m){i=m-1; pv[s]=*u; qu[s]=1+j; qv[s++]=i-j;}
RZ(j1=indexof(p,ind)); jv=AV(j1);
c=0; DO(n, if(s>jv[i])c+=qv[jv[i]];); i=aii(x); xk=i*bp(AT(x)); d=AN(a1);
GA(y1,INT, c*d,2, 0 ); v=AS(y1); v[0]=c; v[1]=d; yv= AV(y1); yu= AV(y);
GA(x1,AT(x),c*i,AR(x),AS(x)); *AS(x1)=c; xv=CAV(x1); xu=CAV(x);
for(i=0;i<n;++i){
k=jv[i];
if(s>k){
c=qu[k]; d=qv[k]; xuu=xu+c*xk; u=yu+c*an;
for(j=0;j<d;++j){
MC(xv,xuu,xk); xv+=xk; xuu+=xk;
DO(h, *yv++=*u++;); DO(r, *yv++=iv[i];); u++; DO(an-1-h, *yv++=*u++;);
}}
iv+=r;
}
if(h){q=grade1(y1); RZ(y1=ifrom(q,y1)); RZ(x1=ifrom(q,x1));}
SPB(zp,i,y1); SPB(zp,x,x1);
R z;
} /* ind{"r w along a sparse axis */
F2(jtfromis){A ind,x,z;B*b;I acr,af,an,ar,*av,k,m,*v,wcr,wf,wn,wr,*ws,wt;P*wp,*zp;
RZ(a&&w);
ar=AR(a); acr=jt->rank?jt->rank[0]:ar; af=ar-acr;
wr=AR(w); wcr=jt->rank?jt->rank[1]:wr; wf=wr-wcr; jt->rank=0;
if(af)R rank2ex(a,w,0L,acr,wcr,jtfromis);
wn=AN(w); ws=AS(w); wt=AT(w);
RZ(ind=pind(wcr?*(ws+wf):1,a));
GA(z,wt,1L,ar+wr-(0<wcr),ws); v=AS(z); ICPY(v+wf,AS(a),ar); if(wcr)ICPY(v+wf+ar,1+wf+ws,wcr-1);
zp=PAV(z); wp=PAV(w); SPB(zp,e,ca(SPA(wp,e)));
RZ(a=ca(SPA(wp,a))); av=AV(a); an=AN(a);
RZ(b=bfi(wr,a,1));
if(b[wf])R fromis1(ind,w,z,wf);
m=wcr; DO(wcr, m-=b[wf+i];); RZ(x=irs2(ind,SPA(wp,x),0L,ar,m,jtifrom));
if(k=ar-1)DO(an, if(av[i]>=wf)av[i]+=k;);
if(AR(z)){SPB(zp,a,a); SPB(zp,x,x); SPB(zp,i,ca(SPA(wp,i))); R z;}
else R AN(x)?reshape(mtv,x):ca(SPA(wp,e));
} /* a{"r w for numeric a and sparse w */
static A jtaaxis(J jt,A w,I wf,A a,I r,I h,I*pp,I*qq,I*rr){A q;B*b,*c,*d;I wr,x,y,z,zr;
wr=AR(w); zr=wr+r-h;
RZ(b=bfi(wr,a,1));
GA(q,B01,zr,1,0); c=BAV(q);
x=y=z=0; d=b; DO(wf, if(*d++)++x;); DO(h, if(*d++)++y;); DO(wr-wf-h, if(*d++)++z;);
*pp=x; *qq=y; *rr=z;
MC(c,b,wf); memset(c+wf,y?C1:C0,r); MC(c+wf+r,b+wf+h,wr-wf-h);
R ifb(zr,c);
}
A jtfrombsn(J jt,A ind,A w,I wf){A a,j1,p,q,x,x1,y,y1,ys,z;C*xu,*xuu,*xv;
I an,c,d,h,i,*iv,j,*jv,k,m,n,pp,*pv,qq,*qv,r,rr,s,*u,*v,wr,*ws,xk,*yu,*yv;P*wp,*zp;
RZ(ind&&w); /* need to handle 1==n */
wr=AR(w); ws=AS(w);
r=AR(ind)-1; v=AS(ind); h=v[r]; n=AN(ind)/h;
RZ(q=odom(2L,r,v)); iv=AV(q)-r;
GA(z,AT(w),1,wr+r-h,0); u=AS(z); ICPY(u,ws,wf); ICPY(u+wf,v,r); ICPY(u+wf+r,ws+wf+h,wr-wf-h);
zp=PAV(z); wp=PAV(w);
SPB(zp,e,ca(SPA(wp,e))); x=SPA(wp,x); y=SPA(wp,i);
a=SPA(wp,a); an=AN(a); SPB(zp,a,aaxis(w,wf,a,r,h,&pp,&qq,&rr));
if( !qq){SPB(zp,i,ca(y)); SPB(zp,x,frombu(ind,x,AR(x)-(wr-wf-rr))); R z;}
if(h>qq){q=nub(over(a,apv(h,wf,1L))); R frombsn(ind,reaxis(grade2(q,q),w),wf);}
if(1<r)RZ(ind=reshape(v2(n,h),ind));
RZ(ys=fromr(indexof(a,apv(h,wf,1L)),y));
RZ(q=eps(ys,ind));
if(!all1(q)){RZ(ys=repeat(q,ys)); RZ(y=repeat(q,y)); RZ(x=repeat(q,x));}
if(wf){q=grade1(ys); RZ(ys=ifrom(q,ys)); RZ(y=ifrom(q,y)); RZ(x=ifrom(q,x));}
m=*AS(y);
GA(p,INT,m,1,0); pv=AV(p);
GA(q,INT,m,1,0); qv=AV(q);
s=0; j=-1; u=AV(ys); v=u+h;
DO(m-1, if(ICMP(u,v,h)){pv[s]=1+j; qv[s++]=i-j; j=i;} u=v; v+=h;); if(m){pv[s]=1+j; qv[s++]=m-1-j;}
RZ(j1=indexof(ifrom(vec(INT,s,pv),ys),ind)); jv=AV(j1);
c=0; DO(n, if(s>jv[i])c+=qv[jv[i]];);
i=aii(x); j=AN(SPA(zp,a)); xk=i*bp(AT(x));
GA(y1,INT, c*j,2, 0 ); v=AS(y1); v[0]=c; v[1]=j; yv= AV(y1); yu= AV(y);
GA(x1,AT(x),c*i,AR(x),AS(x)); *AS(x1)=c; xv=CAV(x1); xu=CAV(x);
for(i=0;i<n;++i){
k=jv[i]; iv+=r;
if(s>k){
c=pv[k]; d=qv[k]; xuu=xu+c*xk; u=yu+c*an;
for(j=0;j<d;++j){
MC(xv,xuu,xk); xv+=xk; xuu+=xk;
DO(pp, *yv++=*u++;); DO(r, *yv++=iv[i];); u+=qq; DO(rr, *yv++=*u++;);
}}}
if(wf){q=grade1(y1); RZ(y1=ifrom(q,y1)); RZ(x1=ifrom(q,x1));}
SPB(zp,i,y1); SPB(zp,x,x1);
R z;
} /* (<"1 ind){w, sparse w and integer array ind */
static A jtfrombs1(J jt,A ind,A w,I wf){A*iv,x,y,z;I id,j,m,n,old,wr,wcr;
RZ(ind&&w);
if(!(BOX&AT(ind))){ASSERT(!AN(ind)||NUMERIC&AT(ind),EVINDEX); RZ(ind=every(ind,0L,jtright1));}
n=AN(ind); iv=AAV(ind); id=(I)ind*ARELATIVE(ind); wr=AR(w); wcr=wr-wf;
ASSERT(1>=AR(ind),EVRANK);
ASSERT(n<=wr-wf,EVLENGTH);
j=n; DO(n, --j; x=AADR(id,iv[j]); if(BOX&AT(x)&&!AR(x)&&(y=AAV0(x),!AN(y)&&1==AR(y)))--n; else break;);
z=w; old=jt->tbase+jt->ttop;
for(j=0;j<n;++j){
x=AADR(id,iv[j]);
if(BOX&AT(x)){
ASSERT(!AR(x),EVINDEX);
x=AAV0(x); m=*(wf+j+AS(w));
if(!AN(x))continue;
RZ(x=less(IX(m),pind(m,x)));
}
RZ(z=irs2(x,z,0L,RMAX,wcr-j,jtfromis)); gc(z,old);
}
R z;
} /* (<ind){"r w, sparse w */
F2(jtfrombs){A ind;I acr,af,ar,wcr,wf,wr;
RZ(a&&w);
ar=AR(a); acr=jt->rank?jt->rank[0]:ar; af=ar-acr;
wr=AR(w); wcr=jt->rank?jt->rank[1]:wr; wf=wr-wcr; jt->rank=0;
ASSERT(!af,EVNONCE);
if(ar){RE(aindex(a,w,wf,&ind)); ASSERT(ind,EVNONCE); R frombsn(ind,w,wf);}
else R frombs1(AAV0(a),w,wf);
} /* a{"r w for boxed a and sparse w */
F2(jtfromsd){A e,x,z;I acr,af,ar,*v,wcr,wf,wr,*ws;P*ap,*zp;
RZ(a&&w);
ar=AR(a); acr=jt->rank?jt->rank[0]:ar; af=ar-acr;
wr=AR(w); wcr=jt->rank?jt->rank[1]:wr; wf=wr-wcr; jt->rank=0;
if(af)R sprank2(a,w,0L,acr,wcr,jtfrom);
ASSERT(AT(w)&B01+INT+FL+CMPX,EVNONCE);
ap=PAV(a); ws=AS(w);
GA(z,STYPE(AT(w)),1L,ar+wr-(0<wcr),ws); zp=PAV(z);
v=AS(z); ICPY(v+wf,AS(a),ar); if(wcr)ICPY(v+wf+ar,1+wf+ws,wcr-1);
RZ(x=irs2(SPA(ap,e),w,0L,0L,wcr,jtifrom)); RZ(e=reshape(mtv,x));
ASSERT(all1(eq(e,x)),EVSPARSE);
SPB(zp,e,e);
SPB(zp,a,wf?plus(sc(wf),SPA(ap,a)):SPA(ap,a));
SPB(zp,i,SPA(ap,i));
if(wf){
RZ(x=irs2(SPA(ap,x),w,0L,RMAX,wcr,jtifrom));
RZ(x=cant2(less(IX(AR(x)),sc(wf)),x));
SPB(zp,x,x);
}else SPB(zp,x,ifrom(SPA(ap,x),w));
R z;
} /* a{"r w, sparse a, dense w */
F2(jtfromss){A e,x,y,z;B*b;I acr,af,ar,c,d,k,m,n,p,*u,*v,wcr,wf,wr,*ws,*yv;P*ap,*wp,*xp,*zp;
RZ(a&&w);
ar=AR(a); acr=jt->rank?jt->rank[0]:ar; af=ar-acr;
wr=AR(w); wcr=jt->rank?jt->rank[1]:wr; wf=wr-wcr; jt->rank=0;
if(af)R sprank2(a,w,0L,acr,wcr,jtfrom);
ASSERT(DTYPE(AT(w))&B01+INT+FL+CMPX,EVNONCE);
ap=PAV(a); wp=PAV(w); ws=AS(w);
GA(z,AT(w),1L,ar+wr-(0<wcr),ws); zp=PAV(z);
v=AS(z); ICPY(v+wf,AS(a),ar); if(wcr)ICPY(v+wf+ar,1+wf+ws,wcr-1);
RZ(x=irs2(SPA(ap,e),w,0L,0L,wcr,jtfrom)); RZ(e=reshape(mtv,x));
ASSERT(all1(denseit(eq(e,x))),EVSPARSE);
SPB(zp,e,e);
x=SPA(ap,a); if(ar>AN(x)){RZ(a=reaxis(IX(ar),a)); ap=PAV(a);}
x=SPA(wp,a); n=AN(x); RZ(b=bfi(wr,x,1));
if(wcr&&!b[wf]){b[wf]=1; ++n; RZ(w=reaxis(ifb(wr,b),w)); wp=PAV(w);}
GA(x,INT,ar+n-!!wcr,1,0); v=AV(x);
DO(wf, if(b[i])*v++=i;); DO(ar, *v++=wf+i;); DO(wcr-1, if(b[i+wf+1])*v++=wf+ar+i;);
SPB(zp,a,x);
RZ(x=irs2(SPA(ap,x),w,0L,RMAX,wcr,jtfrom)); xp=PAV(x);
y=SPA(xp,i); u=AV(y); c=*(1+AS(y)); m=*AS(y); k=0; DO(wf, if(b[i])++k;);
y=SPA(ap,i); v=AV(y); d=*(1+AS(y)); n=c+d-1; p=c-(1+k);
GA(y,INT,m*n,2,0); *AS(y)=m; *(1+AS(y))=n; yv=AV(y);
DO(m, if(k)ICPY(yv,u,k); ICPY(yv+k,v+d*u[k],d); if(p)ICPY(yv+k+d,u+1+k,p); yv+=n; u+=c;);
SPB(zp,i,y);
SPB(zp,x,SPA(xp,x));
R z;
} /* a{"r w, sparse a, sparse w */