#include "matrix.h"

matrix *matrix_alloc (unsigned int M, unsigned int N) {
#ifdef MATRIX_DEBUG
  if ((M <= 0) | (N <= 0)) {
    printf("%s:%i: warning: Alokovana matice musi mit nenulovy pocet radku i sloupcu (M = %i , N = %i)\n", __FILE__, __LINE__, M, N);
    exit(1);
  }
#endif
  matrix * mat = malloc(sizeof (matrix));
  mat->m = M;
  mat->n = N;
  mat->data = (double **) malloc(sizeof (double *) * M);
  unsigned int i;
  for (i = 0; i < M; i++) {
    mat->data[i] = (double *) malloc(sizeof (double) * N);
  }
  return mat;
}


// Uvolni matici MAT z pameti.

void matrix_free (matrix * mat) {
  unsigned int i;
  for (i = 0; i < mat->m; i++) {
    free((* mat).data[i]);
  }
  free((* mat).data);
  free(mat);
}

unsigned long matrix_save_length (matrix *mat) {
  unsigned int m, n, L;
  m = matrix_size_m(mat);
  n = matrix_size_n(mat);
  L = m * n * (2 * 8 + 3) + 2 * (8 + 1) + 1; //printf("L = %i\n");
  return L;
}

char *
matrix_save_save_string (matrix *mat) {
  char buffer[100];
  char *str;
  unsigned int i, j, m, n, L;
  L = matrix_save_length(mat);
  str = malloc(sizeof (char) *L);
  for (i = 0; i < L; i++) {
    str[i] = NULL;
  }
  m = matrix_size_m(mat);
  n = matrix_size_n(mat);
  sprintf(buffer, "%i %i", m, n);
  strcat(str, buffer);
  for (i = 0; i < m; i++) {
    for (j = 0; j < n; j++) {
      sprintf(buffer, " %8.8f", matrix_get(mat, i, j));
      strcat(str, buffer);
    }
  }
  //printf("MATICE\n%s",str);

  return str;
}

int
matrix_save(matrix *mat, char *file)
{
  printf("Saving : %s\n", file);
  FILE * pFile;
  pFile = fopen (file,"w");
  char *str = matrix_save_save_string(mat);
  int out = fputs(str, pFile);
  free(str);
  return out;
}

matrix *
matrix_load_string (char *string) {
  char *str = string;
  unsigned int m = atoi(str);
  while (str[0] != (" ")[0]) {
    str++;
  };
  str++;
  unsigned int n = atoi(str); //printf("%i %i\n", m, n);
  while (str[0] != (" ")[0]) {
    str++;
  };
  str++;
  unsigned int i, j;
  matrix *A = matrix_alloc(m, n);
  double value;
  for (i = 0; i < m; i++) {
    for (j = 0; j < n; j++) {
      value = atof(str); //printf("%f\n", value);
      matrix_set(A, i, j, value);
      if (!((i == m - 1)&(j == n - 1))) {
        while (str[0] != (" ")[0]) {
          str++;
        };
        str++;
      }
    }
  }
  return A;
}

void
matrix_move_m (matrix *mat, int down) {
  unsigned int i, M;
  M = matrix_size_m(mat);
  double **new_data = (double**) malloc(sizeof (double*) *M);
  double **swap;
  down = down % M;
  if (down < 0) {
    down += M;
  }
  for (i = 0; i < down; i++) {
    new_data[i] = mat->data[M - down + i];
  }
  for (i = down; i < M; i++) {
    new_data[i] = mat->data[i - down];
  }
  free(mat->data);
  mat->data = new_data;

}

matrix *
matrix_clone (matrix * mat) {
  unsigned int i, j;
  matrix * clone = matrix_alloc(matrix_size_m(mat), matrix_size_n(mat));
  for (i = 0; i < mat->m; i++) {
    for (j = 0; j < mat->n; j++) {
      matrix_set(clone, i, j, matrix_get(mat, i, j));
    }
  }

  return clone;
}

// Vrati pocet radku matice MAT.

unsigned int
matrix_size_m (matrix * mat) {
  return mat->m;
}


// Vrati pocet sloupcu matice MAT.

unsigned int
matrix_size_n (matrix * mat) {
  return mat->n;
}


// Nastavi matici MAT na indexech I, J hodnotu VALUE.

