///
  /// @primary authors: S.Gorbunov; I.Kisel
  /// @        authors: H.Pabst et al. (Intel); M.Zyzak; I.Kulakov
  ///

#include <arbb.hpp>
using namespace arbb;
#include "fit.h"



#if defined(CONVENTIONAL)
  #include "fit_interface_conventional.h"
#elif defined(SQUARE_ROOT1) || defined(SQUARE_ROOT1S) || defined(SQUARE_ROOT2)
  #include "fit_interface_squareroot.h"
#elif defined(UD_FILTERING)
  #include "fit_interface_ud_filtering.h"
#endif


#pragma warning(disable:4305)

//! global in
TracksCt<FTYPE> gTracksCt;
StationsCt<FTYPE> gStationsCt;

//! global out
FitResults<FTYPE> gResults;

template <typename U, typename Ptr>
void copy_in(dense<U, 1>& dst, Ptr src, unsigned int size)
{
  dense<U> temp; bind(temp, src, size);
  dst = temp;
}

template <typename U, typename Ptr>
void copy_in(dense<U, 2>& dst, Ptr src, unsigned int size1, unsigned int size2)
{
  dense<U, 2> temp; bind(temp, src, size2, size1);//, size2);
  dst = temp;
}

template <typename U, typename Ptr>
void copy_out(dense<U, 1> src, Ptr dst, unsigned int size)
{
  if (!fit_use_range) {
    dense<U, 1> temp; bind(temp, dst, size);
    temp = src;
  } else {
    // TODO: clean up the code when 'range' is fixed
    const_range<U> r = src.read_only_range();
    for (typename const_range<U>::const_iterator I = r.begin(); I != r.end(); I++) {
//      dst[I - r.begin()] = value(*I);
      dst[I - r.begin()] = *I;
    }
  }
}

template<typename T>
HitInfoCt<T>::HitInfoCt( const HitInfo &st, int nStations )
{
    //! global copy in
    copy_in(cos_phi,  st.cos_phi,  nStations);
    copy_in(sin_phi,  st.sin_phi,  nStations);
    copy_in(sigma2,   st.sigma2,   nStations);
    copy_in(sigma216, st.sigma216, nStations);
}
    
template<typename T>
HitXYInfoCt<T>::HitXYInfoCt( const HitXYInfo &st, int nStations )
{
    //! global copy in
  copy_in(C00, st.C00, nStations);
  copy_in(C10, st.C10, nStations);
  copy_in(C11, st.C11, nStations);
}

template<typename T>
StationsCt<T>::StationsCt( const Stations &st ):
  UInfo(st.UInfo, st.nStations),     VInfo(st.VInfo, st.nStations),    XYInfo(st.XYInfo, st.nStations)
{
    //! global copy in
    copy_in(z, st.z, st.nStations);
    copy_in(thick, st.thick, st.nStations);
    copy_in(zhit, st.zhit, st.nStations);
    copy_in(RL, st.RL, st.nStations);
    copy_in(RadThick, st.RadThick, st.nStations);
    copy_in(logRadThick, st.logRadThick, st.nStations);
    copy_in(SyF, st.SyF, st.nStations);
    copy_in(SyL, st.SyL, st.nStations);
    copy_in(mapX, st.mapX[0], st.nStations, MaxNFieldCoeff);
    copy_in(mapY, st.mapY[0], st.nStations, MaxNFieldCoeff);
    copy_in(mapZ, st.mapZ[0], st.nStations, MaxNFieldCoeff);
}

template<typename T>
TracksCt<T>::TracksCt( Tracks& t )
{
    //! global copy in
    copy_in(MC_x, t.MC_x, t.nTracks);
    copy_in(MC_y, t.MC_y, t.nTracks);
    copy_in(MC_z, t.MC_z, t.nTracks);
    copy_in(MC_px, t.MC_px, t.nTracks);
    copy_in(MC_py, t.MC_py, t.nTracks);
    copy_in(MC_pz, t.MC_py, t.nTracks);
    copy_in(MC_q, t.MC_q, t.nTracks);
    copy_in(nHits, t.nHits, t.nTracks);
    copy_in(NDF, t.NDF, t.nTracks);
    copy_in(hitsX, t.hitsX[0], t.nTracks, t.nStations);
    copy_in(hitsY, t.hitsY[0], t.nTracks, t.nStations);
    copy_in(hitsX2, t.hitsX2, t.nStations, t.nTracks);
    copy_in(hitsY2, t.hitsY2, t.nStations, t.nTracks);
    copy_in(hitsIsta, t.hitsIsta[0], t.nTracks, t.nStations);
}


