// calculates \sum_K v_i'*A*v_i
double dosum(Matrix A, Matrix v)
{
  double alpha=0;
  int t = omp_get_max_threads();
  Matrix temp = createMatrix(A->rows,t);
#pragma omp parallel for schedule(static) reduction(+:alpha)
  for( int i=0;i<v->cols;++i ) {
    myMxV(temp->col[omp_get_thread_num()],A,v->col[i]);
    alpha += myinnerproduct(temp->col[omp_get_thread_num()],v->col[i]);
  }
  freeMatrix(temp);

  return alpha;
}

int readTwoOnDimMatrixFromAFile(char *fileName , oneDimMatrix * a , oneDimMatrix * b){
   FILE *fp;
   int numberOfLines = 0;

   fp = openAFileForReading(fileName);
   numberOfLines = readMatrixFromFile(fp, numberOfLines, DETECT_SIZE, a);
   printf("%d  %d %d \n" ,numberOfLines ,  a->nrows , a->ncols );
   numberOfLines= readMatrixFromFile(fp, numberOfLines, DETECT_SIZE , b);
   printf("%d  %d %d \n" ,numberOfLines , b->nrows , b->ncols);
   if(a->nrows != b->nrows)
      inconsistentErrorReport("Different number of rows in each matrix.", 
                              numberOfLines, a->nrows , a->ncols);
   fclose(fp);

   *a = createMatrix(a->nrows , a->ncols);
   *b = createMatrix(b->nrows , b->ncols);
   fp = openAFileForReading(fileName);
   numberOfLines = readMatrixFromFile(fp, numberOfLines, READ_MATRIX, a);
   numberOfLines= readMatrixFromFile(fp, numberOfLines, READ_MATRIX , b);
   fclose(fp);   
   /*   printOneDimMatrix(a);  printf("\n");  printOneDimMatrix(b); */
   return numberOfLines;
}

// creates an identity matrix with numRow rows and numCol columns; works similar to MATLABs eye()
void eye(struct Matrix *eyeMatrix,int numRow,int numCol) {
  createMatrix(eyeMatrix,numRow,numCol);                // create matrix eyeMatrix
  int i,j;
  for (i=0;i<numRow;i++) {                        // iterate through rows
    for (j=0;j<numCol;j++) {                      // iterate through columns
      if (i==j) {
        eyeMatrix->value[i][j]=1.0;                   // if on the diagonal set to 1
      }
      else {
        eyeMatrix->value[i][j]=0.0;                   // else set value to zero
      }
    }
  }
}

int main(int argc, char** argv)
{
  int rank, size;
  init_app(argc, argv, &rank, &size);

  if (argc < 2) {
    printf("usage: %s <N> [L]\n",argv[0]);
    close_app();
    return 1;
  }

  /* the total number of grid points in each spatial direction is (N+1) */
  /* the total number of degrees-of-freedom in each spatial direction is (N-1) */
  int N  = atoi(argv[1]);
  int M  = N-1;
  double L=1.0;
  if (argc > 2)
    L = atof(argv[2]);

  double h = L/N;
  poisson_info_t ctx;
  ctx.A = createPoisson1D(M);

  Vector grid = createVector(M);
  for (int i=0;i<M;++i)
    grid->data[i] = (i+1)*h;

  Matrix u = createMatrix(M, M);
  evalMesh(u->as_vec, grid, grid, poisson_source);
  scaleVector(u->as_vec, h*h);

  double time = WallTime();
  cg(evaluate, u, 1.e-6, &ctx);

  evalMesh2(u->as_vec, grid, grid, exact_solution, -1.0);
  double max = maxNorm(u->as_vec);

  if (rank == 0) {
    printf("elapsed: %f\n", WallTime()-time);
    printf("max: %f\n", max);
  }

  freeMatrix(u);
  freeVector(grid);
  freeMatrix(ctx.A);

  close_app();
  return 0;
}

void printEigen( Matrix *eigen_values, Matrix *eigen_vectors, int num_of_file ) {
	char fn[256];

	Matrix *eigen = createMatrix( eigen_vectors->height + 1, eigen_values->height );
	for ( int row = 0; row < eigen->height; row++ ) {
		eigen->a[row][0] = eigen_values->a[row][row];
		for ( int col = 1; col < eigen->width; col++ ) {
			eigen->a[row][col] = eigen_vectors->a[col-1][row];
		}
	}

	sort( eigen );
	sprintf( fn, "eigen%02d.txt", num_of_file );
	printMatrixToFile( fn, eigen );
}

