text:
/*
 *  aeif_cond_beta_multisynapse.h
 *
 *  This file is part of NEST.
 *
 *  Copyright (C) 2004 The NEST Initiative
 *
 *  NEST 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 2 of the License, or
 *  (at your option) any later version.
 *
 *  NEST 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 NEST.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

#ifndef AEIF_COND_BETA_MULTISYNAPSE_H
#define AEIF_COND_BETA_MULTISYNAPSE_H

// Generated includes:
#include "config.h"

#ifdef HAVE_GSL

// External includes:
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>

// Includes from nestkernel:
#include "archiving_node.h"
#include "connection.h"
#include "event.h"
#include "nest_types.h"
#include "ring_buffer.h"
#include "universal_data_logger.h"

/* BeginDocumentation
 Name: aeif_cond_beta_multisynapse - Conductance based adaptive exponential
                                      integrate-and-fire neuron model according
                                      to Brette and Gerstner (2005) with
                                      multiple synaptic rise time and decay
                                      time constants, and synaptic conductance
                                      modeled by a beta function.

 Description:

 aeif_cond_beta_multisynapse is a conductance-based adaptive exponential
 integrate-and-fire neuron model. It allows an arbitrary number of synaptic
 rise time and decay time constants. Synaptic conductance is modeled by a
 beta function, as described by A. Roth and M.C.W. van Rossum
 in Computational Modeling Methods for Neuroscientists, MIT Press 2013,
 Chapter 6.

 The time constants are supplied by two arrays, "tau_rise" and "tau_decay" for
 the synaptic rise time and decay time, respectively. The synaptic
 reversal potentials are supplied by the array "E_rev". The port numbers
 are automatically assigned in the range from 1 to n_receptors.
 During connection, the ports are selected with the property "receptor_type".

 The membrane potential is given by the following differential equation:
 C dV/dt = -g_L(V-E_L) + g_L*Delta_T*exp((V-V_T)/Delta_T) + I_syn_tot(V, t)
           - w + I_e

 where:
 I_syn_tot(V,t) = \sum_i g_i(t) (V - E_{rev,i}) ,

 the synapse i is excitatory or inhibitory depending on the value of E_{rev,i}
 and the differential equation for the spike-adaptation current w is:

 tau_w * dw/dt = a(V - E_L) - w

 When the neuron fires a spike, the adaptation current w <- w + b.

Parameters:
The following parameters can be set in the status dictionary.

Dynamic state variables:
  V_m        double - Membrane potential in mV
  w          double - Spike-adaptation current in pA.

Membrane Parameters:
  C_m        double - Capacity of the membrane in pF
  t_ref      double - Duration of refractory period in ms.
  V_reset    double - Reset value for V_m after a spike. In mV.
  E_L        double - Leak reversal potential in mV.
  g_L        double - Leak conductance in nS.
  I_e        double - Constant external input current in pA.
  Delta_T    double - Slope factor in mV
  V_th       double - Spike initiation threshold in mV
  V_peak     double - Spike detection threshold in mV.

Adaptation parameters:
  a          double - Subthreshold adaptation in nS.
  b          double - Spike-triggered adaptation in pA.
  tau_w      double - Adaptation time constant in ms

Synaptic parameters
  E_rev      double vector - Reversal potential in mV.
  tau_rise  double vector - Rise time of synaptic conductance in ms (beta
                      function).
  tau_decay double vector - Decay time of synaptic conductance in ms (beta
                      function).

Integration parameters
  gsl_error_tol  double - This parameter controls the admissible error of the
                          GSL integrator. Reduce it if NEST complains about
                          numerical instabilities.

 Examples:

 import nest
 import numpy as np

 neuron = nest.Create('aeif_cond_beta_multisynapse')
 nest.SetStatus(neuron, {"V_peak": 0.0, "a": 4.0, "b":80.5})
 nest.SetStatus(neuron, {'E_rev':[0.0,0.0,0.0,-85.0],
                         'tau_decay':[50.0,20.0,20.0,20.0],
                         'tau_rise':[10.0,10.0,1.0,1.0]})

 spike = nest.Create('spike_generator', params = {'spike_times':
                                                 np.array([10.0])})

 voltmeter = nest.Create('voltmeter', 1, {'withgid': True})

 delays=[1.0, 300.0, 500.0, 700.0]
 w=[1.0, 1.0, 1.0, 1.0]
 for syn in range(4):
     nest.Connect(spike, neuron, syn_spec={'model': 'static_synapse',
                                           'receptor_type': 1 + syn,
                                           'weight': w[syn],
                                           'delay': delays[syn]})

 nest.Connect(voltmeter, neuron)

 nest.Simulate(1000.0)
 dmm = nest.GetStatus(voltmeter)[0]
 Vms = dmm["events"]["V_m"]
 ts = dmm["events"]["times"]
 import pylab
 pylab.figure(2)
 pylab.plot(ts, Vms)
 pylab.show()

 Sends: SpikeEvent

 Receives: SpikeEvent, CurrentEvent, DataLoggingRequest

 author: Bruno Golosio 07/10/2016
 SeeAlso: aeif_cond_alpha_multisynapse
 */

namespace nest
{
/**
 * Function computing right-hand side of ODE for GSL solver.
 * @note Must be declared here so we can befriend it in class.
 * @note Must have C-linkage for passing to GSL. Internally, it is
 *       a first-class C++ function, but cannot be a member function
 *       because of the C-linkage.
 * @note No point in declaring it inline, since it is called
 *       through a function pointer.
 * @param void* Pointer to model neuron instance.
 */
extern "C" int
aeif_cond_beta_multisynapse_dynamics( double, const double*, double*, void* );

/**
 * Conductance based exponential integrate-and-fire neuron model according to
 * Brette and Gerstner
 * (2005) with multiple ports.
 */
class aeif_cond_beta_multisynapse : public Archiving_Node
{

public:
  aeif_cond_beta_multisynapse();
  aeif_cond_beta_multisynapse( const aeif_cond_beta_multisynapse& );
  virtual ~aeif_cond_beta_multisynapse();

  friend int
  aeif_cond_beta_multisynapse_dynamics( double, const double*, double*, void* );

  /**
   * Import sets of overloaded virtual functions.
   * @see Technical Issues / Virtual Functions: Overriding, Overloading, and
   * Hiding
   */
  using Node::handle;
  using Node::handles_test_event;

  port send_test_event( Node&, rport, synindex, bool );

  void handle( SpikeEvent& );
  void handle( CurrentEvent& );
  void handle( DataLoggingRequest& );

  port handles_test_event( SpikeEvent&, rport );
  port handles_test_event( CurrentEvent&, rport );
  port handles_test_event( DataLoggingRequest&, rport );

  void get_status( DictionaryDatum& ) const;
  void set_status( const DictionaryDatum& );

private:
  void init_state_( const Node& proto );
  void init_buffers_();
  void calibrate();
  void update( Time const&, const long, const long );

  // The next two classes need to be friends to access the State_ class/member
  friend class RecordablesMap< aeif_cond_beta_multisynapse >;
  friend class UniversalDataLogger< aeif_cond_beta_multisynapse >;

  // ----------------------------------------------------------------

  /**
   * Independent parameters of the model.
   */
  struct Parameters_
  {
    double V_peak_;  //!< Spike detection threshold in mV
    double V_reset_; //!< Reset Potential in mV
    double t_ref_;   //!< Refractory period in ms

    double g_L;     //!< Leak Conductance in nS
    double C_m;     //!< Membrane Capacitance in pF
    double E_L;     //!< Leak reversal Potential (aka resting potential) in mV
    double Delta_T; //!< Slope faktor in ms.
    double tau_w;   //!< adaptation time-constant in ms.
    double a;       //!< Subthreshold adaptation in nS.
    double b;       //!< Spike-triggered adaptation in pA
    double V_th;    //!< Spike threshold in mV.

    std::vector< double > tau_rise;  //!< Rise time of synaptic conductance
                                     //!< in ms.
    std::vector< double > tau_decay; //!< Decay time of synaptic conductance
                                     //!< in ms.
    std::vector< double > E_rev;     //!< reversal potentials in mV

    double I_e; //!< Intrinsic current in pA.

    double gsl_error_tol; //!< error bound for GSL integrator

    // boolean flag which indicates whether the neuron has connections
    bool has_connections_;

    Parameters_(); //!< Sets default parameter values

    void get( DictionaryDatum& ) const; //!< Store current values in dictionary
    void set( const DictionaryDatum& ); //!< Set values from dictionary

    //! Return the number of receptor ports
    inline size_t
    n_receptors() const
    {
      return E_rev.size();
    }
  };

  // ----------------------------------------------------------------

  /**
   * State variables of the model.
   * @note Copy constructor and assignment operator required because
   *       of C-style arrays.
   */
  struct State_
  {

    /**
     * Enumeration identifying elements in state vector State_::y_.
     * This enum identifies the elements of the vector. It must be public to be
     * accessible from the iteration function. The last two elements of this
     * enum (DG, G) will be repeated
     * n times at the end of the state vector State_::y with n being the number
     * of synapses.
     */
    enum StateVecElems
    {
      V_M = 0,
      W,  // 1
      DG, // 2
      G,  // 3
      STATE_VECTOR_MIN_SIZE
    };

    static const size_t NUMBER_OF_FIXED_STATES_ELEMENTS = 2;        // V_M, W
    static const size_t NUMBER_OF_STATES_ELEMENTS_PER_RECEPTOR = 2; // DG, G

    std::vector< double > y_; //!< neuron state
    int r_;                   //!< number of refractory steps remaining

    State_( const Parameters_& ); //!< Default initialization
    State_( const State_& );
    State_& operator=( const State_& );

    void get( DictionaryDatum& ) const;
    void set( const DictionaryDatum& );

  }; // State_

  // ----------------------------------------------------------------

  /**
   * Buffers of the model.
   */
  struct Buffers_
  {
    Buffers_( aeif_cond_beta_multisynapse& );
    Buffers_( const Buffers_&, aeif_cond_beta_multisynapse& );

    //! Logger for all analog data
    UniversalDataLogger< aeif_cond_beta_multisynapse > logger_;

    /** buffers and sums up incoming spikes/currents */
    std::vector< RingBuffer > spikes_;
    RingBuffer currents_;

    /** GSL ODE stuff */
    gsl_odeiv_step* s_;    //!< stepping function
    gsl_odeiv_control* c_; //!< adaptive stepsize control function
    gsl_odeiv_evolve* e_;  //!< evolution function
    gsl_odeiv_system sys_; //!< struct describing system

