#-------------------------------------------------------------
#
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.
#
#-------------------------------------------------------------

# THIS SCRIPT IMPLEMENTS KNN( K Nearest Neighbor ) ALGORITHM
#
# INPUT   PARAMETERS:
# ---------------------------------------------------------------------------------------------
# NAME    TYPE     DEFAULT     OPTIONAL     MEANING
# ---------------------------------------------------------------------------------------------
# X                     String   ---   N    Location to read the input matrix as features
# T                     String   ---   N    Location to read the input matrix for nearest neighbor
#                                           search
# Y                     String   ---    Y   Location to read the input matrix as target
# Y_T                   String   ---    Y   Location to read the target type of matrix Y whether
#                                           columns in Y are continuous ( =1 ) or
#                                           categorical ( =2 )
# K                     String   ---    Y   Location to write the result of selected k value
# PR                    String   ---    Y   Location to write the predicted value result,
#                                           it's optional when Y is NOT specified.   
# NNR                   String   ---    N   Location to write the nearest neighbor index result for T.   
# trans_continuous      int      0      Y   Option flag for continuous feature transformed to [-1,1]:
#                                           0 = do not transform continuous variable;
#                                           1 = transform continuous variable;
# select_k              int      0      Y   Use k selection algorithm to estimate k
#                                           ( 1 means yes )
# k_min                 int      1      Y   Min k value(  available if select_k = 1 )
# k_max                 int      100    Y   Max k value(  available if select_k = 1 )
# select_feature        int      0      Y   Use feature selection algorithm to select feature 
#                                           ( 1 means yes )
# START_SELECTED        String   ---    Y   Location to read feature selection initinal value
# FEATURE_SELECTED      String   ---    Y   Location to write feature selection result
# feature_max           int      10     Y   Max feature selection
# interval              int      1000   Y   Interval value for K selecting (  available if select_k = 1 )
# feature_importance    int      0      Y   Use feature importance algorithm to estimate each feature
#                                           ( 1 means yes )
# FEATURE_IMPORTANCE_VALUE String --- Y Location to write feature importance result
# k_value               int      5      Y   k value for KNN, ignore if select_k enable
# predict_con_tg        int      0      Y   Continuous  target predict function: mean(=0) or
#                                           median(=1)
# fmt                   String   "text" Y   Matrix output format for MB, ML, usually "text" or "csv"
# ---------------------------------------------------------------------------------------------
# OUTPUT: Matrix PR, Matrix NNR, Matrix FEATURE_IMPORTANCE_VALUE
#
# HOW TO INVOKE THIS SCRIPT - EXAMPLE:
# hadoop jar SystemDS.jar -f knn.dml -nvargs X=%datapath%/train.csv
#                       T=%datapath%/test.csv NNR=%tmpdir%/NNR.csv fmt=csv


print( "start knn algorithm" );
#input data
TrainPath = $X;
TestPath = $T;
CLPath = ifdef( $Y," " );
CLTypePath = ifdef( $Y_T," " );
#output data
NNRPath = $NNR;
PRPath = ifdef( $PR," " );
#parameter
epsilon = ifdef( $epsilon,0 );
distance_type = ifdef( $distance_type,1 );
trans_continuous = ifdef( $trans_continuous,0 );
predict_con_tg = ifdef($predict_con_tg,0);
k_value = ifdef( $k_value,5 );
FMT = ifdef( $fmt,"text" );
#k select parameter
select_k = ifdef( $select_k,0 );
k_min = ifdef( $k_min,1 );
k_max = ifdef( $k_max,100 );
#feature select parameter
select_feature = ifdef( $select_feature,0 );
feature_max = ifdef( $feature_max,10 );
#feature importance parameter
feature_importance = ifdef ($feature_importance,0);


#Read train data
Train = read ( TrainPath );
Test = read ( TestPath );

#data prepare
if( trans_continuous == 1 ){
    print( "Do data prepare to transfer continuous variable to [-1,1] ..." );
    Train = prepareKNNData( Train);
    Test  = prepareKNNData( Test);
}

