/******************************************************************************
 * @internal
 * @file        rCFD_shift.h
 * @copyright   (C) 2023
 *              Particulate Flow Modelling,
 *              Johannes Kepler University, Linz, Austria
 * @website     www.particulate-flow.at
 * @authors     Stefan Pirker
 *****************************************************************************/

#ifndef RCFD_SHIFT
#define RCFD_SHIFT

#include "rCFD_post.h"
#include "rCFD_parallel.h"
#include "rCFD_utilities.h"
#include <stdio.h>

    void shift_read(rCFD_type* rcfd)
    {
#if RP_NODE

#if 1   // local vars
        int     i_state, i_phase, i_frame, i_shift, i_read, i_read_int, i_read_incoming, i_island;
        int     C2C_index;
        int     *i_island_int, *i_island_incoming;
        int     total_number_of_C2Cs_read = 0;
        tmp_C2C_type    tmp_C2Cs;
#endif

        /* A. Allocate local indices */
        {
#if RP_NODE
            /* local island counters */
            {
                i_island_int = (int*)malloc(rcfd->Solver_Dict->number_of_islands*sizeof(int));
                i_island_incoming = (int*)malloc(rcfd->Solver_Dict->number_of_islands*sizeof(int));
            }
#endif
        }

        /* B. Read files into tmp_shifts; on the fly analysis of number_of_shifts and island_offsets, write tmp_C2Cs into C2Cs */
        {
#if RP_NODE
            char    filename[40];
            FILE    *fi = NULL;
            int     i_tmp, c0, node0, c1, node1, my_island_id, i_layer;
            double  tmp, w0;

            i_layer = 0;

            for(i_state = 0; i_state < rcfd->Solver_Dict->number_of_states; i_state++){

                for(i_phase = 0; i_phase < rcfd->Solver_Dict->number_of_phases; i_phase++){

                    /* try to open file[state][phase] */
                    {
                        sprintf(filename,"%s_%d_%d_%d", rcfd->File_Dict->C2C_filename, i_state, i_phase, myid);

                        fi = fopen(filename,"r");
                    }

                    /* read C2Cs per [state, phase, frame] */
                    if(fi){

                        for(i_frame = 0; i_frame < rcfd->Solver_Dict->number_of_frames; i_frame++){

                            /* B.1. initializations per frame */
                            {
                                i_read = 0; i_read_int = 0; i_read_incoming = 0;

                                for(i_island = 0; i_island < rcfd->Solver_Dict->number_of_islands; i_island++){

                                    i_island_int[i_island] = 0;
                                    i_island_incoming[i_island] = 0;
                                }
                            }

                            /* B.2. allocate tmp_C2Cs, read tmp_C2Cs, count indices */
                            {
                                fscanf(fi,"%d \n",&tmp_C2Cs.number_of_shifts);

                                /* in future, we might read this from file */
                                tmp_C2Cs.format = RCFD_SHIFT_C0_N0_C1_N1_W0;

                                /* allocate tmp_C2Cs.shifts */
                                {
                                    tmp_C2Cs.shifts = (C2C_shift_type*)malloc(tmp_C2Cs.number_of_shifts * sizeof(C2C_shift_type));

                                    for(i_shift = 0; i_shift < tmp_C2Cs.number_of_shifts; i_shift++){

                                        tmp_C2Cs.shifts[i_shift].c0 =     -1;
                                        tmp_C2Cs.shifts[i_shift].node0 =  -1;
                                        tmp_C2Cs.shifts[i_shift].c1 =     -1;
                                        tmp_C2Cs.shifts[i_shift].node1 =  -1;
                                        tmp_C2Cs.shifts[i_shift].w0 =     0.0;
                                    }

                                    tmp_C2Cs.island_offsets =       (int*)malloc((rcfd->Solver_Dict->number_of_islands+1) * sizeof(int));
                                    tmp_C2Cs.island_offsets_in =    (int*)malloc((rcfd->Solver_Dict->number_of_islands+1) * sizeof(int));

                                    for(i_island = 0; i_island < rcfd->Solver_Dict->number_of_islands; i_island++){

                                        tmp_C2Cs.island_offsets[i_island] =     0;
                                        tmp_C2Cs.island_offsets_in[i_island] =  0;
                                    }
                                }

                                /* fill tmp_C2Cs and monitor indices */
                                for(i_shift = 0; i_shift < tmp_C2Cs.number_of_shifts; i_shift++){

                                    if((i_shift % rcfd->Solver_Dict->C2C_loading_reduction) == 0){

                                        if(tmp_C2Cs.format == RCFD_SHIFT_C0_N0_C1_N1_W0){

                                            fscanf(fi,"%d %d %d %d %le\n", &c0, &node0, &c1, &node1, &w0);
                                        }

                                        tmp_C2Cs.shifts[i_read].c0 = c0;
                                        tmp_C2Cs.shifts[i_read].node0 = node0;
                                        tmp_C2Cs.shifts[i_read].c1 = c1;
                                        tmp_C2Cs.shifts[i_read].node1 = node1;
                                        tmp_C2Cs.shifts[i_read].w0 = w0;

                                        i_read++;

                                        my_island_id = _C(i_layer).island_id[c1];

                                        if(node0 == node1){

                                            i_read_int++;
                                            i_island_int[my_island_id]++;
                                        }
                                        else{

                                            i_read_incoming++;
                                            i_island_incoming[my_island_id]++;
                                        }
                                    }
                                    else{

                                        if(tmp_C2Cs.format == RCFD_SHIFT_C0_N0_C1_N1_W0){

                                            fscanf(fi,"%d %d %d %d %le\n", &i_tmp, &i_tmp, &i_tmp, &i_tmp, &tmp);
                                        }
                                    }
                                }
                            }

                            /* B.3 set i_read_outgoing, set island indices */
                            {
                                tmp_C2Cs.island_offsets[0] =        0;
                                tmp_C2Cs.island_offsets_in[0] =     0;

                                for(i_island = 0; i_island < rcfd->Solver_Dict->number_of_islands; i_island++){

                                    tmp_C2Cs.island_offsets[(i_island + 1)] =
                                        tmp_C2Cs.island_offsets[i_island] + i_island_int[i_island];

                                    tmp_C2Cs.island_offsets_in[(i_island + 1)] =
                                        tmp_C2Cs.island_offsets_in[i_island] + i_island_incoming[i_island];
                                }
                            }

                            /* B.4 fill tmp_C2Cs into C2Cs.shifts ... */
                            {
                                /* B.4.1. initializations per frame */
                                {
                                    for(i_island = 0; i_island < rcfd->Solver_Dict->number_of_islands; i_island++){

                                        i_island_int[i_island] = 0;
                                        i_island_incoming[i_island] = 0;
                                    }
                                }

                                /* B.4.2. C2C.numbers, allocate and init C2Cs.shifts/shifts_in */
                                {
                                    rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts[0] = i_read_int;

                                    rcfd->C2Cs[i_state][i_phase][i_frame].shifts[0] = (C2C_shift_type*)malloc(i_read_int * sizeof(C2C_shift_type));

                                    for(i_shift = 0; i_shift < rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts[0]; i_shift++){

                                        rcfd->C2Cs[i_state][i_phase][i_frame].shifts[0][i_shift].c0 =        -1;
                                        rcfd->C2Cs[i_state][i_phase][i_frame].shifts[0][i_shift].node0 =     -1;
                                        rcfd->C2Cs[i_state][i_phase][i_frame].shifts[0][i_shift].c1 =        -1;
                                        rcfd->C2Cs[i_state][i_phase][i_frame].shifts[0][i_shift].node1 =     -1;
                                        rcfd->C2Cs[i_state][i_phase][i_frame].shifts[0][i_shift].w0 =        0.0;
                                    }

                                    rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts_in[0] = i_read_incoming;

                                    rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[0] = (C2C_shift_type*)malloc(i_read_incoming * sizeof(C2C_shift_type));

                                    for(i_shift = 0; i_shift < rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts_in[0]; i_shift++){

                                        rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[0][i_shift].c0 =     -1;
                                        rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[0][i_shift].node0 =  -1;
                                        rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[0][i_shift].c1 =     -1;
                                        rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[0][i_shift].node1 =  -1;
                                        rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[0][i_shift].w0 =     0.0;
                                    }

                                    total_number_of_C2Cs_read +=  rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts[0]
                                                                + rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts_in[0];
                                }

                                /* B.4.3. fill C2Cs.shifts/shifts_in by tmp_C2Cs*/
                                {
                                    for(i_shift = 0; i_shift < tmp_C2Cs.number_of_shifts; i_shift++){

                                        c1 = tmp_C2Cs.shifts[i_shift].c1;

                                        my_island_id = _C(i_layer).island_id[c1];

                                        node0 = tmp_C2Cs.shifts[i_shift].node0;
                                        node1 = tmp_C2Cs.shifts[i_shift].node1;

                                        if(node0 == node1){

                                            C2C_index = tmp_C2Cs.island_offsets[my_island_id] + i_island_int[my_island_id];

                                            rcfd->C2Cs[i_state][i_phase][i_frame].shifts[0][C2C_index].c0 =    tmp_C2Cs.shifts[i_shift].c0;
                                            rcfd->C2Cs[i_state][i_phase][i_frame].shifts[0][C2C_index].node0 = tmp_C2Cs.shifts[i_shift].node0;
                                            rcfd->C2Cs[i_state][i_phase][i_frame].shifts[0][C2C_index].c1 =    tmp_C2Cs.shifts[i_shift].c1;
                                            rcfd->C2Cs[i_state][i_phase][i_frame].shifts[0][C2C_index].node1 = tmp_C2Cs.shifts[i_shift].node1;
                                            rcfd->C2Cs[i_state][i_phase][i_frame].shifts[0][C2C_index].w0 =    tmp_C2Cs.shifts[i_shift].w0;

                                            i_island_int[my_island_id]++;
                                        }
                                        else{

                                            C2C_index = tmp_C2Cs.island_offsets_in[my_island_id] + i_island_incoming[my_island_id];

                                            rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[0][C2C_index].c0 =     tmp_C2Cs.shifts[i_shift].c0;
                                            rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[0][C2C_index].node0 =  tmp_C2Cs.shifts[i_shift].node0;
                                            rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[0][C2C_index].c1 =     tmp_C2Cs.shifts[i_shift].c1;
                                            rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[0][C2C_index].node1 =  tmp_C2Cs.shifts[i_shift].node1;
                                            rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[0][C2C_index].w0 =     tmp_C2Cs.shifts[i_shift].w0;

                                            i_island_incoming[my_island_id]++;
                                        }
                                    }

                                    rcfd->C2Cs[i_state][i_phase][i_frame].island_offsets[0][0] = 0;
                                    rcfd->C2Cs[i_state][i_phase][i_frame].island_offsets_in[0][0] = 0;

                                    for(i_island = 0; i_island < rcfd->Solver_Dict->number_of_islands; i_island++){

                                        rcfd->C2Cs[i_state][i_phase][i_frame].island_offsets[0][(i_island + 1)] = tmp_C2Cs.island_offsets[(i_island + 1)];
                                        rcfd->C2Cs[i_state][i_phase][i_frame].island_offsets_in[0][(i_island + 1)] = tmp_C2Cs.island_offsets_in[(i_island + 1)];
                                    }
                                }
                            }

                            /* B.5 free tmp_C2Cs stuff */
                            {
                                free(tmp_C2Cs.shifts);
                                free(tmp_C2Cs.island_offsets);
                                free(tmp_C2Cs.island_offsets_in);
                            }
                        }

                        fclose(fi);
                    }
                }
            }
#endif
        }

        /* C. Free locally vars */
        {
#if RP_NODE
            free(i_island_int);
            free(i_island_incoming);
#endif
        }

        /* D. Set-up C2Cs for upper layers */
        {
#if RP_NODE
            short   debug_this_code = 0;

            if(rcfd->Solver_Dict->number_of_layers > 1){

#if 1       /* local vars */
                int     i_layer, i_upper_shift;
                int     number_of_shifts, number_of_shifts_in_upper_layer, c0, c1, upper_c0, upper_c1;
                double  w0;
#endif

                for(i_state = 0; i_state < rcfd->Solver_Dict->number_of_states; i_state++){

                    for(i_phase = 0; i_phase < rcfd->Solver_Dict->number_of_phases; i_phase++){

                        for(i_frame = 0; i_frame < rcfd->Solver_Dict->number_of_frames; i_frame++){

                            i_layer = 0;

                            while((i_layer + 1) < rcfd->Solver_Dict->number_of_layers){

                                /* D.1 fill shifts of (i_layer + 1) */
                                {
                                    /* TODO remove doublets */
                                    number_of_shifts = rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts[i_layer];
                                    number_of_shifts_in_upper_layer = 0;

                                    for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                                        c0 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts[i_layer][i_shift].c0;
                                        c1 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts[i_layer][i_shift].c1;

                                        upper_c0 = _C(i_layer).parent_cell[c0];
                                        upper_c1 = _C(i_layer).parent_cell[c1];

                                        if(upper_c0 != upper_c1){

                                            number_of_shifts_in_upper_layer++;
                                        }
                                    }

                                    rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts[(i_layer + 1)] = number_of_shifts_in_upper_layer;

                                    /* TODO account for multiple islands */
                                    rcfd->C2Cs[i_state][i_phase][i_frame].island_offsets[(i_layer + 1)][0] = 0;
                                    rcfd->C2Cs[i_state][i_phase][i_frame].island_offsets[(i_layer + 1)][1] = number_of_shifts_in_upper_layer;

                                    rcfd->C2Cs[i_state][i_phase][i_frame].shifts[(i_layer + 1)] = (C2C_shift_type*)malloc(number_of_shifts_in_upper_layer * sizeof(C2C_shift_type));

                                    i_upper_shift = 0;

                                    for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                                        c0 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts[i_layer][i_shift].c0;
                                        c1 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts[i_layer][i_shift].c1;

                                        upper_c0 = _C(i_layer).parent_cell[c0];
                                        upper_c1 = _C(i_layer).parent_cell[c1];

                                        w0 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts[i_layer][i_shift].w0;

                                        if(upper_c0 != upper_c1){

                                            rcfd->C2Cs[i_state][i_phase][i_frame].shifts[(i_layer + 1)][i_upper_shift].node0 = myid;
                                            rcfd->C2Cs[i_state][i_phase][i_frame].shifts[(i_layer + 1)][i_upper_shift].c0 =    upper_c0;
                                            rcfd->C2Cs[i_state][i_phase][i_frame].shifts[(i_layer + 1)][i_upper_shift].node1 = myid;
                                            rcfd->C2Cs[i_state][i_phase][i_frame].shifts[(i_layer + 1)][i_upper_shift].c1 =    upper_c1;
                                            rcfd->C2Cs[i_state][i_phase][i_frame].shifts[(i_layer + 1)][i_upper_shift].w0 =    w0;

                                            i_upper_shift++;
                                        }
                                    }

                                    if((debug_this_code)&&(i_frame == 0)){

                                        report_debug(rcfd, RCFD_REPORT_CATEGORY_READ_DB,
                                                     "[%d] D.2: i_layer %d, number_of_shifts %d/%d",
                                                     myid, i_layer, number_of_shifts, number_of_shifts_in_upper_layer);
                                    }
                                }

                                i_layer++;
                            }
                        }
                    }
                }
            }
