text:
/* matrix/gsl_matrix_uint.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
 *
 * 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.
 */

#ifndef __GSL_MATRIX_UINT_H__
#define __GSL_MATRIX_UINT_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_uint.h>
#include <gsl/gsl_blas_types.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  unsigned int * data;
  gsl_block_uint * block;
  int owner;
} gsl_matrix_uint;

typedef struct
{
  gsl_matrix_uint matrix;
} _gsl_matrix_uint_view;

typedef _gsl_matrix_uint_view gsl_matrix_uint_view;

typedef struct
{
  gsl_matrix_uint matrix;
} _gsl_matrix_uint_const_view;

typedef const _gsl_matrix_uint_const_view gsl_matrix_uint_const_view;

/* Allocation */

GSL_FUN gsl_matrix_uint *
gsl_matrix_uint_alloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_uint *
gsl_matrix_uint_calloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_uint *
gsl_matrix_uint_alloc_from_block (gsl_block_uint * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

GSL_FUN gsl_matrix_uint *
gsl_matrix_uint_alloc_from_matrix (gsl_matrix_uint * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN gsl_vector_uint *
gsl_vector_uint_alloc_row_from_matrix (gsl_matrix_uint * m,
                                        const size_t i);

GSL_FUN gsl_vector_uint *
gsl_vector_uint_alloc_col_from_matrix (gsl_matrix_uint * m,
                                        const size_t j);

GSL_FUN void gsl_matrix_uint_free (gsl_matrix_uint * m);

/* Views */

GSL_FUN _gsl_matrix_uint_view
gsl_matrix_uint_submatrix (gsl_matrix_uint * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_uint_view
gsl_matrix_uint_row (gsl_matrix_uint * m, const size_t i);

GSL_FUN _gsl_vector_uint_view
gsl_matrix_uint_column (gsl_matrix_uint * m, const size_t j);

GSL_FUN _gsl_vector_uint_view
gsl_matrix_uint_diagonal (gsl_matrix_uint * m);

GSL_FUN _gsl_vector_uint_view
gsl_matrix_uint_subdiagonal (gsl_matrix_uint * m, const size_t k);

GSL_FUN _gsl_vector_uint_view
gsl_matrix_uint_superdiagonal (gsl_matrix_uint * m, const size_t k);

GSL_FUN _gsl_vector_uint_view
gsl_matrix_uint_subrow (gsl_matrix_uint * m, const size_t i,
                         const size_t offset, const size_t n);

GSL_FUN _gsl_vector_uint_view
gsl_matrix_uint_subcolumn (gsl_matrix_uint * m, const size_t j,
                            const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_uint_view
gsl_matrix_uint_view_array (unsigned int * base,
                             const size_t n1,
                             const size_t n2);

GSL_FUN _gsl_matrix_uint_view
gsl_matrix_uint_view_array_with_tda (unsigned int * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_FUN _gsl_matrix_uint_view
gsl_matrix_uint_view_vector (gsl_vector_uint * v,
                              const size_t n1,
                              const size_t n2);

GSL_FUN _gsl_matrix_uint_view
gsl_matrix_uint_view_vector_with_tda (gsl_vector_uint * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_FUN _gsl_matrix_uint_const_view
gsl_matrix_uint_const_submatrix (const gsl_matrix_uint * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_uint_const_view
gsl_matrix_uint_const_row (const gsl_matrix_uint * m,
                            const size_t i);

GSL_FUN _gsl_vector_uint_const_view
gsl_matrix_uint_const_column (const gsl_matrix_uint * m,
                               const size_t j);

GSL_FUN _gsl_vector_uint_const_view
gsl_matrix_uint_const_diagonal (const gsl_matrix_uint * m);

GSL_FUN _gsl_vector_uint_const_view
gsl_matrix_uint_const_subdiagonal (const gsl_matrix_uint * m,
                                    const size_t k);

GSL_FUN _gsl_vector_uint_const_view
gsl_matrix_uint_const_superdiagonal (const gsl_matrix_uint * m,
                                      const size_t k);

GSL_FUN _gsl_vector_uint_const_view
gsl_matrix_uint_const_subrow (const gsl_matrix_uint * m, const size_t i,
                               const size_t offset, const size_t n);

GSL_FUN _gsl_vector_uint_const_view
gsl_matrix_uint_const_subcolumn (const gsl_matrix_uint * m, const size_t j,
                                  const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_uint_const_view
gsl_matrix_uint_const_view_array (const unsigned int * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_FUN _gsl_matrix_uint_const_view
gsl_matrix_uint_const_view_array_with_tda (const unsigned int * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_FUN _gsl_matrix_uint_const_view
gsl_matrix_uint_const_view_vector (const gsl_vector_uint * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN _gsl_matrix_uint_const_view
gsl_matrix_uint_const_view_vector_with_tda (const gsl_vector_uint * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_FUN void gsl_matrix_uint_set_zero (gsl_matrix_uint * m);
GSL_FUN void gsl_matrix_uint_set_identity (gsl_matrix_uint * m);
GSL_FUN void gsl_matrix_uint_set_all (gsl_matrix_uint * m, unsigned int x);

GSL_FUN int gsl_matrix_uint_fread (FILE * stream, gsl_matrix_uint * m) ;
GSL_FUN int gsl_matrix_uint_fwrite (FILE * stream, const gsl_matrix_uint * m) ;
GSL_FUN int gsl_matrix_uint_fscanf (FILE * stream, gsl_matrix_uint * m);
GSL_FUN int gsl_matrix_uint_fprintf (FILE * stream, const gsl_matrix_uint * m, const char * format);

GSL_FUN int gsl_matrix_uint_memcpy(gsl_matrix_uint * dest, const gsl_matrix_uint * src);
GSL_FUN int gsl_matrix_uint_swap(gsl_matrix_uint * m1, gsl_matrix_uint * m2);
GSL_FUN int gsl_matrix_uint_tricpy(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix_uint * dest, const gsl_matrix_uint * src);

GSL_FUN int gsl_matrix_uint_swap_rows(gsl_matrix_uint * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_uint_swap_columns(gsl_matrix_uint * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_uint_swap_rowcol(gsl_matrix_uint * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_uint_transpose (gsl_matrix_uint * m);
GSL_FUN int gsl_matrix_uint_transpose_memcpy (gsl_matrix_uint * dest, const gsl_matrix_uint * src);
GSL_FUN int gsl_matrix_uint_transpose_tricpy (CBLAS_UPLO_t Uplo_src, CBLAS_DIAG_t Diag, gsl_matrix_uint * dest, const gsl_matrix_uint * src);

GSL_FUN unsigned int gsl_matrix_uint_max (const gsl_matrix_uint * m);
GSL_FUN unsigned int gsl_matrix_uint_min (const gsl_matrix_uint * m);
GSL_FUN void gsl_matrix_uint_minmax (const gsl_matrix_uint * m, unsigned int * min_out, unsigned int * max_out);

GSL_FUN void gsl_matrix_uint_max_index (const gsl_matrix_uint * m, size_t * imax, size_t *jmax);
GSL_FUN void gsl_matrix_uint_min_index (const gsl_matrix_uint * m, size_t * imin, size_t *jmin);
GSL_FUN void gsl_matrix_uint_minmax_index (const gsl_matrix_uint * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

GSL_FUN int gsl_matrix_uint_equal (const gsl_matrix_uint * a, const gsl_matrix_uint * b);

GSL_FUN int gsl_matrix_uint_isnull (const gsl_matrix_uint * m);
GSL_FUN int gsl_matrix_uint_ispos (const gsl_matrix_uint * m);
GSL_FUN int gsl_matrix_uint_isneg (const gsl_matrix_uint * m);
GSL_FUN int gsl_matrix_uint_isnonneg (const gsl_matrix_uint * m);

GSL_FUN unsigned int gsl_matrix_uint_norm1 (const gsl_matrix_uint * m);

GSL_FUN int gsl_matrix_uint_add (gsl_matrix_uint * a, const gsl_matrix_uint * b);
GSL_FUN int gsl_matrix_uint_sub (gsl_matrix_uint * a, const gsl_matrix_uint * b);
GSL_FUN int gsl_matrix_uint_mul_elements (gsl_matrix_uint * a, const gsl_matrix_uint * b);
GSL_FUN int gsl_matrix_uint_div_elements (gsl_matrix_uint * a, const gsl_matrix_uint * b);
GSL_FUN int gsl_matrix_uint_scale (gsl_matrix_uint * a, const double x);
GSL_FUN int gsl_matrix_uint_scale_rows (gsl_matrix_uint * a, const gsl_vector_uint * x);
GSL_FUN int gsl_matrix_uint_scale_columns (gsl_matrix_uint * a, const gsl_vector_uint * x);
GSL_FUN int gsl_matrix_uint_add_constant (gsl_matrix_uint * a, const double x);
GSL_FUN int gsl_matrix_uint_add_diagonal (gsl_matrix_uint * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_FUN int gsl_matrix_uint_get_row(gsl_vector_uint * v, const gsl_matrix_uint * m, const size_t i);
GSL_FUN int gsl_matrix_uint_get_col(gsl_vector_uint * v, const gsl_matrix_uint * m, const size_t j);
GSL_FUN int gsl_matrix_uint_set_row(gsl_matrix_uint * m, const size_t i, const gsl_vector_uint * v);
GSL_FUN int gsl_matrix_uint_set_col(gsl_matrix_uint * m, const size_t j, const gsl_vector_uint * v);
/***********************************************************************/

/* inline functions if you are using GCC */

GSL_FUN INLINE_DECL unsigned int   gsl_matrix_uint_get(const gsl_matrix_uint * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL void    gsl_matrix_uint_set(gsl_matrix_uint * m, const size_t i, const size_t j, const unsigned int x);
GSL_FUN INLINE_DECL unsigned int * gsl_matrix_uint_ptr(gsl_matrix_uint * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL const unsigned int * gsl_matrix_uint_const_ptr(const gsl_matrix_uint * m, const size_t i, const size_t j);

#ifdef HAVE_INLINE
INLINE_FUN
unsigned int
gsl_matrix_uint_get(const gsl_matrix_uint * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
        }
    }
#endif
  return m->data[i * m->tda + j] ;
}

INLINE_FUN
void
gsl_matrix_uint_set(gsl_matrix_uint * m, const size_t i, const size_t j, const unsigned int x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  m->data[i * m->tda + j] = x ;
}

INLINE_FUN
unsigned int *
gsl_matrix_uint_ptr(gsl_matrix_uint * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (unsigned int *) (m->data + (i * m->tda + j)) ;
}

INLINE_FUN
const unsigned int *
gsl_matrix_uint_const_ptr(const gsl_matrix_uint * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (const unsigned int *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_UINT_H__ */

max_stars_repo_name: mgreter/astrometrylib

max_issues_repo_name: zzpwahaha/VimbaCamJILA

max_forks_repo_name: zzpwahaha/VimbaCamJILA


text:
#pragma once
#include "BooleanNodes.h"
#include "BooleanDAG.h"
#include "NodeVisitor.h"
#include "VarStore.h"
#include <map>
#include "FloatSupport.h"
#include <tuple>
#include <iostream>
#ifndef _NOGSL
#include <gsl/gsl_vector.h>
#else
#include "CustomSolver.h"
#endif
#include "Util.h"
#include "Interface.h"

using namespace std;

class SimpleEvaluator: NodeVisitor
{
    BooleanDAG& bdag;
    map<string, int>& floatCtrls; // Maps float ctrl names to indices with grad vectors
    gsl_vector* ctrls; // ctrl values
    vector<double> distances; // Keeps track of distance metric for boolean nodes
    double MIN_VALUE = 0.001;
    Interface* inputValues;
    int DEFAULT_INPUT = -32;

public:
    SimpleEvaluator(BooleanDAG& bdag_p, map<string, int>& floatCtrls_p);
	virtual ~SimpleEvaluator();

    virtual void visit( SRC_node& node );
    virtual void visit( DST_node& node );
    virtual void visit( CTRL_node& node );
    virtual void visit( PLUS_node& node );
    virtual void visit( TIMES_node& node );
    virtual void visit( ARRACC_node& node );
    virtual void visit( DIV_node& node );
    virtual void visit( MOD_node& node );
    virtual void visit( NEG_node& node );
    virtual void visit( CONST_node& node );
    virtual void visit( LT_node& node );
    virtual void visit( EQ_node& node );
    virtual void visit( AND_node& node );
    virtual void visit( OR_node& node );
    virtual void visit( NOT_node& node );
    virtual void visit( ARRASS_node& node );
    virtual void visit( UFUN_node& node );
    virtual void visit( TUPLE_R_node& node );
    virtual void visit( ASSERT_node& node );

    void setInputs(Interface* inputValues_p);
    void run(const gsl_vector* ctrls_p);

    double getErrorOnConstraint(int nodeid);
    double getSqrtError(bool_node* node);
    double getAssertError(bool_node* node);
    double getBoolCtrlError(bool_node* node);
    double getBoolExprError(bool_node* node);

    double dist(int nodeid) {
        return d(bdag[nodeid]);
    }

    void setvalue(bool_node& bn, double d) {
        distances[bn.id] = d;
    }

    double d(bool_node& bn) {
        return distances[bn.id];
    }

    double d(bool_node* bn) {
        return d(*bn);
    }

    void print() {
        for (int i = 0; i < bdag.size(); i++) {
            cout << bdag[i]->lprint() << endl;
            cout << d(bdag[i]) << endl;
        }
    }

    bool isFloat(bool_node& bn) {
        return (bn.getOtype() == OutType::FLOAT);
    }

    bool isFloat(bool_node* bn) {
        return (bn->getOtype() == OutType::FLOAT);
    }

    int getInputValue(bool_node& bn) {
        if (inputValues->hasValue(bn.id)) {
            int val = inputValues->getValue(bn.id);
            return val;
        } else {
            return DEFAULT_INPUT;
        }
    }

    int getInputValue(bool_node* bn) {
        return getInputValue(*bn);
    }

};

ext: h

max_stars_repo_name: natebragg/sketch-backend

max_issues_repo_name: natebragg/sketch-backend

max_forks_repo_name: natebragg/sketch-backend


text:
#include <geometry/basic_c.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <math.h>
#include <stdio.h>

// author:    Andrew Liew <liew@arch.ethz.ch>
// copyright: Copyright 2018, BLOCK Research Group - ETH Zurich
// license:   MIT License
// email:     liew@arch.ethz.ch

void drx_solver_c(
    double tol,    // Tolerance value.
    int steps,     // Maximum number of steps.
    int summary,   // Print summary at end (1:yes or 0:no).
    int m,         // Number of elements.
    int n,         // Number of nodes.
    int *u,        // Element start node.
    int *v,        // Element end node.
    double *X,     // Nodal co-ordinates.
    double *f0,    // Initial edge forces.
    double *l0,    // Initial edge lengths.
    double *k0,    // Initial edge axial stiffnesses.
    int *ind_c,    // Indices of compression only edges.
    int *ind_t,    // Indices of tension only edges.
    int ind_c_n,   // Length of ind_c.
    int ind_t_n,   // Length of ind_t.
    double *B,     // Constraint conditions Bx, By, Bz.
    double *P,     // Nodal loads Px, Py, Pz.
    double *S,     // Shear forces Sx, Sy, Sz.
    int *rows,     // Rows of Ct.
    int *cols,     // Columns of Ct.
    double *vals,  // Values of Ct.
    int nv,        // Length of rows/cols/vals.
    double *M,     // Mass matrix.
    double factor, // Convergence factor.
    double *V,     // Nodal velocities Vx, Vy, Vz.
    int *inds,     // Indices of beam element start nodes.
    int *indi,     // Indices of beam element intermediate nodes.
    int *indf,     // Indices of beam element finish nodes beams.
    double *EIx,   // Nodal EIx flexural stiffnesses.
    double *EIy,   // Nodal EIy flexural stiffnesses.
    int beams,     // Includes beams:1 or not:0.
    int nb)        // Length of inds/indi/indf.

    {

    int a;
    int b;
    int c;
    int i;
    int j;
    int k;
    int nans[nb];
    int ts;

    double alpha;
    double f[m];
    double fx[m];
    double fy[m];
    double fz[m];
    double frx[n];
    double fry[n];
    double frz[n];
    double kappa;
    double l;
    double La;
    double Lb;
    double Lc;
    double LQn;
    double Lmu;
    double Lc1;
    double Lc2;
    double Mi;
    double Ms;
    double q;
    double res;
    double Rx;
    double Ry;
    double Rz;
    double Rn;
    double Uo;
    double Un;
    double xd;
    double yd;
    double zd;

    gsl_vector *Xs = gsl_vector_alloc(3);
    gsl_vector *Xi = gsl_vector_alloc(3);
    gsl_vector *Xf = gsl_vector_alloc(3);
    gsl_vector *Qa = gsl_vector_alloc(3);
    gsl_vector *Qb = gsl_vector_alloc(3);
    gsl_vector *Qc = gsl_vector_alloc(3);
    gsl_vector *Qn = gsl_vector_alloc(3);
    gsl_vector *mu = gsl_vector_alloc(3);
    gsl_vector *ex = gsl_vector_alloc(3);
    gsl_vector *ey = gsl_vector_alloc(3);
    gsl_vector *ez = gsl_vector_alloc(3);
    gsl_vector *K  = gsl_vector_alloc(3);
    gsl_vector *Kx = gsl_vector_alloc(3);
    gsl_vector *Ky = gsl_vector_alloc(3);
    gsl_vector *Mc = gsl_vector_alloc(3);
    gsl_vector *ua = gsl_vector_alloc(3);
    gsl_vector *ub = gsl_vector_alloc(3);
    gsl_vector *c1 = gsl_vector_alloc(3);
    gsl_vector *c2 = gsl_vector_alloc(3);
    gsl_matrix *Sa = gsl_matrix_alloc(nb, 3);
    gsl_matrix *Sb = gsl_matrix_alloc(nb, 3);

    ts  = 0;
    Uo  = 0.;
    res = 1000. * tol;

    while (ts <= steps && res > tol)
    {
        #pragma omp parallel for private(i,j,k,xd,yd,zd,l,q)
        for (i = 0; i < m; i++)
        {
            j = 3 * v[i];
            k = 3 * u[i];
            xd = X[j + 0] - X[k + 0];
            yd = X[j + 1] - X[k + 1];
            zd = X[j + 2] - X[k + 2];
            l  = gsl_hypot3(xd, yd, zd);
            f[i] = f0[i] + k0[i] * (l - l0[i]);
            q = f[i] / l;
            fx[i] = xd * q;
            fy[i] = yd * q;
            fz[i] = zd * q;
        }

        if (ind_t_n > 0)
        {
            #pragma omp parallel for private(i)
            for (i = 0; i < ind_t_n; i++)
            {
                if (f[i] < 0)
                {
                    fx[i] = 0;
                    fy[i] = 0;
                    fz[i] = 0;
                }
            }
        }

        if (ind_c_n > 0)
        {
            #pragma omp parallel for private(i)
            for (i = 0; i < ind_c_n; i++)
            {
                if (f[i] > 0)
                {
                    fx[i] = 0;
                    fy[i] = 0;
                    fz[i] = 0;
                }
            }
        }

        if (beams)
        {
            #pragma omp parallel for private(i,j)
            for (i = 0; i < n; i++)
            {
                j = i * 3;
                S[j + 0] = 0.;
                S[j + 1] = 0.;
                S[j + 2] = 0.;
            }

            for (i = 0; i < nb; i++)
            {
                a = inds[i] * 3;
                b = indi[i] * 3;
                c = indf[i] * 3;

                vector_from_pointer(&X[a], Xs);
                vector_from_pointer(&X[b], Xi);
                vector_from_pointer(&X[c], Xf);
                subtract_vectors(Xi, Xs, Qa);
                subtract_vectors(Xf, Xi, Qb);
                subtract_vectors(Xf, Xs, Qc);
                cross_vectors(Qa, Qb, Qn);
                subtract_vectors(Xf, Xs, mu);
                scale_vector(mu, 0.5);

                La  = length_vector(Qa);
                Lb  = length_vector(Qb);
                Lc  = length_vector(Qc);
                LQn = length_vector(Qn);
                Lmu = length_vector(mu);
                alpha = acos((gsl_pow_2(La) + gsl_pow_2(Lb) - gsl_pow_2(Lc)) / (2. * La * Lb));
                kappa = 2. * sin(alpha) / Lc;

                gsl_vector_memcpy(ex, Qn);
                gsl_vector_memcpy(ez, mu);
                scale_vector(ex, 1./LQn);
                scale_vector(ez, 1./Lmu);
                cross_vectors(ez, ex, ey);
                gsl_vector_memcpy(K, Qn);
                scale_vector(K, kappa/LQn);
                gsl_vector_memcpy(Kx, ex);
                gsl_vector_memcpy(Ky, ey);
                scale_vector(Kx, dot_vectors(K, ex));
                scale_vector(Ky, dot_vectors(K, ey));
                scale_vector(Kx, EIx[i]);
                scale_vector(Ky, EIy[i]);
                add_vectors(Kx, Ky, Mc);
                cross_vectors(Mc, Qa, ua);
                cross_vectors(Mc, Qb, ub);
                normalize_vector(ua);
                normalize_vector(ub);
                cross_vectors(Qa, ua, c1);
                cross_vectors(Qb, ub, c2);

                Lc1 = length_vector(c1);
                Lc2 = length_vector(c2);
                Ms  = length_vector_squared(Mc);
                scale_vector(ua, Ms * Lc1 / (La * dot_vectors(Mc, c1)));
                scale_vector(ub, Ms * Lc2 / (Lb * dot_vectors(Mc, c2)));

                for (j = 0; j < 3; j++)
                {
                    gsl_matrix_set(Sa, i, j, gsl_vector_get(ua, j));
                    gsl_matrix_set(Sb, i, j, gsl_vector_get(ub, j));
                }
            }

            #pragma omp parallel for private(i,j)
            for (i = 0; i < nb; i++)
            {
                nans[i] = 0;
                for (j = 0; j < 3; j++)
                {
                    if (gsl_isnan(gsl_matrix_get(Sa, i, j)) || gsl_isnan(gsl_matrix_get(Sb, i, j)))
                    {
                        nans[i] = 1;
                        break;
                    }
                }
            }

            for (i = 0; i < nb; i++)
            {
                a = inds[i] * 3;
                b = indi[i] * 3;
                c = indf[i] * 3;

                if (nans[i] == 0)
                {
                    for (j = 0; j < 3; j++)
                    {
                        S[a + j] += gsl_matrix_get(Sa, i, j);
                        S[b + j] -= (gsl_matrix_get(Sa, i, j) + gsl_matrix_get(Sb, i, j));
                        S[c + j] += gsl_matrix_get(Sb, i, j);
                    }
                }
            }
        }

        #pragma omp parallel for private(i)
        for (i = 0; i < n; i++)
        {
            frx[i] = 0;
            fry[i] = 0;
            frz[i] = 0;
        }

        for (i = 0; i < nv; i++)
        {
            frx[rows[i]] += vals[i] * fx[cols[i]];
            fry[rows[i]] += vals[i] * fy[cols[i]];
            frz[rows[i]] += vals[i] * fz[cols[i]];
        }

        Un = 0.;
        Rn = 0.;

        #pragma omp parallel for private(i,j,Rx,Ry,Rz,Mi) reduction(+:Un,Rn)
        for (i = 0; i < n; i++)
        {
            j = 3 * i;
            Rx = (P[j + 0] - S[j + 0] - frx[i]) * B[j + 0];
            Ry = (P[j + 1] - S[j + 1] - fry[i]) * B[j + 1];
            Rz = (P[j + 2] - S[j + 2] - frz[i]) * B[j + 2];
            Rn += gsl_hypot3(Rx, Ry, Rz);
            Mi = M[i] * factor;
            V[j + 0] += Rx / Mi;
            V[j + 1] += Ry / Mi;
            V[j + 2] += Rz / Mi;
            Un += Mi * (gsl_pow_2(V[j + 0]) + gsl_pow_2(V[j + 1]) + gsl_pow_2(V[j + 2]));
        }

        if (Un < Uo)
        {
            #pragma omp parallel for private(i,j)
            for (i = 0; i < n; i++)
            {
                j = 3 * i;
                V[j + 0] = 0.;
                V[j + 1] = 0.;
                V[j + 2] = 0.;
            }
        }

        Uo = Un;

        #pragma omp parallel for private(i,j)
        for (i = 0; i < n; i++)
        {
            j = 3 * i;
            X[j + 0] += V[j + 0];
            X[j + 1] += V[j + 1];
            X[j + 2] += V[j + 2];
        }

        res = Rn / n;
        ts++;

    }

    if (summary == 1)
    {
        printf("Step: %i, Residual: %f\n", ts - 1, res);
    }

}

ext: c

max_stars_repo_name: philianeles/compas

max_issues_repo_name: philianeles/compas

max_forks_repo_name: philianeles/compas


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>
#include <math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

void read_matrix(int** index, int** matrix, double scaling, int N_kw, char* input_fileName)
{
    FILE *fp = fopen(input_fileName, "r");
    fscanf(fp, "%*[^\n]\n");

    for (int ii = 0; ii < N_kw; ii++)
        fscanf(fp, "%d,", &((*index)[ii]));
    fscanf(fp, "%*[^\n]\n");

    int tmp;
    for (int ii = 0; ii < N_kw; ii++)
    {
        for (int jj = 0; jj < N_kw; jj++)
        {
            fscanf(fp, "%d,", &tmp);
            (*matrix)[ii*N_kw + jj] = (int) (tmp * scaling);
        }

        fscanf(fp, "%*[^\n]\n");
    }
    fclose(fp);
}

void pad_matrix(int** matrix_padded, int** matrix, int m, double p, double q, int N_kw, int N_doc)
{
    // Initialising RNG
    const gsl_rng_type * T;
    gsl_rng * r;
    gsl_rng_env_setup();
    T = gsl_rng_default;
    r = gsl_rng_alloc(T);

    // perform padding
    int ii, jj;
    #pragma omp parallel for private(ii, jj)
    for (ii = 0; ii < N_kw; ii++)
    {
        for (jj = 0; jj < N_kw; jj++)
        {
            if (ii == jj)
                (*matrix_padded)[ii*N_kw + jj] = gsl_ran_binomial(r, p, m*(*matrix)[ii*N_kw + jj]) + gsl_ran_binomial(r, q, m*(N_doc - (*matrix)[ii*N_kw + jj]));

            if (ii > jj)
            {
                int N1 = (*matrix)[ii*N_kw + jj];
                int N2 = (*matrix)[ii*N_kw + ii] + (*matrix)[jj*N_kw + jj] - 2*(*matrix)[ii*N_kw + jj];
                int N3 = N_doc - N1 - N2;

                int n1 = gsl_ran_binomial(r, p*p, m*N1);
                int n2 = gsl_ran_binomial(r, p*q, m*N2);
                int n3 = gsl_ran_binomial(r, q*q, m*N3);

                (*matrix_padded)[ii*N_kw + jj] = n1 + n2 + n3;
                (*matrix_padded)[jj*N_kw + ii] = n1 + n2 + n3;
            }
        }
    }
    gsl_rng_free(r);
}

void observe_matrix(gsl_matrix* matrix_obs, int** matrix_padded, int N_kw)
{
    // perform observed count generation
    for (int ii = 0; ii < N_kw; ii++)
        for (int jj = 0; jj < N_kw; jj++)
            gsl_matrix_set(matrix_obs, ii, jj, (double) ((*matrix_padded)[ii*N_kw + jj]));
}

void solution_initial(int** permutation, gsl_matrix* matrix_obs, int** matrix, int m, double p, double q, int N_kw, int N_obs, int N_doc)
{
    int* diff = (int*) malloc(sizeof(int) * N_kw);
    int* index = (int*) malloc(sizeof(int) * N_kw);
    int* occupied = (int*) malloc(sizeof(int) * N_kw);

    for (int ii = 0; ii < N_kw; ii++)
        occupied[ii] = -1;

    for (int ii = 0; ii < N_obs; ii++)
    {
        // reset index
        for (int jj = 0; jj < N_kw; jj++)
            index[jj] = jj;

        // compute differences
        for (int jj = 0; jj < N_kw; jj++)
            diff[jj] = abs(gsl_matrix_get(matrix_obs, ii, ii) - m * p * (*matrix)[jj*N_kw+jj] - m * q * (N_doc - (*matrix)[jj*N_kw+jj]));

        // sorting
        for (int jj = 0; jj < N_kw-1; jj++)
        {
            for (int kk = 0; kk < N_kw-jj-1; kk++)
            {
                if (diff[kk] > diff[kk+1])
                {
                    int temp = diff[kk+1];
                    diff[kk+1] = diff[kk];
                    diff[kk] = temp;

                    temp = index[kk+1];
                    index[kk+1] = index[kk];
                    index[kk] = temp;
                }
            }
        }

        // use the smallest index available
        int idx = 0;
        while (occupied[index[idx]] > 0)
            idx++;

        // update the data structures
        (*permutation)[ii] = index[idx];
        occupied[index[idx]] = 1;
    }

    //for (int ii = 200; ii < 250; ii++)
    //{
    //    printf("%f\n", gsl_matrix_get(matrix_obs, ii, ii) - m * p * (*matrix)[(*permutation)[ii]*N_kw+(*permutation)[ii]] - m * q * (N_doc - (*matrix)[(*permutation)[ii]*N_kw+(*permutation)[ii]]));
    //}
}

void permutation_generation(int* idx1, int* idx2, int** permutation_tmp, int** permutation, int** permutation_inv, gsl_matrix* matrix_obs, int** matrix, int m, double p, double q, int N_kw, int N_obs, int N_doc)
{
    double N1_mean, N1_var, N1_lower, N1_upper, N2_mean, N2_var, N2_lower, N2_upper;
    int check = -1;
    int count = 0;
    *idx1 = rand() % N_obs;
    *idx2 = -1;
    int idx_old = (*permutation)[*idx1];
    int idx_new = 0;

    do {
        idx_new = rand() % N_kw;
        count++;

        if ((*permutation_inv)[idx_new] >= 0)
        {
            *idx2 = (*permutation_inv)[idx_new];

            N1_mean = m * p * (*matrix)[idx_new*N_kw + idx_new] + m * q * (N_doc - (*matrix)[idx_new*N_kw + idx_new]);
            N1_var = (*matrix)[idx_new*N_kw + idx_new] / (double) N_doc * (N_doc - (*matrix)[idx_new*N_kw + idx_new]) * 2.0 * m;
            N1_var += m * p * (1-p) * (*matrix)[idx_new*N_kw + idx_new] + m * q * (1-q) * (N_doc - (*matrix)[idx_new*N_kw + idx_new]);
            N1_var =  3*sqrt(N1_var);

            N1_lower = N1_mean - N1_var;
            N1_upper = N1_mean + N1_var;

            N2_mean = m * p * (*matrix)[idx_old*N_kw + idx_old] + m * q * (N_doc - (*matrix)[idx_old*N_kw + idx_old]);
            N2_var = (*matrix)[idx_old*N_kw + idx_old] / (double) N_doc * (N_doc - (*matrix)[idx_old*N_kw + idx_old]) * 2.0 * m ;
            N2_var += m * p * (1-p) * (*matrix)[idx_old*N_kw + idx_old] + m * q * (1-q)* (N_doc - (*matrix)[idx_old*N_kw + idx_old]);
            N2_var =  3*sqrt(N2_var);

            N2_lower = N2_mean - N2_var;
            N2_upper = N2_mean + N2_var;

            if ((gsl_matrix_get(matrix_obs, *idx1, *idx1) > N1_lower) && (gsl_matrix_get(matrix_obs, *idx1, *idx1) < N1_upper))
                if ((gsl_matrix_get(matrix_obs, *idx2, *idx2) > N2_lower) && (gsl_matrix_get(matrix_obs, *idx2, *idx2) < N2_upper))
                    check = 1;
        }
        else
        {
	    *idx2 = (*permutation_inv)[idx_new];

            N1_mean = m * p * (*matrix)[idx_new*N_kw + idx_new] + m * q * (N_doc - (*matrix)[idx_new*N_kw + idx_new]);
            N1_var = (*matrix)[idx_new*N_kw + idx_new] / (double) N_doc * (N_doc - (*matrix)[idx_new*N_kw + idx_new]) * 2.0 * m;
            N1_var += m * p * (1-p) * (*matrix)[idx_new*N_kw + idx_new] + m * q * (1-q) * (N_doc - (*matrix)[idx_new*N_kw + idx_new]);
            N1_var =  3*sqrt(N1_var);

            N1_lower = N1_mean - N1_var;
            N1_upper = N1_mean + N1_var;

            if ((gsl_matrix_get(matrix_obs, *idx1, *idx1) > N1_lower) && (gsl_matrix_get(matrix_obs, *idx1, *idx1) < N1_upper))
                check = 1;

        }

    } while ((check < 0) && (count < 200));

    if (count == 200)
        *idx2 = *idx1;

    if (*idx1 != *idx2)
    {
        (*permutation_tmp)[*idx1] = idx_new;
        if ((*permutation_inv)[idx_new] >= 0)
        {
            *idx2 = (*permutation_inv)[idx_new];
            (*permutation_tmp)[*idx2] = idx_old;
        }
    }
}

max_stars_repo_name: RethinkingSSE/Attacks-on-SSE

max_issues_repo_name: RethinkingSSE/Attacks-on-SSE

max_forks_repo_name: RethinkingSSE/Attacks-on-SSE


text:
#ifndef __QUANTUM_H__
#define __QUANTUM_H__

#include <iostream>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_poly.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_sort_vector.h>

#endif

void WignerTransform_au(gsl_matrix_complex * dest, gsl_vector_complex * source, int grades, double dx, double dp);
void WignerTransform_SI(gsl_matrix_complex * dest,gsl_vector_complex * source, int grades, double dx, double dp);

max_stars_repo_name: Walter-Feng/HomebrewLib

max_issues_repo_name: Walter-Feng/HomebrewLib

max_forks_repo_name: Walter-Feng/HomebrewLib


text:
#include "jfftw_Wisdom.h"
#include <fftw.h>

/*
 * Class:     jfftw_Wisdom
 * Method:    get
 * Signature: ()Ljava/lang/String;
 */
JNIEXPORT jstring JNICALL Java_jfftw_Wisdom_get( JNIEnv * env, jclass clazz )
{
	char* cwisdom = fftw_export_wisdom_to_string();

	jstring wisdom = (*env)->NewStringUTF( env, cwisdom );
	fftw_free( cwisdom );

	return wisdom;
}
/*
 * Class:     jfftw_Wisdom
 * Method:    add
 * Signature: (Ljava/lang/String;)V
 */
JNIEXPORT void JNICALL Java_jfftw_Wisdom_add( JNIEnv * env, jclass clazz, jstring wisdom )
{
	const char* cwisdom = (*env)->GetStringUTFChars( env, wisdom, NULL );

	fftw_status s = fftw_import_wisdom_from_string( cwisdom );

	(*env)->ReleaseStringUTFChars( env, wisdom, cwisdom );

	if( s == FFTW_FAILURE )
	{
		(*env)->ThrowNew( env, (*env)->FindClass( env, "java/lang/IllegalArgumentException" ), "unable to parse wisdom" );
	}
}
/*
 * Class:     jfftw_Wisdom
 * Method:    clear
 * Signature: ()V
 */
JNIEXPORT void JNICALL Java_jfftw_Wisdom_clear( JNIEnv * env, jclass clazz )
{
	fftw_forget_wisdom();
}

max_stars_repo_name: CSCSI/Triana

max_issues_repo_name: CSCSI/Triana

max_forks_repo_name: CSCSI/Triana


text:
/*
  NAME:
     bovy_det
  PURPOSE:
     returns the determinant of a matrix
  CALLING SEQUENCE:
     bovy_det(gsl_matrix * A)
  INPUT:
     A - matrix
  OUTPUT:
     det(A)
  REVISION HISTORY:
     2008-09-21 - Written Bovy
 */
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>

double bovy_det(gsl_matrix * A){
  gsl_permutation * p = gsl_permutation_alloc (A->size1);
  int signum;
  double det;
  gsl_linalg_LU_decomp(A,p,&signum);
  det= gsl_linalg_LU_det(A,signum);
  gsl_permutation_free(p);

  return det;
}

max_stars_repo_name: surbut/mashr

max_issues_repo_name: surbut/mashr

max_forks_repo_name: surbut/mashr


text:
#ifndef STDAFX__H
#define STDAFX__H

#include <SDKDDKVer.h>

#if !defined(_STL_EXTRA_DISABLED_WARNINGS)
#define _STL_EXTRA_DISABLED_WARNINGS 4061 4324 4365 4514 4571 4582 4583 4623 4625 4626 4710 4774 4820 4987 5026 5027 5039
#endif

#if !defined(_SCL_SECURE_NO_WARNINGS)
#define _SCL_SECURE_NO_WARNINGS 1
#endif

#if !defined(_CRT_SECURE_NO_WARNINGS)
#define _CRT_SECURE_NO_WARNINGS 1
#endif

#if !defined(_SILENCE_CXX17_ITERATOR_BASE_CLASS_DEPRECATION_WARNING)
#define _SILENCE_CXX17_ITERATOR_BASE_CLASS_DEPRECATION_WARNING 1
#endif

#define STRICT
#define NOMINMAX

#pragma warning(disable: 4668) // warning C4668: '%s' is not defined as a preprocessor macro, replacing with '0' for '#if/#elif'
#pragma warning(disable: 4710) // warning C4710: '%s': function not inlined
#pragma warning(disable: 4711) // warning C4711: function '%s' selected for automatic inline expansion
#pragma warning(disable: 4820) // warning C4820: '%s': '%d' bytes padding added after data member '%s'

#pragma warning(push)
#pragma warning(disable: 5039) // warning C5039: '%s': pointer or reference to potentially throwing function passed to extern C function under -EHc. Undefined behavior may occur if this function throws an exception.
#include <Windows.h>
#pragma warning(pop)

#pragma warning(push)
#pragma warning(disable: 4005) // warning C4005: '%s': macro redefinition
#include <Winternl.h>
#include <ntstatus.h>
#pragma warning(pop)

// disable for everything
#pragma warning(disable: 4061) // warning C4061: enumerator '%s' in switch of enum '%s' is not explicitly handled by a case label
#pragma warning(disable: 4324) // warning C4234: structure was padded due to alignment specifier
#pragma warning(disable: 4514) // warning C4514: '%s': unreferenced inline function has been removed
#pragma warning(disable: 4623) // warning C4623: '%s': default constructor was implicitly defined as deleted
#pragma warning(disable: 4625) // warning C4625: '%s': copy constructor was implicitly defined as deleted
#pragma warning(disable: 4626) // warning C4626: '%s': assignment operator was implicitly defined as deleted
#pragma warning(disable: 4710) // warning C4710: '%s': function not inlined
#pragma warning(disable: 4820) // warning C4820: '%s': '%d' bytes padding added after data member '%s'
#pragma warning(disable: 5026) // warning C5026: '%s': move constructor was implicitly defined as deleted
#pragma warning(disable: 5027) // warning C5027: '%s': move assignment operator was implicitly defined as deleted

#pragma warning(disable: 26412) // warning C26412: Do not dereference an invalid pointer (lifetimes rule 1). 'return of %s' was invalidated at line %d by 'no initialization'.
#pragma warning(disable: 26426) // warning C26426: Global initializer calls a non-constexpr function '%s' (i.22: http://go.microsoft.com/fwlink/?linkid=853919).
#pragma warning(disable: 26481) // warning C26481: Don't use pointer arithmetic. Use span instead. (bounds.1: http://go.microsoft.com/fwlink/p/?LinkID=620413)
#pragma warning(disable: 26482) // warning C26482: Only index into arrays using constant expressions (bounds.2: http://go.microsoft.com/fwlink/p/?LinkID=620414).
#pragma warning(disable: 26485) // warning C26485: Expression '%s::`vbtable'': No array to pointer decay. (bounds.3: http://go.microsoft.com/fwlink/p/?LinkID=620415)
#pragma warning(disable: 26490) // warning C26490: Don't use reinterpret_cast. (type.1: http://go.microsoft.com/fwlink/p/?LinkID=620417)
#pragma warning(disable: 26499) // warning C26499: Could not find any lifetime tracking information for '%s'

// disable for standard headers
#pragma warning(push)
#pragma warning(disable: 5039) // warning C5039: '%s': pointer or reference to potentially throwing function passed to extern C function under -EHc. Undefined behavior may occur if this function throws an exception.

#pragma warning(disable: 26400) // warning C26400: Do not assign the result of an allocation or a function call with an owner<T> return value to a raw pointer, use owner<T> instead. (i.11 http://go.microsoft.com/fwlink/?linkid=845474)
#pragma warning(disable: 26401) // warning C26401: Do not delete a raw pointer that is not an owner<T>. (i.11: http://go.microsoft.com/fwlink/?linkid=845474)
#pragma warning(disable: 26408) // warning C26408: Avoid malloc() and free(), prefer the nothrow version of new with delete. (r.10 http://go.microsoft.com/fwlink/?linkid=845483)
#pragma warning(disable: 26409) // warning C26409: Avoid calling new and delete explicitly, use std::make_unique<T> instead. (r.11 http://go.microsoft.com/fwlink/?linkid=845485)
#pragma warning(disable: 26411) // warning C26411: The parameter '%s' is a reference to unique pointer and it is never reassigned or reset, use T* or T& instead. (r.33 http://go.microsoft.com/fwlink/?linkid=845479)
#pragma warning(disable: 26412) // warning C26412: Do not dereference an invalid pointer (lifetimes rule 1). 'return of %s' was invalidated at line %d by 'end of function scope (local lifetimes end)'.
#pragma warning(disable: 26413) // warning C26413: Do not dereference nullptr (lifetimes rule 2). 'nullptr' was pointed to nullptr at line %d.
#pragma warning(disable: 26423) // warning C26423: The allocation was not directly assigned to an owner.
#pragma warning(disable: 26424) // warning C26424: Failing to delete or assign ownership of allocation at line %d.
#pragma warning(disable: 26425) // warning C26425: Assigning '%s' to a static variable.
#pragma warning(disable: 26461) // warning C26461: The reference argument '%s' for function %s can be marked as const. (con.3: https://go.microsoft.com/fwlink/p/?LinkID=786684)
#pragma warning(disable: 26471) // warning C26471: Don't use reinterpret_cast. A cast from void* can use static_cast. (type.1: http://go.microsoft.com/fwlink/p/?LinkID=620417).
#pragma warning(disable: 26481) // warning C26481: Don't use pointer arithmetic. Use span instead. (bounds.1: http://go.microsoft.com/fwlink/p/?LinkID=620413)
#pragma warning(disable: 26482) // warning C26482: Only index into arrays using constant expressions. (bounds.2: http://go.microsoft.com/fwlink/p/?LinkID=620414)
#pragma warning(disable: 26490) // warning C26490: Don't use reinterpret_cast. (type.1: http://go.microsoft.com/fwlink/p/?LinkID=620417)
#pragma warning(disable: 26493) // warning C26493: Don't use C-style casts that would perform a static_cast downcast, const_cast, or reinterpret_cast. (type.4: http://go.microsoft.com/fwlink/p/?LinkID=620420)
#pragma warning(disable: 26494) // warning C26494: Variable '%s' is uninitialized. Always initialize an object. (type.5: http://go.microsoft.com/fwlink/p/?LinkID=620421)
#pragma warning(disable: 26495) // warning C26495: Variable '%s' is uninitialized. Always initialize a member variable. (type.6: http://go.microsoft.com/fwlink/p/?LinkID=620422)
#pragma warning(disable: 26496) // warning C26496: Variable '%s' is assigned only once, mark it as const. (con.4: https://go.microsoft.com/fwlink/p/?LinkID=784969)
#pragma warning(disable: 26497) // warning C26497: This function %s could be marked constexpr if compile-time evaluation is desired. (f.4: https://go.microsoft.com/fwlink/p/?LinkID=784970)

#include <cstddef>
#include <map>
#include <set>
#include <vector>
#include <string>
#include <algorithm>
#include <iostream>
#include <iomanip>
#include <type_traits>
#include <utility>
#include <memory>
#include <tuple>
#include <thread>
#include <cstdlib>

// disable additionally for third-party libraries
#pragma warning(push)
#pragma warning(disable: 4365) // warning C4365: '%s': conversion from '%s' to '%s', signed/unsigned mismatch
#pragma warning(disable: 4371) // warning C4371: '%s': layout of class may have changed from a previous version of the compiler due to better packing of member '%s'
#pragma warning(disable: 4619) // warning C4619: #pragma warning: there is no warning number '%d'
#pragma warning(disable: 5031) // warning C5031: #pragma warning(pop): likely mismatch, popping warning state pushed in different file
#pragma warning(disable: 26429) // warning C26429: Symbol '%s' is never tested for nullness, it can be marked as not_null (f.23: http://go.microsoft.com/fwlink/?linkid=853921).
#pragma warning(disable: 26432) // warning C26432: If you define or delete any default operation in the type '%s', define or delete them all (c.21: http://go.microsoft.com/fwlink/?linkid=853922).
#pragma warning(disable: 26434) // warning C26434: Function '%s' hides a non-virtual function '%s' (c.128: http://go.microsoft.com/fwlink/?linkid=853923).
#pragma warning(disable: 26439) // warning C26439: This kind of function may not throw. Declare it 'noexcept' (f.6: http://go.microsoft.com/fwlink/?linkid=853927).
#pragma warning(disable: 26440) // warning C26440: Function '%s' can be declared 'noexcept' (f.6: http://go.microsoft.com/fwlink/?linkid=853927).
#pragma warning(disable: 26472) // warning C26472: Don't use a static_cast for arithmetic conversions. Use brace initialization, gsl::narrow_cast or gsl::narow (type.1: http://go.microsoft.com/fwlink/p/?LinkID=620417).

#ifndef BOOST_CONFIG_SUPPRESS_OUTDATED_MESSAGE
#define BOOST_CONFIG_SUPPRESS_OUTDATED_MESSAGE
#endif

#include <gsl/gsl>

#include <fmt/format.h>

#pragma warning(pop)

#pragma warning(pop)

#endif

max_stars_repo_name: DrPizza/benchmarking-tools

max_issues_repo_name: DrPizza/benchmarking-tools

max_forks_repo_name: DrPizza/benchmarking-tools


text:
#include <stdio.h>
#include <time.h>
#include <linux/limits.h>
#include <sys/stat.h>
#include <unistd.h>
#include <math.h>
#include "Globals.h"
#include "constitutive_equations.h"
#include "Watershed.h"
#include "manufactured_solution.h"
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>

extern void GetLGLWeights(int N, double* w);

void *xcalloc(int items, int size)
{
	void *ptr = calloc(items, size);
	if (ptr == NULL)
	{
		printf("Unable to allocate memory\n");
		exit(EXIT_FAILURE);
	}
	return ptr;
}

static int firstOutput = 1;
static char folderName[500];
static char ChannelHeightFileName[PATH_MAX + 1];
static char ChannelQFileName[PATH_MAX + 1];
static char KinFileName[PATH_MAX + 1];
static char JuncHeightFileName[PATH_MAX + 1];
static char JuncQFileName[PATH_MAX + 1];
static char FpHeightFileName[PATH_MAX + 1];
static char FpQFileName[PATH_MAX + 1];

void outputFloodplainDataToFile(double currtime, double finalTime, double recordTimeIntervals)
{
	char openFileFormat[1];
	int writeHeader = 0;
	if (firstOutput)
	{
		time_t rawtime;
		time(&rawtime);
		struct tm *timeinfo = localtime(&rawtime);
		sprintf(folderName, "DataOutput_%d_%02d_%02d_%02d_%02d_%02d",timeinfo->tm_year+1900, timeinfo->tm_mon+1, timeinfo->tm_mday, timeinfo->tm_hour, timeinfo->tm_min, timeinfo->tm_sec);
		if (mkdir(folderName, 0777) == -1)
			perror("The following error occured while making the data output directory");

		sprintf(FpHeightFileName, "%s/Floodplains.63", folderName);
		sprintf(FpQFileName, "%s/Floodplains.64", folderName);

		openFileFormat[0] = 'w';
		writeHeader = 1;
		firstOutput = 0;
	}
	else
	{
		openFileFormat[0] = 'a';
	}

	// Output Floodplain data
	FILE *file1 = fopen(FpHeightFileName, openFileFormat);
	FILE *file2 = fopen(FpQFileName, openFileFormat);
	if (!file1)
		printf("Error while opening file1\n");
	if (!file2)
		printf("Error while opening file2\n");
	if (writeHeader)
	{
		int numRecords;
		if (fmod(finalTime, recordTimeIntervals) == 0)
			numRecords = finalTime/recordTimeIntervals + 1;
		else
			numRecords = finalTime/recordTimeIntervals + 2;

		fprintf(file1, "End of Simulation Time = %lf.\nTotal Number of Records = %d\n", finalTime, numRecords);
		fprintf(file1, "Number of Floodplains = %d\n", NumFloodplains);
		fprintf(file2, "End of Simulation Time = %lf.\nTotal Number of Records = %d\n", finalTime, numRecords);
		fprintf(file2, "Number of Floodplains = %d\n", NumFloodplains);
		for (int j = 0; j < NumFloodplains; j++)
		{
			fprintf(file1, "Number of data points for junction %d = %d\n", j, FloodplainList[j]->NumEl);
			fprintf(file2, "Number of data points for junction %d = %d\n", j, FloodplainList[j]->NumEl);
		}
	}
	for (int i=0; i < NumFloodplains; i++)
	{
		fprintf(file1, "# Average H and Zeta at time %.3f for Junction %d\n", currtime, i);
		fprintf(file2, "# Average Qx and Qy at time %.3f for Junction %d\n", currtime, i);

		for (int j=0; j < FloodplainList[i]->NumEl; ++j)
		{
			double avgZeta = 0, avgHeight = 0, avgQx = 0, avgQy =0, avgz =0, zarr[3];
			for(int k=0; k<3; ++k)
			{
				double zeta = FloodplainList[i]->zeta[j][k];
				double Qx = FloodplainList[i]->Qx[j][k];
				double Qy = FloodplainList[i]->Qy[j][k];
				zarr[k]= FloodplainList[i]->NodalZ[j][k];
				double height = zeta + zarr[k];

				avgZeta += zeta;
				avgHeight += height;
				avgQx += Qx;
				avgQy += Qy;
				avgz += zarr[k];
			}
			fprintf(file1, "%d \t %.13f \t %.13f \n", j, avgHeight/3, avgZeta/3);
			fprintf(file2, "%d \t %.13f \t %.13f\n", j, avgQx/3, avgQy/3 );
		}
	}
	fclose(file1);
	fclose(file2);

}

void output2DNodalError(double time)
{
	char *fileName1 = "NodalErrorsH.out";
	char *fileName2 = "NodalErrorsQ.out";
	FILE* ef1 = fopen(fileName1, "w");
	FILE* ef2 = fopen(fileName2, "w");

	int i = 0;
	for (int j=0; j < FloodplainList[i]->NumVerts; ++j)
	{
		double xval = FloodplainList[i]->Vx[j];
		double yval = FloodplainList[i]->Vy[j];
		int num_conn_els = FloodplainList[i]->ElCount[i];
		double nodalZeta = 0;
		double nodalQx = 0;
		double nodalQy = 0;
		for (int k = 0; k < num_conn_els; k++)
		{
			int curr_el = FloodplainList[i]->VtoEl[j][k];
			int loc_v;
			// find the local vertex number for curr_el
			for (int v = 0; v < 3; v++)
			{
				if (j == FloodplainList[i]->EltoVert[curr_el*3+v])
				{
					loc_v = v;
					break;
				}
			}

			int loc_node = FloodplainList[i]->VtoNode[curr_el][loc_v];
			nodalZeta += FloodplainList[i]->zeta[curr_el][loc_node];
			nodalQx += FloodplainList[i]->Qx[curr_el][loc_node];
			nodalQy += FloodplainList[i]->Qy[curr_el][loc_node];

		}
		nodalZeta /= num_conn_els;
		nodalQx /= num_conn_els;
		nodalQy /= num_conn_els;

		double manzeta = getmanH(xval, yval, time);
		double manQx = getQx(xval, yval, time);
		double manQy = getQy(xval, yval, time);

		double er1 = -nodalZeta + manzeta;
		double er2 = -nodalQx + manQx;
		double er3 = -nodalQy + manQy;

		fprintf(ef1, "%d \t %.13f\n", j, er1);
		fprintf(ef2, "%d \t %.13f \t %.13f\n", j, er2, er3);
	}

	fclose(ef1);
	fclose(ef2);

}

void outputFloodplainNodalDataToFile(double currtime, double finalTime, double recordTimeIntervals)
{
	char openFileFormat[1];
	int writeHeader = 0;
	if (firstOutput)
	{
		time_t rawtime;
		time(&rawtime);
		struct tm *timeinfo = localtime(&rawtime);
		sprintf(folderName, "DataOutput_%d_%02d_%02d_%02d_%02d_%02d",timeinfo->tm_year+1900, timeinfo->tm_mon+1, timeinfo->tm_mday, timeinfo->tm_hour, timeinfo->tm_min, timeinfo->tm_sec);
		if (mkdir(folderName, 0777) == -1)
			perror("The following error occured while making the data output directory");

		sprintf(FpHeightFileName, "%s/Floodplains.63", folderName);
		sprintf(FpQFileName, "%s/Floodplains.64", folderName);

		openFileFormat[0] = 'w';
		writeHeader = 1;
		firstOutput = 0;
	}
	else
	{
		openFileFormat[0] = 'a';
	}

	// Output Floodplain data
	FILE *file1 = fopen(FpHeightFileName, openFileFormat);
	FILE *file2 = fopen(FpQFileName, openFileFormat);
	if (!file1)
		printf("Error while opening file1\n");
	if (!file2)
		printf("Error while opening file2\n");
	if (writeHeader)
	{
		int numRecords;
		if (fmod(finalTime, recordTimeIntervals) == 0)
			numRecords = finalTime/recordTimeIntervals + 1;
		else
			numRecords = finalTime/recordTimeIntervals + 2;

		fprintf(file1, "End of Simulation Time = %lf.\nTotal Number of Records = %d\n", finalTime, numRecords);
		fprintf(file1, "Number of Floodplains = %d\n", NumFloodplains);
		fprintf(file2, "End of Simulation Time = %lf.\nTotal Number of Records = %d\n", finalTime, numRecords);
		fprintf(file2, "Number of Floodplains = %d\n", NumFloodplains);
		for (int j = 0; j < NumFloodplains; j++)
		{
			fprintf(file1, "Number of data points for junction %d = %d\n", j, FloodplainList[j]->NumVerts);
			fprintf(file2, "Number of data points for junction %d = %d\n", j, FloodplainList[j]->NumVerts);
		}
	}
	for (int i=0; i < NumFloodplains; i++)
	{
		fprintf(file1, "# Nodal H at time %.3f for Junction %d\n", currtime, i);
		fprintf(file2, "# Nodal Qx and Qy at time %.3f for Junction %d\n", currtime, i);

		for (int j=0; j < FloodplainList[i]->NumVerts; ++j)
		{
			int num_conn_els = FloodplainList[i]->ElCount[i];
			double nodalZeta = 0;
			double nodalQx = 0;
			double nodalQy = 0;
			for (int k = 0; k < num_conn_els; k++)
			{
				int curr_el = FloodplainList[i]->VtoEl[j][k];
				int loc_v;
				// find the local vertex number for curr_el
				for (int v = 0; v < 3; v++)
				{
					if (j == FloodplainList[i]->EltoVert[curr_el*3+v])
					{
						loc_v = v;
						break;
					}
				}

				int loc_node = FloodplainList[i]->VtoNode[curr_el][loc_v];
				nodalZeta += FloodplainList[i]->zeta[curr_el][loc_node];
				nodalQx += FloodplainList[i]->Qx[curr_el][loc_node];
				nodalQy += FloodplainList[i]->Qy[curr_el][loc_node];

			}
			nodalZeta /= num_conn_els;
			nodalQx /= num_conn_els;
			nodalQy /= num_conn_els;

			fprintf(file1, "%d \t %.13f\n", j, nodalZeta);
			fprintf(file2, "%d \t %.13f \t %.13f\n", j, nodalQx, nodalQy );
		}
	}
	fclose(file1);
	fclose(file2);

}
void outputDataToFile(double currtime, double finalTime, double recordTimeIntervals)
{
	char openFileFormat[1];
	int writeHeader = 0;
	if (firstOutput)
	{
		time_t rawtime;
		time(&rawtime);
		struct tm *timeinfo = localtime(&rawtime);
		sprintf(folderName, "DataOutput_%d_%02d_%02d_%02d_%02d_%02d",timeinfo->tm_year+1900, timeinfo->tm_mon+1, timeinfo->tm_mday, timeinfo->tm_hour, timeinfo->tm_min, timeinfo->tm_sec);
		if (mkdir(folderName, 0777) == -1)
			perror("The following error occured while making the data output directory");

		sprintf(ChannelHeightFileName, "%s/Channels.63", folderName);
		sprintf(ChannelQFileName, "%s/Channels.64", folderName);
		sprintf(KinFileName, "%s/KinField.63", folderName);
		sprintf(JuncHeightFileName, "%s/Junctions.63", folderName);
		sprintf(JuncQFileName, "%s/Junctions.64", folderName);
		//sprintf(FpHeightFileName, "%s/Floodplains.63", folderName);
		//sprintf(FpQFileName, "%s/Floodplains.64", folderName);

		// Move junction grid data file to this newly created folder
		/*char NewJunctionGridFileName[PATH_MAX + 1];
		sprintf(NewJunctionGridFileName, "%s/JunctionMesh.14", folderName);
		int ret = rename("./Output/JunctionMesh.14", NewJunctionGridFileName);
		if (ret != 0)
			printf("Error: unable to move the junction grid file\n");
*/
		// Move flowpaths file to this newly created folder
		int ret;
		char NewFlowPathsFileName[PATH_MAX + 1];
		sprintf(NewFlowPathsFileName, "%s/FlowPaths.out", folderName);
		ret = rename("./Output/FlowPaths.out", NewFlowPathsFileName);
		if (ret != 0)
			printf("Error: unable to move the flow paths file\n");

		// Move ChannelNodes.14 to this newly created folder
/*		char NewChannelsGridFileName[PATH_MAX + 1];
		sprintf(NewChannelsGridFileName, "%s/Channels.14", folderName);
		ret = rename("./Output/Channels.14", NewChannelsGridFileName);
		if (ret != 0)
			printf("Error: unable to move channels grid file\n");
*/
		openFileFormat[0] = 'w';
		writeHeader = 1;
		firstOutput = 0;
	}
	else
	{
		openFileFormat[0] = 'a';
	}

		//char cwd[PATH_MAX+1];
	//if (getcwd(cwd, sizeof(cwd))==NULL)
	//	perror("Error while getting current working directory");

	//Output Channels data
	FILE* file1;
	FILE* file2;
	file1 = fopen(ChannelHeightFileName, openFileFormat);
	file2 = fopen(ChannelQFileName, openFileFormat);
	//if (!file1 || !file2)
	//{
	//	printf("Could not open file for writing. Exiting now\n");
	//	exit(1);
	//}
	if (writeHeader)
	{
		int numRecords;
		if (fmod(finalTime, recordTimeIntervals) == 0)
			numRecords = finalTime/recordTimeIntervals+1;
		else
			numRecords = finalTime/recordTimeIntervals + 2;
		fprintf(file1, "End of Simulation Time = %lf.\nTotal Number of Records = %d\n", finalTime, numRecords);
		fprintf(file1, "Number of Channels = %d\n", NumChannels);
		fprintf(file2, "End of Simulation Time = %lf.\nTotal Number of Records = %d\n", finalTime, numRecords);
		fprintf(file2, "NumChannels = %d\n", NumChannels);
		for (int c = 0; c < NumChannels; c++)
		{
			fprintf(file1, "Number of data points for channel %d = %d\n", c, ChannelList[c]->NumNodes);
			fprintf(file2, "Number of data points for channel %d = %d\n", c, ChannelList[c]->NumNodes);
		}
	}
	for (int i=0; i<NumChannels; i++)
	{
		fprintf(file1, "# H and Zeta at time %.3f for Channel %d\n", currtime, i);
		fprintf(file2, "#Q at time %.3f for Channel %d\n", currtime, i);

		int NumNodes = ChannelList[i]->NumNodes;

		for (int j = 0; j<NumNodes; j++)
		{
			double bval = ChannelList[i]->NodalB[j];
			double zval = ChannelList[i]->NodalZ[j];
			double x_val = ChannelList[i]->NodalX[j];
			double y_val = ChannelList[i]->NodalY[j];
			double A = ChannelList[i]->A[j+1];
			double Q = ChannelList[i]->Q[j+1];
			double m1val = ChannelList[i]->Nodalm1[j];
			double m2val = ChannelList[i]->Nodalm2[j];
			double H = getH(A, bval, m1val, m2val);
			double zeta = H - zval;
			fprintf(file1, "%.13f \t %.13f\n", H, zeta);
			fprintf(file2, "%.13f\n", Q);

		}

	}

	fclose(file1);
	fclose(file2);

	// Output Junction data
	file1 = fopen(JuncHeightFileName, openFileFormat);
	file2 = fopen(JuncQFileName, openFileFormat);
	if (!file1)
		printf("Error while opening file1\n");
	if (!file2)
		printf("Error while opening file2\n");
	if (writeHeader)
	{
		int numRecords;
		if (fmod(finalTime, recordTimeIntervals) == 0)
			numRecords = finalTime/recordTimeIntervals + 1;
		else
			numRecords = finalTime/recordTimeIntervals + 2;

		fprintf(file1, "End of Simulation Time = %lf.\nTotal Number of Records = %d\n", finalTime, numRecords);
		fprintf(file1, "Number of Junctions = %d\n", NumJunctions);
		fprintf(file2, "End of Simulation Time = %lf.\nTotal Number of Records = %d\n", finalTime, numRecords);
		fprintf(file2, "NumJunctions = %d\n", NumJunctions);
		for (int j = 0; j < NumJunctions; j++)
		{
			fprintf(file1, "Number of data points for junction %d = %d\n", j, JunctionList[j]->NumEl);
			fprintf(file2, "Number of data points for junction %d = %d\n", j, JunctionList[j]->NumEl);
		}
	}
	for (int i=0; i<NumJunctions; i++)
	{

		fprintf(file1, "# Average H and Zeta at time %.3f for Junction %d\n", currtime, i);
		fprintf(file2, "# Average Qx and Qy at time %.3f for Junction %d\n", currtime, i);

		for (int j=0; j<JunctionList[i]->NumEl; ++j)
		{
			double avgZeta = 0, avgHeight = 0, avgQx = 0, avgQy =0, avgz =0, zarr[3];
			for(int k=0; k<3; ++k)
			{
				double zeta = JunctionList[i]->zeta[j][k];
				double Qx = JunctionList[i]->Qx[j][k];
				double Qy = JunctionList[i]->Qy[j][k];
				zarr[k]= JunctionList[i]->NodalZ[j][k];
				double height = zeta + zarr[k];

				avgZeta += zeta;
				avgHeight += height;
				avgQx += Qx;
				avgQy += Qy;
				avgz += zarr[k];

			}
			fprintf(file1, "%d \t %.13f \t %.13f \n", j, avgHeight/3, avgZeta/3);
			fprintf(file2, "%d \t %.13f \t %.13f\n", j, avgQx/3, avgQy/3 );
		}
	}
	fclose(file1);
	fclose(file2);

	// Output Kinematic Field Data
	// Here we assume that we are working with first order polynomials only
	if (FloodplainList[0]->floodedStatus == 0)
	{
		file1 = fopen(KinFileName, openFileFormat);
		if (writeHeader)
		{
			int numRecords;
			if (fmod(finalTime, recordTimeIntervals) == 0)
				numRecords = finalTime/recordTimeIntervals + 1;
			else
				numRecords = finalTime/recordTimeIntervals + 2;

			fprintf(file1, "End of Simulation Time = %lf.\nTotal Number of Records = %d\n", finalTime, numRecords);
		}
		fprintf(file1,"# Average height at time %3.3f in Kinematic Field\n", currtime);

		int NumEl = FloodplainList[0]->NumEl;
		for (int i = 0; i < NumEl; i++)
		{
			if (KinematicElList[i]->isActive == 1)
			{
				double weq = KinematicElList[i]->weq;
				double height = KinematicElList[i]->A[0]/weq;
				if (KinematicElList[i]->numUpstreamEls > 0)
				{
					int numUpstreamEls = KinematicElList[i]->numUpstreamEls;
					height = KinematicElList[i]->A[0]/weq;
					for (int j = 0; j < numUpstreamEls; j++)
					{
						int el = KinematicElList[i]->upstreamEls[j];
						height += KinematicElList[el]->A[1]/weq;
					}
					height = height/(numUpstreamEls+1);
				}
				//printf("%d %.13f\n", i, height);
				fprintf(file1, "%d %.13f\n", i, height);
			}
			else
				fprintf(file1, "%d %1.13f\n", i, 1e-7);
		}
		fclose(file1);
	}

	//static int previouslyFlooded = 0;
	//static char fpOpenFileFormat[1];
	//int writefpHeader = 0;

	//// Output Floodplain data if the channels have flooded
	//if (FloodplainList[0]->floodedStatus == 1)
	//{
	//	if (previouslyFlooded)
	//	{
	//		fpOpenFileFormat[0] = 'a';
	//	}
	//	else
	//	{
	//		fpOpenFileFormat[0] = 'w';
	//		writefpHeader = 1;
	//	}

	//	file1 = fopen(FpHeightFileName, fpOpenFileFormat);
	//	file2 = fopen(FpQFileName, fpOpenFileFormat);

	//
	//	if (!file1)
	//		printf("Error while opening file1\n");
	//	if (!file2)
	//		printf("Error while opening file2\n");
	//	if (writefpHeader)
	//	{
	//		int numRecords;
	//		if (fmod(finalTime, recordTimeIntervals) == 0)
	//			numRecords = finalTime/recordTimeIntervals + 1;
	//		else
	//			numRecords = finalTime/recordTimeIntervals + 2;

	//		fprintf(file1, "End of Simulation Time = %lf.\nTotal Number of Records = %d\n", finalTime, numRecords);
	//		fprintf(file2, "End of Simulation Time = %lf.\nTotal Number of Records = %d\n", finalTime, numRecords);
	//		fprintf(file1, "Number of data points for floodplain = %d\n", FloodplainList[0]->NumEl);
	//		fprintf(file2, "Number of data points for floodplain = %d\n", FloodplainList[0]->NumEl);
	//	}
	//
	//	fprintf(file1, "# Average H and Zeta at time %.3f\n", currtime);
	//	fprintf(file2, "# Average Qx and Qy at time %.3f\n", currtime);
	//
	//	for (int j=0; j<FloodplainList[0]->NumEl; ++j)
	//	{
	//		double avgZeta = 0, avgHeight = 0, avgQx = 0, avgQy =0, avgz =0, zarr[3];
	//		for(int k=0; k<3; ++k)
	//		{
	//			double zeta = FloodplainList[0]->zeta[j][k];
	//			double Qx = FloodplainList[0]->Qx[j][k];
	//			double Qy = FloodplainList[0]->Qy[j][k];
	//			zarr[k]= FloodplainList[0]->NodalZ[j][k];
	//			double height = zeta + zarr[k];
	//
	//			avgZeta += zeta;
	//			avgHeight += height;
	//			avgQx += Qx;
	//			avgQy += Qy;
	//			avgz += zarr[k];

	//		}
	//		fprintf(file1, "%d \t %.13f \t %.13f \n", j, avgHeight/3, avgZeta/3);
	//		fprintf(file2, "%d \t %.13f \t %.13f\n", j, avgQx/3, avgQy/3 );
	//	}
	//	fclose(file1);
	//	fclose(file2);

	//	previouslyFlooded = 1;

	//}

}

void calculate_l2_error_2d(double time)
{
	int NumEl = FloodplainList[0]->NumEl;
	int Np = FloodplainList[0]->Np;
	double errsqH = 0, errsqQx = 0, errsqQy = 0;

	gsl_vector *diffH = gsl_vector_alloc(Np);
	gsl_vector *diffQx = gsl_vector_alloc(Np);
	gsl_vector *diffQy = gsl_vector_alloc(Np);

	//double max = 0;
	for (int i = 0; i < NumEl; i++)
	{
		for (int j = 0; j < Np; j++)
		{
			double xval = FloodplainList[0]->NodalX[i][j];
			double yval = FloodplainList[0]->NodalY[i][j];
			double exactH = getmanH(xval, yval, time);
			double exactQx = getQx(xval, yval, time);
			double exactQy = getQy(xval, yval, time);
			double compH = FloodplainList[0]->zeta[i][j];
			double compQx =  FloodplainList[0]->Qx[i][j];
			double compQy = FloodplainList[0]->Qy[i][j];
			gsl_vector_set(diffH, j, fabs(exactH-compH));
			gsl_vector_set(diffQx, j, fabs(exactQx-compQx));
			gsl_vector_set(diffQy, j, fabs(exactQy-compQy));
			//gsl_vector_set(diffH, j, RHSZeta[i*Np+j]);
			//gsl_vector_set(diffQx, j,  RHSQx[i*Np+j]);
			//gsl_vector_set(diffQy,j, RHSQy[i*Np+j]);
			//max = fmax(fabs(exactQx - compQx), max);

		}

		gsl_vector *tmpH = gsl_vector_alloc(Np);
		gsl_vector *tmpQx = gsl_vector_alloc(Np);
		gsl_vector *tmpQy = gsl_vector_alloc(Np);
		double jac = FloodplainList[0]->jac[i];
		gsl_blas_dgemv(CblasNoTrans, jac, MassMatrix2D, diffH, 0.0, tmpH);
		gsl_blas_dgemv(CblasNoTrans, jac, MassMatrix2D, diffQx, 0.0, tmpQx);
		gsl_blas_dgemv(CblasNoTrans, jac, MassMatrix2D, diffQy, 0.0, tmpQy);

		double elerrHsq, elerrQxsq, elerrQysq;
		gsl_blas_ddot(diffH, tmpH, &elerrHsq);
		gsl_blas_ddot(diffQx, tmpQx, &elerrQxsq);
		gsl_blas_ddot(diffQy, tmpQy, &elerrQysq);

		errsqH += elerrHsq;
		errsqQx += elerrQxsq;
		errsqQy += elerrQysq;

		gsl_vector_free(tmpH);
		gsl_vector_free(tmpQx);
		gsl_vector_free(tmpQy);
	}

	gsl_vector_free(diffH);
	gsl_vector_free(diffQx);
	gsl_vector_free(diffQy);

	printf("time at error computation = %lf\n", time);
	printf("l2 error in H = %lf \n", sqrt(errsqH));
	printf("l2 error in Qx = %lf \n", sqrt(errsqQx));
	printf("l2 error in Qy = %lf \n", sqrt(errsqQy));
	//printf("maximum difference in Qx = %lf\n", max);

}

void append_to_file(char *fileName, double time, double data)
{
	FILE* file = fopen(fileName, "a");
	fprintf(file, "%lf %lf\n", time,data);
	fclose(file);
}

double calculateTotalWater(struct channel *Chan)
{
	int NumEl = Chan->NumEl;
	int Np = Chan->Np;
	int P = Chan->P;

	double totalWater = 0;
	double LGLWeight[Np];
	GetLGLWeights(P, LGLWeight);
	for (int k = 0; k < NumEl; k++)
	{
		double water = 0;
		double dh = Chan->dh[k];
		for (int i = 0; i < Np; i++)
		{
			double m1val = Chan->Nodalm1[k*Np+i];
			double m2val = Chan->Nodalm2[k*Np+i];
			double bval = Chan->NodalB[k*Np+i];
			double A = Chan->A[k*Np+i+1];
			double h = getH(A, bval, m1val, m2val);
			water += LGLWeight[i]*(h-0.0000001);
		}
		totalWater += 0.5*dh*water;
	}

	return totalWater;

}

double ftTom(double ft)
{
	double meter = ft*0.3048;
	return meter;
}

double mToft (double m)
{
	double ft = m*3.28084;
	return ft;

}

max_stars_repo_name: evalseth/DG-RAIN

max_issues_repo_name: evalseth/DG-RAIN

max_forks_repo_name: evalseth/DG-RAIN


text:
#include <Advent/Advent.h>
#include <DayConfig.h>
#include <Input.h>
#include <cblas.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#include <unicode/ustdio.h>
#include <unicode/ustring.h>

enum
{
  OXY = 0,
  CO2,
};

typedef struct
{
  size_t oxy;
  size_t co2;
} Part2Count;

typedef struct
{
  bool highlight;
  size_t highlightWhich;
} PrintBinArrayOptions;

static void PrintHelp(void);
static void PrintVersion(void);
static void PrintBinArray(const char* array, size_t count, PrintBinArrayOptions options);
static uint64_t BinToI64(const char* bin, size_t count);

int main(int argc, char** argv)
{
  CommandLineCommand command;
  if (Advent_ParseCommandLine(argc, argv, &command))
  {
    switch (command)
    {
    case COMMAND_LINE_COMMAND_PRINT_HELP:
      PrintHelp();
      return 0;
    case COMMAND_LINE_COMMAND_PRINT_HELP_ERROR:
      PrintHelp();
      return 1;
    case COMMAND_LINE_COMMAND_PRINT_VERSION:
      PrintVersion();
      return 0;
    default:
      fprintf(stderr, "Unknown command\n");
      return 2;
    }
  }

  UChar* unicodeInput = malloc((inputText_size + 1) * sizeof(UChar));
  int32_t unicodeInputSize = 0;
  UErrorCode errorCode = U_ZERO_ERROR;
  u_strFromUTF8(unicodeInput, (int32_t)(inputText_size + 1), &unicodeInputSize, (const char*)inputText,
                (int32_t)inputText_size, &errorCode);
  if (errorCode != U_ZERO_ERROR)
  {
    fprintf(stderr, "Unicode error: %s\n", u_errorName(errorCode));
    return 3;
  }

  char** matrix = NULL;
  size_t m = 0;
  size_t n = 0;
  size_t cap = 0;
  size_t rowCap = 0;

  size_t width = 0;

  U_STRING_DECL(delim, "\r\n", 2);
  UChar* saveptr;
  for (UChar* iter = u_strtok_r(unicodeInput, delim, &saveptr); iter != NULL; iter = u_strtok_r(NULL, delim, &saveptr))
  {
    if (m == cap)
    {
      cap = cap ? cap * 2 : 1;
      char** newMatrix = realloc(matrix, cap * sizeof(char*));
      matrix = newMatrix;
      memset(&matrix[m], 0, (cap - m) * sizeof(char*));
    }
    char* row = malloc(rowCap * sizeof(char));
    matrix[m] = row;
    size_t i = 0;
    while (true)
    {
      UChar32 ci;
      U16_NEXT(iter, i, -1, ci);
      if (ci == 0)
      {
        break;
      }
      if (ci == '0' || ci == '1')
      {
        if (i > n)
        {
          if (n == rowCap)
          {
            rowCap = rowCap ? rowCap * 2 : 1;
            for (size_t i = 0; i < m; i++)
            {
              char* row = realloc(matrix[i], rowCap * sizeof(char));
              matrix[i] = row;
            }
          }
          ++n;
        }
        matrix[m][i - 1] = (char)ci;
      }
    }
    ++m;
  }

  char** transposed = malloc(n * sizeof(char*));
  for (size_t i = 0; i < n; ++i)
  {
    transposed[i] = malloc(m * sizeof(char));
    for (size_t j = 0; j < m; ++j)
    {
      transposed[i][j] = matrix[j][i];
    }
  }

  // calculate!
  size_t** counts = malloc(n * sizeof(size_t*));
  uint64_t gamma = 0;
  uint64_t epsilon = 0;
  for (size_t i = 0; i < n; ++i)
  {
    counts[i] = calloc(2, sizeof(size_t));
    for (size_t j = 0; j < m; ++j)
    {
      // NIT: This assumes that the matrix is a rectangle!
      // If there are ANY missing elements, this will likely segfault due to
      // out-of-bounds array access.
      counts[i][transposed[i][j] - '0']++;
    }
    gamma <<= 1;
    gamma |= counts[i][1] > counts[i][0];
    epsilon <<= 1;
    epsilon |= counts[i][0] > counts[i][1];
  }

  printf("%lu\n", gamma * epsilon);

  Part2Count** p2counts = malloc(n * sizeof(Part2Count*));
  for (size_t i = 0; i < n; ++i)
  {
    p2counts[i] = calloc(2, sizeof(Part2Count));
  }
  bool** keep = malloc(m * sizeof(bool*));
  for (size_t i = 0; i < m; ++i)
  {
    keep[i] = malloc(2 * sizeof(bool));
    keep[i][OXY] = true;
    keep[i][CO2] = true;
  }
  size_t oxy = 0;
  size_t co2 = 0;
  for (size_t i = 0; i < n; ++i)
  {
    memset(p2counts[i], 0, 2 * sizeof(Part2Count));
    for (size_t j = 0; j < m; ++j)
    {
      if (keep[j][OXY])
      {
        p2counts[i][transposed[i][j] - '0'].oxy++;
      }
      if (keep[j][CO2])
      {
        p2counts[i][transposed[i][j] - '0'].co2++;
      }
    }
    bool selector = p2counts[i][0].oxy > p2counts[i][1].oxy;
    // filter out 1's for OXY, 0's for CO2
    size_t keptOxy = SIZE_MAX;
    size_t keptCo2 = SIZE_MAX;
    for (size_t j = 0; j < m; ++j)
    {
      if (keep[j][OXY])
      {
        if (transposed[i][j] == (p2counts[i][0].oxy > p2counts[i][1].oxy ? '1' : '0'))
        {
          keep[j][OXY] = false;
        }
        else if (keptOxy == SIZE_MAX)
        {
          keptOxy = j;
        }
      }
      if (keep[j][CO2])
      {
        if (transposed[i][j] == (p2counts[i][0].co2 > p2counts[i][1].co2 ? '0' : '1'))
        {
          keep[j][CO2] = false;
        }
        else if (keptCo2 == SIZE_MAX)
        {
          keptCo2 = j;
        }
      }
    }
    if (keptOxy != SIZE_MAX)
    {
      oxy = keptOxy;
    }
    if (keptCo2 != SIZE_MAX)
    {
      co2 = keptCo2;
    }
  }

  printf("Filtering complete\n");
  printf("OXY: index %lu\n", oxy);
  printf("OXY: ");
  PrintBinArray(matrix[oxy], n, (PrintBinArrayOptions){0});
  printf(" (as array)\n");
  printf("OXY: %lu (as decimal)\n", BinToI64(matrix[oxy], n));
  printf("CO2: index %lu\n", co2);
  printf("CO2: ");
  PrintBinArray(matrix[co2], n, (PrintBinArrayOptions){0});
  printf(" (as array)\n");
  printf("CO2: %lu (as decimal)\n", BinToI64(matrix[co2], n));
  printf("Part 2: %lu\n", BinToI64(matrix[oxy], n) * BinToI64(matrix[co2], n));

  for (size_t i = 0; i < m; ++i)
  {
    free(keep[i]);
  }
  free(keep);
  for (size_t i = 0; i < n; ++i)
  {
    free(p2counts[i]);
  }
  free(p2counts);

  for (size_t i = 0; i < m; ++i)
  {
    free(matrix[i]);
  }
  free(matrix);
  for (size_t i = 0; i < n; ++i)
  {
    free(transposed[i]);
    free(counts[i]);
  }
  free(transposed);
  free(counts);
  free(unicodeInput);
}

static void PrintHelp(void)
{
  PrintVersion();
  printf("Usage: %s [options]\n", ADVENT_DAY_NAME);
  puts("Options:");
  puts(" -h,--help    Print this help message");
  puts(" -V,--version Print the program version");
}

static void PrintVersion(void)
{
  printf("%s v%s\n", ADVENT_DAY_NAME, PROJECT_VERSION);
}

static void PrintBinArray(const char* array, size_t count, PrintBinArrayOptions options)
{
  for (size_t i = 0; i < count; ++i)
  {
    if (options.highlight && options.highlightWhich == i)
    {
      putchar('{');
    }
    putchar(array[i]);
    if (options.highlight && options.highlightWhich == i)
    {
      putchar('}');
    }
  }
}

static uint64_t BinToI64(const char* bin, size_t count)
{
  char* str = malloc(count + 1);
  memcpy(str, bin, count);
  str[count] = '\0';
  uint64_t result = strtoul(str, NULL, 2);
  free(str);
  return result;
}