text:
/* test.c
 *
 * Copyright (C) 2012-2014, 2016 Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <unistd.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_spmatrix.h>

/*
create_random_sparse()
  Create a random sparse matrix with approximately
M*N*density non-zero entries

Inputs: M       - number of rows
        N       - number of columns
        density - sparse density \in [0,1]
                  0 = no non-zero entries
                  1 = all m*n entries are filled
        r       - random number generator

Return: pointer to sparse matrix in triplet format (must be freed by caller)

Notes:
1) non-zero matrix entries are uniformly distributed in [0,1]
*/

static gsl_spmatrix *
create_random_sparse(const size_t M, const size_t N, const double density,
                     const gsl_rng *r)
{
  size_t nnzwanted = (size_t) floor(M * N * GSL_MIN(density, 1.0));
  gsl_spmatrix *m = gsl_spmatrix_alloc_nzmax(M, N,
                                             nnzwanted,
                                             GSL_SPMATRIX_TRIPLET);

  while (gsl_spmatrix_nnz(m) < nnzwanted)
    {
      /* generate a random row and column */
      size_t i = gsl_rng_uniform(r) * M;
      size_t j = gsl_rng_uniform(r) * N;

      /* generate random m_{ij} and add it */
      double x = gsl_rng_uniform(r);
      gsl_spmatrix_set(m, i, j, x);
    }

  return m;
} /* create_random_sparse() */

static gsl_spmatrix *
create_random_sparse_int(const size_t M, const size_t N, const double density,
                         const gsl_rng *r)
{
  const double lower = 1.0;
  const double upper = 10.0;
  size_t nnzwanted = (size_t) floor(M * N * GSL_MIN(density, 1.0));
  gsl_spmatrix *m = gsl_spmatrix_alloc_nzmax(M, N,
                                             nnzwanted,
                                             GSL_SPMATRIX_TRIPLET);

  while (gsl_spmatrix_nnz(m) < nnzwanted)
    {
      /* generate a random row and column */
      size_t i = gsl_rng_uniform(r) * M;
      size_t j = gsl_rng_uniform(r) * N;

      /* generate random m_{ij} and add it */
      int x = (int) (gsl_rng_uniform(r) * (upper - lower) + lower);
      gsl_spmatrix_set(m, i, j, (double) x);
    }

  return m;
}

static void
test_getset(const size_t M, const size_t N,
            const double density, const gsl_rng *r)
{
  int status;
  size_t i, j;

  /* test that a newly allocated matrix is initialized to zero */
  {
    size_t nzmax = (size_t) (0.5 * M * N);
    gsl_spmatrix *mt = gsl_spmatrix_alloc_nzmax(M, N, nzmax, GSL_SPMATRIX_TRIPLET);
    gsl_spmatrix *mccs = gsl_spmatrix_alloc_nzmax(M, N, nzmax, GSL_SPMATRIX_CCS);
    gsl_spmatrix *mcrs = gsl_spmatrix_alloc_nzmax(M, N, nzmax, GSL_SPMATRIX_CRS);

    status = 0;
    for (i = 0; i < M; ++i)
      {
        for (j = 0; j < N; ++j)
          {
            double mtij = gsl_spmatrix_get(mt, i, j);
            double mcij = gsl_spmatrix_get(mccs, i, j);
            double mrij = gsl_spmatrix_get(mcrs, i, j);

            if (mtij != 0.0)
              status = 1;

            if (mcij != 0.0)
              status = 2;

            if (mrij != 0.0)
              status = 3;
          }
      }

    gsl_test(status == 1, "test_getset: triplet M=%zu N=%zu not initialized to zero", M, N);
    gsl_test(status == 2, "test_getset: CCS M=%zu N=%zu not initialized to zero", M, N);
    gsl_test(status == 3, "test_getset: CRS M=%zu N=%zu not initialized to zero", M, N);

    gsl_spmatrix_free(mt);
    gsl_spmatrix_free(mccs);
    gsl_spmatrix_free(mcrs);
  }

  /* test triplet versions of _get and _set */
  {
    const double val = 0.75;
    size_t k = 0;
    gsl_spmatrix *m = gsl_spmatrix_alloc(M, N);

    status = 0;
    for (i = 0; i < M; ++i)
      {
        for (j = 0; j < N; ++j)
          {
            double x = (double) ++k;
            double y;

            gsl_spmatrix_set(m, i, j, x);
            y = gsl_spmatrix_get(m, i, j);
            if (x != y)
              status = 1;
          }
      }

    gsl_test(status, "test_getset: M=%zu N=%zu _get != _set", M, N);

    /* test setting an element to 0 */
    gsl_spmatrix_set(m, 0, 0, 1.0);
    gsl_spmatrix_set(m, 0, 0, 0.0);

    status = gsl_spmatrix_get(m, 0, 0) != 0.0;
    gsl_test(status, "test_getset: M=%zu N=%zu m(0,0) = %f",
             M, N, gsl_spmatrix_get(m, 0, 0));

    /* test gsl_spmatrix_set_zero() */
    gsl_spmatrix_set(m, 0, 0, 1.0);
    gsl_spmatrix_set_zero(m);
    status = gsl_spmatrix_get(m, 0, 0) != 0.0;
    gsl_test(status, "test_getset: M=%zu N=%zu set_zero m(0,0) = %f",
             M, N, gsl_spmatrix_get(m, 0, 0));

    /* resassemble matrix to ensure nz is calculated correctly */
    k = 0;
    for (i = 0; i < M; ++i)
      {
        for (j = 0; j < N; ++j)
          {
            double x = (double) ++k;
            gsl_spmatrix_set(m, i, j, x);
          }
      }

    status = gsl_spmatrix_nnz(m) != M * N;
    gsl_test(status, "test_getset: M=%zu N=%zu set_zero nz = %zu",
             M, N, gsl_spmatrix_nnz(m));

    /* test gsl_spmatrix_ptr() */
    status = 0;
    for (i = 0; i < M; ++i)
      {
        for (j = 0; j < N; ++j)
          {
            double mij = gsl_spmatrix_get(m, i, j);
            double *ptr = gsl_spmatrix_ptr(m, i, j);

            *ptr += val;
            if (gsl_spmatrix_get(m, i, j) != mij + val)
              status = 2;
          }
      }

    gsl_test(status == 2, "test_getset: M=%zu N=%zu triplet ptr", M, N);

    gsl_spmatrix_free(m);
  }

  /* test duplicate values are handled correctly */
  {
    size_t min = GSL_MIN(M, N);
    size_t expected_nnz = min;
    size_t nnz;
    size_t k = 0;
    gsl_spmatrix *m = gsl_spmatrix_alloc(M, N);

    status = 0;
    for (i = 0; i < min; ++i)
      {
        for (j = 0; j < 5; ++j)
          {
            double x = (double) ++k;
            double y;

            gsl_spmatrix_set(m, i, i, x);
            y = gsl_spmatrix_get(m, i, i);
            if (x != y)
              status = 1;
          }
      }

    gsl_test(status, "test_getset: duplicate test M=%zu N=%zu _get != _set", M, N);

    nnz = gsl_spmatrix_nnz(m);
    status = nnz != expected_nnz;
    gsl_test(status, "test_getset: duplicate test M=%zu N=%zu nnz=%zu, expected=%zu",
             M, N, nnz, expected_nnz);

    gsl_spmatrix_free(m);
  }

  /* test CCS version of gsl_spmatrix_get() */
  {
    const double val = 0.75;
    gsl_spmatrix *T = create_random_sparse(M, N, density, r);
    gsl_spmatrix *C = gsl_spmatrix_ccs(T);

    status = 0;
    for (i = 0; i < M; ++i)
      {
        for (j = 0; j < N; ++j)
          {
            double Tij = gsl_spmatrix_get(T, i, j);
            double Cij = gsl_spmatrix_get(C, i, j);
            double *ptr = gsl_spmatrix_ptr(C, i, j);

            if (Tij != Cij)
              status = 1;

            if (ptr)
              {
                *ptr += val;
                Cij = gsl_spmatrix_get(C, i, j);
                if (Tij + val != Cij)
                  status = 2;
              }
          }
      }

    gsl_test(status == 1, "test_getset: M=%zu N=%zu CCS get", M, N);
    gsl_test(status == 2, "test_getset: M=%zu N=%zu CCS ptr", M, N);

    gsl_spmatrix_free(T);
    gsl_spmatrix_free(C);
  }

  /* test CRS version of gsl_spmatrix_get() */
  {
    const double val = 0.75;
    gsl_spmatrix *T = create_random_sparse(M, N, density, r);
    gsl_spmatrix *C = gsl_spmatrix_crs(T);

    status = 0;
    for (i = 0; i < M; ++i)
      {
        for (j = 0; j < N; ++j)
          {
            double Tij = gsl_spmatrix_get(T, i, j);
            double Cij = gsl_spmatrix_get(C, i, j);
            double *ptr = gsl_spmatrix_ptr(C, i, j);

            if (Tij != Cij)
              status = 1;

            if (ptr)
              {
                *ptr += val;
                Cij = gsl_spmatrix_get(C, i, j);
                if (Tij + val != Cij)
                  status = 2;
              }
          }
      }

    gsl_test(status == 1, "test_getset: M=%zu N=%zu CRS get", M, N);
    gsl_test(status == 2, "test_getset: M=%zu N=%zu CRS ptr", M, N);

    gsl_spmatrix_free(T);
    gsl_spmatrix_free(C);
  }
} /* test_getset() */