static ZSolveMatrix Matrix2zsolve(Matrix *M)
{
    int i, j;
    ZSolveMatrix zmatrix;

    zmatrix = createMatrix(M->NbColumns-2, M->NbRows);
    for (i = 0; i < M->NbRows; ++i)
	for (j = 0; j < M->NbColumns-2; ++j) {
	    assert(mpz_cmp_si(M->p[i][1+j], -MAXINT) > 0);
	    assert(mpz_cmp_si(M->p[i][1+j], MAXINT) < 0);
	    zmatrix->Data[i*zmatrix->Width+j] = mpz_get_si(M->p[i][1+j]);
	}

    return zmatrix;
}

Matrix createPoisson1D(int M)
{
  Matrix result = createMatrix(M, M);
  result->data[0][0] = 2.0;
  result->data[1][0] = -1.0;
  for (int i=1;i<M-1;++i) {
    result->data[i][i] = 2.0;
    result->data[i+1][i] = -1.0;
    result->data[i-1][i] = -1.0;
  }
  result->data[M-1][M-1] = 2.0;
  result->data[M-2][M-1] = -1.0;

  return result;
}

SmithWatermanDP::SmithWatermanDP(const char* s1, size_t n1, const char* s2, size_t n2) 
  // recall that the SW matrix has n+1 and m+1 rows and columns respectively
  : DynamicProgramming(n1+1,n2+1)
{
  this->x = s1;
  this->y = s2;
  // default gap penalty is 1 and efault sim matrix is 1 on
  // diagonal and 0 elsewere
  this->gapPenalty = 1;
  this->sim = createDefaultSimilarityMatrix(256);
  // defualt value for backtracking is not enbaled
  this->btEnabled = false;
  this->btMatrix = NULL;
  createMatrix();
}

int main(int argc, char** argv)
{
  int N, K, i, j;
  Matrix A,v;
  double time, sum;

  if (argc < 3) {
    printf("need two parameters, the matrix size and the number of vectors\n");
    return 1;
  }

  N=atoi(argv[1]);
  K=atoi(argv[2]);

  A = createMatrix(N,N);
  // identity matrix
  for (i=0;i<N;++i)
    A->data[i][i] = 1.0;

  v = createMatrix(N,K);

  // fill with column number
  for (i=0;i<K;++i)
    for (j=0;j<N;++j)
      v->data[i][j] = i;

  time = WallTime();
  sum = dosum(A,v);

  printf("sum: %f\n", sum);
  printf("elapsed: %f\n", WallTime()-time);
  freeMatrix(v);
  freeMatrix(A);

  return 0;
}

int main()
{
    FILE * f = fopen("matrix.txt","r");
    readFromAdjMatrix(f);
    printAdjMatrix();
    NodeT** v=(NodeT**)malloc(nrOfVerteces*sizeof(NodeT));
    v=CreateList(adjMatrix);
    int **matrix=(int**)malloc(sizeof(int*)*nrOfVerteces);
    matrix=createMatrix(v);
    printAllList(v);
    printMatrix(matrix);
    bfs(0);
    dfs(0);
    dfsRecurs(0);
    return 0;
}

// sharpening image - apply laplacian filter and
//subtract the resultant image from the original image
struct matrix* imsharpen(struct matrix* image)
{
	struct matrix* filter;
	filter = createMatrix(3,3);

	initMatrix(filter, (void*)laplacianFilter);

	struct matrix* filteredImage = imfilter(image,filter,CORRELATION_OPERATION,MODE_PAD_ZERO);

	struct matrix* sharpenedImage = subtractMatrices(image,filteredImage);

	destroyMatrix(filter);
	destroyMatrix(filteredImage);

	return sharpenedImage;
}

double HMM::backwardAlgorithm(vector<int> words)
{
	int len = words.size();
	double ** backward = createMatrix(len, _maxState);
	computeBackwardMatrix(words, backward, len);

	double p = 0;
	for (int i = 0; i < _maxState; i++)
	{
		p += backward[0][i] * _pObservation[i][words[0]] * initial_probability[i];
	}

	freeMatrix(backward, len, _maxState);

	return p;
}

