#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <string.h>
//#include <conio.h>
#include <time.h>
#include <math.h>

#define NAZEV "C:\\Users\\Omajo\\Desktop\\VUT\\SP_code\\2D.dat"
#define NAZEVVYSTUPU "C:\\Users\\Omajo\\Desktop\\VUT\\SP_code\\CosxplusxVystup1.txt"
#define NAZEVVYSTUPU2 "C:\\Users\\Omajo\\Desktop\\VUT\\SP_code\\CosxplusxVystup2.txt"
#define NAZEVVYSTUPU3 "C:\\Users\\Omajo\\Desktop\\VUT\\SP_code\\CosxplusxVystup3.txt"
#define NAZEVVYSTUPU4 "C:\\Users\\Omajo\\Desktop\\VUT\\SP_code\\CosxplusxVystup4.txt"
#define NAZEVVYSTUPU5 "C:\\Users\\Omajo\\Desktop\\VUT\\SP_code\\CosxplusxVystup5.txt"
#define NAZEVVYSTUPU6 "C:\\Users\\Omajo\\Desktop\\VUT\\SP_code\\CosxplusxVystup6.txt"
#define iterace 1000
#define dtheta 1
#define dph 1
#define npoints 1000
#define step 0.5
#define fi 1.0/1.618033988749894848

