CORE_Test.hpp

#ifndef CORE_Test_HPP
#define CORE_Test_HPP
 
 
template<typename T,class I>
requires functions_type::isArithmeticType<T>
tBoolean CORE_Test::testArray(CORE_Array<T,I>& A) const {
 
 
    std::cout<<"\n\t testing "<<A.getIdentityString()<<"...\n";
    tBoolean ok=true,succeeds=true;
 
    //eps value
    tReal eps=std::numeric_limits<float>::epsilon();
 
    //test the size setting 
    tIndex n=10;
    A.setSize(n);
   
    if (A.getIdentityString().find("Ptr")!=tString::npos) {
        tMemSize memPtrA=sizeof(T*)+sizeof(tIndex)*2+8+sizeof(CORE_Object)+n*sizeof(T);//+8 instead of +1 because of memory alignement
        //Ptr Array
        if (A.getMemorySize()!=memPtrA) {
            std::cout<<"\n\t size of "<<A.getIdentityString()<<" "<<A.getMemorySize()<<": ("<<memPtrA<<")";
            return false;
        }
    } else {
        //val array
        tMemSize memValA=sizeof(std::valarray<T>)+sizeof(CORE_Object)+n*sizeof(T);
        if (A.getMemorySize()!=memValA) {
            std::cout<<"\n\t size of "<<A.getIdentityString()<<" "<<A.getMemorySize()<<": ("<<memValA<<")";
            return false;
        }
    }
    ok=(A.getSize()==10);
    if (!ok) {
        std::cout<<"\t\t\t memory size verification [FAILED] \n";
    } else {
        std::cout<<"\t\t\t memory size verification [OK] \n";
    }
    succeeds=succeeds&& ok;
    
    //test list initializer
    A={1};
    ok=(A.getSize()==1);
    ok=ok && (Equals(A[0],1));
    if (!ok) {
        std::cout<<" A:"<<A<<" ";
        std::cout<<"\t\t\t initialize to {1} [FAILED] \n";
    } else {
        std::cout<<"\t\t\t initialization to {1}  [OK] \n";
    }
    succeeds=succeeds&& ok;
    
    A={1,2,3,4};
    ok=(A.getSize()==4);
    for(auto i=0; const auto&  v: A) {
        ok = ok && Equals(v,i+1);
        i++;
    }
    if (!ok) {
        std::cout<<" A:"<<A<<" ";
        std::cout<<"\t\t\t initialize to {1,2,3,4} [FAILED] \n";
    } else {
        std::cout<<"\t\t\t initialization to {1,2,3,4}  [OK] \n";
    }
    succeeds=succeeds&& ok;
    
 
 
    //accessor tests
    try {
        A.get(4);
        ok=false;
    } catch(std::exception& e) {
        ok=true;
    }
    if (!ok) {
        std::cout<<" A:"<<A<<" ";
        std::cout<<"\t\t\t out bounds access detected [FAILED] \n";
    } else {
        std::cout<<"\t\t\t out bounds access detected [OK] \n";
    }
    succeeds=succeeds&& ok;
 
    //constant initialization & sum tests
    A=2.;
    n=A.getSize();
    T s=0;
    A.sum(s);
    ok=Equals(s,2*n);
    if (!ok) {
        std::cout<<" A:"<<A<<" ";
        std::cout<<"\t\t\t testing sum() & initialize(const T& v)  [FAILED] \n";
    } else {
        std::cout<<"\t\t\t testing sum() & initialize(const T& v) [OK] \n";
    }
    succeeds=succeeds&& ok;
 
    //prod test
    A.prod(s);
    ok=Equals(s,pow(2,n));
    if (!ok) {
        std::cout<<" A:"<<A<<" ";
        std::cout<<"\t\t\t testing prod()  [FAILED] \n";
    } else {
        std::cout<<"\t\t\t testing prod() [OK] \n";
    }
    succeeds=succeeds&& ok;
 
    //+= operator
    A+=1;
    for (tIndex i=0;i<n;i++) ok=ok && (Equals(A[i],3));
    if (!ok) {
        std::cout<<" A:"<<A<<" ";
        std::cout<<"\t\t\t testing +=(const T& v)  [FAILED] \n";
    } else {
        std::cout<<"\t\t\t testing +=(const T& v) [OK] \n";
    }
    succeeds=succeeds&& ok;
 
    //test uniform randomize & copy
    CORE_UniquePointer<CORE_Array<T,I>> pB=I::New();
    if (pB.get()==null) {
        std::cout<<"New Instanced has failed";
        return false;
    }
    CORE_Array<T,I> &B=*pB.get();
    B.setSize(10);
    B.uniformRandomize();
    A=B;
    ok=(A.getSize()==B.getSize());
    auto b=B.cbegin();
    std::for_each(A.cbegin(),A.cend(),
                  [&](const auto & a) {
                      ok= ok && Equals(a,*b);//validate the copy
                      ok =ok && (a>=0) && (a<=1);//validate the unifoem
                      b++;
                  });
 
    if (!ok) {
        std::cout<<" A:"<<A<<" B="<<B<<"  ";
        std::cout<<"\t\t\t testing =(const CORE_Array& v) [FAILED] \n";
    } else {
        std::cout<<"\t\t\t testing =(const CORE_Array& v) [OK] \n";
    }
    succeeds=succeeds&& ok;
 
    //testing Linf distance
    tIndex imax=0;
    tReal d=A.linfDistance(B,imax);
    ok=ok && (fabs(d)<eps);
    if (!ok) {
        std::cout<<" A:"<<A<<" B="<<B<<"  ";
        std::cout<<" linfMax:"<<d<<" imax:"<<imax<<"\n";
        std::cout<<"\t\t\t testing linfDistance [FAILED] \n";
    } else {
        std::cout<<"\t\t\t testing linfDistance [OK] \n";
    }
    succeeds=succeeds&& ok;
    
    //testing L2 distance
    d=A.l2Distance(B);
    ok=ok && (fabs(d)<eps);
    if (!ok) {
        std::cout<<" A:"<<A<<" B="<<B<<"  ";
        std::cout<<" l2:"<<d<<"\n";
        std::cout<<"\t\t\t testing l2Distance [FAILED] \n";
    } else {
        std::cout<<"\t\t\t testing l2Distance [OK] \n";
    }
    succeeds=succeeds&& ok;
 
 
    //testing swap
    
    A.swap(1,2);
    ok = Equals(A[1],B[2]);
    ok= ok && Equals(A[2],B[1]);
    if (!ok) {
        std::cout<<" A:"<<A<<" B="<<B<<"  ";
        std::cout<<" 1<->2 swaping \n";
        std::cout<<"\t\t\t testing swaping element [FAILED] \n";
    } else {
        std::cout<<"\t\t\t testing swaping element [OK] \n";
    }
    succeeds=succeeds&& ok;
    //transform test
    CORE_IntegerValArray C(4);
    C={0,1,2,3};
    auto lambda=[](const tInteger& x) {return x*x;};
    C.transform(lambda);
    ok=(C.getSize()==4);
    tInteger i=0;
    for(const auto Ci : C.getArray()) {
        ok = ok && (Equals(Ci,i*i));
        i++;
    }
    if (!ok) {
        std::cout<<" C:"<<C<<"\n";
        std::cout<<"\t\t\t testing   transform(x^2) [FAILED] \n";
    } else {
        std::cout<<"\t\t\t testing  transform(x^2) [OK] \n";
    }
    succeeds=succeeds && ok;
 
    
    //uniform randomize
    succeeds = uniformRandomTest<T>(A,B) && succeeds;
    
    //copy an integer array to real array
    A.copy(C);
    ok=(A.getSize()==4);
    for (i=0;i<4;i++) {
        ok = ok && (Equals(A[i],i*i));
    }
    if (!ok) {
        std::cout<<" A:"<<A<<" C:"<<C<<"\n";
        std::cout<<"\t\t\t testing  A.copy(C) [FAILED] \n";
    } else {
        std::cout<<"\t\t\t testing  A.copy(C)  [OK] \n";
    }
    succeeds=ok & succeeds;
 
    C=2;
    A=C;
    ok=(A.getSize()==4);
    for (i=0;i<4;i++) {
        ok = ok && (Equals(A[i],2));
    }
    if (!ok) {
        std::cout<<" A:"<<A<<" C:"<<C<<"\n";
        std::cout<<"\t\t\t testing  A=C [FAILED] \n";
    } else {
        std::cout<<"\t\t\t testing  A=C  [OK] \n";
    }
    succeeds=ok & succeeds;
 
 
    
    //performance test:
    
    CORE_Chrono timer;
    
    tIndex N=256*256*256*3;
    //N=16*16*16*3;//for valgrind
    tReal alpha=1,beta=2;
    std::cout<<"\t\t\t testing Y:=beta.Y+alpha X for a size N="<<N<<" randomizing... "<<std::flush;
 
    //A array initialization
    A.setSize(N);
    A.uniformRandomize();
    
    //B array initialization
    B.setSize(N);
    B.uniformRandomize();
    
    std::cout<<" computing Y... "<<std::flush;
    timer.start();
    A.axpy(alpha,B,beta);
    tULLInt duration=timer.stop();
    std::cout<<" time:"<<duration<<" ms \n";
 
    
    return succeeds;
 
}
 