void
matrix_set (matrix * mat, unsigned int i, unsigned int j, double value) {
#ifdef MATRIX_DEBUG
  if ((i < matrix_size_m(mat))&(j < matrix_size_n(mat))) {
#endif
    mat->data[i][j] = value;
#ifdef MATRIX_DEBUG
  } else {
    fprintf(stderr, "%s: %i: error: Pri nastavovani hodnoty matice doslo k chybe. Spatne zadany index M[%i][%i] mimo rozmer matice %i x %i\n", __FILE__, __LINE__, i + 1, j + 1, matrix_size_m(mat), matrix_size_n(mat));
    exit(1);
  }
#endif
}


// Nastavi matici MAT na indexech I, J hodnotu VALUE.

double
matrix_get (matrix * mat, unsigned int i, unsigned int j) {
#ifdef MATRIX_DEBUG
  if ((i < matrix_size_m(mat))&(j < matrix_size_n(mat))) {
#endif
    return mat->data[i][j];
#ifdef MATRIX_DEBUG
  } else {
    fprintf(stderr, "%s: %i: error: Pri ziskavani hodnoty z matice doslo k chybe. Spatne zadany index M[%i][%i] mimo rozmer matice %i x %i\n", __FILE__, __LINE__, i, j, matrix_size_m(mat), matrix_size_n(mat));
    exit(1);
  }
#endif
}

// Vrati  M-ty az N-ty radek a I-ty az J-ty sloupce z matice MAT.

matrix *
matrix_submatrix (matrix *mat, unsigned int m, unsigned int n,
        unsigned int i, unsigned int j) {
#ifdef MATRIX_DEBUG
  if ((i < matrix_size_n(mat))&(j < matrix_size_n(mat))&(i <= j)&
          (m < matrix_size_m(mat))&(n < matrix_size_m(mat))&(m <= n)) {
#endif

    unsigned int M, I, matI;
    matrix *sub = matrix_alloc(n - m + 1, j - i + 1);
    for (M = 0; M < matrix_size_m(sub); M++) {
      matI = i;
      for (I = 0; I < matrix_size_n(sub); I++) {
        matrix_set(sub, M, I, matrix_get(mat, m, matI++));
      }
      m++;
    }
    return sub;
#ifdef MATRIX_DEBUG
  } else {
    printf("%s: %i: error: Pri ziskavani submatice doslo k chybe. Zadane hranice jsou (%i az %i radek) x (%i az %i sloupec) mimo rozmer matice %i x %i\n", __FILE__, __LINE__, m + 1, n + 1, i + 1, j + 1, matrix_size_m(mat), matrix_size_n(mat));
    exit(1);
  }
#endif
}

matrix *
matrix_rand (unsigned int m, unsigned int n) {
  matrix *mat = matrix_alloc(m, n);
  unsigned int i, j;
  double *row;
  for (i = 0; i < matrix_size_m(mat); i++) {
    row = (*mat).data[i];
    for (j = 0; j < matrix_size_n(mat); j++) {
      row[j] = ((double) rand()) / ((double) RAND_MAX + 1);
    }
  }
  return mat;
}

matrix *
matrix_rand_int (unsigned int m, unsigned int n) {
  unsigned int i, j;
  matrix *rand_mat = matrix_rand(m, n);
  matrix *rand_int_mat = matrix_alloc(m, n);
  for (i = 0; i < m; i++) {
    for (j = 0; j < n; j++) {
      matrix_set(rand_int_mat, i, j, floor(matrix_get(rand_mat, i, j) * m * n));
    }
  }
  matrix_free(rand_mat);
  return rand_int_mat;
}

matrix *
matrix_rand_int_no_repeating (unsigned int m, unsigned int n) {
  unsigned int i, j, I, J, R, corect;
  // Matice nahodnych prvku < 0 , 1 )
  matrix *rand_mat = matrix_rand(m, n);
  // matice vystupnich nahodnych intu
  matrix *rand_int_mat = matrix_alloc(m, n);
  // Matice nahodnych intu, kde se mohou vyskytovat i vicektrat stejna cisla
  matrix *RN = matrix_alloc(m, n);

  // Kazdemu prvku z matice RAND_MAT vypocita kolikaty je v poradi nejvyssi.
  double min;
  for (i = 0; i < m; i++) {
    for (j = 0; j < n; j++) {
      min = matrix_get(rand_mat, i, j);
      R = 0;
      for (I = 0; I < m; I++) {
        for (J = 0; J < n; J++) {
          if ((min > matrix_get(rand_mat, I, J))) {
            R++;
          }
        }
      }
      matrix_set(RN, i, j, R);

      // Osetreni jestli uz dana hodnota R je v matici RN jednou nastavena.
      // Pokud ano tak se spocte kolikrat a tohle cislo se odecte od R.
      //
      // Pokud by bylo na 3 pozicich R = 4 tak u tretiho vyskytu by se hodnota
      // prepsala ve vystupni matici na 6, u druheho 5 a u prvniho by zustala 4.
      corect = 0;
      for (I = 0; I < i; I++) {
        for (J = 0; J < j; J++) {
          if (matrix_get(RN, I, J) == R) {
            corect++;
          }
        }
      }
      // Ulozi hodnotu do vystupni matice.
      matrix_set(rand_int_mat, i, j, R + corect);
    }
  }

  // Uvolneni pomocnych promenych
  matrix_free(rand_mat);
  matrix_free(RN);
  // Matice intu bez opakovanis hodnotami < 0 , M * N - 1 )
  return rand_int_mat;

}