n_records = nrow( Train );
n_features = ncol( Train );
s_selected_k = 5;
m_selected_feature = matrix(1,rows=1,cols=n_records);
if( select_k == 1 | select_feature==1 ){
    default_interval = as.integer( ceil( sqrt( n_records ) ) );
    interval = ifdef( $interval,default_interval );
    print( "interval is " + interval );
    #parameter check
    if(  CLPath == " " | CLTypePath == " "){
        stop( "Y or Y_T must be set for k selection algorithm" );
    }
    #parameter re-define
    if( select_k==1 ){
        if(  k_max >= n_records  ){
            k_max = n_records - 1;
            print( "k_max should no greater than number of record, change k_max equal " +
            "( number of record - 1 ) : " + k_max );
        }
        if(  k_max >= interval  ){
            interval = k_max + 1;
            # k_max should equal interval -1, because we drop self when search nn.
            print( "interval should be no less than k_max, change interval equal : " + 
            interval );
        }
        if(  k_max <= 1  ){
            stop( "uncorrect k_max value"  );
        }   
        if(  k_min >= k_max ){
            stop( "k_min >= k_max" );
        }   
    }
    if( select_feature == 1 ){
        if(  k_value >= n_records  ){
            k_value = n_records - 1;
            print( "k_value should be no greater than number of record, change k_value equal " +
            "( number of record - 1 ) : " + k_value );
        }
        #Select feature only
        StartSelectedPath = ifdef( $START_SELECTED," " );
        SelectedPath = ifdef( $FEATURE_SELECTED," " );
        if( StartSelectedPath == " " ){
            m_start_selected_feature = matrix( 0,1,n_features );
        }else{
            m_start_selected_feature = read( StartSelectedPath );
        }       
    }
    CL = read( CLPath );
    CLType = read( CLTypePath );
    type = as.integer( as.scalar( CLType[1,1] ) );
    if( select_k == 1 & select_feature == 1){
        #Combined k and feature selection
        print("Start combined k and feature selection ...");
        [m_selected_feature,s_selected_k] = 
            getSelectedFeatureAndK( Train,CL,type,m_start_selected_feature,
                feature_max,k_min,k_max,interval );
    }else if( select_k == 1 ){
        #Select k only
        print("Start k select ...");
        s_selected_k = getSelectedKBase( Train,CL,type,k_min,k_max,interval );
    }else if( select_feature == 1 ){
        #Select feature only
        print("Start feature selection ... ");
        [m_selected_feature,d_err] = 
            getSelectedFeature( Train,CL,type,m_start_selected_feature,
                feature_max,k_value,interval );
    }
}
if( feature_importance == 1){
    if(  k_value >= n_records  ){
        k_value = n_records - 1;
        print( "k_value should be no greater than number of record, make k_value equal " + 
        "( number of record - 1 ) : " + k_value );
    }
    CL = read( CLPath );
    CLType = read( CLTypePath );
    type = as.integer( as.scalar( CLType[1,1] ) );  
    FeatureImportanceExportPath = ifdef( $FEATURE_IMPORTANCE_VALUE," " );
    if(FeatureImportanceExportPath != " "){
        [m_feature_importance,m_norm_feature_importance] = 
            getFeatureImportance(Train,CL,type,k_value);
        #print( "feature importance is "  );
        #tmp = matrixDisplay( m_feature_importance );
        write ( m_feature_importance,FeatureImportanceExportPath,format=FMT );          
    }
}

# knn search
print("Start knn search ...")
# apply k selection result
if(select_k == 1){
    k_value = s_selected_k
}
if(  k_value >= n_records ){
    k_value = n_records - 1;
    print( "k_value should be no greater than number of record, make k_value equal " + 
        "( number of record - 1 )" );
}
print( "k value is " + k_value );
# apply feature selection result
if( select_feature == 1 ){
    if( n_features > sum( m_selected_feature ) ){
        Train = removeEmpty( target = Train,margin="cols",select = m_selected_feature);
        Test = removeEmpty( target = Test,margin="cols",select = m_selected_feature );
    }
}

NNR_matrix = naiveKNNsearch(Train,Test,k_value)
write( NNR_matrix,NNRPath,format=FMT );

# main part for knn predicted
if( CLPath != " " ){
    CL = read ( CLPath );
    if(  CLTypePath == " " ){
        stop( "Set Y but Miss Y_T for label type" );
    }
    CLType = read( CLTypePath );
    i_cl_type = as.integer( as.scalar( CLType[1,1] ) );
}
CL_matrix = matrix( 0,nrow( Test ),1 );
if( CLPath != " " ){
    parfor( i in 1:nrow( NNR_matrix ) ,check=0){
        NNR_tmp_matrix = matrix( 0,k_value,1 );
        for( j in 1:k_value ){
            NNR_tmp_matrix[j,1] = CL[as.scalar( NNR_matrix[i,j] ),1] ;
        }
        t_cl_value = 0.0 #test data label value
        if( i_cl_type == 2 ){
            t_cl_value = as.scalar( rowIndexMax( t( table( NNR_tmp_matrix,1 ) ) ) )
        }else{
            if ( predict_con_tg == 0){
                t_cl_value = mean( NNR_tmp_matrix );
            } else if(predict_con_tg == 1){
                t_cl_value = median( NNR_tmp_matrix );
            }
        }
        CL_matrix[i,1] = t_cl_value;
    }
}
if( PRPath != " " ){
    write( CL_matrix,PRPath,format=FMT );
}

#naive knn search implement
naiveKNNsearch = function(
    matrix[double] P,
    matrix[double] Q,
    integer K
)return(
    matrix[double] O
){
    num_records = nrow (P);
    num_features = ncol (P);
    num_queries = nrow (Q);
    Qt = t(Q);
    PQt = P %*% Qt;
    P2 = rowSums (P ^ 2);
    D = -2 * PQt + P2;
    if (K == 1) {
        Dt = t(D);
        O = rowIndexMin (Dt);
    } else {
        O = matrix (0, rows = num_queries, cols = K);
        parfor (i in 1:num_queries) {
            D_sorted=order(target=D[,i], by=1, decreasing=FALSE, index.return=TRUE);
            O[i,] = t(D_sorted[1:K,1]);
        }   
    }   
}
#naive knn search for predict value only implement
naiveKNNsearchForPredict = function(
    matrix[double] P,
    matrix[double] Q,
    matrix[double] L,
    integer K
)return(
    matrix[double] OL
){
    num_records = nrow (P);
    num_features = ncol (P);
    num_queries = nrow (Q);
    Qt = t(Q);
    PQt = P %*% Qt;
    P2 = rowSums (P ^ 2);
    D = -2 * PQt + P2;
    if (K == 1) {
        Dt = t(D);
        O = rowIndexMin (Dt);
        OL = matrix (0, rows = num_queries, cols = 1)
        parfor( i in 1:num_queries){
            OL[i,] = L[as.scalar(O[i,1]),1]
        }
    } else {
        OL = matrix (0, rows = num_queries, cols = K);
        parfor (i in 1:num_queries) {
            D_sorted=order(target=cbind(D[,i],L), by=1, decreasing=FALSE, index.return=FALSE);
            OL[i,] = t(D_sorted[1:K,2]);
        }   
    }   
}
# prepareKNNData:
#   Do data prepare - [-1,1] transform for continues variable 
# Argument:
# in_m_data                     input matrix as features
prepareKNNData = function ( 
    matrix[double] in_m_data
     )
return( 
    matrix[double] out_m_data
     )
{
    m_colmax = colMaxs(in_m_data);
    m_colmin = colMins(in_m_data);
    out_m_data = 2*( in_m_data - m_colmin )/( m_colmax - m_colmin ) - 1;
}

getKNeighbor = function (  matrix[double] in_m_data,
        matrix[double] in_m_test_data,
        matrix[double] in_m_cl, 
        integer in_i_k_max   ) 
    return (  matrix[double] out_m_neighbor_value  ) {
    #to naive
    m_search_result = naiveKNNsearchForPredict(in_m_data,in_m_test_data,in_m_cl,in_i_k_max + 1)
    out_m_neighbor_value = m_search_result[ ,2:in_i_k_max + 1]
}

getErr_k = function (  matrix[double] in_m_neighbor_value,
        matrix[double] in_m_cl,
        integer in_i_cl_type , 
        integer in_i_k_min  ) 
    return (  matrix[double] out_m_err  ){
    i_col = ncol( in_m_neighbor_value  );
    i_row  = nrow( in_m_neighbor_value  );  
    
    out_m_err = matrix( 0,i_row,i_col - in_i_k_min + 1  );
    if( in_i_cl_type == 2  ){
        #category
        m_correct = in_m_neighbor_value != in_m_cl[1:i_row,];
    }else{
        #continues
        m_correct = (in_m_neighbor_value - in_m_cl[1:i_row,])^2;#ppred( in_m_neighbor_value,in_m_cl,"-"  );
    }
    parfor( i in 1:i_col-in_i_k_min+1 ,check = 0 ){
        out_m_err[,i] = 
            ( rowSums( m_correct[,1:in_i_k_min + i - 1]  ) / ( in_i_k_min + i - 1  )  );
    }
    #return err for each record and each k ( belong to range 1~max  );
}
/*
eliminateModel = function (  matrix[double] in_j_err, matrix[double] in_j1_err   )
    return(  boolean out_b_inactived   ){
    out_b_inactived = FALSE;
    i_row  = nrow( in_j_err  ); 
    m_err_mean = as.scalar( colSums( in_j_err - in_j1_err  ) / i_row  );
    m_err_variance = as.scalar( 
        sqrt(  colSums( ( in_j_err - in_j1_err - m_err_mean  ) ^ 2  ) / ( i_row-1  )))
    
    #alpha, beta, gamma, delta
    d_gamma = 0.001;
    d_delta = 0.001;
    
    tmp_d_delta = cdf(  
        target = ( -d_gamma - m_err_mean  ) / m_err_variance, dist="t",df=i_row-1  );
    #print( "m_err_mean=" + m_err_mean + ",m_err_variance=" + 
    #   m_err_variance + ",tmp_d_delta=" + tmp_d_delta  );
    if( tmp_d_delta < d_delta  ){
        out_b_inactived = TRUE;
    }

}*/
eliminateModel = function (  double s_err_mean, double s_err_vars, integer i_row  )
    return(  boolean out_b_inactived   ){
    out_b_inactived = FALSE;
    #alpha, beta, gamma, delta
    d_gamma = 0.001;
    d_delta = 0.001;
    
    tmp_d_delta = cdf(  
        target = ( -d_gamma - s_err_mean  ) / s_err_vars, dist="t",df=i_row-1  );
    #print( "m_err_mean=" + m_err_mean + ",m_err_variance=" + 
    #   m_err_variance + ",tmp_d_delta=" + tmp_d_delta  );
    if( tmp_d_delta < d_delta  ){
        out_b_inactived = TRUE;
    }

}
# getSelectedKBase:
#   k selection algorithm with simple KNN algorithm.  
#   Refer to ADD part "6.3 Optimized LOOCV + BRACE"
# Argument:
# in_m_data                     input matrix as features
# in_m_data_target              input matrix as target value
# in_i_is_categorical           1 = category , 0 = continuous
# k_min                         minimun k
# k_max                         maximun k
# interval                      block size for BRACE algorithm
#
# Reture:
# k                             output k value for k selection
getSelectedKBase = function ( 
    matrix[double] in_m_data,
    matrix[double] in_m_data_target,
    integer in_i_is_categorical, # 1 = category , 0 = continuous
    integer k_min,
    integer k_max,
    integer interval  )