#endif
        }

        /* E. set up parallel C2C communication */
        {
#if RP_NODE
            if(rcfd->Solver_Dict->number_of_layers > 1){

                init_upper_layer_parallel_C2Cs(rcfd);
            }

            init_parallel_C2Cs(rcfd);
#endif
        }

        rcfd->Solver->total_number_of_C2Cs_read = PRF_GISUM1(total_number_of_C2Cs_read);
#endif
    }


int  frame_id_normalized_file_exists(rCFD_type* rcfd){
	
	int file_exits = 0;

#if RP_NODE

	if(myid == 0){
		
		char filename[80];

		FILE *fi = NULL;
		
		sprintf(filename, "%s_normalized", rcfd->File_Dict->frame_id_filename);
		
		fi = fopen(filename, "r");
		
		if(fi){
		
			file_exits = 1;
			
			fclose(fi);
		}
	}
	
	file_exits = PRF_GIHIGH1(file_exits);

#endif

	return file_exits;
}

int  frame_id_file_exists(rCFD_type* rcfd){
	
	int file_exits = 0;

#if RP_NODE

	if(myid == 0){
		
		FILE *fi = NULL;
	
		fi = fopen(rcfd->File_Dict->frame_id_filename, "r");
		
		if(fi){
			
			file_exits = 1;
			
			fclose(fi);
		}
	}
	
	file_exits = PRF_GIHIGH1(file_exits);

#endif

	return file_exits;
}

void read_jumps(rCFD_type* rcfd){
#if RP_NODE

	int i_state, i_state2, i_island, i_frame;

	int jump_file_count = 0;

	int frame_count = rcfd->Solver_Dict->number_of_frames;

	char filename[80];

	FILE *fi = NULL;

	// node-0 reads jumps files
	if(myid == 0){
		
		for(i_state = 0; i_state < rcfd->Solver_Dict->number_of_states; i_state++){

			for(i_state2 = 0; i_state2 < rcfd->Solver_Dict->number_of_states; i_state2++){

				for(i_island = 0; i_island < rcfd->Solver_Dict->number_of_islands; i_island++){

					sprintf(filename,"%s_%d_%d_%d", rcfd->File_Dict->Jump_filename, i_state, i_state2, i_island);

					fi = fopen(filename, "r");

					if(fi != NULL){

						for(i_frame = 0; i_frame < frame_count; i_frame++){

							fscanf(fi, "%d \n", &rcfd->Rec->jumps[i_state][i_state2][i_island][i_frame]);
						}

						jump_file_count++;
						
						fclose(fi);
					}
					else{
						
						report_error(rcfd, 
									 RCFD_REPORT_CATEGORY_READ_DB, 
									 "shift_read_rec: jumps file %s doesn't exist", filename);
					}	
				}
			}
		}
	}
	
	// node-0 bcasts jumps
	for(i_state = 0; i_state < rcfd->Solver_Dict->number_of_states; i_state++){

		for(i_state2 = 0; i_state2 < rcfd->Solver_Dict->number_of_states; i_state2++){

			for(i_island = 0; i_island < rcfd->Solver_Dict->number_of_islands; i_island++){
				
				for(i_frame = 0; i_frame < frame_count; i_frame++){

					PRF_BCAST_INT(rcfd->Rec->jumps[i_state][i_state2][i_island], frame_count, 0);
				}
			}
		}
	}

	report_info0(rcfd, 
				 RCFD_REPORT_CATEGORY_READ_DB, 
				 "shift_read_rec: read %d jump files", jump_file_count);
#endif	
}

void read_frame_id_values(rCFD_type* rcfd, FILE *fi){
#if RP_NODE

	if(myid == 0){
			
		int    i_frame, i_id;

		int    frame_count =     rcfd->Solver_Dict->number_of_frames;
		int    frame_id_length = rcfd->Rec_Dict->frame_id_vector_length;

		for(i_frame = 0; i_frame < frame_count; i_frame++){
			
			for(i_id = 0; i_id < frame_id_length; i_id++){
			
				fscanf(fi, "%le ", &rcfd->Rec->frame_id_database[i_frame][i_id]);
			}

			fscanf(fi, "\n");
		}
	}
	//PRF_GSYNC();
#endif	
}

void read_frame_id(rCFD_type* rcfd){
#if RP_NODE

	if(myid == 0){

		FILE    *fi = NULL;

		fi = fopen(rcfd->File_Dict->frame_id_filename,"r");
		
		if(fi){

			read_frame_id_values(rcfd, fi);

			fclose(fi);
			
			report_info(rcfd, RCFD_REPORT_CATEGORY_READ_DB,
					    "read frame_id from %s", 
						rcfd->File_Dict->frame_id_filename);
		}
		else{
		
			report_error(rcfd, RCFD_REPORT_CATEGORY_READ_DB,
					     "read_frame_id: could not open %s", 
						 rcfd->File_Dict->frame_id_filename);
		}   
	}

	PRF_GSYNC();
#endif	
}

void read_frame_id_normalized(rCFD_type* rcfd){
#if RP_NODE

	if(myid == 0){

		char    filename[40];

		FILE    *fi = NULL;
		
		sprintf(filename,"%s_normalized", rcfd->File_Dict->frame_id_filename);

		fi = fopen(filename,"r");
		
		if(!fi){ 
		
			report_error(rcfd, RCFD_REPORT_CATEGORY_READ_DB,
					     "read_frame_id: could not open %s", 
						 filename);
		}
		else{

			// read two header lines - means & vars
			{
				int i_id;
				
				int frame_id_length = rcfd->Rec_Dict->frame_id_vector_length;
				
				for(i_id = 0; i_id < frame_id_length; i_id++){
					
					fscanf(fi, "%le ", &rcfd->Rec->frame_id_means[i_id]);
				}

				fscanf(fi, "\n");

				for(i_id = 0; i_id < frame_id_length; i_id++){
		
					fscanf(fi, "%le ", &rcfd->Rec->frame_id_vars[i_id]);
				}

				fscanf(fi, "\n");
			}

			read_frame_id_values(rcfd, fi);

			fclose(fi);
		}   

	}
#endif	
}

void normalize_frame_id_database(rCFD_type* rcfd){
#if RP_NODE

	if(myid == 0){
		
		int    i_frame, i_id;

		int    frame_count =     rcfd->Solver_Dict->number_of_frames;
		int    frame_id_length = rcfd->Rec_Dict->frame_id_vector_length;
		
		double *frame_id_means = NULL, *frame_id_vars = NULL;

		double **DB_frame_id = rcfd->Rec->frame_id_database;

		// get mean values
		{
			frame_id_means = malloc(frame_id_length * sizeof(double));
			memset(frame_id_means, 0.0, frame_id_length * sizeof(double));

			for(i_frame = 0; i_frame < frame_count; i_frame++){
				
				for(i_id = 0; i_id < frame_id_length; i_id++){
				
					frame_id_means[i_id] += DB_frame_id[i_frame][i_id];
				}
			}
			
			if( frame_count > 0){
				
				for(i_id = 0; i_id < frame_id_length; i_id++){
				
					frame_id_means[i_id] /= (double) frame_count;
				}
			}
		}

		// get variances
		{
			frame_id_vars = malloc(frame_id_length * sizeof(double));
			memset(frame_id_vars, 0.0, frame_id_length * sizeof(double));

			for(i_frame = 0; i_frame < frame_count; i_frame++){
				
				for(i_id = 0; i_id < frame_id_length; i_id++){
				
					frame_id_vars[i_id] += fabs(DB_frame_id[i_frame][i_id] - frame_id_means[i_id]);
				}
			}
			
			if( frame_count > 0){
				
				for(i_id = 0; i_id < frame_id_length; i_id++){
				
					frame_id_vars[i_id] /= (double) frame_count;
				}
			}
		}
		
		// normalize DB_frame_id values
		{
			for(i_frame = 0; i_frame < frame_count; i_frame++){
				
				for(i_id = 0; i_id < frame_id_length; i_id++){
				
					DB_frame_id[i_frame][i_id] -= frame_id_means[i_id];
					
					if(frame_id_vars[i_id] > 0.0){
						
						DB_frame_id[i_frame][i_id] /= frame_id_vars[i_id];
					}
				}
			}
		}
		
		// set Rec->frame_id_means(_vars)
		{
			for(i_id = 0; i_id < frame_id_length; i_id++){
			
				rcfd->Rec->frame_id_means[i_id] = frame_id_means[i_id];
				rcfd->Rec->frame_id_vars [i_id] = frame_id_vars [i_id];
			}			
		}
		
		// free local vars
		{
			free( frame_id_means);
			free( frame_id_vars);
		}

		report_info0(rcfd, 
					 RCFD_REPORT_CATEGORY_READ_DB, 
					 "shift_read_rec: normalized frame_id");					
	}
#endif	
}

void write_normalized_frame_id(rCFD_type* rcfd){
#if RP_NODE

	if(myid == 0){
		
		char filename[80];

		FILE *fi = NULL;
		
		sprintf(filename, "%s_normalized", rcfd->File_Dict->frame_id_filename);
		
		fi = fopen(filename, "w");
		
		if(!fi){
		
			report_error(rcfd, 
			             RCFD_REPORT_CATEGORY_READ_DB, 
						 "shift_read_rec: write_normalized_frame_id: could not open file %s", filename);
			
			return;
		}
   
		int    i_frame, i_id;

		int    frame_count =     rcfd->Solver_Dict->number_of_frames;
		int    frame_id_length = rcfd->Rec_Dict->frame_id_vector_length;

		// two header lines - means & vars
		{
			for(i_id = 0; i_id < frame_id_length; i_id++){
				
				fprintf(fi, "%e ", rcfd->Rec->frame_id_means[i_id]);
			}

			fprintf(fi, "\n");

			for(i_id = 0; i_id < frame_id_length; i_id++){
				
				fprintf(fi, "%e ", rcfd->Rec->frame_id_vars[i_id]);
			}

			fprintf(fi, "\n");
		}
		
		for(i_frame = 0; i_frame < frame_count; i_frame++){
			
			for(i_id = 0; i_id < frame_id_length; i_id++){
			
				fprintf(fi, "%e ", rcfd->Rec->frame_id_database[i_frame][i_id]);
			}

			fprintf(fi, "\n");
		}
		
		fclose(fi);
		
		report_info0(rcfd, 
				 RCFD_REPORT_CATEGORY_READ_DB, 
				 "shift_read_rec: wrote normalized frame_id into %s", filename);					
	}


	//PRF_GSYNC();
#endif
}
	