matrix *
matrix_reshuffle_impossible (matrix *mat, matrix *impossible) {

#ifdef MATRIX_DEBUG
  if ((matrix_size_m(mat) != (matrix_size_m(impossible)))&
          ((matrix_size_n(mat)) != (matrix_size_n(impossible)))) {
    printf("%s:%i: error: Matice musi mit stejne rozmery (%i x %i) != (%i x %i)\n", matrix_size_m(mat), matrix_size_n(mat), matrix_size_m(impossible), matrix_size_n(impossible));
  }
#endif

  unsigned int i, j, I, J, stop, X;

  matrix *reshuffled_mat = matrix_clone(mat);

  double swap;
  X = 0;
  for (i = 0; i < matrix_size_m(reshuffled_mat); i++) {
    for (j = 0; j < matrix_size_n(reshuffled_mat); j++) {

      if (matrix_get(reshuffled_mat, i, j) == matrix_get(impossible, i, j)) {

        stop = 0;
        for (I = 0; I < matrix_size_m(reshuffled_mat); I++) {

          if (stop == 1) {
            break;
          }
          for (J = 0; J < matrix_size_n(reshuffled_mat); J++) {
            //printf("%i %i\n", I, J);
            if ((matrix_get(reshuffled_mat, I, J) != matrix_get(impossible, i, j))&
                    (matrix_get(reshuffled_mat, i, j) != matrix_get(impossible, I, J))) {
              swap = (matrix_get(reshuffled_mat, I, J));
              matrix_set(reshuffled_mat, I, J, matrix_get(reshuffled_mat, i, j));
              matrix_set(reshuffled_mat, i, j, swap);
              stop = 1;
              break;
            }
          }
        }
        if (stop == 0) {
#ifdef MATRIX_DEBUG
          printf("%s:%i: warning: Nepodarilo se preskladat matici !", __FILE__, __LINE__);
#endif
          return NULL;
        }
      }
    }
  }

  return reshuffled_mat;
}

matrix *
matrix_linspace (unsigned int m, unsigned int n, double a, double b, int how) {
  unsigned int i, j;
  double d;
  matrix *linspace;

  if (how == MATRIX_LINSPACE_COLUMN) {

    linspace = matrix_alloc(m, n);

    for (i = 0; i < m; i++) {
      d = (b - a) / (m - 1) * i + a;
      for (j = 0; j < n; j++) {
        matrix_set(linspace, i, j, d);
      }
    }

    return linspace;

  }

  if (how == MATRIX_LINSPACE_ROW) {

    linspace = matrix_alloc(m, n);

    for (j = 0; j < n; j++) {
      d = (b - a) / (n - 1) * j + a;
      for (i = 0; i < m; i++) {
        matrix_set(linspace, i, j, d);
      }
    }

    return linspace;

  }
  return NULL;
}

// Vynuluje matici.

void
matrix_zeros (matrix *mat) {
  unsigned int i, j;
  double *row;
  for (i = 0; i < matrix_size_m(mat); i++) {
    row = (*mat).data[i];
    for (j = 0; j < matrix_size_n(mat); j++) {
      row[j] = 0;
    }
  }
}

// Vrati 1 pokud maji matice A a B stejne hodnoty a 0 pokud ne.

int
matrix_equal (matrix * A, matrix * B) {
  unsigned int ma, na, mb, nb, i, j;
  // Zjisteni rozmeru matic.
  ma = matrix_size_m(A);
  na = matrix_size_n(A);
  mb = matrix_size_m(B);
  nb = matrix_size_n(B);
  // Pokud maji matice A, B ruzne rozmery tak se vypise chyba a vrati se 0.
  if ((ma != mb)&(na != nb)) {
    printf("%s:%s: Matice A a B maji jine rozmery", (const char *) __FILE__, (const char *) __LINE__);
    return 0;
  }

  for (i = 0; i < ma; i++) {
    for (j = 0; j < na; j++) {
      // Pokud A[i][j] != B[i][j] tak funkce vrati 0;
      if (matrix_get(A, i, j) != matrix_get(B, i, j)) {
        return 0;
      }
    }
  }
  return 1;
}