int compare_ints(const void* a, const void* b)
{
	if (*(double*)a > *(double*)b) return 1;
	else if (*(double*)a < *(double*)b) return -1;
	else return 0;
}
int nacti (double **pole1, double **pole2, double **pole3,int *velikost)
{
	int i,size;
	double x;
	
	FILE *data;
	
	if ((data = fopen(NAZEV, "r")) == NULL) {
        fprintf(stderr, "Chyba otevreni souboru %s pro cteni\n", NAZEV);
        return 1;
    }
	
	i=0;
	while(fscanf(data, "%lf", &x) != EOF){
	i++;	
	}
	
	size=i/3;
	*velikost=i/3;
	//polea = (double *) malloc(sizeof(double) * (size_t)size);
	
	if ((*pole1 = (double *) malloc(sizeof(double) * (size_t)size)) == NULL){
        printf("Nedostatek pameti!\n");
        return 1;
    }
	if ((*pole2 = (double *) malloc(sizeof(double) * (size_t)size)) == NULL){
        printf("Nedostatek pameti!\n");
        return 1;
    }
	if ((*pole3 = (double *) malloc(sizeof(double) * (size_t)size)) == NULL){
        printf("Nedostatek pameti!\n");
        return 1;
    }

	rewind (data);
	
	for (i = 0; i <= size; i++) {
        fscanf(data, "%lf %lf %lf", *pole1 + i, *pole2 + i, *pole3 + i);
    }
	
	//printf("%lf %lf %lf \n", **(pole1), **(pole2), **(pole3));
	//printf("%lf %lf %lf \n", **(pole1+size-1), **(pole2+size-1), **(pole3+size-1));
	fclose(data);
	return 0;
}
double interpoluj(double x1, double x2, double *polex1, double *polex2, double *poley, int size)
{
	int i,nalezeno;
	nalezeno=0;
	double A[3];
	double B[3];
	double C[3];
	double D[3];
	double X[3];
	
	for(i = 0; i < 3; i++)	//inicializace bodu
	{
		A[i]=0;
		B[i]=0;
		C[i]=0;
		D[i]=0;
		X[i]=0;
	}
	
	X[0]=x1;
	X[1]=x2;
	
	if(x1<x2){
		for(i = 0; i <= size; i++){
			if(X[0] == *(polex1+i) && X[1] == *(polex2+i) ){		//Prvni test zda se hodnota v pozadovanem bode nenachazi v souboru hodnot. 
					X[2]=*(poley+i);								//Pokud ano pouze tuto hodnotu vratim.
					nalezeno=1;												
					break;
			}
		}
		if(nalezeno==0){
			for(i = 0; i <= size; i++){								//Pokud se bod v souboru nenachazi je treba nalezt nejblizsi body pravidelne
				if(*(polex1+i) <= X[0] && *(polex1+i) > A[0]){		//prostorove mrizky, aby bylo mozne interpolovat.
						A[0] = *(polex1+i);
				}
				if(*(polex2+i) <= X[1] && *(polex2+i) > A[1]){
						A[1] = *(polex2+i);
				}
			}
			D[0]=A[0];
			D[1]=*(polex2+size-1);
			B[1]=A[1];
			B[0]=*(polex1+size-1);
			
			for(i = 0; i <= size; i++){								
				if(*(polex1+i) >= X[0] && *(polex1+i) < B[0]){		
						B[0] = *(polex1+i);
				}
			}	
			
			for(i = 0; i <= size; i++){								
				if(*(polex2+i) >= X[1] && *(polex2+i) < D[1]){		
						D[1] = *(polex2+i);
				}
			}
			C[0]=B[0];
			C[1]=D[1];
		}
		
		for(i = 0; i <= size; i++){									//prirazeni funkcni hodnoty bodum
			if(A[0] == *(polex1+i) && A[1] == *(polex2+i) ){ 
					A[2]=*(poley+i);								
			}
			if(B[0] == *(polex1+i) && B[1] == *(polex2+i) ){ 
					B[2]=*(poley+i);								
			}
			if(C[0] == *(polex1+i) && C[1] == *(polex2+i) ){ 
					C[2]=*(poley+i);								
			}
			if(D[0] == *(polex1+i) && D[1] == *(polex2+i) ){ 
					D[2]=*(poley+i);								
			}
		}													
	}
	else if(x1>x2){
		for(i = 0; i <= size; i++){
			if(X[1] == *(polex1+i) && X[0] == *(polex2+i) ){		//Prvni test zda se hodnota v pozadovanem bode nenachazi v souboru hodnot. 
					X[2]=*(poley+i);								//Pokud ano pouze tuto hodnotu vratim.
					nalezeno=1;												
					break;
			}
		}
		if(nalezeno==0){
			for(i = 0; i <= size; i++){								//Pokud se bod v souboru nenachazi je treba nalezt nejblizsi body pravidelne
				if(*(polex1+i) <= X[1] && *(polex1+i) > A[0]){		//prostorove mrizky, aby bylo mozne interpolovat.
						A[0] = *(polex1+i);
				}
				if(*(polex2+i) <= X[0] && *(polex2+i) > A[1]){
						A[1] = *(polex2+i);
				}
			}
			D[0]=A[0];
			D[1]=*(polex2+size-1);
			B[1]=A[1];
			B[0]=*(polex1+size-1);
			
			for(i = 0; i <= size; i++){								
				if(*(polex1+i) >= X[1] && *(polex1+i) < B[0]){		
						B[0] = *(polex1+i);
				}
			}	
			
			for(i = 0; i <= size; i++){								
				if(*(polex2+i) >= X[0] && *(polex2+i) < D[1]){		
						D[1] = *(polex2+i);
				}
			}
			C[0]=B[0];
			C[1]=D[1];
		}	
		for(i = 0; i <= size; i++){									//prirazeni funkcni hodnoty bodum
			if(A[0] == *(polex1+i) && A[1] == *(polex2+i) ){ 
					A[2]=*(poley+i);								
			}
			if(B[0] == *(polex1+i) && B[1] == *(polex2+i) ){ 
					B[2]=*(poley+i);								
			}
			if(C[0] == *(polex1+i) && C[1] == *(polex2+i) ){ 
					C[2]=*(poley+i);								
			}
			if(D[0] == *(polex1+i) && D[1] == *(polex2+i) ){ 
					D[2]=*(poley+i);								
			}
		}	
	}
	else{
		for(i = 0; i <= size; i++){
			if(X[0] - *(polex1+i)  <= 0.005 && X[1] - *(polex2+i)  <= 0.005 ){		//Pokud x1=x2, tak hodnotu naleznu s presnosti na 0.005 primo ze zadanych hodnot.
					X[2]=*(poley+i);												
					nalezeno=1;												
					break;
			}
		}		
	}
	
	if(nalezeno==0){	//Pokud nebyla nalezena pozadovana funkcni hodnota primo z hodnot, tak interpoluji.
	//printf("Interpoluji\n");
					if(x1<x2){
						X[2]=(((A[2])/((B[0]-A[0])*(D[1]-A[1])))*(B[0]-X[0])*(D[1]-X[1])) + 
						(((B[2])/((B[0]-A[0])*(D[1]-A[1])))*(X[0]-A[0])*(D[1]-X[1]))+
						(((D[2])/((B[0]-A[0])*(D[1]-A[1])))*(B[0]-X[0])*(X[1]-A[1])) + 
						(((C[2])/((B[0]-A[0])*(D[1]-A[1])))*(X[0]-A[0])*(X[1]-A[1]));
					}
					else if(x1>x2){
						X[2]=(((A[2])/((B[0]-A[0])*(D[1]-A[1])))*(B[0]-X[1])*(D[1]-X[0])) + 
						(((B[2])/((B[0]-A[0])*(D[1]-A[1])))*(X[1]-A[0])*(D[1]-X[0]))+
						(((D[2])/((B[0]-A[0])*(D[1]-A[1])))*(B[0]-X[1])*(X[0]-A[1])) + 
						(((C[2])/((B[0]-A[0])*(D[1]-A[1])))*(X[1]-A[0])*(X[0]-A[1]));
					}
	}

	
	/*printf("Bod X[%lf %lf %lf]\n", X[0], X[1], X[2]);
	printf("Bod A[%lf %lf %lf]\n", A[0], A[1], A[2]);
	printf("Bod B[%lf %lf %lf]\n", B[0], B[1], B[2]);
	printf("Bod C[%lf %lf %lf]\n", C[0], C[1], C[2]);
	printf("Bod D[%lf %lf %lf]\n", D[0], D[1], D[2]);
	*/
	return X[2];
}
double nahodnehledani(double *polex1, double *polex2,double *poley,double *bodmin, int size)
{
	double min,x,y,z;
	int i;
	min=*(polex1+size-1);
	
	srand((unsigned int) time(NULL));  
	
	for (i = 0; i < iterace; i++){ 
		
		x = *polex1+((double)rand() / (double)RAND_MAX)*(*(polex1+size-1)-*polex1);
		y = *polex2+((double)rand() / (double)RAND_MAX)*(*(polex2+size-1)-*polex2);
	
		z=interpoluj(x,y,polex1 ,polex2 ,poley ,size);
		
        if (z<min)
		{
		min=z;
		}
	}
	
	*bodmin = x;
	*(bodmin+1) = y;
	return min;
}
double evolucnistrategie(double *polex1, double *polex2,double *poley, double *bodmin, int size)
{
	double min,x,x1,y,y1,z,smodch,krok,uspech,pravidlo,c1,c2;
	smodch = (*(polex1+size-1) - *polex1) / 3;
	int i;
	min=*(polex1+size-1);
	krok = 0;
	uspech = 0;
	pravidlo=0.2;
	c1=0.82;
	c2=1.22;
	
	srand((unsigned int) time(NULL)); 
	
	x = *polex1+((double)rand() / (double)RAND_MAX)*(*(polex1+size-1)-*polex1);
	y = *polex2+((double)rand() / (double)RAND_MAX)*(*(polex2+size-1)-*polex2);
	
	for (i = 0; i < iterace; i++){ 
		
		while(1){
		x1 = x + smodch * (-1 + ((double)rand() / (double)RAND_MAX) * 2);
		y1 = y + smodch * (-1 + ((double)rand() / (double)RAND_MAX) * 2);
		
		if(x1>*polex1 && x1<*(polex1+size-1) && y1>*polex2 && y1<*(polex1+size-1)){
			break;
		}
		}
		z=interpoluj(x1,y1,polex1 ,polex2 ,poley ,size);
		krok++;
		
        if (z<min){
			uspech++;
			min=z;
			x=x1;
			y=y1;
		}
		
		if((uspech/krok)<pravidlo){
			smodch=smodch*c1;
		}
		else if((uspech/krok)>pravidlo){
			smodch=smodch*c2;
		}
	}
	
	*bodmin = x;
	*(bodmin+1) = y;	
	return min;
}
double vzdalenost(double *boda, double *bodb, double theta, double ph, int size)
{
	double mojevzdalenost;
	int i;
	
	mojevzdalenost = 0;
		
	for (i = 0; i < size; i++){ 
		mojevzdalenost = mojevzdalenost + (theta*(pow(fabs(*(boda+i)-*(bodb+i)),ph)));
	}

	return mojevzdalenost;
}
int maticeR(double *polebodu, double matrixR[100][100], int pocetbodu, int rad)
{
	int i,j,k;
	//double a,b,c,d;
	double boda[rad];
	double bodb[rad];
	
	for (i = 0; i < pocetbodu; i++) {				//Vypocet matice R
		for (j = 0; j < pocetbodu; j++){
			
			if(i==j){
				matrixR[i][j] = 1;
			}
			else if(i<j){
				for (k = 0; k < rad; k++){
					boda[k]=*(polebodu+(i*(rad+1))+k);
					bodb[k]=*(polebodu+(j*(rad+1))+k);
					/*if(k==0){
						a=*(polebodu+(i*rad)+k);
						c=*(polebodu+(j*rad)+k);
					}
					if(k==1){
						b=*(polebodu+(i*rad)+k);
						d=*(polebodu+(j*rad)+k);
					}*/
				}
				matrixR[i][j] = exp(-(vzdalenost(&boda[0], &bodb[0], dtheta, dph, rad)));
				matrixR[j][i] = exp(-(vzdalenost(&boda[0], &bodb[0], dtheta, dph, rad)));
			}
		}
    }
	
	return 0;
}
double LUdeterminant(double A[100][100], int rad)
{
	double L[100][100], U[100][100], sum;
    int i = 0, j = 0, k = 0, n = 0;

	n = rad;
	sum = 1;
    
    for (i = 0; i < n; i++)
    {
        for (j = 0; j < n; j++)
        {
            if (j < i)
                L[j][i] = 0.0;
            else
            {
                L[j][i] = A[j][i];
                for (k = 0; k < i; k++)
                {
                    L[j][i] = L[j][i] - L[j][k] * U[k][i];
                }
            }
        }
        for (j = 0; j < n; j++)
        {
            if (j < i)
                U[i][j] = 0.0;
            else if (j == i)
                U[i][j] = 1.0;
            else
            {
                U[i][j] = A[i][j] / L[i][i];
                for (k = 0; k < i; k++)
                {
                    U[i][j] = U[i][j] - ((L[i][k] * U[k][j]) / L[i][i]);
                }
            }
        }
    }
								// Determinant je soucin determinantu matic L a U, coz jsou matice v trojuhelnikovem tvaru.
    							// Jejich deterinant je soucin prvku na diagonale.
        for(j=0; j<n; j++)
		{
				sum = sum * L[j][j]*U[j][j];
		}

	
	if(sum != sum){
		return 0;
	}
	else{
		return sum;
	}	
}
double determinant(double a[100][100], int k)
{
 double s = 1, det = 0, b[100][100];
 int i, j, m, n, c;
 if (k == 1) {
	return (a[0][0]);
 } 
 else {
	det = 0;
	for (c = 0; c < k; c++) {
		m = 0;
		n = 0;
		for (i = 0; i < k; i++) {
			for (j = 0; j < k; j++) {
				b[i][j] = 0;
				if (i != 0 && j != c) {
					b[m][n] = a[i][j];
					if (n < (k - 2))
					n++;
				else {
					n = 0;
					m++;
				}
				}
			}
		}

		det = det + s * (a[0][c] * determinant(b, k - 1));
		s = -1 * s;
	}
 }
 return (det);
}
double Determinant(double **a,int n)
{
    int i,j,j1,j2 ;                    // general loop and matrix subscripts
    double det = 0 ;                   // init determinant
    double **m = NULL ;                // pointer to pointers to implement 2d
                                       // square array

    if (n < 1)    {   }                // error condition, should never get here

    else if (n == 1) {                 // should not get here
        det = a[0][0] ;
        }

    else if (n == 2)  {                // basic 2X2 sub-matrix determinate
                                       // definition. When n==2, this ends the
        det = a[0][0] * a[1][1] - a[1][0] * a[0][1] ;// the recursion series
        }


                                       // recursion continues, solve next sub-matrix
    else {                             // solve the next minor by building a
                                       // sub matrix
        det = 0 ;                      // initialize determinant of sub-matrix

                                           // for each column in sub-matrix
        for (j1 = 0 ; j1 < n ; j1++) {
                                           // get space for the pointer list
            m = (double **) malloc((n-1)* sizeof(double *)) ;

            for (i = 0 ; i < n-1 ; i++)
                m[i] = (double *) malloc((n-1)* sizeof(double)) ;
                       //     i[0][1][2][3]  first malloc
                       //  m -> +  +  +  +   space for 4 pointers
                       //       |  |  |  |          j  second malloc
                       //       |  |  |  +-> _ _ _ [0] pointers to
                       //       |  |  +----> _ _ _ [1] and memory for
                       //       |  +-------> _ a _ [2] 4 doubles
                       //       +----------> _ _ _ [3]
                       //
                       //                   a[1][2]
                      // build sub-matrix with minor elements excluded
            for (i = 1 ; i < n ; i++) {
                j2 = 0 ;               // start at first sum-matrix column position
                                       // loop to copy source matrix less one column
                for (j = 0 ; j < n ; j++) {
                    if (j == j1) continue ; // don't copy the minor column element

                    m[i-1][j2] = a[i][j] ;  // copy source element into new sub-matrix
                                            // i-1 because new sub-matrix is one row
                                            // (and column) smaller with excluded minors
                    j2++ ;                  // move to next sub-matrix column position
                    }
                }

            det += pow(-1.0,1.0 + j1 + 1.0) * a[0][j1] * Determinant(m,n-1) ;
                                            // sum x raised to y power
                                            // recursively get determinant of next
                                            // sub-matrix which is now one
                                            // row & column smaller

            for (i = 0 ; i < n-1 ; i++) free(m[i]) ;// free the storage allocated to
                                            // to this minor's set of pointers
            free(m) ;                       // free the storage for the original
                                            // pointer to pointer
        }
    }
    return(det) ;
}
int trans(double num[100][100], double fac[100][100], double invMatrix[100][100], int r) 
{
 int i, j;
 double b[100][100], inv[100][100], d;
 
 for (i = 0; i < r; i++) {
	for (j = 0; j < r; j++) {
   b[i][j] = fac[j][i];
  }
 }
 
 d = LUdeterminant(num, r);
 //printf("ftrans determinant = %20.15lf\n", d);
 inv[i][j] = 0.0;
 for (i = 0; i < r; i++) {
  for (j = 0; j < r; j++) {
   inv[i][j] = b[i][j] / d;
  }
 }
 
 //printf("\nTHE INVERSE OF THE MATRIX:\n");
 for (i = 0; i < r; i++) {
  for (j = 0; j < r; j++) {
   //printf("\t%20.15lf", inv[i][j]);
   invMatrix[i][j] = inv[i][j];
  }
  //printf("\n");
 }
 return 0;
}
int cofactors(double num[100][100], double invMatrix[100][100], int f) 
{
	double b[100][100], fac[100][100];
	int p, q, m, n, i, j;
	
	for (q = 0; q < f; q++) {
		for (p = 0; p < f; p++) {
			m = 0;
			n = 0;
			for (i = 0; i < f; i++) {
				for (j = 0; j < f; j++) {
					b[i][j] = 0.0;
					if (i != q && j != p) {
						b[m][n] = num[i][j];
						if (n < (f - 2))
							n++;
						else {
						n = 0;
						m++;
						}
					}
				}
			}
			fac[q][p] = pow(-1, q + p) * LUdeterminant(b, f - 1);
			//printf("%d determinant = %20.15lf\n", f-1, determinant(b, f - 1));
			//printf("%d LUdeterminant = %20.15lf\n", f-1, LUdeterminant(b, f - 1));
			}
	}	
	trans(num,fac,invMatrix,f);
	return 0;
}
void multiply (double matice1[100][100], double matice2[100][100], double maticevys[100][100], int pocr1, int pocs1, int pocr2, int pocs2)
{
	int i, j, k;
	double sum=0;
	if (pocs1 != pocr2){
    printf("Matice nelze nasobit\n");
	}
	else{
		for (i = 0; i < pocr1; i++) {
			for (j = 0; j < pocs2; j++) {
				for (k = 0; k < pocr2; k++) {
					sum = sum + matice1[i][k]*matice2[k][j];
				}
 
				maticevys[i][j] = sum;
				sum = 0;
			}
		}
	}
}
void fjednotkovaMatice(double matrixR[100][100], int rad)
{
	int i, j;
	
	for (i = 0; i < rad; i++) {
		for (j = 0; j < rad; j++) {
			if(i==j){
				matrixR[i][j] = 1;
			}
			else{
				matrixR[i][j] = 0;
				
			}
		}
	}
}
double mi(double *poleznamychbodu, double invMatrixR[100][100], int size, int rad)
{
	int i;
	double vectorof1[100][100], vectorofy[100][100], resultmatrix[100][100], resultmatrix2[100][100], citatel, jmenovatel;
	for (i = 0; i < size; i++){
		vectorof1[0][i] = 1;
	}
	multiply(vectorof1,invMatrixR,resultmatrix,1,size,size,size);		//= 1´R
	
	for (i = 0; i < size; i++){
		vectorofy[i][0] = *(poleznamychbodu+((i*(rad+1))+rad));
	}
	multiply(resultmatrix,vectorofy,resultmatrix2,1,size,size,1);		//= 1'Ry
	citatel = resultmatrix2[0][0];
	multiply(vectorof1,invMatrixR,resultmatrix,1,size,size,size);		//= 1´R
	for (i = 0; i < size; i++){
		vectorof1[i][0] = 1;
	}
	multiply(resultmatrix,vectorof1,resultmatrix2,1,size,size,1);		//= 1´R1
	jmenovatel = resultmatrix2[0][0];
	return citatel/jmenovatel;										//u = (1'Ry)/(1'R1)
}
double predictor(double *odhadovanybod, double *poleznamychbodu, double invMatrixR[100][100], int size, int rad)
{
	int i, j;
	double u, vektorr[100][100],vectorofy[100][100], resultmatrix[100][100], resultmatrix2[100][100], boda[rad], bodb[rad];
	u = mi(poleznamychbodu,invMatrixR,size,rad);
	for (j = 0; j < rad; j++){
		boda[j] = *(odhadovanybod+j);
	}
	for (i = 0; i < size; i++){
		for (j = 0; j < rad; j++){
			bodb[j] = *(poleznamychbodu+(i*(rad+1))+j);
		}
		vektorr[0][i] = exp(-vzdalenost(&boda[0],&bodb[0],dtheta,dph,rad));		//vypocet korelacniho vektoru r'
	}
	for (j = 0; j < size; j++){
		vectorofy[j][0] = *(poleznamychbodu+(j*(rad+1))+rad) - u;	// = (y-1u)
	}
	multiply(vektorr,invMatrixR,resultmatrix,1,size,size,size);		// = r'R^-1
	multiply(resultmatrix,vectorofy,resultmatrix2,1,size,size,1);	// = r'R^-1(y-1u)
	return resultmatrix2[0][0]+u; 									// y^(x*) = u + r'R^-1(y-1u)
}
double sigma(double *poleznamychbodu, double invMatrixR[100][100], int size, int rad)
{
	int i;
	double vectorofy[100][100], resultmatrix[100][100], resultmatrix2[100][100], u;
	u = mi(poleznamychbodu,invMatrixR,size,rad);
	for (i = 0; i < size; i++){
		vectorofy[0][i] = (*(poleznamychbodu+((i*(rad+1))+rad))) - u;		// = (y - 1u)'
	}
	multiply(vectorofy,invMatrixR,resultmatrix,1,size,size,size);			// = (y - 1u)'R^-1
	for (i = 0; i < size; i++){
		vectorofy[i][0] = (*(poleznamychbodu+((i*(rad+1))+rad))) - u;		// = (y - 1u)
	}
	multiply(resultmatrix,vectorofy,resultmatrix2,1,size,size,1);			// = (y - 1u)'R^-1 (y - 1u)
	
	return resultmatrix2[0][0]/size;										// = (y - 1u)'R^-1 (y - 1u) / n
}
double meansquarederror(double *odhadovanybod, double *poleznamychbodu, double invMatrixR[100][100], int size, int rad)
{
	int i, j;
	double sig, mezivysledek, citatel, jmenovatel, vektorr[100][100], vektorr2[100][100], vectorof1[100][100], 
	vectorof12[100][100], resultmatrix[100][100], resultmatrix2[100][100], boda[rad], bodb[rad];
	sig = sigma(poleznamychbodu,invMatrixR,size,rad);
	for (j = 0; j < rad; j++){
		boda[j] = *(odhadovanybod+j);
	}
	for (i = 0; i < size; i++){
		for (j = 0; j < rad; j++){
			bodb[j] = *(poleznamychbodu+(i*(rad+1))+j);
		}
		vektorr[0][i] = exp(-vzdalenost(&boda[0],&bodb[0],dtheta,dph,rad));		//vypocet korelacniho vektoru r'
		vektorr2[i][0] = exp(-vzdalenost(&boda[0],&bodb[0],dtheta,dph,rad));		//vypocet korelacniho vektoru r
	}
	for (j = 0; j < size; j++){
		vectorof1[0][j] = 1;										//vektor 1'
		vectorof12[j][0] = 1;										//vektor 1
	}
	multiply(vektorr,invMatrixR,resultmatrix,1,size,size,size);		// = r'R^-1
	multiply(resultmatrix,vektorr2,resultmatrix2,1,size,size,1);	// = r'R^-1 r
	mezivysledek = 1 - resultmatrix2[0][0];							// = 1 - r'R^-1 r
	multiply(vectorof1,invMatrixR,resultmatrix,1,size,size,size);	// = 1'R^-1
	multiply(resultmatrix,vektorr2,resultmatrix2,1,size,size,1);	// = 1'R^-1 r
	citatel = pow(1 - resultmatrix2[0][0],2);						// = (1 - 1'R^-1 r)^2
	multiply(vectorof1,invMatrixR,resultmatrix,1,size,size,size);	// = 1'R^-1
	multiply(resultmatrix,vectorof12,resultmatrix2,1,size,size,1);	// = 1'R^-1 1
	jmenovatel = resultmatrix2[0][0];
	return fabs(sig * (mezivysledek + (citatel/jmenovatel)));						// = sigma^2 (1 - r'R^-1 r + ((1 - 1'R^-1 r)^2)/(1'R^-1 1))
}
double expectedimprovement(double *odhadovanybod, double *poleznamychbodu, double invMatrixR[100][100], int size, int rad)
{
	int i, j;
	double PDF, CDF, odhady, MSE , sig, u, x, polebodu[100], boda[rad];
	for (j = 0; j < rad; j++){
		boda[j] = *(odhadovanybod+j);
	}
	
	for (i = 0; i < size; i++){
		polebodu[i]= *(poleznamychbodu+((i*(rad+1))+rad));
	}
	//u = mi(poleznamychbodu,invMatrixR,size,rad);
	qsort(&polebodu[0],size,sizeof(double),compare_ints);
	//sig = polebodu[0];
	odhady = predictor(odhadovanybod,poleznamychbodu,invMatrixR,size,rad);
	MSE = meansquarederror(odhadovanybod,poleznamychbodu,invMatrixR,size,rad);
	x = (polebodu[0]-odhady)/sqrt(MSE);
	PDF = (1.0/(sqrt(2*M_PI)))*(exp(-((pow(x,2))/2)));
	CDF = 0.5*(1+erf(x/sqrt(2)));
	//PDF = (1/(sqrt(1)*(sqrt(2*M_PI))))*(exp(-((pow(x,2))/(2*1))));
	//CDF = 0.5*(1+erf((x)/(sqrt(1)*sqrt(2))));
	u = ((polebodu[0]-odhady)*CDF+(sqrt(MSE)*PDF));
	if(u != u){
		return 0;
	}
	else{
		return u;
	}
}
double expectedimprovementBB(double *poleznamychbodu, double invMatrixR[100][100], double yprediktor,double s2chyba, int size, int rad)
{
	int i;
	double PDF, CDF, odhady, MSE , u, x, polebodu[100];
	
	for (i = 0; i < size; i++){
		polebodu[i]= *(poleznamychbodu+((i*(rad+1))+rad));
	}
	qsort(&polebodu[0],size,sizeof(double),compare_ints);
	
	odhady = yprediktor;
	MSE = s2chyba;
	x = (polebodu[0]-odhady)/sqrt(MSE);
	PDF = (1.0/(sqrt(2*M_PI)))*(exp(-((pow(x,2))/2)));
	CDF = 0.5*(1+erf(x/sqrt(2)));

	u = ((polebodu[0]-odhady)*CDF+(sqrt(MSE)*PDF));
	if(u != u){
		return 0;
	}
	else{
		return u;
	}
}
double polynom(double x)
{
return 5.95121023505635 + 0.887754463839314*x - 7.70721913502845*pow(x,2) + 13.8332014875226*pow(x,3) - 13.3092977965037*pow(x,4) + 7.36278255763972*pow(x,5) - 2.49262822525486*pow(x,6) + 0.537433122512059*pow(x,7) - 
0.0737094124205644*pow(x,8) + 0.00586951381472542*pow(x,9) - 0.000136327435548531*pow(x,10) - 0.0000238489156809096*pow(x,11) + 0.00000296118529300681*pow(x,12) - 0.000000160787297342575*pow(x,13) + 
0.00000000451685442306565*pow(x,14) - 0.0000000000533679890011691*pow(x,15);
}
void BBoptimalizace2D(double *vzorky,double *interval, int size, int rad)
{
	double mojeMaticeR[100][100], mojeInvMaticeR[100][100], A[4], B[4], C[4], D[4], miny=0,maxs2=0, intervaly[100], maxzlepseni=0, pozbod[3], odhady;
	int i,j,jemnost=15;
	//Vypocet matice R a inverzni matice R-1
	maticeR(vzorky,mojeMaticeR,size,rad);
	cofactors(mojeMaticeR,mojeInvMaticeR,size);
	
	//Urceni hranic jednotlivych intervalu
	A[0]=*interval;
	A[1]=((*(interval+1)-*interval)*fi)+*interval;
	A[2]=*(interval+2);
	A[3]=((*(interval+3)-*(interval+2))*fi)+*(interval+2);
	
	B[0]=((*(interval+1)-*interval)*fi)+*interval;
	B[1]=*(interval+1);
	B[2]=*(interval+2);
	B[3]=((*(interval+3)-*(interval+2))*fi)+*(interval+2);
	
	C[0]=((*(interval+1)-*interval)*fi)+*interval;
	C[1]=*(interval+1);
	C[2]=((*(interval+3)-*(interval+2))*fi)+*(interval+2);
	C[3]=*(interval+3);
	
	D[0]=*interval;
	D[1]=((*(interval+1)-*interval)*fi)+*interval;
	D[2]=((*(interval+3)-*(interval+2))*fi)+*(interval+2);
	D[3]=*(interval+3);
	
	//test intervalu A
	for (i = 0; i < jemnost; i++) {	
		for (j = 0; j < jemnost; j++) {
		pozbod[0] = A[0]+(((double)i / (jemnost-1))*(A[1]-A[0]));
			pozbod[1] = A[2]+(((double)j / (jemnost-1))*(A[3]-A[2]));
			odhady = predictor(&pozbod[0],vzorky,mojeInvMaticeR,size,2);
			if(i==0 && j==0){
				miny=odhady;
			}
			if(miny>odhady){
				miny=odhady;
			}
		}
	}
	for (i = 0; i < jemnost; i++) {	
		for (j = 0; j < jemnost; j++) {
			pozbod[0] = A[0]+(((double)i / (jemnost-1))*(A[1]-A[0]));
			pozbod[1] = A[2]+(((double)j / (jemnost-1))*(A[3]-A[2]));
			odhady = meansquarederror(&pozbod[0],vzorky,mojeInvMaticeR,size,2);
			if(i==0 && j==0){
				maxs2=odhady;
			}
			if(maxs2<odhady){
				maxs2=odhady;
			}
		}
	}
	for (i = 0; i < 4; i++){		//zmensuji interval na interval A a v dalsi casti porovnam zda B,C nebo D nejsou vhodnejsi
		*(interval+i)=A[i];
	}
	
	maxzlepseni=expectedimprovementBB(vzorky,mojeInvMaticeR,miny,maxs2,size,rad);//vypocet jake maximalni zlepseni muze nastat v tomto intervalu
	
		for (i = 0; i < jemnost; i++) {	
		for (j = 0; j < jemnost; j++) {
			pozbod[0] = B[0]+(((double)i / (jemnost-1))*(B[1]-B[0]));
			pozbod[1] = B[2]+(((double)j / (jemnost-1))*(B[3]-A[2]));
			odhady = predictor(&pozbod[0],vzorky,mojeInvMaticeR,size,2);
			if(i==0 && j==0){
				miny=odhady;
			}
			if(miny>odhady){
				miny=odhady;
			}
		}
	}
	for (i = 0; i < jemnost; i++) {	
		for (j = 0; j < jemnost; j++) {
			pozbod[0] = B[0]+(((double)i / (jemnost-1))*(B[1]-B[0]));
			pozbod[1] = B[2]+(((double)j / (jemnost-1))*(B[3]-B[2]));
			odhady = meansquarederror(&pozbod[0],vzorky,mojeInvMaticeR,size,2);
			if(i==0 && j==0){
				maxs2=odhady;
			}
			if(maxs2<odhady){
				maxs2=odhady;
			}
		}
	}
	if(expectedimprovementBB(vzorky,mojeInvMaticeR,miny,maxs2,size,rad)>maxzlepseni){		//pokud je pravdepobnost zlepseni v tomto intervalu vetsi tak jej vyberu jako vhodnejsi
		for (i = 0; i < 4; i++){		//vybiram vychozi interval jako B
		*(interval+i)=B[i];
		}
		maxzlepseni=expectedimprovementBB(vzorky,mojeInvMaticeR,miny,maxs2,size,rad);															
	}
	
	for (i = 0; i < jemnost; i++) {	
		for (j = 0; j < jemnost; j++) {
			pozbod[0] = C[0]+(((double)i / (jemnost-1))*(C[1]-C[0]));
			pozbod[1] = C[2]+(((double)j / (jemnost-1))*(C[3]-C[2]));
			odhady = predictor(&pozbod[0],vzorky,mojeInvMaticeR,size,2);
			if(i==0 && j==0){
				miny=odhady;
			}
			if(miny>odhady){
				miny=odhady;
			}
		}
	}
	for (i = 0; i < jemnost; i++) {	
		for (j = 0; j < jemnost; j++) {
			pozbod[0] = C[0]+(((double)i / (jemnost-1))*(C[1]-C[0]));
			pozbod[1] = C[2]+(((double)j / (jemnost-1))*(C[3]-C[2]));
			odhady = meansquarederror(&pozbod[0],vzorky,mojeInvMaticeR,size,2);
			if(i==0 && j==0){
				maxs2=odhady;
			}
			if(maxs2<odhady){
				maxs2=odhady;
			}
		}
	}
	if(expectedimprovementBB(vzorky,mojeInvMaticeR,miny,maxs2,size,rad)>maxzlepseni){		//pokud je pravdepobnost zlepseni v tomto intervalu vetsi tak jej vyberu jako vhodnejsi
		for (i = 0; i < 4; i++){		//vybiram vychozi interval jako C
		*(interval+i)=C[i];
		}
		maxzlepseni=expectedimprovementBB(vzorky,mojeInvMaticeR,miny,maxs2,size,rad);															
	}
	for (i = 0; i < jemnost; i++) {	
		for (j = 0; j < jemnost; j++) {
			pozbod[0] = D[0]+(((double)i / (jemnost-1))*(D[1]-D[0]));
			pozbod[1] = D[2]+(((double)j / (jemnost-1))*(D[3]-D[2]));
			odhady = predictor(&pozbod[0],vzorky,mojeInvMaticeR,size,2);
			if(i==0 && j==0){
				miny=odhady;
			}
			if(miny>odhady){
				miny=odhady;
			}
		}
	}
	for (i = 0; i < jemnost; i++) {	
		for (j = 0; j < jemnost; j++) {
			pozbod[0] = D[0]+(((double)i / (jemnost-1))*(D[1]-D[0]));
			pozbod[1] = D[2]+(((double)j / (jemnost-1))*(D[3]-D[2]));
			odhady = meansquarederror(&pozbod[0],vzorky,mojeInvMaticeR,size,2);
			if(i==0 && j==0){
				maxs2=odhady;
			}
			if(maxs2<odhady){
				maxs2=odhady;
			}
		}
	}
	if(expectedimprovementBB(vzorky,mojeInvMaticeR,miny,maxs2,size,rad)>maxzlepseni){		//pokud je pravdepobnost zlepseni v tomto intervalu vetsi tak jej vyberu jako vhodnejsi
		for (i = 0; i < 4; i++){		//vybiram vychozi interval jako D
		*(interval+i)=D[i];
		}
		maxzlepseni=expectedimprovementBB(vzorky,mojeInvMaticeR,miny,maxs2,size,rad);															
	}
	
}