template<typename T,class I>
tBoolean CORE_Test::uniformRandomTest(CORE_Array<T,I>& A,CORE_Array<T,I>& B) const {
     
    //uniform randomize
      
    //area of a disk defined by (x^2+y^2)=a^2  : Pi.a^2
    //area of the smallest square containing the previous circle : 4a^2
    //for arbitrary point of the square has a proprobility P= Pi.a^2/4.a^2=Pi/4.
    //so PI= 4 . P
 
    CORE_Chrono timer;
 
    tIndex p,N=128*128*128;
    //N=10;//for debug
    //N=16*16*16;//for valgrind
    A.setSize(2*N);
    
    timer.start();
    A.uniformRandomize();
  
    const T* vA=A.getValues();
    const T* eA=vA+2*N;
    T d=0,dk=0;
    p=0;
    while (vA!=eA) {
        d=0;
        
        dk=(*vA);
        dk*=dk;
        d+=dk;
        vA++;
 
        dk=(*vA);
        dk*=dk;
        d+=dk;
        vA++;
        
        p+=(d<=1);
        
    }
  
    //the probality is p=s/n
    //so pi ~4. p
    tReal Pi=4.*((tReal)p)/((tReal)N);
    //duration in ms
    tULLInt duration=timer.stop();
    tReal err=fabs(1-(Pi/M_PI));
    tReal errMax=1.e-3;
    std::cout<<"\t\t\t Pi:"<<std::setprecision(12)<<Pi<<" M_PI:"<<M_PI<<" err="<<err<<" (errMax="<<errMax<<") p="<<p<<" N="<<N<<" in "<<duration<<"ms ";
    tBoolean ok=(err<errMax);
    if (!ok) {
        std::cout<<" [FAILED] \n";
    } else {
        std::cout<<" [OK] \n";
    }
    return ok;
}
 
 
 