static void
test_memcpy(const size_t M, const size_t N,
            const double density, const gsl_rng *r)
{
  int status;

  {
    gsl_spmatrix *A = create_random_sparse(M, N, density, r);
    gsl_spmatrix *A_ccs = gsl_spmatrix_ccs(A);
    gsl_spmatrix *A_crs = gsl_spmatrix_crs(A);
    gsl_spmatrix *B_t, *B_ccs, *B_crs;

    B_t = gsl_spmatrix_alloc(M, N);
    gsl_spmatrix_memcpy(B_t, A);

    status = gsl_spmatrix_equal(A, B_t) != 1;
    gsl_test(status, "test_memcpy: _memcpy M=%zu N=%zu triplet format", M, N);

    B_ccs = gsl_spmatrix_alloc_nzmax(M, N, A_ccs->nzmax, GSL_SPMATRIX_CCS);
    B_crs = gsl_spmatrix_alloc_nzmax(M, N, A_ccs->nzmax, GSL_SPMATRIX_CRS);

    gsl_spmatrix_memcpy(B_ccs, A_ccs);
    gsl_spmatrix_memcpy(B_crs, A_crs);

    status = gsl_spmatrix_equal(A_ccs, B_ccs) != 1;
    gsl_test(status, "test_memcpy: _memcpy M=%zu N=%zu CCS", M, N);

    status = gsl_spmatrix_equal(A_crs, B_crs) != 1;
    gsl_test(status, "test_memcpy: _memcpy M=%zu N=%zu CRS", M, N);

    gsl_spmatrix_free(A);
    gsl_spmatrix_free(A_ccs);
    gsl_spmatrix_free(A_crs);
    gsl_spmatrix_free(B_t);
    gsl_spmatrix_free(B_ccs);
    gsl_spmatrix_free(B_crs);
  }

  /* test transpose_memcpy */
  {
    gsl_spmatrix *A = create_random_sparse(M, N, density, r);
    gsl_spmatrix *AT = gsl_spmatrix_alloc(N, M);
    gsl_spmatrix *B = gsl_spmatrix_ccs(A);
    gsl_spmatrix *BT = gsl_spmatrix_alloc_nzmax(N, M, 1, GSL_SPMATRIX_CCS);
    gsl_spmatrix *C = gsl_spmatrix_crs(A);
    gsl_spmatrix *CT = gsl_spmatrix_alloc_nzmax(N, M, 1, GSL_SPMATRIX_CRS);
    size_t i, j;

    gsl_spmatrix_transpose_memcpy(AT, A);
    gsl_spmatrix_transpose_memcpy(BT, B);
    gsl_spmatrix_transpose_memcpy(CT, C);

    status = 0;
    for (i = 0; i < M; ++i)
      {
        for (j = 0; j < N; ++j)
          {
            double Aij = gsl_spmatrix_get(A, i, j);
            double ATji = gsl_spmatrix_get(AT, j, i);

            if (Aij != ATji)
              status = 1;
          }
      }

    gsl_test(status, "test_memcpy: _transpose_memcpy M=%zu N=%zu triplet format", M, N);

    status = 0;
    for (i = 0; i < M; ++i)
      {
        for (j = 0; j < N; ++j)
          {
            double Aij = gsl_spmatrix_get(A, i, j);
            double Bij = gsl_spmatrix_get(B, i, j);
            double BTji = gsl_spmatrix_get(BT, j, i);

            if ((Bij != BTji) || (Aij != Bij))
              status = 1;
          }
      }

    gsl_test(status, "test_memcpy: _transpose_memcpy M=%zu N=%zu CCS format", M, N);

    status = 0;
    for (i = 0; i < M; ++i)
      {
        for (j = 0; j < N; ++j)
          {
            double Aij = gsl_spmatrix_get(A, i, j);
            double Cij = gsl_spmatrix_get(C, i, j);
            double CTji = gsl_spmatrix_get(CT, j, i);

            if ((Cij != CTji) || (Aij != Cij))
              status = 1;
          }
      }

    gsl_test(status, "test_memcpy: _transpose_memcpy M=%zu N=%zu CRS format", M, N);

    gsl_spmatrix_free(A);
    gsl_spmatrix_free(AT);
    gsl_spmatrix_free(B);
    gsl_spmatrix_free(BT);
    gsl_spmatrix_free(C);
    gsl_spmatrix_free(CT);
  }
} /* test_memcpy() */

static void
test_transpose(const size_t M, const size_t N,
               const double density, const gsl_rng *r)
{
  int status;
  gsl_spmatrix *A = create_random_sparse(M, N, density, r);
  gsl_spmatrix *AT = gsl_spmatrix_alloc_nzmax(M, N, A->nz, A->sptype);
  gsl_spmatrix *AT2 = gsl_spmatrix_alloc_nzmax(M, N, A->nz, A->sptype);
  gsl_spmatrix *AT2_ccs, *AT2_crs;
  size_t i, j;

  /* test triplet transpose */

  gsl_spmatrix_memcpy(AT, A);
  gsl_spmatrix_memcpy(AT2, A);

  gsl_spmatrix_transpose(AT);
  gsl_spmatrix_transpose2(AT2);

  status = 0;
  for (i = 0; i < M; ++i)
    {
      for (j = 0; j < N; ++j)
        {
          double Aij = gsl_spmatrix_get(A, i, j);
          double ATji = gsl_spmatrix_get(AT, j, i);
          double AT2ji = gsl_spmatrix_get(AT2, j, i);

          if (Aij != ATji)
            status = 1;

          if (Aij != AT2ji)
            status = 2;
        }
    }

  gsl_test(status == 1, "test_transpose: transpose M=%zu N=%zu triplet format",
           M, N);
  gsl_test(status == 2, "test_transpose: transpose2 M=%zu N=%zu triplet format",
           M, N);

  /* test CCS transpose */

  AT2_ccs = gsl_spmatrix_ccs(A);
  gsl_spmatrix_transpose2(AT2_ccs);

  status = 0;
  for (i = 0; i < M; ++i)
    {
      for (j = 0; j < N; ++j)
        {
          double Aij = gsl_spmatrix_get(A, i, j);
          double AT2ji = gsl_spmatrix_get(AT2_ccs, j, i);

          if (Aij != AT2ji)
            status = 2;
        }
    }

  gsl_test(status == 2, "test_transpose: transpose2 M=%zu N=%zu CCS format",
           M, N);

  /* test CRS transpose */

  AT2_crs = gsl_spmatrix_crs(A);
  gsl_spmatrix_transpose2(AT2_crs);

  status = 0;
  for (i = 0; i < M; ++i)
    {
      for (j = 0; j < N; ++j)
        {
          double Aij = gsl_spmatrix_get(A, i, j);
          double AT2ji = gsl_spmatrix_get(AT2_crs, j, i);

          if (Aij != AT2ji)
            status = 2;
        }
    }

  gsl_test(status == 2, "test_transpose: transpose2 M=%zu N=%zu CRS format",
           M, N);

  gsl_spmatrix_free(A);
  gsl_spmatrix_free(AT);
  gsl_spmatrix_free(AT2);
  gsl_spmatrix_free(AT2_ccs);
  gsl_spmatrix_free(AT2_crs);
}