return( 
    integer k
      )
{
    b_continue_loop = TRUE;
    i_iter = 1;
    i_record  = nrow( in_m_data  );

    m_active_flag = matrix( 0,1,k_max-k_min+1  );
    i_active_model_number = k_max - k_min + 1;
    
    m_iter_err_sum = matrix( 0,1,k_max-k_min+1);
    m_iter_err_sum_squared = matrix( 0,1,k_max-k_min+1);
    while( b_continue_loop  ){
        # 1.build k-d tree? , or use hash method
        # 2.search data to get k_max nearest neighbor
        i_process_item = i_iter*interval;
        if(  i_process_item >= i_record  ){
            i_process_item = i_record;
            b_continue_loop = FALSE;
        }
        i_process_begin_item = (i_iter - 1)*interval + 1
		#i_process_begin_item = 1
        i_process_end_item = i_process_item
        
        m_neighbor_value = getKNeighbor( in_m_data,in_m_data[i_process_begin_item:i_process_end_item,],in_m_data_target,k_max  );
        # 3.get matrix of err from k_min to k_max 
        m_err = getErr_k( m_neighbor_value,in_m_data_target[i_process_begin_item:i_process_end_item,],in_i_is_categorical,k_min  );

        # 4.check this matrix to drop unnessary record
        m_active_flag_tmp = matrix( 0,1,ncol( m_err  )  );
        
        s_rows_number = i_process_item#nrow(m_err);
            
        m_iter_err_sum = colSums(m_err) + m_iter_err_sum
        m_iter_err_sum_squared = colSums(m_err ^ 2) + m_iter_err_sum_squared
        
        m_err_mean = - outer(t(m_iter_err_sum), m_iter_err_sum , "-") / s_rows_number
        m_err_vars = ( m_err_mean ^2 * s_rows_number -
            2 * m_err_mean * m_iter_err_sum  + m_iter_err_sum_squared) / (s_rows_number-1)
        m_err_vars = sqrt(m_err_vars)
        
        parfor( i in 1:ncol( m_err  )  ){
            #m_err_i = m_err - m_err[,i];
            #m_err_i_mean = colMeans(m_err_i);
            #m_err_i_vars = colVars(m_err_i);
            parfor( j in 1:ncol( m_err  ),check=0  ){
                b_execute_block = TRUE;
                if( j==i  ) b_execute_block = FALSE;
                if( as.scalar( m_active_flag_tmp[1,i]  ) == 1  )  b_execute_block = FALSE; 
                #i has dropped, ignore this case 
                if( as.scalar( m_active_flag_tmp[1,j]  ) == 1  )  b_execute_block = FALSE; 
                #j has dropped, ignore this case
                if( b_execute_block  ){
                    b_flag = eliminateModel( as.scalar(m_err_mean[i,j]),as.scalar(m_err_vars[i,j]),s_rows_number);
                    if(  b_flag == TRUE  ){
                        m_active_flag_tmp[1,i] = 1;
                    }
                }
            }
        }
        
        #m_tmp = ppred( colSums( m_active_flag_tmp  ),1,">="  );
        m_active_flag =  m_active_flag + m_active_flag_tmp >= 1;
        i_active_model_number = -sum( m_active_flag - 1  );
        
        # 5.break while check
        if( i_active_model_number <= 1  ){
            b_continue_loop = FALSE;
        }
        i_iter = i_iter + 1;
        print("i_iter" + i_iter)
    }

    k = 0;
    if(  i_active_model_number == 0   ){
        print( "All k kick out, use min of range " + k_min  );
        k = k_min;
    }else if( i_active_model_number == 1  ){
        k = k_min + as.integer( as.scalar( rowIndexMin( m_active_flag  )  )  ) - 1;
        print( "Get k, which value is " + k  );
    }else{
        m_err_for_order = 
            cbind( t( m_iter_err_sum ), matrix( seq( k_min,k_max,1  ),k_max-k_min+1,1  )  );
        m_err_for_order = removeEmpty( 
            target=m_err_for_order * t( m_active_flag == 0 ),margin="rows"  );
        for( i in 1:nrow( m_err_for_order  )  ){
            print( "k:" + as.scalar( m_err_for_order[i,2]  ) + 
                ",err:" + as.scalar( m_err_for_order[i,1]  )  );
        }
        m_err_order = order( target=m_err_for_order,by=1,decreasing=FALSE,index.return=FALSE  );
        k = as.integer( as.scalar( m_err_order[1,2]  )  );
        print( "Get minimum LOOCV error, which value is " + k  );
    }
}

getErr = function (  matrix[double] in_m_neighbor_value,
        matrix[double] in_m_cl,
        integer in_i_cl_type ) 
    return (  matrix[double] out_m_err ){
    i_col = ncol( in_m_neighbor_value );
    i_row  = nrow( in_m_neighbor_value );   
    
    if( in_i_cl_type == 2 ){
        #category
        m_correct = in_m_neighbor_value != in_m_cl[1:i_row,];
    }else{
        #continues
        m_correct = (in_m_neighbor_value - in_m_cl[1:i_row,])^2;#ppred( in_m_neighbor_value,in_m_cl,"-" );
    }
    out_m_err = ( rowSums( m_correct[,1:i_col] )/( i_col ) );
}