// Zkontroluje jestli ma matice nulove rozmery.

int
matrix_not_null (matrix * mat) {
  if ((mat->m == 0)&(mat->n == 0)) {
    return 1;
  }
  return 0;
}

// Vypise na obrazovku matici.

void
matrix_print (matrix * mat) {
  unsigned int m = matrix_size_m(mat);
  unsigned int n = matrix_size_n(mat);
  unsigned int i, j;
  for (i = 0; i < m; i++) {
    for (j = 0; j < n; j++) {
      printf("%.5f ", matrix_get(mat, i, j));
    }
    printf("\n");
  }
}

void
matrix_print_info (matrix * mat) {
  printf("matrix %i x %i\n", mat->m, mat->n);
}

matrix *
matrix_transpose (matrix *mat) {
  matrix *trans = matrix_alloc(matrix_size_n(mat), matrix_size_m(mat));
  unsigned int i, j;
  for (i = 0; i < matrix_size_m(trans); i++) {
    for (j = 0; j < matrix_size_n(trans); j++) {
      matrix_set(trans, i, j, matrix_get(mat, j, i));
    }
  }
  return trans;
}

double
matrix_forbenius_norm (matrix *mat) {
  unsigned int i, j;
  double N = 0;
  for (i = 0; i < matrix_size_m(mat); i++) {
    for (j = 0; j < matrix_size_n(mat); j++) {
      N += pow(matrix_get(mat, i, j), 2);
    }
  }
  return sqrt(N);
}

matrix *
matrix_multiply (matrix *A, matrix *B) {
  unsigned int i, j, k;

#ifdef MATRIX_DEBUG
  if ((matrix_size_m(A) != matrix_size_n(B)) | (matrix_size_n(A) != matrix_size_m(B))) {
    printf("%s:%i: Matice nelze nasobit, protoze maji spatne rozmery A (%i x %i) , B (%i , %i)\n", __FILE__, __LINE__, matrix_size_m(A), matrix_size_n(A), matrix_size_m(B), matrix_size_n(B));
    return NULL;
  }
#endif

  matrix *C = matrix_alloc(matrix_size_m(A), matrix_size_n(B));

  double prod;
  for (i = 0; i < matrix_size_m(C); i++) {
    for (j = 0; j < matrix_size_n(C); j++) {

      prod = 0;
      for (k = 0; k < matrix_size_n(A); k++) {
        prod += matrix_get(A, i, k) * matrix_get(B, k, j);
      }
      
      matrix_set(C, i, j, prod);
    }

  }
  return C;
}

matrix *
matrix_multiply_constant (matrix *A, double c) {
  unsigned int i, j;
  matrix *B = matrix_alloc(matrix_size_m(A), matrix_size_n(A));

  for (i = 0; i < matrix_size_m(A); i++) {
    for (j = 0; j < matrix_size_n(A); j++) {
      matrix_set(B, i, j, matrix_get(A, i, j) * c);
    }
  }

  return B;
}

matrix *
matrix_adition (matrix *A, matrix *B) {
  unsigned int i, j;
  matrix *C = matrix_alloc(matrix_size_m(A), matrix_size_n(A));

  for (i = 0; i < matrix_size_m(A); i++) {
    for (j = 0; j < matrix_size_n(A); j++) {
      matrix_set(C, i, j, matrix_get(A, i, j) + matrix_get(B, i, j));
    }
  }

  return C;
}

matrix *
matrix_substraction (matrix *A, matrix *B) {
  unsigned int i, j;
  matrix *C = matrix_alloc(matrix_size_m(A), matrix_size_n(A));

  for (i = 0; i < matrix_size_m(A); i++) {
    for (j = 0; j < matrix_size_n(A); j++) {
      matrix_set(C, i, j, matrix_get(A, i, j) - matrix_get(B, i, j));
    }
  }

  return C;
}

double
matrix_sum (matrix *mat) {
  unsigned int i, j;
  double sum = 0; // Inicializace souctu
  // Nyni sectu vsechny prvky
  for (i = 0; i < matrix_size_m(mat); i++) {
    for (j = 0; j < matrix_size_n(mat); j++) {
      sum += matrix_get(mat, i, j); // Pricteni prvku MAT[I][J]
    }
  }
  return sum;
}