static void
test_ops(const size_t M, const size_t N,
         const double density, const gsl_rng *r)
{
  size_t i, j;
  int status;

  /* test gsl_spmatrix_add */
  {
    gsl_spmatrix *A = create_random_sparse(M, N, density, r);
    gsl_spmatrix *B = create_random_sparse(M, N, density, r);

    gsl_spmatrix *A_ccs = gsl_spmatrix_ccs(A);
    gsl_spmatrix *B_ccs = gsl_spmatrix_ccs(B);
    gsl_spmatrix *C_ccs = gsl_spmatrix_alloc_nzmax(M, N, 1, GSL_SPMATRIX_CCS);

    gsl_spmatrix *A_crs = gsl_spmatrix_crs(A);
    gsl_spmatrix *B_crs = gsl_spmatrix_crs(B);
    gsl_spmatrix *C_crs = gsl_spmatrix_alloc_nzmax(M, N, 1, GSL_SPMATRIX_CRS);

    gsl_spmatrix_add(C_ccs, A_ccs, B_ccs);
    gsl_spmatrix_add(C_crs, A_crs, B_crs);

    status = 0;
    for (i = 0; i < M; ++i)
      {
        for (j = 0; j < N; ++j)
          {
            double aij, bij, cij;

            aij = gsl_spmatrix_get(A_ccs, i, j);
            bij = gsl_spmatrix_get(B_ccs, i, j);
            cij = gsl_spmatrix_get(C_ccs, i, j);
            if (aij + bij != cij)
              status = 1;

            aij = gsl_spmatrix_get(A_crs, i, j);
            bij = gsl_spmatrix_get(B_crs, i, j);
            cij = gsl_spmatrix_get(C_crs, i, j);
            if (aij + bij != cij)
              status = 2;
          }
      }

    gsl_test(status == 1, "test_ops: add M=%zu N=%zu CCS", M, N);
    gsl_test(status == 2, "test_ops: add M=%zu N=%zu CRS", M, N);

    gsl_spmatrix_free(A);
    gsl_spmatrix_free(B);
    gsl_spmatrix_free(A_ccs);
    gsl_spmatrix_free(B_ccs);
    gsl_spmatrix_free(C_ccs);
    gsl_spmatrix_free(A_crs);
    gsl_spmatrix_free(B_crs);
    gsl_spmatrix_free(C_crs);
  }
} /* test_ops() */

static void
test_io_ascii(const size_t M, const size_t N,
              const double density, const gsl_rng *r)
{
  int status;
  gsl_spmatrix *A = create_random_sparse_int(M, N, density, r);

  char filename[] = "test.dat";

  /* test triplet I/O */
  {
    FILE *f = fopen(filename, "w");

    gsl_spmatrix_fprintf(f, A, "%g");

    fclose(f);
  }

  {
    FILE *f = fopen(filename, "r");
    gsl_spmatrix *B = gsl_spmatrix_fscanf(f);

    status = gsl_spmatrix_equal(A, B) != 1;
    gsl_test(status, "test_io_ascii: fprintf/fscanf M=%zu N=%zu triplet format", M, N);

    fclose(f);
    gsl_spmatrix_free(B);
  }

  /* test CCS I/O */
  {
    FILE *f = fopen(filename, "w");
    gsl_spmatrix *A_ccs = gsl_spmatrix_ccs(A);

    gsl_spmatrix_fprintf(f, A_ccs, "%g");

    fclose(f);
    gsl_spmatrix_free(A_ccs);
  }

  {
    FILE *f = fopen(filename, "r");
    gsl_spmatrix *B = gsl_spmatrix_fscanf(f);

    status = gsl_spmatrix_equal(A, B) != 1;
    gsl_test(status, "test_io_ascii: fprintf/fscanf M=%zu N=%zu CCS format", M, N);

    fclose(f);
    gsl_spmatrix_free(B);
  }

  /* test CRS I/O */
  {
    FILE *f = fopen(filename, "w");
    gsl_spmatrix *A_crs = gsl_spmatrix_crs(A);

    gsl_spmatrix_fprintf(f, A_crs, "%g");

    fclose(f);
    gsl_spmatrix_free(A_crs);
  }

  {
    FILE *f = fopen(filename, "r");
    gsl_spmatrix *B = gsl_spmatrix_fscanf(f);

    status = gsl_spmatrix_equal(A, B) != 1;
    gsl_test(status, "test_io_ascii: fprintf/fscanf M=%zu N=%zu CRS format", M, N);

    fclose(f);
    gsl_spmatrix_free(B);
  }

  unlink(filename);

  gsl_spmatrix_free(A);
}

static void
test_io_binary(const size_t M, const size_t N,
               const double density, const gsl_rng *r)
{
  int status;
  gsl_spmatrix *A = create_random_sparse(M, N, density, r);
  gsl_spmatrix *A_ccs, *A_crs;

  char filename[] = "test.dat";

  /* test triplet I/O */
  {
    FILE *f = fopen(filename, "wb");

    gsl_spmatrix_fwrite(f, A);

    fclose(f);
  }

  {
    FILE *f = fopen(filename, "rb");
    gsl_spmatrix *B = gsl_spmatrix_alloc_nzmax(M, N, A->nz, A->sptype);

    gsl_spmatrix_fread(f, B);

    status = gsl_spmatrix_equal(A, B) != 1;
    gsl_test(status, "test_io_binary: fwrite/fread M=%zu N=%zu triplet format", M, N);

    fclose(f);
    gsl_spmatrix_free(B);
  }

  /* test CCS I/O */
  A_ccs = gsl_spmatrix_ccs(A);

  {
    FILE *f = fopen(filename, "wb");

    gsl_spmatrix_fwrite(f, A_ccs);

    fclose(f);
  }

  {
    FILE *f = fopen(filename, "rb");
    gsl_spmatrix *B = gsl_spmatrix_alloc_nzmax(M, N, A->nz, GSL_SPMATRIX_CCS);

    gsl_spmatrix_fread(f, B);

    status = gsl_spmatrix_equal(A_ccs, B) != 1;
    gsl_test(status, "test_io_binary: fwrite/fread M=%zu N=%zu CCS format", M, N);

    fclose(f);
    gsl_spmatrix_free(B);
  }

  /* test CRS I/O */
  A_crs = gsl_spmatrix_crs(A);

  {
    FILE *f = fopen(filename, "wb");

    gsl_spmatrix_fwrite(f, A_crs);

    fclose(f);
  }

  {
    FILE *f = fopen(filename, "rb");
    gsl_spmatrix *B = gsl_spmatrix_alloc_nzmax(M, N, A->nz, GSL_SPMATRIX_CRS);

    gsl_spmatrix_fread(f, B);

    status = gsl_spmatrix_equal(A_crs, B) != 1;
    gsl_test(status, "test_io_binary: fwrite/fread M=%zu N=%zu CRS format", M, N);

    fclose(f);
    gsl_spmatrix_free(B);
  }

  unlink(filename);

  gsl_spmatrix_free(A);
  gsl_spmatrix_free(A_ccs);
  gsl_spmatrix_free(A_crs);
}

int
main()
{
  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);

  test_memcpy(10, 10, 0.2, r);
  test_memcpy(10, 15, 0.3, r);
  test_memcpy(53, 213, 0.4, r);
  test_memcpy(920, 2, 0.2, r);
  test_memcpy(2, 920, 0.3, r);

  test_getset(20, 20, 0.3, r);
  test_getset(30, 20, 0.3, r);
  test_getset(15, 210, 0.3, r);

  test_transpose(50, 50, 0.5, r);
  test_transpose(10, 40, 0.3, r);
  test_transpose(40, 10, 0.3, r);
  test_transpose(57, 13, 0.2, r);

  test_ops(20, 20, 0.2, r);
  test_ops(50, 20, 0.3, r);
  test_ops(20, 50, 0.3, r);
  test_ops(76, 43, 0.4, r);

  test_io_ascii(30, 30, 0.3, r);
  test_io_ascii(20, 10, 0.2, r);
  test_io_ascii(10, 20, 0.2, r);
  test_io_ascii(34, 78, 0.3, r);

  test_io_binary(50, 50, 0.3, r);
  test_io_binary(25, 10, 0.2, r);
  test_io_binary(10, 25, 0.2, r);
  test_io_binary(101, 253, 0.3, r);

  gsl_rng_free(r);

  exit (gsl_test_summary());
} /* main() */

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/*
 *
 *
 *
 */

#include <inttypes.h>
#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <cblas.h>
#include <string.h> // memcpy
#include <math.h>
#include <stdarg.h>

#include "cgraph.h"
#include "cg_operation.h"
#include "cg_types.h"
#include "cg_variables.h"
#include "cg_errors.h"
#include "cg_constants.h"
#include "cg_enums.h"
#include "cg_factory.h"
#include "cg_math.h"
#include "cg_ops.h"