void shift_read_rec(rCFD_type* rcfd)
{
#if RP_NODE

	if( rcfd->Rec_Dict->format <= RCFD_REC_REPLAY){
		
		if(rcfd->Rec->jumps != NULL){

			read_jumps(rcfd);
		}
		else{
			
			report_error(rcfd, 
			             RCFD_REPORT_CATEGORY_READ_DB, 
						 "shift_read_rec: jumps variable doesn't exist");
		}			
	}
	else if( rcfd->Rec_Dict->format == RCFD_REC_FRAME_ID){

		if( rcfd->Rec_Dict->frame_id_normalize_on){
			
			if( frame_id_normalized_file_exists(rcfd)){

				read_frame_id_normalized(rcfd);
			}
			else if (frame_id_file_exists(rcfd)){

				read_frame_id(rcfd);

				normalize_frame_id_database(rcfd);

				write_normalized_frame_id(rcfd);
			}
			else{

			report_error(rcfd, 
						 RCFD_REPORT_CATEGORY_READ_DB, 
						 "shift_read_rec: frame_id file doesn't exist");
		}
		}
		else{	// no frame_id_normalization
			
			if (frame_id_file_exists(rcfd)){
				
				read_frame_id(rcfd);
			}
			else{

				report_error(rcfd, 
						 RCFD_REPORT_CATEGORY_READ_DB, 
						 "shift_read_rec: frame_id file doesn't exist");
			}
		}
	}

	PRF_GSYNC();
#endif		
}

    void shift_local_data(rCFD_type* rcfd, const int i_state, const int i_layer, const int i_phase, const int i_island, const int i_frame)
    {
#if RP_NODE

#if 1   // local vars
        int     i_shift, i_data;
        int     loop_offset0, loop_offset1, c0, c1;
        double  w0, data0;
#endif

        // parallel mapping to corona
        for(i_data = 0; i_data < rcfd->Phase_Dict[i_phase].number_of_data; i_data++){

            parallel_communicate_data(rcfd, i_layer, i_phase, i_data, map_host_format);
        }

        loop_offset0    = rcfd->C2Cs[i_state][i_phase][i_frame].island_offsets[i_layer][i_island];
        loop_offset1    = rcfd->C2Cs[i_state][i_phase][i_frame].island_offsets[i_layer][(i_island + 1)];

        for(i_shift = loop_offset0; i_shift < loop_offset1; i_shift++){

            c0 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts[i_layer][i_shift].c0;
            c1 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts[i_layer][i_shift].c1;
            w0 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts[i_layer][i_shift].w0;

            if(w0 > 0.0){

                for(i_data = 0; i_data < rcfd->Phase_Dict[i_phase].number_of_data; i_data++){

                    data0 = _C(i_layer).data[i_phase][c0][i_data];

                    _C(i_layer).data_shift[i_phase][c1][i_data] =
                        (w0 * data0 + _C(i_layer).weight_after_shift[c1] * _C(i_layer).data_shift[i_phase][c1][i_data])/
                        (w0         + _C(i_layer).weight_after_shift[c1]);
                }

                _C(i_layer).weight_after_shift[c1] += w0;
            }
        }
#endif
    }

    void shift_data(rCFD_type* rcfd, const int i_phase)
    {
#if RP_NODE
        int i_state, i_layer, i_island, i_frame;
        i_state = rcfd->Solver->current_state;
        i_layer = rcfd->Solver->current_layer;
        i_frame = rcfd->Rec->global_frame;

        if(i_layer > 0)
            return; // TODO: implement for layer>0

        for(i_island = 0; i_island < rcfd->Solver_Dict->number_of_islands; i_island++){

#ifdef RCFD_LIB_ISLANDS
            i_frame = rcfd->Rec->global_frame_of_island[i_island];
#endif
            shift_local_data(rcfd, i_state, i_layer, i_phase, i_island, i_frame);

            parallel_shift_data(rcfd, i_state, i_layer, i_phase, i_island, i_frame);

            parallel_monitor_data_fluxes(rcfd, i_state, i_layer, i_phase, i_island, i_frame);
        }
#endif
    }

    void shift_fill(rCFD_type* rcfd, const int i_phase)
    {
#if RP_NODE

#if 1   // local vars
        int     i_layer, i_frame, i_face, i_data, i_cell;
        int     number_of_unhit_cells, number_of_fill_loops, c0, c1, c;
#endif

        i_layer = rcfd->Solver->current_layer;
        if(i_layer > 0)
            return; // TODO: implement for layer>0

        i_frame = rcfd->Rec->global_frame;

        if(rcfd->Solver_Dict->shifts_from_displacement_fields_on){   // waterproof C2C patterns

            number_of_unhit_cells = 0;
        }
        else{

            number_of_unhit_cells = 1;
        }

        number_of_fill_loops = 0;

        // fill holes with data_shift
        while((number_of_unhit_cells > 0) && (number_of_fill_loops < rcfd->Solver_Dict->max_fill_loops)){

            number_of_fill_loops++;

            /* fill holes by face swaps */
            for(i_face = 0; i_face < _FACE_DICT(i_layer).number_of_int_faces; i_face++){

                c0 = _F(i_layer).c0[i_face];
                c1 = _F(i_layer).c1[i_face];

                // swap c0, c1
                if((_C(i_layer).weight_after_shift[c0] > 0.)&&(_C(i_layer).weight_after_shift[c1] == 0.0)){

                    c = c0; c0 = c1; c1 = c;
                }

                if((_C(i_layer).weight_after_shift[c1] > 0.)&&(_C(i_layer).weight_after_shift[c0] == 0.0)){

                    // fill c0 with c1 data
                    for(i_data = 0; i_data < rcfd->Phase_Dict[i_phase].number_of_data; i_data++){

                        _C(i_layer).data_swap[i_phase][c0][i_data] =  (
                                                                          _C(i_layer).weight_after_shift[c1]
                                                                        * _C(i_layer).data_shift[i_phase][c1][i_data]
                                                                        + _C(i_layer).weight_after_swap[c0]
                                                                        * _C(i_layer).data_swap[i_phase][c0][i_data]
                                                                      )
                                                                      /
                                                                      (
                                                                          _C(i_layer).weight_after_shift[c1]
                                                                        + _C(i_layer).weight_after_swap[c0]
                                                                      );
                    }

                    _C(i_layer).weight_after_swap[c0] += _C(i_layer).weight_after_shift[c1];
                }
            }

            number_of_unhit_cells = 0;

            for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_int_cells; i_cell++){

                if(_C(i_layer).weight_after_swap[i_cell] > 0.0){

                    for(i_data = 0; i_data < rcfd->Phase_Dict[i_phase].number_of_data; i_data++){

                        _C(i_layer).data_shift[i_phase][i_cell][i_data] = _C(i_layer).data_swap[i_phase][i_cell][i_data];
                    }

                    _C(i_layer).weight_after_shift[i_cell] = _C(i_layer).weight_after_swap[i_cell];
                    _C(i_layer).weight_after_swap[i_cell] = 0.0;
                }

#ifdef RCFD_LIB_ISLANDS

                i_frame = rcfd->Rec->global_frame_of_island[_C(i_layer).island_id[i_cell]];
#endif // RCFD_LIB_ISLANDS

                if((_C(i_layer).weight_after_shift[i_cell] == 0.0) && (_C(i_layer).vof[i_frame][i_cell][i_phase] > 0.0)){

                    number_of_unhit_cells++;
                }
            }
        }

        // data = data_shift
        for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_int_cells; i_cell++){

            if((rcfd->Solver_Dict->shifts_from_displacement_fields_on) || (_C(i_layer).weight_after_shift[i_cell] > 0.0)){

                for(i_data = 0; i_data < rcfd->Phase_Dict[i_phase].number_of_data; i_data++){

                    _C(i_layer).data[i_phase][i_cell][i_data] = _C(i_layer).data_shift[i_phase][i_cell][i_data];
                }
            }
        }

        // parallel mapping to corona
        for(i_data = 0; i_data < rcfd->Phase_Dict[i_phase].number_of_data; i_data++){

            parallel_communicate_data(rcfd, i_layer, i_phase, i_data, map_host_format);
        }