matrix *
matrix_sum_horizontaly (matrix *mat) {
  unsigned int i, j;
  matrix *out = matrix_alloc(matrix_size_m(mat), 1);
  // Nyni sectu vsechny prvky
  for (i = 0; i < matrix_size_m(mat); i++) {
    double sum = 0; // Inicializace souctu radku
    for (j = 0; j < matrix_size_n(mat); j++) {
      sum += matrix_get(mat, i, j); // Pricteni prvku MAT[I][J]
    }
    matrix_set(out, i, 0, sum);
  }
  return out;
}

matrix *
matrix_sum_verticaly (matrix *mat) {
  unsigned int i, j;
  matrix *out = matrix_alloc(1, matrix_size_n(mat));
  // Nyni sectu vsechny prvky
  for (j = 0; j < matrix_size_n(mat); j++) {
    double sum = 0; // Inicializace souctu radku
    for (i = 0; i < matrix_size_m(mat); i++) {
      sum += matrix_get(mat, i, j); // Pricteni prvku MAT[I][J]
    }
    matrix_set(out, 0, j, sum);
  }
  return out;
}

double
matrix_mean (matrix *mat) {
  return matrix_sum(mat) / ((double) (matrix_size_m(mat) * matrix_size_n(mat)));
}


matrix *matrix_row_mean(matrix *mat)
{
  unsigned int i, j;
  double mean;
  matrix *out = matrix_alloc(matrix_size_m(mat), matrix_size_n(mat));
  for (i = 0; i<matrix_size_m(mat); i++) {
    mean = 0;
    for (j = 0; j<matrix_size_n(mat); j++) {
      mean += matrix_get(mat, i, j);
    }
    matrix_set(out, i, 0, mean/((double) matrix_size_n(mat)));
  }
  return out;
}

matrix *matrix_row_stdev(matrix *mat)
{
  unsigned int i, j;
  matrix *mean = matrix_row_mean(mat);
  double stdev;
  matrix *out = matrix_alloc(matrix_size_m(mat), matrix_size_n(mat));
  for (i = 0; i<matrix_size_m(mat); i++) {
    stdev = 0;
    for (j = 0; j<matrix_size_n(mat); j++) {
      stdev += pow(matrix_get(mat, i, j)-matrix_get(mean, i, 0), 2);
    }
    matrix_set(out, i, 0, sqrt(stdev)/((double) matrix_size_n(mat)));
  }
  matrix_free(mean);
  return out;
}

matrix *matrix_row_max(matrix *mat)
{
  unsigned int i, j;
  double max, foo;
  matrix *out = matrix_alloc(matrix_size_m(mat), matrix_size_n(mat));
  for (i = 0; i<matrix_size_m(mat); i++) {
    max = matrix_get(mat, i, 0);
    for (j = 1; j<matrix_size_n(mat); j++) {
      foo = matrix_get(mat, i, j);
      if (foo>max) {
        max = foo;
      }
    }
    matrix_set(out, i, 0, max);
  }
  return out;
}