void cg_assert(int cond, const char * rawcond, const char * fmt, ...)
{
	if(cond)
		return;

	char temp[1024];
	va_list vl;
	va_start(vl, fmt);
	vsprintf(temp, fmt, vl);
	va_end(vl);
	fprintf(stdout, "Fatal error, assertion failed: %s\n", rawcond);
	fprintf(stdout, temp);
	fprintf(stdout, "\n");

	exit(-1);
}

#define CHECK_RESULT(node) \
if(node->error != NULL){\
	return node;\
}

/*
 * Works only with constant types
 */
CGNode* copyNode(CGNode* node){
	CGNode* n = calloc(1, sizeof(CGNode));
	if(node->type != CGNT_CONSTANT){
		fprintf(stderr, "Call to copyNode with a non-constant node.\n... This should not happen, but who knows these days.\n");
		fprintf(stderr, "copy node with nodetype = %d\n", node->type);
		exit(-1);
	}

	n->type = CGNT_CONSTANT;
	n->constant = calloc(1, sizeof(CGPConstant));
	n->constant->type = node->constant->type;
	n->result = NULL;
	n->diff = NULL;

	switch(node->constant->type){
		case CGVT_DOUBLE: {
			CGDouble* d = calloc(1, sizeof(CGDouble));
			d->value = ((CGDouble*)node->constant->value)->value;

			n->constant->value = d;
			break;
		}

		case CGVT_VECTOR: {
			CGVector* V = calloc(1, sizeof(CGVector));
			CGVector* src = (CGVector*)node->constant->value;

			V->len = src->len;
			if(src->data != NULL) {
				V->data = calloc(V->len, sizeof(cg_float));
				memcpy(V->data, src->data, V->len * sizeof(cg_float));
			}

#ifdef CG_USE_OPENCL
	        V->buf = NULL;
	        V->loc = CG_DATALOC_HOST_MEM;
#endif
			n->constant->value = V;
			break;
		}

		case CGVT_MATRIX: {
			CGMatrix* M = calloc(1, sizeof(CGMatrix));
			CGMatrix* src = (CGMatrix*)node->constant->value;

			uint64_t size = src->cols * src->rows;

			M->rows = src->rows;
			M->cols = src->cols;

			if(src->data != NULL) {
				M->data = calloc(size, sizeof(cg_float));
				memcpy(M->data, src->data, size*sizeof(cg_float));
			}

#ifdef CG_USE_OPENCL
			M->buf = NULL;
			M->loc = CG_DATALOC_HOST_MEM;
#endif
			n->constant->value = M;
			break;
		}
	}

	return n;
}

// @Deprecated @NotUsed
void* copyNodeValue(CGNode* node){
	if(node->type != CGNT_CONSTANT){
		fprintf(stderr, "Call to copyNodeValue with a non-constant node.\n... This should not happen, but who knows these days.");
		exit(-1);
	}

	switch(node->constant->type){
		case CGVT_DOUBLE: {
			CGDouble* d = calloc(1, sizeof(CGDouble));
			d->value = ((CGDouble*)node->constant->value)->value;

			return d;
		}

		case CGVT_VECTOR: {
			CGVector* V = calloc(1, sizeof(CGVector));
			CGVector* src = (CGVector*)node->constant->value;

			V->len = src->len;
			V->data = calloc(V->len, sizeof(cg_float));

			memcpy(V->data, src->data, V->len*sizeof(cg_float));
			return V;
		}

		case CGVT_MATRIX: {
			CGMatrix* M = calloc(1, sizeof(CGMatrix));
			CGMatrix* src = (CGMatrix*)node->constant->value;

			uint64_t size = src->cols * src->rows;

			M->rows = src->rows;
			M->cols = src->cols;
			M->data = calloc(size, sizeof(cg_float));

			memcpy(M->data, src->data, size*sizeof(cg_float));
			return M;
		}
	}
}

/*
 * @Deprecated
 */
void* copyRNodeValue(CGResultNode* node){
	switch(node->type){
		case CGVT_DOUBLE: {
			CGDouble* d = calloc(1, sizeof(CGDouble));
			d->value = ((CGDouble*)node->value)->value;

			return d;
		}

		case CGVT_VECTOR: {
			CGVector* V = calloc(1, sizeof(CGVector));
			CGVector* src = (CGVector*)node->value;

			V->len = src->len;
			V->data = calloc(V->len, sizeof(cg_float));

			memcpy(V->data, src->data, V->len*sizeof(cg_float));
			return V;
		}

		case CGVT_MATRIX: {
			CGMatrix* M = calloc(1, sizeof(CGMatrix));
			CGMatrix* src = (CGMatrix*)node->value;

			uint64_t size = src->cols * src->rows;

			M->rows = src->rows;
			M->cols = src->cols;
			M->data = calloc(size, sizeof(cg_float));

			memcpy(M->data, src->data, size*sizeof(cg_float));
			return M;
		}
	}
}

CGResultNode* copyResultNode(CGResultNode* node){
	CGResultNode* res = calloc(1, sizeof(CGResultNode));

	switch(node->type){
		case CGVT_DOUBLE: {
			CGDouble* d = calloc(1, sizeof(CGDouble));
			CGDouble* d2 = (CGDouble*)node->value;
			d->value = d2->value;

			res->value = d;
			break;
		}

		case CGVT_VECTOR: {
			CGVector* V = calloc(1, sizeof(CGVector));
			CGVector* src = (CGVector*)node->value;

			V->len = src->len;
			V->data = calloc(V->len, sizeof(cg_float));

#ifdef CG_USE_OPENCL
			V->buf = src->buf;
#endif
			memcpy(V->data, src->data, V->len*sizeof(cg_float));

			res->value = V;
			break;
		}

		case CGVT_MATRIX: {
			CGMatrix* M = calloc(1, sizeof(CGMatrix));
			CGMatrix* src = (CGMatrix*)node->value;

			uint64_t size = src->cols * src->rows;

			M->rows = src->rows;
			M->cols = src->cols;
			M->data = calloc(size, sizeof(cg_float));
#ifdef CG_USE_OPENCL
            M->buf = src->buf;
#endif

			memcpy(M->data, src->data, size*sizeof(cg_float));

			res->value = M;
			break;
		}
	}
	res->type = node->type;

	return res;
}

CGMatrix* vectorToMatrix(CGVector* v){
	CGMatrix* m = calloc(1, sizeof(CGMatrix));
	m->rows = v->len;
	m->cols = 1;
	m->data = v->data;

	return m;
}

/*
 * Computational Graph traversing
 */