    // IntergrationStep_ should be reset with the neuron on ResetNetwork,
    // but remain unchanged during calibration. Since it is initialized with
    // step_, and the resolution cannot change after nodes have been created,
    // it is safe to place both here.
    double step_;            //!< simulation step size in ms
    double IntegrationStep_; //!< current integration time step,
                             //!< updated by solver

    /**
     * Input current injected by CurrentEvent.
     * This variable is used to transport the current applied into the
     * _dynamics function computing the derivative of the state vector.
     * It must be a part of Buffers_, since it is initialized once before
     * the first simulation, but not modified before later Simulate calls.
     */
    double I_stim_;
  };

  // ----------------------------------------------------------------

  /**
   * Internal variables of the model.
   */
  struct Variables_
  {

    /** initial value to normalise synaptic conductance */
    std::vector< double > g0_;

    /**
     * Threshold detection for spike events: P.V_peak if Delta_T > 0.,
     * P.V_th if Delta_T == 0.
     */
    double V_peak;

    unsigned int refractory_counts_;
  };

  // Access functions for UniversalDataLogger -------------------------------

  //! Read out state vector elements, used by UniversalDataLogger
  template < State_::StateVecElems elem >
  double
  get_y_elem_() const
  {
    return S_.y_[ elem ];
  }

  // Data members -----------------------------------------------------------

  /**
   * @defgroup aeif_cond_beta_multisynapse
   * Instances of private data structures for the different types
   * of data pertaining to the model.
   * @note The order of definitions is important for speed.
   * @{
   */
  Parameters_ P_;
  State_ S_;
  Variables_ V_;
  Buffers_ B_;
  /** @} */

  //! Mapping of recordables names to access functions
  static RecordablesMap< aeif_cond_beta_multisynapse > recordablesMap_;
};

inline port
aeif_cond_beta_multisynapse::send_test_event( Node& target,
  rport receptor_type,
  synindex,
  bool )
{
  SpikeEvent e;
  e.set_sender( *this );

  return target.handles_test_event( e, receptor_type );
}

inline port
aeif_cond_beta_multisynapse::handles_test_event( CurrentEvent&,
  rport receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
  return 0;
}

inline port
aeif_cond_beta_multisynapse::handles_test_event( DataLoggingRequest& dlr,
  rport receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
  return B_.logger_.connect_logging_device( dlr, recordablesMap_ );
}

inline void
aeif_cond_beta_multisynapse::get_status( DictionaryDatum& d ) const
{
  P_.get( d );
  S_.get( d );
  Archiving_Node::get_status( d );

  ( *d )[ names::recordables ] = recordablesMap_.get_list();
}

inline void
aeif_cond_beta_multisynapse::set_status( const DictionaryDatum& d )
{
  Parameters_ ptmp = P_; // temporary copy in case of errors
  ptmp.set( d );         // throws if BadProperty
  State_ stmp = S_;      // temporary copy in case of errors
  stmp.set( d );         // throws if BadProperty

  // We now know that (ptmp, stmp) are consistent. We do not
  // write them back to (P_, S_) before we are also sure that
  // the properties to be set in the parent class are internally
  // consistent.
  Archiving_Node::set_status( d );

  // if we get here, temporaries contain consistent set of properties
  P_ = ptmp;
  S_ = stmp;
}

} // namespace

#endif // HAVE_GSL
#endif // AEIF_COND_BETA_MULTISYNAPSE_H //

max_stars_repo_name: OpenHEC/SNN-simulator-on-PYNQcluster

max_issues_repo_name: zlchai/SNN-simulator-on-PYNQcluster

max_forks_repo_name: OpenHEC/SNN-simulator-on-PYNQcluster


text:
/*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
**                                                                                **
** This file forms part of the Underworld geophysics modelling application.         **
**                                                                                **
** For full license and copyright information, please refer to the LICENSE.md file  **
** located at the project root, or contact the authors.                             **
**                                                                                **
**~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*/
#include <petsc.h>
#include <petscvec.h>
#include <petscmat.h>

#include "common-driver-utils.h"

// Calculate {r} = {F} - [K]{u} - [G]{p}
// Calculate || {r} ||_2
double BSSCR_StokesMomentumResidual( Mat K, Mat G, Vec F, Vec u, Vec p )
{
	Vec                     uStar;
	//      PetscReal       negOne = -1.0;
	PetscReal               f_norm;
	PetscReal               r1_norm;

	VecNorm( F, NORM_2, &f_norm );          // u_norm = || {uStar} ||_2

	VecDuplicate( u, &uStar );
	MatMult( K, u, uStar );                                 // {uStar} = [K]{u}
	MatMultAdd( G, p, uStar, uStar );               // {uStar} = {uStar} + [G]{p}
	VecAYPX( uStar, -1.0, F );                      // {uStar} = {F} - {uStar}
	VecNorm( uStar, NORM_2, &r1_norm );             // r_norm = || {uStar} ||_2
	Stg_VecDestroy(&uStar );
	/*
	printf("%s \n", __func__ );
	printf("\t||f - Ku - Gp||       = %g \n", r1_norm );
	printf("\t||f - Ku - Gp||/||f|| = %g \n", r1_norm/f_norm );
	*/

	//      return ( (double)(r1_norm/f_norm) );
	return ( (double)(r1_norm) );
}

// Calculate {r} = {H} - [G]^T{u} - [C]{p}
// Calculate || {r} ||_2
double BSSCR_StokesContinuityResidual( Mat G, Mat C, Vec H, Vec u, Vec p )
{
	Vec                     pStar;
	//      PetscReal       negOne = -1.0;
	PetscReal               r2_norm;
	PetscReal               u_norm;

	VecNorm( u, NORM_2, &u_norm );          // u_norm = || {uStar} ||_2

	VecDuplicate( H, &pStar );
	MatMultTranspose( G, u, pStar );                // {pStar} = [G]^T{u}
	if( C != PETSC_NULL ) {
		MatMultAdd( C, p, pStar, pStar );	/* {pStar} = {pStar} + [C] {p} */
	}
	VecAYPX( pStar, -1.0, H );                              // {pStar} = {H} - {pStar}

	VecNorm( pStar, NORM_2, &r2_norm );             // norm = || {pStar} ||_2
	Stg_VecDestroy(&pStar );

	/*
	printf("%s \n", __func__ );
	printf("\t||h - Du||       = %g \n", r2_norm );
	printf("\t||h - Du||/||u|| = %g \n", r2_norm/u_norm );
	*/

	//      return ( (double)(r2_norm/u_norm) );
	return ( (double)(r2_norm) );
}

ext: c


text:
#pragma once

#define EXPORT __declspec(dllexport)

#define WIN32_LEAN_AND_MEAN

// STL
#include <array>
#include <cassert>
#include <cstdint>
#include <deque>
#include <exception>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <stack>
#include <stdexcept>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>

// Windows
#include <Windows.h>
#include <VersionHelpers.h>
#include <WTypes.h>

// DirectX 11
#include <d3d11_1.h>
#include "d3d8types.h"
#include <d3dcompiler.h>
#include <DirectXMath.h>

// GSL
#include <gsl/span>

// DirectXTK
#include <wrl/client.h>
#include <SimpleMath.h>

// Local
#include "cbuffers.h"
#include "CBufferWriter.h"
#include "d3d8to11.hpp"
#include "d3d8to11_base.h"
#include "d3d8to11_device.h"
#include "d3d8to11_index_buffer.h"
#include "d3d8to11_resource.h"
#include "d3d8to11_surface.h"
#include "d3d8to11_texture.h"
#include "d3d8to11_vertex_buffer.h"
#include "d3d8types.hpp"
#include "defs.h"
#include <dirty_t.h>
#include "hash_combine.h"
#include "int_multiple.h"
#include "Light.h"
#include "Material.h"
#include "Shader.h"
#include "simple_math.h"
#include "SimpleMath.h"
#include "typedefs.h"
#include "Unknown.h"
#include "dynarray.h"

ext: h


text:
#pragma once

#include "common.h"

#include <string>
#include <vector>
#include <array>
#include <map>
#include <gsl\gsl>
#include <type_traits>

namespace multiformats {

    enum base_t {
        dynamic_base = -1,
        identity = 0,
        base256,
        //base1,        // '1'
        base2,        // '01'
        base8,        // '01234567'
        base10,       // '0123456789'
        base16,       // '0123456789abcdef'
        BASE16,       // '0123456789ABCDEF'
        base32,       // 'abcdefghijklmnopqrstuvwxyz234567' - rfc4648 no padding
        BASE32,       // 'ABCDEFGHIJKLMNOPQRSTUVWXYZ234567' - rfc4648 no padding
        base32pad,    // 'abcdefghijklmnopqrstuvwxyz234567=' - rfc4648 with padding
        BASE32pad,    // 'ABCDEFGHIJKLMNOPQRSTUVWXYZ234567=' - rfc4648 with padding
        base32hex,    // '0123456789abcdefghijklmnopqrstuv' - rfc4648 no padding - highest char
        BASE32hex,    // '0123456789ABCDEFGHIJKLMNOPQRSTUV' - rfc4648 no padding - highest char
        base32hexpad, // '0123456789abcdefghijklmnopqrstuv=' - rfc4648 with padding
        BASE32hexpad, // '0123456789ABCDEFGHIJKLMNOPQRSTUV=' - rfc4648 with padding
        base32z,      // 'ybndrfg8ejkmcpqxot1uwisza345h769' - z - base - 32 - used by Tahoe - LAFS - highest letter
        base58flickr, // '123456789abcdefghijkmnopqrstuvwxyzABCDEFGHJKLMNPQRSTUVWXYZ' - highest letter
        base58btc,    // '123456789ABCDEFGHJKLMNPQRSTUVWXYZabcdefghijkmnopqrstuvwxyz' - highest letter
        base64,       // 'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/' - rfc4648 no padding
        base64pad,    // 'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/=' - rfc4648 with padding - MIME encoding
        base64url,    // 'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789-_' - rfc4648 no padding
        base64urlpad, // 'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789-_=' - rfc4648 with padding
    };

    namespace details {

        typedef buffer_t(*CodecFunc)(bufferview_t, const char*);
        inline buffer_t codec_noimpl(bufferview_t /*data*/, const char* /*digits*/) { return {}; }
        buffer_t encode_base0(bufferview_t data, const char* digits);
        buffer_t encode_base2(bufferview_t data, const char* digits);
        buffer_t encode_base8(bufferview_t data, const char* digits);
        //buffer_t encode_base10(bufferview_t data, const char* digits);
        buffer_t encode_base16(bufferview_t data, const char* digits);
        buffer_t encode_base32(bufferview_t data, const char* digits);
        //buffer_t encode_base58(bufferview_t data, const char* digits);
        buffer_t encode_base64(bufferview_t data, const char* digits);

