int i;
real [int] tab(10), tab1(10); // 2 array of 10 real 
complex [int] ctab(10), ctab1(10); // 2 array of 10 complex 
//string [int] stab(10);   pas code ... 


//  real [int] tab2; //  bug 
tab = 1;  //  set all the array to 1
tab[1]=2;
ctab = 1+2i;  //  set all the array to 1+2i
ctab[1]=2;
cout << tab[1] << " " << tab[9] << " size of tab = " << tab.n << endl;
cout << ctab[1] << " " << ctab[9] << " size of ctab = " << ctab.n << endl;

tab1=tab;
tab=tab+tab1;
tab=2*tab+tab1*5;
tab1=2*tab-tab1*5;
tab+=tab;
cout << tab << endl;
cout << tab[1] << " " << tab[9] <<  endl;

ctab1=ctab;
ctab=ctab+ctab1;
ctab=2*ctab+ctab1*5;
ctab1=2*ctab-ctab1*5;
ctab+=ctab;
cout << ctab << endl;
cout << ctab[1] << " " << ctab[9] <<  endl;


real [string] map; //  a dynamique array
 string[string] smap;
for (i=0;i<10;i=i+1)
  {
    tab[i] = i*i;
    cout << i << " " << tab[i] << "\n";
  };

map["1"]=2.0;
map[2]=3.0; //  2 is automaticaly cast to the string "2"

cout << " map[\"1\"] = " << map["1"] << " == 2.0 ; "<< endl;
cout << " map[2] = " << map[2] << " == 3.0 "<< endl;
assert( abs(map["1"]-2.0)<1.e-6);
assert(abs(map[2]-3.0)<1e-6);



real [int] tab2=[1,2,3,3.14];
int  [int] itab2=[1,2,3,5];

cout << tab2 << endl;
cout << itab2 << endl;
tab2.resize(10);
for (int i=4;i<10;i++) tab2[i]=i;
cout << "tab2 = " << tab2 << endl;
tab2 /= 2; // bug  before v2.0-3
cout << "tab2 = " << tab2 << endl;
tab2 *= 2;
cout << "tab2 = " << tab2 << endl;


real [int,int] mat(10,10),mmat(10,10);
mat=0;
for(int i=0;i<mat.n;i++)
  for(int j=0;j<mat.m;j++)
    mat(i,j)=i+100*(j+1);
mmat=mat;
cout << mmat << endl;
mat.resize(15,15);
for(int i=10;i<mat.n;i++)
  for(int j=0;j<mat.m;j++)
    mat(i,j)=i+100*(j+1);
for(int i=0;i<mat.n;i++)
  for(int j=10;j<mat.m;j++)
    mat(i,j)=i+100*(j+1);
cout << mat << endl; // wrong
//mat(0:9,0:9)=mmat;  todo 
//cout << mat << endl; // wrong

// array of mesh 
mesh[int] aTh(10);
aTh[1]= square(2,2);
plot(aTh[1]);

// add 1 sep 2005 FH 
real[int]  a(5),b(5),c(5),d(5);
a = 1;
b = 2;
c = 3;
a[2]=0;
d = ( a ? b : c ); // i = 0, n-1  d[i] = a[i] ? b[i] : c[i] , 
cout << d << endl;
cout << "==   2       2       3       2       2 \n";
d = (a ? 1 : 10);
cout << " (a ? 1 : 10) "<< d << endl;
d = (a ? b : -1);
cout << " (a ? b : -1 ) "<< d << endl;
d = (a ? -2 : c);
cout << " (a ? -2 : c) " <<  d << endl;
d =   1./d; 
cout << " 1/ d == (a ? -2 : c) " <<  d << endl;
d =   a/3.; 
cout << "a/3 ==  " <<  d << endl;