void fitTracksCtEntry(i32 vSize)
{
    /* to alleviate global vector overhead, we currently use their local copies 
       during the computation */
    FitResults<FTYPE> results(vSize);
    StationsCt<FTYPE> vStationsCt = gStationsCt;
    TracksCt<FTYPE> vTracksCt = gTracksCt;
    FieldRegionCt<FTYPE> magFieldCt(vSize);
    
#if defined(CONVENTIONAL)
    FitC< FTYPE, FTYPE > fitter;
    FitC< FTYPE, dense<FTYPE> > fitter2;
#elif defined(SQUARE_ROOT1) || defined(SQUARE_ROOT1S) || defined(SQUARE_ROOT2)
    FitSR< FTYPE, FTYPE > fitter;
    FitSR< FTYPE, dense<FTYPE> > fitter2;
#elif defined(UD_FILTERING)
    FitUD< FTYPE, FTYPE > fitter;
    FitUD< FTYPE, dense<FTYPE> > fitter2;
#endif
    
    usize j;
    _for( j = 0, j < Ntimes, j ++ ){
        fitter.Fit( vTracksCt, vStationsCt, magFieldCt, results );
    }_end_for;

    fitter2.Extrapolate( results, vTracksCt.MC_z, results.VecT[4], magFieldCt, fill<FTYPE>(0,results.length()) );

    /* write back the values of local copies to global vectors */
    gResults = results;
}

void fitTracksCt( ftype (*T)[6], ftype (*C)[15] )
{
    int nt = vTracks.nTracks;

    //! reversed T (to parallel with Ct)
    ARBB_CPP_ALIGN(ftype *TBuf[6]);
    for( int i = 0; i < 6; i ++ ){
        TBuf[i] = ( ftype * ) arbb::aligned_malloc( sizeof( ftype ) * nt );
    }
    //! reversed C
    ARBB_CPP_ALIGN(ftype *CBuf[15]);
    for( int i = 0; i < 15; i ++ ){
        CBuf[i] = ( ftype * ) arbb::aligned_malloc( sizeof( ftype ) * nt );
    }

    //! ======================= Ct ==============================
    //! global copy in
    gTracksCt = TracksCt<FTYPE>( vTracks );
    gStationsCt = StationsCt<FTYPE>( vStations );

    i32 vSize(nt);

    call(fitTracksCtEntry)(vSize);
    //fitTracksCtEntry(vSize);
    
    //! global copy out
    for( int i = 0; i < 6; i ++ ){
      copy_out(gResults.VecT[i], TBuf[i], nt);
    }
#if defined(CONVENTIONAL)
    const CovV<FTYPE>& tmp = gResults.VecC;
#else
    CovVConventional<FTYPE> tmp = gResults.VecC;
#endif
    copy_out(tmp.C00, CBuf[0], nt);
    copy_out(tmp.C10, CBuf[1], nt);
    copy_out(tmp.C11, CBuf[2], nt);
    copy_out(tmp.C20, CBuf[3], nt);
    copy_out(tmp.C21, CBuf[4], nt);
    copy_out(tmp.C22, CBuf[5], nt);
    copy_out(tmp.C30, CBuf[6], nt);
    copy_out(tmp.C31, CBuf[7], nt);
    copy_out(tmp.C32, CBuf[8], nt);
    copy_out(tmp.C33, CBuf[9], nt);
    copy_out(tmp.C40, CBuf[10], nt);
    copy_out(tmp.C41, CBuf[11], nt);
    copy_out(tmp.C42, CBuf[12], nt);
    copy_out(tmp.C43, CBuf[13], nt);
    copy_out(tmp.C44, CBuf[14], nt);
    //! ======================= Ct ==============================

    //! transpose results
    for( int i = 0; i < nt; i ++ ){
        for( int j = 0; j < 6; j ++ ){
            T[i][j] = TBuf[j][i];
        }
        for( int j = 0; j < 15; j ++ ){
            C[i][j] = CBuf[j][i];
        }
    }

    //! release memory
    for( int i = 0; i < 6; i ++ ){
        arbb::aligned_free( TBuf[i] );
    }
    for( int i = 0; i < 15; i ++ ){
        arbb::aligned_free( CBuf[i] );
    }
}