        buffer_t decode_base0(bufferview_t data, const char* digits);
        buffer_t decode_base2(bufferview_t data, const char* digits);
        buffer_t decode_base8(bufferview_t data, const char* digits);
        //buffer_t decode_base10(bufferview_t data, const char* digits);
        buffer_t decode_base16(bufferview_t data, const char* digits);
        buffer_t decode_base32(bufferview_t data, const char* digits);
        //buffer_t decode_base58(bufferview_t data, const char* digits);
        buffer_t decode_base64(bufferview_t data, const char* digits);

        template <base_t _FromBase, base_t _ToBase>
        buffer_t convert_base(bufferview_t from, const char*);

        struct baseimpl {
            base_t    key;
            const char* name;
            char        code;
            int         radix;
            CodecFunc   encode;
            CodecFunc   decode;
            const char* digits;
        };

        constexpr baseimpl _BaseTable[] = {
            { dynamic_base, "dynamic_base",       -1,  -1, codec_noimpl                       , codec_noimpl                       , nullptr },
            { base256,      "base256",        0, 256, codec_noimpl                       , codec_noimpl                       , nullptr },
            { identity,     "identity",       0,   0, encode_base0                       , decode_base0                       , nullptr },
            { base2,        "base2",        '0',   2, encode_base2                       , decode_base2                       , "01" },
            { base8,        "base8",        '7',   8, encode_base8                       , decode_base8                       , "01234567" },
            { base10,       "base10",       '9',  10, convert_base<base256, base10>      , convert_base<base10, base256>      , "0123456789" },
            { base16,       "base16",       'f',  16, encode_base16                      , decode_base16                      , "0123456789abcdef" },
            { BASE16,       "BASE16",       'F',  16, encode_base16                      , decode_base16                      , "0123456789ABCDEF" },
            { base32,       "base32",       'b',  32, encode_base32                      , decode_base32                      , "abcdefghijklmnopqrstuvwxyz234567" },
            { BASE32,       "BASE32",       'B',  32, encode_base32                      , decode_base32                      , "ABCDEFGHIJKLMNOPQRSTUVWXYZ234567" },
            { base32pad,    "base32pad",    'c',  32, encode_base32                      , decode_base32                      , "abcdefghijklmnopqrstuvwxyz234567=" },
            { BASE32pad,    "BASE32pad",    'C',  32, encode_base32                      , decode_base32                      , "ABCDEFGHIJKLMNOPQRSTUVWXYZ234567=" },
            { base32hex,    "base32hex",    'v',  32, encode_base32                      , decode_base32                      , "0123456789abcdefghijklmnopqrstuv" },
            { BASE32hex,    "BASE32hex",    'V',  32, encode_base32                      , decode_base32                      , "0123456789ABCDEFGHIJKLMNOPQRSTUV" },
            { base32hexpad, "base32hexpad", 't',  32, encode_base32                      , decode_base32                      , "0123456789abcdefghijklmnopqrstuv=" },
            { BASE32hexpad, "BASE32hexpad", 'T',  32, encode_base32                      , decode_base32                      , "0123456789ABCDEFGHIJKLMNOPQRSTUV=" },
            { base32z,      "base32z",      'h',  32, encode_base32                      , decode_base32                      , "ybndrfg8ejkmcpqxot1uwisza345h769" },
            { base58flickr, "base58flickr", 'Z',  58, convert_base<base256, base58flickr>, convert_base<base58flickr, base256>, "123456789abcdefghijkmnopqrstuvwxyzABCDEFGHJKLMNPQRSTUVWXYZ" },
            { base58btc,    "base58btc",    'z',  58, convert_base<base256, base58btc>   , convert_base<base58btc, base256>   , "123456789ABCDEFGHJKLMNPQRSTUVWXYZabcdefghijkmnopqrstuvwxyz" },
            { base64,       "base64",       'm',  64, encode_base64                      , decode_base64                      , "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/" },
            { base64pad,    "base64pad",    'M',  64, encode_base64                      , decode_base64                      , "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/=" },
            { base64url,    "base64url",    'u',  64, encode_base64                      , decode_base64                      , "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789-_" },
            { base64urlpad, "base64urlpad", 'U',  64, encode_base64                      , decode_base64                      , "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789-_=" },
        };

        constexpr int find_baseimpl(base_t code) {
            for (auto i = 0; i < _countof(_BaseTable); i++)
                if (_BaseTable[i].key == code) return i;
            return 0;
        }
        constexpr base_t find_basecode(byte_t code) {
            for (auto i = 0; i < _countof(_BaseTable); i++)
                if (_BaseTable[i].code == code) return _BaseTable[i].key;
            return dynamic_base;
        }

        inline byte_t from_digit(byte_t digit, const char* begin, const char* end)
        {
            auto ptr = std::find(begin, end, digit);
            if (ptr == end) throw std::out_of_range("Provided digit is not in alphabet");
            return gsl::narrow<byte_t>(ptr - begin);
        }

        template <base_t _Base>
        byte_t from_digit(byte_t digit) {
            constexpr auto _Index = details::find_baseimpl(_Base);
            constexpr auto _Impl = details::_BaseTable[_Index];

            static_assert(_Index > 0, "base not implemented");

            const auto begin = _Impl.digits;
            const auto end = _Impl.digits + _Impl.radix;
            return from_digit(digit, begin, end);
        }
        template<>
        inline byte_t from_digit<base256>(byte_t digit) {
            return digit;
        }

        template <base_t _Base>
        byte_t to_digit(byte_t value) {
            constexpr auto _Index = details::find_baseimpl(_Base);
            constexpr auto _Impl = details::_BaseTable[_Index];

            static_assert(_Index > 0, "base not implemented");
            if (value >= _Impl.radix) throw std::out_of_range(std::string("Provided value is not in accepted range of ") + _Impl.name + ".");

            return _Impl.digits[value];
        }
        template<>
        inline byte_t to_digit<base256>(byte_t value) {
            return value;
        }

        template <base_t _FromBase, base_t _ToBase>
        buffer_t convert_base(bufferview_t from, const char*)
        {
            constexpr auto _FromIndex = details::find_baseimpl(_FromBase);
            constexpr auto _ToIndex = details::find_baseimpl(_ToBase);
            constexpr auto _FromImpl = details::_BaseTable[_FromIndex];
            constexpr auto _ToImpl = details::_BaseTable[_ToIndex];

            static_assert(_FromIndex > 0 && _ToIndex > 0, "base not implemented");

            // Skip & count leading zeroes.
            auto first = std::find_if_not(std::begin(from), std::end(from), [](auto v) { return v == 0; });
            auto last = std::end(from);

            const auto leadingZeroes = first - std::begin(from);
            const auto dataSize = last - first;

            // Compute the max size of the encoded string, maybe shorter (log(256) / log(10), rounded up).
            const auto tmpSize = size_t(dataSize * log((double)_FromImpl.radix) / log((double)_ToImpl.radix)) + 1;
            auto tmp = buffer_t(tmpSize, 0);

            // Process the data
            auto encodedLen = 0;
            while (first != last)
            {
                auto carry = static_cast<uint32_t>(from_digit<_FromBase>(*first++));

                auto i = 0;
                for (auto it = tmp.rbegin(); it != tmp.rend() && (carry || i < encodedLen); i++, it++)
                {
                    carry += _FromImpl.radix * (*it);
                    *it = gsl::narrow<byte_t>(carry % _ToImpl.radix);
                    carry /= _ToImpl.radix;
                }

                encodedLen = i;
            }

            // Skip leading zeroes in encoded buffer (keep at least one zero if null)
            auto tmpFirst = std::find_if_not(std::begin(tmp), std::end(tmp), [](auto v) { return v == 0; });
            if (tmpFirst == std::end(tmp)) tmpFirst--;

            // translate with digits
            const auto encodedSize = leadingZeroes + (std::end(tmp) - tmpFirst);
            auto encoded = buffer_t(encodedSize, to_digit<_ToBase>(0));
            for (auto i = leadingZeroes; i < encodedSize; ++i)
                encoded[i] = to_digit<_ToBase>(*tmpFirst++);

            return encoded;
        }

        template <base_t _Base = dynamic_base, int _Index = find_baseimpl(_Base)>
        class basecode_type
        {
        public:
            static_assert(_Base != dynamic_base, "dynamic_base is not allowed here");
            static_assert(_Index > 0, "The base_t is not implemented");

            constexpr basecode_type(base_t base) { Expects(base == _Base); }

            constexpr base_t      base()   const { return _BaseTable[_Index].key; }
            constexpr const char* name()   const { return _BaseTable[_Index].name; }
            constexpr uint32_t    code()   const { return _BaseTable[_Index].code; }
            constexpr int32_t     radix()  const { return _BaseTable[_Index].radix; }
            constexpr const char* digits() const { return _BaseTable[_Index].digits; }
            constexpr CodecFunc   encode() const { return _BaseTable[_Index].encode; }
            constexpr CodecFunc   decode() const { return _BaseTable[_Index].decode; }
            constexpr int         index()  const { return _Index; }
        };
        template <>
        class basecode_type<dynamic_base>
        {
        public:
            explicit constexpr basecode_type(base_t base) : _index(find_baseimpl(base)) { Expects(_index > 0); }

            constexpr base_t      base()   const { return _BaseTable[_index].key; }
            constexpr const char* name()   const { return _BaseTable[_index].name; }
            constexpr uint32_t    code()   const { return _BaseTable[_index].code; }
            constexpr int32_t     radix()  const { return _BaseTable[_index].radix; }
            constexpr const char* digits() const { return _BaseTable[_index].digits; }
            constexpr CodecFunc   encode() const { return _BaseTable[_index].encode; }
            constexpr CodecFunc   decode() const { return _BaseTable[_index].decode; }
            constexpr int         index()  const { return _index; }
        private:
            const int _index;
        };
    }

    //
    // The class `encoded_string<base_t>` is used to strongly type strings that are encoded in a base known at compile-time.
    // The specialization `encoded_string<dynamic_base>` is used to type encoded strings for which the base is not known at compile-time.
    //
    template <base_t _Base = dynamic_base>
    class encoded_string
    {
    public:
        // Construct empty encoded_string<>
        encoded_string() : _base(_Base) {}
        encoded_string(base_t code) : _base(code) {}

        // Construct by copying _Right
        //    - from string
        encoded_string(stringview_t _Right) : _string(to_string(_Right)), _base(_Base) {}
        encoded_string(base_t code, stringview_t _Right) : _string(to_string(_Right)), _base(code) {}

        encoded_string(const char* _Right) : encoded_string(gsl::ensure_z(_Right)) {}
        encoded_string(base_t code, const char* _Right) : encoded_string(code, gsl::ensure_z(_Right)) {}

        //    - from encoded_string<>
        encoded_string(const encoded_string<_Base>& _Right) : _string(_Right._string), _base(_Right.base()) {}
        template <base_t _RightBase>
        encoded_string(const encoded_string<_RightBase>& _Right) : _string(_Right._string), _base(_Right.base())
        {
            static_assert(_RightBase == _Base || _RightBase == dynamic_base || _Base == dynamic_base, "Mismatch between codes.");
        }

        // Construct by moving _Right
        //    - from string
        encoded_string(string_t&& _Right) : _string(std::move(_Right)), _base(_Base) {}
        encoded_string(base_t code, string_t&& _Right) : _string(std::move(_Right)), _base(code) {}

        //    - from encoded_string<>
        encoded_string(encoded_string<_Base>&& _Right) : _string(std::move(_Right._string)), _base(_Right.base()) {}
        template <base_t _RightBase>
        encoded_string(encoded_string<_RightBase>&& _Right) : _string(std::move(_Right._string)), _base(_Right.base())
        {
            static_assert(_RightBase == _Base || _RightBase == dynamic_base || _Base == dynamic_base, "Mismatch between codes.");
        }

        // Assign by copying _Right
        // Note: an encoded_string<> is immutable once initialized
        //    - from string
        encoded_string<_Base>& operator =(stringview_t _Right)
        {
            Expects(_string.empty() && !_Right.empty());
            _string = to_string(_Right);
            return *this;
        }
        encoded_string<_Base>& operator =(const char* _Right)
        {
            return operator+(gsl::ensure_z(_Right));
        }

        //    - from encoded_string<>
        encoded_string<_Base>& operator =(const encoded_string<_Base>& _Right)
        {
            Expects(_Right.base() == base());
            Expects(_string.empty() && !_Right.empty());
            _string = _Right._string;
            return *this;
        }
        template <base_t _RightBase>
        encoded_string<_Base>& operator =(const encoded_string<_RightBase>& _Right)
        {
            static_assert(_RightBase == _Base || _RightBase == dynamic_base || _Base == dynamic_base, "Mismatch between codes.");
            Expects(_Right.base() == base());
            Expects(_string.empty() && !_Right.empty());
            _string = _Right._string;
            return *this;
        }

        // Assign by moving _Right
        // Note: an encoded_string<> is immutable once initialized
        //    - from string
        encoded_string<_Base>& operator =(string_t&& _Right)
        {
            Expects(_string.empty() && !_Right.empty());
            _string = std::move(_Right);
            return *this;
        }

        //    - from encoded_string<>
        encoded_string<_Base>& operator =(encoded_string<_Base>&& _Right)
        {
            Expects(_Right.base() == base());
            Expects(_string.empty() && !_Right.empty());
            _string = std::move(_Right._string);
            return *this;
        }

        template <base_t _RightBase>
        encoded_string<_Base>& operator =(encoded_string<_RightBase>&& _Right)
        {
            static_assert(_RightBase == _Base || _RightBase == dynamic_base || _Base == dynamic_base, "Mismatch between codes.");
            Expects(_Right.base() == base());
            Expects(_string.empty() && !_Right.empty());
            _string = std::move(_Right._string);
            return *this;
        }

        base_t base() const { return _base.base(); }
        const string_t& str() const { return _string; }

        bool empty() const { return _string.empty(); }

    private:
        const details::basecode_type<_Base> _base;
        string_t _string;

        template <base_t _FriendBase> friend class encoded_string;
    };

    // operator+
    template <base_t _Base>
    string_t operator+(stringview_t _Left, const encoded_string<_Base>& _Right) { return (_Left + _Right.str()); }
    template <base_t _Base>
    string_t operator+(const encoded_string<_Base>& _Left, stringview_t _Right) { return (_Left.str() + _Right); }

    template <base_t _Base>
    string_t operator+(const char* _Left, const encoded_string<_Base>& _Right) { return (_Left + _Right.str()); }
    template <base_t _Base>
    string_t operator+(const encoded_string<_Base>& _Left, const char* _Right) { return (_Left.str() + _Right); }

    // operator==
    template <base_t _BaseLeft, base_t _BaseRight>
    bool operator==(const encoded_string<_BaseLeft>& _Left, const encoded_string<_BaseRight>& _Right)
    {
        return (_Left.base() == _Right.base()) && (_Left.str() == _Right.str());
    }

    template <base_t _Base>
    bool operator==(const std::string& _Left, const encoded_string<_Base>& _Right) { return (_Left == _Right.str()); }
    template <base_t _Base>
    bool operator==(const encoded_string<_Base>& _Left, const std::string& _Right) { return (_Left.str() == _Right); }

    // operator!=
    template <base_t _BaseLeft, base_t _BaseRight>
    bool operator!=(const encoded_string<_BaseLeft>& _Left, const encoded_string<_BaseRight>& _Right) { return !(_Left == _Right); }
    template <base_t _Base>
    bool operator!=(const std::string& _Left, const encoded_string<_Base>& _Right) { return !(_Left == _Right); }
    template <base_t _Base>
    bool operator!=(const encoded_string<_Base>& _Left, const std::string& _Right) { return !(_Left == _Right); }

    // operator<
    template <base_t _BaseLeft, base_t _BaseRight>
    bool operator<(const encoded_string<_BaseLeft>& _Left, const encoded_string<_BaseRight>& _Right)
    {
        if (_Left.base() < _Right.base()) return true;
        if (_Left.base() > _Right.base()) return false;

        return (_Left.str() < _Right.str());
    }
    template <base_t _Base>
    bool operator<(const std::string& _Left, const encoded_string<_Base>& _Right) { return (_Left < _Right.str()); }
    template <base_t _Base>
    bool operator<(const encoded_string<_Base>& _Left, const std::string& _Right) { return (_Left.str() < _Right); }

    //
    template <base_t _Base>
    std::ostream& operator<< (std::ostream& os, const encoded_string<_Base>& _Right) { return os << _Right.str(); }

    //
    // Encode a buffer/string into a _Base encoded_string
    //
    template <base_t _Base>
    encoded_string<_Base> encode(bufferview_t data)
    {
        constexpr auto _Index = details::find_baseimpl(_Base);
        constexpr auto _Impl = details::_BaseTable[_Index];

        static_assert(_Base != dynamic_base, "encode<_Base>(bufferview_t) is not allowed for _Base == basecode::dynamic_base");
        static_assert(_Index > 0, "encode<_Base>(bufferview_t) is not implementeed for this _Base");

        Expects(!data.empty());
        return as_string(_Impl.encode(data, _Impl.digits));
    }

    inline encoded_string<> encode(base_t base, bufferview_t data)
    {
        Expects(!data.empty());

        const auto index = details::find_baseimpl(base);
        const auto impl = details::_BaseTable[index];
        Expects(index > 0);

        return { base, as_string(impl.encode(data, impl.digits)) };
    }

    template <base_t _Base>
    encoded_string<_Base> encode(const std::string& string)                { return encode<_Base>(as_buffer(string)); }
    inline encoded_string<> encode(base_t code, const std::string& string) { return encode(code, as_buffer(string)); }

    //
    // Decode a _Base encoded_string into a buffer
    //
    template <base_t _Base>
    buffer_t decode(const encoded_string<_Base>& string)
    {
        constexpr auto _Index = details::find_baseimpl(_Base);
        constexpr auto _Impl = details::_BaseTable[_Index];

        static_assert(_Base != dynamic_base, "decode<_Base>(encoded_string) is not allowed for _Base == basecode::dynamic_base");
        static_assert(_Index > 0, "decode<_Base>(encoded_string) is not implementeed for this _Base");
        Expects(!string.empty());

        return _Impl.decode(as_buffer(string.str()), _Impl.digits);
    }

    inline buffer_t decode(const encoded_string<>& string)
    {
        Expects(!string.empty());

        const auto index = details::find_baseimpl(string.base());
        const auto impl = details::_BaseTable[index];
        Expects(index > 0);

        return impl.decode(as_buffer(string.str()), impl.digits);
    }

    inline stringview_t decode(base_t base, stringview_t src, buffer_t& dst)
    {
        Expects(!src.empty());

        const auto index = details::find_baseimpl(base);
        const auto impl = details::_BaseTable[index];
        Expects(index > 0);

        auto pos = to_string(src).find_first_not_of(impl.digits);
        if (pos == string_t::npos) pos = src.size();

        dst += decode({ base, src.first(pos) });

        return src.last(src.size() - pos);
    }
    inline buffer_t decode(base_t base, stringview_t src)
    {
        auto dst = buffer_t{};
        decode(base, src, dst);
        return dst;
    }
    inline buffer_t decode(base_t base, const char* src)
    {
        return decode(base, gsl::ensure_z(src));
    }

    //
    // Self-identifying base-encoded string
    //
    template <base_t _Base = dynamic_base>
    using multibase = std::string;

    template <base_t _Base>
    multibase<_Base> encode_multibase(bufferview_t data) {
        constexpr auto _Index = details::find_baseimpl(_Base);
        constexpr auto _Impl = details::_BaseTable[_Index];

        static_assert(_Base != dynamic_base, "make_multibase<_Base>(bufferview_t) is not allowed for _Base == basecode::dynamic_base");
        static_assert(_Index > 0, "make_multibase<_Base>(bufferview_t) is not implementeed for this _Base");

        return _Impl.code + encode<_Base>(data).str();
    }

    inline multibase<> encode_multibase(base_t base, bufferview_t data) {
        const auto index = details::find_baseimpl(base);
        const auto impl = details::_BaseTable[index];
        Expects(index > 0);
        return impl.code + encode(base, data).str();
    }

    template <base_t _Base>
    multibase<_Base> encode_multibase(const std::string& string)                { return encode_multibase<_Base>(as_buffer(string)); }
    inline multibase<> encode_multibase(base_t code, const std::string& string) { return encode_multibase(code, as_buffer(string)); }

    template <base_t _Base>
    encoded_string<_Base> decode_multibase(multibase<_Base> mb)        { return mb.substr(1); }
    inline encoded_string<> decode_multibase(const std::string& mb)    { return { details::find_basecode(mb[0]), mb.substr(1) }; }
}

inline multiformats::encoded_string<multiformats::base16> operator "" _16(const char* s, std::size_t)
{ return multiformats::encoded_string<multiformats::base16>(s); }

inline multiformats::encoded_string<multiformats::base64url> operator "" _64url(const char* s, std::size_t)
{
    return multiformats::encoded_string<multiformats::base64url>(s);
}

max_stars_repo_name: cedrou/cpp-multiformats

max_issues_repo_name: cedrou/multiformats

max_forks_repo_name: cedrou/multiformats


text:
/* filter/gsl_filter.h
 *
 * Copyright (C) 2018 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.
 */

#ifndef __GSL_FILTER_H__
#define __GSL_FILTER_H__

#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_movstat.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

/* end point handling methods */
typedef enum
{
  GSL_FILTER_END_PADZERO = GSL_MOVSTAT_END_PADZERO,
  GSL_FILTER_END_PADVALUE = GSL_MOVSTAT_END_PADVALUE,
  GSL_FILTER_END_TRUNCATE = GSL_MOVSTAT_END_TRUNCATE
} gsl_filter_end_t;

/* robust scale estimates */
typedef enum
{
  GSL_FILTER_SCALE_MAD, /* median absolute deviation */
  GSL_FILTER_SCALE_IQR, /* interquartile range */
  GSL_FILTER_SCALE_SN,  /* S_n scale statistic */
  GSL_FILTER_SCALE_QN   /* Q_n scale statistic */
} gsl_filter_scale_t;

/* workspace for Gaussian filter */
typedef struct
{
  size_t K;        /* window size */
  double *kernel;  /* Gaussian kernel, size K */
  gsl_movstat_workspace *movstat_workspace_p;
} gsl_filter_gaussian_workspace;

gsl_filter_gaussian_workspace *gsl_filter_gaussian_alloc(const size_t K);
void gsl_filter_gaussian_free(gsl_filter_gaussian_workspace * w);
int gsl_filter_gaussian(const gsl_filter_end_t endtype, const double alpha, const size_t order, const gsl_vector * x,
                        gsl_vector * y, gsl_filter_gaussian_workspace * w);
int gsl_filter_gaussian_kernel(const double alpha, const size_t order, const int normalize, gsl_vector * kernel);

/* workspace for standard median filter */
typedef struct
{
  gsl_movstat_workspace *movstat_workspace_p;
} gsl_filter_median_workspace;

gsl_filter_median_workspace *gsl_filter_median_alloc(const size_t K);
void gsl_filter_median_free(gsl_filter_median_workspace * w);
int gsl_filter_median(const gsl_filter_end_t endtype, const gsl_vector * x, gsl_vector * y, gsl_filter_median_workspace * w);

/* workspace for recursive median filter */
typedef struct
{
  size_t H;                            /* window half-length (K / 2) */
  size_t K;                            /* window size */
  void *state;                         /* workspace for min/max accumulator */
  double *window;                      /* array holding first window */
  const gsl_movstat_accum * minmaxacc; /* minimum/maximum accumulator */
  gsl_movstat_workspace *movstat_workspace_p;
} gsl_filter_rmedian_workspace;

gsl_filter_rmedian_workspace *gsl_filter_rmedian_alloc(const size_t K);
void gsl_filter_rmedian_free(gsl_filter_rmedian_workspace * w);
int gsl_filter_rmedian(const gsl_filter_end_t, const gsl_vector * x, gsl_vector * y, gsl_filter_rmedian_workspace * w);
int gsl_filter_rmedian2(const gsl_vector * x, gsl_vector * y, gsl_filter_rmedian_workspace * w);

typedef struct
{
  gsl_movstat_workspace *movstat_workspace_p;
} gsl_filter_impulse_workspace;

gsl_filter_impulse_workspace *gsl_filter_impulse_alloc(const size_t K);
void gsl_filter_impulse_free(gsl_filter_impulse_workspace * w);
int gsl_filter_impulse(const gsl_filter_end_t endtype, const gsl_filter_scale_t scale_type, const double t,
                       const gsl_vector * x, gsl_vector * y, gsl_vector * xmedian, gsl_vector * xsigma, size_t * noutlier,
                       gsl_vector_int * ioutlier, gsl_filter_impulse_workspace * w);

__END_DECLS

#endif /* __GSL_FILTER_H__ */


text:
#include "matrix.h"
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <math.h>

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>

/*************************** Complex Number *************************/

/*
typedef struct {
    double dat[2];
} Complex;
*/

Complex ComplexPolar(double r, double theta) {
    Complex z = {r * cos(theta), r * sin(theta)};
    return z;
}

void ComplexPrint(Complex z) {
    printf("%f + %fi\n", z.re, z.im);
}

Complex ComplexAdd(Complex z, Complex w) {
    Complex x = {z.re + w.re, z.im + w.im};
    return x;
}

Complex ComplexMul(Complex z, Complex w) {
    Complex x = {z.re*w.re - z.im*w.im, z.re*w.im + w.re*z.im};
    return x;
}

Complex ComplexMulReal(Complex z, double r) {
    Complex x = {z.re*r, z.im*r};
    return x;
}

double ComplexMag(Complex z) {
    return sqrt(ComplexMagSquare(z));
}

double ComplexMagSquare(Complex z) {
    return z.re*z.re + z.im*z.im;
}

/*************************** Bit Vector *************************/

/*
struct BitVector {
    unsigned int size;
    unsigned char* data;
};
*/

// alloc and init to 0
struct BitVector* newBitVector(unsigned int size) {
    assert(size > 0);

    struct BitVector* vec = (struct BitVector*)malloc(sizeof(struct BitVector));
    vec->size = size;

    vec->data = (unsigned char*)gsl_vector_alloc(vec->size);
    gsl_vector_set_zero((gsl_vector*)vec->data);

    return vec;
}

// alloc and init randomly
struct BitVector* newBitVectorRandom(unsigned int size) {
    assert(size > 0);

    struct BitVector* vec = (struct BitVector*)malloc(sizeof(struct BitVector));
    vec->size = size;
    vec->data = (unsigned char*)gsl_vector_alloc(vec->size);

    BitVectorSetRandom(vec);

    return vec;
}

void BitVectorSetRandom(struct BitVector* vec) {
    for (unsigned int i = 0; i < vec->size; i++) {
        BitVectorSet(vec, i, rand());
    }
}

void BitVectorSetZero(struct BitVector* vec) {
    gsl_vector_set_zero((gsl_vector*)vec->data);
}

void BitVectorFree(struct BitVector* vec) {
    gsl_vector_free((gsl_vector*)vec->data);
    free(vec);
}

// copy vec2 to vec1
void BitVectorCopy(struct BitVector* vec1, struct BitVector* vec2) {
    assert(vec1->size == vec2->size);

    gsl_vector_memcpy((gsl_vector*)vec1->data, (gsl_vector*)vec2->data);

}

unsigned int BitVectorGet(struct BitVector* vec, unsigned int loc) {
    assert(vec->size > loc);

    return (int)gsl_vector_get((gsl_vector*)vec->data, loc) % 2;
};

void BitVectorSet(struct BitVector* vec, unsigned int loc, unsigned int value) {
    assert(vec->size > loc);

    gsl_vector_set((gsl_vector*)vec->data, loc, (double) (value % 2));
}

void BitVectorFlip(struct BitVector* vec, unsigned int loc) {
    assert(vec->size > loc);

    BitVectorSet(vec, loc, 1-BitVectorGet(vec, loc));
}

int BitVectorSame(struct BitVector* vec1, struct BitVector* vec2) {
    assert(vec1->size == vec2->size);

    for (unsigned int i = 0; i < vec1->size; i++) {
        if (BitVectorGet(vec1, i) != BitVectorGet(vec2, i)) return 0;
    }
    return 1;
}

void BitVectorPrint(struct BitVector* vec) {
    printf("[");
    for (unsigned int i = 0; i < vec->size; i++) {
        printf("%d", BitVectorGet(vec, i));
    }
    printf("]\n");
}

unsigned int BitVectorInner(struct BitVector* vec1, struct BitVector* vec2) {
    assert(vec1->size == vec2->size);

    double inner;
    gsl_blas_ddot((gsl_vector*)vec1->data, (gsl_vector*)vec2->data, &inner);

    return (unsigned int)inner;
}

 // vec1 <- vec1 ^ vec2
void BitVectorXorSet(struct BitVector* vec1, struct BitVector* vec2) {
    assert(vec1->size == vec2->size);

    for (unsigned int i = 0; i < vec1->size; i++) {
        BitVectorSet(vec1, i, BitVectorGet(vec1,i) + BitVectorGet(vec2,i) % 2);
    }
}

// helper
void BitVectorMod2(struct BitVector* vec) {
    for (unsigned int i = 0; i < vec->size; i++) {
        BitVectorSet(vec, i, BitVectorGet(vec, i));
    }
}

/*************************** Bit Matrix *************************/

/*
struct BitMatrix {
    unsigned int rows;
    unsigned int cols;
    unsigned char* data;
};
*/

// alloc and init to 0
struct BitMatrix* newBitMatrixZero(unsigned int rows, unsigned int cols) {
    assert(rows > 0);
    assert(cols > 0);

    struct BitMatrix* mat = (struct BitMatrix*)malloc(sizeof(struct BitMatrix));
    mat->rows = rows;
    mat->cols = cols;

    mat->data = (unsigned char*)gsl_matrix_alloc(mat->rows, mat->cols);
    gsl_matrix_set_zero((gsl_matrix*)mat->data);

    return mat;
}

// alloc and init to identity
struct BitMatrix* newBitMatrixIdentity(unsigned int rows) {
    struct BitMatrix* mat = newBitMatrixZero(rows, rows);

    gsl_matrix_set_identity((gsl_matrix*)mat->data);

    return mat;
}

// alloc and init randomly
struct BitMatrix* newBitMatrixRandom(unsigned int rows, unsigned int cols) {
    assert(rows > 0);
    assert(cols > 0);

    struct BitMatrix* mat = newBitMatrixZero(rows, rows);
    BitMatrixSetRandom(mat);
    return mat;
}

void BitMatrixFree(struct BitMatrix* mat) {
    gsl_matrix_free((gsl_matrix*)mat->data);
    free(mat);
}

void BitMatrixCopy(struct BitMatrix* mat1, struct BitMatrix* mat2) {
    assert(mat1->rows == mat2->rows);
    assert(mat1->cols == mat2->cols);

    gsl_matrix_memcpy((gsl_matrix*)mat1->data, (gsl_matrix*)mat2->data);
}

void BitMatrixSetZero(struct BitMatrix* mat) {
    gsl_matrix_set_zero((gsl_matrix*)mat->data);
}

void BitMatrixSetRandom(struct BitMatrix* mat) {
    for (unsigned int i = 0; i < mat->rows; i++) {
        for (unsigned int j = 0; j < mat->cols; j++) {
            BitMatrixSet(mat, i, j, rand());
        }
    }
}

void BitMatrixSetIdentity(struct BitMatrix* mat) {
    gsl_matrix_set_identity((gsl_matrix*)mat->data);
}

void BitMatrixPrint(struct BitMatrix* mat) {
    printf("[");
    for (unsigned int i = 0; i < mat->rows; i++) {
        if (i != 0) printf("\n ");
        printf("[");
        for (unsigned int j = 0; j < mat->cols; j++) {
            printf("%d", BitMatrixGet(mat, i, j));
        }
        printf("]");
    }
    printf("]\n");

}

unsigned int BitMatrixGet(struct BitMatrix* mat, unsigned int row, unsigned int col) {
    assert(mat->cols > col);
    assert(mat->rows > row);

    return (int)gsl_matrix_get((gsl_matrix*)mat->data, row, col) % 2;
}

void BitMatrixSet(struct BitMatrix* mat, unsigned int row, unsigned int col, unsigned int value) {
    assert(mat->cols > col);
    assert(mat->rows > row);

    gsl_matrix_set((gsl_matrix*)mat->data, row, col, (double)(value % 2));
}

void BitMatrixFlip(struct BitMatrix* mat, unsigned int row, unsigned int col) {
    assert(mat->cols > col);
    assert(mat->rows > row);

    BitMatrixSet(mat, row, col, 1-BitMatrixGet(mat, row, col));
}

int BitMatrixSame(struct BitMatrix* mat1, struct BitMatrix* mat2) {
    assert(mat1->rows == mat2->rows);
    assert(mat1->cols == mat2->cols);

    for (unsigned int i = 0; i < mat1->rows; i++) {
        for (unsigned int j = 0; j < mat1->cols; j++) {
            if (BitMatrixGet(mat1, i, j) != BitMatrixGet(mat2, i, j)) return 0;
        }
    }
    return 1;
}

void BitMatrixXorSet(struct BitMatrix* mat1, struct BitMatrix* mat2) { // mat1 <- mat1 ^ mat2
    assert(mat1->cols == mat2->cols);
    assert(mat1->rows == mat2->cols);

    for (unsigned int i = 0; i < mat1->rows; i++) {
        for (unsigned int j = 0; j < mat1->cols; j++) {
            BitMatrixSet(mat1, i, j, BitMatrixGet(mat1,i,j) + BitMatrixGet(mat2,i,j) % 2);
        }
    }
}

void BitMatrixSetCol(struct BitMatrix* mat, struct BitVector* vec, unsigned int col) {
    assert(vec->size == mat->rows);
    assert(col < mat->cols);

    gsl_matrix_set_col((gsl_matrix*)mat->data, col, (gsl_vector*)vec->data);
}

void BitMatrixSetRow(struct BitMatrix* mat, struct BitVector* vec, unsigned int row) {
    assert(vec->size == mat->cols);
    assert(row < mat->rows);

    gsl_matrix_set_row((gsl_matrix*)mat->data, row, (gsl_vector*)vec->data);
}

struct BitVector* BitMatrixGetCol(struct BitMatrix* mat, unsigned int col) {
    assert(col < mat->cols);

    struct BitVector* vec = newBitVector(mat->rows);
    gsl_matrix_get_col((gsl_vector*)vec->data, (gsl_matrix*)mat->data, col);

    return vec;
}

struct BitVector* BitMatrixGetRow(struct BitMatrix* mat, unsigned int row) {
    assert(row < mat->rows);

    struct BitVector* vec = newBitVector(mat->cols);

    gsl_matrix_get_row((gsl_vector*)vec->data, (gsl_matrix*)mat->data, row);

    return vec;
}

void BitMatrixSwapCols(struct BitMatrix* mat, unsigned int col1, unsigned int col2) {
    gsl_matrix_swap_columns((gsl_matrix*)mat->data, col1, col2);
}

void BitMatrixSwapRows(struct BitMatrix* mat, unsigned int row1, unsigned int row2) {
    gsl_matrix_swap_rows((gsl_matrix*)mat->data, row1, row2);

}

// helper
void BitMatrixMod2(struct BitMatrix* mat) {
    for (unsigned int i = 0; i < mat->rows; i++) {
        for (unsigned int j = 0; j < mat->cols; j++) {
            BitMatrixSet(mat, i, j, BitMatrixGet(mat,i,j));
        }
    }
}

struct BitMatrix* BitMatrixTranspose(struct BitMatrix* mat) {
    struct BitMatrix* matT = newBitMatrixZero(mat->cols, mat->rows);

    gsl_matrix_transpose_memcpy((gsl_matrix*)matT->data, (gsl_matrix*)mat->data);

    return matT;
}

void BitMatrixTransposeSet(struct BitMatrix* mat) {
    struct BitMatrix* matT = BitMatrixTranspose(mat);

    mat->rows = matT->rows;
    mat->cols = matT->cols;
    BitMatrixCopy(matT, mat);

    BitMatrixFree(matT);
}

struct BitMatrix* BitMatrixMulMatrix(struct BitMatrix* mat1, struct BitMatrix* mat2) {
    assert(mat1->cols == mat2->rows);

    struct BitMatrix* matOut = newBitMatrixZero(mat1->rows, mat2->cols);
    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, (gsl_matrix*)mat1->data, (gsl_matrix*)mat2->data, 0, (gsl_matrix*)matOut->data);

    BitMatrixMod2(matOut);

    return matOut;
}

// mat2 <- mat1*mat2;
void BitMatrixMulMatrixLeft(struct BitMatrix* mat1, struct BitMatrix* mat2) {
    struct BitMatrix* matTmp = BitMatrixMulMatrix(mat1, mat2);
    BitMatrixCopy(mat2, matTmp);
    BitMatrixFree(matTmp);
}

// mat1 <- mat1*mat2;
void BitMatrixMulMatrixRight(struct BitMatrix* mat1, struct BitMatrix* mat2) {
    struct BitMatrix* matTmp = BitMatrixMulMatrix(mat1, mat2);
    BitMatrixCopy(mat1, matTmp);
    BitMatrixFree(matTmp);
}

struct BitVector* BitMatrixMulVector(struct BitMatrix* mat, struct BitVector* vec) {
    assert(mat->cols == vec->size);

    struct BitVector* vecOut = newBitVector(mat->rows);
    gsl_blas_dgemv(CblasNoTrans, 1, (gsl_matrix*)mat->data, (gsl_vector*)vec->data, 0, (gsl_vector*)vecOut->data);

    BitVectorMod2(vecOut);

    return vecOut;
}

// vec <- mat*vec
void BitMatrixMulVectorSet(struct BitMatrix* mat, struct BitVector* vec) {
    assert(mat->cols == vec->size);
    assert(mat->cols == mat->rows);

    gsl_blas_dgemv(CblasNoTrans, 1, (gsl_matrix*)mat->data, (gsl_vector*)vec->data, 0, (gsl_vector*)vec->data);
    BitVectorMod2(vec);
}

unsigned int BitMatrixRank(struct BitMatrix* mat) {
    struct BitMatrix* null = newBitMatrixIdentity(mat->cols);
    unsigned int rank = mat->cols; // rank of the null space

    struct BitMatrix* good = newBitMatrixZero(mat->cols, mat->cols);
    struct BitMatrix* bad = newBitMatrixZero(mat->cols, mat->cols);
    unsigned int numgood, numbad;

    for (unsigned int row = 0; row < mat->rows; row++) {
        struct BitVector* rowVec = BitMatrixGetRow(mat, row);
        numgood = 0;
        numbad = 0;

        for (unsigned int i = 0; i < rank; i++) {
            struct BitVector* nullVec = BitMatrixGetRow(null, i);

            if (BitVectorInner(nullVec, rowVec) % 2 == 0) {
                BitMatrixSetRow(good, nullVec, numgood);
                numgood += 1;
            } else {
                BitMatrixSetRow(bad, nullVec, numbad);
                numbad += 1;
            }
        }

        // copy good into null
        for (rank = 0; rank < numgood; rank++) {
            struct BitVector* goodVec = BitMatrixGetRow(good, rank);
            BitMatrixSetRow(null, goodVec, rank);
        }

        if (numbad > 1) {
            struct BitVector* flipVec = BitMatrixGetRow(bad, 0);

            for (; rank+1 < numbad+numgood; rank++) {
                struct BitVector* badVec = BitMatrixGetRow(bad, 1+rank-numgood);
                BitVectorXorSet(badVec, flipVec);
                BitMatrixSetRow(null, badVec, rank);
            }
        }
    }

    BitMatrixFree(null);
    BitMatrixFree(good);
    BitMatrixFree(bad);

    rank = mat->cols - rank; // rank of null space to rank of mat
    return rank;
}

max_stars_repo_name: patrickrall/CircuitSimulator

max_issues_repo_name: patrickrall/CircuitSimulator

max_forks_repo_name: patrickrall/CircuitSimulator


text:
#include <assert.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <float.h>
#include "nn.h"
#ifdef __APPLE__
#include <Accelerate/Accelerate.h>
#else
#include <cblas.h>
#endif

#define AVERAGE_FISHER 10

/**
 * A function for generating random numbers according to a N(mu, sigma) Gaussian
 * distribution.
 */
double randn (double mu, double sigma)
{
  double U1, U2, W, mult;
  static double X1, X2;
  static int call = 0;

  if (call == 1) {
      call = !call;
      return (mu + sigma * (double) X2);
  }
  do {
      U1 = -1 + ((double) rand () / RAND_MAX) * 2;
      U2 = -1 + ((double) rand () / RAND_MAX) * 2;
      W = pow (U1, 2) + pow (U2, 2);
  } while (W >= 1 || W == 0);

  mult = sqrt ((-2 * log (W)) / W);
  X1 = U1 * mult;
  X2 = U2 * mult;

  call = !call;

  return (mu + sigma * (double) X1);
}

/**
 * Shorthand to generate normally distributed random numbers.
 */
double rand_normal()
{
  return randn(0.0, 1.0);
}

/**
 * @brief Shuffle the index array.
 */
void shuffle_index(size_t num_pattern, size_t random_idx[num_pattern])
{
  size_t p, np, op;
  for (p = 0; p < num_pattern; ++p) {
    random_idx[p] = p;
  }
  for (p = 0; p < num_pattern; ++p) {
    np = p + ((double)rand()/((double)RAND_MAX+1)) * (num_pattern - 1 - p);
    op = random_idx[p];
    random_idx[p] = random_idx[np];
    random_idx[np] = op;
  }
}

/**
 * This function fills the input and target vector with the list of training
 * examples from automaton.
 */
void fill_input_target(size_t size, double* input, uint8_t* target,
                       uint8_t* automaton, int offset, int states)
{
  size_t index;
  int counter;
  int side = 2 * offset + 1;
  int num_input = states * (side * side - 1);
  uint8_t val;

  for (size_t i = 0; i < size; ++i) {
    for (size_t j = 0; j < size; ++j) {
      index = i * size + j;
      counter = 1;

      /* Add bias in the main vector */
      input[index * (num_input + 1)] = 1.0;

      for (int a = -offset; a < offset + 1; ++a) {
        for (int b = -offset; b < offset + 1; ++b) {
          if (a != 0 || b != 0) {  /* Don't take index i,j */
            val = automaton[((i + a + size) % size) * size
                            + ((j + b + size) % size)];
            for (uint8_t s = 0; s < states; ++s) {
              input[index * (num_input + 1) + counter] = (val == s) ? 1.: 0.;
              counter++;
            }
          }
        }
      }

      target[index] = automaton[i * size + j];
    }
  }
}

void init_weights(int num_hidden, int num_input, int num_output,
                  double* delta_w_ih, double* weight_ih,
                  double* delta_w_ho, double* weight_ho)
{
  int i, j, k;
  /* Initialize weights input -> hidden */
  for (j = 0; j < num_hidden; ++j) {
    delta_w_ih[j] = 0.0;
    weight_ih[j] = 0.0;
    for (i = 1; i < num_input + 1; ++i) {
      delta_w_ih[i * num_hidden + j] = 0.0;
      weight_ih[i * num_hidden + j] = rand_normal() *
        sqrt(1 / (double)(num_input));
    }
  }

  /* Initialize weights hidden -> output */
  for (k = 0; k < num_output; ++k) {
    delta_w_ho[k] = 0.0;
    weight_ho[k] = 0.0;
    for (j = 1; j < num_hidden + 1; ++j) {
      delta_w_ho[j * num_output + k] = 0.0;
      weight_ho[j * num_output + k] = rand_normal() *
        sqrt(1 / (double)(num_hidden));
    }
  }
}

void forward(double* input, double* output,
             int num_hidden, int num_pattern,
             int num_input, int num_output,
             double* hidden,
             double* hidden_bias,
             double* weight_ih,
             double* weight_ho)
{
  int j, k, p;
  double max_out, agg_out = 0.0;

  /* Compute hidden activations */
  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
              num_pattern, num_hidden, num_input + 1, 1.0,
              input, num_input + 1, weight_ih, num_hidden,
              0.0, hidden, num_hidden);

  /* ReLU non-linearity */
  for (p = 0; p < num_pattern; ++p) {
    hidden_bias[p * (num_hidden + 1)] = 1.0;
    for (j = 1; j < num_hidden + 1; ++j) {
      hidden_bias[p * (num_hidden + 1) + j] =
        (hidden[p * num_hidden + j - 1] > 0.0) ?
        hidden[p * num_hidden + j - 1]: 0.0;
    }
  }

  /* Compute output unit activations */
  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
              num_pattern, num_output, num_hidden + 1, 1.0,
              hidden_bias, num_hidden + 1, weight_ho, num_output,
              0.0, output, num_output);

  /* Compute softmax of output */
  for (p = 0; p < num_pattern; ++p) {
    max_out = output[p * num_output];
    for (k = 1; k < num_output; ++k) {
      if (output[p * num_output + k] > max_out) {
        max_out = output[p * num_output + k];
      }
    }
    agg_out = 0.0;
    for (k = 0; k < num_output; ++k) {
      output[p * num_output + k] = exp(output[p * num_output + k] - max_out);
      agg_out += output[p * num_output + k];
    }
    for (k = 0; k < num_output; ++k) {
      output[p * num_output + k] /= agg_out;
    }
  }
}

void update_weights(int num_input, int num_hidden, int num_output,
                    double eta,
                    double* batch_error, double reg, double alpha,
                    double* weight_ih, double* delta_w_ih,
                    double* delta_w_ih_prev,
                    double* weight_ho, double* delta_w_ho,
                    double* delta_w_ho_prev)
{
  int i, j, k;
  /* Update weights with average gradient */
  for (i = 0; i < num_input + 1; ++i) {
    for (j = 0; j < num_hidden; ++j) {
      if (reg > 0.) {
        *batch_error += 0.5 * reg *
          weight_ih[i * num_hidden + j] * weight_ih[i * num_hidden + j] ;
        weight_ih[i * num_hidden + j] += (1 + alpha) *
          reg * weight_ih[i * num_hidden + j];
      }

      weight_ih[i * num_hidden + j] -= (1 + alpha) *
        (eta * delta_w_ih[i * num_hidden + j]);

      if (alpha > 0.) {
        weight_ih[i * num_hidden + j] += (-alpha) *
          delta_w_ih_prev[i * num_hidden + j];
      }
    }
  }
  for (j = 0; j < num_hidden + 1; ++j) {
    for (k = 0; k < num_output; ++k) {
      if (reg > 0.) {
        *batch_error += 0.5 * reg *
          weight_ho[j * num_output + k] * weight_ho[j * num_output + k];
        weight_ho[j * num_output + k] += (1 + alpha) *
        reg * weight_ho[j * num_output + k];
      }

      weight_ho[j * num_output + k] -= (1 + alpha) *
        (eta * delta_w_ho[j * num_output + k]);

      if (alpha > 0.) {
        weight_ho[j * num_output + k] += (-alpha) *
          delta_w_ho_prev[j * num_output + k];
      }
    }
  }
}

void compute_batch_gradients(int base_index, double alpha,
                             int batch_size, int num_input, int num_output,
                             int num_hidden, size_t* random_idx,
                             double* output,
                             uint8_t* target, double* delta_output,
                             double* delta_w_ho, double* hidden_bias,
                             double* weight_ho, double* delta_h,
                             double* delta_w_ih, double* input,
                             double* delta_w_ih_prev,
                             double* delta_w_ho_prev,
                             network_opts_t* opts)
{
  int i, j, k, p;

  /* Gradient initialization (Nesterov momentum) */
  if (opts->optim_type == NESTEROV) {
    memcpy(delta_w_ih_prev, delta_w_ih,
           sizeof(double) * num_hidden * (num_input + 1));
    for (i = 0; i < num_input + 1; ++i) {
      for (j = 0; j < num_hidden; ++j) {
        delta_w_ih[i * num_hidden + j] *= alpha;
      }
    }

    memcpy(delta_w_ho_prev, delta_w_ho,
           sizeof(double) * num_output * (num_hidden + 1));
    for (j = 0; j < num_hidden + 1; ++j) {
      for (k = 0; k < num_output; ++k) {
        delta_w_ho[j * num_output + k] *= alpha;
      }
    }
  }

  for (int b = 0; b < batch_size; ++b) {
    p = random_idx[(base_index + b)];

    /* Backpropagation */
    for (k = 0; k < num_output; ++k) {
      /* Output gradients */
      delta_output[b * num_output + k] =
        (output[b * num_output + k] - ((k == target[p])? 1.0: 0.0));
      delta_output[b * num_output + k] /= (double) batch_size;
    }
  }

  /* Dot product hidden_bias x delta_output stored in delta_w_ho */
  cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans,
              num_hidden + 1, num_output, batch_size,
              1 , hidden_bias, num_hidden + 1,
              delta_output, num_output,
              (opts->optim_type == NESTEROV)? 1.0: 0.0,
              delta_w_ho, num_output);

  /* Dot product weight_ho x delta_output (we don't count the bias in
     weight_ho) stored in delta_h */
  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans,
              batch_size, num_hidden, num_output,
              1.0, delta_output, num_output,
              &weight_ho[num_output], num_output,
              0.0, delta_h, num_hidden);

  /* Hidden layer non linearity gradients */
  for (int b = 0; b < batch_size; ++b) {
    for (j = 0; j < num_hidden; ++j) {
      delta_h[b * num_hidden + j] *=
        ((hidden_bias[b * (num_hidden + 1) + j + 1] > 0) ? 1.0: 0.0);
    }
  }

  cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans,
              num_input + 1, num_hidden, batch_size,
              1 , input, num_input + 1,
              delta_h, num_hidden,
              (opts->optim_type == NESTEROV)? 1.0: 0.0,
              delta_w_ih, num_hidden);

}

double compute_loss(int base_index, int batch_size,
                    size_t* random_idx,
                    int num_output, double* output,
                    uint8_t* target)
{
  double batch_error = 0.0;
  int p;
  for (int b = 0; b < batch_size; ++b) {
    p = random_idx[(base_index + b)];
    /* Compute loss */
    batch_error +=
      - log((output[b * num_output + target[p]] > 0) ?
            output[b * num_output + target[p]]: DBL_MIN);
  }
  batch_error /= (double) batch_size;
  return batch_error;
}

double compute_fisher(int states, int neighbors,
                      int num_pattern, int num_output,
                      int num_hidden, int num_input,
                      double* input,
                      double* weight_ih, double* weight_ho)
{
  double fisher_information = 0.0;
  int base, index, delta;

  /* Build the secondary input that will be perturbed */
  double* perturbed_input =
    (double*) malloc(sizeof(double) * num_pattern * (num_input + 1));

  /* Allocate placeholders in the function to make it more
     adaptable to the input dataset. */
  double* output =
    (double*) malloc(sizeof(double) * num_pattern * num_output);
  double* output_pert =
    (double*) malloc(sizeof(double) * num_pattern * num_output);
  double* hidden =
    (double*) malloc(sizeof(double) * num_pattern * num_hidden);
  double* hidden_bias =
    (double*) malloc(sizeof(double) * num_pattern * (num_hidden + 1));

  /* Compute the output of the network on the base dataset */
  forward(input, output, num_hidden, num_pattern, num_input, num_output,
          hidden, hidden_bias, weight_ih, weight_ho);

  for (int n = 0; n < AVERAGE_FISHER; ++n) {
    for (int neighbor_n = 0; neighbor_n < neighbors; ++neighbor_n) {

      /* Perturb second input */
      memcpy(perturbed_input, input,
            sizeof(double) * num_pattern * (num_input + 1));
      for (int i = 0; i < num_pattern; ++i) {
        /* Choose index of the input cell to perturb */
        index = i * (num_input + 1) + neighbor_n * states;

        for (int t = 0; t < states; ++t) {
          base = 1 + (rand() % (states - 1));
          perturbed_input[index + ((base + t) % states)] =
            input[index + t];
        }
      }
      /* Compute the output of the network on the perturbed dataset */
      forward(perturbed_input, output_pert, num_hidden, num_pattern,
              num_input, num_output, hidden, hidden_bias, weight_ih, weight_ho);

      for (int i = 0; i < num_pattern; ++i) {
        for (int j = 0; j < num_output; ++j) {
          if (output_pert[i * num_output + j] > 0
              && output[i * num_output + j] > 0) {
            delta = log(output[i * num_output + j]
                        / output_pert[i * num_output + j]);
            fisher_information += output[i * num_output + j]
              * delta * delta;
          }
        }
      }

    }
  }

  free(output);
  free(hidden);
  free(hidden_bias);
  free(perturbed_input);
  free(output_pert);

  return fisher_information / (AVERAGE_FISHER * num_pattern);
}

double compute_error(int num_pattern, int num_output,
                     int num_hidden, int num_input,
                     double* input, uint8_t* target,
                     double* weight_ih, double* weight_ho)
{
  double test_error = 0.0;
  double val;

  /* Allocate placeholders in the function to make it more
     adaptable to the input dataset. */
  double* output =
    (double*) malloc(sizeof(double) * num_pattern * num_output);
  double* hidden =
    (double*) malloc(sizeof(double) * num_pattern * num_hidden);
  double* hidden_bias =
    (double*) malloc(sizeof(double) * num_pattern * (num_hidden + 1));

  /* Compute the output of the network */
  forward(input, output, num_hidden, num_pattern, num_input, num_output,
          hidden, hidden_bias, weight_ih, weight_ho);

  /* Compute loss */
  for (int p = 0; p < num_pattern; ++p) {
    val = output[p * num_output + target[p]];
    test_error += - log((val > 0) ? val: DBL_MIN);
  }
  test_error /= num_pattern;

  free(output);
  free(hidden);
  free(hidden_bias);

  return test_error;
}

void train_nn_on_automaton(size_t size, int states,
                           uint8_t* train_automaton,
                           uint8_t** test_automata,
                           int n_tests,
                           network_opts_t* opts,
                           network_result_t* res)
{
  if (opts->verbosity >= 1) {
    fprintf(stdout, "\nProcessing with options: h%i r%i e%i\n",
            opts->num_hid, opts->offset, opts->max_epoch);
  }

  size_t num_pattern = size * size;
  int side = 2 * opts->offset + 1;
  int num_input = states * (side * side - 1);
  int num_hidden = opts->num_hid;
  int num_output = states;

  double batch_error, error, eta = 1, alpha = 0.9;
  if (opts->optim_type != NESTEROV) {
    alpha = 0.;
  }

  const int batch_size = 8;
  /* Regression parameter */
  double reg = 0.;

  /* ====== Network and training variables declaration ====== */

  /* num_pattern x (num_input + 1) array that holds all the training set */
  double* base_input =
    (double *) malloc(num_pattern * (num_input + 1) * sizeof(double));
  /* Array that holds the training labels */
  uint8_t* target = (uint8_t *) malloc(num_pattern * sizeof(uint8_t));
  /* Fill those arrays with the automaton's content */
  fill_input_target(size, base_input, target,
                    train_automaton, opts->offset, states);

  /* Arrays for the test data */
  double* test_input =
    (double *) malloc(num_pattern * (num_input + 1) * sizeof(double));
  uint8_t* test_target = (uint8_t *) malloc(num_pattern * sizeof(uint8_t));

  /* Weights of the network */
  double weight_ih[(num_input + 1) * num_hidden];
  double weight_ho[(num_hidden + 1) * num_output];

  /* Gradients of the output and hidden layer placeholders */
  double delta_output[batch_size * num_output];
  double delta_h[batch_size * num_hidden];

  /* Weight gradients and previous gradients for Nesterov momentum */
  double delta_w_ih[(num_input + 1) * num_hidden];
  double delta_w_ho[(num_input + 1) * num_hidden];
  double* delta_w_ih_prev = NULL;
  double* delta_w_ho_prev = NULL;

  /* Allocate only if being used later */
  if (opts->optim_type == NESTEROV) {
    delta_w_ih_prev = malloc(sizeof(double) * (num_input + 1) * num_hidden);
    delta_w_ho_prev =  malloc(sizeof(double) * (num_hidden + 1) * num_output);
  }

  /* Allocate the arrays that will hold data for each batch */
  double input[(batch_size) * (num_input + 1)];
  double hidden[batch_size * num_hidden];
  double hidden_bias[batch_size * (num_hidden + 1)];
  double output[batch_size * num_output];

  int epoch;
  size_t random_idx[num_pattern];
  double test_error = 0.0;
  double test_var = 0.0;

  /* ===================== End declarations ====================== */

  /* Initialize the weights of the network */
  init_weights(num_hidden, num_input, num_output,
               delta_w_ih, weight_ih,
               delta_w_ho, weight_ho);

  for (epoch = 0; epoch < opts->max_epoch; ++epoch) {
    /* Initialize error */
    error = 0.0;

    /* Learning rate decay */
    if (epoch > 0 && epoch%10 == 0 && opts->decay == DECAY) {
      eta -= .5 * eta;
    }

    /* Random ordering of input patterns done at each epoch */
    shuffle_index(num_pattern, random_idx);

    /* Loop through every batch in the dataset */
    for (size_t s = 0; s < num_pattern; s += batch_size) {
      for (int b = 0; b < batch_size; ++b) {
        if (s + b >= num_pattern) {
          break;
        }
        /* Copy batch elements to the input array for processing */
        memcpy(input + b * (num_input + 1),
               base_input + (random_idx[(s + b)] *
                             (num_input + 1)),
               sizeof(double) * (num_input + 1));
      }

      /* Forward pass */
      forward(input, output, num_hidden, batch_size, num_input, num_output,
              hidden, hidden_bias, weight_ih, weight_ho);

      batch_error = compute_loss(s, batch_size, random_idx, num_output,
                                 output, target);

      /* Compute the gradients for each weight matrix */
      compute_batch_gradients(s, alpha, batch_size, num_input, num_output,
                              num_hidden,
                              random_idx, output, target,
                              delta_output, delta_w_ho, hidden_bias, weight_ho,
                              delta_h, delta_w_ih, input,
                              delta_w_ih_prev, delta_w_ho_prev, opts);

      /* Update the weights of the network */
      update_weights(num_input, num_hidden, num_output, eta, &batch_error, reg,
                     alpha, weight_ih, delta_w_ih, delta_w_ih_prev, weight_ho,
                     delta_w_ho, delta_w_ho_prev);

      /* Error is only updated here because it might have been changed by
         `update_weights` */
      error += batch_error * batch_size;
    }
    error /= (double)(num_pattern);

    if (opts-> verbosity >= 1 && epoch%5 == 0) {
      fprintf(stdout, "\nEpoch %d: Error = %f", epoch, error);
    }
  }

  /* Compute error on the training set */
  error = compute_error(num_pattern, num_output, num_hidden,
                        num_input, base_input, target,
                        weight_ih, weight_ho);

  if (opts->verbosity >= 1) {
    fprintf(stdout, "\nTrain error: %f", error);
  }

  /* Compute Fisher information if the flag requires it */
  if (opts->fisher == FISHER) {
    res->fisher_info = compute_fisher(states, (side*side - 1), num_pattern,
                                      num_output, num_hidden, num_input,
                                      base_input, weight_ih, weight_ho);
  }

  /* Was an array of states to test on provided ? */
  if (test_automata != NULL) {
    double test_errors[n_tests];
    for (int i = 0; i < n_tests; ++i) {
      /* Fill the placeholders with test data */
      fill_input_target(size, test_input, test_target,
                        test_automata[i], opts->offset, states);

      /* Compute error on the test set */
      test_errors[i] = compute_error(num_pattern, num_output, num_hidden,
                                     num_input, test_input, test_target,
                                     weight_ih, weight_ho);
      test_error += error / test_errors[i];
    }

    /* Average score across all tested states */
    test_error /= n_tests;

    /* Compute variance of scores */
    for (int i = 0; i < n_tests; ++i) {
      test_var += (error / test_errors[i] - test_error)
        * (error / test_errors[i] - test_error);
    }
    test_var /= n_tests;
    test_var = sqrt(test_var);

    if (opts->verbosity >= 1) {
      fprintf(stdout, "\tTest error: %f\tVar: %f\tRatio: %f",
              error / test_error, test_var, test_error);
    }
  }

  /* Log results */
  if (opts->verbosity >= 1 && opts->fisher== FISHER) {
      fprintf(stdout, "\tFisher: %f\n", res->fisher_info);
  } else if (opts->verbosity >= 1){
      fprintf(stdout, "\n");
  }

  /* Output data */
  res->train_error = error;
  res->test_error = error / test_error;
  res->error_var = test_var;

  /* Cleanup allocated arrays */
  free(target);
  free(base_input);
  free(test_input);
  free(test_target);
  free(delta_w_ih_prev);
  free(delta_w_ho_prev);
}


text:
////////////////////////////////////////////////////////////
//
// Copyright (c) 2018 Sir Ementaler
//
// 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.
//
////////////////////////////////////////////////////////////

#ifndef PICTOLEV_BINARY_SERIALIZATION_H
#define PICTOLEV_BINARY_SERIALIZATION_H

#include <algorithm>
#include <iterator>
#include <limits>
#include <ostream>
#include <type_traits>
#include <gsl/gsl_util>

enum class endian {
#ifdef _WIN32
	little = 0,
	big = 1,
	native = little
#else
	little = __ORDER_LITTLE_ENDIAN__,
	big = __ORDER_BIG_ENDIAN__,
	native = __BYTE_ORDER__
#endif
};

template<std::size_t N>
void write_buffer(std::ostream& stream, const char (&buffer)[N]) {
	stream.write(buffer, N);
}

template<endian Endian = endian::native, typename T>
std::enable_if_t<std::is_integral_v<T> && (sizeof(T) == 1)>
write_binary(std::ostream& stream, T value) {
	stream.put(value);
}

template<endian Endian, typename T>
std::enable_if_t<std::is_integral_v<T> && (sizeof(T) > 1)>
write_binary(std::ostream& stream, T value) {
	static_assert(Endian == endian::little || Endian == endian::big);
	std::make_unsigned_t<T> u_value = value;
	const auto generator = [&u_value]() {
		const char result = gsl::narrow_cast<unsigned char>(u_value);
		u_value >>= std::numeric_limits<unsigned char>::digits;
		return result;
	};
	char buffer[sizeof(T)] {};
	if constexpr (Endian == endian::little)
		std::generate(std::begin(buffer), std::end(buffer), generator);
	else
		std::generate(std::rbegin(buffer), std::rend(buffer), generator);
	write_buffer(stream, buffer);
}

#endif

max_stars_repo_name: SirEmentaler/PicToLev

max_issues_repo_name: SirEmentaler/PicToLev

max_forks_repo_name: SirEmentaler/PicToLev