ivp_multidimensional_interp.cxx

hl2 source code. Do not use it illegal.

#endif
	success = IVP_OK;
    } else {
	//more than 1 input vector needed for the interpolation
	if ((nr_of_vectors_involved < 2) || (nr_of_vectors_involved > nr_occupied)) {
	    //in the previous interpolation run one vector was sufficient (but not now)
	    nr_of_vectors_involved = 2;
	}
	
	while (nr_of_vectors_involved <= nr_occupied) {
	    //initialize the scratchboard
	    for (int i=0; i<(nr_of_vectors_involved-1); i++) {
		scratchboard_A[i]->set(previous_inputs[i+1]);
		scratchboard_A[i]->subtract(previous_inputs[0]);
	    }
	    scratchboard_A[nr_of_vectors_involved-1]->set(new_input_difference); //now the matrix [A|b] is ready for LINFIT

	    //run the linfit algorithm
	    success = linfit(nr_of_elements_input, nr_of_vectors_involved-1, scratchboard_A, interpolation_weights, &res);

#ifdef WITH_DEBUG_OUTPUT
	    printf("nr_of_vectors_involved = %d, nr_occupied = %d, res = %1.4e\n", nr_of_vectors_involved, nr_occupied, res);
#endif
	    
	    if (success) {

		interpolation_weigh_first = 1.0f;
		
		if (res > max_res) {
		    success = IVP_FAULT;
#ifdef WITH_DEBUG_OUTPUT
		    printf("res is over limit!\n");
		    nr_res_over_limit[nr_of_vectors_involved-1]++;
#endif
		} else {
		
		    //ensure that the interpolation_weights are in the interval [-0.2f; 1.2f]
		    for (int ix=0; ix<(nr_of_vectors_involved-1); ix++) {
			if (IVP_Inline_Math::fabsd(interpolation_weights[ix] - 0.5f) > 0.7f) {
			    success = IVP_FAULT;
#ifdef WITH_DEBUG_OUTPUT
			    printf("interpolation_weigh is over limit!\n");
			    nr_int_weight_over_limit[nr_of_vectors_involved-1]++;
#endif
			    break;
			}
			interpolation_weigh_first -= interpolation_weights[ix];
		    }

		    if ((IVP_Inline_Math::fabsd(interpolation_weigh_first - 0.5f) > 0.7f) && success) {
			success = IVP_FAULT;
#ifdef WITH_DEBUG_OUTPUT
			printf("interpolation_weigh_first is over limit!\n");
			nr_int_weight_over_limit[nr_of_vectors_involved-1]++;
#endif
		    }
		}

#ifdef WITH_DEBUG_OUTPUT
		printf("int_w[0] = %e ", interpolation_weigh_first);
		for (int g=0; g<nr_of_vectors_involved-1; g++) {
		    printf("int_w[%d] = %e ", g+1, interpolation_weights[g]);
		}
		printf("\n");
#endif
		
		if (success) {
#ifdef DEBUG
		  nr_of_success[nr_of_vectors_involved-1]++;
#endif
		  break;
		}
	    } else {
		//linfit could not solve the interpolation
#ifdef DEBUG
		nr_of_linfit_failure[nr_of_vectors_involved-1]++;
#endif
	    }
	    nr_of_vectors_involved++;
	}
    }
    
    	
    if (success) {  //new_input was successfully interpolated, the weights of the old input vectors involved can be found in interpolation_weights
	
	//clear the space for the new solution
	for (int j=0; j<nr_of_elements_solution; j++) {
	    previous_solutions[scratch_area_index]->set(j, 0.0f);
	}

	
	//calc weight statistics for every vector involved and
	//calc the output as result of the linear combination of the solution vectors
	for (int i=1; i<nr_of_vectors_involved; i++) {
	    previous_inputs[i]->weight_statistic = previous_inputs[i]->weight_statistic * influence_of_old_weight + interpolation_weights[i-1];

	    for (int k=0; k<nr_of_elements_solution; k++) {
		previous_solutions[scratch_area_index]->set(k, previous_solutions[scratch_area_index]->element[k] + interpolation_weights[i-1] * previous_solutions[i]->element[k]);
	    }

	}

	//calc weigh statistics for the first vector and add first vector's share to the output
	previous_inputs[0]->weight_statistic = (previous_inputs[0]->weight_statistic * influence_of_old_weight) + interpolation_weigh_first;
	for (int k=0; k<nr_of_elements_solution; k++) {
	    previous_solutions[scratch_area_index]->set(k, previous_solutions[scratch_area_index]->element[k] + interpolation_weigh_first * previous_solutions[0]->element[k]);
	}
	
	output->set(previous_solutions[scratch_area_index]);

#ifdef WITH_DEBUG_OUTPUT
	{
	    IVP_MI_Vector *test_vector;
	    IVP_MI_VECTOR_ALLOCA(test_vector, nr_of_elements_input);
	    for (int i=1; i<nr_of_vectors_involved; i++) {
		for (int k=0; k<nr_of_elements_input; k++) {
		    test_vector->set(k, test_vector->element[k] + interpolation_weights[i-1] * previous_inputs[i]->element[k]);
		}
	    }

	    for (int k=0; k<nr_of_elements_input; k++) {
		test_vector->set(k, test_vector->element[k] + interpolation_weigh_first * previous_inputs[0]->element[k]);
	    }

	    printf("interpolated_output_vector:\n");
	    output->print();
	}
#endif
	//sort the vectors (i.e. index[]) in declining order of their weight statistics
	sort_vectors(nr_occupied);
    }
	