matrix *matrix_inverse(matrix *A)
{
  int m, n, i, j, k;
  m = matrix_size_m(A);
  n = matrix_size_n(A);
  double c;

  // Kontrola ctvercove matice
  if (m != n) {
    fprintf(stderr, "%s: %i: ERROR: Inverzni matici lze vypocitat jenom ze ctvercove matice. Vstupni matice ma pritom rozmery (%i x %i).\n", __FILE__, __LINE__, m, n);
    exit(1);
  }

  // Vytvoreni matice [A I]
  matrix *AI = matrix_alloc(m, 2*m);
  for (i = 0; i<m; i++) {

    for (j = 0; j<m; j++) {
      matrix_set(AI, i, j, matrix_get(A, i, j));
    }

    for (j = m; j<2*m; j++) {
      if (j-m == i) {
        matrix_set(AI, i, j, 1);
      } else {
        matrix_set(AI, i, j, 0);
      }
    }

  }

  // Gauss-jordanova eliminacni metoda
  for (i = 0; i<m-1; i++) { // I - radek ktery bude prictem k J-tem radku

    for (j = 1+i; j<m; j++) { // J - upravovany radek
      c = matrix_get(AI, j, i)/matrix_get(AI, i, i);

      for (k = i; k<m; k++) { // K - upravovany sloupec
        matrix_set(AI, j, k, -matrix_get(AI, j, k)+c*matrix_get(AI, i, k));
      }

      for (k = m; k<2*m; k++) { // K - upravovany sloupec
        matrix_set(AI, j, k, -matrix_get(AI, j, k)+c*matrix_get(AI, i, k));
      }

    }

    c = matrix_get(AI, i, i);
    for (k = i; k<2*m; k++) { // K - upravovany sloupec
      matrix_set(AI, i, k, matrix_get(AI, i, k)/c);
    }

  }

  i = m-1;
  c = matrix_get(AI, i, i);
  for (k = i; k<2*m; k++) { // K - upravovany sloupec
    matrix_set(AI, i, k, matrix_get(AI, i, k)/c);
  }

  for (i = m-1; i>0; i--) { // I - radek ktery bude prictem k J-tem radku

    for (j = i-1; j>=0; j--) { // J - upravovany radek
      c = matrix_get(AI, j, i)/matrix_get(AI, i, i);

      matrix_set(AI, j, i, matrix_get(AI, j, i)-c*matrix_get(AI, i, i));

      for (k = m; k<2*m; k++) { // K - upravovany sloupec
        matrix_set(AI, j, k, matrix_get(AI, j, k)-c*matrix_get(AI, i, k));
      }

    }

  }

  

  matrix * iA = matrix_alloc(m, m);

  // Ulozeni inverzni matice
  for (i = 0; i<m; i++) {

    for (j = 0; j<m; j++) {

      matrix_set(iA, i, j, matrix_get(AI, i, j+m));

    }

  }

  matrix_free(AI);
  return iA;
}

matrix *matrix_levinson_durbin(matrix *R /* Musi byt symetricka matice (Töplitzova) */, unsigned int P)
{
  double km, Em, sum, ajm;
  unsigned int i, j;
  matrix *a;
  
  a = matrix_alloc(matrix_size_m(R), matrix_size_m(R));
  matrix_zeros(a);
  
  Em = matrix_get(R, 0, 0); // E(0) = R(0)
  km = -matrix_get(R, 1, 0)/matrix_get(R, 0, 0); // k(1) = -R(1)/R(0)
  matrix_set(a, 1, 1, km); // a(1) = k(1)
  Em = Em*(1-pow(km, 2)); // E(1) = E(0)*(1-k(1)^2)

  for (i = 2; i<=P-1; i++) {

    sum = 0;
    for (j = 1; j<=i-1; j++) {
      sum += matrix_get(a, j, i-1)*matrix_get(R, i-j, 0);
    }
    km = -(matrix_get(R, i, 0)+sum)/Em; // km = - ( R(m) + sum( a(j, m-1) * R(m-j) ) ) / E(m-1)
    matrix_set(a, i, i, km); // a(m, m) = k(m)

    for (j = 0; j<=i-1; j++) { // j = 1, 2, ... m-1
      ajm = matrix_get(a, j, i-1)+km*matrix_get(a, i-j, i-1); // a(j, m) = a(j, m-1) + k(m) * a(m-j, m-1)
      matrix_set(a, j, i, ajm);
    }

    Em = Em * (1-pow(km, 2)); // E(m) = E(m-1) * (1 - k(m)^2 )

  }

  return a;
}

matrix *
matrix_cat (matrix *A, matrix *B, unsigned int how) {
  unsigned int ma, na, mb, nb, m_new, n_new, i, j;
  matrix *NEW;
  // Rozmery
  ma = matrix_size_m(A);
  na = matrix_size_n(A);
  mb = matrix_size_m(B);
  nb = matrix_size_n(B);

  if (how == MATRIX_CAT_HORIZONTALLY) {
#ifdef MATRIX_DEBUG
    if (ma == mb) {
#endif
      // Nove rozmery matice
      m_new = ma;
      n_new = na + nb;
      // Alokace nove matice.
      NEW = matrix_alloc(m_new, n_new);

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

        // Naplneni nove matice daty ze struktury A
        for (j = 0; j < na; j++) {
          matrix_set(NEW, i, j, matrix_get(A, i, j));
        }

        // Naplneni nove matice daty ze struktury B
        for (j = na; j < n_new; j++) {
          matrix_set(NEW, i, j, matrix_get(B, i, j - na));
        }
      }

      return NEW;

#ifdef MATRIX_DEBUG
    } else {
      printf("%s:%i: error: Matice nelze hotizontalne spojit protoze ma kazda jiny pocet radku (MA = %i a MB = %i).\n", __FILE__, __LINE__, ma, mb);
      return NULL;
    }
#endif
  }

  if (how == MATRIX_CAT_VERTICALLY) {
#ifdef MATRIX_DEBUG
    if (na == nb) {
#endif
      // Nove rozmery matice
      m_new = ma + mb;
      n_new = na;
      // Alokace nove matice.
      NEW = matrix_alloc(m_new, n_new);

      for (j = 0; j < n_new; j++) {
        // Naplneni nove matice daty ze struktury A
        for (i = 0; i < ma; i++) {
          matrix_set(NEW, i, j, matrix_get(A, i, j));
        }

        // Naplneni nove matice daty ze struktury B
        for (i = ma; i < m_new; i++) {
          matrix_set(NEW, i, j, matrix_get(B, i - ma, j));
        }
      }

      return NEW;

#ifdef MATRIX_DEBUG
    } else {
      printf("%s:%i: error: Matice nelze vertikalne spojit protoze ma kazda jiny pocet sloupu (NA = %i a NB = %i).\n", __FILE__, __LINE__, na, nb);
      return NULL;
    }
#endif
  }