CGResultNode* processUnaryOperation(CGraph* graph, CGUnaryOperationType type, CGNode* uhs, CGNode* parentNode){
	CGVarType uhsType = CGVT_DOUBLE;
	void* uhsValue = NULL;
	CGResultNode* newres = NULL;

	CGResultNode* lhsResult = computeCGNode(graph, uhs);
	CHECK_RESULT(lhsResult)
	uhsType = lhsResult->type;
	uhsValue = lhsResult->value;

	switch(type){
		case CGUOT_EXP:{
			if(uhsType == CGVT_DOUBLE){
				newres = expD((CGDouble*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if(uhsType == CGVT_VECTOR){
				newres = expV((CGVector*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if(uhsType == CGVT_MATRIX){
				newres = expM((CGMatrix*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}
			break;
		}

		case CGUOT_LOG:{
			if(uhsType == CGVT_DOUBLE){
				newres = logD((CGDouble*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if(uhsType == CGVT_VECTOR){
				newres = logV((CGVector*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if(uhsType == CGVT_MATRIX){
				newres = logM((CGMatrix*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}
			break;
		}

		case CGUOT_MINUS:{
			if(uhsType == CGVT_DOUBLE){
				CGDouble* rhs = calloc(1, sizeof(CGDouble));
				rhs->value = -1;

				CGResultNode* res = mulDD((CGDouble*)uhsValue, rhs, graph, parentNode);

				freeDoubleValue(rhs);

				parentNode->result = res;
				return res;
			}

			if(uhsType == CGVT_VECTOR){
				CGDouble* lhs = calloc(1, sizeof(CGDouble));
				lhs->value = -1;

				CGResultNode* res = mulDV(lhs, (CGVector*)uhsValue, graph, parentNode);

				free(lhs);
				parentNode->result = res;
				return res;
			}

			if(uhsType == CGVT_MATRIX){
				CGDouble* lhs = calloc(1, sizeof(CGDouble));
				lhs->value = -1;
				CGResultNode* res = mulDM(lhs, (CGMatrix*)uhsValue, graph, parentNode);

				free(lhs);

				parentNode->result = res;
				return res;
			}
			break;
		}

		case CGUOT_SIN:{
			if(uhsType == CGVT_DOUBLE){
				newres = sinD((CGDouble*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if(uhsType == CGVT_VECTOR){
				newres = sinV((CGVector*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if(uhsType == CGVT_MATRIX){
				newres = sinM((CGMatrix*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}
			break;
		}

		case CGUOT_COS:{
			if(uhsType == CGVT_DOUBLE){
				newres = cosD((CGDouble*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if(uhsType == CGVT_VECTOR){
				newres = cosV((CGVector*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if(uhsType == CGVT_MATRIX){
				newres = cosM((CGMatrix*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}
			break;
		}

		case CGUOT_TAN:{
			if(uhsType == CGVT_DOUBLE){
				newres = tanD((CGDouble*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if(uhsType == CGVT_VECTOR){
				newres = tanV((CGVector*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if(uhsType == CGVT_MATRIX){
				newres = tanM((CGMatrix*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}
			break;
		}

		case CGUOT_TANH:{
			if(uhsType == CGVT_DOUBLE){
				newres = tanhD((CGDouble*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if(uhsType == CGVT_VECTOR){
				newres = tanhV((CGVector*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if(uhsType == CGVT_MATRIX){
				newres = tanhM((CGMatrix*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}
			break;
		}

		case CGUOT_INV:{
			char msg[MAX_ERR_FMT_LEN];
			snprintf(msg, MAX_ERR_FMT_LEN, "Operation `%s` is not implemented/supported", getUnaryOperationTypeString(type));
			newres = returnResultError(graph, CGET_OPERATION_NOT_IMPLEMENTED, parentNode, msg);
		}

		case CGUOT_TRANSPOSE:{

			if(uhsType == CGVT_DOUBLE){
				newres = transposeD((CGDouble*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if(uhsType == CGVT_VECTOR){
				newres = transposeV((CGVector*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if(uhsType == CGVT_MATRIX){
				newres = transposeM((CGMatrix*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			break;
		}

		case CGUOT_RELU:{
			if(uhsType == CGVT_DOUBLE){
				newres = reluD((CGDouble*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if(uhsType == CGVT_VECTOR){
				newres = reluV((CGVector*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if(uhsType == CGVT_MATRIX){
				newres = reluM((CGMatrix*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			break;
		}

        case CGUOT_SOFTPLUS:{
			if(uhsType == CGVT_DOUBLE){
				newres = softplusD((CGDouble*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if(uhsType == CGVT_VECTOR){
				newres = softplusV((CGVector*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if(uhsType == CGVT_MATRIX){
				newres = softplusM((CGMatrix*)uhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			break;
        }
	}
	char msg[MAX_ERR_FMT_LEN];
	snprintf(msg, MAX_ERR_FMT_LEN, "Operation [%s %s] cannot be applied", getVariableTypeString(uhsType), getUnaryOperationTypeString(type));
	newres = returnResultError(graph, CGET_INCOMPATIBLE_ARGS_EXCEPTION, parentNode, msg);
	return newres;
}

CGResultNode* processBinaryOperation(CGraph* graph, CGBinaryOperationType type, CGNode* lhs, CGNode* rhs, CGNode* parentNode){
	CGVarType lhsType = CGVT_DOUBLE;
	CGVarType rhsType = CGVT_DOUBLE;
	CGResultNode* newres = NULL;
	void* lhsValue = NULL;
	void* rhsValue = NULL;

	CGResultNode* lhsResult = computeCGNode(graph, lhs);
	CHECK_RESULT(lhsResult)
	lhsType = lhsResult->type;
	lhsValue = lhsResult->value;

	CGResultNode* rhsResult = computeCGNode(graph, rhs);
	CHECK_RESULT(rhsResult)
	rhsType = rhsResult->type;
	rhsValue = rhsResult->value;

	switch(type){
		case CGBOT_ADD:{
			if((lhsType == CGVT_DOUBLE) && (rhsType == CGVT_DOUBLE)){
				newres = addDD((CGDouble*)lhsValue, (CGDouble*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_VECTOR) && (rhsType == CGVT_DOUBLE)){
				newres = addVD((CGVector*)lhsValue, (CGDouble*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_DOUBLE) && (rhsType == CGVT_VECTOR)){
				newres = addVD((CGVector*)rhsValue, (CGDouble*)lhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_MATRIX) && (rhsType == CGVT_DOUBLE)){
				newres = addMD((CGMatrix*)lhsValue, (CGDouble*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_DOUBLE) && (rhsType == CGVT_MATRIX)){
				newres = addMD((CGMatrix*)rhsValue, (CGDouble*)lhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_VECTOR) && (rhsType == CGVT_VECTOR)){
				newres = addVV((CGVector*)lhsValue, (CGVector*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_MATRIX) && (rhsType == CGVT_VECTOR)){
				newres = addMV((CGMatrix*)lhsValue, (CGVector*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_VECTOR) && (rhsType == CGVT_MATRIX)){
				newres = addMV((CGMatrix*)rhsValue, (CGVector*)lhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_MATRIX) && (rhsType == CGVT_MATRIX)){
				newres = addMM((CGMatrix*)lhsValue, (CGMatrix*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}
			break;
		}

		case CGBOT_SUB:{
			if((lhsType == CGVT_DOUBLE) && (rhsType == CGVT_DOUBLE)){
				newres = subDD((CGDouble*)lhsValue, (CGDouble*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_VECTOR) && (rhsType == CGVT_DOUBLE)){
				newres = subVD((CGVector*)lhsValue, (CGDouble*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_DOUBLE) && (rhsType == CGVT_VECTOR)){
				newres = subDV((CGDouble*)lhsValue, (CGVector*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_MATRIX) && (rhsType == CGVT_DOUBLE)){
				newres = subMD((CGMatrix*)lhsValue, (CGDouble*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_DOUBLE) && (rhsType == CGVT_MATRIX)){
				newres = subDM((CGDouble*)lhsValue, (CGMatrix*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_VECTOR) && (rhsType == CGVT_VECTOR)){
				newres = subVV((CGVector*)lhsValue, (CGVector*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_MATRIX) && (rhsType == CGVT_MATRIX)){
				newres = subMM((CGMatrix*)lhsValue, (CGMatrix*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_MATRIX) && (rhsType == CGVT_VECTOR)){
				newres = subMV((CGMatrix*)lhsValue, (CGVector*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_VECTOR) && (rhsType == CGVT_MATRIX)){
				newres = subVM((CGVector*)lhsValue, (CGMatrix*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}
			break;
		}

		case CGBOT_DIV:{
			if((lhsType == CGVT_DOUBLE) && (rhsType == CGVT_DOUBLE)){
				newres = divDD((CGDouble*)lhsValue, (CGDouble*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_VECTOR) && (rhsType == CGVT_DOUBLE)){
				newres = divVD((CGVector*)lhsValue, (CGDouble*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_DOUBLE) && (rhsType == CGVT_VECTOR)){
				newres = divDV((CGDouble*)lhsValue, (CGVector*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_VECTOR) && (rhsType == CGVT_VECTOR)){
				newres = divVV((CGVector*)lhsValue, (CGVector*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_MATRIX) && (rhsType == CGVT_DOUBLE)){
				newres = divMD((CGMatrix*)lhsValue, (CGDouble*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_MATRIX) && (rhsType == CGVT_VECTOR)){
				newres = divMV((CGMatrix*)lhsValue, (CGVector*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_DOUBLE) && (rhsType == CGVT_MATRIX)){
				newres = divDM((CGDouble*)lhsValue, (CGMatrix*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			break;
		}

		case CGBOT_MULT:{
			// TODO: update
			if((lhsType == CGVT_MATRIX) && (rhsType == CGVT_MATRIX)){
				newres = mulMM((CGMatrix*)lhsValue, (CGMatrix*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			// TODO: update
			if((lhsType == CGVT_MATRIX) && (rhsType == CGVT_VECTOR)){
				newres = mulMV((CGMatrix*)lhsValue, (CGVector*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			// TODO: update
			if((lhsType == CGVT_VECTOR) && (rhsType == CGVT_MATRIX)){
				newres = mulMV((CGMatrix*)rhsValue, (CGVector*)lhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_DOUBLE) && (rhsType == CGVT_DOUBLE)){
				newres = mulDD((CGDouble*)lhsValue, (CGDouble*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_DOUBLE) && (rhsType == CGVT_VECTOR)){
				newres = mulDV((CGDouble*)lhsValue, (CGVector*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_VECTOR) && (rhsType == CGVT_DOUBLE)){
				newres = mulDV((CGDouble*)rhsValue, (CGVector*)lhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_DOUBLE) && (rhsType == CGVT_MATRIX)){
				newres = mulDM((CGDouble*)lhsValue, (CGMatrix*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_MATRIX) && (rhsType == CGVT_DOUBLE)){
				newres = mulDM((CGDouble*)rhsValue, (CGMatrix*)lhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_VECTOR) && (rhsType == CGVT_VECTOR)){
				newres = crossVV((CGVector*)rhsValue, (CGVector*)lhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}
			break;
		}

		case CGBOT_POW:{
			if((lhsType == CGVT_DOUBLE) && (rhsType == CGVT_DOUBLE)){
				newres = powDD((CGDouble*)lhsValue, (CGDouble*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_VECTOR) && (rhsType == CGVT_DOUBLE)){
				newres = powVD((CGVector*)lhsValue, (CGDouble*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_MATRIX) && (rhsType == CGVT_DOUBLE)){
				newres = powMD((CGMatrix*)lhsValue, (CGDouble*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}
			break;
		}

		case CGBOT_DOT: {
			/*
			 * The following are the same as MUL
			 */
			if((lhsType == CGVT_DOUBLE) && (rhsType == CGVT_DOUBLE)){
				newres = mulDD((CGDouble*)lhsValue, (CGDouble*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_DOUBLE) && (rhsType == CGVT_VECTOR)){
				newres = mulDV((CGDouble*)lhsValue, (CGVector*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_VECTOR) && (rhsType == CGVT_DOUBLE)){
				newres = mulDV((CGDouble*)rhsValue, (CGVector*)lhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_DOUBLE) && (rhsType == CGVT_MATRIX)){
				newres = mulDM((CGDouble*)lhsValue, (CGMatrix*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_MATRIX) && (rhsType == CGVT_DOUBLE)){
				newres = mulDM((CGDouble*)rhsValue, (CGMatrix*)lhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			/* here starts dot specific impl */

			if((lhsType == CGVT_MATRIX) && (rhsType == CGVT_MATRIX)){
				newres = dotMM((CGMatrix*)lhsValue, (CGMatrix*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_MATRIX) && (rhsType == CGVT_VECTOR)){
				newres = dotMV((CGMatrix*)lhsValue, (CGVector*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			if((lhsType == CGVT_VECTOR) && (rhsType == CGVT_MATRIX)){
				newres = dotVM((CGVector*)lhsValue, (CGMatrix*)rhsValue, graph, parentNode);
				//newres = dotMM(vectorToMatrix((CGVector*)lhsValue), (CGMatrix*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}

			/*
			if((lhsType == CGVT_VECTOR) && (rhsType == CGVT_MATRIX)){
				newres = mulMV((CGMatrix*)rhsValue, (CGVector*)lhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}
			*/

			if((lhsType == CGVT_VECTOR) && (rhsType == CGVT_VECTOR)){
				newres = dotVV((CGVector*)lhsValue, (CGVector*)rhsValue, graph, parentNode);
				parentNode->result = newres;
				return newres;
			}
			break;
		}

		case CGBOT_TMULT:{
			char msg[MAX_ERR_FMT_LEN];
			snprintf(msg, MAX_ERR_FMT_LEN, "Operation `%s` is not implemented/supported", getBinaryOperationTypeString(type));
			return returnResultError(graph, CGET_OPERATION_NOT_IMPLEMENTED, parentNode, msg);
		}
	}

	char msg[MAX_ERR_FMT_LEN];
	snprintf(msg, MAX_ERR_FMT_LEN, "Operation [%s %s %s] cannot be applied", getVariableTypeString(lhsType), getBinaryOperationTypeString(type), getVariableTypeString(rhsType));
	return returnResultError(graph, CGET_INCOMPATIBLE_ARGS_EXCEPTION, parentNode, msg);
}

CGResultNode* computeRawNode(CGNode* node){
	CGraph* tmp_G = makeGraph("temp_g");

	tmp_G->root = node;
	storeNodesInGraph(tmp_G, node);

	CGResultNode* res = computeGraph(tmp_G);
	res = copyResultNode(res);
	freeGraph(tmp_G);
	free(tmp_G);

	return res;
}

CGResultNode* computeCGNode(CGraph* graph, CGNode* node){
	CGResultNode* result = NULL;

	if(node->result != NULL){
		return node->result;
	}

	switch(node->type){
		case CGNT_CONSTANT:{
			result = constantNodeToResultNodeCopy(node);
			break;
		}

		case CGNT_VARIABLE:{
			if(graph == NULL)
			{
				char msg[MAX_ERR_FMT_LEN];
				snprintf(msg, MAX_ERR_FMT_LEN, "Cannot compute variable`%s` without the graph instance", node->var->name);
				return returnResultError(graph, CGET_NO_GRAPH_INSTANCE, node, msg);
			}

			CGNode* constantNode = *map_get(&graph->vars, node->var->name);
			if(constantNode == NULL)
			{
				char msg[MAX_ERR_FMT_LEN];
				snprintf(msg, MAX_ERR_FMT_LEN, "No variable `%s` was found in graph `%s`", node->var->name, graph!=NULL?graph->name:"[anonymous]");
				return returnResultError(graph, CGET_VARIABLE_DOES_NOT_EXIST, node, msg);
			}

			CGResultNode* rnode = computeCGNode(graph, constantNode);
			CHECK_RESULT(rnode)
			constantNode->result = rnode;
			node->result = copyResultNode(rnode);

			result = node->result;
			break;
		}
		case CGNT_BINARY_OPERATION:
			result = processBinaryOperation(graph, node->bop->type, node->bop->lhs, node->bop->rhs, node);
			break;
		case CGNT_UNARY_OPERATION:
			result = processUnaryOperation(graph, node->uop->type, node->uop->uhs, node);
			break;

		/*
		 * TODO: add this test to unittest
		 */
		case CGNT_AXIS_BOUND_OPERATION:
		{
            //CGResultNode* res = computeCGNode(graph, node->axop->uhs);

			switch(node->axop->type){
				case CGABOT_SUM:{
					CGResultNode* newres = computeCGNode(graph, node->axop->uhs);
					CHECK_RESULT(newres)

					if(newres->type == CGVT_DOUBLE){
						result = sumD((CGDouble*)newres->value, graph, node);
					}

					if(newres->type == CGVT_VECTOR){
						result = sumV((CGVector*)newres->value, graph, node);
					}

					if(newres->type == CGVT_MATRIX){
						result = sumM((CGMatrix*)newres->value, graph, node, node->axop->axis);
					}

					break;
				}

				case CGABOT_MAX:{
					result = max(node, graph);
					break;
				}

				case CGABOT_MIN:{
					result = min(node, graph);
					break;
				}

				case CGABOT_MEAN:{
					result = mean(node, graph);
					break;
				}

				case CGABOT_SOFTMAX:{
					char msg[MAX_ERR_FMT_LEN];
					snprintf(msg, MAX_ERR_FMT_LEN, "Operation [CGABOT_SOFTMAX] is not implemented/supported");
					return returnResultError(graph, CGET_OPERATION_NOT_IMPLEMENTED, node, msg);
					//CGResultNode* res = computeCGNode(graph, node->axop->uhs);
					//break;
				}

				case CGABOT_ARGMAX:{
				    result = argmax(node, graph);
                    break;
				}

                case CGABOT_ARGMIN:{
                    result = argmin(node, graph);
                    break;
                }

				default:
				    break;
			}

			break;
		}
		case CGNT_GRAPH:{
			result = computeGraph(node->graph);
			break;
			/*
			char msg[MAX_ERR_FMT_LEN];
			snprintf(msg, MAX_ERR_FMT_LEN, "Operation [GRAPH] is not implemented/supported");
			return returnResultError(graph, CGET_OPERATION_NOT_IMPLEMENTED, node, msg);
			*/
		}

		case CGNT_CROSS_ENTROPY_LOSS_FUNC:
		{
			CGResultNode* X_res = computeCGNode(graph, node->crossEntropyLoss->x);
			CGNode* X = resultNodeToConstantNodeCopy(X_res);

			CGNode* X_val = softmax_node(X);

			CGResultNode* x_res = computeRawNode(X_val);

			CGResultNode* y_res = computeCGNode(graph, node->crossEntropyLoss->y);

			result = crossEntropy(x_res, y_res, node->crossEntropyLoss->num_classes);

			freeResultNode(x_res);
			free(x_res);
			break;
		}
	}

	node->result = reduceDim(result);

	if(node->diff == NULL)
		switch(node->result->type){
			case CGVT_DOUBLE:
				node->diff = makeZeroDoubleConstantNode();
				break;
			case CGVT_VECTOR:{
				CGVector* v = (CGVector*)node->result ->value;
				node->diff = makeZeroVectorConstantNode(v->len);
				break;
			}

			case CGVT_MATRIX:{
				CGMatrix* v = (CGMatrix*)node->result ->value;
				node->diff = makeZeroMatrixConstantNode(v->rows, v->cols);
				break;
			}
		}

	return node->result;
}

CGResultNode* reduceDim(CGResultNode* result){
    //return result;

	switch(result->type){
		case CGVT_DOUBLE:{
			return result;
		}

		case CGVT_VECTOR:{
			CGVector* vec = (CGVector*)result->value;
			if (vec->len > 1)
				return result;

#ifdef CG_USE_OPENCL
			copyDataToHost(result);
#endif

			CGDouble* d = calloc(1, sizeof(CGDouble));
			d->value = vec->data[0];

			freeVectorValue(result->value);
			free(result->value);

			result->type = CGVT_DOUBLE;
			result->value = d;

			return result;
		}

		case CGVT_MATRIX:{
			CGMatrix* mat = (CGMatrix*)result->value;

			if((mat->rows>1) &&(mat->cols>1))
				return result;

			if((mat->rows == 1) && (mat->cols == 1)){
                copyDataToHost(result);

				CGDouble* d = calloc(1, sizeof(CGDouble));
				d->value = mat->data[0];

				freeMatrixValue(result->value);
				free(result->value);

				result->type = CGVT_DOUBLE;
				result->value = d;

				return result;
			}

			if(mat->rows == 1){
				CGVector* vec = calloc(1, sizeof(CGVector));
				vec->len = mat->cols;
				vec->data = mat->data;
				//TODO: Memory leak here
#ifdef CG_USE_OPENCL
				vec->buf = mat->buf;
#endif

				//freeMatrixValue(result->value);
				free(result->value);

				result->type = CGVT_VECTOR;
				result->value = vec;

				return result;
			}

			return result;
		}
	}
}

CGResultNode* computeGraph(CGraph* graph){
	resetGraphResultNodes(graph, graph->root);
	CGResultNode* res = computeCGNode(graph, graph->root);
	copyDataToHost(res);

	return res;
}

CGResultNode* computeGraphNode(CGraph* graph, CGNode* node){
	storeNodesInGraph(graph, node);
	resetGraphResultNodes(graph, node);
	return computeCGNode(graph, node);
}

void storeNodesInGraph(CGraph* graph, CGNode* node){
	int idx = -1;

	vec_find(&graph->nodes, node, idx);

	if(idx != -1)
		return;

	vec_push(&graph->nodes, node);

	switch(node->type){
		case CGNT_CONSTANT:
			break;
		case CGNT_VARIABLE:
			break;
		case CGNT_BINARY_OPERATION:
			storeNodesInGraph(graph, node->bop->lhs);
			storeNodesInGraph(graph, node->bop->rhs);
			break;
		case CGNT_UNARY_OPERATION:
			storeNodesInGraph(graph, node->uop->uhs);
			break;
		case CGNT_AXIS_BOUND_OPERATION:
			storeNodesInGraph(graph, node->axop->uhs);
			break;
		case CGNT_GRAPH:
			storeNodesInGraph(graph, node->graph->root);
			break;
		case CGNT_CROSS_ENTROPY_LOSS_FUNC:
			storeNodesInGraph(graph, node->crossEntropyLoss->x);
			storeNodesInGraph(graph, node->crossEntropyLoss->y);
			break;
	}
}

void resetGraphResultNodes(CGraph* graph, CGNode* node){
	if(node->result != NULL){
		freeResultNode(node->result);
		free(node->result);
		node->result = NULL;
	}

	if(node->diff != NULL){
		freeNode(graph, node->diff);
		free(node->diff);
		node->diff = NULL;
	}

	switch(node->type){
		case CGNT_CONSTANT:
			break;
		case CGNT_VARIABLE:{
			CGNode* var = graphGetVar(graph, node->var->name);
			if(var != NULL)
				resetGraphResultNodes(graph, var);
			break;
		}
		case CGNT_BINARY_OPERATION:
			resetGraphResultNodes(graph, node->bop->lhs);
			resetGraphResultNodes(graph, node->bop->rhs);
			break;
		case CGNT_UNARY_OPERATION:
			resetGraphResultNodes(graph, node->uop->uhs);
			break;
		case CGNT_AXIS_BOUND_OPERATION:
			resetGraphResultNodes(graph, node->axop->uhs);
			break;
		case CGNT_GRAPH:
			resetGraphResultNodes(graph, node->graph->root);
			break;
		case CGNT_CROSS_ENTROPY_LOSS_FUNC:
			resetGraphResultNodes(graph, node->crossEntropyLoss->x);
			resetGraphResultNodes(graph, node->crossEntropyLoss->y);
			break;
	}
}

void graphSetVar_lua(CGraph* graph, const char* name, CGNode* value){
	CGNode** old = map_get(&graph->vars, name);
	if(old != NULL){
		//freeNode(graph, *old);
		map_remove(&graph->vars, name);
	}

	int res = map_set(&graph->vars, name, value);
}

ext: c

max_stars_repo_name: Enehcruon/cgraph

max_issues_repo_name: Enehcruon/cgraph

max_forks_repo_name: praisethemoon/ccgraph


text:
#ifndef LIBTENSOR_CBLAS_H_H
#define LIBTENSOR_CBLAS_H_H

extern "C" {  // Fixes older cblas.h versions without extern "C"
#include <cblas.h>
}

#endif  // LIBTENSOR_CBLAS_H_H

ext: h

max_stars_repo_name: maxscheurer/libtensor

max_issues_repo_name: maxscheurer/libtensor

max_forks_repo_name: maxscheurer/libtensor


text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#pragma once

#ifndef painterpath_d2d84280_c369_48f2_99f8_28652ad13baa_h
#define painterpath_d2d84280_c369_48f2_99f8_28652ad13baa_h

#include <ariel/config.h>
#include <gslib/std.h>

__ariel_begin__

class painter_path;
typedef list<painter_path> painter_paths;

class ariel_export painter_linestrip
{
public:
    typedef vector<vec2> points;
    typedef points::iterator ptiter;
    typedef points::const_iterator ptciter;

protected:
    bool                _closed;
    points              _pts;

public:
    painter_linestrip();
    ~painter_linestrip() {}
    int get_size() const { return (int)_pts.size(); }
    void get_bound_rect(rectf& rc) const;
    vec2& get_point(int i) { return _pts.at(i); }
    const vec2& get_point(int i) const { return _pts.at(i); }
    vec2& get_last_point() { return _pts.back(); }
    const vec2& get_last_point() const { return _pts.back(); }
    void add_point(const vec2& pt) { _pts.push_back(pt); }
    void set_point(int i, const vec2& pt) { _pts.at(i) = pt; }
    void reserve(int s) { _pts.reserve(s); }
    void clear() { _pts.clear(); }
    void swap(painter_linestrip& another);
    void finish();
    void transform(const mat3& m);
    vec2* expand(int size);
    void expand_to(int size) { _pts.resize(size); }
    void reverse();
    void offset(const painter_linestrip& org, float off);
    int point_inside(const vec2& pt) const;
    bool is_closed() const { return _closed; }
    void set_closed(bool c) { _closed = c; }
    bool is_clockwise() const;
    bool is_convex() const;
    bool is_convex(int i) const;
    void tracing() const;
    void trace_segments() const;
    void trace_mel(float z) const;

protected:
    enum
    {
        st_clockwise_mask = 1,
        st_convex_mask = 2,
    };
    mutable uint        _tst_table;
    mutable bool        _is_clock_wise;
    mutable bool        _is_convex;

protected:
    bool is_clockwise_inited() const { return (_tst_table & st_clockwise_mask) != 0; }
    void set_clockwise_inited() const { _tst_table |= st_clockwise_mask; }
    bool is_convex_inited() const { return (_tst_table & st_convex_mask) != 0; }
    void set_convex_inited() const { _tst_table |= st_convex_mask; }
};

typedef list<painter_linestrip> linestrips;
typedef vector<painter_linestrip*> linestripvec;
typedef linestrips painter_linestrips;

/* create a random access view for linestrips */
ariel_export extern void append_linestrips_rav(linestripvec& rav, linestrips& src);
ariel_export extern void create_linestrips_rav(linestripvec& rav, linestrips& src);

class ariel_export painter_path
{
public:
    enum tag
    {
        pt_moveto,
        pt_lineto,
        pt_quadto,
        pt_cubicto,
    };

    class __gs_novtable ariel_export node abstract
    {
    public:
        virtual ~node() {}
        virtual tag get_tag() const = 0;
        virtual const vec2& get_point() const = 0;
        virtual void set_point(const vec2&) = 0;
        virtual void interpolate(painter_linestrip& c, const node* last, float step_len) const = 0;

    public:
        template<class _c>
        _c* as_node() { return static_cast<_c*>(this); }
        template<class _c>
        const _c* as_const_node() const { return static_cast<const _c*>(this); }
    };

    template<tag _tag>
    class node_tpl:
        public node
    {
    public:
        node_tpl() {}
        node_tpl(float x, float y): _pt(x, y) {}
        node_tpl(const vec2& pt): _pt(pt) {}

    public:
        tag get_tag() const override { return _tag; }
        const vec2& get_point() const override { return _pt; }
        void set_point(const vec2& pt) override { _pt = pt; }
        void interpolate(painter_linestrip& c, const node*, float) const override { c.add_point(_pt); }

    protected:
        vec2            _pt;
    };

    typedef node_tpl<pt_moveto> move_to_node;
    typedef node_tpl<pt_lineto> line_to_node;

    class ariel_export quad_to_node:
        public node_tpl<pt_quadto>
    {
    public:
        quad_to_node() {}
        quad_to_node(const vec2& p1, const vec2& p2): node_tpl(p2) { _c = p1; }
        void set_control(const vec2& c) { _c = c; }
        const vec2& get_control() const { return _c; }
        void interpolate(painter_linestrip& c, const node* last, float step_len) const override;

    protected:
        vec2            _c;
    };

    class ariel_export cubic_to_node:
        public node_tpl<pt_cubicto>
    {
    public:
        cubic_to_node() {}
        cubic_to_node(const vec2& p1, const vec2& p2, const vec2 p3): node_tpl(p3) { _c[0] = p1, _c[1] = p2; }
        void set_control1(const vec2& c) { _c[0] = c; }
        void set_control2(const vec2& c) { _c[1] = c; }
        const vec2& get_control1() const { return _c[0]; }
        const vec2& get_control2() const { return _c[1]; }
        void interpolate(painter_linestrip& c, const node* last, float step_len) const override;

    protected:
        vec2            _c[2];
    };

    typedef vector<node*> node_list;
    typedef node_list::iterator iterator;
    typedef node_list::const_iterator const_iterator;
    typedef vector<int> indices;

protected:
    node_list           _nodelist;

public:
    painter_path() {}
    painter_path(const painter_path& path) { duplicate(path); }
    ~painter_path() { destroy(); }
    bool empty() const { return _nodelist.empty(); }
    int size() const { return (int)_nodelist.size(); }
    void resize(int len);
    void destroy();
    void duplicate(const painter_path& path);
    void add_path(const painter_path& path);
    void add_rect(const rectf& rc);
    void swap(painter_path& path);
    iterator begin() { return _nodelist.begin(); }
    iterator end() { return _nodelist.end(); }
    const_iterator begin() const { return _nodelist.begin(); }
    const_iterator end() const { return _nodelist.end(); }
    node* get_node(int i) { return _nodelist.at(i); }
    const node* get_node(int i) const { return _nodelist.at(i); }
    void close_path();
    void close_sub_path();
    void get_boundary_box(rectf& rc) const;
    void move_to(float x, float y) { move_to(vec2(x, y)); }
    void line_to(float x, float y) { line_to(vec2(x, y)); }
    void arc_to(float x, float y, float r) { arc_to(vec2(x, y), r); }
    void rarc_to(float x, float y, float r) { rarc_to(vec2(x, y), r); }
    void quad_to(float x1, float y1, float x2, float y2) { quad_to(vec2(x1, y1), vec2(x2, y2)); }
    void cubic_to(float x1, float y1, float x2, float y2, float x3, float y3) { cubic_to(vec2(x1, y1), vec2(x2, y2), vec2(x3, y3)); }
    void move_to(const vec2& pt) { _nodelist.push_back(new move_to_node(pt)); }
    void line_to(const vec2& pt) { _nodelist.push_back(new line_to_node(pt)); }
    void quad_to(const vec2& p1, const vec2& p2) { _nodelist.push_back(new quad_to_node(p1, p2)); }
    void cubic_to(const vec2& p1, const vec2& p2, const vec2& p3) { _nodelist.push_back(new cubic_to_node(p1, p2, p3)); }
    void arc_to(const vec2& p1, const vec2& p2, float r);
    void arc_to(const vec2& pt, float r);
    void rarc_to(const vec2& pt, float r);
    void transform(const mat3& m);
    void get_linestrips(linestrips& c, float step_len = -1.f) const;
    int get_control_contour(painter_linestrip& ls, int start) const;
    int get_sub_path(painter_path& sp, int start) const;
    void to_sub_paths(painter_paths& paths) const;
    bool is_clockwise() const;
    bool is_convex() const;
    void simplify(painter_path& path) const;
    void reverse();
    void tracing() const;
    void tracing_segments() const;
};

typedef painter_path::node painter_node;

struct ariel_export path_info
{
    const painter_node*     node[2];

    path_info() { node[0] = node[1] = 0; }
    path_info(const painter_node* n1, const painter_node* n2) { node[0] = n1, node[1] = n2; }
    int get_order() const;
    int get_point_count() const;
    bool get_points(vec2 pt[], int cnt) const;
    void tracing() const;
};

enum curve_type
{
    ct_quad,
    ct_cubic,
};

struct __gs_novtable ariel_export curve_splitter abstract
{
    vec2                fixedpt;
    float               ratio;
    curve_splitter*     child[2];
    curve_splitter*     parent;

public:
    curve_splitter();
    virtual ~curve_splitter();
    virtual curve_type get_type() const = 0;
    virtual int get_point_count() const = 0;
    virtual bool get_points(vec2 p[], int count) const = 0;
    virtual void split(float t) = 0;
    virtual void interpolate(vec2& p, float t) const = 0;
    virtual float reparameterize(const vec2& p) const = 0;
    virtual void tracing() const = 0;
    bool is_leaf() const;
};

struct ariel_export curve_splitter_quad:
    public curve_splitter
{
    vec2                cp[3];
    vec3                para[2];

public:
    curve_splitter_quad(const vec2 p[3]);
    curve_type get_type() const override { return ct_quad; }
    int get_point_count() const override { return 3; }
    bool get_points(vec2 p[], int count) const override;
    void split(float t) override;
    void interpolate(vec2& p, float t) const override;
    float reparameterize(const vec2& p) const override;
    void tracing() const override;
};

struct ariel_export curve_splitter_cubic:
    public curve_splitter
{
    vec2                cp[4];
    vec4                para[2];

public:
    curve_splitter_cubic(const vec2 p[4]);
    curve_type get_type() const override { return ct_cubic; }
    int get_point_count() const override { return 4; }
    bool get_points(vec2 p[], int count) const override;
    void split(float t) override;
    void interpolate(vec2& p, float t) const override;
    float reparameterize(const vec2& p) const override;
    void tracing() const override;
};

struct ariel_export curve_helper
{
    static curve_splitter* create_splitter(const path_info* pnf);
    static curve_splitter* create_next_splitter(vec2& p, curve_splitter* cs, float t);
    static curve_splitter* query_splitter(curve_splitter* cs, float t);
    static curve_splitter* query_splitter(curve_splitter* cs, const vec2& p);
    static curve_splitter* query_splitter(curve_splitter* cs, const vec2& p1, const vec2& p2);
};

struct ariel_export painter_helper
{
    enum
    {
        merge_straight_line = 0x01,
        fix_loop = 0x02,
        fix_inflection = 0x04,
        reduce_straight_curve = 0x08,
        reduce_short_line = 0x10,
    };

    static void transform(painter_path& output, const painter_path& path, uint mask);
    static void close_sub_paths(painter_path& output, const painter_path& path);
};

__ariel_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib