use Array:all;
use Math:all;
use SimplePrint:all;
use String:{tochar};

#define APEXRAVEL(y) (_reshape_([prod(_shape_(y))],y))

inline int[+] indraxis0(int i, int[+]M)
{ /* M[scalar;;;;] */
 junk = print(tochar("psel i,M are:"));
 junk = junk +print(i);
 junk = junk +print(M);
 z = psel([i],M);
 junk = junk+ print(tochar("z=psel([i],M) is:"));
 junk = junk +print(z);
 
 return(z,junk);
}

inline int[+] indraxis0(int[+] i, int[+]M)
{ /* M[i;;;;] i may be matrix */
 cellshape=_drop_SxV_(1,_shape_(M));
 cell=genarray(cellshape,0);
 z = with(. <= iv <= .)
	genarray(_shape_(i),psel([i[iv]],M),cell);
 return(z);
}

inline int[+] indsaxis0(int i, int[+]M, int[*] Y)
{ /* M[scalar;;;;] <- Y */
 z = modarray(M,i,Y);
 return(z);
}

inline int[+] indrrank1(int[*] i, int[+]V)
{ /* Index reference V[any] */
 cellshape=_drop_SxV_(1,_shape_(V));
 cell=genarray(cellshape,0);
 z = with(0*_shape_(i) <= [iv] < _shape_(i))
	genarray(_shape_(i),psel([i[iv]],V),cell);

 return(z);
}

inline int[+] indsaxis0(int[*] i, int[+]M, int[+]Y)
{ /* M[i;;;;] <- Y     i may be matrix */
  /* Each scalar in i replaces one item (subarray) in M */
 ri = APEXRAVEL(i);
 for(k=0; k< _shape_(ri)[0]; k++){
	iv = [ri[k]];
	M = modarray(M,iv,psel(iv,Y));
 }
 return(M);
}

inline int[*] psel( int[.] idx, int[*] array)
{
 new_shape = _drop_SxV_( 
	_sel_( [0], _shape_(idx)), _shape_(array));
  res = with( . <= iv <= .) {
          new_idx = _cat_VxV_( idx, iv);

        } genarray( new_shape, _sel_(new_idx, array));

  return( res);

}


int main()
{
  shp=[14,5];
  v = with (. <= [iv] <= . ) 
	genarray([prod(shp)],20+iv);
  z =  print(tochar("v is:"));
  z= z+print(v);
  m2 = _reshape_(shp,v);
  z =  z+print(tochar("m2 is:"));
  z= z+ print(m2);
  i = _reshape_([3, 2], [1, 3, 2, 8, 12, 0]);
  vec = [1,3,2,4,9,12];
  new = with (. <= iv <= .)
	genarray(_shape_(vec),m2[vec[iv]]+1000);
  z =  z+print(tochar("new is:"));
  z=z+print(new);

  k = _reshape_([2,3,4], [2,1,0, 2,1,0, 2,1,0, 2,1,0, 
                        2,1,0, 2,1,0, 2,1,0, 2,1,0]);
  z =  z+print(tochar("k is:"));
  z = z+ print(k);

/* V1[v2] <- v3 */
 for(i=0; i < _shape_(vec)[0]; i++){
  	z =  z+print(tochar("for v,iv are:"));
	z=z+print(v);
	iv= [vec[i]];
	z=z+print(iv);
	v = modarray(v,iv,i);
 }
 z =  z+print(tochar("after for v,m2 are:"));
 z=z +print(v);
 z=z +print(m2);

/*  M1[vec;] <- M2  E.g., mdiv pivot */
 v1,junk = indraxis0(vec,m2)+10000;
 m3 = indsaxis0(vec,m2,v1);
 z =  z+junk+print(tochar("m3 is:"));
 z=z +print(m3);

/* M1[;vec] <- M2 */

m22 = _reshape_([2,2], [0,1,1,1]);
k2 = indrrank1([0,1,2,3],k);
m4,junk= indraxis0(m22,k2);
 z = z+junk+ print(tochar("m4 is:"));
z=z +print(m4);

return(z);

}