#endif // RP_NODE
    }

    void shift_add_to_list(rCFD_type* rcfd, const int i_layer, const int i_phase, int c0, int node0, double w0, int c1, int node1)
    {
#if RP_NODE

#if 1
        int     shift_count, max_number_of_shifts;
#endif

        shift_count = rcfd->Tracer->monitor_counter[i_phase];

        max_number_of_shifts = _CELL_DICT(i_layer).number_of_int_cells;

        if((shift_count >= 0) && (shift_count < max_number_of_shifts)){

            rcfd->Tracer->shifts[i_phase][shift_count].c0 =       c0;
            rcfd->Tracer->shifts[i_phase][shift_count].node0 =    node0;
            rcfd->Tracer->shifts[i_phase][shift_count].w0 =       w0;
            rcfd->Tracer->shifts[i_phase][shift_count].c1 =       c1;
            rcfd->Tracer->shifts[i_phase][shift_count].node1 =    node1;

            rcfd->Tracer->monitor_counter[i_phase]++;
        }
        else{

            report_warning(rcfd, RCFD_REPORT_CATEGORY_CREATE_SHIFT,
                           "[%d] shift_add_to_list: shift_count %d out of range",
                           myid, shift_count);
        }
#endif
    }

    int shift_from_displacement_fields(rCFD_type* rcfd, int i_layer, int i_phase)
    {
#if RP_NODE

#define LOCAL_DEBUG 0

#if 1
        int     i_dim, i_cell, i_cell2, i_nnb, i_node, i_shift;
        int     cal_1st_cell_count, cal_next_neigh_cell_count, cal_all_cell_count, uncal_cell_count, cal_node_cell_count, tmp_count = 0, global_uncal_cell_count;
        int     cal_cell_found, c0, node0, c1, node1, c0_nnb_min_dist, c0_nnb, c0_min_dist, index_in_global_list;
        int     *uncal_cell_c0 = NULL, *uncal_cell_node0 = NULL, *uncal_cell_c1 = NULL, *uncal_cell_node1 = NULL, *list_of_int = NULL;
        int     *global_uncal_cell_c0 = NULL, *global_uncal_cell_node0 = NULL, *global_uncal_cell_c1 = NULL, *global_uncal_cell_node1 = NULL;
        double  dist_sq, min_dist_sq, w0;
        double  x_cal[3];
        double  *uncal_cell_dist = NULL, *uncal_cell_w0 = NULL, *cell_length = NULL, *list_of_double = NULL, *list_of_dist = NULL;
        double  *global_uncal_cell_w0 = NULL, *uncal_cell_displacement = NULL, *global_uncal_cell_displacement = NULL, *global_uncal_cell_dist = NULL;

#endif

        // allocate and init local vars
        {
            uncal_cell_c0 =    (int*)malloc(_CELL_DICT(i_layer).number_of_int_cells * sizeof(int));
            uncal_cell_node0 = (int*)malloc(_CELL_DICT(i_layer).number_of_int_cells * sizeof(int));
            uncal_cell_c1 =    (int*)malloc(_CELL_DICT(i_layer).number_of_int_cells * sizeof(int));
            uncal_cell_node1 = (int*)malloc(_CELL_DICT(i_layer).number_of_int_cells * sizeof(int));

            uncal_cell_w0 =   (double*)malloc(_CELL_DICT(i_layer).number_of_int_cells * sizeof(double));
            uncal_cell_dist = (double*)malloc(_CELL_DICT(i_layer).number_of_int_cells * sizeof(double));

            uncal_cell_displacement = (double*)malloc((3 * _CELL_DICT(i_layer).number_of_int_cells) * sizeof(double));

            cal_1st_cell_count = 0;
            cal_next_neigh_cell_count = 0;
            cal_all_cell_count =  0;
            uncal_cell_count =   0;

            rcfd->Tracer->monitor_counter[i_phase] = 0;

            cell_length = (double*)malloc(_CELL_DICT(i_layer).number_of_int_cells * sizeof(double));

            for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_int_cells; i_cell++){

                cell_length[i_cell] = pow((3.0 / (4.0 * 3.14) * _C(i_layer).volume[i_cell]), (1./3.));   // radius of equivalent sphere
            }
        }

        // look for nearest {c0 <> x_cal(c1)} pairing on myid
        for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_int_cells; i_cell++){

            if(_C(i_layer).vof[0][i_cell][i_phase] > rcfd->Solver_Dict->vof_cutoff_limit){

                // init loop vars
                {
                    cal_cell_found = 0;

                    for(i_dim = 0; i_dim < 3; i_dim++){

                        x_cal[i_dim] = _C(i_layer).cfd_displ[i_phase][i_cell][i_dim];
                    }

                    c1 = i_cell;
                }

                // check nnb cells of _C(i_layer).cal_cell_id
                if(!cal_cell_found){

                    c0_nnb_min_dist = c1;

                    min_dist_sq = cell2point_distance_sq(rcfd, i_layer, c0_nnb_min_dist, x_cal);

                    for(i_nnb = 0; i_nnb < _C(i_layer).number_of_next_neigh_cells[c1]; i_nnb++){

                        c0_nnb = _C(i_layer).next_neigh_cell[c1][i_nnb];

                        if(_C(i_layer).vof[0][c0_nnb][i_phase] > rcfd->Solver_Dict->vof_cutoff_limit){

                            dist_sq = cell2point_distance_sq(rcfd, i_layer, c0_nnb, x_cal);

                            if(c0_nnb < _CELL_DICT(i_layer).number_of_int_cells){    // interior cell

                                if(dist_sq < min_dist_sq){

                                    c0_nnb_min_dist = c0_nnb;
                                    min_dist_sq = dist_sq;
                                }
                            }
                        }
                    }

                    if(min_dist_sq < cell_length[c0_nnb_min_dist] * cell_length[c0_nnb_min_dist]){

                        w0 = 1.0;

                        shift_add_to_list(rcfd, i_layer, i_phase, c0_nnb_min_dist, myid, w0, c1, myid);

                        cal_cell_found = 1;

                        cal_next_neigh_cell_count++;
                    }
                }

                // check all cells of myid
                if(!cal_cell_found){

                    min_dist_sq = rcfd->Calibration_Dict->great;

                    c0_min_dist = -1;

                    for(i_cell2 = 0; i_cell2 < _CELL_DICT(i_layer).number_of_int_cells; i_cell2++){

                        if(_C(i_layer).vof[0][i_cell2][i_phase] > rcfd->Solver_Dict->vof_cutoff_limit){

                            dist_sq = cell2point_distance_sq(rcfd, i_layer, i_cell2, x_cal);

                            if(dist_sq < min_dist_sq){

                                min_dist_sq = dist_sq;
                                c0_min_dist = i_cell2;
                            }
                        }
                    }

                    w0 = 1.0;

                    if((c0_min_dist != -1) && (min_dist_sq < cell_length[c0_min_dist] * cell_length[c0_min_dist])){

                        shift_add_to_list(rcfd, i_layer, i_phase, c0_min_dist, myid, w0, c1, myid);

                        cal_cell_found = 1;

                        cal_all_cell_count++;
                    }
                    else{

                        uncal_cell_c0[uncal_cell_count] =     c0_min_dist;
                        uncal_cell_node0[uncal_cell_count] =  myid;
                        uncal_cell_w0[uncal_cell_count] =     w0;
                        uncal_cell_c1[uncal_cell_count] =     c1;
                        uncal_cell_node1[uncal_cell_count] =  myid;

                        uncal_cell_dist[uncal_cell_count] =   sqrt(min_dist_sq);

                        for(i_dim = 0; i_dim < 3; i_dim++){

                            uncal_cell_displacement[(3 * uncal_cell_count + i_dim)] = x_cal[i_dim];
                        }

                        uncal_cell_count++;
                    }
                }
            }
        }

        // look for nearest {c0 <> x_cal(c1)} pairing on other nodes
        {
            // provide global list for all nodes
            {
                if(myid > 0){

                    PRF_CSEND_INT(node_zero, &uncal_cell_count, 1, myid);
                }

                if(myid == 0){

                    global_uncal_cell_count = uncal_cell_count;

                    for(i_node = 1; i_node < compute_node_count; i_node++){

                        PRF_CRECV_INT(i_node, &tmp_count, 1, i_node);

                        global_uncal_cell_count += tmp_count;
                    }

                    for(i_node = 1; i_node < compute_node_count; i_node++){

                        PRF_CSEND_INT(i_node, &global_uncal_cell_count, 1, myid);
                    }
                }

                if(myid > 0){

                    PRF_CRECV_INT(node_zero, &global_uncal_cell_count, 1, node_zero);
                }

                global_uncal_cell_c0 =    (int*)malloc(global_uncal_cell_count * sizeof(int));
                global_uncal_cell_node0 = (int*)malloc(global_uncal_cell_count * sizeof(int));
                global_uncal_cell_c1 =    (int*)malloc(global_uncal_cell_count * sizeof(int));
                global_uncal_cell_node1 = (int*)malloc(global_uncal_cell_count * sizeof(int));

                global_uncal_cell_w0 =   (double*)malloc(global_uncal_cell_count * sizeof(double));
                global_uncal_cell_dist = (double*)malloc(global_uncal_cell_count * sizeof(double));

                global_uncal_cell_displacement =   (double*)malloc((3 * global_uncal_cell_count) * sizeof(double));

                if(myid > 0){

                    PRF_CSEND_INT(node_zero, &uncal_cell_count, 1, myid);

                    PRF_CSEND_INT(node_zero, uncal_cell_c0, uncal_cell_count, myid);
                    PRF_CSEND_INT(node_zero, uncal_cell_c1, uncal_cell_count, myid);
                    PRF_CSEND_INT(node_zero, uncal_cell_node0, uncal_cell_count, myid);
                    PRF_CSEND_INT(node_zero, uncal_cell_node1, uncal_cell_count, myid);
                    PRF_CSEND_REAL(node_zero, uncal_cell_w0, uncal_cell_count, myid);
                    PRF_CSEND_REAL(node_zero, uncal_cell_dist, uncal_cell_count, myid);

                    PRF_CSEND_REAL(node_zero, uncal_cell_displacement, (3 * uncal_cell_count), myid);

                    PRF_CRECV_INT(node_zero, global_uncal_cell_c0, global_uncal_cell_count, node_zero);
                    PRF_CRECV_INT(node_zero, global_uncal_cell_c1, global_uncal_cell_count, node_zero);
                    PRF_CRECV_INT(node_zero, global_uncal_cell_node0, global_uncal_cell_count, node_zero);
                    PRF_CRECV_INT(node_zero, global_uncal_cell_node1, global_uncal_cell_count, node_zero);
                    PRF_CRECV_REAL(node_zero, global_uncal_cell_w0, global_uncal_cell_count, node_zero);
                    PRF_CRECV_REAL(node_zero, global_uncal_cell_dist, global_uncal_cell_count, node_zero);

                    PRF_CRECV_REAL(node_zero, global_uncal_cell_displacement, (3 * global_uncal_cell_count), node_zero);
                }
                else if(myid == 0){

                    for(i_shift = 0; i_shift < uncal_cell_count; i_shift++){

                        global_uncal_cell_c0[i_shift] =    uncal_cell_c0[i_shift];
                        global_uncal_cell_c1[i_shift] =    uncal_cell_c1[i_shift];
                        global_uncal_cell_node0[i_shift] = uncal_cell_node0[i_shift];
                        global_uncal_cell_node1[i_shift] = uncal_cell_node1[i_shift];

                        global_uncal_cell_w0[i_shift] =    uncal_cell_w0[i_shift];
                        global_uncal_cell_dist[i_shift] =  uncal_cell_dist[i_shift];

                        for(i_dim = 0; i_dim < 3; i_dim++){

                            global_uncal_cell_displacement[(3 * i_shift + i_dim)] = uncal_cell_displacement[(3 * i_shift + i_dim)];
                        }
                    }

                    index_in_global_list = i_shift;

                    for(i_node = 1; i_node < compute_node_count; i_node++){

                        PRF_CRECV_INT(i_node, &tmp_count, 1, i_node);

                        if(tmp_count > 0){

                            PRF_CRECV_INT(i_node, &global_uncal_cell_c0[index_in_global_list], tmp_count, i_node);
                            PRF_CRECV_INT(i_node, &global_uncal_cell_c1[index_in_global_list], tmp_count, i_node);
                            PRF_CRECV_INT(i_node, &global_uncal_cell_node0[index_in_global_list], tmp_count, i_node);
                            PRF_CRECV_INT(i_node, &global_uncal_cell_node1[index_in_global_list], tmp_count, i_node);
                            PRF_CRECV_REAL(i_node, &global_uncal_cell_w0[index_in_global_list], tmp_count, i_node);
                            PRF_CRECV_REAL(i_node, &global_uncal_cell_dist[index_in_global_list], tmp_count, i_node);

                            PRF_CRECV_REAL(i_node, &global_uncal_cell_displacement[3 * index_in_global_list], (3 * tmp_count), i_node);

                            index_in_global_list += tmp_count;
                        }
                    }

                    for(i_node = 1; i_node < compute_node_count; i_node++){

                        PRF_CSEND_INT(i_node, global_uncal_cell_c0, global_uncal_cell_count, myid);
                        PRF_CSEND_INT(i_node, global_uncal_cell_c1, global_uncal_cell_count, myid);
                        PRF_CSEND_INT(i_node, global_uncal_cell_node0, global_uncal_cell_count, myid);
                        PRF_CSEND_INT(i_node, global_uncal_cell_node1, global_uncal_cell_count, myid);
                        PRF_CSEND_REAL(i_node, global_uncal_cell_w0, global_uncal_cell_count, myid);
                        PRF_CSEND_REAL(i_node, global_uncal_cell_dist, global_uncal_cell_count, myid);

                        PRF_CSEND_REAL(i_node, global_uncal_cell_displacement, (3 * global_uncal_cell_count), myid);
                    }
                }
            }

            // loop all cells and try to find min_dist
            {
                cal_node_cell_count = 0;

                for(i_shift = 0; i_shift < global_uncal_cell_count; i_shift++){

                    if(global_uncal_cell_node1[i_shift] != myid){   // inter-node shift

                        for(i_dim = 0; i_dim < 3; i_dim++){

                            x_cal[i_dim] = global_uncal_cell_displacement[(3 * i_shift + i_dim)];
                        }

                        if((x_cal[0] > _C(i_layer).x_min_node[0]) && (x_cal[0] < _C(i_layer).x_max_node[0]) &&
                           (x_cal[1] > _C(i_layer).x_min_node[1]) && (x_cal[1] < _C(i_layer).x_max_node[1]) &&
                           (x_cal[2] > _C(i_layer).x_min_node[2]) && (x_cal[2] < _C(i_layer).x_max_node[2])){

                            min_dist_sq = global_uncal_cell_dist[i_shift] * global_uncal_cell_dist[i_shift];

                            c0_min_dist = -1;

                            for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_int_cells; i_cell++){

                                if(_C(i_layer).vof[0][i_cell][i_phase] > rcfd->Solver_Dict->vof_cutoff_limit){

                                    dist_sq = cell2point_distance_sq(rcfd, i_layer, i_cell, x_cal);

                                    if(dist_sq < min_dist_sq){

                                        min_dist_sq = dist_sq;
                                        c0_min_dist = i_cell;
                                    }
                                }
                            }

                            if(c0_min_dist != -1){

                                global_uncal_cell_dist[i_shift] =  sqrt(min_dist_sq);
                                global_uncal_cell_c0[i_shift] =    c0_min_dist;
                                global_uncal_cell_node0[i_shift] = myid;
                                //global_uncal_cell_w0[i_shift] =    _C(i_layer).volume[c0_min_dist];
                                global_uncal_cell_w0[i_shift] = 1.0;

                                cal_node_cell_count++;
                            }
                        }
                    }
                }
            }

            // provide adapted list to all nodes
            {
                if(myid == 0){

                    list_of_dist =   (double*)malloc(global_uncal_cell_count * sizeof(double));
                    list_of_double = (double*)malloc(global_uncal_cell_count * sizeof(double)); // w0
                    list_of_int =    (int*)malloc(global_uncal_cell_count * sizeof(int));       // c0

                    for(i_node = 1; i_node < compute_node_count; i_node++){

                        PRF_CRECV_REAL(i_node, list_of_dist, global_uncal_cell_count, i_node);
                        PRF_CRECV_INT(i_node, list_of_int, global_uncal_cell_count, i_node);
                        PRF_CRECV_REAL(i_node, list_of_double, global_uncal_cell_count, i_node);

                        for(i_shift = 0; i_shift < global_uncal_cell_count; i_shift++){

                            if(list_of_dist[i_shift] < global_uncal_cell_dist[i_shift]){

                                global_uncal_cell_c0[i_shift] = list_of_int[i_shift];
                                global_uncal_cell_node0[i_shift] = i_node;
                                global_uncal_cell_w0[i_shift] = list_of_double[i_shift];
                                global_uncal_cell_dist[i_shift] = list_of_dist[i_shift];
                            }
                        }
                    }

                    free(list_of_dist);
                    free(list_of_double);
                    free(list_of_int);

                    for(i_node = 1; i_node < compute_node_count; i_node++){

                        PRF_CSEND_INT(i_node,  global_uncal_cell_c0,    global_uncal_cell_count, myid);
                        PRF_CSEND_INT(i_node,  global_uncal_cell_node0, global_uncal_cell_count, myid);
                        PRF_CSEND_REAL(i_node, global_uncal_cell_w0,    global_uncal_cell_count, myid);
                        PRF_CSEND_INT(i_node,  global_uncal_cell_c1,    global_uncal_cell_count, myid);
                        PRF_CSEND_INT(i_node,  global_uncal_cell_node1, global_uncal_cell_count, myid);
                    }
                }

                if(myid> 0){

                    PRF_CSEND_REAL(node_zero, global_uncal_cell_dist, global_uncal_cell_count, myid);
                    PRF_CSEND_INT(node_zero,  global_uncal_cell_c0,   global_uncal_cell_count, myid);
                    PRF_CSEND_REAL(node_zero, global_uncal_cell_w0,   global_uncal_cell_count, myid);

                    PRF_CRECV_INT(node_zero,  global_uncal_cell_c0,    global_uncal_cell_count, node_zero);
                    PRF_CRECV_INT(node_zero,  global_uncal_cell_node0, global_uncal_cell_count, node_zero);
                    PRF_CRECV_REAL(node_zero, global_uncal_cell_w0,    global_uncal_cell_count, node_zero);
                    PRF_CRECV_INT(node_zero,  global_uncal_cell_c1,    global_uncal_cell_count, node_zero);
                    PRF_CRECV_INT(node_zero,  global_uncal_cell_node1, global_uncal_cell_count, node_zero);
                }
            }

            // add shifts
            {
                for(i_shift = 0; i_shift < global_uncal_cell_count; i_shift++){

                    if(global_uncal_cell_node1[i_shift] == myid){

                        c0 =    global_uncal_cell_c0[i_shift];
                        node0 = global_uncal_cell_node0[i_shift];
                        w0 =    global_uncal_cell_w0[i_shift];
                        c1 =    global_uncal_cell_c1[i_shift];
                        node1 = myid;

                        shift_add_to_list(rcfd, i_layer, i_phase, c0, node0, w0, c1, node1);
                    }
                }
            }
        }

        // free local vars
        {
            free(uncal_cell_c0);
            free(uncal_cell_node0);
            free(uncal_cell_w0);
            free(uncal_cell_c1);
            free(uncal_cell_node1);

            free(uncal_cell_dist);
            free(uncal_cell_displacement);

            free(cell_length);

            free(global_uncal_cell_c0);
            free(global_uncal_cell_node0);
            free(global_uncal_cell_w0);
            free(global_uncal_cell_c1);
            free(global_uncal_cell_node1);

            free(global_uncal_cell_dist);
            free(global_uncal_cell_displacement);
        }

        if(LOCAL_DEBUG){

            report_debug(rcfd, RCFD_REPORT_CATEGORY_CREATE_SHIFT,
                         "[%d] shift_from_displacement_fields: 1st_cell %d next_neigh_cell %d all_cell %d no_cell %d node_cell %d",
                         myid, cal_1st_cell_count, cal_next_neigh_cell_count, cal_all_cell_count, uncal_cell_count, cal_node_cell_count);
        }

        return rcfd->Tracer->monitor_counter[i_phase];