#ifdef MATRIX_DEBUG
  printf("%s:%i: error: Neznama hodnota '%i' vstupni promene HOW\n", __FILE__, __LINE__, how);
#endif
  return NULL;
}

matrix *
matrix_cat_list (list *L, unsigned int how) {
  unsigned int I, J, i, j, ma, na, mb, nb, M, N;
  matrix *A, *B, *OUT;

  if (list_get_size(L) == 0) {
#ifdef MATRIX_DEBUG
    printf("%s:%i: Seznam je prazdny.\n", __FILE__, __LINE__);
#endif
    return NULL;
  }

  if (how == MATRIX_CAT_HORIZONTALLY) {

    N = 0; // inicializace noveho poctu sloupcu

    // Zjisteni rozmeru prvniho prvku
    list_first(L);
    A = element_get_matrix(list_get_curent(L));
    ma = matrix_size_m(A);
    N = matrix_size_n(A);

    // Test jestli matice v seznamu maji stejny pocet radku.
    while (list_next(L)) {
      // Zjisteni rozmeru aktualniho prvku.
      B = element_get_matrix(list_get_curent(L));
      mb = matrix_size_m(B);
      nb = matrix_size_n(B);
      if (ma != mb) {
#ifdef MATRIX_DEBUG
        printf("%s:%i: Matice nelze spojit horizontalne, protoze nemaji stejny pocet radku.\n", __FILE__, __LINE__);
#endif
        return NULL;
      }

      // Ulozeni rozmeru tohoto prvku do dalsiho cyklu
      ma = mb;
      na = nb;
      // Inkrementace poctu sloupcu.
      N += nb;

    }

    M = ma;
    //printf("%i x %i\n", M, N);
    // Alokace vystupni matice.
    OUT = matrix_alloc(M, N);

    // naplneni struktury OUT daty.
    J = 0;
    N = 0;
    list_first(L);
    do {
      A = element_get_matrix(list_get_curent(L));
      N += matrix_size_n(A);
      j = 0;
      // bude se nastavovat J-ty sloupec
      while (J < N) {

        // Bude se nastavovat I-ty radek V j-tem sloupci.
        for (i = 0; i < M; i++) {
          matrix_set(OUT, i, J, matrix_get(A, i, j));
        }
        j++;
        J++;

      }
    } while (list_next(L));

    return OUT;

  }

  if (how == MATRIX_CAT_VERTICALLY) {

    M = 0; // inicializace noveho poctu radku

    // Zjisteni rozmeru prvniho prvku
    list_first(L);
    A = element_get_matrix(list_get_curent(L));
    M = matrix_size_m(A);
    na = matrix_size_n(A);

    // Test jestli struktury v seznamu maji stejny pocet radku.
    while (list_next(L)) {
      // Zjisteni rozmeru aktualniho prvku.
      B = element_get_matrix(list_get_curent(L));
      mb = matrix_size_m(B);
      nb = matrix_size_n(B);
      if (na != nb) {
#ifdef MATRIX_DEBUG
        printf("%s:%i: Matice nelze spojit horizontalne, protoze nemaji stejny pocet radku.\n", __FILE__, __LINE__);
#endif
        return NULL;
      }

      // Ulozeni rozmeru tohoto prvku do dalsiho cyklu
      ma = mb;
      na = nb;
      // Inkrementace poctu radku.
      M += mb;

    }

    N = na;

    // Alokace vystupni matice.
    OUT = matrix_alloc(M, N);

    // naplneni struktury OUT daty.
    I = 0;
    M = 0;
    list_first(L);
    do {
      A = element_get_matrix(list_get_curent(L));
      M += matrix_size_m(A);
      i = 0;
      // bude se nastavovat I-ty radek
      while (I < M) {

        // Bude se nastavovat J-ty sloupec V j-tem sloupci.
        for (j = 0; j < N; j++) {
          matrix_set(OUT, I, j, matrix_get(A, i, j));
        }
        i++;
        I++;
      }
    } while (list_next(L));

    //cell_print_info(OUT);

    return OUT;

  }

}