CvMatND* FeatPyramid::padArray (CvMatND *mat, int dimPad[3], float val)
{
  int dims[3];

  if (mat == NULL)
  {
    for (int i = 0; i < 2; i++)
      dims[i] = dimPad[i] * 2;

    dims[2] = 32 + (dimPad[2] * 2);
  }

  else
  {
    for (int i = 0; i < 3; i++)
      dims[i] = mat->dim[i].size + (dimPad[i] * 2);
  }

  // New bigger matrix is created
  assert (dims[0] > 0);
  assert (dims[1] > 0);
  assert (dims[2] > 0);

  //CvMatND *padMat = createNDMatrix (3, dims, CV_64FC1);
	CvMatND *padMat;
	createMatrix (3, dims, CV_64FC1, &padMat);
  assert (padMat != NULL);

  // Firstable, all the elements are setted to val
  for ( int i = 0; i < padMat->dim[0].size; i++ )
    for ( int j = 0; j < padMat->dim[1].size; j++ )
      for ( int k = 0; k < padMat->dim[2].size; k++ )
        cvSetReal3D (padMat, i, j, k, val);

  if (mat != NULL)
  {
    // Original values are setted
    for ( int i = dimPad[0]; i < padMat->dim[0].size - dimPad[0]; i++ )
      for ( int j = dimPad[1]; j < padMat->dim[1].size - dimPad[1]; j++ )
        for ( int k = dimPad[2]; k < padMat->dim[2].size - dimPad[2]; k++ )
          cvSetReal3D (padMat, i, j, k, cvGetReal3D (mat, i-dimPad[0], j-dimPad[1], k-dimPad[2]));
  }

  return padMat;
}

NS_CORE InternalGearRatios DefKNuton::solveNuton( const Jacobi& jacobian, System& system )
{
	const double eps = 0.001;
	const int maxIterCount = 100;
	double norm;
	int iterCount = 0;
	bool notFinded = false;

	do
	{
		// создаем матрицу с уравнениями
		auto matrix = createMatrix( jacobian, system );
		// создаем вектор правых частей уравнений
		auto rightParts = createRightParts( system );
		// решаем систему уравнений и получаем дельту
		auto next = MatrixOperations::solveGaus( matrix, rightParts );
		if ( next.size() == 0 )
		{
			next = MatrixOperations::solveKramer( matrix, rightParts );
		}
		// если решений нет - прерываемся, иначе вычисляем начальные условия для следующей итерации
		if ( next.size() == 0 || ++iterCount > maxIterCount )
		{
			notFinded = true;
		}
		else
		{
			// вычисляем норму ( максимальную дельту )
			norm = calcNorm( next );
			// вычисляем начальные значения для следующей итерации
			for ( size_t i = 0; i < next.size(); i++ )
			{
				system.getUnknownVariables()[i].setValue( system.getUnknownVariables()[i].getValue() + next[i] );
			}
		}
	} while ( !notFinded && norm >= eps );

	NS_CORE InternalGearRatios ans;
	if ( !notFinded )
	{
		ans = geInternalGearRatioValuesFromSystem( system );
	}

	return ans;
}

int main(int argc, char** argv)
{
  int i, j, N, flag;
  Vector grid;
  Matrix b, e;
  double time, sum, h;
  int rank, size;
  int mpi_top_coords[2];
  int mpi_top_sizes[2];

  init_app(argc, argv, &rank, &size);

  N=atoi(argv[1]);

  // setup topology
  mpi_top_sizes[0] = mpi_top_sizes[1] = 0;
  MPI_Dims_create(size, 2, mpi_top_sizes);
  int periodic[2] = {0, 0};
  MPI_Comm comm;
  MPI_Cart_create(MPI_COMM_WORLD, 2, mpi_top_sizes, periodic, 0, &comm);
  MPI_Cart_coords(comm, rank, 2, mpi_top_coords);

  int* size1;
  int* displ1;
  int* size2;
  int* displ2;
  splitVector(N, mpi_top_sizes[0], &size1, &displ1);
  splitVector(N, mpi_top_sizes[1], &size2, &displ2);

  b = createMatrix(size1[mpi_top_coords[0]], size2[mpi_top_coords[1]]);
  for (j=0;j<b->cols;++j)
    for(i=0;i<b->rows;++i)
      b->data[j][i] = (j+displ2[mpi_top_coords[1]])*N+1+(i+displ1[mpi_top_coords[0]]);
  b->glob_rows = N;
  b->glob_cols = N;
  b->as_vec->comm = &comm;

  saveMatrixMPI(b, "meh.asc");

  freeMatrix(b);
  MPI_Comm_free(&comm);

  close_app();
  return 0;
}