#undef LOCAL_DEBUG

#else

        return -1;
#endif
    }

/*****************************************************************************/
/************************** start shift_calibrate() **************************/
/*****************************************************************************/
#if RP_NODE
/**
 * @brief Find the closest neighbor of c0 to a point.
 *
 * @param rcfd          the rcfd struct
 * @param c0            the starting point of the search for the closest cell
 * @param x_cal         the point for which the closest cell should be found.
 * @param i_layer       the current layer
 * @param i_phase       the current phase (necessary for vof_cutoff_limit)
 * @param inclued_c0    flag sets if c0 should also be considered to obtain the closest cell to x_cal
 * @return the cell id of the closest neighbor to x_cal
 */
int find_closest_neighbor(rCFD_type *rcfd, int c0, const double *x_cal, int i_layer, int i_phase, int include_c0)
{
    double dist_sq, min_dist_sq;
    int i_nnb, c0_nnb, min_dist_id;

    // initialize the minimal distance with c0 or the first neighbor
    if(include_c0){

        min_dist_sq = cell2point_distance_sq(rcfd, i_layer, c0, x_cal);
        min_dist_id = c0;
        i_nnb = 0;
    }
    else{

        min_dist_sq = cell2point_distance_sq(rcfd, i_layer, _C(i_layer).next_neigh_cell[c0][0], x_cal);
        min_dist_id = _C(i_layer).next_neigh_cell[c0][0];
        i_nnb = 1;
    }

    // loop over all neighbors of c0 to find the one closest to x_cal
    for(; i_nnb < _C(i_layer).number_of_next_neigh_cells[c0]; i_nnb++){

        c0_nnb = _C(i_layer).next_neigh_cell[c0][i_nnb];

        if(_C(i_layer).vof[0][c0_nnb][i_phase] > rcfd->Solver_Dict->vof_cutoff_limit){

            dist_sq = cell2point_distance_sq(rcfd, i_layer, c0_nnb, x_cal);

            if(dist_sq < min_dist_sq){

                min_dist_sq = dist_sq;
                min_dist_id = c0_nnb;
            }
        }
    }

    return min_dist_id;
}

/**
 * @brief Find the closest cell to a point by iteratively looking through
 *        c0 and its neighboring cells.
 *
 * @param rcfd      the rcfd struct
 * @param c0        the starting cell of the search for the closest cell
 * @param x_cal     the point for which the closest cell should be found.
 * @param i_layer   the current layer
 * @param i_phase   the current phase (necessary for vof_cutoff_limit)
 * @return the cell id of the closest cell to x_cal
 */
int find_closest_cell(rCFD_type *rcfd, int c0, const double *x_cal, int i_layer, int i_phase)
{
    int min_dist_id;
    int c0_old = -1;

    // as long as a a neighboring cell has been found that is closer, check the neighbors of the
    // new closest cell for even closer neighbors.
    while(c0 != c0_old){

        c0_old = c0;

        min_dist_id = find_closest_neighbor(rcfd, c0, x_cal, i_layer, i_phase, 1);

        // if interior cell, adapt c0
        if(min_dist_id < _CELL_DICT(i_layer).number_of_int_cells){

            c0 = min_dist_id;
        }
    }

    return c0;
}

/**
 * @brief Finds all neighbors that could possibly used for interpolation of vec_c0_cal.
 *        It does so by checking if the scalar product of vec_c0_cal and the calclulated
 *        C2C vectors is positive.
 *
 * @param c0                    origin cell from which the neighbors are iterated
 * @param vec_c0_cal            the vector that should be interpolated via C2C vectors
 * @param cell_centers          used to calculate C2C vectors
 * @param next_neigh_cells      The neighboring cells to be iterated.
 * @param len_next_neigh_cells  length of next_neigh_cells
 * @param aligned_neigh_list    list to be filled with aligned neighbors
 * @param len_aligned_neigh_list length of aligned_neigh_list
 * @param cell_vofs             contains vof of all cells (necessary for vof_cutoff)
 * @param vof_cutoff            used to ignore shifts from cells with too little vof's
 * @param i_phase               the current phase (necessary for vof_cutoff)
 * @return number of aligned neighbor cells or -1 if there were more than len_aligned_neigh_list
 */
int find_aligned_neighbors(int c0, const double *vec_c0_cal, double **cell_centers,
                           const int *next_neigh_cells, int len_next_neigh_cells,
                           int *aligned_neigh_list, int len_aligned_neigh_list,
                           double **cell_vofs, int vof_cutoff,
                           int i_phase)
{
    int n_align_neigh_cells = 0;
    int i_nb, i_nnb, i_dim, c0_nb;

    double vec_c0_nb[3];
    double vec_align;

    // init aligned_neigh_list with -1 == no aligned neighbors
    for(i_nb = 0; i_nb < len_aligned_neigh_list; i_nb++){

        aligned_neigh_list[i_nb] = -1;
    }

    // iterate through all neighbors of c0...
    for(i_nnb = 0; i_nnb < len_next_neigh_cells; i_nnb++){

        c0_nb = next_neigh_cells[i_nnb];

        if(cell_vofs[c0_nb][i_phase] > vof_cutoff){

            // ...calculate the C2C vector from c0 to c0_nb via the cell-centers...
            for(i_dim = 0; i_dim < 3; i_dim++){

                vec_c0_nb[i_dim] = cell_centers[c0_nb][i_dim] - cell_centers[c0][i_dim];
            }

            // ... and check if the scalar projection on vec_c0_cal is positive.
            // (if it is not, you immediately know that the neighbor cell can only
            // extrapolate vec_c0_cal and not interpolate as it is pointing away
            // from x_cal.)
            vec_align = scalar_projection(vec_c0_nb, vec_c0_cal);

            if(vec_align > 0.0){

                // if it is, add it to the list of possible aligned cells.
                if(n_align_neigh_cells < len_aligned_neigh_list){

                    aligned_neigh_list[n_align_neigh_cells] = c0_nb;

                    n_align_neigh_cells++;
                }
                else{

                    return -1;
                }
            }
        }
    }

    return n_align_neigh_cells;
}

/**
 * @brief Calculate the C2C shift vectors between new_shift[0] and
 *        new_shift[1,2,3] and writes them to M
 *
 * @param new_shifts    The list of shifts from which the C2C shift vectors are calculated.
 * @param M             The shift matrix to be set up.
 * @param cell_centers The cell positions for the calculation of the C2C shift vectors.
 */
void setup_shift_matrix_3d(const C2C_shift_type *new_shifts, double **M, double **cell_centers)
{
    int i_dim;

    for(i_dim = 0; i_dim < 3; i_dim++){

        M[i_dim][0] = cell_centers[new_shifts[1].c0][i_dim] - cell_centers[new_shifts[0].c0][i_dim];
        M[i_dim][1] = cell_centers[new_shifts[2].c0][i_dim] - cell_centers[new_shifts[0].c0][i_dim];
        M[i_dim][2] = cell_centers[new_shifts[3].c0][i_dim] - cell_centers[new_shifts[0].c0][i_dim];
    }
}

/**
 * @brief Calculate the C2C shift vectors between new_shift[0] and
 *        new_shift[1,2] and the cross product of new_shift[1,2] and writes them to M.
 *
 * @param new_shifts    The list of shifts from which the C2C shift vectors are calculated.
 * @param M             The shift matrix to be set up.
 * @param cell_centers The cell positions for the calculation of the C2C shift vectors.
 */
void setup_shift_matrix_2d(const C2C_shift_type *new_shifts, double **M, double **cell_centers)
{
    int i_dim;
    double a[3], b[3], axb[3];

    for(i_dim = 0; i_dim < 3; i_dim++){

        a[i_dim] = cell_centers[new_shifts[1].c0][i_dim] - cell_centers[new_shifts[0].c0][i_dim];
        b[i_dim] = cell_centers[new_shifts[2].c0][i_dim] - cell_centers[new_shifts[0].c0][i_dim];
    }
    vector_cross_product(a, b, axb);

    for(i_dim = 0; i_dim < 3; i_dim++){

        M[i_dim][0] =   a[i_dim];
        M[i_dim][1] =   b[i_dim];
        M[i_dim][2] = axb[i_dim];
    }
}

/**
 * @brief Advance a set of 3 indices so that all possible combinations of
 *        indices will be reached for a given length of a list.
 *
 * @param indices   the indices to be advanced.
 * @param len_list  the length of the list the indices point at.
 */
void advance_indices(int *indices, int len_list)
{
    if(indices[2] < (len_list - 2)){

        indices[2]++;
    }
    else if(indices[1] < (len_list - 3)){

        indices[1]++;
        indices[2] = indices[1] + 1;
    }
    else{

        indices[0]++;
        indices[1] = indices[0] + 1;
        indices[2] = indices[0] + 2;
    }
}

/**
 * @brief Try to interpolate the vector vec_c0_cal by using a combination of cell centers.
 *
 * @param vec_c0_cal             The vector to be interpolated.
 * @param new_shifts             The C2C shifts that should be filled with correct shifts and weights
 * @param aligned_neigh_list     The list of cells that can be used for the interpolation.
 * @param len_aligned_neigh_list The length of aligned_neigh_list.
 * @param cell_centers           The cell positions for the calculation of the C2C shift vectors.
 * @return 1 if interpolation is successfull and 0 if it fails.
 */