matrix *matrix_boolean(matrix *mat)
{
  matrix *mbool = matrix_alloc(matrix_size_m(mat),matrix_size_n(mat));
  unsigned int i, j;
  for (i = 0; i<matrix_size_m(mbool); i++) {
    for (j = 0; j<matrix_size_n(mbool); j++) {
      if (matrix_get(mat, i, j)==0) {
        matrix_set(mbool, i, j, 0);
      } else {
        matrix_set(mbool, i, j, 1);
      }
    }
  }
  return mbool;
}

matrix *matrix_neg(matrix *mat)
{
  matrix *mbool = matrix_alloc(matrix_size_m(mat),matrix_size_n(mat));
  unsigned int i, j;
  for (i = 0; i<matrix_size_m(mbool); i++) {
    for (j = 0; j<matrix_size_n(mbool); j++) {
      if (matrix_get(mat, i, j)==0) {
        matrix_set(mbool, i, j, 1);
      } else {
        matrix_set(mbool, i, j, 0);
      }
    }
  }
  return mbool;
}

matrix *matrix_and(matrix *A, matrix *B)
{
  if (matrix_size_m(A)!=matrix_size_m(B) |  matrix_size_n(A)!=matrix_size_n(B)) {
    printf("%s: %i: error: Matice A ma jiny rozmer nez matice B.\n", __FILE__, __LINE__, matrix_size_m(A), matrix_size_m(A), matrix_size_m(B), matrix_size_m(B));
    exit(1);
  }

  matrix *mbool = matrix_alloc(matrix_size_m(A),matrix_size_n(A));
  double aij, bij;
  unsigned int i, j;
  for (i = 0; i<matrix_size_m(A); i++) {
    for (j = 0; j<matrix_size_n(A); j++) {
      aij = matrix_get(A, i, j);
      bij = matrix_get(B, i, j);
      if (aij==bij & aij==1) {
        matrix_set(mbool, i, j, 1);
      } else {
        matrix_set(mbool, i, j, 0);
      }
    }
  }
  return mbool;
}

matrix *matrix_or(matrix *A, matrix *B)
{
  if (matrix_size_m(A)!=matrix_size_m(B) |  matrix_size_n(A)!=matrix_size_n(B)) {
    printf("%s: %i: error: Matice A () ma jiny rozmer nez matice B ().\n", __FILE__, __LINE__, matrix_size_m(A), matrix_size_m(A), matrix_size_m(B), matrix_size_m(B));
    exit(1);
  }

  matrix *mbool = matrix_alloc(matrix_size_m(A),matrix_size_n(A));
  double aij, bij;
  unsigned int i, j;
  for (i = 0; i<matrix_size_m(A); i++) {
    for (j = 0; j<matrix_size_n(A); j++) {
      aij = matrix_get(A, i, j);
      bij = matrix_get(B, i, j);
      if (aij==1 | bij==1) {
        matrix_set(mbool, i, j, 1);
      } else {
        matrix_set(mbool, i, j, 0);
      }
    }
  }
  return mbool;
}

matrix *matrix_biimp(matrix *A, matrix *B)
{
  if (matrix_size_m(A)!=matrix_size_m(B) |  matrix_size_n(A)!=matrix_size_n(B)) {
    printf("%s: %i: error: Matice A () ma jiny rozmer nez matice B ().\n", __FILE__, __LINE__, matrix_size_m(A), matrix_size_m(A), matrix_size_m(B), matrix_size_m(B));
    exit(1);
  }

  matrix *mbool = matrix_alloc(matrix_size_m(A),matrix_size_n(A));
  double aij, bij;
  unsigned int i, j;
  for (i = 0; i<matrix_size_m(A); i++) {
    for (j = 0; j<matrix_size_n(A); j++) {
      aij = matrix_get(A, i, j);
      bij = matrix_get(B, i, j);
      if (aij==bij) {
        matrix_set(mbool, i, j, 1);
      } else {
        matrix_set(mbool, i, j, 0);
      }
    }
  }
  return mbool;
}