Matrix* processThroughReceptiveFields(Matrix *ts, const double *centers, const Constants *c) {
    Matrix *ts_out = createMatrix(c->M, ts->ncol);
//    printf("%d %d %d", c->M, ts->nrow, c->M % ts->nrow);
    assert(c->M % ts->nrow == 0);
    
    size_t neurons_per_ts = c->M / ts->nrow;
    size_t i, ri;
    for(i=0, ri=0; (i < c->M) && (ri < ts->nrow); i+=neurons_per_ts, ri++) {
        for(size_t ni=i; ni<(i+neurons_per_ts); ni++) {
            for(size_t vi=0; vi < ts->ncol; vi++) {
                double val = getMatrixElement(ts, ri, vi);                
                double val_field = c->preproc->gain * exp( - ((centers[ni] - val) * (centers[ni] - val))/c->preproc->sigma );
                setMatrixElement(ts_out, ni, vi, val_field);
            }
        }
    }
    return(ts_out);
}

// calculates \sum_K v_i'*A*v_i
double dosum(Matrix A, Matrix v)
{
  double alpha=0.0;
  Matrix temp;
  int i, t;

  t = getMaxThreads();

  temp = createMatrix(A->rows, t);
#pragma omp parallel for schedule(static) reduction(+:alpha)
  for(i=0;i<v->cols;++i) {
    MxV(temp->col[getCurrentThread()],A,v->col[i], 1.0, 0.0, 'N');
    alpha += dotproduct(temp->col[getCurrentThread()],v->col[i]);
  }
  freeMatrix(temp);

  return alpha;
}

void Image::readImage(const QString & path)
{
    QFile file(path);

    if (file.open(QIODevice::ReadOnly))
    {
        QTextStream stream(&file);
        QString line = "";

        for(int i=0;i<3;i++)
        {
            line = stream.readLine();
            if(!line.startsWith("#"))
            {
                if(i==0)
                        format = line;
                if(i==1)
                {
                    QStringList list = line.split(" ");
                    columns=list[0].toInt();
                    rows=list[1].toInt();
                    if(QString::compare(format, "P3") == 0 || QString::compare(format, "P6") == 0)
                        width=columns*3;
                    else
                        width=columns;
                }
                if(i==2)
                    levels = line.toInt();
            }
            else
                i--;
        }

        createLut();
        createMatrix();

        for (int idx=0; idx<rows; idx++ )
        {
            for (int idy=0; idy<width; idy++ )
            {
                stream>>matrix[idx][idy];
            }
        }
    }

calc_data do_cucl(int k, int n, double m, double g, QVector<double> trend_m, QVector<double> trend_s,
                  double level)
{
    QVector<QVector<double>> matrix = createMatrix(k,n,m,g,trend_m,trend_s);
    QVector<double> y;
    for(auto sample : matrix)
    {
        y.push_back(sample_mean(sample));
    }
    QVector<double> x;
    for(int i = 0; i<k;i++)
    {
        x.push_back(i);
    }

    QVector<bool> u = U(y);
    QVector<bool> l = L(y);

    int a = Asum(u,l);
    int b = Bsum(u,l);

    //Тренд средних
    double F = f(k);
    double t_stat_mean = a/F;

    double L = l_const(k);
    double t_stat_s = (b-F*F)/L;
    double t_cr = t_test_cr(k,level);
    qDebug() << "a" << a;
    qDebug() << "b" << b;
    qDebug() << "t mean" << t_stat_mean;
    qDebug() << "t s" << t_stat_s;
    qDebug() << t_cr;
    calc_data data;
    data.x=x;
    data.y=y;
    data.t_A=t_stat_mean;
    data.t_B=t_stat_s;
    data.t_cr=t_cr;
    data.A=a;
    data.B=b;
    return data;
}

Matrix *multiply( Matrix *a, Matrix *b ) {
	Matrix *product = createMatrix( a->height, b->width );

	if ( a->width != b->height ) {
		printf("計算出来ない行列の組みわせです．\n");
		exit(1);
	}

	for ( int row = 0; row < a->height; row++ ) {
		for ( int col = 0; col < b->width; col++ ) {
			product->a[ row ][ col ] = 0;
			for ( int k = 0; k < a->width; k++ ) {
				product->a[ row ][ col ] += a->a[ row ][ k ] * b->a[ k ][ col ];
			}
		}
	}

	return product;
}