int interpolate_vec_c0_cal(const double *vec_c0_cal, C2C_shift_type *new_shifts,
                           const int *aligned_neigh_list, int len_aligned_neigh_list,
                           double **cell_centers)
{
    // allocate and init M, M_inv
    int nbytes = 3 * 3 * sizeof(double);

    int i_neighbors[3]= {0,1,2}; // indices to iterate through all combinations of neighbors
    double det;

    double **M       = malloc_r_2d(3, 3);
    double **M_inv   = malloc_r_2d(3, 3);

    memset(M[0],    0, nbytes);
    memset(M_inv[0],0, nbytes);

    // loop over indices of aligned nb cells
    while(i_neighbors[2] < (len_aligned_neigh_list - 3)){

        // set vectors to triple nb cells
        new_shifts[1].c0 = aligned_neigh_list[i_neighbors[0]];
        new_shifts[2].c0 = aligned_neigh_list[i_neighbors[1]];
        new_shifts[3].c0 = aligned_neigh_list[i_neighbors[2]];

        setup_shift_matrix_3d(new_shifts, M, cell_centers);

        // calculate individual weights
        det = inverse_3x3(M, M_inv);

        if(det != 0.0){

            new_shifts[1].w0 = M_inv[0][0]*vec_c0_cal[0] + M_inv[0][1]*vec_c0_cal[1] + M_inv[0][2]*vec_c0_cal[2];
            new_shifts[2].w0 = M_inv[1][0]*vec_c0_cal[0] + M_inv[1][1]*vec_c0_cal[1] + M_inv[1][2]*vec_c0_cal[2];
            new_shifts[3].w0 = M_inv[2][0]*vec_c0_cal[0] + M_inv[2][1]*vec_c0_cal[1] + M_inv[2][2]*vec_c0_cal[2];

            new_shifts[0].w0 = 1.0 - (new_shifts[1].w0 + new_shifts[2].w0 + new_shifts[3].w0);
        }
        else{

            new_shifts[0].w0 = -1.0;
            new_shifts[1].w0 = -1.0;
            new_shifts[2].w0 = -1.0;
            new_shifts[3].w0 = -1.0;
        }

        // check if all pos - return sucessfully if true
        if( (new_shifts[0].w0 > 0.0) &&
            (new_shifts[1].w0 > 0.0) &&
            (new_shifts[2].w0 > 0.0) &&
            (new_shifts[3].w0 > 0.0)){

            free_r_2d(M);
            free_r_2d(M_inv);

            return 1;
        }

        // advance indices i_neighbors
        advance_indices(i_neighbors, len_aligned_neigh_list);
    }

    // free M, M_inv and return fail
    free_r_2d(M);
    free_r_2d(M_inv);

    new_shifts[0].w0 = -1.0;
    new_shifts[1].w0 = -1.0;
    new_shifts[2].w0 = -1.0;
    new_shifts[3].w0 = -1.0;

    return 0;
}

double calc_projection_accuracy(const double *vec_c0_cal,
                                    const C2C_shift_type *new_shifts, int len_new_shifts,
                                    double **cell_centers)
{
    int i_dim, i_shift;
    double displ_pos[3] = {cell_centers[new_shifts[0].c0][0] + vec_c0_cal[0],
                           cell_centers[new_shifts[0].c0][1] + vec_c0_cal[1],
                           cell_centers[new_shifts[0].c0][2] + vec_c0_cal[2]};
    double proj_pos[3] = {0.0, 0.0, 0.0};
    double weight;

    //assume the weights sum up to 1, otherwise we would have to divide by weight
    for(i_shift = 0; i_shift < len_new_shifts; i_shift++){

        for(i_dim = 0; i_dim < 3; i_dim++){

            proj_pos[i_dim] += cell_centers[new_shifts[i_shift].c0][i_dim] * new_shifts[i_shift].w0;
        }
    }

    weight =  (displ_pos[0] - proj_pos[0])*(displ_pos[0] - proj_pos[0])
            + (displ_pos[1] - proj_pos[1])*(displ_pos[1] - proj_pos[1])
            + (displ_pos[2] - proj_pos[2])*(displ_pos[2] - proj_pos[2]);

    return weight;
}

/**
* @brief Projects vec_c0_cal onto a 2D plane, 1D line or 0D point - whichever
 *        gives the shortest projection distance. Never Fails, because you can
 *        always project onto a 0D point.
 *
 * @param vec_c0_cal             The vector to be interpolated.
 * @param min_new_shifts         The C2C shifts that should be filled with correct shifts and weights.
 * @param aligned_neigh_list     The list of cells that can be used for the interpolation.
 * @param len_aligned_neigh_list The length of aligned_neigh_list.
 * @param cell_centers           The cell positions for the calculation of the C2C shift vectors.
 */
void project_vec_c0_cal(const double *vec_c0_cal, C2C_shift_type *min_new_shifts,
                        const int *aligned_neigh_list, int len_aligned_neigh_list,
                        double **cell_centers)
{
    int nbytes = 3 * 3 * sizeof(double);

    int i_neighbors[3]= {0,1,len_aligned_neigh_list}; // indices to iterate through all tuple-combinations of neighbors
    double det;

    double **M       = malloc_r_2d(3, 3);
    double **M_inv   = malloc_r_2d(3, 3);

    memset(M[0],    0, nbytes);
    memset(M_inv[0],0, nbytes);

    double min_proj_dist = -1.0, current_proj_dist;

    C2C_shift_type new_shifts[3] = {min_new_shifts[0], min_new_shifts[1], min_new_shifts[2]};

    // FIRST: Try 2D-Interpolation (only use i_neighbors[0,1] and their cross-product for interpolation)
    // loop over indices of aligned nb cells

    while(i_neighbors[2] < (len_aligned_neigh_list + 1) && len_aligned_neigh_list > 1){

        // set vectors to tuple nb cells
        new_shifts[1].c0 = aligned_neigh_list[i_neighbors[0]];
        new_shifts[2].c0 = aligned_neigh_list[i_neighbors[1]];

        setup_shift_matrix_2d(new_shifts, M, cell_centers);

        // calculate individual weights
        det = inverse_3x3(M, M_inv);

        if(det != 0.0){

            new_shifts[1].w0 = M_inv[0][0]*vec_c0_cal[0] + M_inv[0][1]*vec_c0_cal[1] + M_inv[0][2]*vec_c0_cal[2];
            new_shifts[2].w0 = M_inv[1][0]*vec_c0_cal[0] + M_inv[1][1]*vec_c0_cal[1] + M_inv[1][2]*vec_c0_cal[2];

            new_shifts[0].w0 = 1.0 - (new_shifts[1].w0 + new_shifts[2].w0);
        }
        else{

            new_shifts[0].w0 = -1.0;
            new_shifts[1].w0 = -1.0;
            new_shifts[2].w0 = -1.0;
        }

        // check if all weights pos and if projection distance is smaller - if it is:
        // reset min_new_shifts to better projection
        if( (new_shifts[0].w0 > 0.0) &&
            (new_shifts[1].w0 > 0.0) &&
            (new_shifts[2].w0 > 0.0)){

            current_proj_dist = calc_projection_accuracy(vec_c0_cal, new_shifts, 3, cell_centers);

            if(current_proj_dist < min_proj_dist || min_proj_dist < 0.0){

                min_proj_dist = current_proj_dist;
                min_new_shifts[0] = new_shifts[0];
                min_new_shifts[1] = new_shifts[1];
                min_new_shifts[2] = new_shifts[2];
                //return;
            }
        }

        // advance indices i_neighbors
        advance_indices(i_neighbors, 0);
    }

    // SECOND: Try 1D Interpolation (by using the closest cell and i_neighbors[0])

    int i_dim;
    // vector between cell centers of closest cell and i_neighbors[0]
    double c2c_vector[3], c2c_vec_length;

    new_shifts[2].w0 = -1.0;

    for(i_neighbors[0] = 0; i_neighbors[0] < len_aligned_neigh_list; i_neighbors[0]++){

        new_shifts[1].c0 = aligned_neigh_list[i_neighbors[0]];
        c2c_vec_length = 0.0;

        //calc c2c_vector and it's length
        for(i_dim = 0; i_dim < 3; i_dim++){
            c2c_vector[i_dim] =   cell_centers[new_shifts[1].c0][i_dim]
                                - cell_centers[new_shifts[0].c0][i_dim];
            c2c_vec_length += c2c_vector[i_dim] * c2c_vector[i_dim];
        }
        c2c_vec_length = sqrt(c2c_vec_length);

        // weight is calculated by:
        // c2c_vector . vec_c0_cal / |c2c_vector|^2 = cos(theta) * |vec_c0_cal| / |c2c_vector|
        // which can be seen as the length of vec_c0_cal when you project it onto c2c_vector,
        // diviced by the length of c2c_vector itself.
        // In other words: It is the length ratio between the projected vec_c0_cal and c2c_vector.
        new_shifts[1].w0 = scalar_projection(vec_c0_cal, c2c_vector) / c2c_vec_length;

        // check if all weights pos and if projection distance is smaller - if it is:
        // reset min_new_shifts to better projection
        if(new_shifts[1].w0 <= 1.0 && new_shifts[1].w0 >= 0.0){

            new_shifts[0].w0 = 1.0 - new_shifts[1].w0;

            current_proj_dist = calc_projection_accuracy(vec_c0_cal, new_shifts, 2, cell_centers);

            if(current_proj_dist < min_proj_dist || min_proj_dist < 0.0){

                min_proj_dist = current_proj_dist;
                min_new_shifts[0] = new_shifts[0];
                min_new_shifts[1] = new_shifts[1];
                min_new_shifts[2] = new_shifts[2];
            }
        }
    }

    new_shifts[0].w0 = 1.0;
    current_proj_dist = calc_projection_accuracy(vec_c0_cal, new_shifts, 1, cell_centers);

    // if 2D and 1D interpolation failed, fall back to 0D interpolation
    // by using closest cell_center
    if(current_proj_dist < min_proj_dist || min_proj_dist < 0.0){

        min_new_shifts[0].w0 = 1.0;
        min_new_shifts[1].w0 = -1.0;
        min_new_shifts[2].w0 = -1.0;
    }

    // free M, M_inv
    free_r_2d(M);
    free_r_2d(M_inv);
}

/**
 * @brief Multiply the weights of the given shifts with the vof of the
 *        corresponding c0 cell and re-normalize the resulting weigths
 *
 * @param new_shifts        list of C2C shifts to be re-weighted
 * @param len_new_shifts    length of new_shifts
 * @param vofs              list of vofs to be used as factors for the new weight.
 * @param i_phase           the phase for which the vof should be used.
 */
void shift_multiply_weights_with_vof(C2C_shift_type *new_shifts, int len_new_shifts, double **vofs, int i_phase)
{
    int i_shift;
    double sum_of_weights = 0.0;

    for(i_shift=0; i_shift < len_new_shifts; i_shift++){

        if(new_shifts[i_shift].w0 > 0.0){

            new_shifts[i_shift].w0 *= vofs[new_shifts[i_shift].c0][i_phase];
            sum_of_weights += new_shifts[i_shift].w0;
        }
    }

    if(sum_of_weights > 0.0){

        double inv_sum_of_weights = 1.0/sum_of_weights;

        for(i_shift=0; i_shift < len_new_shifts; i_shift++){

            if(new_shifts[i_shift].w0 > 0.0){

                new_shifts[i_shift].w0 *= inv_sum_of_weights;
            }
        }
    }
}

/**
 * @brief Copy the content of one shift to another shift
 *
 * @param shift     C2C_shift_type to be overwritten
 * @param new_shift C2C_shift_type to be copied
 */
void shiftcpy(C2C_shift_type *destination, const C2C_shift_type *source)
{
    // shallow copy
    *destination = *source;
}

/**
 * @brief prune new_shifts and adds them to the list of shifts. Overwrites
 *        prim_shift and shifts.
 *
 * @param prim_shift        pointer to the primary shift to be replaced by a new shift
 * @param shifts            list of shifts where the new shifts should be added
 * @param new_shifts        list of new shifts that should be added
 * @param new_shifts_max    maximum number of shifts that should be added
 * @param new_shifts_len    length of new shifts
 * @param pruning_on        flag sets wether new shifts should be pruned
 * @param min_shifts_after_pruning sets how many shifts should remain at the very least
 * @param pruning_weight    sets the minimal weight fraction for pruning shifts
 * @return either returns the number of new shifts or the negative number of shifts that could not be added anymore
 */
int prune_add_shifts(C2C_shift_type *prim_shift, C2C_shift_type *shifts,
                     C2C_shift_type *new_shifts, int new_shifts_max, int new_shifts_len,
                     int pruning_on, int min_shifts_after_pruning, double pruning_weight)
{
    int new_shifts_count = 0;
    double sum_of_weights = 0.0;
    double *weights = (double*)malloc(new_shifts_len * sizeof(double));

    int *sorted_ids;
    int i;

    // calc sum_of_weights for comparison with pruning_weight and weights for argsort
    for(i = 0; i < new_shifts_len; i++){

        if(new_shifts[i].w0 > 0.0){

            sum_of_weights += new_shifts[i].w0;
        }

        weights[i] = new_shifts[i].w0;
    }

    //sort shifts with their weights.
    sorted_ids = selectionArgsort(weights, new_shifts_len, 0);

    free(weights);

    for(i = 0; i < new_shifts_len; i++){

        // break if all valid shifts have been added
        if(new_shifts[sorted_ids[i]].w0 < 0.0)
            break;

        // break if pruning is on and all desired shifts have been added
        if( pruning_on &&
            (i >= min_shifts_after_pruning) &&
            (new_shifts[sorted_ids[i]].w0 / sum_of_weights) < pruning_weight)
            break;

        // if no more shifts can be added, returns the negative number of missing shifts.
        if(new_shifts_count >= new_shifts_max){

            free(sorted_ids);

            if(i == 0)
                return new_shifts_count - new_shifts_len;
            else
                return new_shifts_count - new_shifts_len + 1;
        }

        // add new_shifts[sorted_ids[i]] either at the primary shift location or the new_shifts location
        if(i == 0){

            shiftcpy(prim_shift, &new_shifts[sorted_ids[i]]);
        }
        else{

            shiftcpy(&shifts[new_shifts_count], &new_shifts[sorted_ids[i]]);
            new_shifts_count++;
        }
    }

    free(sorted_ids);

    return new_shifts_count;
}