# getSelectedFeature:
#   feature selection algorithm.  
#   Refer to ADD part "7.1 FS+BRACE"
# Argument:
# in_m_data                     input matrix as features
# in_m_data_target              input matrix as target value
# in_i_is_categorical           1 = category , 0 = continuous
# in_m_init_selected            S.user, initial selected feature which use specified
# in_i_max_select               J, max feature selected
# k_value                       k
# interval                      block size for BRACE algorithm
#
# Reture:
# out_m_selected_feature        output matrix for feature selection
# out_d_min_LOOCV               output err

getSelectedFeature = function ( 
    matrix[double] in_m_data,
    matrix[double] in_m_data_target,
    integer in_i_is_categorical, # 1 = category , 0 = continuous
    matrix[double] in_m_init_selected,
    integer in_i_max_select,
    integer k_value,
    integer interval )
return( 
    matrix[double] out_m_selected_feature,
    double out_d_min_LOOCV
     )
{
    i_n_record  = nrow( in_m_data );
    i_n_column = ncol( in_m_data );
    m_main_selected_flag = in_m_init_selected;
    b_no_feature_selected = TRUE;
    if( sum( in_m_init_selected ) >= 1 ){
        b_no_feature_selected = FALSE;
    }
    d_max_err_value = ( max( in_m_data_target ) - min( in_m_data_target ) ) * 100;
    b_continue_main_loop = TRUE;    #level 1 while loop flag
    d_min_LOOCV = Inf;
    while( b_continue_main_loop ){
        m_feature_selected_flag = m_main_selected_flag;
        m_this_model_selected_flag = TRUE;
        i_index_min_LOOCV = -1; # flag for which model win in BRACE algorith
        b_selected_morethan_one = FALSE;
        b_continue_loop = TRUE; #level 2 while loop flag
        i_iter = 1;
        m_iter_err_sum = matrix( 0,1,i_n_column+1 );
        m_iter_err_sum_squared = matrix( 0,1,i_n_column+1 );
        while( b_continue_loop ){
            #print( "*************" + i_iter + "*************" );
            i_process_item = i_iter*interval;
            if(  i_process_item >= i_n_record ){
                i_process_item = i_n_record;
                b_continue_loop = FALSE;
            }
            i_process_begin_item = (i_iter - 1)*interval + 1
            #i_process_begin_item = 1
            i_process_end_item = i_process_item
            m_err = matrix( d_max_err_value,i_process_end_item - i_process_begin_item + 1,i_n_column+1 );
            if( b_no_feature_selected == TRUE ){
                parfor( i in 1:i_n_column ,check=0){
                    if( as.scalar( m_feature_selected_flag[1,i] ) != 1 ){
                        m_tmp_process_data = in_m_data[,i];
                        m_neighbor_value = getKNeighbor( 
                            m_tmp_process_data,m_tmp_process_data[i_process_begin_item:i_process_end_item,],
                            in_m_data_target,k_value );                            
                        m_tmp_err = getErr( 
                            m_neighbor_value,in_m_data_target[i_process_begin_item:i_process_end_item,],in_i_is_categorical );
                        m_err[,i] = m_tmp_err;
                    }
                }
            }else{
                #Use m_main_selected_flag but not m_feature_selected_flag,
                # m_main_selected_flag: which feature are init selected
                # m_feature_selected_flag: which feature are dropped & init selected
                m_tmp_data = removeEmpty( target=in_m_data ,margin="cols", select = m_main_selected_flag);
                if( m_this_model_selected_flag == TRUE ){
                    m_neighbor_value = getKNeighbor( 
                        m_tmp_data,m_tmp_data[i_process_begin_item:i_process_end_item,],in_m_data_target, k_value );
                    
                    m_tmp_err = getErr( m_neighbor_value,in_m_data_target[i_process_begin_item:i_process_end_item,],in_i_is_categorical );
                    m_err[,i_n_column+1] = m_tmp_err;
                }
                parfor( i in 1:i_n_column ,check=0){
                    if( as.scalar( m_feature_selected_flag[1,i] ) != 1 ){
                        m_tmp_process_data = cbind( m_tmp_data,in_m_data[,i] );
                        m_neighbor_value = getKNeighbor( 
                            m_tmp_process_data,m_tmp_process_data[i_process_begin_item:i_process_end_item,],in_m_data_target,k_value );
                        m_tmp_err = getErr( 
                            m_neighbor_value,in_m_data_target[i_process_begin_item:i_process_end_item,],in_i_is_categorical );
                        m_err[,i] = m_tmp_err;
                    }
                }
            }
            if( m_this_model_selected_flag == TRUE ){
                m_active_flag_tmp = cbind( m_feature_selected_flag,matrix( 0,1,1 ) );
            }else{
                m_active_flag_tmp = cbind( m_feature_selected_flag,matrix( 1,1,1 ) );
            }
            
            s_rows_number = i_process_item#nrow(m_err);
            
            m_iter_err_sum = colSums(m_err) + m_iter_err_sum
            m_iter_err_sum_squared = colSums(m_err ^ 2) + m_iter_err_sum_squared
            
            m_err_mean = - outer(t(m_iter_err_sum), m_iter_err_sum , "-") / s_rows_number
            m_err_vars = ( m_err_mean ^2 * s_rows_number -
                    2 * m_err_mean * m_iter_err_sum  + m_iter_err_sum_squared) / (s_rows_number-1)
            m_err_vars = sqrt(m_err_vars)
            
            parfor( i in 1:ncol( m_err ) ){
                #m_err_i_mean = (m_iter_err_sum - m_iter_err_sum[,i])/s_rows_number
                #m_err_i_vars = (m_iter_err_sum_squared + m_err_i_mean ^2 * s_rows_number -
                #   2 * m_iter_err_sum * m_err_i_mean) / (s_rows_number-1)
                #m_err_i_vars = sqrt(m_err_i_vars)
                #m_err_i_mean = m_err_mean[i,]
                #m_err_i_vars = m_err_vars[i,]
                parfor( j in 1:ncol( m_err ) ,check=0){
                    b_execute_block = TRUE;
                    if( j==i ) b_execute_block = FALSE;
                    if( as.scalar( m_active_flag_tmp[1,i] ) == 1 )  b_execute_block = FALSE; 
                    #i has dropped, ignore this case 
                    if( as.scalar( m_active_flag_tmp[1,j] ) == 1 )  b_execute_block = FALSE; 
                    #j has dropped, ignore this case
                    if( b_execute_block ){
                        #print( "i=" + i + ",j=" + j );
                        b_flag = eliminateModel( as.scalar(m_err_mean[i,j]),as.scalar(m_err_vars[i,j]),s_rows_number);
                        if(  b_flag == TRUE ){
                            m_active_flag_tmp[1,i] = 1;
                        }
                    }
                }
            }
            #We mark bit to 1 for selected feature before current loop, 
            #and mark bit to 1 also for dropped feature in current loop
                        
            if( sum( m_active_flag_tmp != 1 ) > 1 ){
                b_selected_morethan_one = TRUE;
            }
            m_col_sums_err = m_iter_err_sum #colSums( m_err );
            i_index_min_LOOCV = as.scalar( rowIndexMin( m_col_sums_err ) );
            d_min_LOOCV = as.scalar( m_col_sums_err[1,i_index_min_LOOCV] );
            i_index_min_LOOCV = i_index_min_LOOCV%% ( i_n_column+1 )
            
            if( sum( m_active_flag_tmp != 1 ) <= 1 ){
                b_continue_loop = FALSE;
            }
            if( as.scalar( m_active_flag_tmp[1,i_n_column+1] ) == 1 ){
                m_this_model_selected_flag = FALSE;
            }
            m_feature_selected_flag = m_active_flag_tmp[,1:i_n_column];
            i_iter = i_iter + 1;            
        }
        #select current model, jump out.
        if( i_index_min_LOOCV == 0 ){
            b_continue_main_loop = FALSE;
            print( "Select Current model" );
        }else{
            print( "select feature " + i_index_min_LOOCV + ", change bit value to 1" );
            m_main_selected_flag[1,i_index_min_LOOCV] = 1;
            b_no_feature_selected = FALSE;
        }
        if( sum( m_main_selected_flag - in_m_init_selected ) >= in_i_max_select ){
            #select more than 10
            b_continue_main_loop = FALSE;
        }
        if( sum( m_main_selected_flag ) == i_n_column ){
            #all selected
            b_continue_main_loop = FALSE;
        }
    }
    out_m_selected_feature = m_main_selected_flag;
    out_d_min_LOOCV = d_min_LOOCV;
}