// Liste all array operator and method 
// ok in version 2.0-3 
// they exist in the 3 type of array
// real[int] , long[int], complex[int] 
// ---------------------------------
int N=5;
{
real[int] a(N),b(N),c(N);
a =1; 
 a(0:4:2) = 2;
a(3:4) = 4;
cout <<" a = " << a << endl;
b = a+ a;
cout <<" b = a+a : " << b << endl;
b += a;
cout <<" b += a : " << b << endl;
b += 2*a;
cout <<" b += 2*a : " << b << endl;
b /= 2;
cout <<" b /= 2 : " << b << endl;
b .*= a; // same b = b .* a 
cout << "b*=a; b =" << b << endl;
b ./= a; // same b = b ./ a 
cout << "b/=a; b =" << b << endl;
c = a+b;
cout << " c =a+b : c=" << c << endl;
c = 2*a+4*b;
cout << " c =2*a+4b : c= " << c << endl;

c = a+4*b;
cout << " c =a+4b : c= " << c << endl;
c = -a+4*b;
cout << " c =-a+4b : c= " << c << endl;
c = -a-4*b;
cout << " c =-a-4b : c= " << c << endl;
c = -a-b;
cout << " c =-a-b : c= " << c << endl;

c = a .* b; 
cout << " c =a.*b  : c= " << c << endl;
c = a ./ b; 
cout << " c =a./b  : c= " << c << endl;
c = 2 * b;
cout << " c =2*b   : c= " << c << endl;
c =  b*2 ;
cout << " c =b*2   : c= " << c << endl;

/*  this operator do not exist 
c =  b/2 ;
cout << " c =b/2   : c= " << c << endl;
*/

// ---- the  methods --
cout << " ||a||_1     = " <<  a.l1     << endl;
cout << " ||a||_2     = " <<  a.l2     << endl;
cout << " ||a||_infty = " <<  a.linfty << endl;
cout << " sum a_i     = " <<  a.sum    << endl;
cout << " max a_i     = " <<  a.max  << " a[  " << a.imax << " ] = " << a[a.imax]     << endl;
cout << " imax a_i     = " <<  a.imax  << " a[  " << a.imax << " ] = " << a[a.imax]     << endl;
cout << " min a_i     = " <<  a.min  << " a[  " << a.imin << " ] = " << a[a.imin]    << endl;
cout << " imin a_i     = " <<  a.imin  << " a[  " << a.imin << " ] = " << a[a.imin]    << endl;
cout << " a'*a        = " <<  (a'*a)   << endl;
cout << " a quantile 0.2 =  " <<  a.quantile(0.2) << endl;

// array mapping  after version 2.3
int[int] I=[2,3,4,-1,3];
b=c=-3; 
b= a(I); // for( i=0;i<b.n;i++) if(I[i] >=0)  b[i]=a[I[i]];
c(I)= a; // for( i=0;i<I.n;i++) if(I[i] >=0)  C(I[i])=a[i];
cout << " b = a(I) : " << b << "\n  c(I) = a " << c << endl;
c(I) += a;// for( i=0;i<I.n;i++) if(I[i] >=0)  C(I[i])+=a[i];
cout << " b = a(I) : " << b << "\n  c(I) = a " << c << endl;   


}
{
// bidimensionnal  array
  int N=3,M=4;

  real[int,int] A(N,M),RA(N,M);
  real[int]  b(N),c(M);
  b=[1,2,3];
  c=[4,5,6,7];

  complex[int,int]  C(N,M);
  complex[int]  cb=[1,2,3],cc=[10i,20i,30i,40i];


  b=[1,2,3];

  int [int] I=[2,0,1];
  int [int] J=[2,0,1,3];

  A=1; // set the all matrix
  A(2,:) = 4; //  the full line 2
  A(:,1) = 5; //  the full column 1
  A(0:N-1,2) = 2; // set the column 2
  A(1,0:2) = 3; // set the line 1 from 0 to 2

  cout << " A = " << A << endl;
  // outer product 
  C  =  cb*cc';
  C +=  3*cb*cc';
  C -=  5i*cb*cc';
  cout << " C = " << C << endl;
  // the way to transform a array to a sparce matrix
  matrix B;
  B = A; 
  B=A(I,J); // B(i,j)= A(I(i),J(j))
  B=A(I^-1,J^-1);  // B(I(i),J(j))= A(i,j)  

  A = 2.*b*c'; // outer product  
  cout << " A = " << A << endl;
  cout << " A(1,2) " << A(1,2)  << endl;
  B = b*c'; // outer product  B(i,j)  = b(i)*c(j) 
  B = b*c'; // outer product  B(i,j)  = b(i)*c(j) 
  B = (2*b*c')(I,J); //   outer product  B(i,j)  = b(I(i))*c(J(j)) 
  B = (3.*b*c')(I^-1,J^-1); // outer product  B(I(i),J(j))  = b(i)*c(j) 
  cout << "B = (3.*b*c')(I^-1,J^-1) =  " << B << endl;  
  cout << " b =" <<  b << endl;
  b = exp(b) ;
  cout << " exp(b) =" <<  b << endl;
  cb += complex(10.)*cc(0:2);
  cout << " cb =" <<  cb << endl;
  cb = exp(cb) ;
  cout << " exp(cb) =" <<  cb << endl;
  cout << " exp(cb).re =" <<  cb.re << endl;
  cout << " exp(cb).im =" <<  cb.im << endl;
  cb.im = 0.;
  cout << cb << endl;
  b += cb.re + cb.im;   //  do not work to do
  cout << " b = " << endl; 
  {ofstream FA("A.txt");
  FA << A << endl;
  }
  {ifstream FA("A.txt");
  FA >> RA ;
  cout << RA << endl; 
  RA -= A;
  cout << " RA =  00 ??? " << RA.linfty  << endl; 
  assert( RA.linfty < 1e-12);
  }

}
{
// sort array : 
real[int] a=[3,5,7,9,0];
real[int] b(a);
int[int] p=[0,1,2,3,4];

b=a;
cout << " a =" <<a << endl;
sort(b,p);
cout << " b.sort=" <<  b << endl;
cout << " b = " << b << endl;
cout << " p =  " << p << endl;
b=a;
b(0:5:2).sort;
cout << b << endl; 
cout << b(0:5:2)  << endl; 
cout << " quantile(0.2)  = " <<  b.quantile(0.2) << endl; 
p[0] = 1000;
p[1] = 888;
cout << "before    p(1:3).sort  " << p  << endl; 
p(1:3).sort;
cout << "after   p(1:3).sort " << p  << endl; 
p.sort; // version 3.19 ##
cout << "after   p.sort " << p  << endl; 

macro AA [ [1,2],[3,4]] //
macro BB [ [1,2],[5,10]] //
 cout << AA[1][0] << " 3  " << endl;
 cout << trace( AA * BB' + AA - BB' ) <<endl;
}
//  version 3.2  mai 2009
//  like math lab. and scilab
{
int[int] tt(2:10); //  2,3,4,5,6,7,8,9,10  
int[int] t1(2:3:10); // 2,5,8,
cout << " tt(2:10)= " << tt << endl;
cout << " t1(2:3:10)= " << t1 << endl;
tt=1:2:5;
cout << " 1.:2:5 =>  " << tt << endl;
}

{
real[int] tt(2:10); //  2,3,4,5,6,7,8,9,10  
real[int] t1(2.:3:10.); // 2,5,8,
cout << " tt(2:10) = = " << tt << endl;
cout << " t1(2.:3:10.)= " << t1 << endl;
tt=1.:0.5:3.999;
cout << " 1.:0.5:3.999 =>  " << tt << endl;
}
{
complex[int] tt(2.+0i:10.+0i); //  2,3,4,5,6,7,8,9,10  
complex[int] t1(2.:3.:10.); // 2,5,8,
cout << " tt(2.+0i:10.+0i)= " << tt << endl;
cout << " t1(2.:3.:10.);= " << t1 << endl;
cout << " tt.re real part array   " << tt.re << endl ; 
 //  the real part array of the complex array \index{re}\index{array!re}
cout << " tt.im imag part array   " << tt.im << endl ; 
//  the imag part array of the complex array \index{im}\index{array!im}


}
{

real [int] tab(10), tab1(10); // 2 array of 10 real
//real [int] tab2;    //  bug array with no size
tab = 1.03;                //  set all the array to 1.03
tab[1]=2.15;
cout << tab[1] << " " << tab[9] << " size of tab = "
     << tab.n << " min: " << tab.min << "  max:" << tab.max
     << " sum : "   << tab.sum <<   endl; //
tab.resize(12); //  change the size of array tab
  // to 12 with preserving first value
tab(10:11)=3.14; //  set unset value
cout <<" resize tab: " <<  tab << endl;
real [string] tt;
tt["+"]=1.5;
cout<<tt["a"]<<"  "<<tt["+"]<<endl;
real[int]  a(5),b(5),c(5),d(5);
a = 1;
b = 2;
c = 3;
a[2]=0;
d = ( a ? b : c ); // for i = 0, n-1  : d[i] = a[i] ? b[i] : c[i] ,
cout << " d = ( a ? b : c )  is " << d << endl;
d = ( a ? 1 : c );// for i = 0, n-1: d[i] = a[i] ? 1 : c[i] ,   (v2.23-1)
d = ( a ? b : 0 );// for i = 0, n-1: d[i] = a[i] ? b[i] : 0 ,   (v2.23-1)
d = ( a ? 1 : 0 );// for i = 0, n-1: d[i] = a[i] ? 0 : 1 ,     \hfill(v2.23-1)
 tab.sort ; //  sort the array tab  (version 2.18) 
cout << " tab (after sort) "  << tab << endl;
int[int] ii(0:d.n-1); // set array ii to 0,1, ..., d.n-1
d=-1:-5; // set d to  -1,-2, .. -5
sort(d,ii); // sort array d and ii in parallele
cout << " d " << d << "\n ii = " << ii << endl;

}

//  version 3.8-1
for(int i=0;i<ARGV.n;++i)
  {
    cout << ARGV[i] << endl;
  }
// Array of FE.
{
mesh Th=square(2,2);
fespace Vh(Th,P1);// scalar 
fespace Wh(Th,[P1,P1]);// vector
Vh[int] vh(10);
Wh[int] [uh,wh](10); 
vh[2]=x;
[uh[2],wh[2]]=[y,x];
[uh[2],wh[2]]=[y,x]; // set 3rd value 
cout << wh.n << endl;
cout << vh.n <<endl;
vh.resize(20);
cout << vh.n <<endl;
}



{ // version 3.43-3  init of map ...
  cout << " new stuff " << endl; 
  int[string] a=["2",1,"12",4];
  string[string] b= ["-1",1,"13","qsdqdq"];
  cout << " a " << a  << endl;
   cout << " b " << b  << endl; 
// the new kind of loop: 
  for [i,v:a] cout << " a: " << i << " " << v << endl;

  real[int,int] t=[[1,2],[2,3],[3,4],[5,6]]; 
  cout << t << endl;
  cout << t(1:3,0:1)' << endl; 

}