/**
 * @brief Count the number of valid shifts (with a weight larger than zero)
 *        contained in the given list of shifts.
 *
 * @param new_shifts list of C2C shifts
 * @param len_new_shifts length of new_shifts
 * @return Number of valid shifts (-1 as one shift is old)
 */
int valid_new_shifts(C2C_shift_type *new_shifts, int len_new_shifts)
{
    int i, valid_shifts = 0;

    for(i = 0; i < len_new_shifts; i++){

        if(new_shifts[i].w0 > 0.0){

            valid_shifts++;
        }
    }

    // one shift was old
    return valid_shifts - 1;
}
#endif

int shift_calibrate(rCFD_type *rcfd, int i_layer, int i_phase, C2C_shift_type *shifts, double **displ, int number_of_shifts)
{
#if RP_NODE
    /* known limitations
        if format == RCFD_INTERP_NEAREST_ALIGNED_NEIGHBOR_CELLS:
            in some cases three c0-nb cells are not enough to reach cal_x (we might need to use more nb cells)
    */
#define LOCAL_DEBUG 0

    int     i_shift, i_shift_new, i_dim, n_align_neigh_cells;
    int     success_interpolation, unaligned_count, aligned_wo_success_count, aligned_w_success_count, pruning_count;
    int     *aligned_neigh_list = NULL;
    double  vec_c0_cal[3];
    int     len_new_shifts;
    C2C_shift_type *new_shifts = NULL;

    switch(rcfd->Calibration_Dict->format){

        case RCFD_INTERP_NEAREST_ALIGNED_NEIGHBOR_CELLS:
        {
            // init global counters, allocate aligned_neigh_list
            i_shift_new = 0;

            aligned_neigh_list = (int*)malloc(rcfd->Calibration_Dict->max_number_of_aligned_neigh_cells * sizeof(int));

            unaligned_count = 0;
            aligned_wo_success_count = 0;
            aligned_w_success_count = 0;
            pruning_count = 0;

            // the interpolation method will take up to 4 shifts
            len_new_shifts = 4;
            new_shifts = (C2C_shift_type*)malloc(len_new_shifts * sizeof(C2C_shift_type));

            // loop primary shifts and add new shifts by calibration
            for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                // set loop vars
                {
                    int i;
                    C2C_shift_type init_shift = {.c0 = -1,
                                                 .c1 = shifts[i_shift].c1,
                                                 .node0 = shifts[i_shift].node0,
                                                 .node1 = shifts[i_shift].node1,
                                                 .w0 = -1};

                    // init new_shifts - first shift is primary shift
                    shiftcpy(&new_shifts[0], &shifts[i_shift]);

                    for(i = 1; i < len_new_shifts; i++){

                        shiftcpy(&new_shifts[i], &init_shift);
                    }

                    success_interpolation = 0;
                }

                new_shifts[0].c0 = find_closest_cell(rcfd, new_shifts[0].c0, displ[i_shift], i_layer, i_phase);

                // vec x_c0 to displ[i_shift]
                for(i_dim = 0; i_dim < 3; i_dim++){

                    vec_c0_cal[i_dim] = displ[i_shift][i_dim] - _C(i_layer).x[new_shifts[0].c0][i_dim];
                }

                // create list of aligned nb cells
                // aligned cells are cells that can be used to interpolate displ[i_shift],
                // unaligned cells can only extrapolate displ[i_shift].
                n_align_neigh_cells = find_aligned_neighbors(new_shifts[0].c0,
                                                             vec_c0_cal,
                                                             _C(i_layer).x,
                                                             _C(i_layer).next_neigh_cell[new_shifts[0].c0],
                                                             _C(i_layer).number_of_next_neigh_cells[new_shifts[0].c0],
                                                             aligned_neigh_list,
                                                             rcfd->Calibration_Dict->max_number_of_aligned_neigh_cells,
                                                             _C(i_layer).vof[0],
                                                             rcfd->Solver_Dict->vof_cutoff_limit,
                                                             i_phase);

                // if there were more aligned neighbor cells than what fit inside the list, send a warning...
                if(n_align_neigh_cells == -1){

                    report_warning(rcfd, RCFD_REPORT_CATEGORY_CREATE_SHIFT,
                                   "[%d] shift_calibrate: new_shifts[0].c0 = %d has more than %d aligned neighbours",
                                   myid, new_shifts[0].c0, rcfd->Calibration_Dict->max_number_of_aligned_neigh_cells);

                    n_align_neigh_cells = rcfd->Calibration_Dict->max_number_of_aligned_neigh_cells;
                }

                // if three or more aligned cells have been found, try to interpolate.
                if(n_align_neigh_cells >= 3){

                    success_interpolation = interpolate_vec_c0_cal(vec_c0_cal, new_shifts, aligned_neigh_list, n_align_neigh_cells, _C(i_layer).x);

                    // update debug information:
                    if(success_interpolation)
                        aligned_w_success_count++;
                    else
                        aligned_wo_success_count++;
                }
                else{

                    unaligned_count++;
                }

                if(!success_interpolation){ // 2D or 1D projection - whichever gives the smallest projection error

                    project_vec_c0_cal(vec_c0_cal, new_shifts, aligned_neigh_list, n_align_neigh_cells, _C(i_layer).x);
                }

                // adapt the weights of all valid shifts with their vof's
                if((rcfd->Solver_Dict->number_of_phases > 1) && (rcfd->Calibration_Dict->multiply_shift_weights_with_vof_on)){

                    shift_multiply_weights_with_vof(new_shifts, len_new_shifts, _C(i_layer).vof[0], i_phase);
                }

                // prune and add shifts
                int added_shifts = prune_add_shifts(&shifts[i_shift],
                                                    &shifts[number_of_shifts + i_shift_new],
                                                    new_shifts,
                                                    (3 * number_of_shifts - i_shift_new),
                                                    len_new_shifts,
                                                    rcfd->Calibration_Dict->calibration_pruning_on,
                                                    rcfd->Calibration_Dict->calibration_min_shifts_after_pruning,
                                                    rcfd->Calibration_Dict->calibration_pruning_weight);

                pruning_count += valid_new_shifts(new_shifts, len_new_shifts) - added_shifts;

                if(added_shifts < 0){

                    i_shift_new = 3 * number_of_shifts; // set to max number of new shifts
                    report_warning(rcfd, RCFD_REPORT_CATEGORY_CREATE_SHIFT,
                                   "[%d] shift_calibrate: failed to add new shifts (list is full)", myid);
                }
                else{

                    i_shift_new += added_shifts;
                }
            }

            // free local var
            free(new_shifts);
            free(aligned_neigh_list);

            if(LOCAL_DEBUG){

                report_debug(rcfd, RCFD_REPORT_CATEGORY_CREATE_SHIFT,
                             "[%d] shift_calibrate: total new shifts %d "
                             "( = unaligned %d + aligned_w/o_sucess %d + 3 * aligned_success %d - pruning_count %d)",
                             myid, i_shift_new,
                             unaligned_count, aligned_wo_success_count, aligned_w_success_count, pruning_count);
            }

            return i_shift_new;
        }
    }
#undef LOCAL_DEBUG

#endif // RP_NODE

    return 0;
}

    int shift_write(rCFD_type* rcfd, int i_state, int i_frame, int i_phase, int i_layer)
    {
#if RP_NODE

#define LOCAL_DEBUG 0

#if 1   /* local vars */
        int     i_shift;
        int     total_number_of_stored_C2Cs;
        int     c0, c1, node0, node1;
        double  w0;
        FILE    *f_out = NULL;
        char    file_name[80];
#endif

        /* open f_out */
        {
            sprintf(file_name, "%s_%d_%d_%d", rcfd->File_Dict->C2C_filename, i_state, i_phase, myid);

            if(rcfd->Tracer->frame_counter == 0){

                f_out = fopen(file_name, "w");
            }
            else{

                f_out = fopen(file_name, "a");
            }

            if(f_out == NULL){

                report_error(rcfd, RCFD_REPORT_CATEGORY_WRITE_SHIFT, "[%d] shift_write: f_out == NULL", myid);
                return -1;
            }
        }

        /* write C2C shifts */
        {
            if(rcfd->Tracer->monitor_counter[i_phase] > 0){

                /* write valid c2c shifts, update _C(i_layer).weight_after_shift */
                {
                    //fprintf(f_out,"%d \n", (rcfd->Tracer->monitor_counter[i_phase] - number_of_lock_cells));
                    fprintf(f_out,"%d \n", rcfd->Tracer->monitor_counter[i_phase]);

                    for(i_shift = 0; i_shift < rcfd->Tracer->monitor_counter[i_phase]; i_shift++){

                        c0 =    rcfd->Tracer->shifts[i_phase][i_shift].c0;
                        node0 = rcfd->Tracer->shifts[i_phase][i_shift].node0;
                        c1 =    rcfd->Tracer->shifts[i_phase][i_shift].c1;
                        node1 = rcfd->Tracer->shifts[i_phase][i_shift].node1;
                        w0 =    rcfd->Tracer->shifts[i_phase][i_shift].w0;

                        fprintf(f_out,"%d %d %d %d %e\n", c0, node0, c1, node1, w0);

#if 0   // old code

                        if(rcfd->Solver_Dict->interpolation_on){

                            fprintf(f_out,"%d %d %d %d %e\n", c0, node0, c1, node1, w0);
                        }
                        else if(_C(i_layer).hit_by_other_cell[c1]){

                            fprintf(f_out,"%d %d %d %d %e\n", c0, node0, c1, node1, w0);

                            _C(i_layer).weight_after_shift[c1] += w0; // set per phase
                        }
#endif
                    }
                }

                total_number_of_stored_C2Cs = rcfd->Tracer->monitor_counter[i_phase];

                rcfd->Tracer->monitor_counter[i_phase] = 0;
            }
            else{

                fprintf(f_out,"%d \n", rcfd->Tracer->monitor_counter[i_phase]);

                total_number_of_stored_C2Cs = 0;
            }

            fclose(f_out);
        }

        return total_number_of_stored_C2Cs;

#undef LOCAL_DEBUG

#else
        return -1;
#endif
    }


/**
 * @brief node-zero writes frame_id values into file
 * @param rcfd pointer to previously created rCFD instance
 * @retval -1 Error
 * @retval 0 OK 
 */
int shift_write_frame_id(rCFD_type* rcfd)
{
#if RP_NODE

    if(myid == 0){

        int     i_id;
        
        FILE    *f_out = NULL;
        char    file_name[80];

        /* open f_out */
        {
            sprintf(file_name, "%s", rcfd->File_Dict->frame_id_filename);

            if(rcfd->Tracer->frame_counter == 0){

                f_out = fopen(file_name, "w");
            }
            else{

                f_out = fopen(file_name, "a");
            }

            if(f_out == NULL){

                report_error(rcfd, RCFD_REPORT_CATEGORY_WRITE_SHIFT, 
                             "shift_write_frame_id: f_out == NULL");
                return -1;
            }
        }

        /* write C2C shifts */
        {
            for(i_id = 0; i_id < rcfd->Rec_Dict->frame_id_vector_length; i_id++){

                fprintf(f_out,"%e ", rcfd->Rec->frame_id[i_id]);
            }
            
            fprintf(f_out,"\n");
            
            fclose(f_out);
        }
    }
    
    return 0;
#else
    return -1;
#endif
}

/*****************************************************************************/
/************************** start shift_reduce_DB  ***************************/
/*****************************************************************************/

typedef struct seq_struct
{
    short   *valid;
    
    int     *min_diff_id, *min_2nd_diff_id;
    
    double  sum_of_norm;
    
    double  *mean_norm, *sum_min_diff;

} Seq_type;

void read_norms(rCFD_type *rcfd, double** norms, int frame_count, int int_cell_count)
{
#if RP_NODE
    int  i_frame, i_cell, i_tmp;
    int  my_state = rcfd->Solver->current_state;
    char filename[80]; 
    FILE *fi = NULL;

    sprintf(filename, "%s_%d_%d", rcfd->File_Dict->Norm_filename, my_state, myid);
    fi = fopen(filename, "r");

    if(!fi){    

        report_error(rcfd, RCFD_REPORT_CATEGORY_WRITE_SHIFT, 
                     "shift_reduce_database: read_norms fi == NULL");

        return;
    }

    for(i_frame = 0; i_frame < frame_count; i_frame++){
        
        fscanf(fi, "%d\n", &i_tmp);
        
        if(i_tmp != int_cell_count){
            
            report_error(rcfd, RCFD_REPORT_CATEGORY_WRITE_SHIFT, 
                     "shift_reduce_database: read_norms: i_tmp != int_cell_count");

            return;
        }
        
        for(i_cell = 0; i_cell < int_cell_count; i_cell++){
            
            fscanf(fi, "%le\n", &norms[i_frame][i_cell]);
        }
    }

    fclose(fi);
#endif  
}