# getSelectedFeatureAndK:
#   Combine k and feature selection algorithm.  
#   Refer to ADD part "8.Combined k and feature selection"
# Argument:
# in_m_data                     input matrix as features
# in_m_data_target              input matrix as target value
# in_i_is_categorical           1 = category , 0 = continuous
# in_m_init_selected            S.user, initial selected feature which use specified
# in_i_max_select               J, max feature selected
# k_min                         minimun k
# k_max                         maximun k
# interval                      block size for BRACE algorithm
#
# Reture:
# out_m_selected_feature        output matrix for feature selection
# out_i_selected_k              output k value for k selection
getSelectedFeatureAndK = function ( 
    matrix[double] in_m_data,
    matrix[double] in_m_data_target,
    integer in_i_is_categorical, # 1 = category , 0 = continuous
    matrix[double] in_m_init_selected,
    integer in_i_max_select,
    integer k_min,
    integer k_max,
    integer interval )
return( 
    matrix[double] out_m_selected_feature,
    integer out_i_selected_k
     )
{
    m_err = matrix( 0,1,k_max-k_min+1 );
    m_feature = matrix( 0,k_max-k_min+1,ncol( in_m_data ) );
    
    #Step 1. For each k in [k_min,k_max] ( k_min has default value 1, k_max has default value 100 )
    #in parallel select relevant features using FS+BRACE or schemata search described in Section 7.
    parfor( i in k_min:k_max,check=0 ){
        #print( " ************** " + i  + " ****************" );
        [m_selected_feature,d_err] = 
            getSelectedFeature( in_m_data,in_m_data_target,in_i_is_categorical
            ,in_m_init_selected,in_i_max_select,i,interval );
        
        m_err[1,i] = d_err;
        m_feature[i,] = m_selected_feature;
        #print( "err:" + d_err );
    }
    
    #Step 2. Output the combination of features and k with the smallest LOOCV error.
    i_min_err_index = as.integer( as.scalar( rowIndexMin( m_err ) ) );
    out_m_selected_feature = m_feature[i_min_err_index,];
    out_i_selected_k = i_min_err_index + k_min - 1;
    #print( "min err k is " + out_i_selected_k );
}