template<typename T,typename K,K D,class S,class I>
tBoolean  CORE_Test::testField(CORE_Field<T,K,D,S,I>& A) const {
    tBoolean ok=true,succeeds=true;
 
    std::cout<<"\n\t testing "<<A.getIdentityString()<<"...\n";
    tDimension k;
    
    //init tests
    A={1,2,3,4,5,6,7,8,9,10};
    ok=(A.getElementsNumber()==3) &&  (A.getSize()==9);
    
    tInteger i=1;
    for (const auto &Ai : A) {
        ok = ok && ( Equals(Ai,i));
        i++;
    }
    if (!ok) {
        std::cout<<"error in copy initializer \n";
        std::cout<<"A size : "<<A.getElementsNumber()<<" array size:"<<A.getSize()<<"\n";
        std::cout<<"A="<<A<<"\n";
        std::cout<<"\t\t\t testing A={1,2,3,4,5,6,7,8,9} [FAILED] \n";
    } else {
        std::cout<<"\t\t\t testing A={1,2,3,4,5,6,7,8,9}  [OK] \n";
    }
    
    succeeds=succeeds&& ok;
    
    //test the setting normalized values
    std::array<tReal,3> V={1,1,1};
    CORE_RealPtrArray::Normalize(V);
    A.initialize(V);
    ok=(A.getElementsNumber()==3);
    for (i=0;i<3;i++) {
        for (k=0;k<3;k++) {
            ok=ok && Equals(A(i,k),1./sqrt(3.));
        }
    }
    if (!ok) {
        std::cout<<A<<"\n";
        std::cout<<" error in setting value to normalized of ("<<V[0]<<","<<V[1]<<","<<V[2]<<")\n";
        std::cout<<"\t\t\t testing  initialize(std::array<D>) [FAILED] \n";
    } else {
        std::cout<<"\t\t\t testing initialize(std::array<D>)  [OK] \n";
    }
    succeeds=succeeds&& ok;
    
    //nan test
    ok=!(A.isNANContained());
    A(1,0)=sqrt(-1);
    ok=A.isNANContained() && ok;
    if (!ok) {
        std::cout<<" A:"<<A<<" ";
        std::cout<<"\t\t\t testing nan detection  [FAILED] \n";
    } else {
        std::cout<<"\t\t\t testing nan detection [OK] \n";
    }
  
    succeeds=succeeds&& ok;
 
    //+= operator
    A(1,0)=0;
    A+=1;
    ok=ok && Equals(A(1,0),1);
    if (!ok) {
        std::cout<<" A:"<<A<<" ";
        std::cout<<"\t\t\t testing +=1 [FAILED] \n";
    } else {
        std::cout<<"\t\t\t testing +=1 [OK] \n";
    }
    succeeds=succeeds&& ok;
 
 
    
    //random test  
    //volume defined by (x^2+y^2)=a^2  : (4/3).Pi.a^3
    //volume of the smallest square containing the previous sphere : 8a^3
    //for arbitrary point of the square has a proprobility P= (4/3).Pi.a^3/(8.a^3)=Pi/6.
    //so PI= 6 . P
 
    
    //tIndex N=10;
    tIndex N=10000000;
    //N=16*16*16;//for valgrind
    CORE_Chrono timer;
    A.setElementsNumber(N);
    ok =(A.getSize()==N*D);
    if (!ok) {
        std::cout<<"error in setting size/elements of A \n";
        return false;
    }
    CORE_RealValArray nD;
 
    timer.start();
    A.uniformRandomize(0,1);
    //std::cout<<"A:"<<A<<"\n";
    A.mod2(nD);
    //std::cout<<"norm2:"<<nD<<"\n";
    //X:=(X^2+Y^2+Z^2)<=1)
    nD.transform([](const tReal& r) { return (r<=1);});
    //std::cout<<"r<=1:"<<nD<<"\n";
    T p=1;
    nD.sum(p);
    //the probality is p=s/n
     //so pi ~6. p
    tReal Pi=6.*((tReal)p)/((tReal)N);
    //duration in ms
    tULLInt duration=timer.stop();
    tReal err=fabs(1-(Pi/M_PI));
    tReal errMax=1.e-3;
    std::cout<<"\t\t\t Pi:"<<Pi<<" err="<<err<<" (errrMax="<<errMax<<") with  PI="<<M_PI<<" with "<<N<<" random  values computed in "<<duration<<"ms ";
    ok=(err<errMax);
    if (!ok) {
        std::cout<<" [FAILED] \n";
    } else {
        std::cout<<" [OK] \n";
    }
    succeeds=succeeds&& ok;
 
 
    //normalized 
    A.setSize(4);
    A.uniformRandomize(-1,1);
    A.normalize();
    ok=(A.getElementsNumber()==1);
    if (!ok) {
        std::cout<<" error in getElementsNumber() & setSize() \n";
        return false;
    }
    tReal r;
    ok=true;
    for (i=0;i<A.getElementsNumber();i++) {
        r=0;
        for (k=0;k<D;k++) {
            r+=A(i,k)*A(i,k);
        }
        r=sqrt(r);
        ok=(fabs(r-1)<1.e-12);
        if (!ok) break;
    }
    if (!ok) {
        std::cout<<"\n error in normalization of A at index "<<i<<" A=("<<A(i,0)<<","<<A(i,1)<<","<<A(i,2)<<") norm:"<<r<<"\n";
        std::cout<<"\t\t\t testing field normalization [FAILED] \n";
    } else {
        std::cout<<"\t\t\t testing field normalization [OK] \n";
    }
    succeeds=succeeds&& ok;
 
    //verify the ApplyDArray method
    //transform X as a circle of radius 1 to a radius 2 and center (1,1,1)
    //F(X)=alpha X + C
    A.setElementsNumber(4);
    A.uniformRandomize(-1,1);
    A.normalize();
    ok=(A.getSize()==4*D);
    if (!ok) {
        std::cout<<"error in setting setElementsNumber() & getSize() \n";
        return false;
    }
    std::array<tReal,D> C= {1,1,1};
    r=2;
    auto similitude=[&r,&C](auto& Xi) {
                        auto iCid=C.cbegin();//iterator to C(0)
                        std::for_each_n(&Xi,D,//iterator on d in [0,3[
                                        [&r,&iCid](auto &Xid){//X(i,d)
                                            Xid=Xid*r+(*iCid);
                                            iCid++;//&C(i,d)
                                        });
                        return Xi;
                    };
    A.elementsTransform(similitude);
    ok=true;
    for (i=0;i<A.getElementsNumber();i++) {
        r=0;
        for (k=0;k<D;k++) {
            r+=(A(i,k)-C[k])*(A(i,k)-C[k]);
        }
        r=sqrt(r);
        ok=(fabs(r-2)<1.e-12);
        if (!ok) break;
    }
    
    if (!ok) {
        std::cout<<"\n error in similitude of A at index "<<i<<" A=("<<A(i,0)<<","<<A(i,1)<<","<<A(i,2)<<") |A-C(1,1,1)|:"<<r<<"\n";
        std::cout<<"\t\t\t testing similitude transformation [FAILED] \n";
    } else {
        std::cout<<"\t\t\t testing similitude transformation [OK] \n";
    }
     
    succeeds=succeeds&& ok;
 
 
    //scalar product test
    //A={1,2,3,4,5,6};
    A.copy({1,2,3,4,5,6});
    //test uniform randomize & copy
    CORE_UniquePointer<CORE_Field<T,K,D,S,I>> pB=I::New();
    if (pB.get()==null) {
        std::cout<<"New Instanced has failed";
        return false;
    }
    CORE_Field<T,K,D,S,I> &B=*dynamic_cast<CORE_Field<tReal,K,D,S,I>*>(pB.get());
    B={1,1,1,1,1,1};
    T s;
    A.scalarProduct(B,s);
    ok=Equals(s,21);    
    if (!ok) {
        std::cout<<"\n error <A,B>="<<s<<" A="<<A<<" B="<<B<<"\n";
        std::cout<<"\t\t\t testing scalar product [FAILED] \n";
    } else {
        std::cout<<"\t\t\t testing scalar product [OK] \n";
    }
    succeeds=succeeds&& ok;
    
    
    B.mod2(nD);
    ok=(nD.getSize()==2);
    ok=ok && (nD[0]==3);
    ok=ok && (nD[1]==3);
    if (!ok) {
       std::cout<<"\n error |B|="<<nD<<" B="<<B<<" \n";
        std::cout<<"\t\t\t testing norm2 [FAILED] \n";
    } else {
        std::cout<<"\t\t\t testing norm2 [OK] \n";
    }
    succeeds=succeeds&& ok;
 
    //performance test:
      
    //performance test:
    N=256*256*256;
    //N=16*16*16;//for valgrind
    tReal alpha=1,beta=2;
    std::cout<<"\t\t\t testing Y:=beta.Y+alpha X for a size N="<<N<<" ";
    
    //A field
    A.setElementsNumber(N);
    A.uniformRandomize(0,1);
    B.setElementsNumber(N);
    B.uniformRandomize(0,1);
    
    timer.start();
    A.axpy(alpha,B,beta);
    duration=timer.stop();
    std::cout<<" time:"<<duration<<" ms \n";
    
    
    return succeeds;
}
#endif