void set_single_seq_diff(rCFD_type *rcfd, Seq_type *seq, double **norms, int i_seq, int seq_count, int seq_length, int int_cell_count)
{
#if RP_NODE

    int i_seq2, i_in_seq, i_frame, i_frame2, i_cell;
    
    double dist_seq2seq;
    
    double seq_min_diff =     seq->sum_of_norm;
    double seq_min_2nd_diff = seq->sum_of_norm;

    for(i_seq2 = 0; i_seq2 < seq_count; i_seq2++){
    
        if((i_seq != i_seq2) && (seq->valid[i_seq2])){

            dist_seq2seq = 0.0;

            for(i_in_seq = 0; i_in_seq < seq_length; i_in_seq++){

                i_frame  = i_seq  * seq_length + i_in_seq;
                i_frame2 = i_seq2 * seq_length + i_in_seq;

                for(i_cell = 0; i_cell < int_cell_count; i_cell++){

                    if(fabs(norms[i_frame][i_cell]) > seq->mean_norm[i_frame]){

                        dist_seq2seq += fabs(norms[i_frame][i_cell]-norms[i_frame2][i_cell]);
                    }
                }
            }
            
            dist_seq2seq = PRF_GRSUM1(dist_seq2seq);

            if(seq_min_diff > dist_seq2seq){

                seq_min_diff =             dist_seq2seq;
                seq->min_diff_id [i_seq] = i_seq2;
                seq->sum_min_diff[i_seq] = seq_min_diff + seq_min_2nd_diff;
            }
            else if(seq_min_2nd_diff > dist_seq2seq){

                seq_min_2nd_diff            = dist_seq2seq;
                seq->min_2nd_diff_id[i_seq] = i_seq2;
                seq->sum_min_diff   [i_seq] = seq_min_diff + dist_seq2seq;
            }
        }
    }

#endif
}

void init_seq_diff(rCFD_type *rcfd, Seq_type *seq, double **norms, int seq_count, int seq_length, int int_cell_count)
{
#if RP_NODE
    int i_seq, i_in_seq, i_frame, i_cell;
    
    int global_int_cell_count;
    
    double mean_norm;

    // set seq->mean_norm[i_seq]
    {
        seq->sum_of_norm = 0.0;
        mean_norm =   0.0;

        for(i_seq = 0; i_seq < seq_count; i_seq++){

            for(i_in_seq = 0; i_in_seq < seq_length; i_in_seq++){

                i_frame  = i_seq  * seq_length + i_in_seq;

                for(i_cell = 0; i_cell < int_cell_count; i_cell++){

                    mean_norm +=        fabs(norms[i_frame][i_cell]);
                    seq->sum_of_norm += fabs(norms[i_frame][i_cell]);
                }
            }

            mean_norm =             PRF_GRSUM1(mean_norm);
            global_int_cell_count = PRF_GISUM1(int_cell_count);
            
            if((seq_length * global_int_cell_count) > 0){

                mean_norm /= (double)(seq_length * global_int_cell_count);
            }

            seq->mean_norm[i_seq] = mean_norm;
        }
        
        seq->sum_of_norm = PRF_GRSUM1(seq->sum_of_norm);

        report_debug0(rcfd, RCFD_REPORT_CATEGORY_WRITE_SHIFT, 
                      "set_seq_diff: seq->mean_norm[10] %e", 
                      seq->mean_norm[10]);
        
    }
    
    // set seq->diff
    {
        for(i_seq = 0; i_seq < seq_count; i_seq++){
            
            set_single_seq_diff(rcfd, seq, norms, i_seq, seq_count, seq_length, int_cell_count);
        }
    }
#endif
}

int get_obsolet_seq_id(rCFD_type *rcfd, Seq_type *seq, int seq_count)
{
#if RP_NODE
    int    i_seq;

    int    min_diff_id = -1;

    double min_diff = 1e10; 
    
    for(i_seq = 0; i_seq < seq_count; i_seq++){
        
        if((seq->valid[i_seq]) &&
           (seq->sum_min_diff[i_seq] <= min_diff)){

            min_diff = seq->sum_min_diff[i_seq];

            min_diff_id = i_seq;
        }
    }
    
    return min_diff_id;
#else
    return -1;
#endif
}

void adjust_seq_diff(rCFD_type *rcfd, Seq_type *seq, double **norms, int obsolet_seq_id, int seq_count, int seq_length, int int_cell_count)
{
#if RP_NODE
    int i_seq;

    int obsolet_neigh_id =     seq->min_diff_id    [obsolet_seq_id];
    int obsolet_2nd_neigh_id = seq->min_2nd_diff_id[obsolet_seq_id];

    for(i_seq = 0; i_seq < seq_count; i_seq++){

        if((seq->valid[i_seq]) &&
           ((seq->min_diff_id[i_seq] ==     obsolet_neigh_id) ||
            (seq->min_2nd_diff_id[i_seq] == obsolet_neigh_id) ||
            (seq->min_diff_id[i_seq] ==     obsolet_2nd_neigh_id) ||
            (seq->min_2nd_diff_id[i_seq] == obsolet_2nd_neigh_id))){

            set_single_seq_diff(rcfd, seq, norms, i_seq, seq_count, seq_length, int_cell_count);
        }
    }   

#endif  
}

void rewrite_c2c(rCFD_type *rcfd, Seq_type *seq, int seq_count, int seq_length)
{
#if RP_NODE

    int i_seq, i_in_seq, i_shift;

    int c0, n0, c1, n1, shift_count;

    double w0;

    char filename_in[80], filename_out[80];
    
    FILE *fi = NULL, *fo = NULL;

    shift_count = PRF_GISUM1(myid);     // just for synchronization

    report_debug0(rcfd, RCFD_REPORT_CATEGORY_WRITE_SHIFT, 
                  "shift_reduce_database: rewrite_c2c: start for all nodes");   

    sprintf(filename_in,  "./data/c2c_raw/c2c_0_0_%d", myid);   // TODO better file management
    sprintf(filename_out, "./data/c2c/c2c_0_0_%d",     myid);
    
    fi = fopen(filename_in,  "r");
    fo = fopen(filename_out, "w");
    
    if(!fi || !fo){

        report_error(rcfd, RCFD_REPORT_CATEGORY_WRITE_SHIFT, 
                     "shift_reduce_database: rewrite_c2c fi || fo == NULL");

        return;
    }
    
    for(i_seq = 0; i_seq < seq_count; i_seq++){
        
        for(i_in_seq = 0; i_in_seq < seq_length; i_in_seq++){
            
            fscanf(fi, "%d\n", &shift_count);
            
            if(seq->valid[i_seq]) fprintf(fo, "%d\n", shift_count);
            
            for(i_shift = 0; i_shift < shift_count; i_shift++){

                fscanf(fi, "%d %d %d %d %le\n", &c0, &n0, &c1, &n1, &w0);
                
                if(seq->valid[i_seq]) fprintf(fo,  "%d %d %d %d %e\n", c0, n0, c1, n1, w0);
            }
        }
    }

    fclose(fi);
    fclose(fo);

#endif
}

void rewrite_frame_id(rCFD_type *rcfd, Seq_type *seq, int seq_count, int seq_length, int frame_id_length)
{
#if RP_NODE

    if(myid == 0){
    
        int i_seq, i_in_seq, i_id;

        double tmp;

        FILE *fi = NULL, *fo = NULL;

        fi = fopen("./data/c2c_raw/c2c_id", "r");
        fo = fopen("./data/c2c/c2c_id", "w");
        
        if(!fi || !fo){

            report_error(rcfd, RCFD_REPORT_CATEGORY_WRITE_SHIFT, 
                         "shift_reduce_database: rewrite_c2c_id fi || fo == NULL");

            return;
        }           
        
        for(i_seq = 0; i_seq < seq_count; i_seq++){
            
            for(i_in_seq = 0; i_in_seq < seq_length; i_in_seq++){
                
                for(i_id = 0; i_id < frame_id_length; i_id++){
                    
                    fscanf(fi, "%le ", &tmp);
                    
                    if(seq->valid[i_seq]) fprintf(fo, "%e ", tmp);
                }
                
                fscanf(fi, "\n");

                if(seq->valid[i_seq]) fprintf(fo, "\n");
            }
        }

        fclose(fi);
        fclose(fo);
    }
#endif
}

int shift_reduce_database(rCFD_type *rcfd, int target_seq_count)
{
#if RP_NODE

    if(!rcfd) return -1;

    int frame_count = rcfd->Solver_Dict->number_of_frames;
    int seq_length =  rcfd->Rec_Dict->frame_id_sequence_length;
    int seq_count =   (int)floor((double)frame_count / (double)seq_length);

    if(seq_count <= target_seq_count)  return -1;

#if 1   // local vars

    int obsolet_seq_id, reduced_seq_count;

    Seq_type *seq = NULL;

    int base_layer =     0;
    int int_cell_count = _CELL_DICT(base_layer).number_of_int_cells;

    double  **norms = NULL; 

#endif

    // alloc, load norms and alloc, set seq
    {
        report_debug0(rcfd, RCFD_REPORT_CATEGORY_WRITE_SHIFT, 
                      "shift_reduce_database: before alloc");   

        // alloc and load norms
        {
            norms = malloc_r_2d(frame_count, int_cell_count);

            memset(norms[0], 0.0, (frame_count * int_cell_count * sizeof(double)));
            
            read_norms(rcfd, norms, frame_count, int_cell_count);
        }
        
        // alloc and set seq
        {
            seq = malloc(sizeof(Seq_type));
            
            seq->valid =           malloc(seq_count * sizeof(short));
            seq->min_diff_id =     malloc(seq_count * sizeof(int));
            seq->min_2nd_diff_id = malloc(seq_count * sizeof(int));
            seq->mean_norm =       malloc(seq_count * sizeof(double));
            seq->sum_min_diff =    malloc(seq_count * sizeof(double));
            
            memset(seq->valid,            1, seq_count * sizeof(short));
            memset(seq->min_diff_id,     -1, seq_count * sizeof(int));
            memset(seq->min_2nd_diff_id, -1, seq_count * sizeof(int));
            memset(seq->mean_norm,      0.0, seq_count * sizeof(double));
            memset(seq->sum_min_diff,   0.0, seq_count * sizeof(double));
            
            init_seq_diff(rcfd, seq, norms, seq_count, seq_length, int_cell_count);
        }

        report_debug0(rcfd, RCFD_REPORT_CATEGORY_WRITE_SHIFT, 
                      "shift_reduce_database: after alloc");
    }
    
    // look for obsolet sequences until target_seq_count is reached
    {
        reduced_seq_count = seq_count;
        
        while(reduced_seq_count > target_seq_count){
        
            obsolet_seq_id = get_obsolet_seq_id(rcfd, seq, seq_count);
            
            report_debug0(rcfd, RCFD_REPORT_CATEGORY_WRITE_SHIFT, 
                          "shift_reduce_database: obsolet_seq_id %d", 
                          obsolet_seq_id);
                          
            seq->valid[obsolet_seq_id] = 0;
            
            adjust_seq_diff(rcfd, seq, norms, obsolet_seq_id, seq_count, seq_length, int_cell_count);
            
            reduced_seq_count--;
        }

        report_debug0(rcfd, RCFD_REPORT_CATEGORY_WRITE_SHIFT, 
                      "shift_reduce_database: after adjust_seq_diff");  
        
    }
    
    // rewrite c2c and frame_id files
    {
        // rename c2c to c2c_raw and create new c2c
        // TODO create c2c_red375 and sync nodes
        // create_new_c2c_folder(rcfd)
        if(myid == 0){
        
            int status = 0;
            
            // TODO check if ./data/c2c_raw already exists
            status |= rename("./data/c2c", "./data/c2c_raw");
            status |= mkdir ("./data/c2c", 0777);
            
            report_debug0(rcfd, RCFD_REPORT_CATEGORY_WRITE_SHIFT, 
                          "shift_reduce_database: c2c, c2c_raw: status %d", status);    
        }   
        
        rewrite_c2c(rcfd, seq, seq_count, seq_length);
        
        if(rcfd->Rec_Dict->format == RCFD_REC_FRAME_ID){
            
            int frame_id_length = rcfd->Rec_Dict->frame_id_vector_length;
            
            rewrite_frame_id(rcfd, seq, seq_count, seq_length, frame_id_length);
        }
        
        report_debug0(rcfd, RCFD_REPORT_CATEGORY_WRITE_SHIFT, 
                      "shift_reduce_database: after rewrite");          
    }
    
    // free local vars
    {
        free_r_2d(norms);
        
        norms = NULL;
        
        if(seq != NULL){

            free(seq->valid);

            free(seq->min_diff_id);
            free(seq->min_2nd_diff_id);
            
            free(seq->mean_norm);
            free(seq->sum_min_diff);
            
            free(seq);
            
            seq = NULL;
        }
        
        report_debug0(rcfd, RCFD_REPORT_CATEGORY_WRITE_SHIFT, 
                      "shift_reduce_database: after free");
    }
    
    return reduced_seq_count;
#else
    return -1;
#endif
}

#endif // RCFD_SHIFT