getFeatureImportance = function ( 
    matrix[double] in_m_data,
    matrix[double] in_m_data_target,
    integer in_i_is_categorical, # 1 = category , 0 = continuous
    integer k_value)
return( 
    matrix[double] out_m_feature_importance,
    matrix[double] out_m_norm_feature_importance
     )
{
    n_feature = ncol(in_m_data)
    n_record = nrow(in_m_data)
    if(n_feature <= 1){
        stop("can't estimate feature importance when ncol = 1")
    }
    m_err = matrix( 0,n_record,n_feature);
    parfor(i_feature in 1:n_feature, check=0){
        m_feature_select = matrix(1,1,n_feature)
        m_feature_select[1,i_feature] = 0;
        m_in_tmp_data = removeEmpty(target=in_m_data,margin="cols", select= m_feature_select)
        m_neighbor_value = getKNeighbor( m_in_tmp_data,m_in_tmp_data,in_m_data_target,k_value);
        m_tmp_err = getErr( m_neighbor_value,in_m_data_target ,in_i_is_categorical );
        m_err[,i_feature] = m_tmp_err
    }
    out_m_feature_importance = colSums( m_err );
    out_m_norm_feature_importance = 
        out_m_feature_importance / as.scalar(rowSums(out_m_feature_importance))
}