#ifdef WITH_DEBUG_OUTPUT
    for (int y=0; y<nr_occupied; y++) {
	printf("timestamp = %f  ", previous_inputs[y]->time_stamp);
	previous_inputs[y]->print();
    }
#endif
    
    return(success);
}
#endif


void IVP_Multidimensional_Interpolator::add_new_input_solution_combination_conventional(const IVP_MI_Vector *new_input, const IVP_MI_Vector *new_solution) {

    //move all vectors one position further, the new vectors will be copied to the first position
    IVP_MI_Vector *tmp_vec_input = previous_inputs[nr_of_vectors-1];
    IVP_MI_Vector *tmp_vec_solution = previous_solutions[nr_of_vectors-1];
    for (int i=nr_of_vectors-2; i>=0; i--) {
	previous_inputs[i+1] = previous_inputs[i];
	previous_solutions[i+1] = previous_solutions[i];
    }
    previous_inputs[0] = tmp_vec_input;
    previous_solutions[0] = tmp_vec_solution;
    
    //copy new vectors into first position
    previous_inputs[0]->set(new_input);
    previous_solutions[0]->set(new_solution);

    //debugging!!!!!!!!!!!!!!!!!!!!!!!!!
    previous_inputs[0]->set_time_stamp(new_input->time_stamp);
    previous_solutions[0]->set_time_stamp(new_solution->time_stamp);
    

    //give the new input vector a weigh that ensures it will be positioned very high after a sort
    previous_inputs[0]->weight_statistic = 1.0f / (1.0f - influence_of_old_weight);

    if (nr_occupied < nr_of_vectors)
	nr_occupied++;
}


void IVP_Multidimensional_Interpolator::add_new_input_solution_combination_stochastic(const IVP_MI_Vector *new_input, const IVP_MI_Vector *new_solution) {

    //define the location of the new pair of values by stochastic means
    int solution_slot = (counter_vector_replacement-- * 101) % nr_of_vectors;
    
    previous_inputs[solution_slot]->set(new_input);
    previous_solutions[solution_slot]->set(new_solution);

    //debugging!!!!!!!!!!
    previous_inputs[solution_slot]->set_time_stamp(new_input->time_stamp);
    previous_solutions[solution_slot]->set_time_stamp(new_solution->time_stamp);

    //give the new input vector a weigh that ensures it will be positioned very high after a sort
    previous_inputs[solution_slot]->weight_statistic = 1.0f / (1.0f - influence_of_old_weight);

    sort_vectors(nr_occupied);
    
    //update the counter
    if (counter_vector_replacement < 0)
	counter_vector_replacement = initial_value_vector_replacement;
}



/************************
 * Constructor
 ************************/
IVP_Multidimensional_Interpolator::IVP_Multidimensional_Interpolator(int nr_of_vectors_, int nr_of_elements_input_, int nr_of_elements_solution_) {

    int i;
    nr_of_vectors = nr_of_vectors_;
    nr_of_elements_input = nr_of_elements_input_;
    nr_of_elements_solution = nr_of_elements_solution_;
    //max_tries_nr_of_vectors_involved = max_tries_nr_of_vectors_involved_;

    //initialize counter_tries_nr_of_vectors_involved
    counter_tries_nr_of_vectors_involved = 0;

    //initialize counter_vector_replacement
    initial_value_vector_replacement = 50;
    counter_vector_replacement = initial_value_vector_replacement;

    //previous_inputs = new IVP_MI_Vector[nr_of_vectors](nr_of_elements_input);  //contains vectors of input data from previous runs of a solver
    previous_inputs = (IVP_MI_Vector **)p_malloc( sizeof(void*) * (nr_of_vectors) );
    for (i=0; i<nr_of_vectors; i++) {
	previous_inputs[i] = IVP_MI_Vector::malloc_mi_vector(nr_of_elements_input);
    }
    
    //IMPORTANT: The vectors have to be linearly unrelated
    previous_solutions = (IVP_MI_Vector **)p_malloc( sizeof(void*) * (nr_of_vectors+1) );
    for (i=0; i<=nr_of_vectors; i++) {
	previous_solutions[i] = IVP_MI_Vector::malloc_mi_vector(nr_of_elements_solution);
    }
    
    influence_of_old_weight = 0.8f;
	
    nr_occupied = 0;
    nr_of_vectors_involved = 2;

    scratch_area_index = nr_of_vectors;
    MI_eps = 1.0e-14f;
    //max_res = 1.0e-2f;

    //begin debug variables!!!!!!!!!!!!
#ifdef DEBUG
    nr_one_vector_sufficient = 0;
    nr_res_over_limit = new int[nr_of_vectors];
    nr_int_weight_over_limit = new int[nr_of_vectors];
    nr_of_linfit_failure = new int[nr_of_vectors];
    nr_of_success = new int[nr_of_vectors];
    for (i=0; i<nr_of_vectors; i++) {
	nr_res_over_limit[i] = 0;
	nr_int_weight_over_limit[i] = 0;
	nr_of_linfit_failure[i] = 0;
	nr_of_success[i] = 0;
    }
#endif
    //end debug variables!!!!!!!!!!!!!
}